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


SUBROUTINE microph_nem(nx,ny,nz, dtbig1,                                & 1,27
           t,p,qv,qc,qr,qi,qs,qh,raing,prcrate,                         &
           rhostr,pbar,ptbar,ppi,j3,j3inv,                              &
           rho,rsat,rliq,rice,x,t1,                                     &
           change,rate,maxmelt,rp_nuc,                                  &
           needi,needr,needs,tvq,tem15,tem16)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Paul Schultz's (FSL) Ice microphysics scheme.  Contains
!  equations for water vapor, cloud water, rain water, ice,
!  snow, and hail.  The scheme is based on:
!
!  Schultz, P., 1995: An Explicit Cloud Physics Parameterization
!  for Operational Numerical Weather Prediction.  Monthly Weather
!  Review, 123, 3331-3343.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Jason J. Levit
!  10/01/1996.
!
!  MODIFICATION HISTORY:
!
!  02/19/97 (J. Levit)
!  A major bug was found in the diffusional crystal growth
!  equation.  The conversion parameter 'c2p' should have been
!  'v2p'.
!
!  01/27/97 (J. Levit)
!  Added a variable 'delt', which is equal to 2.0*dtbig1, to fix
!  a bug associated with the temporal finite differencing in some of
!  the subroutines.  Also, uncommented out the code which allows
!  riming to occur.
!  Cleaned up some more of the code and documentation.
!
!  12/09/1996 (J. Levit)
!  Added some extra code at the end of the subroutine to
!  account for generation of negative values in the precipitation
!  fields.  The terminal velocity (tqv) was set equal to
!  zero before calculation, and a DO LOOP was added to
!  take only the positive values of the precipitation fields
!  once the fallout was calculated.
!
!  11/15/1996 (J. Levit)
!  Completed conversion from RAMS code into ARPS format.
!  Cleaned up code and documentation.
!
!  07/10/1997 (Fanyou Kong)
!  Include MPDCD advection option (sadvopt = 5)
!
!-----------------------------------------------------------------------
!
!  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 of pressure
!    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:
!
!    tem1     Temporary work array.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  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 :: t (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
  REAL :: p (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 :: ppi   (nx,ny,nz)     ! Exner function.
  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 ).
!
!-----------------------------------------------------------------------
!
!  Temporary arrays
!
!-----------------------------------------------------------------------
!
  REAL :: rho(nx,ny,nz)           ! Base state air density (kg/m**3)
  REAL :: rsat    (nx,ny,nz)      ! Temporary work array
  REAL :: rliq    (nx,ny,nz)      ! Temporary work array
  REAL :: rice    (nx,ny,nz)      ! Temporary work array
  REAL :: x       (nx,ny,nz)      ! Temporary work array
  REAL :: t1      (nx,ny,nz)      ! Temporary work array
  REAL :: change  (nx,ny,nz)      ! Temporary work array
  REAL :: rate    (nx,ny,nz)      ! Temporary work array
  REAL :: maxmelt (nx,ny,nz)      ! Temporary work array
  REAL :: rp_nuc  (nx,ny,nz)      ! Temporary work array
  REAL :: needi   (nx,ny,nz)      ! Temporary work array
  REAL :: needr   (nx,ny,nz)      ! Temporary work array
  REAL :: needs   (nx,ny,nz)      ! Temporary work array
  REAL :: tvq     (nx,ny,nz)      ! Temporary work array
  REAL :: tem15   (nx,ny,nz)      ! Temporary work array
  REAL :: tem16   (nx,ny,nz)      ! Temporary work array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
!
  INTEGER :: i,j,k                ! Temporary variable
  REAL :: denwater,tema,tem       ! Temporary variable
  REAL :: delt                    ! Temporary variable
  REAL :: eff                     ! Temporary variable
  REAL :: err_num,err_den,dr      ! Temporary variable
  INTEGER :: lvlq                 ! Temporayr variable
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'nemcst.inc'
  INCLUDE 'globcst.inc'
  INCLUDE 'phycst.inc'
!
!-----------------------------------------------------------------------
!
!  Function f_qvsatl, f_qvsati, f_desdtl, and f_desdti, and inline
!  directive for Cray PVP
!
!-----------------------------------------------------------------------
!
  REAL :: f_qvsatl, f_qvsati, f_desdtl, f_desdti

!fpp$ expand (f_qvsatl)
!fpp$ expand (f_qvsati)
!fpp$ expand (f_desdtl)
!fpp$ expand (f_desdti)
!dir$ inline always f_qvsatl, f_qvsati, f_desdtl, f_desdti
!*$*  inline routine (f_qvsatl, f_qvsati, f_desdtl, f_desdti)

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  lvlq=3

!C    IF( sadvopt.ge.1.and.sadvopt.le.3) THEN ! Leapfrog scheme
  IF( sadvopt /= 4) THEN                  ! Leapfrog scheme
    delt = 2*dtbig1
  ELSE                                    ! Forward scheme
    delt = dtbig1
  END IF

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        qv(i,j,k,lvlq)=MAX(0.0,qv(i,j,k,lvlq))
        qc(i,j,k,lvlq)=MAX(0.0,qc(i,j,k,lvlq))
        qr(i,j,k,lvlq)=MAX(0.0,qr(i,j,k,lvlq))
        qi(i,j,k,lvlq)=MAX(0.0,qi(i,j,k,lvlq))
        qs(i,j,k,lvlq)=MAX(0.0,qs(i,j,k,lvlq))
        qh(i,j,k,lvlq)=MAX(0.0,qh(i,j,k,lvlq))
        t(i,j,k,lvlq)=(t(i,j,k,lvlq)+ptbar(i,j,k))*ppi(i,j,k)
        p(i,j,k,lvlq)=(p(i,j,k,lvlq)+pbar(i,j,k))
        rho(i,j,k)=rhostr(i,j,k)*j3inv(i,j,k)
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Nothing in this routine can change the quantities til, which is
!  the temperature with corrections for latent heating by ice and
!  liquid, and rtot.  til is similar to theta-il of Tripoli and
!  Cotton (1981 MWR).  These are checked at the end.
!
!-----------------------------------------------------------------------
!

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        rliq(i,j,k) = qc(i,j,k,lvlq) + qr(i,j,k,lvlq)
        rice(i,j,k) = qi(i,j,k,lvlq) + qs(i,j,k,lvlq) + qh(i,j,k,lvlq)
      END DO
    END DO
  END DO

! rtot0 = qv(i,j,k,lvlq) + rliq(i,j,k) + rice(i,j,k)
! til0 = t  - (Lvl*rliq(i,j,k)+Lvi*rice(i,j,k))/cp

!
!-----------------------------------------------------------------------
!
!  Condensation and evaporation of liquid.  Cloud liquid is assumed
!  to occur instantaneously.  Evaporation of rain is as in Dudhia
!  and Moncrieff (JAS 89).  No condensational growth of rain.
!
!  The evaporation process eats up liquid before ice, and small
!  particles before large, so the order is cloud liquid, rain,
!  pristine crystals, snow, and finally ice.
!
!-----------------------------------------------------------------------
!
! write (*,*) "ZUWEN subsatopt/rhsat", subsatopt, rhsat 

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        rsat(i,j,k) = rhsat*f_qvsatl( p(i,j,k,lvlq), t(i,j,k,lvlq) )

        IF ( .NOT. (rliq(i,j,k) == 0. .AND.                             &
                      qv(i,j,k,lvlq) < rsat(i,j,k)) ) THEN

! Call SatAdjL (qv(i,j,k,lvlq), p(i,j,k,lvlq), t(i,j,k,lvlq), rsat(i,j,k), t1)

          t1(i,j,k) = t(i,j,k,lvlq)
          rsat(i,j,k) = rhsat*f_qvsatl( p(i,j,k,lvlq), t1(i,j,k) )
          dr = qv(i,j,k,lvlq) - rsat(i,j,k)
          err_num = t1(i,j,k) - t(i,j,k,lvlq) - lvl/cp * dr
          err_den = 1. + lvl/cp * rsat(i,j,k)      & ! 1+Lv1/cp*d(qvs)/dt
                   * (rddrv+(1.-rddrv)*rsat(i,j,k))                     &
                       * f_desdtl(t1(i,j,k))

          t1(i,j,k) = t1(i,j,k) - err_num/err_den
          rsat(i,j,k) = rhsat*f_qvsatl( p(i,j,k,lvlq), t1(i,j,k) )
          dr = qv(i,j,k,lvlq) - rsat(i,j,k)
          err_num = t1(i,j,k) - t(i,j,k,lvlq) - lvl/cp * dr
          err_den = 1. + lvl/cp * rsat(i,j,k)      & ! 1+Lv1/cp*d(qvs)/dt
                   * (rddrv+(1.-rddrv)*rsat(i,j,k))                     &
                       * f_desdtl(t1(i,j,k))

          t1(i,j,k) = t1(i,j,k) - err_num/err_den
          rsat(i,j,k) = rhsat*f_qvsatl( p(i,j,k,lvlq), t1(i,j,k) )
          dr = qv(i,j,k,lvlq) - rsat(i,j,k)
          err_num = t1(i,j,k) - t(i,j,k,lvlq) - lvl/cp * dr
          err_den = 1. + lvl/cp * rsat(i,j,k)      & ! 1+Lv1/cp*d(qvs)/dt
                   * (rddrv+(1.-rddrv)*rsat(i,j,k))                     &
                       * f_desdtl(t1(i,j,k))

          t1(i,j,k) = t1(i,j,k) - err_num/err_den
          rsat(i,j,k) = rhsat*f_qvsatl( p(i,j,k,lvlq), t1(i,j,k) )
          dr = qv(i,j,k,lvlq) - rsat(i,j,k)

          qc(i,j,k,lvlq) = qc(i,j,k,lvlq) + (qv(i,j,k,lvlq)-rsat(i,j,k))

          IF (qc(i,j,k,lvlq) > 0.) THEN

            qv(i,j,k,lvlq) = rsat(i,j,k)
            t(i,j,k,lvlq) = t1(i,j,k)

          ELSE ! evaporation of cloud liquid, then rain

            qv(i,j,k,lvlq) = rsat(i,j,k) + qc(i,j,k,lvlq)
            qc(i,j,k,lvlq) = 0.

            IF (qr(i,j,k,lvlq) > 0. .AND. qv(i,j,k,lvlq) < rsat(i,j,k)) THEN
              rate(i,j,k) = r2v * (rsat(i,j,k)-qv(i,j,k,lvlq))
              change(i,j,k) = rate(i,j,k) * delt
              IF (change(i,j,k) > (rsat(i,j,k)-qv(i,j,k,lvlq)))         &
                  change(i,j,k)=(rsat(i,j,k)-qv(i,j,k,lvlq))
              IF (change(i,j,k) > qr(i,j,k,lvlq)) THEN
                qv(i,j,k,lvlq) = qv(i,j,k,lvlq) + qr(i,j,k,lvlq)
                qr(i,j,k,lvlq) = 0.
              ELSE
                qv(i,j,k,lvlq) = qv(i,j,k,lvlq) + change(i,j,k)
                qr(i,j,k,lvlq) = qr(i,j,k,lvlq) - change(i,j,k)
              END IF
            END IF

            x(i,j,k) = qr(i,j,k,lvlq) - rliq(i,j,k)
            t(i,j,k,lvlq) = t(i,j,k,lvlq) + x(i,j,k)*lvl/cp

          END IF

          rliq(i,j,k) = qc(i,j,k,lvlq) + qr(i,j,k,lvlq)

        END IF

      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Ice growth and evaporation.  Ice nucleation is more a function of
!  supersaturation than temperature (the same is true of diffusional
!  ice crystal growth).  In the presence of cloud water, ice
!  supersaturation is greatest at -12C, and actually decreases at
!  very cold temperatures.
!
!  Ice nucleation generates ice mass much slower than diffusional
!  growth, so we don't compute that function it if there's any cloud
!  ice present.  Also, no nucleation in the absence of cloud water.
!
!  The diffusion from vapor to pristine crystals is proportional to
!  the vapor excess and is a function of the crystal mass already
!  there.  In the presence of cloud water, the excess is .17 g/kg
!   at 1000 mb and .85 g/kg at 200 mb.  This is a supersaturation of
!  12.4% at -12 C, the temperature at which the difference is
!  greatest.
!
!  Diffusional growth of snow is not allowed at this time.  The
!  water mass will get there anyway via crystal growth and
!  collection.
!
!  If the vapor transfer should be independent of pressure, the
!  equation should be rate(i,j,k) = V2P * (rho*(rv-rsat)) * (rho*rp)]
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        rsat(i,j,k) = rhsat*f_qvsati( p(i,j,k,lvlq), t(i,j,k,lvlq) )
        IF (.NOT.(rice(i,j,k) == 0. .AND. qv(i,j,k,lvlq) < rsat(i,j,k))) THEN

