!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE ADJUVW                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE adjuvw( nx,ny,nz, u,v,w,wcont,ptprt,ptbar,                   & 3,8
           zp,j1,j2,j3,aj3z,mapfct,rhostr,wcloud,                       &
           wndadj,obropt,obrzero,cldwopt,                               &
           rhostru,rhostrv,rhostrw,rhs,divh,tem1,tem2,phi)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Diagnose w and wcont from u and v fields. Adjust u and v to ensure
!  anelastic mass conservation at each grid point.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  3/27/1995
!
!  MODIFICATION HISTORY:
!
!  4/19/1995 (M. Xue)
!  Redefined phi as a 2-D array, and added it into the argument list.
!
!  9/25/1996 (K. Brewster)
!  Added optional processing to set wcont=0 at bottom of
!  Rayleigh damping layer.
!
!  2/16/1998 (K. Brewster)
!  Added processing of wcloud, an estimate of w from the cloud
!  analysis.
!
!  7/27/1998 (M. Xue)
!  Map factor properly included. WCTOW and WCONTRA are now called
!  to convert between w and wcont. WCONTRA required that crdtrns
!  be defined, which is added to INITADAS.
!
!  2001-06-01 (Gene Bassett)
!  Fixed use of obropt > 9.  If subroutine was called more than once,
!  subsequent calls would use obropt-10 instead of obropt.
!
!-----------------------------------------------------------------------
!
!  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 vertical
!
!    u        x component of velocity at a given time level (m/s)
!    v        y component of velocity at a given time level (m/s)
!    w        z component of velocity at a given time level (m/s)
!    ptprt    perturbation potential temperature (K)
!    ptbar    base state potential temperature (K)
!    zp       model physical height
!    j1       Coordinate transformation Jacobian -d(zp)/dx
!    j2       Coordinate transformation Jacobian -d(zp)/dy
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!    mapfct   map factor at scalar, u and v points.
!    rhostr   Base state density rhobar times j3 (kg/m**3)
!    wcloud   W from cloud analysis, if available
!
!    wndadj   Wind adjustment option.
!    obropt   O'Brien adjustment option.  Determines
!             distribution of mean divergence error used to
!             enforce upper boundary condition of w=0
!             = 1  Linearly in computational z. (default)
!             = 2  Linearly in physical z.
!             = 3  Linearly in potential temperature.
!    cldwopt  Option to replace w with w from cloud analysis
!
!  OUTPUT:
!
!    u        x component of velocity at a given time level (m/s)
!    v        y component of velocity at a given time level (m/s)
!    w        Vertical component of velocity in Cartesian
!             coordinates at a given time level (m/s)
!    wcont    Contravariant vertical velocity (m/s)
!
!  WORK ARRAYS:
!
!    rhostru  Temporary work array.
!    rhostrv  Temporary work array.
!    rhostrw  Temporary work array.
!    rhs      Temporary work array.
!    divh     Temporary work array storing horizontal mass divergence.
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    phi      Temporary work array.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE

  INTEGER :: nx,ny,nz
