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


SUBROUTINE microph(nx,ny,nz, dtbig1,                                    & 1,4
           ptprt,pprt,qv,qc,qr,qi,qs,qh,raing,prcrate,                  &
           rhostr,pbar,ptbar,ppi,j3,j3inv,                              &
           rhobar,p,lathvt,tem1,tem2,tem3)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate and apply the microphysical contributions to the water
!  and temperature fields. The current version implements the Kessler
!  warm rain microphysics parameterization scheme. Ice processes are
!  not included in this version.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/1991.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/28/92 (K. Brewster)
!  Further facelift.
!
!  11/08/98 (K. Brewster)
!  Added calculation of total p and pi for use in revap and satadj
!
!-----------------------------------------------------------------------
!
!  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
!
!    dtbig1   Big time step size (s)
!    ptprt    Perturbation potential temperature at all time levels (K)
!    pprt     Perturbation pressure at all time levels (Pascal)
!    qv       Water vapor specific humidity at all time levels (kg/kg)
!    qc       Cloud water mixing ratio at all time levels (kg/kg)
!    qr       Rainwater mixing ratio at all time levels (kg/kg)
!    qi       Cloud ice mixing ratio at all time levels (kg/kg)
!    qs       Snow mixing ratio at all time levels (kg/kg)
!    qh       Hail mixing ratio at all time levels (kg/kg)
!    raing    Accumulated grid-scale rainfall (mm)
!
!    rhostr   Base state air density times j3 (kg/m**3)
!    pbar     Base state pressure (Pascal)
!    ptbar    Base state potential temperature (K)
!    ppi      Exner function at tfuture
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!
!  OUTPUT:
!
!    ptprt    Perturbation potential temperature at time tfuture (K)
!    pprt     Perturbation pressure at time tfuture (Pascal)
!    qv       Water vapor specific humidity at time tfuture (kg/kg)
!    qc       Cloud water mixing ratio at time tfuture (kg/kg)
!    qr       Rainwater mixing ratio at time tfuture (kg/kg)
!    qi       Cloud ice mixing ratio at time tfuture (kg/kg)
!    qs       Snow mixing ratio at time tfuture (kg/kg)
!    qh       Hail mixing ratio at time tfuture (kg/kg)
!    raing    Accumulated grid-scale rainfall (mm)
!    prcrate  Precipitation rate (kg/(m**2*s))
!
!  WORK ARRAYS:
!
!    rhobar   Base state air density (kg/m**3)
!    lathvt   Temperature-dependent latent heat of evaporation
!    p        Total pressure at tfuture
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INCLUDE 'timelvls.inc'
!
  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions
!
  REAL :: dtbig1               ! Big time step size (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)  ! Rainwater 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 :: raing (nx,ny)        ! Accumulated grid-scale rainfall (mm)
  REAL :: prcrate(nx,ny)       ! Precipitation rate (kg/(m**2*s))


  REAL :: rhostr(nx,ny,nz)     ! Base state air density times j3 (kg/m**3)
  REAL :: pbar  (nx,ny,nz)     ! Base state pressure (Pascal).
  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! d( zp )/d( z ).
  REAL :: ppi   (nx,ny,nz)     ! Exner function at tfuture
  REAL :: j3inv (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! d( zp )/d( z ).
!
!-----------------------------------------------------------------------
!
!  Temporary arrays
!
!-----------------------------------------------------------------------
!
  REAL :: p     (nx,ny,nz)    ! Total pressure at tfuture (Pa)
  REAL :: rhobar(nx,ny,nz)    ! Base state air density (kg/m**3)
  REAL :: lathvt(nx,ny,nz)    ! Temperature dependent latent heat
                              ! of evaporation. Local array

  REAL :: tem1(nx,ny,nz)      ! Temporary work array
  REAL :: tem2(nx,ny,nz)      ! Temporary work array
  REAL :: tem3(nx,ny,nz)      ! Temporary work array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
!
  INTEGER :: i,j,k,INDEX
  REAL :: tk
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!-----------------------------------------------------------------------
!
!  Since rhobar is used most often in the microphysics parameterizations,
!  we calculate and save rhobar for later use.
!
!-----------------------------------------------------------------------
!
  DO k = 1,nz-1
    DO j = 1,ny-1
      DO i = 1,nx-1
        rhobar(i,j,k) = rhostr(i,j,k)*j3inv(i,j,k)
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  To remove negative mixing ratios, which result from computational
!  inaccuracies in the advection process, we set all negative mixing
!  ratios to zero. This is an artificial adjustment and, as a result,
!  total water will not be conserved. The adjustment can be averted
!  by enhancing the numerical accuracy of subsequent model versions.
!
!-----------------------------------------------------------------------
!
  DO k = 1,nz
    DO j = 1,ny
      DO i = 1,nx
        qv(i,j,k,tfuture) = MAX(0.0,qv(i,j,k,tfuture))
        qc(i,j,k,tfuture) = MAX(0.0,qc(i,j,k,tfuture))
        qr(i,j,k,tfuture) = MAX(0.0,qr(i,j,k,tfuture))
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Use the look up table to calculate the temperature dependent
!  latent heat of evaporation lathvt.
!  Also calculate Exner function.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1

        p(i,j,k) = pprt(i,j,k,tfuture)+pbar(i,j,k)
        tk = (ptprt(i,j,k,tfuture)+ptbar(i,j,k)) * ppi(i,j,k)
        INDEX = MAX( MIN(151, nint(tk-172.15)), 1)
        lathvt(i,j,k) = latheatv(INDEX)

      END DO
    END DO
  END DO

  IF( mphyopt /= 0 .OR. cnvctopt == 2) THEN

!
!
!-----------------------------------------------------------------------
!
!  Calculate the autoconversion of cloud water to rainwater
!  and the accretion (collection) of cloud droplets by rain drops.
!
!-----------------------------------------------------------------------
!

    CALL autocac(nx,ny,nz, dtbig1, qc,qr)
!
!
!-----------------------------------------------------------------------
!
!  Rainwater evaporates when the air is subsaturated. Under these
!  conditions, water vapor specific humidity increases at the expense
!  of rainwater and the air is subsequently cooled due to evaporation.
!
!  The evaporation rate is subject to certain specified limits.
!  Rainwater, water vapor and temperature fields are then adjusted.
!
!-----------------------------------------------------------------------
!

    CALL revap(nx,ny,nz, dtbig1,ptprt,qv,qr,                            &
               rhobar,ptbar,p,ppi,lathvt, tem1,tem2)

!
!-----------------------------------------------------------------------
!
!                        RAINWATER SEDIMENTATION
!
!  Calculate the rain fallout rate and apply it to the rain
!  water field.
!
!  Since the rainwater fallout is an advective process,
!  negative values of rainwater can be generated by the finite
!  difference scheme. If this occurs, we set the negative water
!  content to zero.
!
!-----------------------------------------------------------------------
!
!call test_check_in("Bqrfall")
    CALL qrfall (nx,ny,nz,dtbig1,                                       &
                 qr,rhobar,j3,j3inv,raing,prcrate,                      &
                 tem1,tem2,tem3)
!call test_check_in("Aqrfall")
!
    DO k = 1,nz
      DO j = 1,ny
        DO i = 1,nx
          qr(i,j,k,tfuture) = MAX(0.0,qr(i,j,k,tfuture))
        END DO
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!                     SATURATION ADJUSTMENT
!             (CONDENSATION AND EVAPORATION BETWEEN QV AND QC)
!
!  Under the conditions of supersaturation, condensation occurs
!  which converts water vapor to cloud water.
!  .
!
!  If the air is subsaturated, cloud water evaporates and becomes
!  water vapor until either the air is saturated or the cloud water
!  has completely evaporated.
!
!  During these processes, the air temperature changes.
!
!  Adjust the potential temperature, water vapor and cloud water
!  accordingly.
!
!-----------------------------------------------------------------------
!
!
  END IF

  CALL satadj(nx,ny,nz, ptprt,qv,qc,ptbar,                              &
              p,ppi,lathvt,  tem1,tem2,tem3)

  DO k = 1,nz
    DO j = 1,ny
      DO i = 1,nx
        qv(i,j,k,tfuture) = MAX(0.0,qv(i,j,k,tfuture))
        qc(i,j,k,tfuture) = MAX(0.0,qc(i,j,k,tfuture))
        qr(i,j,k,tfuture) = MAX(0.0,qr(i,j,k,tfuture))
      END DO
    END DO
  END DO

!
  RETURN
END SUBROUTINE microph
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE AUTOCAC                    ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE autocac(nx,ny,nz, dtbig1, qc,qr) 1

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the autoconversion of cloud water to rainwater and the
!  accretion (collection) of cloud droplets by rain drops.
!  First calculate the conversion and accretion rates, then impose
!  limits on the amount of conversion and accretion. Adjust the
!  cloud and rainwater fields accordingly.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/1991.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/28/92 (K. Brewster)
!  Further facelift.
!
!  6/10/94 (M Xue & A. Sathye)
!  Power calculation for the accretion rate is replaced by a
!  lookup table function.
!
!  07/10/97 (Fanyou Kong - CMRP)
!  Include MPDCD advection option (sadvopt = 5)
!
!  11/18/98 (Keith Brewster)
!  Changed pibar to ppi (full pi).
!
!-----------------------------------------------------------------------
!
!  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
!
!    qc       Cloud water mixing ratio at all time levels (kg/kg)
!    qr       Rainwater mixing ratio at all time levels (kg/kg)
!
!  OUTPUT:
!
!    qc       Cloud water mixing ratio at time tfuture (kg/kg)
!    qr       Rainwater mixing ratio at time tfuture (kg/kg)
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  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.

  PARAMETER (nt=3, tpast=1, tpresent=2, tfuture=3)
!
  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  REAL :: dtbig1               ! Big time step size (s)

  REAL :: qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (kg/kg)
  REAL :: qr    (nx,ny,nz,nt)  ! Rainwater mixing ratio (kg/kg)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
!
  INTEGER :: i,j,k,lvlq
  REAL :: deltat,qcplus,qrplus,ar,arcrdt,cr
  REAL :: temp, f1, f2, rstep
  INTEGER :: INDEX
!
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!

  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  lvlq = tfuture
!
!-----------------------------------------------------------------------
!
!  When water and temperture fields are advected using leapfrog scheme:
!  (if using forward scheme, deltat = dtbig1.)
!
!-----------------------------------------------------------------------
!
!C    IF( sadvopt.ge.1.and.sadvopt.le.3.or.sadvopt.eq.5) THEN ! Leapfrog scheme
  IF( sadvopt /= 4) THEN                  ! Leapfrog scheme
    deltat = 2*dtbig1
  ELSE                                    ! Forward scheme
    deltat = dtbig1
  END IF
!
!-----------------------------------------------------------------------
!
!  Autoconversion rate of cloud water to rainwater:
!
!      AR = autort * (QC - autotr)
!
!  autort : Autoconversion rate parameter (defined in globcst.inc)
!  autotr:  Threshold value for the autoconversion process to begin.
!
!
!  Accretion (collection) rate of cloud water by rainwater
!
!      CR = accrrt * QC * (QR ** 0.875)
!
!  accrrt:  Accretion rate multiplier (defined in globcst.inc)
!
!-----------------------------------------------------------------------

  rstep = 1.0/50.0E-6

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

!
!-----------------------------------------------------------------------
!
!  Take only the positive values of qc and qr.
!
!-----------------------------------------------------------------------
!
        qcplus = MAX(0.0, qc(i,j,k,lvlq) )
        qrplus = MAX(0.0, qr(i,j,k,lvlq) )
!
!-----------------------------------------------------------------------
!
!  The autoconversion rate:
!
!-----------------------------------------------------------------------
!
        ar = autort *( qcplus - autotr)
!
!-----------------------------------------------------------------------
!
!  Autoconversion only occurs when qc > autotr
!  i.e., AR calculated above should be >0.0
!
!-----------------------------------------------------------------------
!
        ar = MAX(0.0, ar)

        arcrdt = MIN( ar*deltat, qcplus )
        qc(i,j,k,lvlq) = qc(i,j,k,lvlq) - arcrdt
        qr(i,j,k,lvlq) = qr(i,j,k,lvlq) + arcrdt
        qcplus = MAX(0.0, qc(i,j,k,lvlq) )
        qrplus = MAX(0.0, qr(i,j,k,lvlq) )
!
!-----------------------------------------------------------------------
!
!  Accretion rate:
!
!      CR = accrrt * QC * (QR ** 0.875)
!
!  Here a linear interpolation from a lookup table data is used to
!  evaluate power function qr ** 0.875
!  for which it is assumed that 0.0=< qr =< 0.05.
!
!-----------------------------------------------------------------------
!

        temp = MIN(0.05, qrplus)*rstep
        INDEX = INT(temp)

        f1 = pwr875(INDEX)
        f2 = pwr875(INDEX + 1)

        cr = accrrt * qcplus * (f1 + ((f2 - f1) * (temp - INDEX) ))
!
!-----------------------------------------------------------------------
!
!  Calculate the total conversion from cloud water to rainwater due
!  to autoconversion and accretion. This should not exceed the
!  amount of cloud water present:
!
!-----------------------------------------------------------------------
!
!    arcrdt = min( (ar + cr)*deltat, qcplus )

        arcrdt = MIN( cr*deltat, qcplus )
!
!-----------------------------------------------------------------------
!
!  Update qc and qr to account for loss in cloud water
!  and gain in rainwater.
!
!-----------------------------------------------------------------------

        qc(i,j,k,lvlq) = qc(i,j,k,lvlq) - arcrdt

        qr(i,j,k,lvlq) = qr(i,j,k,lvlq) + arcrdt


      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE autocac
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE REVAP                      ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE revap(nx,ny,nz, dtbig1,ptprt,qv,qr,                          & 1,1
           rhobar,ptbar,p,ppi, lathvt, qvs, tem1)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the rainwater evaporation rate. Apply these changes
!  to QR and QV and then adjust the potential temperature PTPRT
!  due to the evaporative cooling.
!
!  Rainwater evaporates when the air is subsaturated. Under these
!  conditions, water vapor specific humidity increases at the expense
!  of rainwater and the air is subsequently cooled due to evaporation.
!
!  The evaporation rate is subject to certain specified limits.
!  Rainwater, water vapor and temperature fields are then adjusted.
!
!-----------------------------------------------------------------------

!
!  AUTHOR: Ming Xue
!  2/29/1992.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/28/92 (K. Brewster)
!  Further facelift.
!
!  6/10/94 (M Xue & A. Sathye)
!  Power calculations are replaced by lookup table functions.
!
!  11/08/98 (K. Brewster)
!  Changed calculation of qvs to be based on total pi and p, not pbar
!
!
!-----------------------------------------------------------------------
!
!  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
!
!    ptprt    Perturbation potential temperature at all time levels (K)
!    qv       Water vapor specific humidity at all time levels (kg/kg)
!    qr       Rainwater mixing ratio at all time levels (kg/kg)
!
!    ptbar    Base state potential temperature (K)
!    p        Total pressure at tfuture
!    ppi      Exner function of pressure
!    lathvt   Temperature-dependent latent heat of evaporation
!
!  OUTPUT:
!
!    ptprt    Perturbation potential temperature at time tfuture (K)
!    qv       Water vapor specific humidity at time tfuture (kg/kg)
!    qr       Rainwater mixing ratio at time tfuture (kg/kg)
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  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.

  PARAMETER (nt=3, tpast=1, tpresent=2, tfuture=3)
!
  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  REAL :: dtbig1               ! Big time step size (s)

  REAL :: ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
  REAL :: qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (kg/kg)
  REAL :: qr    (nx,ny,nz,nt)  ! Rainwater mixing ratio (kg/kg)

  REAL :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)
  REAL :: rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