! Call SatAdjI (qv(i,j,k,lvlq), p, t, rsat(i,j,k), t1(i,j,k))

          t1(i,j,k) = t(i,j,k,lvlq)
          rsat(i,j,k) = rhsat*f_qvsati( p(i,j,k,lvlq), t1(i,j,k) )
          dr = qv(i,j,k,lvlq) - rsat(i,j,k)
          err_num = t1(i,j,k) - t(i,j,k,lvlq) - lvi/cp * dr
          err_den = 1. + lvl/cp * rsat(i,j,k)      & ! 1+Lv1/cp*d(qvs)/dt
                   * (rddrv+(1.-rddrv)*rsat(i,j,k))                     &
                       * f_desdti(t1(i,j,k))

          t1(i,j,k) = t1(i,j,k) - err_num/err_den
          rsat(i,j,k) = rhsat*f_qvsati( p(i,j,k,lvlq), t1(i,j,k) )
          dr = qv(i,j,k,lvlq) - rsat(i,j,k)
          err_num = t1(i,j,k) - t(i,j,k,lvlq) - lvi/cp * dr
          err_den = 1. + lvl/cp * rsat(i,j,k)      & ! 1+Lv1/cp*d(qvs)/dt
                   * (rddrv+(1.-rddrv)*rsat(i,j,k))                     &
                       * f_desdti(t1(i,j,k))

          t1(i,j,k) = t1(i,j,k) - err_num/err_den
          rsat(i,j,k) = rhsat*f_qvsati( p(i,j,k,lvlq), t1(i,j,k) )
          dr = qv(i,j,k,lvlq) - rsat(i,j,k)
          err_num = t1(i,j,k) - t(i,j,k,lvlq) - lvi/cp * dr
          err_den = 1. + lvl/cp * rsat(i,j,k)      & ! 1+Lv1/cp*d(qvs)/dt
                   * (rddrv+(1.-rddrv)*rsat(i,j,k))                     &
                       * f_desdti(t1(i,j,k))

          t1(i,j,k) = t1(i,j,k) - err_num/err_den
          rsat(i,j,k) = rhsat*f_qvsati( p(i,j,k,lvlq), t1(i,j,k) )
          dr = qv(i,j,k,lvlq) - rsat(i,j,k)

          IF (qv(i,j,k,lvlq) > rsat(i,j,k)) THEN ! ice growth

            IF (qc(i,j,k,lvlq) > 0. .AND. qi(i,j,k,lvlq) < 1E-6) THEN

