subroutine update(curlev,nx,ny,nz,p)
use hpfprocessors
use param
integer :: curlev,nx,ny,nz
real, dimension(:,:,:,:) :: p

!HPF$ DISTRIBUTE p(*,BLOCK,*,*) onto PROC


real, parameter :: cx=1.0,cy=1.0
real  sx,sy,sx2,sy2,sxy,sx1,sy1,sxy1


 sx=cx*dtbydr
 sy=cy*dtbydr
 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)

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

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

!p(2,:,:,:)=p(1,:,:,:)

return
end
