!
!################################################################
!################################################################
!#####                                                      #####
!#####               SUBROUTINE ASSIMVEL                    #####
!#####                                                      #####
!#####               Copyright (c) 1993                     #####
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!#####                                                      #####
!################################################################
!################################################################
!


SUBROUTINE assimvel(nx,ny,nz,x,y,z,zp,                                  & 1,21
           u,v,w,wcont,ptprt,pprt,qv,qc,qr,qi,qs,qh,                    &
           ubar,vbar,ptbar,pbar,rhostr,qvbar,j1,j2,j3,                  &
           tem1,tem2,tem3,tem4,tem5,tem6,tem7,tem8,tem9,                &
           tem10,tem11,tem12,tem13)
!
!---------------------------------------------------------------------
!
!  PURPOSE:
!
!  This subroutine solves a 3-D elliptic equation with
!  variable coeficients:
!
!  c1*d2q/dx2 + c2*d2q/dy2 + c3*d2q/dz2 - c4*d2q/dxdy - c5*d2q/dxdz
!  - c6*d2q/dydz - c7*dq/dx - c8*dq/dy - c9*dq/dz = 'Known Forcing'
!
!  subject to the Dirichlet boundary condition:
!
!  q = 0 on the x & y boundaries, and
!
!  subject to the z-boundary condition:
!
!  dq/dz = (z/x**2+y**2)*[2r(vro - vrm) + (xdq/dx + ydq/dy)]
!  where:    vro - observed radial velocity
!            vrm - model derived radial velocity
!
!  The solution is obtained by applying a direct SOR technique
!  (see Haltiner and Williams "Numerical Prediction and Dynamical
!  Weather Prediction", Ch. 5 pg 159.)
!
!---------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Steve Lazarus and Alan Shapiro
!  9/29/93
!
!  MODIFICATION HISTORY:
!
!  Limin Zhao (10/18/95)
!  Modified the code to obtain the 'lamda2' and 'lamda1' simultaneously
!  by solving the two related equations interatively. It helps to
!  eliminate the sensitivity of the solution to the different
!  numerical schemes.
!
!  Alan Shapiro and Limin Zhao (10/19/95)
!  Put 'lamda1' in new staggered points and modified the related code.
!  The above two modifications work together to ensure the divergence
!  of the wind field be small enough after the variational adjustment.
!
!  Limin Zhao and Alan Shapiro(11/10/95)
!  Put a 'patch' on the top boundary to avoid the singularity problem.
!  Inside the 'patch', the zero gradient B.C. is used to replace the
!  non-zero gradient B.C.
!
!  Limin Zhao (11/16/95)
!  Set an option for scaling 'lamda1'. When the coefficient 'scoef'
!  is set to 1, no scaling is performed.
!
!  Limin Zhao (1/11/96)
!  Strip out the hole-filling and time interpolation part, and make
!  the program mainly to the velocity adjustment, either use Vr or the
!  variational adjustment.
!
!  Alan Shapiro and Limin Zhao (03/25/96)
!  The boundary values in the patch are calculated linearly from
!  its latent boundary values.
!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    nx       Number of grid points in the x-direction (east/west)
!    ny       Number of grid points in the y-direction (north/south)
!    nz       Number of grid points in the z-direction (vertical)
!
!    x        x coordinate of grid points in physical/comp. space (m)
!    y        y coordinate of grid points in physical/comp. space (m)
!    z        z coordinate of grid points in computational space (m)
!    zp       Vertical coordinate of grid points in physical space (m)
!
!    u        x component of velocity (m/s)
!    v        y component of velocity (m/s)
!    w        Vertical component of velocity in Cartesian
!             coordinates (m/s)
!    wcont    Contravariant vertical velocity (m/s)
!
!    ptprt    Perturbation potential temperature (K)
!    pprt     Perturbation pressure (Pascal)
!    qv       Water vapor specific humidity (kg/kg)
!    qc       Cloud water mixing ratio (kg/kg)
!    qr       Rainwater mixing ratio (kg/kg)
!    qi       Cloud ice mixing ratio (kg/kg)
!    qs       Snow mixing ratio (kg/kg)
!    qh       Hail mixing ratio (kg/kg)
!
!    ubar     Base state x velocity component (m/s)
!    vbar     Base state y velocity component (m/s)
!    ptbar    Base state potential temperature (K)
!    pbar     Base state pressure (Pascal)
!    rhostr   Base state air density (kg/m**3) times j3
!    qvbar    Base state water vapor specific humidity (kg/kg)
!
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!
!    tem1     Observed radial velocity
!    tem5     Model or retrieved radial velocity
!
!---------------------------------------------------------------------
!
!  OUTPUT:
!
!    tem7     Solution of the Poisson equation (q above).
!
!---------------------------------------------------------------------
!
!  WORK ARRAYS:
!
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!    tem6     Temporary work array.
!    tem8     Temporary work array.
!    tem9     Temporary work array.
!
!    (These arrays are defined and used locally (i.e. inside this
!    subroutine), they may also be passed into routines called by
!    this one. Exiting the call to this subroutine, these temporary
!    work arrays may be used for other purposes, and therefore their
!    contents may be overwritten. Please examine the usage of work
!    arrays before you alter the code.)
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE     ! Force explicit declarations
!
!-----------------------------------------------------------------------
!
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'      ! Model physical constants
  INCLUDE 'globcst.inc'     ! Model control constants
  INCLUDE 'assim.inc'       ! Assim/Retr control parameters
  INCLUDE 'bndry.inc'       ! Boundary condition parameters
  INCLUDE 'grid.inc'
!
!-----------------------------------------------------------------------
!
  INTEGER :: nt                ! The no. of t-levels of t-dependent arrays
  INTEGER :: tpast             ! Index of time level for the past time.
  INTEGER :: tpresent          ! Index of time level for the present time.
  INTEGER :: tfuture           ! Index of time level for the future time.

  INTEGER :: in,jn,ir0,jr0,ipe,ips,jps,jpe,ipatch
  REAL :: acof,bcof

  PARAMETER (nt=3, tpast=1, tpresent=2, tfuture=3)

  INTEGER :: nx, ny, nz        ! Number of grid points in x, y and z directions

  INTEGER :: grdbas

  REAL :: x     (nx)           ! The x-coord. of the physical and
                               ! computational grid. Defined at u-point.
  REAL :: y     (ny)           ! The y-coord. of the physical and
                               ! computational grid. Defined at v-point.
  REAL :: z     (nz)           ! The z-coord. of the computational grid.
                               ! Defined at w-point on the staggered grid.
  REAL :: zp    (nx,ny,nz)     ! The physical height coordinate defined at
                               ! w-point of the staggered grid.

  REAL :: u     (nx,ny,nz,nt)  ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz,nt)  ! Total v-velocity (m/s)
  REAL :: w     (nx,ny,nz,nt)  ! Total w-velocity (m/s)
  REAL :: wcont (nx,ny,nz)     ! Contravariant vertical velocity (m/s)

  REAL :: ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
  REAL :: pprt  (nx,ny,nz,nt)  ! Perturbation pressure (Pascal)
  REAL :: qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg)
  REAL :: qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (kg/kg)

  REAL :: qr    (nx,ny,nz,nt)  ! Rain water mixing ratio (kg/kg)
  REAL :: qi    (nx,ny,nz,nt)  ! Cloud ice mixing ratio (kg/kg)
  REAL :: qs    (nx,ny,nz,nt)  ! Snow mixing ratio (kg/kg)
  REAL :: qh    (nx,ny,nz,nt)  ! Hail mixing ratio (kg/kg)

  REAL :: ubar  (nx,ny,nz)     ! Base state u-velocity (m/s)
  REAL :: vbar  (nx,ny,nz)     ! Base state v-velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
  REAL :: rhostr(nx,ny,nz)     ! Base state air density (kg/m**3) times j3
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                               ! (kg/kg)

  REAL :: j1    (nx,ny,nz)     ! Coordinate transformation Jacobian -d(zp)/d(x)
  REAL :: j2    (nx,ny,nz)     ! Coordinate transformation Jacobian -d(zp)/d(y)
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian  d(zp)/d(z)

  REAL :: tem1(nx,ny,nz)
  REAL :: tem2(nx,ny,nz)
  REAL :: tem3(nx,ny,nz)
  REAL :: tem4(nx,ny,nz)
  REAL :: tem5(nx,ny,nz)
  REAL :: tem6(nx,ny,nz)
  REAL :: tem7(nx,ny,nz)
  REAL :: tem8(nx,ny,nz)
  REAL :: tem9(nx,ny,nz)
  REAL :: tem10(nx,ny,nz)
  REAL :: tem11(nx,ny,nz)
  REAL :: tem12(nx,ny,nz)
  REAL :: tem13(nx,ny,nz)