!    Call Nucleate_pristine (qv(i,j,k,lvlq), rsat(i,j,k),
!  :       t(i,j,k,lvlq), rho(i,j,k), rp_nuc(i,j,k))

              rp_nuc(i,j,k) = 0.
              IF (.NOT.(t(i,j,k,lvlq) > 268.)) THEN

                IF (rp_nuc(i,j,k) > qi(i,j,k,lvlq)) THEN
                  qv(i,j,k,lvlq) = qv(i,j,k,lvlq) + qi(i,j,k,lvlq)
                  qi(i,j,k,lvlq) = rp_nuc(i,j,k)
                  qv(i,j,k,lvlq) = qv(i,j,k,lvlq) - rp_nuc(i,j,k)
                END IF

              ELSE

                rp_nuc(i,j,k) = 1E3 * EXP(-.639+12.96*(qv(i,j,k,lvlq)   &
                                          / rsat(i,j,k)-1.))            &
                              * pmas / rho(i,j,k)

                IF (rp_nuc(i,j,k) > (qv(i,j,k,lvlq)-rsat(i,j,k))/2.)    &
                    rp_nuc(i,j,k) = (qv(i,j,k,lvlq)-rsat(i,j,k))/2.

                IF (rp_nuc(i,j,k) > qi(i,j,k,lvlq)) THEN
                  qv(i,j,k,lvlq) = qv(i,j,k,lvlq) + qi(i,j,k,lvlq)
                  qi(i,j,k,lvlq) = rp_nuc(i,j,k)
                  qv(i,j,k,lvlq) = qv(i,j,k,lvlq) - rp_nuc(i,j,k)
                END IF
              END IF
            END IF

            IF (qi(i,j,k,lvlq) > 0. .AND. qv(i,j,k,lvlq) > rsat(i,j,k)) THEN
              rate(i,j,k) = v2p*(qv(i,j,k,lvlq)-rsat(i,j,k)) *          &
                            (qi(i,j,k,lvlq)*rho(i,j,k))
              change(i,j,k) = rate(i,j,k) * delt

              IF (change(i,j,k) > (qv(i,j,k,lvlq)-rsat(i,j,k)))         &
                  change(i,j,k)=(qv(i,j,k,lvlq)-rsat(i,j,k))

              qi(i,j,k,lvlq) = qi(i,j,k,lvlq) + change(i,j,k)
              qv(i,j,k,lvlq) = qv(i,j,k,lvlq) - change(i,j,k)
            END IF

          ELSE ! ice evaporation

