c     domain-dcmp/heat/heat-2d_y
c     1-D domain decomposition in y-direction
c
      include "mpif.h"
      parameter (MAX_X=400,MAX_Y=200)
      integer my_rank,nproc,ierr,k,ny1,ny2
      real*8 h1,h2,r1,r2,ua(MAX_X,MAX_Y),ub(MAX_X,MAX_Y),dt,t
      real*8 va(MAX_X,MAX_Y),vb(MAX_X,MAX_Y),yl
c     
      call init(nproc,my_rank)
c
c     r=dt/h^2=0.25
c     fine grid nx2=8, r2=dt/h^2=0.25
      r2=0.25
      ny2=8
      ntoty2=ny2*nproc+2
      ntotx2=ntoty2
      yl=1.0
      h2=yl/(ntoty2-1)
      dt=r2*h2*h2

c     coarse grid nx1=4
      ny1=4
      ntoty1=ny1*nproc+2
      ntotx1=ntoty1
      yl=1.0
      h1=yl/(ntoty1-1)
      r1=dt/(h1*h1)
c
      call init_cn(nproc,my_rank,MAX_X,ny1,ntotx1,h1,ua)
      call init_cn(nproc,my_rank,MAX_X,ny2,ntotx2,h2,va)
c
c     4 time levels
      do k=1,4
         t=dt*dble(k-1)
c
c        coarse grid
         call bndry_cn(nproc,my_rank,MAX_X,ny1,ntotx1,ua,t)
         call comput(nproc,my_rank,MAX_X,ny1,ntotx1,r1,ua,ub)
         call exchange(nproc,my_rank,MAX_X,ny1,ntotx1,ub)
         call gather(k,nproc,my_rank,MAX_X,ny1,ntotx1,ub)
         do i=1,ntotx1
            do j=1,ny1+2
               ua(i,j)=ub(i,j)
            enddo
         enddo
c
c        fine grid
         call bndry_cn(nproc,my_rank,MAX_X,ny2,ntotx2,va,t)
         call comput(nproc,my_rank,MAX_X,ny2,ntotx2,r2,va,vb)
         call exchange(nproc,my_rank,MAX_X,ny2,ntotx2,vb)
         call gather(k,nproc,my_rank,MAX_X,ny2,ntotx2,vb)
         do i=1,ntotx2
            do j=1,ny2+2
               va(i,j)=vb(i,j)
            enddo
         enddo
      enddo
      call MPI_Finalize(ierr)
      stop
      end
c
c--------------------------------------------------------
c     Initial Condition
c--------------------------------------------------------
      subroutine init_cn(nproc,my_rank,MAX_X,ny,ntotx,h,u1)
      real*8 u1(MAX_X,1),h
      real*8 x(100),y(100),pi
      integer nproc,my_rank,ny
      pi=dacos(-1.d0)
      do j=1,ny+2
         y(j)=my_rank*dble(ny)*h+dble(j-1)*h
      enddo
      do i=1,ntotx
         x(i)=h*(i-1)
      enddo
c
c     initial condition
      do i=1,ntotx
         do j=1,ny+2
            u1(i,j)=dsin(pi*x(i))*dsin(pi*y(j))
         enddo
      enddo
      return
      end
c
c----------------------------------------------------------
c     Boundary Condition
c----------------------------------------------------------
      subroutine bndry_cn(nproc,my_rank,MAX_X,ny,ntotx,u,t)
      integer nproc,my_rank,ny,ntotx
      real*8 u(MAX_X,1),t
c
c     bndry u(x,0)=sin(t)
c
      if(my_rank.eq.0) then
         do i=1,ntotx
            u(i,1)=dsin(t)
         enddo
      endif

c     bndry u(x,1)=sin(t)
      if(my_rank.eq.nproc-1) then
         do i=1,ntotx
            u(i,ny+2)=dsin(t)
         enddo
      endif

c     bndry u(0,y)=sin(t)
      do j=1,ny+2
         u(1,j)=dsin(t)
      enddo