!
!-----------------------------------------------------------------------
!
!  Misc. local variables:
!
!-----------------------------------------------------------------------
!
  INTEGER :: icount ! Parameter input flag
  SAVE icount
  DATA icount /1/

  INTEGER :: i, j, k, ir, jr
  INTEGER :: istor,jstor,kstor
  INTEGER :: nchout    ! I/O Channel of unformatted binary restart data
  INTEGER :: it,ix
  INTEGER :: itmax     ! Maximum number of iterations
  INTEGER :: istat
  INTEGER :: isrc
  INTEGER :: tim
  INTEGER :: idummy
  INTEGER :: itag      ! Counter indicating input file date/time
  INTEGER :: count
  INTEGER :: mdpt

  REAL :: t1           ! Time of first data file
  REAL :: t2           ! Time of second data file
  REAL :: t3           ! Time of third data file

  REAL :: assimtim(100)! Time of (all) input data files

  REAL :: xs           ! Scalar coordinate x-component
  REAL :: ys           ! Scalar coordinate y-component
  REAL :: zs           ! Scalar coordinate z-component
  REAL :: xs1          ! Scalar coordinate x-component
  REAL :: ys1          ! Scalar coordinate y-component
  REAL :: zs1          ! Scalar coordinate z-component

  REAL :: xmove        ! (x,y) position of radar relative to
  REAL :: ymove        ! moving grid (grid origin at (0,0,0))

  REAL :: znz2         ! = z(nz-2) - zshift
  REAL :: znz1         ! = z(nz-1) - zshift
  REAL :: z2           ! = z(2)    - zshift

  REAL :: alpha        ! Over-relaxation coefficient
  REAL :: tol          ! An error tolerance, proportional to grid spacing
                       ! and max velocity error
  REAL :: veltol       ! Maximum tolerable error in velocity
  REAL :: newit        ! New iterate of q at a point
  REAL :: stor         ! Stores the new value of q temporarily
  REAL :: resid        ! the residual of the Poisson equation
  REAL :: diff         ! Difference between the updated q and old q at a point
  REAL :: test         ! Value of the maximum difference (qnew-qold) in grid
  REAL :: ptol         ! Error tolerance for hole-filling
  REAL :: anorm
  REAL :: anormf
  REAL :: test1
  REAL :: iterr

  REAL :: b1,b2,b3,b4  ! Coefficients of the Poisson equation.
  REAL :: b5,b6,b7,b8
  REAL :: b9
  REAL :: c1,c2,c3,c4  ! Coefficient combinations of b1...b9 in the Poisson
  REAL :: c5,c6,c7,c8  ! equation. The coefficients are a function of grid
  REAL :: c9,c10       ! spacing and location.
  REAL :: tm1,tm2,tm3
  REAL :: tm4,tm5,tm6
  REAL :: tm7,tm8,tm9
  REAL :: scoef
!ccxxx   You might need change veltol sometime
!   parameter(itmax=100000,veltol=0.01,scoef=1.0)  ! original
!   parameter(itmax=900000,veltol=0.01,scoef=1.0)  ! original
  PARAMETER(itmax=900000, veltol=0.5,scoef=1.0)
!   parameter(itmax=2, veltol=0.05,scoef=1.0)  ! only for test
  REAL :: rad          ! Radius of sphere in Cartesian coord
  REAL :: radh         ! Radius of circle in hz. plane
  REAL :: radinv       ! Inverse radius magnitude

  REAL :: rdx2         ! Reciprocal of (dx)**2
  REAL :: rdy2         ! Reciprocal of (dy)**2
  REAL :: rdz2         ! Reciprocal of (dz)**2

  REAL :: term1,term2  ! Analytic terms representing the RHS of the
  REAL :: term3,term4  ! Poisson eqn.
  REAL :: term5

  REAL :: vphi,vtheta  ! Model tangential and polar velocities
  REAL :: rcos         ! Squareroot of x**2+y**2

  INTEGER :: newebc    ! East boundary condition parameter.
  INTEGER :: newwbc    ! West boundary condition parameter.
  INTEGER :: newnbc    ! North boundary condition parameter.
  INTEGER :: newsbc    ! South boundary condition parameter.

  REAL :: vro,vrm   ! Used for diagnostics/testing
  REAL :: vradj
  REAL :: fpat      ! temp variable for forcing "patch"
  REAL :: maxdiv    ! Maximum divergence (diagnostic calculation)
  INTEGER :: imx,jmx   ! Gridpoint location of max divp
  INTEGER :: kmx
  INTEGER :: imx1,jmx1,kmx1
  REAL :: maxw
  REAL :: rrr       ! used for (x**2+y**2+z**2)**0.5