!
!-----------------------------------------------------------------------
!
!    Pristine crystals.  Might make this instantaneous.
!
!-----------------------------------------------------------------------
!
            IF (qi(i,j,k,lvlq) > 0. .AND. qv(i,j,k,lvlq) < rsat(i,j,k)) THEN

              rate(i,j,k) = p2v * (rsat(i,j,k)-qv(i,j,k,lvlq))
              change(i,j,k) = rate(i,j,k) * delt

              IF (change(i,j,k) > (rsat(i,j,k)-qv(i,j,k,lvlq)))         &
                  change(i,j,k)=(rsat(i,j,k)-qv(i,j,k,lvlq))

              IF (change(i,j,k) > qi(i,j,k,lvlq)) THEN
                qv(i,j,k,lvlq) = qv(i,j,k,lvlq) + qi(i,j,k,lvlq)
                qi(i,j,k,lvlq) = 0.
              ELSE
                qv(i,j,k,lvlq) = qv(i,j,k,lvlq) + change(i,j,k)
                qi(i,j,k,lvlq) = qi(i,j,k,lvlq) - change(i,j,k)
              END IF

            END IF

!
!-----------------------------------------------------------------------
!
!    Then snow.
!
!-----------------------------------------------------------------------
!
            IF (qs(i,j,k,lvlq) > 0. .AND. qv(i,j,k,lvlq) < rsat(i,j,k)) THEN

              rate(i,j,k) = s2v * (rsat(i,j,k)-qv(i,j,k,lvlq))
              change(i,j,k) = rate(i,j,k) * delt

              IF (change(i,j,k) > (rsat(i,j,k)-qv(i,j,k,lvlq)))         &
                  change(i,j,k)=(rsat(i,j,k)-qv(i,j,k,lvlq))

              IF (change(i,j,k) > qs(i,j,k,lvlq)) THEN
                qv(i,j,k,lvlq) = qv(i,j,k,lvlq) + qs(i,j,k,lvlq)
                qs(i,j,k,lvlq) = 0.
              ELSE
                qv(i,j,k,lvlq) = qv(i,j,k,lvlq) + change(i,j,k)
                qs(i,j,k,lvlq) = qs(i,j,k,lvlq) - change(i,j,k)
              END IF
            END IF

!
!-----------------------------------------------------------------------
!
!    And finally ice.  It might be argued that graupel and hail
!    can be water-coated and thus evaporate wrt liquid
!    saturation (i.e., faster).
!
!-----------------------------------------------------------------------
!
            IF (qh(i,j,k,lvlq) > 0. .AND. qv(i,j,k,lvlq) < rsat(i,j,k)) THEN

              rate(i,j,k) = i2v * (rsat(i,j,k)-qv(i,j,k,lvlq))
              change(i,j,k) = rate(i,j,k) * delt

              IF (change(i,j,k) > (rsat(i,j,k)-qv(i,j,k,lvlq)))         &
                  change(i,j,k)=(rsat(i,j,k)-qv(i,j,k,lvlq))

              IF (change(i,j,k) > qh(i,j,k,lvlq)) THEN
                qv(i,j,k,lvlq) = qv(i,j,k,lvlq) + qh(i,j,k,lvlq)
                qh(i,j,k,lvlq) = 0.
              ELSE
                qv(i,j,k,lvlq) = qv(i,j,k,lvlq) + change(i,j,k)
                qh(i,j,k,lvlq) = qh(i,j,k,lvlq) - change(i,j,k)
              END IF
            END IF

          END IF

          x(i,j,k) = (qi(i,j,k,lvlq) + qs(i,j,k,lvlq) +                 &
                     qh(i,j,k,lvlq)) - rice(i,j,k)
          t(i,j,k,lvlq) = t(i,j,k,lvlq) + x(i,j,k)*lvi/cp
          rice(i,j,k) = qi(i,j,k,lvlq) +                                &
                        qs(i,j,k,lvlq) + qh(i,j,k,lvlq)

        END IF

      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  MELTING AND FREEZING.
