subroutine update(a,t,dr,nx,ny)
use parameters
real, dimension(:,:,:) :: a
real t
integer nx,ny
real, dimension(0:2) :: dr
!HPF$ PROCESSORS P(PX,PY)
!HPF$ DISTRIBUTE *(*,BLOCK,BLOCK) onto *P :: a

real,parameter :: cx=0.75,cy=0.75
real sx,sy,sx2,sy2,sxy,sx1,sy1,sxy1
real coef1,coef2,coef3,coef4,coef5,coef6

 sx=cx*dr(0)/dr(1)
 sy=cy*dr(0)/dr(2)
 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


do i=1,500
a(2,:,:)=cshift(a(2,:,:),shift=1,dim=2)
FORALL(i=1:nx,j=1:ny) a(2,i,j)=exp(sin(a(1,i,j)**2))
FORALL(i=1:nx,j=1:ny) a(2,i,j)=a(2,i,j)**3
enddo

t=t+dr(0)
a(2,:,:)=a(1,:,:)

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

!     inner boundary
      a(1,1,:)=a(1,2,:)
      a(1,:,1)=sin(0.0628*t)


return
end