!
  REAL :: p     (nx,ny,nz)     ! Total pressure at tfuture.
  REAL :: ppi   (nx,ny,nz)     ! Exner function of pressure.
  REAL :: lathvt(nx,ny,nz)     ! Temperature dependent latent heat

!-----------------------------------------------------------------------
!
!  Temporary arrays
!
!-----------------------------------------------------------------------
!
  REAL :: qvs   (nx,ny,nz)      ! Saturation specific humidity
                                ! with respect to liquid water.
  REAL :: tem1  (nx,ny,nz)      ! Temperature (K)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,lvlq,lvlt
  REAL :: deltat,qrplus,qvplus,c,er,erdt
  REAL :: interp, temp, f1, f2, rstep
  INTEGER :: INDEX
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

  lvlq = tfuture
  lvlt = tfuture
!
!-----------------------------------------------------------------------
!
!  The time increment when water and temperature fields are advected
!  using the leapfrog scheme (if using forward scheme, deltat =
!  dtbig1).
!
!-----------------------------------------------------------------------
!
!C    IF( sadvopt.ge.1.and.sadvopt.le.3) THEN ! Leapfrog scheme
  IF( sadvopt /= 4) THEN                  ! Leapfrog scheme
    deltat = 2*dtbig1
  ELSE                                    ! Forward scheme
    deltat = dtbig1
  END IF