!  Melting first.  Cloud ice melts before snow, but they both melt
!  immediately.  Cloud ice melts into cloud liquid.  Snow melts into
!  rain.  Ice also melts into rain, but not immediately.  In all
!  cases, the amount of melting is limited to the amount it takes to
!  cool the parcel to the freezing point.  It works out to about
!  3 g/kg per centigrade degree.  Start by calculating
!  the maximum amount of melting possible in this time step.
!
!-----------------------------------------------------------------------
!

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

        IF (t(i,j,k,lvlq) > 273.1) THEN

          maxmelt(i,j,k) = (t(i,j,k,lvlq)-273.1) * cp/lli

          IF (qi(i,j,k,lvlq) > 0.) THEN
            IF (maxmelt(i,j,k) < qi(i,j,k,lvlq)) THEN
              qc(i,j,k,lvlq) = qc(i,j,k,lvlq) + maxmelt(i,j,k)
              qi(i,j,k,lvlq) = qi(i,j,k,lvlq) - maxmelt(i,j,k)
              GO TO 501
            ELSE
              maxmelt(i,j,k) = maxmelt(i,j,k) - qi(i,j,k,lvlq)
              qc(i,j,k,lvlq) = qc(i,j,k,lvlq) + qi(i,j,k,lvlq)
              qi(i,j,k,lvlq) = 0.
            END IF
          END IF

          IF (qs(i,j,k,lvlq) > 0.) THEN
            IF (maxmelt(i,j,k) < qs(i,j,k,lvlq)) THEN
              qr(i,j,k,lvlq) = qr(i,j,k,lvlq) + maxmelt(i,j,k)
              qs(i,j,k,lvlq) = qs(i,j,k,lvlq) - maxmelt(i,j,k)
              GO TO 501
            ELSE
              maxmelt(i,j,k) = maxmelt(i,j,k) - qs(i,j,k,lvlq)
              qr(i,j,k,lvlq) = qr(i,j,k,lvlq) + qs(i,j,k,lvlq)
              qs(i,j,k,lvlq) = 0.
            END IF
          END IF

          IF (qh(i,j,k,lvlq) > 0.) THEN
            rate(i,j,k) = i2r * (t(i,j,k,lvlq)-273.1)/5.
            change(i,j,k) = rate(i,j,k)*delt
            IF (change(i,j,k) > maxmelt(i,j,k)) change(i,j,k)=maxmelt(i,j,k)
            IF (change(i,j,k) < qh(i,j,k,lvlq)) THEN
              qr(i,j,k,lvlq) = qr(i,j,k,lvlq) + change(i,j,k)
              qh(i,j,k,lvlq) = qh(i,j,k,lvlq) - change(i,j,k)
            ELSE
              qr(i,j,k,lvlq) = qr(i,j,k,lvlq) + qh(i,j,k,lvlq)
              qh(i,j,k,lvlq) = 0.
            END IF
          END IF

          501           CONTINUE

!
!-----------------------------------------------------------------------
!
!  Now freezing.  First the cloud liquid, then the rain.
!
!-----------------------------------------------------------------------
!
        ELSE
!
!-----------------------------------------------------------------------
!
!  Cloud water stays supercooled well below freezing, but how much?
!  This is just freezing because it's cold.  There is very little
!  cloud water below -25C, and almost none observed below -40C.
!  Ramp parabolically from -20 to -40.
!
!-----------------------------------------------------------------------
!
          IF (qc(i,j,k,lvlq) > 0. .AND. t(i,j,k,lvlq) < 253.) THEN
            rate(i,j,k) = c2p * ((253.-t(i,j,k,lvlq))/20.)**2
            change(i,j,k) = rate(i,j,k) * delt
            IF (change(i,j,k) > qc(i,j,k,lvlq)) THEN
              qi(i,j,k,lvlq) = qi(i,j,k,lvlq) + qc(i,j,k,lvlq)
              qc(i,j,k,lvlq) = 0.
            ELSE
              qi(i,j,k,lvlq) = qi(i,j,k,lvlq) + change(i,j,k)
              qc(i,j,k,lvlq) = qc(i,j,k,lvlq) - change(i,j,k)
            END IF
          END IF
!
!-----------------------------------------------------------------------
!
!  Rain freezing into ice; parabolic function similar to C2P.
!  Based loosely on Fig.1 from Cotton (MWR 72b).
!
!-----------------------------------------------------------------------
!
          IF (qr(i,j,k,lvlq) > 0. .AND. t(i,j,k,lvlq) < 267.) THEN
            rate(i,j,k) = r2i * ((267.-t(i,j,k,lvlq))/14.)**2
            change(i,j,k) = rate(i,j,k) * delt
            IF (change(i,j,k) > qr(i,j,k,lvlq)) THEN
              qh(i,j,k,lvlq) = qh(i,j,k,lvlq) + qr(i,j,k,lvlq)
              qr(i,j,k,lvlq) = 0.
            ELSE
              qh(i,j,k,lvlq) = qh(i,j,k,lvlq) + change(i,j,k)
              qr(i,j,k,lvlq) = qr(i,j,k,lvlq) - change(i,j,k)
            END IF
          END IF

        END IF