!
!-----------------------------------------------------------------------
!
!  Routines called:
!
!-----------------------------------------------------------------------
!
  EXTERNAL aamult
  EXTERNAL assimrd
  EXTERNAL retrint
  EXTERNAL avgx
  EXTERNAL avgy
  EXTERNAL avgz
  EXTERNAL dtadump
  EXTERNAL gtdmpfn
  EXTERNAL bcu
  EXTERNAL bcv
  EXTERNAL rhouvw
  EXTERNAL lbcw
  EXTERNAL vbcw
  EXTERNAL smth
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code ...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!
!-----------------------------------------------------------------------
!
!
!  The following is a brief summary of the assimilation flags
!  available within this routine.  For a more detailed review
!  of the options see ASSIMCON.F.
!
!  ivar    Variational adjustment flag:
!
!          = 0, NO variational adjustment this model time step
!          = 1, Perform variational adjustment this model time step
!
!  insrt   Insertion flag:
!
!          = 0, Do NOT insert velocities this time step
!          = 1, Insert velocities this time step
!
!  A matrix describing the assimilation options herein:
!
!       ivar   insrt   Comments
!       --------------------------------------------------------
!        0       0     Nothing - Exit assimilation mode
!        0       1     Direct insertion only
!        1       0     Variational adjustment only
!        1       1     Variational adjustment + direct insertion
!
!  NOTE: If irecov=1, then the ingest and time interpolation of
!       3 input data files occurs herein.
!
!-----------------------------------------------------------------------
!
  WRITE(6,*)'code in assimvel'

  IF ( insrt == 0.AND.ivar == 0.) RETURN
!
!-----------------------------------------------------------------------
!
!  Has switch been set for a direct insertion only? If yes, ingest
!  data.
!
!  NOTE: Insertion of data is done here ONLY if there is no
!        variational adjustment (ivopt=0). If ivopt=1, then the
!        adjusted data is inserted at the end of this routine.
!
!-----------------------------------------------------------------------
!
  tim = tpresent

  umove=0.0
  vmove=0.0
  xmove= xshift - umove*(curtim-assimtim(1))
  ymove= yshift - vmove*(curtim-assimtim(1))
!
!-------------------------------------------------------------------------
!
!  Interpolate input 'model' velocity components to scalar grid point.
!  Where:
!              tem2 = scalar u component
!              tem3 = scalar v component
!              tem4 = scalar w component
!
!-------------------------------------------------------------------------
!
  DO k = 1,nz
    DO j = 1,ny
      DO i = 1,nx
        tem7(i,j,k) = u(i,j,k,tim)
        tem8(i,j,k) = v(i,j,k,tim)
        tem9(i,j,k) = w(i,j,k,tim)
      END DO
    END DO
  END DO

!

  CALL avgx(tem7, 0, nx,ny,nz,                                          &
            1,nx-1,1,ny-1,1,nz-1,tem2)
  CALL avgy(tem8, 0, nx,ny,nz,                                          &
              1,nx-1,1,ny-1,1,nz-1,tem3)
  CALL avgz(tem9, 0, nx,ny,nz,                                          &
              1,nx-1,1,ny-1,1,nz-1,tem4)

!
!----------------------------------------------------------------------
!
!    Calculate the first guess radial velocity component using
!    the hole-filled and blended velocity data.
!
!----------------------------------------------------------------------
!

  DO k = 1, nz-1
    DO j = 1, ny-1
      DO i = 1, nx-1


        xs = 0.5*(x(i)+x(i+1)) - xmove
        ys = 0.5*(y(j)+y(j+1)) - ymove
        zs = 0.5*(z(k)+z(k+1)) - zshift

        rad         = SQRT(xs**2+ys**2+zs**2)
        tem5(i,j,k) = (tem2(i,j,k)*xs+tem3(i,j,k)*ys+                   &
                        tem4(i,j,k)*zs)/rad

      END DO
    END DO
  END DO
!
!----------------------------------------------------------------------
!
!  Save rhobar in tem9
!
!----------------------------------------------------------------------
!
  DO k=1,nz
    DO j=1,ny
      DO i=1,nx
        tem9(i,j,k) = rhostr(i,j,k)/j3(i,j,k)
      END DO
    END DO
  END DO
!
!----------------------------------------------------------------------
!
!    Check whether the variational adjusment is needed.
!
!----------------------------------------------------------------------
!
  IF ( insrt == 1.AND.ivar == 0) THEN

    WRITE(6,'(/5x,a)') 'Performing a direct insertion only'
!
!-----------------------------------------------------------------------
!
!  Calculate the model tangential velocity, vphi, and model polar
!  velocity, vtheta, at a scalar point.  The adjusted winds ua,va,
!  and wa, are constructed (at a scalar point) as follows:
!
!     tem6 = ua = Vr*x/rad + vphi*y/rcos - vtheta*x*z/rad*rcos
!     tem7 = va = Vr*y/rad - vphi*x/rcos - vtheta*y*z/rad*rcos
!     tem8 = wa = Vr*z/rad + vtheta*rcos/rad
!
!  Where,
!     rcos = sqrt(x**2+y**2)
!     rad  = sqrt(x**2+y**2+z**2)
!     tem1 = Vr (model (OSSE) or observed radial velocity)
!
!  For more details on this insertion technique, see the write-up
!  "Direct insertion of single-Doppler radial winds into a numerical
!  model" by Alan Shapiro.
!
!-----------------------------------------------------------------------
!

    DO k = 1, nz-1
      DO j = 1, ny-1
        DO i = 1, nx-1

          xs    = 0.5*(x(i)+x(i+1)) - xmove
          ys    = 0.5*(y(j)+y(j+1)) - ymove
          zs    = 0.5*(z(k)+z(k+1)) - zshift

          rad   = SQRT(xs**2+ys**2+zs**2)

          rcos  = SQRT(xs**2+ys**2)

          vphi  = (1./rcos)*(tem2(i,j,k)*ys - tem3(i,j,k)*xs)

          vtheta= (1./(rad*rcos))*(-tem2(i,j,k)*xs*zs                   &
                                   -tem3(i,j,k)*ys*zs                   &
                                   +tem4(i,j,k)*(rcos**2))

!-----------------------------------------------------
!        Comment out overwrite of obs radial component
!        01-03-99 SSW for AMS99
!
!        tem6(i,j,k) = tem1(i,j,k)*xs/rad + vphi*ys/rcos
!    :                  - vtheta*xs*zs/(rad*rcos)
!        tem7(i,j,k) = tem1(i,j,k)*ys/rad - vphi*xs/rcos
!    :                  - vtheta*ys*zs/(rad*rcos)
!        tem8(i,j,k) = tem1(i,j,k)*zs/rad + vtheta*rcos/rad
!
          tem6(i,j,k) = tem2(i,j,k)
          tem7(i,j,k) = tem3(i,j,k)
          tem8(i,j,k) = tem4(i,j,k)
