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

!nx=size(a,dim=2)
!ny=size(a,dim=3)


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

!write(6,*) 'start update'
t=t+dr(0)
a(2,:,:)=a(1,:,:)

!write(6,*) 'new=old: done'

!     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 

!write(6,*) 'interior: done'

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

!write(6,*) 'inner boundary: done'

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

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

!write(6,*) 'outer boundaries: done'

!    a(1,nx,ny)=sxy1*(                                          &
!      -a(1,nx-1,ny-1)*(1.0-sx-sy)                               &
!       -a(1,nx-1,ny)*(1.0-sx+sy)-a(2,nx,ny-1)*(1.0+sx-sy)          &
!      +a(2,nx,ny)*(1-sx-sy)+a(2,nx-1,ny-1)*(1.0+sx+sy)               &
!      +a(2,nx,ny-1)*(1.0-sx+sy)+a(2,nx-1,ny)*(1+sx-sy))


write(6,*) 'time',t,' value=',a(1,2,2)

return
end