!
!-----------------------------------------------------------------------
!
!  Temperature after freezing or melting.
!  This should give the same answer.
!  x = (qi(i,j,k,lvlq) + qs(i,j,k,lvlq) + qh(i,j,k,lvlq)) -
!    :     rice(i,j,k)
!  t = t + x*Lli/cp
!
!-----------------------------------------------------------------------
!
        x(i,j,k) = (qc(i,j,k,lvlq) + qr(i,j,k,lvlq)) - rliq(i,j,k)
        t(i,j,k,lvlq) = t(i,j,k,lvlq) - x(i,j,k)*lli/cp

        rliq(i,j,k) = qc(i,j,k,lvlq) + qr(i,j,k,lvlq)
        rice(i,j,k) = qi(i,j,k,lvlq) + qs(i,j,k,lvlq) + qh(i,j,k,lvlq)

      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  COLLECTION.  These processes are determined by spacing between
!  particles regardless of how much gas is also in the volume.  So
!  we'll first convert to specific contents, and then later back to
!  mixing ratios.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        qc(i,j,k,lvlq) = qc(i,j,k,lvlq) * rho(i,j,k)
        qi(i,j,k,lvlq) = qi(i,j,k,lvlq) * rho(i,j,k)
        qr(i,j,k,lvlq) = qr(i,j,k,lvlq) * rho(i,j,k)
        qs(i,j,k,lvlq) = qs(i,j,k,lvlq) * rho(i,j,k)
        qh(i,j,k,lvlq) = qh(i,j,k,lvlq) * rho(i,j,k)
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Autoconversion.  As soon as you build up enough cloud matter, it
!  starts converting to rain or snow.  This "nucleates" the
!  collection process, which is nonlinear.  The nucleated amount
!  determines how long before rapid collection occurs.
!  If some precipitate is already present, the autoconv procedure
!  just makes sure there's enough.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        IF (qc(i,j,k,lvlq) > qcmin) THEN
          change(i,j,k) = qc(i,j,k,lvlq) - qcmin
          IF (change(i,j,k) > qr(i,j,k,lvlq)) THEN
            change(i,j,k) = change(i,j,k) - qr(i,j,k,lvlq)
            qr(i,j,k,lvlq) = qr(i,j,k,lvlq) + change(i,j,k)
            qc(i,j,k,lvlq) = qc(i,j,k,lvlq) - change(i,j,k)
          END IF
        END IF
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        IF (qi(i,j,k,lvlq) > qpmin) THEN
          change(i,j,k) = qi(i,j,k,lvlq) - qpmin
          IF (change(i,j,k) > qs(i,j,k,lvlq)) THEN
            change(i,j,k) = change(i,j,k) - qs(i,j,k,lvlq)
            qs(i,j,k,lvlq) = qs(i,j,k,lvlq) + change(i,j,k)
            qi(i,j,k,lvlq) = qi(i,j,k,lvlq) - change(i,j,k)
          END IF
        END IF
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  There can be a three-way competition for cloud water among the
!  rain, snow, and ice, so we distribute it into the three
!  categories.  The "need" variables are the amount of cloud water
!  that process would use up in a time step if it didn't have to
!  compete.  The collection formula for rain is very nearly the same
!  as Soong and Ogura (1973); the other functions are based on that.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        IF (qc(i,j,k,lvlq) > 0.) THEN
          rate(i,j,k) = c2r * qc(i,j,k,lvlq) * qr(i,j,k,lvlq)
          needr(i,j,k) = rate(i,j,k)*delt
          rate(i,j,k) = c2s * qc(i,j,k,lvlq) * qs(i,j,k,lvlq)
          needs(i,j,k) = rate(i,j,k)*delt
          rate(i,j,k) = c2i * qc(i,j,k,lvlq) * qh(i,j,k,lvlq)
          needi(i,j,k) = rate(i,j,k)*delt
          tem15(i,j,k) = needr(i,j,k) + needs(i,j,k) + needi(i,j,k)
          IF (tem15(i,j,k) > qc(i,j,k,lvlq)) THEN
            needr(i,j,k) = needr(i,j,k) * qc(i,j,k,lvlq) / tem15(i,j,k)
            needs(i,j,k) = needs(i,j,k) * qc(i,j,k,lvlq) / tem15(i,j,k)
            needi(i,j,k) = needi(i,j,k) * qc(i,j,k,lvlq) / tem15(i,j,k)
          END IF

!
!-----------------------------------------------------------------------
!
!  The riming process nucleates a little cloud ice.  Until a better
!  number comes along we'll say 1% of the collected liquid, both for
!  snow and graupel.
!
!-----------------------------------------------------------------------
!
          qi(i,j,k,lvlq) = qi(i,j,k,lvlq) + .01*needs(i,j,k)
          needs(i,j,k) = needs(i,j,k) - .01*needs(i,j,k)
          qi(i,j,k,lvlq) = qi(i,j,k,lvlq) + .01*needi(i,j,k)
          needi(i,j,k) = needi(i,j,k) - .01*needi(i,j,k)
!
!-----------------------------------------------------------------------
!
!  For simplicity (i.e., until a better way to do it comes along),
!  we assume that the result of snow riming is to convert the
!  cloud water, but not the snow, to the graupel category.
!
!-----------------------------------------------------------------------
!
          qr(i,j,k,lvlq) = qr(i,j,k,lvlq) + needr(i,j,k)
          qh(i,j,k,lvlq) = qh(i,j,k,lvlq) + needi(i,j,k) + needs(i,j,k)
          qc(i,j,k,lvlq) = qc(i,j,k,lvlq) - needr(i,j,k) -              &
                           needs(i,j,k) - needi(i,j,k)

        END IF

      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Unlike the collection of cloud water, which can be complete, we