!
!-----------------------------------------------------------------------
!
!   Calculate qvs, the saturation specific humidity, using
!   enhanced Teten's formula
!
!-----------------------------------------------------------------------

  DO k = 1,nz-1
    DO j = 1,ny-1
      DO i = 1,nx-1
        tem1(i,j,k) = (ptprt(i,j,k,lvlt)+ptbar(i,j,k)) * ppi(i,j,k)
      END DO
    END DO
  END DO

  CALL getqvs(nx,ny,nz, 1,nx-1,1,ny-1,1,nz-1, p,tem1, qvs)

! write (*,*) "ZUWEN subsatopt/rhsat", subsatopt, rhsat 
  IF (subsatopt /= 0) THEN 

    DO k = 1,nz-1
      DO j = 1,ny-1
        DO i = 1,nx-1
          qvs(i,j,k) = rhsat*qvs(i,j,k)   ! adjust qvs for sub-saturation
        END DO
      END DO
    END DO

  END IF 

!
!-----------------------------------------------------------------------
!
!  Calculate the amount of rainwater evaporation and adjust
!  the rainwater, water vapor and temperature fields.
!
!-----------------------------------------------------------------------
!
  rstep = 1.0/50.0E-6

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

        qrplus = MAX(0.0, qr(i,j,k,lvlq) )
        qvplus = MAX(0.0, qv(i,j,k,lvlq) )