!
  INTEGER :: wndadj
  INTEGER :: obropt
  REAL :: obrzero
  INTEGER :: cldwopt

  REAL :: u     (nx,ny,nz)     ! Total u-velocity (m/s)
  REAL :: v     (nx,ny,nz)     ! Total v-velocity (m/s)
  REAL :: w     (nx,ny,nz)     ! Total w-velocity (m/s)
  REAL :: ptprt (nx,ny,nz)     ! Perturbation potential temperature (K)
  REAL :: wcont (nx,ny,nz)     ! Contravariant vertical velocity (m/s)
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)

  REAL :: mapfct(nx,ny,8)
  REAL :: rhostr(nx,ny,nz)     ! Base state density rhobar times j3.
  REAL :: wcloud(nx,ny,nz)

  REAL :: zp(nx,ny,nz)

  REAL :: j1    (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! - d( zp )/d( x ).
  REAL :: j2    (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! - d( zp )/d( y ).
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! d( zp )/d( z ).
  REAL :: aj3z  (nx,ny,nz)     ! Coordinate transformation Jacobian defined
                               ! as d( zp )/d( z ) AVERAGED IN THE Z-DIR.

  REAL :: rhostru(nx,ny,nz)    ! Work array
  REAL :: rhostrv(nx,ny,nz)    ! Work array
  REAL :: rhostrw(nx,ny,nz)    ! Work array
  REAL :: rhs    (nx,ny,nz)    ! Work array
  REAL :: divh   (nx,ny,nz)    ! Work array
  REAL :: tem1   (nx,ny,nz)    ! Work array
  REAL :: tem2   (nx,ny,nz)    ! Work array

  REAL :: phi  (0:nx,0:ny)     ! Work array
!
!-----------------------------------------------------------------------
!
!  Control parameters
!
!-----------------------------------------------------------------------
!
  LOGICAL :: chkmas
  PARAMETER (chkmas=.true.)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,ktop
  INTEGER :: iterations,imax,imin,jmax,jmin,kmax,kmin
  REAL :: dxx,dyy,absphi,delphi,rij,dphi
  REAL :: divmax, divmin, wcmax, wcmin, dz2, dth2, zk, zk2, depth2
  REAL :: pi, alpha, rdxx, rdyy, tem
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'        ! Grid parameters
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  CALL rhouvw(nx,ny,nz,rhostr,rhostru,rhostrv,rhostrw)

  IF(wndadj == 1) THEN

    DO k=1,nz
      DO j=1,ny-1
        DO i=1,nx-1
          wcont(i,j,k)=0.
        END DO
      END DO
    END DO

    CALL wctow(nx,ny,nz,u,v,wcont,mapfct,                               &
               j1,j2,j3,aj3z,rhostr,rhostru,rhostrv,rhostrw,1,1,        &
               w,                                                       &
               tem1,tem2)

  ELSE IF (wndadj == 2 .OR. wndadj == 3 ) THEN

    w = 0.0
    wcont = 0.0
 
!-----------------------------------------------------------------------
!
!  Storing m**2 difx(ustr/m)+dify(vstr/m) in divh.
!
!-----------------------------------------------------------------------
!
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          divh(i,j,k)=((u(i+1,j,k)*rhostru(i+1,j,k)*mapfct(i+1,j,5)     &
                       -u(i  ,j,k)*rhostru(i  ,j,k)*mapfct(i  ,j,5))    &
                       *dxinv                                           &
                      +(v(i,j+1,k)*rhostrv(i,j+1,k)*mapfct(i,j+1,6)     &
                       -v(i,j  ,k)*rhostrv(i,j  ,k)*mapfct(i,j  ,6))    &
                       *dyinv) *mapfct(i,j,1)*mapfct(i,j,1)
          tem1(i,j,k)=zp(i,j,k)
        END DO
      END DO
    END DO

    DO k=3,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          tem2(i,j,k)=0.5*(ptprt(i,j,k)  +ptbar(i,j,k) +                &
                           ptprt(i,j,k-1)+ptbar(i,j,k-1))
        END DO
      END DO
    END DO

    DO j=1,ny-1
      DO i=1,nx-1
        wcont(i,j,2)=0.
      END DO
    END DO

    DO k=3,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          wcont(i,j,k)= wcont(i,j,k-1) - dz*divh(i,j,k-1)
        END DO
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  wcont is now carrying wcont * rhostrw
!
!-----------------------------------------------------------------------
!

    WRITE(6,'(a,i4)') 'obropt = ',obropt

    IF(obropt > 10) THEN
      ktop=0
      DO j=2,ny-1
        DO i=2,nx-1
          DO k=nz-2,3,-1
            IF(zp(i,j,k) < obrzero) EXIT
          END DO
!          271       CONTINUE
          ktop=ktop+k
        END DO
      END DO

      ktop=nint(FLOAT(ktop)/FLOAT((nx-2)*(ny-2)))
      ktop=MIN((ktop+1),nz-1)

      WRITE(6,'(a,i5,a,f9.0,a)')                                        &
          '   Found k-index of bottom of Rayleigh damping zone:',       &
          ktop,' obrzero=',obrzero,' m.'

      DO j=1,ny-1
        DO i=1,nx-1
          DO k=ktop+1,nz
            wcont(i,j,k)=0.
          END DO
        END DO
      END DO


    ELSE

      ktop=nz-1

    END IF

    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,1)=1./((tem1(i,j,ktop)-tem1(i,j,2))*                   &
                        (tem1(i,j,ktop)-tem1(i,j,2)))
        tem2(i,j,2)=1.5*tem2(i,j,3)-0.5*tem2(i,j,4)
        tem2(i,j,1)=1./((tem2(i,j,ktop)-tem2(i,j,2))*                   &
                      (tem2(i,j,ktop)-tem2(i,j,2)))
      END DO
    END DO

!
    WRITE(6,'(a)')                                                      &
        ' Before w adjust, rho*wcont(top):'

    CALL a3dmax(wcont,1,nx,1,nx-1,                                      &
              1,ny,1,ny-1,1,nz,ktop,ktop,                               &
              wcmax,wcmin, imax,jmax,kmax, imin,jmin,kmin)

    WRITE(6,'(2(1x,a,g12.4,3(a,i3)))')                                  &
        'wcmin =',wcmin,' at i=',imin,', j=',jmin,', k=',kmin,          &
        'wcmax =',wcmax,' at i=',imax,', j=',jmax,', k=',kmax

    IF( obropt == 1 .or. obropt == 11) THEN
!
!-----------------------------------------------------------------------
!
!  For obropt = 1,
!  we apply correction to wcont * rhostrw by assuming that
!  difz ( delta(wcont * rhostr) ) = A k, where A is
!  a constant. Therefore  delta( wcont*rhostr ) = A k**2/2.
!
!-----------------------------------------------------------------------
!
      WRITE(6,'(a)')                                                    &
          'Horizontal divergence correction distributed linearly in k.'

      depth2 = 1./((ktop-2)*(ktop-2)*dz*dz)

      DO k=3,ktop
        zk2 = ((k-2)*dz)*((k-2)*dz)
        DO j=1,ny-1
          DO i=1,nx-1
            wcont(i,j,k)=wcont(i,j,k)-                                  &
                         wcont(i,j,ktop)*(zk2*depth2)
          END DO
        END DO
      END DO

    ELSE IF( obropt == 2 .or. obropt == 12) THEN
!
!-----------------------------------------------------------------------
!
!  For obropt = 2,
!  we apply correction to wcont * rhostrw by assuming that
!  difz ( delta(wcont * rhostr) ) = A zs, where A is
!  a constant. Therefore  delta( wcont*rhostr ) = A zp**2/2.
!
!-----------------------------------------------------------------------
!
      WRITE(6,'(a)')                                                    &
          'Horizontal divergence correction distributed linearly in z.'

      DO k=3,ktop
        DO j=1,ny-1
          DO i=1,nx-1
            dz2=(tem1(i,j,k)-tem1(i,j,2))*(tem1(i,j,k)-tem1(i,j,2))
            wcont(i,j,k)=wcont(i,j,k)-                                  &
                         wcont(i,j,ktop)*(dz2*tem1(i,j,1))
          END DO
        END DO
      END DO

    ELSE IF( obropt == 3 .or. obropt == 13) THEN
!
!-----------------------------------------------------------------------
!
!  For obropt = 3,
!  We apply correction to wcont * rhostrw by assuming that
!  difz ( delta(wcont * rhostr) ) = A theta, where A is
!  a constant. Therefore  delta( wcont*rhostr ) = A theta**2/2.
!
!-----------------------------------------------------------------------
!
      WRITE(6,'(a)')                                                    &
          'Horizontal div. correction distributed linearly in theta.'

      DO k=3,ktop
        DO j=1,ny-1
          DO i=1,nx-1
            dth2=(tem2(i,j,k)-tem2(i,j,2))*(tem2(i,j,k)-tem2(i,j,2))
            wcont(i,j,k)=wcont(i,j,k)-                                  &
                       wcont(i,j,ktop)*(dth2*tem2(i,j,1))
          END DO
        END DO
      END DO
    ELSE
      WRITE(6,'(a)')                                                    &
          'No OBrien vertical velocity/mean divergence adjustment'
    END IF

    WRITE(6,'(a)')                                                      &
        ' After wcont adjustment, rho*wcont(top):'

    CALL a3dmax(wcont,1,nx,1,nx-1,                                      &
              1,ny,1,ny-1,1,nz,ktop,ktop,                               &
              wcmax,wcmin, imax,jmax,kmax, imin,jmin,kmin)

    WRITE(6,'(2(1x,a,g12.4,3(a,i3)))')                                  &
        'wcmin =',wcmin,' at i=',imin,', j=',jmin,', k=',kmin,          &
        'wcmax =',wcmax,' at i=',imax,', j=',jmax,', k=',kmax
!
!-----------------------------------------------------------------------
!
!  Obtain wcont from wcont*rhostrw.
!
!-----------------------------------------------------------------------
!
    DO k=2,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          wcont(i,j,k)=wcont(i,j,k)/rhostrw(i,j,k)
        END DO
      END DO
    END DO
    DO j=1,ny-1
      DO i=1,nx-1
        wcont(i,j,1 )=-wcont(i,j,   3)
        wcont(i,j,nz)=-wcont(i,j,nz-2)
      END DO
    END DO

!
!-----------------------------------------------------------------------
!
!  Apply w from cloud analysis, if desired.
!  Reset wcont
!
!-----------------------------------------------------------------------
!
    IF(cldwopt > 0) THEN

      DO k=1,nz-1
        DO j=1,ny-1
          DO i=1,nx-1
            w(i,j,k)=AMAX1(          END DO
        END DO
      END DO

      CALL wcontra(nx,ny,nz,u,v,w,mapfct,j1,j2,j3,aj3z,                 &
           rhostr,rhostru,rhostrv,rhostrw,wcont,tem1,tem2)


    END IF

    IF(wndadj == 3) THEN

      IF( nx*ny*nz < (nx+1)*(ny+1) ) THEN
        PRINT*,'Work array PHI in ADJUVW too small. Job stopped.'
        STOP
      END IF

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=1,nx-1
            tem1(i,j,k) = divh(i,j,k) +                                 &
                       (wcont(i,j,k+1)*rhostrw(i,j,k+1)-                &
                        wcont(i,j,k  )*rhostrw(i,j,k  ))*dzinv
          END DO
        END DO
      END DO

      PRINT*,' Div max/min before u-v adjustment.'

      CALL a3dmax(tem1,1,nx,1,nx-1,                                     &
              1,ny,1,ny-1,1,nz,2,nz-2,                                  &
              divmax,divmin, imax,jmax,kmax, imin,jmin,kmin)

      WRITE(6,'(2(1x,a,g25.12,3(a,i3)))')                               &
          'divmin =',divmin,' at i=',imin,', j=',jmin,', k=',kmin,      &
          'divmax =',divmax,' at i=',imax,', j=',jmax,', k=',kmax


      dxx = dx*dx
      dyy = dy*dy

      DO k=2,nz-2
        DO j=1,ny-1
          DO i=1,nx-1
            rhs(i,j,k) = - tem1(i,j,k)/(mapfct(i,j,1)**2)
          END DO
        END DO
      END DO

      pi = 4.0*ATAN(1.0)
      alpha = 4*(0.5 - pi/(2*SQRT(2.0))*                                &
          SQRT(1.0/(nx+1)**2+1.0/(ny+1)**2) )

      rdxx = 1.0/dxx
      rdyy = 1.0/dyy
      tem = alpha*dxx*dyy/(2*(dxx+dyy))

      PRINT*,'Over-relaxation coeff = ', alpha

      DO j=0,ny
        DO i=0,nx
          phi(i,j) = 0.0
        END DO
      END DO

      DO k=2,nz-2

        zk = (k-1.5)*dz
        IF(k > 2) THEN
          DO j=1,ny-1
            DO i=1,nx-1
              phi(i,j) = phi(i,j)*zk/(zk-dz)
            END DO
          END DO
        END IF
!
!-----------------------------------------------------------------------
!
!  Solve Possion equation for the potential of the corrective
!  horizontal velocity. SOR method is used.
!  phi = zero is assumed at i=0 and nx, j=0 and ny.
!
!-----------------------------------------------------------------------
!
        DO iterations = 1,1000

          absphi = 0.0
          delphi = 0.0

          DO j=1,ny-1
            DO i=1,nx-1
              rij = (phi(i+1,j)-2*phi(i,j)+phi(i-1,j))*rdxx+            &
                    (phi(i,j+1)-2*phi(i,j)+phi(i,j-1))*rdyy-            &
                    rhs(i,j,k)
              dphi = tem*rij

              phi(i,j) = phi(i,j) + dphi
              absphi = ABS(phi(i,j))+absphi
              delphi = ABS(  dphi  )+delphi
            END DO
          END DO

          IF( delphi/(absphi+1.0E-30 ) < 1.0E-5 ) EXIT

        END DO
!        710       CONTINUE

        PRINT*,'Iterations ended at', iterations,                       &
               ' for k=,',k,' error=',                                  &
                 delphi/(absphi+1.0E-30 )

!
!-----------------------------------------------------------------------
!
!  Apply the adjustment to u and v.
!
!-----------------------------------------------------------------------
!
        DO j=1,ny-1
          DO i=1,nx
            u(i,j,k)=u(i,j,k)+(phi(i,j)-phi(i-1,j))*mapfct(i,j,2)       &
                     /(rhostru(i,j,k)*dx)
          END DO
        END DO

        DO j=1,ny
          DO i=1,nx-1
            v(i,j,k)=v(i,j,k)+(phi(i,j)-phi(i,j-1))*mapfct(i,j,3)       &
                     /(rhostrv(i,j,k)*dy)
          END DO
        END DO

      END DO

      IF( chkmas ) THEN  ! Check mass continuity after adjustment

        DO k=2,nz-2
          DO j=1,ny-1
            DO i=1,nx-1
              tem1(i,j,k)=((u(i+1,j,k)*rhostru(i+1,j,k)*mapfct(i+1,j,5) &
                           -u(i  ,j,k)*rhostru(i  ,j,k)*mapfct(i  ,j,5)) &
                           *dxinv                                       &
                          +(v(i,j+1,k)*rhostrv(i,j+1,k)*mapfct(i,j+1,6) &
                           -v(i,j  ,k)*rhostrv(i,j  ,k)*mapfct(i,j  ,6)) &
                           *dyinv) *mapfct(i,j,1)*mapfct(i,j,1)         &
                           +(wcont(i,j,k+1)*rhostrw(i,j,k+1)            &
                            -wcont(i,j,k  )*rhostrw(i,j,k  ))*dzinv
            END DO
          END DO
        END DO

        PRINT*,' Div max/min after u-v adjustment.'

        CALL a3dmax(tem1,1,nx,1,nx-1,                                   &
                1,ny,1,ny-1,1,nz,2,nz-2,                                &
                divmax,divmin, imax,jmax,kmax, imin,jmin,kmin)

        WRITE(6,'(2(1x,a,g25.12,3(a,i3)))')                             &
            'divmin =',divmin,' at i=',imin,', j=',jmin,', k=',kmin,    &
            'divmax =',divmax,' at i=',imax,', j=',jmax,', k=',kmax

      END IF
    END IF

    CALL wctow(nx,ny,nz,u,v,wcont,mapfct,                               &
               j1,j2,j3,aj3z,rhostr,rhostru,rhostrv,rhostrw,1,1,        &
               w,                                                       &
               tem1,tem2)

  END IF

  RETURN
END SUBROUTINE adjuvw