subroutine update(phi_1,phi_2,phit_1,phit_2,phix_1,phix_2, &
                  phiy_1,phiy_2,phiz_1,phiz_2,dt,dx,dy,dz,nx,ny,nz)
!--------------------------------------------------------------------
!user function called from user supplied update driver: gridstep
!VERSION for SW
!--------------------------------------------------------------------

use types
real(SP), dimension(:,:,:) :: phi_1,phi_2,phit_1,phit_2,phix_1,phix_2
real(SP), dimension(:,:,:) :: phiy_1,phiy_2,phiz_1,phiz_2
real(SP)                   :: dt,dx,dy,dz
integer(I4B)               :: nx,ny,nz

real(SP)                   :: dtbydx,dtbydy,dtbydz,temp
real(SP)                   :: vist,visx,visy,visz,vis_coeff
integer(I4B)               :: i,j,k,ii,jj,kk

write(56,*) 'start update'

dtbydx = dt / dx 
dtbydy = dt / dy 
dtbydz = dt / dz 

vis_coeff = 1.0/100.0

do i=2,nx-1
do j=2,ny-1
do k=2,nz-1

!  CALCULATE VISCOSITY TERMS
   vist = (phix_1(i+1,j,k) - 2.0 * phix_1(i,j,k) + phix_1(i-1,j,k) ) &
        + (phiy_1(i,j+1,k) - 2.0 * phiy_1(i,j,k) + phiy_1(i,j-1,k) ) &
        + (phiz_1(i,j,k+1) - 2.0 * phiz_1(i,j,k) + phiz_1(i,j,k-1) )
   vist = vist * vis_coeff

   visx = phit_1(i+1,j,k) - 2.0 * phit_1(i,j,k) + phit_1(i-1,j,k)
   visx = visx * vis_coeff

   visy = phit_1(i,j+1,k) - 2.0 * phit_1(i,j,k) + phit_1(i,j-1,k)
   visy = visy * vis_coeff

   visz = phit_1(i,j,k+1) - 2.0 * phit_1(i,j,k) + phit_1(i,j,k-1)
   visz = visz * vis_coeff

!  temp = 0.0
!  do ii=i-1,i+1
!  do jj=j-1,j+1
!  do kk=k-1,k+1
!     temp = temp + phi_1(ii,jj,kk)
!  enddo
!  enddo
!  enddo
!  temp = temp/27.0

!  PHI update
phi_2(i,j,k) = phi_2(i,j,k) + 2.0 * dt * phit_1(i,j,k)

!  PHIT update
phit_2(i,j,k) = phit_2(i,j,k) + &
   (dtbydx) * (phix_1(i+1,j,k) - phix_1(i-1,j,k)) + &
   (dtbydy) * (phiy_1(i,j+1,k) - phiy_1(i,j-1,k)) + &
   (dtbydz) * (phiz_1(i,j,k+1) - phiz_1(i,j,k-1)) + vist

!  PHIX update
phix_2(i,j,k) = phix_2(i,j,k) + &
   (dtbydx) * (phit_1(i+1,j,k) - phit_1(i-1,j,k)) + visx

!  PHIY update
phiy_2(i,j,k) = phiy_2(i,j,k) + &
   (dtbydy) * (phit_1(i,j+1,k) - phit_1(i,j-1,k)) + visy

!  PHIZ update
phiz_2(i,j,k) = phiz_2(i,j,k) + &
   (dtbydz) * (phit_1(i,j,k+1) - phit_1(i,j,k-1)) + visz

enddo
enddo
enddo

!  phi swap: parallel implementation will need to be made better 

do i=2,nx-1
do j=2,ny-1
do k=2,nz-1
   temp = phi_1(i,j,k)
   phi_1(i,j,k) = phi_2(i,j,k)
   phi_2(i,j,k) = temp
   temp = phit_1(i,j,k)
   phit_1(i,j,k) = phit_2(i,j,k)
   phit_2(i,j,k) = temp
   temp = phix_1(i,j,k)
   phix_1(i,j,k) = phix_2(i,j,k)
   phix_2(i,j,k) = temp
   temp = phiy_1(i,j,k)
   phiy_1(i,j,k) = phiy_2(i,j,k)
   phiy_2(i,j,k) = temp
   temp = phiz_1(i,j,k)
   phiz_1(i,j,k) = phiz_2(i,j,k)
   phiz_2(i,j,k) = temp


enddo
enddo
enddo

return
end







