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

real(SP),parameter :: cx=1.0,cy=1.0
real(SP) sx,sy,sx2,sy2,sxy,sx1,sy1,sxy1
real(SP) coef1,coef2,coef3,coef4,coef5,coef6
write(56,*) '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)
! 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

write(56,*) dr,dt,t
write(56,*) nx,ny
write(56,*) 'AOLD'
write(56,*) aold(5,:), '*'

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

write(56,*) 'ANEW'
write(56,*) anew(5,:), '*'

return
end