!  don't want to zero out the pristine crystals.  For one thing,
!  some are so tiny they won't get collected, but also, we want to
!  be able to produce more condensate if it still supersaturated wrt
!  ice, which won't happen if there's zero cloud ice.  We'll leave
!  behind the equivalent of 100 per liter.  At 1E-11 kg per crystal
!  and 1000 liters per m**3, that would be 1E-6 kg/m**3 (1 mg/m**3).
!  The temperature-dependent efficiency follows Lin et al (JCAM 83).
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        IF (qi(i,j,k,lvlq) > 0. .AND. qs(i,j,k,lvlq) > 0.) THEN
          eff = 1. - (273.1-t(i,j,k,lvlq))/50.
          rate(i,j,k) = p2s * eff * qi(i,j,k,lvlq) * qs(i,j,k,lvlq)
          change(i,j,k) = rate(i,j,k) * delt
          IF (change(i,j,k) > qi(i,j,k,lvlq)) THEN
            qs(i,j,k,lvlq) = qs(i,j,k,lvlq) + qi(i,j,k,lvlq) - 1E-6
            qi(i,j,k,lvlq) = 1E-6
          ELSE
            qs(i,j,k,lvlq) = qs(i,j,k,lvlq) + change(i,j,k)
            qi(i,j,k,lvlq) = qi(i,j,k,lvlq) - change(i,j,k)
          END IF
        END IF
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!   Convert back to mixing ratios.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        qc(i,j,k,lvlq) = qc(i,j,k,lvlq) / rho(i,j,k)
        qi(i,j,k,lvlq) = qi(i,j,k,lvlq) / rho(i,j,k)
        qr(i,j,k,lvlq) = qr(i,j,k,lvlq) / rho(i,j,k)
        qs(i,j,k,lvlq) = qs(i,j,k,lvlq) / rho(i,j,k)
        qh(i,j,k,lvlq) = qh(i,j,k,lvlq) / rho(i,j,k)

!
!-----------------------------------------------------------------------
!
!  Calculate temperature after phase change(i,j,k)s resulting from
!  collection.
!
!-----------------------------------------------------------------------
!
        x(i,j,k) = (qc(i,j,k,lvlq) + qr(i,j,k,lvlq)) - rliq(i,j,k)
        t(i,j,k,lvlq) = t(i,j,k,lvlq) - x(i,j,k)*lli/cp
        rliq(i,j,k) = qc(i,j,k,lvlq) + qr(i,j,k,lvlq)
        rice(i,j,k) = qi(i,j,k,lvlq) + qs(i,j,k,lvlq) + qh(i,j,k,lvlq)
      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Now comes the final saturation adjustment, if necessary.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        rsat(i,j,k) = rhsat*f_qvsatl( p(i,j,k,lvlq), t(i,j,k,lvlq) )

        IF (.NOT.(qc(i,j,k,lvlq) == 0. .AND. qv(i,j,k,lvlq) < rsat(i,j,k))) THEN

! Call SatAdjL (qv(i,j,k,lvlq), p, t, rsat(i,j,k), t1)

          t1(i,j,k) = t(i,j,k,lvlq)
          dr = qv(i,j,k,lvlq) - rsat(i,j,k)
          err_num = t1(i,j,k) - t(i,j,k,lvlq) - lvl/cp * dr
          err_den = 1. + lvl/cp * rsat(i,j,k) * f_desdtl( t1(i,j,k) )

          t1(i,j,k) = t1(i,j,k) - err_num/err_den
          rsat(i,j,k) = rhsat*f_qvsatl( p(i,j,k,lvlq), t1(i,j,k) )
          dr = qv(i,j,k,lvlq) - rsat(i,j,k)
          err_num = t1(i,j,k) - t(i,j,k,lvlq) - lvl/cp * dr
          err_den = 1. + lvl/cp * rsat(i,j,k) * f_desdtl( t1(i,j,k) )

          t1(i,j,k) = t1(i,j,k) - err_num/err_den
          rsat(i,j,k) = rhsat*f_qvsatl( p(i,j,k,lvlq), t1(i,j,k) )
          dr = qv(i,j,k,lvlq) - rsat(i,j,k)
          err_num = t1(i,j,k) - t(i,j,k,lvlq) - lvl/cp * dr
          err_den = 1. + lvl/cp * rsat(i,j,k) * f_desdtl( t1(i,j,k) )

          t1(i,j,k) = t1(i,j,k) - err_num/err_den
          rsat(i,j,k) = rhsat*f_qvsatl( p(i,j,k,lvlq), t1(i,j,k) )
          dr = qv(i,j,k,lvlq) - rsat(i,j,k)

          qc(i,j,k,lvlq) = qc(i,j,k,lvlq) + (qv(i,j,k,lvlq)-rsat(i,j,k))

          IF (qc(i,j,k,lvlq) > 0.) THEN
            qv(i,j,k,lvlq) = rsat(i,j,k)
          ELSE
            qv(i,j,k,lvlq) = rsat(i,j,k) + qc(i,j,k,lvlq)
            qc(i,j,k,lvlq) = 0.
          END IF

          x(i,j,k) = qc(i,j,k,lvlq) + qr(i,j,k,lvlq) - rliq(i,j,k)
          t(i,j,k,lvlq) = t(i,j,k,lvlq) + x(i,j,k)*lvl/cp
          rliq(i,j,k) = qc(i,j,k,lvlq) + qr(i,j,k,lvlq)

        END IF

      END DO
    END DO
  END DO