!-----------------------------------------------------
!
        END DO
      END DO
    END DO

!
!-----------------------------------------------------------------------
!
!    Reconstruct the model velocity field on the Arakawa C-grid. u,v,
!    and w are stored in tem2,tem3, and tem4 respectively until the
!    boundary conditions are applied at the end of this routine.
!
!-----------------------------------------------------------------------
!

    DO k = 2, nz-1
      DO j = 2, ny-1
        DO i = 2, nx-1
          tem2(i,j,k) = .5*(tem6(i-1,j,k)+tem6(i,j,k))
          tem3(i,j,k) = .5*(tem7(i,j-1,k)+tem7(i,j,k))
          tem4(i,j,k) = .5*(tem8(i,j,k-1)+tem8(i,j,k))
        END DO
      END DO
    END DO

    GO TO 270                          ! Apply boundary conditions
!
!----------------------------------------------------------------------
!
!  Perform a variational adjustment
!
!----------------------------------------------------------------------
!
  ELSE

    PRINT *,'Performing a variational adjustment'

!
!----------------------------------------------------------------------
!
!  Redefine the following arrays such that:
!
!            tem1 = tem1*rhostr   tem4 = tem4*rhostr
!            tem2 = tem2*rhostr   tem5 = tem5*rhostr
!            tem3 = tem3*rhostr
!
!----------------------------------------------------------------------
!
    CALL aamult(tem1,tem9,nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,tem1)

    CALL avgsu(tem9,nx,ny,nz, 1,ny-1, 1,nz-1, tem6, tem10)
    CALL aamult(
    CALL avgsv(tem9,nx,ny,nz, 1,nx-1, 1,nz-1, tem6, tem10)
    CALL aamult(
    CALL avgsw(tem9,nx,ny,nz, 1,nx-1, 1,ny-1, tem6)
    CALL aamult(
    CALL aamult(tem5,tem9,nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,tem5)

    CALL avgsw(tem9,nx,ny,nz, 1,nx-1, 1,ny-1, tem12)  ! store rhobar at
                                                      ! w points.

    WRITE(6,*) 'after redefine'
!
!----------------------------------------------------------------------
!
!  Set the optimum over-relaxation coefficient, alpha.
!  WARNING: In order to avoid solution 'blow-up', alpha should be
!           near 1.
!
!----------------------------------------------------------------------
!
!     alpha=0.25    ! original value; under relaxation, slow convergence
!    alpha=1.00    ! under relaxation, slow convergence
    alpha=0.75    ! under relaxation, slow convergence
!     alpha=0.15    ! under relaxation, slow convergence
!     alpha=1.5    ! SOR; optimum -> 2 for large dimensions.

    DO k = 1,nz
      DO j = 1,ny
        DO i = 1,nx
          tem6(i,j,k)  = 0.0
          tem7(i,j,k)  = 0.0
          tem8(i,j,k)  = 0.0
          tem10(i,j,k) = 0.0
          tem11(i,j,k) = 0.0
        END DO
      END DO
    END DO

    dxinv     = 1./dx
    dyinv     = 1./dy
    dzinv     = 1./dz

    rdx2 = dxinv*dxinv        ! reciprocal of (dx)**2
    rdy2 = dyinv*dyinv        ! reciprocal of (dy)**2
    rdz2 = dzinv*dzinv        ! reciprocal of (dz)**2
!
!----------------------------------------------------------------------
!
!    Start the interation loop.
!
!----------------------------------------------------------------------
!
    WRITE(6,*) 'start iteration'

    DO it = 1,itmax
!
!----------------------------------------------------------------------
!
!  Compute the right hand side of the Poisson equation and store it
!  in tem6.  This will be handled externally once the code is
!  completed and tested.
!
!----------------------------------------------------------------------
!
      DO k = 2, nz-2
        DO j = 2, ny-2
          DO i = 2, nx-2
!
!
!       xs1       = x(i+1) - xmove
!       ys1       = y(j+1) - ymove
!       zs1       = z(k+1) - zshift
!
!       xs       = x(i) - xmove
!       ys       = y(j) - ymove
!       zs       = z(k) - zshift
!
!       term1 = dxinv*(tem2(i+1,j,k)-tem2(i,j,k))
!  :          + dyinv*(tem3(i,j+1,k)-tem3(i,j,k))
!  :          + dzinv*(tem4(i,j,k+1)-tem4(i,j,k))
!
!
!       term2 = dxinv*(xs1*(tem8(i+1,j,k)+tem8(i+1,j+1,k)
!  :          +             tem8(i+1,j,k+1)+tem8(i+1,j+1,k+1))
!  :          -         xs*(tem8(i,j,k)+tem8(i,j+1,k)
!  :          +             tem8(i,j,k+1)+tem8(i,j+1,k+1)))
!  :          + dyinv*(ys1*(tem8(i,j+1,k)+tem8(i+1,j+1,k)
!  :          +             tem8(i,j+1,k+1)+tem8(i+1,j+1,k+1))
!  :          -         ys*(tem8(i,j,k)+tem8(i+1,j,k)
!  :          +             tem8(i,j,k+1)+tem8(i+1,j,k+1)))
!  :          + dzinv*(zs1*(tem8(i,j,k+1)+tem8(i+1,j,k+1)
!  :          +             tem8(i,j+1,k+1)+tem8(i+1,j+1,k+1))
!  :          -         zs*(tem8(i,j,k)+tem8(i+1,j,k)
!  :          +             tem8(i,j+1,k)+tem8(i+1,j+1,k)))
!
!       term2       = term2/scoef
!
!       tem6(i,j,k) = -2.0*term1 + 0.25*term2
!
! **** new formulation
!
            xs       = 0.5 * ( x(i) + x(i+1) ) - xmove
            ys       = 0.5 * ( y(j) + y(j+1) ) - ymove
            zs       = 0.5 * ( z(k) + z(k+1) ) - zshift
            rrr      = SQRT(xs*xs + ys*ys + zs*zs)

            term1 = dxinv*(tem2(i+1,j,k)-tem2(i,j,k))                   &
                  + dyinv*(tem3(i,j+1,k)-tem3(i,j,k))                   &
                  + dzinv*(tem4(i,j,k+1)-tem4(i,j,k))

            term2 = xs*dxinv*( (tem8(i+1,j,k)+tem8(i+1,j+1,k)           &
                  +             tem8(i+1,j,k+1)+tem8(i+1,j+1,k+1))      &
                  -            (tem8(i,j,k)+tem8(i,j+1,k)               &
                  +             tem8(i,j,k+1)+tem8(i,j+1,k+1)) )        &
                           + ys*dyinv*( (tem8(i,j+1,k)+tem8(i+1,j+1,k)  &
                           +             tem8(i,j+1,k+1)+tem8(i+1,j+1,k+1)) &
                           -            (tem8(i,j,k)+tem8(i+1,j,k)      &
                  +             tem8(i,j,k+1)+tem8(i+1,j,k+1)) )        &
                           + zs*dzinv*( (tem8(i,j,k+1)+tem8(i+1,j,k+1)  &
                           +             tem8(i,j+1,k+1)+tem8(i+1,j+1,k+1)) &
                           -            (tem8(i,j,k)+tem8(i+1,j,k)      &
                           +             tem8(i,j+1,k)+tem8(i+1,j+1,k)) )

            term2 = 0.25 * term2 + 2.0 * 0.125 *                        &
                     ( tem8(i,  j,  k  )+tem8(i+1,j,  k  ) +            &
                       tem8(i,  j+1,k  )+tem8(i+1,j+1,k  ) +            &
                       tem8(i,  j,  k+1)+tem8(i+1,j,  k+1) +            &
                       tem8(i,  j+1,k+1)+tem8(i+1,j+1,k+1) )

            term2 = term2/rrr

            tem6(i,j,k) = - 2.0 * term1 + term2

          END DO
        END DO
      END DO
!
!----------------------------------------------------------------------
!
!  Solve the elliptic equation via the successive overrelaxation
!  method.  Check for convergence. The following coefficients
!  are defined as:
!
!  b1 = rdx2
!  b2 = rdy2
!  b3 = rdz2
!
!----------------------------------------------------------------------
!
      b1          = rdx2
      b2          = rdy2
      b3          = rdz2
      c1          = -2.*(b1 + b2 + b3)
      c2          = b1
      c3          = b1
      c4          = b2
      c5          = b2
      c6          = b3
      c7          = b3

      test   = 0.0
      DO k = 2, nz-2
        DO j = 2, ny-2
          DO i = 2, nx-2

            tm1         = c2*tem7(i+1,j,k)
            tm2         = c3*tem7(i-1,j,k)
            tm3         = c4*tem7(i,j+1,k)
            tm4         = c5*tem7(i,j-1,k)
            tm5         = c6*tem7(i,j,k+1)
            tm6         = c7*tem7(i,j,k-1)

            resid  = c1*tem7(i,j,k) + tm1 + tm2 + tm3 + tm4             &
                                    + tm5 + tm6 - tem6(i,j,k)

            newit  = tem7(i,j,k) - alpha*resid/c1   ! SOR method

!       diff   = abs(newit      - tem7(i,j,k))
!       diff = abs(newit - tem7(i,j,k)) / tem9(i,j,k)

            diff = ABS(resid/c1) / tem9(i,j,k)   ! alpha excluded;
                                                 ! normakized residual

            IF(diff > test) THEN
              test = diff
              imx = i
              jmx = j
              kmx = k
            END IF

!
!       tem7(i,j,k) = newit

            tem11(i,j,k) = newit

            tem10(i,j,k) = resid          !outputs for test

          END DO
        END DO
      END DO

!      write(6,*) 'after relaxation'

!ccxxx
!ccxxx Test simutenous relaxation
!ccxxx
      DO k = 2, nz-2
        DO j = 2, ny-2
          DO i = 2, nx-2

            tem7(i,j,k) = tem11(i,j,k)

          END DO
        END DO
      END DO
!
!----------------------------------------------------------------------
!
!  Set the boundary conditions on Xl, Xr, Yl, & Yr
!
!----------------------------------------------------------------------
!
      DO k = 2, nz-2          !XL,XR B.C.
        DO j = 2, ny-2
          tem7(1,j,k)    = - tem7(2,j,k)
          tem7(nx-1,j,k) = - tem7(nx-2,j,k)
        END DO
      END DO

      DO k = 2,nz-2          !YL,YR B.C
        DO i = 2,nx-2
          tem7(i,1,k)    = - tem7(i,2,k)
          tem7(i,ny-1,k) = - tem7(i,ny-2,k)
        END DO
      END DO
!
!----------------------------------------------------------------------
!
!  Set the top and bottom boundary conditions.
!
!----------------------------------------------------------------------
!
      z2   = z(2)    - zshift
      znz1 = z(nz-1) - zshift
      DO i=2,nx-2
        DO j=2,ny-2

!       xs = x(i) - xmove
!       ys = y(j) - ymove
!
!       radh = sqrt(xs*xs + ys*ys)
!
!       tem7(i,j,1)    = tem7(i,j,2)          !'lamda2' grad B.C.
!  : -                   0.25*dz*z2*(tem8(i,j,2)+tem8(i,j+1,2)
!  : +                   tem8(i+1,j,2)+tem8(i+1,j+1,2))/scoef
!
!       tem7(i,j,nz-1) = tem7(i,j,nz-2)
!  : +                   0.25*dz*znz1*(tem8(i,j,nz-1)+tem8(i,j+1,nz-1)
!  : +                   tem8(i+1,j,nz-1)+tem8(i+1,j+1,nz-1))/scoef
!
          xs = 0.5 * ( x(i) + x(i+1) ) - xmove
          ys = 0.5 * ( y(j) + y(j+1) ) - ymove

          rrr = SQRT(xs*xs + ys*ys + z2*z2)

          tem7(i,j,1)    = tem7(i,j,2)          & !'lamda2' grad B.C.
          -                   0.25*dz*z2*(tem8(i,j,2)+tem8(i,j+1,2)     &
              +                   tem8(i+1,j,2)+tem8(i+1,j+1,2)) / rrr  &
              +                   2.0 * dz * tem4(i,j,2)

          rrr = SQRT(xs*xs + ys*ys + znz1*znz1)

          tem7(i,j,nz-1) = tem7(i,j,nz-2)                               &
              +                   0.25*dz*znz1*(tem8(i,j,nz-1)+tem8(i,j+1,nz-1) &
              +                   tem8(i+1,j,nz-1)+tem8(i+1,j+1,nz-1)) / rrr &
              -                   2.0 * dz * tem4(i,j,nz-1)


        END DO
      END DO

!ccxxx
!ccxxx  Patch: you might like to change it.
!ccxxx

      ir0 = xmove/dx + 1
      jr0 = ymove/dy + 1
      ipatch =5
      ips=ir0-ipatch
      ipe=ir0+ipatch
      jps=jr0-ipatch
      jpe=jr0+ipatch

      IF((ir0-ipatch) < 2) ips=2
      IF((jr0-ipatch) < 2) jps=2
      IF((ir0+ipatch) > (nx-2)) ipe=nx-2
      IF((jr0+ipatch) > (ny-2)) jpe=ny-2

      DO in = ips,ipe
        DO jn = jps,jpe


          acof = (jn-jps)/(jpe-jps)
          tem7(in,jn,nz-1) = (tem7(in,jpe,nz-1)-tem7(in,jps,nz-1))*acof &
                            +  tem7(in,jps,nz-1)

          acof = (in-ips)/(ipe-ips)
          tem7(in,jn,nz-1) = 0.5*(tem7(in,jn,nz-1) + ((tem7(ipe,jn,nz-1) &
                         - tem7(ips,jn,nz-1))*acof) + tem7(ips,jn,nz-1))

!    acof = (in-ips)/(ipe-ips)
!    bcof = (jn-jps)/(jpe-jps)
!    tem7(in,jn,nz-1) = (1.0-acof)*(1.0-bcof)*tem7(ips,jps,nz-1)
!  :                  + (1.0-acof)*     bcof *tem7(ips,jpe,nz-1)
!  :                  +      acof *     bcof *tem7(ipe,jpe,nz-1)
!  :                  +      acof *(1.0-bcof)*tem7(ipe,jps,nz-1)
!
!    tem7(in,jn,nz-1) = 0.5 * (tem7(in,jn,nz-1) + tem7(in,jn,nz-2))
!
!    tem7(in,jn,nz-1) = tem7(in,jn,nz-2)     ! zeor gradient

        END DO
      END DO


!   ir0 = xmove/dx + 1
!   jr0 = ymove/dy + 1
!   ipatch = 10             ! within a radius of 20 km
!   ips=ir0-ipatch
!   ipe=ir0+ipatch
!   jps=jr0-ipatch
!   jpe=jr0+ipatch

!   if((ir0-ipatch).lt.2) ips=2
!   if((jr0-ipatch).lt.2) jps=2
!   if((ir0+ipatch).gt.(nx-2)) ipe=nx-2
!   if((jr0+ipatch).gt.(ny-2)) jpe=ny-2
!
!   DO in = ips,ipe
!   DO jn = jps,jpe
!     tem7(in,jn,1) = 0.5 * ( tem7(in,jn,1) +
!    :          0.25*(tem7(in-1,jn,1)+tem7(in+1,jn,1)+
!    :                tem7(in,jn-1,1)+tem7(in,jn+1,1)) )
!     tem7(in,jn,nz-1) = 0.5 * ( tem7(in,jn,nz-1) +
!    :          0.25*(tem7(in-1,jn,nz-1)+tem7(in+1,jn,nz-1)+
!    :                tem7(in,jn-1,nz-1)+tem7(in,jn+1,nz-1)) )
!   END DO
!   END DO

!   DO  k = 2, nz-2
!   DO in = ips,ipe
!   DO jn = jps,jpe
!     tem7(in,jn,k) = 0.5 * ( tem7(in,jn,k) +
!    :               (tem7(in-1,jn,k)+tem7(in+1,jn,k)+
!    :                tem7(in,jn-1,k)+tem7(in,jn+1,k)+
!    :                tem7(in,jn,k-1)+tem7(in,jn,k+1))/6.0 )
!   END DO
!   END DO
!   END DO

      DO i=2,nx-2
        tem7(i,1,1)       = - tem7(i,2,1)
        tem7(i,1,nz-1)    = - tem7(i,2,nz-1)
        tem7(i,ny-1,1)    = - tem7(i,ny-2,1)
        tem7(i,ny-1,nz-1) = - tem7(i,ny-2,nz-1)
      END DO

      DO j=2,ny-2
        tem7(1,j,1)       = - tem7(2,j,1)
        tem7(1,j,nz-1)    = - tem7(2,j,nz-1)
        tem7(nx-1,j,1)    = - tem7(nx-2,j,1)
        tem7(nx-1,j,nz-1) = - tem7(nx-2,j,nz-1)
      END DO
!
!----------------------------------------------------------------------
!
!    Set the corner points.
!
!----------------------------------------------------------------------
!
      DO k = 1,nz-1
        tem7(1,1,k)       = - tem7(2,2,k)
        tem7(1,ny-1,k)    = - tem7(2,ny-2,k)
        tem7(nx-1,1,k)    = - tem7(nx-2,2,k)
        tem7(nx-1,ny-1,k) = - tem7(nx-2,ny-2,k)
      END DO
!
!----------------------------------------------------------------------
!
!   Determine the first Lagrange Multiplier, Lamda1. Lamda1 is
!   assumed valid at integer points(i,j,k). The result is
!   stored in tem8.
!
!----------------------------------------------------------------------
!
      test1=0.0
      DO k = 2, nz-1
        DO j = 2, ny-1
          DO i = 2, nx-1

            xs       = x(i) - xmove
            ys       = y(j) - ymove
            zs       = z(k) - zshift

            rad        = SQRT(xs**2+ys**2+zs**2)
            radinv     = 1./(rad**2)


            IF( rad < 1.0 ) THEN
              rad = 1.0
              radinv     = 1./(rad**2)
            END IF

            term1 = dxinv*xs*((tem7(i,j-1,k-1)-tem7(i-1,j-1,k-1))       &
                +                    (tem7(i,j-1,k)-tem7(i-1,j-1,k))    &
                +                    (tem7(i,j,k-1)-tem7(i-1,j,k-1))    &
                +                    (tem7(i,j,k)-tem7(i-1,j,k)))

            term2 = dyinv*ys*((tem7(i-1,j,k-1)-tem7(i-1,j-1,k-1))       &
                +                    (tem7(i,j,k-1)-tem7(i,j-1,k-1))    &
                +                    (tem7(i-1,j,k)-tem7(i-1,j-1,k))    &
                +                    (tem7(i,j,k)-tem7(i,j-1,k)))

            term3 = dzinv*zs*((tem7(i,j-1,k)-tem7(i,j-1,k-1))           &
                +                    (tem7(i-1,j-1,k)-tem7(i-1,j-1,k-1)) &
                +                    (tem7(i-1,j,k)-tem7(i-1,j,k-1))    &
                +                    (tem7(i,j,k)-tem7(i,j,k-1)))

            term4 = 0.125*((tem5(i,j,k)+tem5(i-1,j,k)                   &
                +                  tem5(i-1,j-1,k)+tem5(i,j-1,k)        &
                +                  tem5(i,j,k-1)+tem5(i-1,j,k-1)        &
                +                  tem5(i-1,j-1,k-1)+tem5(i,j-1,k-1))   &
                -                 (tem1(i,j,k)+tem1(i-1,j,k)            &
                +                  tem1(i-1,j-1,k)+tem1(i,j-1,k)        &
                +                  tem1(i,j,k-1)+tem1(i-1,j,k-1)        &
                +                  tem1(i-1,j-1,k-1)+tem1(i,j-1,k-1)))

!        term5  = 2.0*term4/rad + 0.25*radinv*(term1+term2+term3)

            term5  = 2.0*term4 + 0.25*(term1+term2+term3)/rad

            diff = ABS(term5-tem8(i,j,k)) / tem12(i,j,k)
            IF(diff > test1) THEN
              test1 = diff
              imx1 = i
              jmx1 = j
              kmx1 = k
            END IF

            tem8(i,j,k) = scoef*term5

          END DO
        END DO
      END DO

      DO k = 2, nz-1
        DO j = 2, ny-1
          DO i = 2, nx-1
            xs = x(i) - xmove
            ys = y(j) - ymove
            zs = z(k) - zshift
            rad = SQRT(xs**2+ys**2+zs**2)

            IF( rad < 1.0 ) THEN
              tem8(i,j,k) = 0.0
              radinv = 0.0

              IF(i > 2) THEN
                tem8(i,j,k) = tem8(i,j,k) + tem8(i-1,j,k)
                radinv = radinv + 1.
              END IF

              IF(i < nx-1) THEN
                tem8(i,j,k) = tem8(i,j,k) + tem8(i+1,j,k)
                radinv = radinv + 1.
              END IF

              IF(j > 2) THEN
                tem8(i,j,k) = tem8(i,j,k) +  tem8(i,j-1,k)
                radinv = radinv + 1.
              END IF

              IF(j < ny-1) THEN
                tem8(i,j,k) = tem8(i,j,k) +  tem8(i,j+1,k)
                radinv = radinv + 1.
              END IF

              IF(k > 2) THEN
                tem8(i,j,k) = tem8(i,j,k) +  tem8(i,j,k-1)
                radinv = radinv + 1.
              END IF

              IF(k < nz-1) THEN
                tem8(i,j,k) = tem8(i,j,k) +  tem8(i,j,k+1)
                radinv = radinv + 1.
              END IF

              tem8(i,j,k) = tem8(i,j,k) / radinv

              WRITE(6,*) 'Do patching for lamda1 at i,j,k=',i,j,k
            END IF

          END DO
        END DO
      END DO

!
!----------------------------------------------------------------------
!
!  Test for solution convergence.
!
!----------------------------------------------------------------------
!
!     tol = sqrt(dx**2+dy**2+dz**2)*veltol
      tol = MIN( MIN(dx, dy), dz) * veltol

!     ix = mod(it,1)
      ix = MOD(it,50)
      IF(it == 50) THEN
        WRITE(6,*) ' Max diff at iteration: imx,jmx,kmx,test,test1'
      END IF

      IF(ix == 0) THEN
        WRITE(6,'(i9,2x,2i4,i3,e16.6,2x,2i4,i3,e16.6)')                 &
            it,imx,jmx,kmx,test, imx1,jmx1,kmx1,test1
      END IF

!     IF(test.le.tol.and.it.gt.2) THEN
      IF(test <= tol .AND. test1 <= veltol .AND. it > 2) THEN

        WRITE(6,'(/5x,a,i6,/5x,a,2f12.6)')                              &
            'The SOR technique converged. The # of iterations was ',    &
            it,' The tolerance level was set at: ',veltol,tol
        WRITE(6,'(i9,2x,2i4,i3,e16.6,2x,2i4,i3,e16.6)')                 &
            it,imx,jmx,kmx,test, imx1,jmx1,kmx1,test1

        GO TO 258

      END IF

    END DO

    STOP


    258     CONTINUE

    WRITE(6,*) 'patch: xmove,ymove,ir0,jr0,ipatch=',                    &
                       xmove,ymove,ir0,jr0,ipatch
    WRITE(6,*) 'patch: ips,ipe,jps,jpe=',ips,ipe,jps,jpe

    2201    CONTINUE

!
!----------------------------------------------------------------------
!
!   Determine the adjusted velocity, ua, store results in tem2:
!
!           2*rhobar*(Ua-Um) + lamda1*x - dlamda2/dx = 0
!
!   NOTE:  rhobar (tem9) is assumed horizontally homogenous.
!
!----------------------------------------------------------------------
!
    DO k = 1,nz-1
      DO j = 1,ny-1
        DO i = 1,nx-1
          tem9(i,j,k)= rhostr(i,j,k)/j3(i,j,k)
        END DO
      END DO
    END DO

    CALL avgsu(tem9,nx,ny,nz, 1,ny-1, 1,nz-1, tem6, tem10)

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

!       xs       = x(i) - xmove
!
!       tem2(i,j,k)  =  u(i,j,k,tim) + (.5/tem6(i,j,k))
!  :                    *(dxinv*(tem7(i,j,k)-tem7(i-1,j,k))
!  :                 -  0.25*xs*(tem8(i,j,k)+tem8(i,j+1,k)
!  :                 +  tem8(i,j,k+1)+tem8(i,j+1,k+1)))
!
          xs  =         x(i)            - xmove
          ys  = 0.5 * ( y(j) + y(j+1) ) - ymove
          zs  = 0.5 * ( z(k) + z(k+1) ) - zshift
          rrr = SQRT(xs*xs + ys*ys + zs*zs)

          tem2(i,j,k)  =  u(i,j,k,tim) + (.5/tem6(i,j,k))               &
                       * ( dxinv*(tem7(i,j,k)-tem7(i-1,j,k))            &
                         - 0.25*(xs/rrr)*(tem8(i,j,k)+tem8(i,j+1,k)     &
                           +tem8(i,j,k+1)+tem8(i,j+1,k+1)) )

        END DO
      END DO
    END DO

!
!----------------------------------------------------------------------
!
!   Determine the adjusted velocity, va, store results in tem3:
!
!           2*rhobar*(Va-Vm) + lamda1*y - dlamda2/dy = 0
!
!   NOTE:  rhobar (tem9) is assumed horizontally homogenous.
!
!----------------------------------------------------------------------
!
    CALL avgsv(tem9,nx,ny,nz, 1,nx-1, 1,nz-1, tem6, tem10)

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

!        ys       = y(j) - ymove
!
!        tem3(i,j,k)  =  v(i,j,k,tim) + (.5/tem6(i,j,k))
!  :                     *(dyinv*(tem7(i,j,k)-tem7(i,j-1,k))
!  :                  -  0.25*ys*(tem8(i,j,k)+tem8(i+1,j,k)
!  :                  +  tem8(i,j,k+1)+tem8(i+1,j,k+1)))

          xs  = 0.5 * ( x(i) + x(i+1) ) - xmove
          ys  =         y(j)            - ymove
          zs  = 0.5 * ( z(k) + z(k+1) ) - zshift
          rrr = SQRT(xs*xs + ys*ys + zs*zs)

          tem3(i,j,k)  =  v(i,j,k,tim) + (.5/tem6(i,j,k))               &
                       * ( dyinv*(tem7(i,j,k)-tem7(i,j-1,k))            &
                         - 0.25*(ys/rrr)*(tem8(i,j,k)+tem8(i+1,j,k)     &
                           +tem8(i,j,k+1)+tem8(i+1,j,k+1)) )

        END DO
      END DO
    END DO

!
!----------------------------------------------------------------------
!
!   Determine the adjusted velocity, wa, store results in tem4:
!
!           2*rhobar*(Wa-Wm) + lamda1*z - dlamda2/dz = 0
!
!----------------------------------------------------------------------
!
    CALL avgsw(tem9,nx,ny,nz, 1,nx-1, 1,ny-1, tem6)

    WRITE(6,*) 'Maximum vertical velocity after adjustment'
    DO k = 2, nz-1
      maxw=0
      DO j = 2, ny-2
        DO i = 2, nx-2

!       zs       = z(k) - zshift
!
!        tem4(i,j,k)  =  w(i,j,k,tim) + (0.5/tem6(i,j,k))
!  :                     *(dzinv*(tem7(i,j,k)-tem7(i,j,k-1))
!  :                  -  0.25*zs*(tem8(i,j,k)+tem8(i+1,j,k)
!  :                  +  tem8(i,j+1,k)+tem8(i+1,j+1,k)))

          xs  = 0.5 * ( x(i) + x(i+1) ) - xmove
          ys  = 0.5 * ( y(j) + y(j+1) ) - ymove
          zs  =         z(k)            - zshift
          rrr = SQRT(xs*xs + ys*ys + zs*zs)

          tem4(i,j,k)  =  w(i,j,k,tim) + (0.5/tem6(i,j,k))              &
                       * ( dzinv*(tem7(i,j,k)-tem7(i,j,k-1))            &
                         - 0.25*(zs/rrr)*(tem8(i,j,k)+tem8(i+1,j,k)     &
                           +tem8(i,j+1,k)+tem8(i+1,j+1,k)) )

          IF (ABS(tem4(i,j,k)) > maxw) THEN
            maxw = ABS(tem4(i,j,k))
            imx=i
            jmx=j
          END IF
        END DO
      END DO

      WRITE(6,*) '  imx,jmx,k,maxw: ',imx,jmx,k,maxw

    END DO

  END IF

!
!-----------------------------------------------------------------------
!
!  Exit adjustment/or direct insertion
!
!-----------------------------------------------------------------------
!

  270   CONTINUE

!
!-----------------------------------------------------------------------
!
!  Compare the model radial velocity and observed radial velocity
!
!-----------------------------------------------------------------------
!
  CALL avgx(tem2, 0, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,tem6)
  CALL avgy(tem3, 0, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,tem7)
  CALL avgz(tem4, 0, nx,ny,nz,1,nx-1,1,ny-1,1,nz-1,tem8)

  WRITE(6,*) 'Check vr & print out if abs(vradj-vro) > 1.5 m/s'

  DO k=2,nz-2
    DO j=2,ny-2
      DO i=2,nx-2
        xs    = 0.5*(x(i)+x(i+1)) - xmove
        ys    = 0.5*(y(j)+y(j+1)) - ymove
        zs    = 0.5*(z(k)+z(k+1)) - zshift

        rad   = SQRT(xs**2+ys**2+zs**2)
        vradj = (tem6(i,j,k)*xs + tem7(i,j,k)*ys                        &
                     + tem8(i,j,k)*zs)/rad
        vro   = tem1(i,j,k)/tem9(i,j,k)
        vrm   = tem5(i,j,k)/tem9(i,j,k)

!      IF(abs(vradj-vro).gt.1.5) THEN
!        write(6,*) '  i,j,k,vradj,vro,vrm: ',i,j,k,vradj,vro,vrm
!      ENDIF

      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  NOTE: This next section is designed to determine whether or not
!        mass conservation has been satisfied.  It is not an essential
!        part of the assimilation package.
!
!       Stagger rhostar(stored in tem9).
!
!-----------------------------------------------------------------------
!
  CALL avgsu(tem9,nx,ny,nz, 1,ny-1, 1,nz-1, tem6, tem10)
  CALL avgsv(tem9,nx,ny,nz, 1,nx-1, 1,nz-1, tem7, tem10)
  CALL avgsw(tem9,nx,ny,nz, 1,nx-1, 1,ny-1, tem8)
!
!-----------------------------------------------------------------------
!
!  Calculate del dot rho V using the adjusted velocities.
!
!-----------------------------------------------------------------------
!

  WRITE(6,*) 'background divergence'

  DO k=2,nz-2
    maxdiv = 0.0
    DO j=4,ny-4
      DO i=4,nx-4
        tem1(i,j,k)=                                                    &
            (u(i+1,j,k,tim)*tem6(i+1,j,k)-u(i,j,k,tim)*tem6(i,j,k))*dxinv &
            +(v(i,j+1,k,tim)*tem7(i,j+1,k)-v(i,j,k,tim)*tem7(i,j,k))*dyinv &
            +(w(i,j,k+1,tim)*tem8(i,j,k+1)-w(i,j,k,tim)*tem8(i,j,k))*dzinv
        IF (ABS(tem1(i,j,k)) > maxdiv) THEN
          maxdiv = ABS(tem1(i,j,k))
          imx = i
          jmx = j
        END IF
      END DO
    END DO
    maxdiv=tem1(imx,jmx,k)
    PRINT *,'  maxdiv,imx,jmx,k = ',maxdiv,imx,jmx,k
  END DO

  WRITE(6,*) 'adjusted divergence'

  DO k=2,nz-2
    maxdiv = 0.0
    DO j=4,ny-4
      DO i=4,nx-4
        tem1(i,j,k)=                                                    &
            (tem2(i+1,j,k)*tem6(i+1,j,k)-tem2(i,j,k)*tem6(i,j,k))*dxinv &
            +(tem3(i,j+1,k)*tem7(i,j+1,k)-tem3(i,j,k)*tem7(i,j,k))*dyinv &
            +(tem4(i,j,k+1)*tem8(i,j,k+1)-tem4(i,j,k)*tem8(i,j,k))*dzinv
        IF (ABS(tem1(i,j,k)) > maxdiv) THEN
          maxdiv = ABS(tem1(i,j,k))
          imx = i
          jmx = j
        END IF
      END DO
    END DO
    maxdiv=tem1(imx,jmx,k)
    PRINT *,'  maxdiv,imx,jmx,k = ',maxdiv,imx,jmx,k
  END DO

  RETURN
END SUBROUTINE assimvel