!
!-----------------------------------------------------------------------
!
!  The formula for the calculation of the coefficient for the
!  rainwater evaporation rate
!
!    c = 1.6 + 30.3922 * ( rhobar(i,j,k)*qrplus )**0.2046.
!
!  Here a linear interpolation from a lookup table data is used to
!  evaluate power function (rhobar*qr) ** 0.2046
!  for which it is assumed that 0.0=< rhobar*qr =< 0.05.
!
!-----------------------------------------------------------------------
!
        temp = MIN (0.05, qrplus * rhobar(i,j,k)) * rstep
        INDEX = INT(temp)

        f1 = pwr2046(INDEX)
        f2 = pwr2046(INDEX + 1)

        interp= f1 + ((f2 - f1) * (temp - INDEX) )

        c = 1.6 + 30.3922 * interp
!
!-----------------------------------------------------------------------
!
!  The formula for calculation of the rainwater evaporation rate
!
!     (rhobar(i,j,k)*qrplus)**0.525
!
!  Here a linear interpolation from a lookup table data is used to
!  evaluate power function (rhobar*qr) ** 0.525
!  for which it is assumed that 0.0=< rhobar*qr =< 0.05.
!
!-----------------------------------------------------------------------
!

        f1 = pwr525(INDEX)
        f2 = pwr525(INDEX + 1)

        interp= f1 + ((f2 - f1) * (temp - INDEX) )

        er = c * (1.0- qvplus/qvs(i,j,k) )*                             &
                interp/( (2.03E4 + 9.584E6/(p(i,j,k)*qvs(i,j,k)))       &
                *rhobar(i,j,k) )
!
!-----------------------------------------------------------------------
!
!  The amount of rain evaporation is limited by the available
!  rainwater:
!
!-----------------------------------------------------------------------
!
        erdt = er * deltat
        erdt = MIN( qrplus, MAX(0.0, erdt) )
!
!-----------------------------------------------------------------------
!
!  Adjust qr, qv and ptprt accordingly.
!
!-----------------------------------------------------------------------
!
        qr(i,j,k,lvlq) = qr(i,j,k,lvlq) - erdt

        qv(i,j,k,lvlq) = qv(i,j,k,lvlq) + erdt

        ptprt(i,j,k,lvlt) = ptprt(i,j,k,lvlt)                           &
              - lathvt(i,j,k)*erdt/( ppi(i,j,k)*cp )

      END DO
    END DO
  END DO


  RETURN
END SUBROUTINE revap
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE QRFALL                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE qrfall(nx,ny,nz,dtbig1,                                      & 1,1
           qr,rhobar,j3,j3inv,raing, prcrate,                           &
           draing,vtr,tem1)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the fall-out rate of rainwater and apply this effect to
!  the rainwater field, qr. Since the leapfrog-centered scheme is
!  used, the rainwater flux divergence (i.e., the advection by the
!  terminal fallout speed of rain) is calculated at time ""tpresent".
!  This subroutine is called by the warm rain package.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/1991.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/28/92 (K. Brewster)
!  Further facelift.
!
!  6/10/94 (M Xue & A. Sathye)
!  Power calculations are replaced by lookup table functions.
!
!  5/8/1995 (M. Xue)
!  Time-splitting upstream advection used for rainwater sedimentation
!
!  10/20/1996 (M. Xue)
!  Created a general routine for hydrometeor sedimentation called
!  QFALLOUT, which is now also called by ice microphysics routine.
!  Precipitation rate array prcrate is now correctly calculated
!  for split-step integration.
!
!-----------------------------------------------------------------------
!
!  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
!
!    dtbig1   Big time step size (s)
!    qr       Rainwater mixing ratio at all time levels (kg/kg)
!
!    rhobar   Base state air density (kg/m**3)
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!    raing    Accumulated grid-scale rainfall (mm)
!
!  OUTPUT:
!
!    qr       Rainwater mixing ratio at time tfuture (kg/kg)
!    raing    Accumulated grid-scale rainfall (mm)
!    prcrate  Precipitation rate
!
!  WORK ARRAYS:
!
!    draing   Gridscale rainfall (mm)
!    vtr      terminal velocity
!    tem1     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
!
  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.

  PARAMETER (nt=3, tpast=1, tpresent=2, tfuture=3)
