      program transport
      parameter(nx=25,ny=25)
      real, dimension(nx,ny) :: Gnew,G
       
! initialize G      
      Gnew=0.0
!      Gnew(:,1)=5.0
!      Gnew(1,:)=10.0

      write(11,1001) nx,ny
      write(11,1002) Gnew

 1001 Format(2I5)
 1002 Format(E15.5)

      cx=0.5
      cy=0.5
      dt=1.0
      dx=1.0
      dy=1.0
 
      sx=cx*dt/dx
      sy=cy*dt/dy
      sx2=sx**2
      sy2=sy**2
      sxy=sx*sy
      sx1=1.0/(1.0+sx)
      sy1=1.0/(1.0+sy)
      sxy1=1.0/(1.0+sx+sy)

      coef1=-sy/2.0+sxy/2.0+sy2
      coef2= sy/2.0-sxy/2.0+sy2
      coef3=sxy/2.0+sy/2.0

      coef4=-sx/2.0+sxy/2.0+sx2
      coef5= sx/2.0-sxy/2.0+sx2
      coef6=sxy/2.0+sx/2.0

! evolution
      nsteps=100

      do ns=1,nsteps

      G=Gnew

!     interior
      Forall(i=2:nx-1,j=2:ny-1) Gnew(i,j)=                         &
       G(i,j)                                                      &
       -(G(i+1,j)-G(i-1,j))*sx/2.0-(G(i,j+1)-G(i,j-1))*sy/2.0      &      
       +(G(i+1,j)-2.0*G(i,j)+G(i-1,j))*sx2/2.0                     &
       +(G(i,j+1)-2.0*G(i,j)+G(i,j-1))*sy2/2.0                     &
       +(G(i+1,j+1)-G(i+1,j-1)-G(i-1,j+1)+G(i-1,j-1))*sxy/4.0 

!     inner boundary
      Gnew(1,:)=Gnew(2,:)
      Gnew(:,1)=ns*0.1*sin(0.314*ns)
       

!     outer boundaries
       Forall(j=2:ny-1) Gnew(nx,j)=sx1*(                           &
        -Gnew(nx-1,j)*(1-sx) + G(nx-1,j)*(1+sx)                    &
        +G(nx,j+1)*coef1 + G(nx,j-1)*coef2                         & 
        +(G(nx-1,j-1)-G(nx-1,j+1))*coef3) 

       Forall(i=2:nx-1) Gnew(i,ny)=sy1*(                           &
        -Gnew(i,ny-1)*(1-sy) + G(i,ny-1)*(1+sy)                    &
        +G(i+1,ny)*coef4 + G(i-1,ny)*coef4                         & 
        +(G(i-1,ny-1)-G(i+1,ny-1))*coef6) 

       Gnew(nx,ny)=sxy1*(                                          &
        -Gnew(nx-1,ny-1)*(1.0-sx-sy)                               &
        -Gnew(nx-1,ny)*(1.0-sx+sy)-G(nx,ny-1)*(1.0+sx-sy)          &
        +G(nx,ny)*(1-sx-sy)+G(nx-1,ny-1)*(1.0+sx+sy)               &
        +G(nx,ny-1)*(1.0-sx+sy)+G(nx-1,ny)*(1+sx-sy))

      write(11,1001) nx,ny
      write(11,1002) Gnew

  
       enddo 
    
     end 




