subroutine update(anew,aold,t,dt,dr,nx,ny)
use parameters
!--------------------------------------------------------------------
!user function called from user supplied update driver: gridstep
!VERSION for TRANSPORT demo
!--------------------------------------------------------------------
real, dimension(:,:,:)   :: anew,aold
real                     :: t,dt
integer                  :: nx,ny
real, dimension(2)       :: dr

!HPF$ PROCESSORS P(PX,PY)
!HPF$ DISTRIBUTE (BLOCK,BLOCK,*) ONTO *P :: anew,aold

real, parameter :: cx=1.0,cy=1.0
real  sx,sy,sx2,sy2,sxy,sx1,sy1,sxy1
!write(6,*) 'start update'

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

!write(6,*) 'TIME=',t,nx,ny
!write(6,*) aold(5,1:10,1), '*'

do i=1,100
FORALL(i=1:nx,j=1:ny) anew(i,j,1)=exp(sin(aold(i,j,1)**2))
FORALL(i=1:nx,j=1:ny) anew(i,j,1)=anew(i,j,1)**3
enddo

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

!write(6,*) anew(5,1:10,1), '*'

return
end