c     bndry u(1,y)=sin(t)
      do j=1,ny+2
         u(ntotx,j)=dsin(t)
      enddo
      return
      end
c
c-----------------------------------------------------------
c     Exchangeing the first and the last u to the neighbors
c-----------------------------------------------------------
      subroutine exchange(nproc,my_rank,MAX_X,ny,ntotx,u)
      include "mpif.h"
      real*8 u(MAX_X,1)
      integer status(MPI_STATUS_SIZE)
      integer nproc,my_rank,itag,nproc1,b_nbr,t_nbr
      nproc1=nproc-1
c
c     sending the first u to my_rank-1 for 0< my_rank < nproc-1
c
      if((my_rank-1).lt.0) then
         b_nbr=MPI_PROC_NULL
      else
         b_nbr=my_rank-1
      endif
      if((my_rank+1).gt.nproc1) then
         t_nbr=MPI_PROC_NULL
      else
         t_nbr=my_rank+1
      endif
      itag=20
c
c     send u(i,2) to bottom neighbor and receive u(i,n+2) from top neighbor
c
      call MPI_Sendrecv(u(1,2),ntotx,MPI_DOUBLE_PRECISION,b_nbr,
     *     itag,u(1,ny+2),ntotx,MPI_DOUBLE_PRECISION,t_nbr,itag,
     *     MPI_COMM_WORLD,status,ierr)
c     
c     sending the last u to my_rank+1 for 0< my_rank < nproc-1
c     
      itag=10
c
c     send u(i,ny+1) to top neighbor and receive u(i,1) from bottom neighbor
c
      call MPI_Sendrecv(u(1,ny+1),ntotx,MPI_DOUBLE_PRECISION,t_nbr,
     *     itag,u(1,1),ntotx,MPI_DOUBLE_PRECISION,b_nbr,itag,
     *     MPI_COMM_WORLD,status,ierr)
      return
      end
c
c--------------------------------------------------------------
c     Computing u for the next time step
c--------------------------------------------------------------
      subroutine comput(nproc,my_rank,MAX_X,ny,ntotx,r,u1,u2)
      real*8 u1(MAX_X,1),u2(MAX_X,1),r
      integer nproc,my_rank
      do i=1,ntotx
         do j=1,ny+2
            u2(i,j)=u1(i,j)
         enddo
      enddo
      do i=2,ntotx-1
         do j=2,ny+1
            u2(i,j)=r*(u1(i-1,j)-2.*u1(i,j)+u1(i+1,j)+
     *           u1(i,j-1)-2.*u1(i,j)+u1(i,j+1))+u1(i,j)
         enddo
      enddo
      return
      end
c
      subroutine gather(k,nproc,my_rank,MAX_X,ny,ntotx,u)
      integer k,iroot,nproc,my_rank,ny,ntotx
      real*8 u(MAX_X,1),uu(200),array(400)
      include 'mpif.h'
c     collocting u of each proc to global array
      iroot=0
      do i=2,ntotx-1
         do j=1,ny+2
            uu(j)=u(i,j)
         enddo
         call MPI_Gather(uu(2),ny,MPI_DOUBLE_PRECISION,array,ny,
     *        MPI_DOUBLE_PRECISION,iroot,MPI_COMM_WORLD,ierr)
         if(my_rank.eq.0) then
            write(*,*) 'k=',k,' i=',i,' ',(array(j),j=1,ny*nproc)
         endif
      enddo
      return
      end
c-----------------------------------------------------------
c     Initializing MPI
c-----------------------------------------------------------
      subroutine init(nproc,my_rank)

      include "mpif.h"
      integer my_rank,nproc
      character*30 nodename
      integer nchar
      call MPI_init(ierr)
      call MPI_Comm_rank(MPI_COMM_WORLD,my_rank,ierr)
      call MPI_Comm_size(MPI_COMM_WORLD,nproc,ierr)
      call MPI_Get_processor_name(nodename,nchar,ierr)
      write(*,*) 'my_rank=',my_rank,' processor=',nodename
      return
      end