!
!-----------------------------------------------------------------------
!
!  Conservation checks.
!
!-----------------------------------------------------------------------
!

! rtot1(i,j,k) = qv(i,j,k,lvlq) + rliq(i,j,k) + rice(i,j,k)
! If (abs(rtot0-rtot1(i,j,k)).gt..000001) then
!  print*, 'rtot check', rtot0, rtot1(i,j,k)
!  print*, qv(i,j,k,lvlq), qc(i,j,k,lvlq), qr(i,j,k,lvlq), qi(i,j,k,lvlq), qs(i,j,k,lvlq), qh(i,j,k,lvlq)
! End if

! til1 = t - (Lvl*rliq(i,j,k)+Lvi*rice(i,j,k))/cp
! If (abs(til1-til0).gt..0001) then
!  pqh(i,j,k,lvlq)nt*, 'til check', til0, til1
! End if

!
!-----------------------------------------------------------------------
!
!  These terminal velocity formulations are similar in form to Ogura
!  and Takahashi (1971).  The curves for rain and snow were tweaked
!  until they matched the curves on page 241 of Pielke's book (they
!  are very, very close).  The curve for ice is based on the curve
!  for rain; the only difference is in the exponent.  The effect is
!  that small values of ice, presumed to be heavily-rimed snow, fall
!  slower than rain of the same concentration, but higher values,
!  presumed to be big graupel or hail, fall faster than rain.  The
!  transition is at 1 g/kg.  The terminal velocity for pristine
!  crystals does not have a dependency on mixing ratio like the
!  others, because the others incorporate(i,j,k) the assumption that
!  higher mixing ratios imply bigger particles, which fall faster,
!  but that's not the case for pristine crystals.
!
!-----------------------------------------------------------------------
!
  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        t(i,j,k,lvlq)=(t(i,j,k,lvlq)/ppi(i,j,k))-ptbar(i,j,k)
        p(i,j,k,lvlq)=p(i,j,k,lvlq)-pbar(i,j,k)
        qv(i,j,k,lvlq)=MAX(0.0,qv(i,j,k,lvlq))
        qc(i,j,k,lvlq)=MAX(0.0,qc(i,j,k,lvlq))
        qr(i,j,k,lvlq)=MAX(0.0,qr(i,j,k,lvlq))
        qi(i,j,k,lvlq)=MAX(0.0,qi(i,j,k,lvlq))
        qs(i,j,k,lvlq)=MAX(0.0,qs(i,j,k,lvlq))
        qh(i,j,k,lvlq)=MAX(0.0,qh(i,j,k,lvlq))
      END DO
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tvq(i,j,k)=0.0
        IF (qi(i,j,k,lvlq) > 0.) tvq(i,j,k) = 0.5 * SQRT(1./rho(i,j,k))
      END DO
    END DO
  END DO

  CALL qfallout(nx,ny,nz,dtbig1,qi,rho,j3,j3inv,tem15(1,1,1),           &
                tvq,tem16)

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tvq(i,j,k)=0.0
        IF (qr(i,j,k,lvlq) > 0.) tvq(i,j,k) =                           &
            5.5 * (1000.*qr(i,j,k,lvlq))**.125 * SQRT(1./rho(i,j,k))
      END DO
    END DO
  END DO

  CALL qfallout(nx,ny,nz,dtbig1,qr,rho,j3,j3inv,tem15(1,1,2),           &
                tvq,tem16)

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tvq(i,j,k)=0.0
        IF (qs(i,j,k,lvlq) > 0.) tvq(i,j,k) =                           &
            2.0 * (1000.*qs(i,j,k,lvlq))**.100 * SQRT(1./rho(i,j,k))
      END DO
    END DO
  END DO

  CALL qfallout(nx,ny,nz,dtbig1,qs,rho,j3,j3inv,tem15(1,1,3),           &
                tvq,tem16)

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        tvq(i,j,k)=0.0
        IF (qh(i,j,k,lvlq) > 0.) tvq(i,j,k) = 5.5 *                     &
            (1000.*qh(i,j,k,lvlq))**.333 * SQRT(1./rho(i,j,k))
      END DO
    END DO
  END DO

  CALL qfallout(nx,ny,nz,dtbig1,qh,rho,j3,j3inv,tem15(1,1,4),           &
                tvq,tem16)

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

  DO j=1,ny-1
    DO i=1,nx-1
      tema=tem15(i,j,2)+tem15(i,j,3)+tem15(i,j,4)
      raing(i,j)=raing(i,j)+tema
      prcrate(i,j) = tema*tem
    END DO
  END DO

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        qr(i,j,k,lvlq)=MAX(0.0,qr(i,j,k,lvlq))
        qi(i,j,k,lvlq)=MAX(0.0,qi(i,j,k,lvlq))
        qs(i,j,k,lvlq)=MAX(0.0,qs(i,j,k,lvlq))
        qh(i,j,k,lvlq)=MAX(0.0,qh(i,j,k,lvlq))
      END DO
    END DO
  END DO

  RETURN
END SUBROUTINE microph_nem