!
  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  REAL :: dtbig1               ! Big time step size (s)
  REAL :: qr    (nx,ny,nz,nt)  ! Rainwater mixing ratio (kg/kg)
  REAL :: rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! defined as d( zp )/d( z ).
  REAL :: j3inv (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! defined as d( zp )/d( z ).
  REAL :: raing (nx,ny)        ! Accumulated grid-scale rainfall (mm)
  REAL :: prcrate(nx,ny)       ! Precipitation rate (kg/(m**2*s))
!
!-----------------------------------------------------------------------
!
!  Temporary arrays
!
!-----------------------------------------------------------------------
!
  REAL :: draing(nx,ny)       ! Gridscale rainfall (mm)
  REAL :: vtr   (nx,ny,nz)    ! terminal velocify
  REAL :: tem1  (nx,ny,nz)    ! work array

!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,lvlq
  REAL :: temp, interp, f1, f2, rstep
  INTEGER :: INDEX
  REAL :: denwater,tem
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  lvlq = tfuture
!
!-----------------------------------------------------------------------
!
!  The formula for terminal fallout speed of rainwater:
!
!  Vtr = 36.34 * ((0.001 * rhobar * qr) ** 0.1364) * (rho0 /rhobar)**0.5
!  Vtr is positive downwards.
!
!  A linear interpolation from a lookup table data is used to
!  evaluate power function (0.001 * rhobar * qr) ** 0.1364
!  for which it is assumed that 0.0=< rhobar * qr =< 0.05.
!
!
!-----------------------------------------------------------------------
!
  rstep = 1.0/50.0E-9

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

        temp = MIN (0.00005, MAX(0.0,0.001 * rhobar(i,j,k)              &
              * qr(i,j,k,lvlq-1)))*rstep
        INDEX = INT(temp)

        f1 = pwr1364(INDEX)
        f2 = pwr1364(INDEX+1)

        interp= f1 + ((f2 - f1) * (temp - INDEX) )

        vtr(i,j,k) = 36.34 *  interp * SQRT(rho0/rhobar(i,j,k))

      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate the rainwater fallout term and the precipitation rate.
!
!-----------------------------------------------------------------------
!

!call test_check_in("Bqrfallout")
  CALL qfallout(nx,ny,nz,dtbig1,qr,rhobar,j3,j3inv,draing,vtr,          &
                tem1)
!call test_check_in("Aqrfallout")

  denwater = 1000.0   ! Density of liquid water (kg/m**3)
  tem = denwater/(1000.0*dtbig)

  DO j=1,ny-1
    DO i=1,nx-1
      raing(i,j)=raing(i,j)+draing(i,j)
      prcrate(i,j) = draing(i,j)*tem
    END DO
  END DO

  RETURN
END SUBROUTINE qrfall
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE QFALLOUT                   ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE qfallout(nx,ny,nz,dtbig1,q,rhobar,j3,j3inv,draing,vtr,       & 6,3
           rhovq)

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the fall-out of hydrometeor given it's fallout velocity.
!  Upstream advection with split time steps is used here.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/1991.
!
!  MODIFICATION HISTORY:
!
!  5/8/1995 (M. Xue)
!  Time-splitting upstream advection used for rainwater sedimentation
!
!  12/6/1996 (Yuhe Liu)
!  Moved the initialization of draing to the beginning of executable
!  statements. This solved the problem that the returned draing
!  became undefined when no rain sedimentation needed.
!
!-----------------------------------------------------------------------
!
!  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
!
!    dtbig1   Big time step size (s)
!    q        Hydrometeor mixing ratio at all time levels (kg/kg)
!
!    rhobar   Base state air density (kg/m**3)
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!    vtr      terminal velocity (m/s)
!
!  OUTPUT:
!
!    q        Rainwater mixing ratio at time tfuture (kg/kg)
!    draing   Grid-scale rainfall (mm)
!
!  WORK ARRAYS:
!
!    rhovq    rhobar*terminal velocity*q
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  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.

  PARAMETER (nt=3, tpast=1, tpresent=2, tfuture=3)
!
  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  REAL :: dtbig1               ! Big time step size (s)
  REAL :: q     (nx,ny,nz,nt)  ! Rainwater mixing ratio (kg/kg)
  REAL :: rhobar(nx,ny,nz)     ! Base state air density (kg/m**3)
  REAL :: j3    (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! defined as d( zp )/d( z ).
  REAL :: j3inv (nx,ny,nz)     ! Coordinate transformation Jacobian
                               ! defined as d( zp )/d( z ).
  REAL :: draing (nx,ny)       ! Grid-scale rainfall (mm)
  REAL :: vtr   (nx,ny,nz)     ! terminal velocify
!
!-----------------------------------------------------------------------
!
!  Temporary arrays
!
!-----------------------------------------------------------------------
!
  REAL :: rhovq (nx,ny,nz)    ! rhobar*terminal velocity*q

!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,nrnstp,irnstp
  REAL :: denwater,tem, vtrdzmax,dtrnf,dtrnf0
  REAL :: deltat
  REAL :: temmin
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'grid.inc'          ! Grid & map parameters.
  INCLUDE 'phycst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  CALL acct_interrupt(qfall_acct)

  DO j=1,ny-1
    DO i=1,nx-1
      draing(i,j)=0.0
    END DO
  END DO

  vtrdzmax = 0.0
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        vtrdzmax=MAX(vtrdzmax,vtr(i,j,k)*(j3inv(i,j,k)*dzinv))
      END DO
    END DO
  END DO

!wdt 2001-04-17 gmb: slows execution down significantly:
!  ! for bit-for-bit MP accuracy:   !wdt - don't use!
  temmin = 0
  CALL mpmax0(vtrdzmax,temmin)

  IF( vtrdzmax == 0.0 ) RETURN  ! vtr=0, no rain sedimentation.

  IF( sadvopt /= 4) THEN                  ! Leapfrog scheme
    deltat = 2*dtbig1
  ELSE                                    ! Forward scheme
    deltat = dtbig1
  END IF

  dtrnf = MIN(deltat, 0.5/MAX(1.0E-20,vtrdzmax))
                         ! 0.5 is the max Courant number

  dtrnf0 = dtrnf
  nrnstp = nint ( deltat/dtrnf )
  IF(nrnstp == 0) nrnstp=1

  dtrnf  = deltat/nrnstp
  IF (dtrnf > dtrnf0) THEN
    nrnstp = nrnstp + 1
    dtrnf  = deltat/nrnstp
  END IF

  DO irnstp = 1, nrnstp
!
!-----------------------------------------------------------------------
!
!  The density weighted rainwater fallout flux:
!  rhovq = rhobar * vtr * q
!
!-----------------------------------------------------------------------
!
    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          rhovq(i,j,k) = rhobar(i,j,k)*vtr(i,j,k)*q(i,j,k,tfuture)
        END DO
      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  Upstream-foreward advection for the rainwater sedimentation
!
!-----------------------------------------------------------------------
!
    DO k=2,nz-2
      DO j=1,ny-1
        DO i=1,nx-1
          q(i,j,k,tfuture) = q(i,j,k,tfuture)+ dtrnf *                  &
              ((rhovq(i,j,k+1)-rhovq(i,j,k))/(dz*rhobar(i,j,k)*j3(i,j,k)))
        END DO
      END DO
    END DO
!call test_dump2(q(1,1,1,tfuture),"XXXqfallout_q",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(rhovq,"XXXqfallout_rhovq",nx,ny,nz,1,nx-1,2,ny-2,2,nz-1)
!call test_check_in("qfallout1")

!-----------------------------------------------------------------------
!
!  Accumulated precipitation.  The depth of water
!  accumulated per timestep is given by:
!
!  depth (m) = terminal velocity * q * dt * (rhobar/denwater)
!
!  This equation is derived by temporally integrating the vertical
!  flux of rain through the lowest model level.  The total is given
!  for the entire computational domain.
!
!-----------------------------------------------------------------------

    denwater = 1000.0   ! Density of liquid water.


    tem = dtbig1/deltat *dtrnf/denwater *1000.0
                          ! dtbig1/deltat=0.5 for leapfrog scheme

    DO j=1,ny-1
      DO i=1,nx-1
        draing(i,j)=draing(i,j)+rhovq(i,j,2) * tem
      END DO
    END DO

  END DO

  CALL acct_stop_inter

  RETURN
END SUBROUTINE qfallout
!
!##################################################################
!##################################################################
!######                                                      ######
!######                SUBROUTINE SATADJ                     ######
!######                                                      ######
!######                     Developed by                     ######
!######     Center for Analysis and Prediction of Storms     ######
!######                University of Oklahoma                ######
!######                                                      ######
!##################################################################
!##################################################################
!


SUBROUTINE satadj(nx,ny,nz, ptprt,qv,qc,ptbar,                          & 1,2
           p, ppi, lathvt,  t,qvs, tem1 )

!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Perform the saturation adjustment (condensation and evaporation
!  between qv and qc). Under the conditions of supersaturation
!  (qv > qvs) water vapor, qv, is condensed into cloud water, qc.
!  When the air is subsaturated, (qv <= qvs) cloud water is then
!  evaporated. Adjustments are then made to the potential
!  temperature, water vapor and cloud water fields.

!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  2/29/1992.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/28/92 (K. Brewster)
!  Further facelift.
!
!  11/05/98 (K. Brewster)
!  Added "time constant" to dqv calculation.
!
!  11/08/98 (K. Brewster)
!  Changed calculation of qvs to be based on total pi and p, not pbar
!
!-----------------------------------------------------------------------
!
!  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
!
!    ptprt    Perturbation potential temperature at all time levels(K)
!    pprt     Perturbation pressure at all time levels (Pascal)
!    qv       Water vapor specific humidity at all time levels (kg/kg)
!    qc       Cloud water mixing ratio at all time levels (kg/kg)
!
!    ptbar    Base state potential temperature (K)
!
!  OUTPUT:
!
!    ptprt    Perturbation potential temperature at time tfuture (K)
!    qv       Water vapor specific humidity at time tfuture (kg/kg)
!    qc       Cloud water mixing ratio at time tfuture (kg/kg)
!
!  WORK ARRAYS:
!
!    qvs      Saturation specific humidity (kg/kg)
!    t        Temperature (K)
!
!   (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
!
  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.

  PARAMETER (nt=3, tpast=1, tpresent=2, tfuture=3)
!
  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions

  REAL :: ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
  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 :: ptbar (nx,ny,nz)     ! Base state potential temperature (K)

  REAL :: lathvt(nx,ny,nz)     ! Temperature dependent latent heat of
                               ! vaporization (kg/(m*s**2).
!
!-----------------------------------------------------------------------
!
!  Temporary arrays
!
!-----------------------------------------------------------------------
!
  REAL :: p     (nx,ny,nz)     ! Total pressure at tfuture (Pa)
  REAL :: ppi   (nx,ny,nz)     ! Exner function
  REAL :: t     (nx,ny,nz)     ! Temperature (K)
  REAL :: qvs   (nx,ny,nz)     ! Saturation specific humidity (kg/kg)
  REAL :: tem1  (nx,ny,nz)     ! Temporary arrays
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,lvlq,lvlt
  REAL :: dqv,qvplus,qcplus
  REAL :: dqvcst
  PARAMETER (dqvcst=1.0)
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Function f_desdt and inline directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_desdt

!fpp$ expand (f_desdt)
!dir$ inline always f_desdt
!*$*  inline routine (f_desdt)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  lvlq = tfuture
  lvlt = tfuture
!
!-----------------------------------------------------------------------
!
!  Calculate the temperature.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        t(i,j,k) = (ptprt(i,j,k,lvlt)+ptbar(i,j,k)) * ppi(i,j,k)
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate the saturation specific humidity.
!
!-----------------------------------------------------------------------
!
  CALL getqvs(nx,ny,nz, 1,nx-1,1,ny-1,1,nz-1, p, t, qvs)

! write (*,*) "ZUWEN subsatopt/rhsat", subsatopt, rhsat 
  IF (subsatopt /= 0) THEN 

    DO k=1,nz-1
      DO j=1,ny-1
        DO i=1,nx-1
          qvs(i,j,k) = rhsat*qvs(i,j,k)  ! adjust qvs for sub-saturation
        END DO
      END DO
    END DO

  END IF 

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tem1(i,j,k) = f_desdt( t(i,j,k) )     ! d(es)/dt/es
        tem1(i,j,k) = (rddrv+(1.-rddrv)*qvs(i,j,k))/rddrv * tem1(i,j,k)
                                              ! d(qvs)/dt/qvs
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
!
!-----------------------------------------------------------------------
!
!  Calculate the amount of evaporation (when subsaturated, dqv>0)
!  or condensation (when supersaturated, dqv<0).
!
!-----------------------------------------------------------------------
!
        qvplus = MAX(0.0,qv(i,j,k,lvlq))
        qcplus = MAX(0.0,qc(i,j,k,lvlq))

        dqv = (qvs(i,j,k)-qvplus)                                       &
            / (1.0+tem1(i,j,k)*qvs(i,j,k)*lathvt(i,j,k)/cp)
!
!-----------------------------------------------------------------------
!
!  If evaporation occurs (dqv>0), the amount should not exceed
!  the available cloud water.
!
!  Make evaporation and condensation happen gradually.
!
!-----------------------------------------------------------------------
!
        dqv = dqvcst*MIN( dqv, qcplus )
!
!-----------------------------------------------------------------------
!
!  Adjust qv, qc and ptprt accordingly.
!
!-----------------------------------------------------------------------
!
        qv(i,j,k,lvlq) = qv(i,j,k,lvlq) + dqv

        qc(i,j,k,lvlq) = qc(i,j,k,lvlq) - dqv

        ptprt(i,j,k,lvlt) = ptprt(i,j,k,lvlt) -                         &
              dqv*lathvt(i,j,k)/(ppi(i,j,k)*cp)


      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE satadj