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


SUBROUTINE microph_ice(nx,ny,nz, dtbig1,                                & 1,5
           ptprt,pprt,qv,qc,qr,qi,qs,qh,raing,prcrate,                  &
           rhostr,pbar,ptbar,qvbar,ppi,j3,j3inv,                        &
           rhobar,cgsrhobar,cgspres,tem1,tem2,tem3,tem4,tem5,           &
           tem6,tem7,tem8,tem9,tem10,tem11,tem12)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate and apply the microphysical contributions to the water,
!  ice and temperature fields, using an ice microphysics parameterization
!  scheme.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/1991.
!
!  MODIFICATION HISTORY:
!
!  5/02/92 (M. Xue)
!  Added full documentation.
!
!  5/28/92 (K. Brewster)
!  Further facelift.
!
!  9/10/92 (M. Xue)
!  Negative water contents are not zeroed out in this version.
!
!  09/15/94 (M. Xue)
!  Modified to adapt Tao's ice code.
!
!  10/20/94 (Yuhe Liu)
!  Fixed a bug in the argument list of calling subroutine ICECVT.
!
!  03/05/97 (Fanyou Kong -- CMRP)
!  Modify the code to apply all three time levels in calculate
!  production terms in Tao scheme (icecvt)
!
!  07/10/97 (Fanyou Kong - CMRP)
!  Include MPDCD advection option (sadvopt = 5)
!
!  11/18/98 (Keith Brewster)
!  Changed pibar to ppi (full pi).
!
!  08/27/02 (Dan Weber)
!  Added option for using Ferrier (1994) fall velocity coefficients 
!  and replaces the Lin constants when fallopt=2.
!
!-----------------------------------------------------------------------
!
!  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   The large time step size for this call.
!
!    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)
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!    exner
!
!  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)
!    cgsrhobarBase state density in cgs unit.
!    cgspres  Base state pressure in cgs unit.
!    tem1     Temporary work array.
!    tem2     Temporary work array.
!    tem3     Temporary work array.
!    tem4     Temporary work array.
!    tem5     Temporary work array.
!    tem6     Temporary work array.
!    tem7     Temporary work array.
!    tem8     Temporary work array.
!    tem9     Temporary work array.
!    tem10    Temporary work array.
!    tem11    Temporary work array.
!    tem12    Temporary work array.
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INCLUDE 'timelvls.inc'

  INTEGER :: irsh         ! Flag identifying hydrometer fields
                          ! 0 = rain; 1 = snow; 2 = hail
!
  INTEGER :: nx,ny,nz          ! Number of grid points in 3 directions
  REAL :: dtbig1               ! The large time step size for this call.
!
  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 :: j3inv (nx,ny,nz)     ! Coordinate transformation Jacobian defined as
                               ! d( zp )/d( z ).
  REAL :: qvbar (nx,ny,nz)     ! Base state water vapor specific humidity
                               ! (kg/kg)
  REAL :: ppi   (nx,ny,nz)     ! Exner function.
!
!-----------------------------------------------------------------------
!
!  Temporary arrays
!
!-----------------------------------------------------------------------
!
  REAL :: rhobar(nx,ny,nz)     ! Temporary work array
  REAL :: cgsrhobar(nx,ny,nz)  ! Temporary work array
  REAL :: cgspres(nx,ny,nz)    ! Temporary work 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
  REAL :: tem4  (nx,ny,nz)     ! Temporary work array
  REAL :: tem5  (nx,ny,nz)     ! Temporary work array
  REAL :: tem6  (nx,ny,nz)     ! Temporary work array
  REAL :: tem7  (nx,ny,nz)     ! Temporary work array
  REAL :: tem8  (nx,ny,nz)     ! Temporary work array
  REAL :: tem9  (nx,ny,nz)     ! Temporary work array
  REAL :: tem10 (nx,ny,nz)     ! Temporary work array
  REAL :: tem11 (nx,ny,nz)     ! Temporary work array
  REAL :: tem12 (nx,ny,nz)     ! Temporary work array
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k,lvlq
!  real tbar
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
!
!-----------------------------------------------------------------------
!
!  Following constants are defined in subroutine STCSTICE
!
!  (Suggest to move the following COMMON block to a inlcude file.)
!
!-----------------------------------------------------------------------
!
  REAL :: tnw,tns,tng,roqr,roqs,roqg
  REAL :: c76,c358,c172,c409,c218,c580,c610,c149,c879,c141
  REAL :: zrc,zgc,zsc,vrc,vgc,vsc
  REAL :: ag,bg,as,bs,aww,bww,bgh,bgq,bsh,bsq,bwh,bwq
  REAL :: alv,alf,als,t0,t00,avc,afc,asc,rn1,bnd1,rn2,bnd2,             &
       rn3,rn4,rn5,rn6,rn7,rn8,rn9,rn10,rn101,rn10a,rn11,rn11a,         &
       rn12,rn12a(31),rn12b(31),rn13(31),rn14,rn15,rn15a,rn16,rn17,     &
       rn17a,rn17b,rn17c,rn18,rn18a,rn19,rn19a,rn19b,rn20,rn20a,        &
       rn20b,bnd3,rn21,rn22,rn23,rn23a,rn23b,rn25,rn25a(31),rn30a,      &
       rn30b,rn30c,rn31,beta,rn32

  COMMON/size/ tnw,tns,tng,roqr,roqs,roqg
  COMMON/cont/ c76,c358,c172,c409,c218,c580,c610,c149,c879,c141
  COMMON/bterv/ zrc,zgc,zsc,vrc,vgc,vsc
  COMMON/b3cs/ ag,bg,as,bs,aww,bww,bgh,bgq,bsh,bsq,bwh,bwq
  COMMON/bsnw/ alv,alf,als,t0,t00,avc,afc,asc,rn1,bnd1,rn2,bnd2,        &
       rn3,rn4,rn5,rn6,rn7,rn8,rn9,rn10,rn101,rn10a,rn11,rn11a,         &
       rn12,rn12a,rn12b,rn13,rn14,rn15,rn15a,rn16,rn17,                 &
       rn17a,rn17b,rn17c,rn18,rn18a,rn19,rn19a,rn19b,rn20,rn20a,        &
       rn20b,bnd3,rn21,rn22,rn23,rn23a,rn23b,rn25,rn25a,rn30a,          &
       rn30b,rn30c,rn31,beta,rn32

  INTEGER :: constset
  SAVE constset
  DATA constset /0/
  REAL :: denwater,tem,tema, deltat
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  lvlq = tfuture
!
!-----------------------------------------------------------------------
!
!  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.
!
!-----------------------------------------------------------------------
!
  IF (constset == 0) THEN

    CALL stcstice
    constset = 1

  END IF

  DO k = 1,nz-1
    DO j = 1,ny-1
      DO i = 1,nx-1
        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))
        qi(i,j,k,tfuture) = MAX(0.0,qi(i,j,k,tfuture))
        qs(i,j,k,tfuture) = MAX(0.0,qs(i,j,k,tfuture))
        qh(i,j,k,tfuture) = MAX(0.0,qh(i,j,k,tfuture))
      END DO
    END DO
  END DO
!
  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)
        cgsrhobar(i,j,k) = rhobar(i,j,k) * 0.001
        cgspres  (i,j,k) = (pprt(i,j,k,tfuture)+pbar(i,j,k))*10.0
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Call ice parameterization routine that does conversions among
!  water and ice categories.
!
!-----------------------------------------------------------------------
!
!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

!call test_dump2(qr(1,1,1,tpast),"XXXBicecvt_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXBicecvt_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXBicecvt_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_check_in("Bicecvt")
  CALL icecvt(nx, ny, nz, deltat,                                       &
              ptprt,qv,qc,qr,qi,qs,qh,                                  &
              cgsrhobar, ptbar,qvbar,cgspres, ppi,                      &
              tem1(1,1,1), tem1(1,1,2), tem1(1,1,3), tem1(1,1,4),       &
              tem2(1,1,1), tem2(1,1,2), tem2(1,1,3), tem2(1,1,4),       &
              tem3(1,1,1), tem3(1,1,2), tem3(1,1,3), tem3(1,1,4),       &
              tem4(1,1,1), tem4(1,1,2), tem4(1,1,3), tem4(1,1,4),       &
              tem5(1,1,1), tem5(1,1,2), tem5(1,1,3), tem5(1,1,4),       &
              tem6(1,1,1), tem6(1,1,2), tem6(1,1,3), tem6(1,1,4),       &
              tem7(1,1,1), tem7(1,1,2), tem7(1,1,3), tem7(1,1,4),       &
              tem8(1,1,1), tem8(1,1,2), tem8(1,1,3), tem8(1,1,4),       &
              tem9(1,1,1), tem9(1,1,2), tem9(1,1,3), tem9(1,1,4),       &
              tem10(1,1,1),tem10(1,1,2),tem10(1,1,3),tem10(1,1,4),      &
              tem11(1,1,1),tem11(1,1,2),tem11(1,1,3),tem11(1,1,4),      &
              tem12(1,1,1),tem12(1,1,2))
!call test_dump2(qr(1,1,1,tpast),"XXXAicecvt_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXAicecvt_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXAicecvt_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qv(1,1,1,tpast),"XXXAicecvt_qv1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qv(1,1,1,tpresent),"XXXAicecvt_qv2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qv(1,1,1,tfuture),"XXXAicecvt_qv3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qc(1,1,1,tpast),"XXXAicecvt_qc1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qc(1,1,1,tpresent),"XXXAicecvt_qc2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qc(1,1,1,tfuture),"XXXAicecvt_qc3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qi(1,1,1,tpast),"XXXAicecvt_qi1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qi(1,1,1,tpresent),"XXXAicecvt_qi2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qi(1,1,1,tfuture),"XXXAicecvt_qi3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qs(1,1,1,tpast),"XXXAicecvt_qs1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qs(1,1,1,tpresent),"XXXAicecvt_qs2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qs(1,1,1,tfuture),"XXXAicecvt_qs3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qh(1,1,1,tpast),"XXXAicecvt_qh1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qh(1,1,1,tpresent),"XXXAicecvt_qh2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qh(1,1,1,tfuture),"XXXAicecvt_qh3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(tem1,"XXXAicecvt_tem1",nx,ny,4,1,nx-1,2,ny-2,1,4)
!call test_dump2(tem2,"XXXAicecvt_tem2",nx,ny,4,1,nx-1,2,ny-2,1,4)
!call test_dump2(tem3,"XXXAicecvt_tem3",nx,ny,4,1,nx-1,2,ny-2,1,4)
!call test_dump2(tem4,"XXXAicecvt_tem4",nx,ny,4,1,nx-1,2,ny-2,1,4)
!call test_dump2(tem5,"XXXAicecvt_tem5",nx,ny,4,1,nx-1,2,ny-2,1,4)
!call test_dump2(tem6,"XXXAicecvt_tem6",nx,ny,4,1,nx-1,2,ny-2,1,4)
!call test_dump2(tem7,"XXXAicecvt_tem7",nx,ny,4,1,nx-1,2,ny-2,1,4)
!call test_dump2(tem8,"XXXAicecvt_tem8",nx,ny,4,1,nx-1,2,ny-2,1,4)
!call test_dump2(tem9,"XXXAicecvt_tem9",nx,ny,4,1,nx-1,2,ny-2,1,4)
!call test_dump2(tem10,"XXXAicecvt_tem10",nx,ny,4,1,nx-1,2,ny-2,1,4)
!call test_dump2(tem11,"XXXAicecvt_tem11",nx,ny,4,1,nx-1,2,ny-2,1,4)
!call test_dump2(tem12,"XXXAicecvt_tem12",nx,ny,4,1,nx-1,2,ny-2,1,4)
!call test_dump2(ptprt,"XXXAicecvt_ptprt",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(cgsrhobar,"XXXAicecvt_cgsrhobar",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(ptbar,"XXXAicecvt_ptbar",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qvbar,"XXXAicecvt_qvbar",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(cgspres,"XXXAicecvt_cgspres",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(ppi,"XXXAicecvt_ppi",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_check_in("Aicecvt")
!
!-----------------------------------------------------------------------
!
!  Hydrometeor sedimentation
!
!-----------------------------------------------------------------------
!
  irsh = 0  ! for rainwater
!call test_check_in("Bqhfall_qr")
  CALL qhfall (irsh,dtbig1,nx,ny,nz, qr, rhobar,j3,j3inv,               &
               tem4(1,1,1), tem1,tem2)
!call test_dump2(qr(1,1,1,tpast),"XXXAqhfall_qr1",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tpresent),"XXXAqhfall_qr2",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_dump2(qr(1,1,1,tfuture),"XXXAqhfall_qr3",nx,ny,nz,1,nx-1,2,ny-2,2,nz-2)
!call test_check_in("Aqhfall_qr")

  irsh = 1  ! for snow
!call test_check_in("Bqhfall_qs")
  CALL qhfall (irsh,dtbig1,nx,ny,nz, qs, rhobar,j3,j3inv,               &
               tem4(1,1,2), tem1,tem2)
!call test_check_in("Aqhfall_qs")

  irsh = 2  ! for hail
!call test_check_in("Bqhfall_qh")
  CALL qhfall (irsh,dtbig1,nx,ny,nz, qh, rhobar,j3,j3inv,               &
               tem4(1,1,3), tem1,tem2)
!call test_check_in("Aqhfall_qh")

  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
      tema = tem4(i,j,1)+tem4(i,j,2)+tem4(i,j,3)
      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
        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))
        qi(i,j,k,tfuture) = MAX(0.0,qi(i,j,k,tfuture))
        qs(i,j,k,tfuture) = MAX(0.0,qs(i,j,k,tfuture))
        qh(i,j,k,tfuture) = MAX(0.0,qh(i,j,k,tfuture))
      END DO
    END DO
  END DO
!
  RETURN
END SUBROUTINE microph_ice

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


SUBROUTINE qhfall(irsh,dtbig1,nx,ny,nz,q,rhobar,j3,j3inv,               & 3,2
           draing, vtr3d,tem1)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate the fall-out of hydrometers.
!  This subroutine is called by the ice microphysics package.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Ming Xue
!  10/10/1991. Written based on QRFALL used by warmrain microphysics.
!
!  MODIFICATION HISTORY:
!
!  10/20/1996 (M. Xue)
!  This routine now calls a standard routine QFALLOUT.
!  Precipitation rate array prcrate is now correctly calculated
!  for split-step integration.

!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!    irsh     Flag identifying hydrometer fields
!    dtbig1   The large time step size for this call.
!    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
!
!    q        Hydrometer mixing ratio at all time levels (kg/kg)
!
!    rhobar   Base state air density (kg/m**3)
!    j3       Coordinate transformation Jacobian  d(zp)/dz
!    j3inv    1/j3.
!
!  OUTPUT:
!
!    q        Hydrometer mixing ratio at time tfuture (kg/kg)
!    draing   Grid-scale rainfall (mm)
!
!  WORK ARRAYS:
!
!    vtr3d    Work array
!    tem1     Work array
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  Variable Declarations.
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INCLUDE 'timelvls.inc'
!
  INTEGER :: irsh         ! Flag identifying hydrometer fields
                          ! 0 = rain, 1 = snow, 2 = hail
  REAL :: dtbig1          ! The large time step size for this call.
  INTEGER :: nx,ny,nz     ! Number of grid points in 3 directions

  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)     ! 1/j3.
  REAL :: draing (nx,ny)       ! Grid-scale rainfall (mm)
!
!-----------------------------------------------------------------------
!
!  Temporary arrays
!
!-----------------------------------------------------------------------
!
  REAL :: vtr3d(nx,ny,nz)
  REAL :: tem1 (nx,ny,nz)
!
!-----------------------------------------------------------------------
!
!  Misc. local variables
!
!-----------------------------------------------------------------------
!
  INTEGER :: i,j,k
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'globcst.inc'
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  DO k = 1, nz-1
    DO j = 1, ny-1
      DO i = 1, nx-1
        tem1(i,j,k) = rhobar(i,j,k) * 0.001
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate terminal falling speed for q, and store in tem1
!
!-----------------------------------------------------------------------
!
  CALL terv(nx,ny,nz,irsh, tem1, q(1,1,1,tpresent), vtr3d)

!-----------------------------------------------------------------------
!
!  vtr3d is defined at the w point and is in cgs unit.
!
!-----------------------------------------------------------------------

  DO k=1,nz-1
    DO j=1,ny-1
      DO i=1,nx-1
        vtr3d(i,j,k)=vtr3d(i,j,k)*0.01
      END DO
    END DO
  END DO
!
!-----------------------------------------------------------------------
!
!  Calculate the rainwater fallout term and the precipitation rate.
!
!-----------------------------------------------------------------------
!

!call test_check_in("Bqfall_in_qhfall_micro_ice3d")
  CALL qfallout(nx,ny,nz,dtbig1,q, rhobar,j3,j3inv, draing,vtr3d,       &
                tem1)
!call test_check_in("Aqfall_in_qhfall_micro_ice3d")

  RETURN
END SUBROUTINE qhfall


!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE ICECVT                    ######
!######                                                      ######
!######                     Developed by                     ######
!######                                                      ######
!######    Goddard Cumulus Ensemble Modeling Group, NASA     ######
!######                                                      ######
!######     Center for Analysis and Prediction of Storms     ######
!######               University of Oklahoma                 ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE icecvt(nx,ny,nz, d2t,                                        & 1
           ptprt,qv,qc,qr,qi,qs,qh,                                     &
           rhobar,ptbar,qvbar,pres,ppi,                                 &
           wgacr, scv, tca, dwv, zr, vr, zs, vs,                        &
           zg, vg, psaut, psaci,psacw, qsacw,praci,piacr,               &
           praut, pracw, psfw, psfi, dgacs, dgacw, dgaci,               &
           dgacr, pgacs, wgacs, qgacw,wgaci,qgacr,pgaut,                &
           pracs, psacr, qsacr, pgfr, tair, tairc,                      &
           tair3, tairc3,                                               &
           pr, pg, ps, dlt2, dlt3, rtair, tem2d1, tem2d2)
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Calculate and apply the microphysical contributions to the water/ice
!  and temperature fields. It uses Y.L. Lin's ice microphysics
!  parameterization scheme and the code is adopted from W. G. Tao's
!  ice package.
!
!-----------------------------------------------------------------------
!
!  REFERENCES:
!
!  Lin et. al. (1983)
!  TAo and Simpson (1989, JAS)
!
!-----------------------------------------------------------------------
!
!  AUTHOR: W.G. Tao, Goddard Cumulus Ensemble Modeling Group, NASA
!
!  MODIFICATION HISTORY:
!
!  12/29/93 (A. Sathye)
!  Adopted the original code and converted 2-D into 3-D.
!
!  01/05/94 (Y. Liu)
!  Formatting the code in accordance with the ARPS coding format.
!  Continue the conversion of 2-D to 3-D.
!
!  8/9/94 (MX)
!  A statement using unset variable r4f corrected.
!
!  10/14/94 (V. Wong, A. Sathye, Y. Liu, & X. Song)
!  Walked through the code and fixed a few bugs.
!
!  03/01/96 (Vince Wong)
!  Modified to prevent qs divided by zero.
!
!  10/23/1996 (M. Xue)
!  Altered codes in loops 150 and 200 to avoid trunction errors that
!  were causing the program to bomb.
!
!  1/23/1997 (J. Zong, M. Xue, V. Wong)
!  Corrected bugs in loop 150 and 200 introduced in modification
!  on 10/23/1996.
!
!  2/11/1997 (J. Zong)
!  Set a positive lower limit to qi when calculating deposition of
!  qi to avoid arithmetic exception.
!
!  2/26/1997 (J. Zong, M. Xue and Yuhe Liu)
!  Fractional power calculations are replaced by lookup table functions.
!
!  3/05/1997 (Fanyou Kong)
!  Modify the time levels used to calculate microphysical source and
!  sink terms, to make the scheme more physically rational;
!
!  07/10/97 (Fanyou Kong - CMRP)
!  Change sqrt(sqrt(zr**11)) to zr**2.75 to allow the code
!  executable on DEC Alpha system
!
!  5/19/1998 (M. Xue)
!  Changed reference density at surface from rhobar(1,1,2) to rhobar0.
!
!  11/18/98 (Keith Brewster)
!  Changed pibar to ppi (full pi).
!
!  July 2001 
!  Correction to rainwater evaporation term.  See ARPSSUPPORT e-mails
!  by Vince Wong and Eric Kemp.
!
!  Fall 2002
!  Added fallopt option to specify which hydrometeor fall velocity
!  formulations -- original Lin-Tao or Ferrier -- should be used.
!  ***NOTE*** Additional work is required for some of the transfer
!  terms when the Ferrier fall speeds are used.
!
!  13 May 2003 (Eric Kemp)
!  Bug fix to Praut term.  Added rhobar to denominator following
!  Orville and Kopp, JAS, 1977, equation 17.  The exclusion of the
!  rhobar in the Lin et al., JAM, 1983 manuscript is a typo.  This
!  error was detected by Jianzhong Wang.
!
!-----------------------------------------------------------------------
!
!  ORIGINAL COMMENTS:
!
!  lin et al (83) ice phase microphysical processes
!  modified and coded by goddard cumulus ensemble modeling group
!  tao, simpson and mccumber's saturation technique (mwr, 1989)
!  tao and simpson (jas, 1989; tao, 1993)
!
!  d2t - leapfrog time step
!  dpt - deviation potential temperature field
!  dqv - deviation water vapor field
!  qcl - cloud water field
!  qrn - rain field
!  qci - cloud ice field
!  qcs - snow field
!  qcg - hail field
!  rho - air density
!  ta1 - base air potential temperature at the level
!  qa1 - base water vapor at the level
!  p0 - air pressure
!  ppi - exner function
!
!  nx: dimension in x-direction
!  nz: dimension in z-direction
!  note: physical domain extends from k=2 to k=nz-1 and i=2 to i=nx-1
!
!  c.g.s units
!
!-----------------------------------------------------------------------
!
  IMPLICIT NONE
!
  INCLUDE 'timelvls.inc'
!
  INTEGER :: nx, ny, nz

  REAL :: ptprt (nx,ny,nz,nt)  ! Perturbation potential temperature (K)
  REAL :: qv    (nx,ny,nz,nt)  ! Water vapor specific humidity (g/g)
  REAL :: qc    (nx,ny,nz,nt)  ! Cloud water mixing ratio (g/g)
  REAL :: qr    (nx,ny,nz,nt)  ! Rainwater mixing ratio (g/g)
  REAL :: qi    (nx,ny,nz,nt)  ! Cloud ice mixing ratio (g/g)
  REAL :: qs    (nx,ny,nz,nt)  ! Snow mixing ratio (g/g)
  REAL :: qh    (nx,ny,nz,nt)  ! Hail mixing ratio (g/g)

  REAL :: rhobar(nx,ny,nz)  ! Base state air density (g/cm**3)
  REAL :: pres  (nx,ny,nz)  ! Pressure ((g cm/s**2)/cm**2)
  REAL :: ptbar (nx,ny,nz)  ! Base state potential temperature (K)
  REAL :: qvbar (nx,ny,nz)  ! Base state air density (g/g)
  REAL :: ppi   (nx,ny,nz)  ! Exner function
!
!-----------------------------------------------------------------------
!
!  Include files:
!
!-----------------------------------------------------------------------
!
  INCLUDE 'phycst.inc'
  INCLUDE 'globcst.inc'
!
!-----------------------------------------------------------------------
!
!  Local work arrays for production terms etc.
!
!-----------------------------------------------------------------------
!

  REAL :: wgacr(nx,ny)
  REAL :: psaut(nx,ny)
  REAL :: psaci(nx,ny)
  REAL :: psacw(nx,ny)
  REAL :: qsacw(nx,ny)
  REAL :: praci(nx,ny)
  REAL :: piacr(nx,ny)
  REAL :: praut(nx,ny)
  REAL :: pracw(nx,ny)
  REAL :: psfw (nx,ny)
  REAL :: psfi (nx,ny)
  REAL :: dgacs(nx,ny)
  REAL :: dgacw(nx,ny)
  REAL :: dgaci(nx,ny)
  REAL :: dgacr(nx,ny)
  REAL :: pgacs(nx,ny)
  REAL :: wgacs(nx,ny)
  REAL :: qgacw(nx,ny)
  REAL :: wgaci(nx,ny)
  REAL :: qgacr(nx,ny)
  REAL :: pgaut(nx,ny)
  REAL :: pracs(nx,ny)
  REAL :: psacr(nx,ny)
  REAL :: qsacr(nx,ny)
  REAL :: pgfr (nx,ny)

  REAL :: tair (nx,ny)
  REAL :: tairc(nx,ny)
  REAL :: tair3 (nx,ny)
  REAL :: tairc3(nx,ny)
  REAL :: pr   (nx,ny)
  REAL :: pg   (nx,ny)
  REAL :: ps   (nx,ny)
  REAL :: dlt2 (nx,ny)
  REAL :: dlt3 (nx,ny)
  REAL :: rtair(nx,ny)

  REAL :: scv  (nx,ny)
  REAL :: tca  (nx,ny)
  REAL :: dwv  (nx,ny)
  REAL :: zr   (nx,ny)
  REAL :: vr   (nx,ny)
  REAL :: zs   (nx,ny)
  REAL :: vs   (nx,ny)
  REAL :: zg   (nx,ny)
  REAL :: vg   (nx,ny)

  REAL :: tem2d1(nx,ny)
  REAL :: tem2d2(nx,ny)
!
!-----------------------------------------------------------------------
!
!  Following constants are defined and set in subroutine STCSTICE
!
!  (Suggest to move the following COMMON block to a inlcude file.)
!
!-----------------------------------------------------------------------
!
  REAL :: tnw,tns,tng,roqr,roqs,roqg
  REAL :: c76,c358,c172,c409,c218,c580,c610,c149,c879,c141
  REAL :: zrc,zgc,zsc,vrc,vgc,vsc
  REAL :: ag,bg,as,bs,aww,bww,bgh,bgq,bsh,bsq,bwh,bwq
  REAL :: alv,alf,als,t0,t00,avc,afc,asc,rn1,bnd1,rn2,bnd2,             &
      rn3,rn4,rn5,rn6,rn7,rn8,rn9,rn10,rn101,rn10a,rn11,rn11a,          &
      rn12,rn12a(31),rn12b(31),rn13(31),rn14,rn15,rn15a,rn16,rn17,      &
      rn17a,rn17b,rn17c,rn18,rn18a,rn19,rn19a,rn19b,rn20,rn20a,rn20b,   &
      bnd3,rn21,rn22,rn23,rn23a,rn23b,rn25,rn25a(31),rn30a,rn30b,       &
      rn30c,rn31,beta,rn32

  COMMON/size/ tnw,tns,tng,roqr,roqs,roqg
  COMMON/cont/ c76,c358,c172,c409,c218,c580,c610,c149,c879,c141
  COMMON/bterv/ zrc,zgc,zsc,vrc,vgc,vsc
  COMMON/b3cs/ ag,bg,as,bs,aww,bww,bgh,bgq,bsh,bsq,bwh,bwq
  COMMON/bsnw/ alv,alf,als,t0,t00,avc,afc,asc,rn1,bnd1,rn2,bnd2,        &
      rn3,rn4,rn5,rn6,rn7,rn8,rn9,rn10,rn101,rn10a,rn11,rn11a,          &
      rn12,rn12a,rn12b,rn13,rn14,rn15,rn15a,rn16,rn17,                  &
      rn17a,rn17b,rn17c,rn18,rn18a,rn19,rn19a,rn19b,rn20,rn20a,rn20b,   &
      bnd3,rn21,rn22,rn23,rn23a,rn23b,rn25,rn25a,rn30a,rn30b,           &
      rn30c,rn31,beta,rn32

  REAL :: r19bt,r20t,r23t,betah,r10t,r11at,ee2,a1,a2,ee1,bs6,x3,        &
       cpcgs,rt0,d2t,x1,x2,bw3,bsh5,bgh5,bw6,bs3,bg3,bwh5
!
!-----------------------------------------------------------------------
!
!  Above constants defined in subroutine STCSTICE
!
!  The followings are locally used variables
!
!-----------------------------------------------------------------------
!
  REAL :: dep
  REAL :: dd
  REAL :: dd1
  REAL :: dm
  REAL :: ern
  REAL :: rsub1
  REAL :: cnd
  REAL :: pgwet
  REAL :: egs
  REAL :: esi
  REAL :: qsi
  REAL :: ssi
  REAL :: qsw
  REAL :: pihom
  REAL :: ssw
  REAL :: pidw
  REAL :: pimlt
  REAL :: psmlt
  REAL :: pgmlt
  REAL :: psdep
  REAL :: pssub
  REAL :: pgsub
  REAL :: pint
  REAL :: pidep
  REAL :: prn
  REAL :: psn
  REAL :: dlt1

  REAL :: y1
  REAL :: y2
  REAL :: y3
  REAL :: y4
  REAL :: y5

  INTEGER :: it, i,j,k
  REAL :: t1, t2, t3, t4
  REAL :: tem

  REAL :: temp0, temp, interp, f1, f2, rstep, scvstep
  INTEGER :: INDEX
  REAL :: rhobar0
!
!-----------------------------------------------------------------------
!
!  Define an inline function cvmgp(x1,x2,x3).
!
!  cvmgp=x1 for x3>=0.0
!  cvmgp=x2 for x3< 0.0
!
!-----------------------------------------------------------------------
!
  REAL :: eps
  REAL :: cvmgp
  cvmgp(x1,x2,x3)= x1*(0.5+SIGN(0.5,x3))+x2*(0.5-SIGN(0.5,x3))
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!

!
!-----------------------------------------------------------------------
!
!  Two water and three classes of ice-phase.
!
!-----------------------------------------------------------------------
!
  it=1

  cpcgs=1.003E7    ! Missing in Tao's original code
                   ! Since ARPS has defined cp, better to use that.

  rt0=1./(t0-t00)

  bw3=bww+3.
  bs3=bs+3.
  bg3=bg+3.
  bwh5=2.5+bwh
  bsh5=2.5+bsh
  bgh5=2.5+bgh
  bw6=bww+6.
  bs6=bs+6.
  betah=.5*beta
  r10t=rn10*d2t
  r11at=rn11a*d2t
  r19bt=rn19b*d2t
  r20t=-rn20*d2t
  r23t=-rn23*d2t

  rhobar0 = 1.275*0.001     ! Air density in cgs at 0 C and 1000 mb
  t1 = SQRT(rhobar0)

  t2 = 1.e5 * zrc
  t3 = 1.e5 * zsc
  t4 = 1.e5 * zgc
!
!-----------------------------------------------------------------------
!
!  Begin k loop for vertical dimension
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  The inverse of the increment of lookup tables for rhobar*q is held
!  by rstep, and that for (rhobar/mu/psi**2) by scvstep, where mu is
!  dynamic viscosity of air and psi is diffusivity of water vapor in
!  air.
!
!-----------------------------------------------------------------------
!
  rstep   = 1.0 / 50.0E-10
  scvstep = 1.0 / 0.3

  DO k=2,nz-2

    DO j = 1, ny-1
      DO i = 1, nx-1
        tem2d1(i,j) = SQRT(rhobar(i,j,k))
        tem2d2(i,j) = t1 / tem2d1(i,j)
      END DO
    END DO

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

        qc(i,j,k,2) = MAX(0.0, qc(i,j,k,2))
        qr(i,j,k,2) = MAX(0.0, qr(i,j,k,2))
        qi(i,j,k,2) = MAX(0.0, qi(i,j,k,2))
        qs(i,j,k,2) = MAX(0.0, qs(i,j,k,2))
        qh(i,j,k,2) = MAX(0.0, qh(i,j,k,2))

        tair(i,j) = (ptprt(i,j,k,2)+ptbar(i,j,k))*ppi(i,j,k)
        tairc(i,j) = tair(i,j)-t0
!
!-----------------------------------------------------------------------
!
!  Compute zr,zs,zg,vr,vs,vg
!
!-----------------------------------------------------------------------
!
        zr(i,j) = t2/SQRT(tem2d1(i,j))
        zs(i,j) = t3/SQRT(tem2d1(i,j))
        zg(i,j) = t4/SQRT(tem2d1(i,j))
        vr(i,j) = 0.0
        vs(i,j) = 0.0
        vg(i,j) = 0.0

!
!      for rain we have a more complicated modification
!       vtr =  4854*sqrt(rho0/rho)* [gamma(5)*lambda**(4)/
!                                   gamma(4)*(lambda+f)**(5)
!       

!        IF (qr(i,j,k,2) > 1.e-20) THEN
        IF (qr(i,j,k,2) > 1.e-16) THEN ! EMK
          dd = rhobar(i,j,k)*qr(i,j,k,2)
          y1 = SQRT(SQRT(dd))
          zr(i,j) = zrc/y1

          IF (fallopt == 2) THEN  ! Ferrier with updated coeff.
       
!            vr(i,j)=vrc*tem2d2(i,j)*(zr(i,j)**4)/((zr(i,j)+1.95)**5)

            temp  = min( 50.0e-6, max(0.0,dd) ) * rstep
            index = int(temp)
            f1 = pwrlam195ratio(index)
            f2 = pwrlam195ratio(index+1)
           vr(i,j)=vrc*tem2d2(i,j)*( f1 + (f2-f1)*(temp-index) )

          ELSE            ! original Lin scheme...

            temp  = min( 50.0e-6, max(0.0,dd) ) * rstep
            index = int(temp)
            f1 = pwr2(index)
            f2 = pwr2(index+1)
            vr(i,j) = vrc * tem2d2(i,j) * ( f1 + (f2-f1)*(temp-index) )

          END IF

        END IF

!        IF (qs(i,j,k,2) > 1.e-20) THEN
        IF (qs(i,j,k,2) > 1.e-16) THEN ! EMK
          dd = rhobar(i,j,k)*qs(i,j,k,2)
          y1 = SQRT(SQRT(dd))
          zs(i,j) = zsc/y1
          IF (fallopt == 2) THEN  ! Ferrier scheme with new coeff.
!            vs(i,j) = vsc*tem2d2(i,j)*(y1**.42)
!            vs(i,j) = vsc*tem2d2(i,j)*(dd**0.105) ! EMK

            temp  = min( 50.0e-6, max(0.0,dd) ) * rstep
            index = int(temp)
            f1 = pwr105(index)
            f2 = pwr105(index+1)
            vs(i,j) = vsc*tem2d2(i,j)*( f1 + (f2-f1)*(temp-index) ) ! EMK

          ELSE                  !  Lin scheme..
!            vs(i,j) = vsc*tem2d2(i,j)*sqrt(sqrt(y1))

            ! EMK...Adapted code from subroutine TERV
            temp  = min( 50.0e-6, max(0.0,dd) ) * rstep
            index = int(temp)
            f1 = pwr0625(index)
            f2 = pwr0625(index+1)
            vs(i,j) = vsc * tem2d2(i,j) * ( f1 + (f2-f1)*(temp-index) )

          END IF
        END IF

!  note zs is lambda for snow and vs is mass weighted mean fall
!  velocity for snow.

!        IF (qh(i,j,k,2) > 1.e-20) THEN
        IF (qh(i,j,k,2) > 1.e-16) THEN
          dd = rhobar(i,j,k)*qh(i,j,k,2)
          y1 = SQRT(SQRT(dd))
          zg(i,j) = zgc/y1

          IF (fallopt == 2) THEN  ! Ferrier scheme with new coeff.
!            vg(i,j) = vgc*tem2d2(i,j)*(y1**.6384)
!            vg(i,j) = vgc*tem2d2(i,j)*(dd**.1596) ! EMK

            temp  = min( 50.0e-6, max(0.0,dd) ) * rstep
            index = int(temp)
            f1 = pwr1596(index)
            f2 = pwr1596(index+1)
            vg(i,j) = vgc*tem2d2(i,j)*( f1 + (f2-f1)*(temp-index) )

          ELSE                  !  Lin scheme..
!            vg(i,j) = vgc/tem2d1(i,j)*sqrt(y1)

            ! EMK...Adapted code from subroutine TERV
            temp  = min( 50.0e-6, max(0.0,dd) ) * rstep
            index = int(temp)
            f1 = pwr0625(index)
            f2 = pwr0625(index+1)
            interp = f1 + (f2 - f1) * (temp - index)
            vg(i,j) = vgc / sqrt(rhobar(i,j,k)) * interp * interp

          END IF
        END IF
!
!-----------------------------------------------------------------------
!
!  y1  : dynamic viscosity of air (u)
!  dwv : diffusivity of water vapor in air (pi)
!  tca : thermal conductivity of air (ka)
!  y2  : kinetic viscosity (v)
!
!-----------------------------------------------------------------------
!
        y1 = c149*SQRT(tair(i,j)**3)/(tair(i,j)+120.)

        temp  = MIN( 150.0, MAX( 0.0, (tair(i,j)-173.15) ) )
        INDEX = INT(temp)
        f1 = pwr81(INDEX)
        f2 = pwr81(INDEX+1)
        interp = f1 + (f2-f1)*(temp-INDEX)

        dwv(i,j) = c879/pres(i,j,k)*tair(i,j) * interp
        tca(i,j) = c141*y1

        temp = MIN( 3000.0, rhobar(i,j,k)/(y1*dwv(i,j)**2) )            &
                * scvstep
        INDEX = INT( temp )
        f1 = pwr1666(INDEX)
        f2 = pwr1666(INDEX+1)

        scv(i,j) = f1 + (f2-f1)*(temp-INDEX)

      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!   1 * psaut : autoconversion of qi to qs
!   3 * psaci : accretion of qi to qs
!   4 * psacw : accretion of qc by qs (riming) (qsacw for psmlt)
!   5 * praci : accretion of qi by qr
!   6 * piacr : accretion of qr or qg by qi
!
!-----------------------------------------------------------------------
!
    DO j = 1,ny-1
      DO i = 1,nx-1

        psaut(i,j) = 0.0
        psaci(i,j) = 0.0
        praci(i,j) = 0.0
        piacr(i,j) = 0.0
        psacw(i,j) = 0.0
        qsacw(i,j) = 0.0

        tem = 1.0/zs(i,j)
        dd = tem**3 * SQRT(SQRT(tem))  ! dd= (1/zs)**3.25

        IF(qr(i,j,k,2) > 1.0E-20) THEN
          temp0 = rhobar(i,j,k)*qr(i,j,k,2)
        ELSE
          temp0 = rhobar(i,j,k)*1.0E-20
        END IF
        temp  = MIN( 50.0E-6, temp0 ) * rstep
        INDEX = INT(temp)
        f1 = pwr2(INDEX)
        f2 = pwr2(INDEX+1)
        interp = f1 + ( (f2 - f1) * (temp - INDEX) )

        IF (tair(i,j) < t0) THEN
          esi = EXP(.025*tairc(i,j))
          psaut(i,j) = AMAX1(rn1*esi*(qi(i,j,k,2)-bnd1) ,0.0)
          psaci(i,j) = rn3*tem2d2(i,j)*esi*qi(i,j,k,2)*dd
          psacw(i,j) = rn4*tem2d2(i,j)*qc(i,j,k,2)*dd
          praci(i,j) = rn5*tem2d2(i,j)*qi(i,j,k,2)*interp/zr(i,j)**3
          piacr(i,j) = rn6*tem2d2(i,j)*qi(i,j,k,2)*interp               &
                          *(zr(i,j))**(-6)
!                          *(1.0/zr(i,j))**6
        ELSE
          qsacw(i,j) = rn4*tem2d2(i,j)*qc(i,j,k,2)*dd
        END IF
!
!-----------------------------------------------------------------------
!
!   praut : autoconversion of qc to qr                         (21)
!   pracw : accretion of qc by qr                              (22)
!
!-----------------------------------------------------------------------
!
        pracw(i,j) = rn22*tem2d2(i,j)*qc(i,j,k,2)*interp/zr(i,j)**3
        praut(i,j) = 0.0
        y1 = qc(i,j,k,2)-bnd3

! Original
!        IF (y1 > 0.0) praut(i,j) = rhobar(i,j,k)*y1*y1/(1.2E-4+rn21/y1)
! EMK 12 May 2003 Density correction...see equation 27 in 
! Orville and Kopp, JAS, 1977.
        IF (y1 > 0.0) praut(i,j) = &
          rhobar(i,j,k)*y1*y1/(1.2E-4+rn21/(rhobar(i,j,k)*y1))
!
!-----------------------------------------------------------------------
!
!   psfw : bergeron processes for qs (koening, 1971)        (12)
!   psfi : bergeron processes for qs                        (13)
!
!-----------------------------------------------------------------------
!
        psfw(i,j) = 0.0
        psfi(i,j) = 0.0

        IF (tair(i,j) < t0) THEN
          y1 = AMAX1( AMIN1(tairc(i,j), -1.), -31.)
          it = INT(ABS(y1))
          y1 = rn12a(it)
          y2 = rn12b(it)
          psfw(i,j) = AMAX1(d2t*y1*(y2+rn12*rhobar(i,j,k)               &
                     *qc(i,j,k,2))*qi(i,j,k,2),0.0)
          psfi(i,j) = rn13(it)*qi(i,j,k,2)
        END IF
!
!-----------------------------------------------------------------------
!
!  qg = qg+min(pgdry,pgwet)
!  pgacs : accretion of qs by qg (dgacs,wgacs: dry and wet)    (9)
!  dgacw : accretion of qc by qg (qgacw for pgmlt)            (14)
!  dgacr : accretion of qr to qg (qgacr for pgmlt)            (16)
!
!-----------------------------------------------------------------------
!
        ee1 = 1.
        ee2 = 0.09
        egs = ee1*EXP(ee2*tairc(i,j))

        IF (tair(i,j) >= t0) egs = 1.0
        y1 = ABS(vg(i,j)-vs(i,j))
        y2 = zs(i,j)*zg(i,j)
        y3 = 5./y2
        y4 = .08*y3*y3
        y5 = .05*y3*y4
        dd = y1*(y3/zs(i,j)**5+y4/zs(i,j)**3+y5/zs(i,j))
        pgacs(i,j) = rn9/rhobar(i,j,k)*egs*dd
        dgacs(i,j) = pgacs(i,j)
        wgacs(i,j) = rn9/rhobar(i,j,k)*dd
        y1 = 1.0 / ( zg(i,j)**3 * SQRT(zg(i,j)) )
        dgacw(i,j) = AMAX1(rn14/tem2d1(i,j)*qc(i,j,k,2)*y1, 0.0)
        qgacw(i,j) = dgacw(i,j)
        y1 = ABS(vg(i,j)-vr(i,j))
        y2 = zr(i,j)*zg(i,j)
        y3 = 5./y2
        y4 = .08*y3*y3
        y5 = .05*y3*y4
        dd = rn16/rhobar(i,j,k)                                         &
                *y1*(y3/zr(i,j)**5+y4/zr(i,j)**3                        &
                +y5/zr(i,j))
        dgacr(i,j) = AMAX1(dd, 0.0)
        qgacr(i,j) = dgacr(i,j)

        IF (tair(i,j) >= t0) THEN
          dgacs(i,j) = 0.0
          wgacs(i,j) = 0.0
          dgacw(i,j) = 0.0
          dgacr(i,j) = 0.0
        ELSE
          pgacs(i,j) = 0.0
          qgacw(i,j) = 0.0
          qgacr(i,j) = 0.0
        END IF
!
!-----------------------------------------------------------------------
!
!  dgaci : accretion of qi by qg (wgaci for wet growth)       (15)
!  pgwet : wet growth of qg                                   (17)
!
!-----------------------------------------------------------------------
!
        dgaci(i,j) = 0.0
        wgaci(i,j) = 0.0
        pgwet = 0.0

        IF (tair(i,j) < t0) THEN

          y1 = qi(i,j,k,2)/( zg(i,j)**3 * SQRT(zg(i,j)) )
          dgaci(i,j) = rn15/tem2d1(i,j)*y1
          wgaci(i,j) = rn15a/tem2d1(i,j)*y1
          y1 = 1./(alf+rn17c*tairc(i,j))

          IF(qh(i,j,k,2) > 1.0E-20) THEN
            temp0 = rhobar(i,j,k)*qh(i,j,k,2)
          ELSE
            temp0 = rhobar(i,j,k)*1.0E-20
          END IF
          temp  = MIN( 50.0E-6, temp0 ) * rstep
          INDEX = INT(temp)
          f1 = pwr0625(INDEX)
          f2 = pwr0625(INDEX+1)
          interp = f1 + (f2 - f1) * (temp - INDEX)

          y3 = .78/zg(i,j)**2+rn17a/SQRT(tem2d1(i,j))                   &
               *scv(i,j)/( zg(i,j)**3 * interp )
          y4 = rhobar(i,j,k)*alv*dwv(i,j)                               &
               *(3.799052E3/pres(i,j,k)-(qv(i,j,k,2)))                  &
               -tca(i,j)*tairc(i,j)
          dd = y1*(rn17/rhobar(i,j,k)*y4*y3                             &
               +(wgaci(i,j)+wgacs(i,j))*(alf+rn17b*tairc(i,j)))
          pgwet = AMAX1(dd, 0.0)

        END IF
!
!-----------------------------------------------------------------------
!
!  Shed process (wgacr = pgwet-dgacw-wgaci-wgacs)
!
!-----------------------------------------------------------------------
!

        wgacr(i,j) = pgwet-dgacw(i,j)-wgaci(i,j)-wgacs(i,j)
        y2 = dgacw(i,j)+dgaci(i,j)+dgacr(i,j)+dgacs(i,j)

        IF (pgwet >= y2) THEN
          wgacr(i,j) = 0.0
          wgaci(i,j) = 0.0
          wgacs(i,j) = 0.0
        ELSE
          dgacr(i,j) = 0.0
          dgaci(i,j) = 0.0
          dgacs(i,j) = 0.0
        END IF

      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  Handling the negative cloud water (qc)
!  Handling the negative cloud ice (qi)
!
!-----------------------------------------------------------------------
!
    eps = 1.0E-30

    DO j = 1,ny-1
      DO i = 1,nx-1
        y1 = (psacw(i,j)+praut(i,j)+pracw(i,j)+psfw(i,j)+               &
              dgacw(i,j)+qsacw(i,j)+qgacw(i,j))*d2t

        IF (qc(i,j,k,3) < y1) THEN
          a1 = qc(i,j,k,3)/(y1+eps)
          psacw(i,j) = psacw(i,j)*a1
          praut(i,j) = praut(i,j)*a1
          pracw(i,j) = pracw(i,j)*a1
          psfw (i,j) = psfw (i,j)*a1
          dgacw(i,j) = dgacw(i,j)*a1
          qsacw(i,j) = qsacw(i,j)*a1
          qgacw(i,j) = qgacw(i,j)*a1
          qc(i,j,k,3)  = 0.0
        ELSE
          qc(i,j,k,3)  = qc(i,j,k,3)-y1
        END IF

        y2 = (psaut(i,j)+psaci(i,j)+praci(i,j)+psfi(i,j)+               &
              dgaci(i,j)+wgaci(i,j))*d2t

        IF (y2 > qi(i,j,k,3)) THEN
          a2 = qi(i,j,k,3)/(y2+eps)
          psaut(i,j) = psaut(i,j)*a2
          psaci(i,j) = psaci(i,j)*a2
          praci(i,j) = praci(i,j)*a2
          psfi(i,j) = psfi(i,j)*a2
          dgaci(i,j) = dgaci(i,j)*a2
          wgaci(i,j) = wgaci(i,j)*a2
          qi(i,j,k,3) = 0.0
        ELSE
          qi(i,j,k,3) = qi(i,j,k,3)-y2
        END IF

        dlt3(i,j) = 0.0
        dlt2(i,j) = 0.0

        IF (tair(i,j) < t0) THEN

          IF (qr(i,j,k,2) < 1.e-4) THEN
            dlt3(i,j) = 1.0
            dlt2(i,j) = 1.0
          END IF

          IF (qs(i,j,k,2) >= 1.e-4) dlt2(i,j) = 0.0

        END IF

        pr(i,j) = (qsacw(i,j)+praut(i,j)+pracw(i,j)+                    &
                   qgacw(i,j))*d2t
        ps(i,j) = (psaut(i,j)+psaci(i,j)+psacw(i,j)+                    &
                   psfw(i,j)+psfi(i,j)+dlt3(i,j)*praci(i,j))*d2t
        pg(i,j) = ((1.-dlt3(i,j))*praci(i,j)+dgaci(i,j)+                &
                   wgaci(i,j)+dgacw(i,j))*d2t

      END DO
    END DO


!
!-----------------------------------------------------------------------
!
!  pracs : accretion of qs by qr                               (7)
!  psacr : accretion of qr by qs (qsacr for psmlt)             (8)
!
!-----------------------------------------------------------------------
!

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

        y1 = ABS(vr(i,j)-vs(i,j))
        y2 = zr(i,j)*zs(i,j)
        y3 = 5./y2
        y4 = .08*y3*y3
        y5 = .05*y3*y4
        pracs(i,j) = rn7/rhobar(i,j,k)                                  &
                   *y1*(y3/zs(i,j)**5+y4/zs(i,j)**3                     &
                   +y5/zs(i,j))
        psacr(i,j) = rn8/rhobar(i,j,k)                                  &
                   *y1*(y3/zr(i,j)**5+y4/zr(i,j)**3                     &
                   +y5/zr(i,j))
        qsacr(i,j) = psacr(i,j)

        IF (tair(i,j) >= t0) THEN
          pracs(i,j) = 0.0
          psacr(i,j) = 0.0
        ELSE
          qsacr(i,j) = 0.0
        END IF
!
!-----------------------------------------------------------------------
!
!  pgaut : autoconversion of qs to qg                          (2)
!  pgfr  : freezing of qr to qg                               (18)
!
!-----------------------------------------------------------------------
!
        pgaut(i,j) = 0.0
        pgfr (i,j) = 0.0

        IF (tair(i,j) < t0) THEN
          y2 = EXP(rn18a*(t0-tair(i,j)))
          pgfr(i,j) = AMAX1(rn18/rhobar(i,j,k)*(y2-1.)                  &
                       *(zr(i,j))**(-7), 0.0)
!                       *(1.0/zr(i,j))**7, 0.0)
        END IF

      END DO
    END DO

!-----------------------------------------------------------------------
!
!  Handling the negative rain water (qr)
!  Handling the negative snow (qs)
!
!-----------------------------------------------------------------------

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

        y1 = (piacr(i,j)+dgacr(i,j)+wgacr(i,j)+psacr(i,j)+              &
              pgfr(i,j))*d2t

        tem = qr(i,j,k,3)+pr(i,j)
        IF (tem < y1) THEN
          a1 = tem/(y1+eps)
          piacr(i,j) = piacr(i,j)*a1
          dgacr(i,j) = dgacr(i,j)*a1
          wgacr(i,j) = wgacr(i,j)*a1
          pgfr (i,j) = pgfr (i,j)*a1
          psacr(i,j) = psacr(i,j)*a1
          qr(i,j,k,3)  = 0.0
        ELSE
          qr(i,j,k,3)  = qr(i,j,k,3)+pr(i,j)-y1
        END IF

        prn = d2t*((1.-dlt3(i,j))*piacr(i,j)+dgacr(i,j)+                &
              wgacr(i,j)+(1.-dlt2(i,j))*psacr(i,j)+pgfr(i,j))
        ps(i,j) = ps(i,j)+d2t*(dlt3(i,j)*piacr(i,j)+                    &
               dlt2(i,j)* psacr(i,j))
        pracs(i,j) = (1.-dlt2(i,j))*pracs(i,j)
        psn = d2t*(pgacs(i,j)+dgacs(i,j)+wgacs(i,j)+pgaut(i,j)+         &
               pracs(i,j))

        tem = qs(i,j,k,3)+ps(i,j)
        qs(i,j,k,3) = qs(i,j,k,3)+ps(i,j)-psn

        IF (qs(i,j,k,3) < 0.0) THEN

          IF ( psn > 0.0 ) THEN
            a2 = tem/psn
            pgacs(i,j) = pgacs(i,j)*a2
            dgacs(i,j) = dgacs(i,j)*a2
            wgacs(i,j) = wgacs(i,j)*a2
            pgaut(i,j) = pgaut(i,j)*a2
            pracs(i,j) = pracs(i,j)*a2
            psn = psn*a2
          ELSE
            psn = psn + qs(i,j,k,3)
          END IF

          qs(i,j,k,3) = 0.0

        END IF

        y2 = d2t*(psacw(i,j)+psfw(i,j)+dgacw(i,j)+piacr(i,j)+           &
              dgacr(i,j)+wgacr(i,j)+psacr(i,j)+pgfr(i,j))
        ptprt(i,j,k,3) = ptprt(i,j,k,3)+afc/ppi(i,j,k)*y2
        qh   (i,j,k,3) = qh(i,j,k,3)+pg(i,j)+prn+psn

      END DO
    END DO


!
!-----------------------------------------------------------------------
!
!  psmlt : melting of qs                                      (11)
!  pgmlt : melting of qg to qr                                (19)
!
!-----------------------------------------------------------------------
!

!   write (*,*) "ZUWEN subsatopt/rhsat", subsatopt, rhsat 

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

        psmlt = 0.0
        pgmlt = 0.0
!C        tair(i,j) = (ptprt(i,j,k,2)+ptbar(i,j,k))*ppi(i,j,k)

        IF (tair(i,j) >= t0) THEN

!C          tairc(i,j) = tair(i,j)-t0
          y1 = tca(i,j)*tairc(i,j)-rhobar(i,j,k)*alv                    &
               *dwv(i,j)*(3.799052E3/pres(i,j,k)-(qv(i,j,k,2)))

          IF(qs(i,j,k,2) > 1.0E-20) THEN
            temp0 = rhobar(i,j,k)*qs(i,j,k,2)
          ELSE
            temp0 = rhobar(i,j,k)*1.0E-20
          END IF
          temp  = MIN( 50.0E-6, temp0 ) * rstep
          INDEX = INT(temp)
          f1 = pwr15625(INDEX)
          f2 = pwr15625(INDEX+1)
          interp = f1 + (f2 - f1) * (temp - INDEX)

          y2 = .78/zs(i,j)**2+rn101*SQRT(tem2d2(i,j))                   &
               *scv(i,j)*interp/zs(i,j)**2
          dd = rn11/rhobar(i,j,k)*d2t*y1*y2                             &
               +r11at*tairc(i,j)*(qsacw(i,j)+qsacr(i,j))
          psmlt = AMAX1(0.0, AMIN1(dd, qs(i,j,k,3)))

          IF(qh(i,j,k,2) > 1.0E-20) THEN
            temp0 = rhobar(i,j,k)*qh(i,j,k,2)
          ELSE
            temp0 = rhobar(i,j,k)*1.0E-20
          END IF
          temp  = MIN( 50.0E-6, temp0 ) * rstep
          INDEX = INT(temp)
          f1 = pwr0625(INDEX)
          f2 = pwr0625(INDEX+1)
          interp = f1 + (f2 - f1) * (temp - INDEX)

          y3  = .78/zg(i,j)**2+rn19a/SQRT(tem2d1(i,j))                  &
               *scv(i,j)/( zg(i,j)**3 * interp )
          dd1 = rn19/rhobar(i,j,k)*d2t*y1*y3                            &
                +r19bt*tairc(i,j)*(qgacw(i,j)+qgacr(i,j))
          pgmlt = AMAX1(0.0, AMIN1(dd1, qh(i,j,k,3)))
          ptprt(i,j,k,3) = ptprt(i,j,k,3)-afc/ppi(i,j,k)*(psmlt+pgmlt)
          qr(i,j,k,3) = qr(i,j,k,3)+psmlt+pgmlt
          qs(i,j,k,3) = qs(i,j,k,3)-psmlt
          qh(i,j,k,3) = qh(i,j,k,3)-pgmlt

        END IF
!
!-----------------------------------------------------------------------
!
!  pihom : homogeneous freezing of qc to qi (t < t00)         (24)
!  pidw : deposition growth of qc to qi ( t0 < t < =  t00)    (25)
!  pimlt : melting of qi to qc (t > =  t0)                    (26)
!
!-----------------------------------------------------------------------

!C        IF (qc(i,j,k,2).le.1.e-20) qc(i,j,k,2) = 0.0
!C        IF (qi(i,j,k,2).le.1.e-20) qi(i,j,k,2) = 0.0

        tair3(i,j) = (ptprt(i,j,k,3)+ptbar(i,j,k))*ppi(i,j,k)
        pihom = cvmgp(0.,qc(i,j,k,3),tair3(i,j)-t00)
        pimlt = cvmgp(qi(i,j,k,3),0.,tair3(i,j)-t0)
        pidw = 0.0

        IF (tair(i,j) < t0 .AND. tair(i,j) > t00) THEN
!C          tairc(i,j) = tair(i,j)-t0
          y1 = AMAX1( AMIN1(tairc(i,j), -1.), -31.)
          it = INT(ABS(y1))
          y2 = rn25a(it)
          y3 = EXP(.5*ABS(tairc(i,j)))
          pidw = AMIN1(rn25/rhobar(i,j,k)*d2t*y2*                       &
               y3, qc(i,j,k,3))
        END IF

        y1 = pihom-pimlt+pidw
        ptprt(i,j,k,3) = ptprt(i,j,k,3)+afc/ppi(i,j,k)*y1
        qc(i,j,k,3) = qc(i,j,k,3)-y1
        qi(i,j,k,3) = qi(i,j,k,3)+y1
!
!-----------------------------------------------------------------------
!
!  pint  : initiation of qi                                   (31)
!  pidep : deposition of qi                                   (32)
!
!  Note:   rhsat is the minimum RH for condensation to occur. 
!          It is used to adjust saturated specific humidity, 
!          qsw and qsi, in the following code. Zuwen, 04/2002
!
!-----------------------------------------------------------------------
!
        pint = 0.0
!C        tair(i,j) = (ptprt(i,j,k,2)+ptbar(i,j,k))*ppi(i,j,k)

        IF (tair(i,j) < t0) THEN
!C          tairc(i,j) = tair(i,j)-t0
          dd = rn31/rhobar(i,j,k)*EXP(beta*tairc(i,j))
          rtair(i,j) = 1./(tair(i,j)-c76)
          y2 = EXP(c218-c580*rtair(i,j))
          qsi = rhsat*(3.799052E3/pres(i,j,k))*y2
          esi = c610*y2
          ssi = (qv(i,j,k,2))/qsi-1.
          dm = AMAX1( (qv(i,j,k,3)-qsi), 0.)
          rsub1 = c580*asc*qsi*rtair(i,j)*rtair(i,j)
          pint = AMIN1(dd,dm)
          y1 = 1./tair(i,j)
          y2 = EXP(betah*tairc(i,j))
          IF(qi(i,j,k,2) > 1.0E-30) THEN
            y3 = SQRT(qi(i,j,k,2))
          ELSE
            y3 = 1.0E-15
          END IF
          dd = y1*(rn30a*y1-rn30b)+rn30c*                               &
               tair(i,j)/esi
          pidep = AMAX1(rn32*d2t/tem2d1(i,j)                            &
               *ssi*y2*y3/dd, 0.0)
          pint = pint+pidep
          dep = dm/(1.+rsub1)
          pint = AMIN1(pint,dep)
          ptprt(i,j,k,3) = ptprt(i,j,k,3)+asc/ppi(i,j,k)*pint
          qv(i,j,k,3) = qv(i,j,k,3)-pint
          qi(i,j,k,3) = qi(i,j,k,3)+pint
        END IF

      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  Tao et al (1989) saturation technique
!
!-----------------------------------------------------------------------
!
    DO j = 1,ny-1
      DO i = 1,nx-1

        tair3(i,j) = (ptprt(i,j,k,3)+ptbar(i,j,k))*ppi(i,j,k)
        cnd = rt0*(tair3(i,j)-t00)
        dep = rt0*(t0-tair3(i,j))
        y1 = 1./(tair3(i,j)-c358)
        y2 = 1./(tair3(i,j)-c76)
        qsw = rhsat*(3.799052E3/pres(i,j,k))*EXP(c172-c409*y1)
        qsi = rhsat*(3.799052E3/pres(i,j,k))*EXP(c218-c580*y2)
        dd = c409*ppi(i,j,k)*y1*y1
        dd1 = c580*ppi(i,j,k)*y2*y2

        IF (qc(i,j,k,3) <= 1.e-20) qc(i,j,k,3) = 1.e-20
        IF (qi(i,j,k,3) <= 1.e-20) qi(i,j,k,3) = 1.e-20

        IF (tair3(i,j) >= t0) THEN
          dep = 0.0
          cnd = 1.0
          qi(i,j,k,3) = 0.0
        END IF

        IF (tair3(i,j) <= t00) THEN
          cnd = 0.0
          dep = 1.0
          qc(i,j,k,3) = 0.0
        END IF

        y5 = avc/ppi(i,j,k)*cnd+asc/ppi(i,j,k)*dep
        y1 = qc(i,j,k,3)*qsw/(qc(i,j,k,3)+qi(i,j,k,3))
        y2 = qi(i,j,k,3)*qsi/(qc(i,j,k,3)+qi(i,j,k,3))
        y4 = dd*y1+dd1*y2
        dm = qv(i,j,k,3)-y1-y2
        rsub1 = dm/(1.+y4*y5)
        cnd = cnd*rsub1
        dep = dep*rsub1

        IF (qc(i,j,k,3) <= 1.e-20) qc(i,j,k,3) = 0.
        IF (qi(i,j,k,3) <= 1.e-20) qi(i,j,k,3) = 0.
!
!-----------------------------------------------------------------------
!
!  Condensation or evaporation of qc
!
!-----------------------------------------------------------------------
!
        cnd = AMAX1(-qc(i,j,k,3),cnd)
!
!-----------------------------------------------------------------------
!
!  Deposition or sublimation of qi
!
!-----------------------------------------------------------------------
!
        dep = AMAX1(-qi(i,j,k,3),dep)
        ptprt(i,j,k,3) = ptprt(i,j,k,3)+avc/ppi(i,j,k)*                 &
                     cnd+asc/ppi(i,j,k)*dep
        qv(i,j,k,3) = qv(i,j,k,3)-cnd-dep
        qc(i,j,k,3) = qc(i,j,k,3)+cnd
        qi(i,j,k,3) = qi(i,j,k,3)+dep

      END DO
    END DO
!
!-----------------------------------------------------------------------
!
!  psdep : deposition or sublimation of qs                    (10)
!  pgsub : sublimation of qg                                  (20)
!
!-----------------------------------------------------------------------
!
    DO j = 1,ny-1
      DO i = 1,nx-1

        psdep = 0.0
        pssub = 0.0
        pgsub = 0.0
!C        tair(i,j) = (ptprt(i,j,k,2)+ptbar(i,j,k))*ppi(i,j,k)

        IF (tair(i,j) < t0) THEN

          IF (qs(i,j,k,2) < 1.e-20) qs(i,j,k,2) = 0.0
          IF (qh(i,j,k,2) < 1.e-20) qh(i,j,k,2) = 0.0

          rtair(i,j) = 1./(tair(i,j)-c76)
          qsi = rhsat*(3.799052E3/pres(i,j,k))*EXP(c218-c580*rtair(i,j))
          ssi = (qv(i,j,k,2))/qsi-1.
          y1 = rn10a*rhobar(i,j,k)/(tca(i,j)*tair(i,j)**2)              &
                +1./(dwv(i,j)*qsi)

          IF(qs(i,j,k,2) > 1.0E-20) THEN
            temp0 = rhobar(i,j,k)*qs(i,j,k,2)
          ELSE
            temp0 = rhobar(i,j,k)*1.0E-20
          END IF
          temp  = MIN( 50.0E-6, temp0 ) * rstep
          INDEX = INT(temp)
          f1 = pwr15625(INDEX)
          f2 = pwr15625(INDEX+1)
          interp = f1 + (f2 - f1) * (temp - INDEX)

          y2 = .78/zs(i,j)**2+rn101*SQRT(tem2d2(i,j))                   &
               *scv(i,j)*interp/zs(i,j)**2
          psdep = r10t*ssi*y2/y1
          pssub = psdep
          psdep = AMAX1(psdep, 0.)
          pssub = AMAX1(-qs(i,j,k,3), AMIN1(pssub, 0.))

          IF(qh(i,j,k,2) > 1.0E-20) THEN
            temp0 = rhobar(i,j,k)*qh(i,j,k,2)
          ELSE
            temp0 = rhobar(i,j,k)*1.0E-20
          END IF
          temp  = MIN( 50.0E-6, temp0 ) * rstep
          INDEX = INT(temp)
          f1 = pwr0625(INDEX)
          f2 = pwr0625(INDEX+1)
          interp = f1 + (f2 - f1) * (temp - INDEX)

          y2 = .78/zg(i,j)**2+rn20b/SQRT(tem2d1(i,j))                   &
                *scv(i,j)/( zg(i,j)**3 * interp )
          pgsub = r20t*ssi*y2/y1
          dm = qv(i,j,k,2)-qsi
          rsub1 = c580*asc*qsi*rtair(i,j)*rtair(i,j)
!
!-----------------------------------------------------------------------
!
!  Deposition or sublimation of qs
!
!-----------------------------------------------------------------------
!
          y1 = dm/(1.+rsub1)
          psdep = AMIN1(psdep,AMAX1(y1,0.))
          y2 = AMIN1(y1,0.)
          pssub = AMAX1(pssub,y2)
!
!-----------------------------------------------------------------------
!
!  Sublimation of qg
!
!-----------------------------------------------------------------------
!
          dd = AMAX1((-y2-qs(i,j,k,3)), 0.)
          pgsub = AMIN1(dd, qh(i,j,k,3), AMAX1(pgsub,0.))
          dlt1 = cvmgp(1.,0.,qc(i,j,k,2)+qi(i,j,k,2)-1.e-8)
          psdep = dlt1*psdep
          pssub = (1.-dlt1)*pssub
          pgsub = (1.-dlt1)*pgsub
          ptprt(i,j,k,3) = ptprt(i,j,k,3)+asc/ppi(i,j,k)                &
                *(psdep+ pssub-pgsub)
          qv(i,j,k,3) = qv(i,j,k,3)+pgsub-pssub-psdep
          qs(i,j,k,3) = qs(i,j,k,3)+psdep+pssub
          qh(i,j,k,3) = qh(i,j,k,3)-pgsub

        END IF

!
!-----------------------------------------------------------------------
!
!* 23 * ern : evaporation of qr (subsaturation)                    (23**
!
!-----------------------------------------------------------------------
!
        ern = 0.0

        IF (qr(i,j,k,2) > 1.e-20) THEN

!C          tair(i,j) = (ptprt(i,j,k,2)+ptbar(i,j,k))*ppi(i,j,k)
          rtair(i,j) = 1./(tair(i,j)-c358)
          qsw = rhsat*(3.799052E3/pres(i,j,k))*EXP(c172-c409*rtair(i,j))
          ssw = (qv(i,j,k,2))/qsw-1.0
          dm = qv(i,j,k,3)-qsw
          rsub1 = c409*avc*qsw*rtair(i,j)*rtair(i,j)
          dd1 = AMAX1(-dm/(1.+rsub1), 0.0)
!wdt update: (see arpssupport emails by Vince Wong & Eric Kemp
!"Re: micro_ice3d.f (ARPS Version 4.5.2.)" 9 Jul 2001)
          y1 = .78/zr(i,j)**2+rn23a*SQRT(tem2d2(i,j))                   &
                  *scv(i,j)/zr(i,j)**2.9
          y2 = rn23b*rhobar(i,j,k)/(tca(i,j)*tair(i,j)**2)              &
                +1./(dwv(i,j)*qsw)
          ern = r23t*ssw*y1/y2
          ern = AMIN1(dd1,qr(i,j,k,3),AMAX1(ern,0.))
          ptprt(i,j,k,3) = ptprt(i,j,k,3)-avc/ppi(i,j,k)*ern
          qv(i,j,k,3) = qv(i,j,k,3)+ern
          qr(i,j,k,3) = qr(i,j,k,3)-ern

        END IF

        qc(i,j,k,3) = MAX(0., qc(i,j,k,3))
        qs(i,j,k,3) = MAX(0., qs(i,j,k,3))
        qi(i,j,k,3) = MAX(0., qi(i,j,k,3))
        qh(i,j,k,3) = MAX(0., qh(i,j,k,3))
        qr(i,j,k,3) = MAX(0., qr(i,j,k,3))

      END DO
    END DO

  END DO

  RETURN
END SUBROUTINE icecvt
!
!
!##################################################################
!##################################################################
!######                                                      ######
!######                 SUBROUTINE TERV                      ######
!######                                                      ######
!######                     Developed by                     ######
!######                                                      ######
!######    Goddard Cumulus Ensemble Modeling Group, NASA     ######
!######                                                      ######
!######     Center for Analysis and Prediction of Storms     ######
!######               University of Oklahoma                 ######
!######                                                      ######
!##################################################################
!##################################################################
!
!


SUBROUTINE terv(nx,ny,nz,irsg,rhobar,q,vtr) 1
!
!-----------------------------------------------------------------------
!
!  PURPOSE:
!
!  Compute the terminal velocities (vtr) of rain water qr, snow qs and
!  hail qg.
!
!-----------------------------------------------------------------------
!
!  AUTHOR: Goddard Cumulus Ensemble Modeling Group, NASA
!
!  MODIFICATION HISTORY:
!
!  10/20/1996 (M. Xue)
!  Rewritten for ARPS.
!
!  02/24/1997 (J. Zong, M. Xue and Yuhe Liu)
!  Power calculations are replaced by lookup table functions for
!  terminal velocity.
!
!  08/27/2002 (D. Weber)
!  Added fallopt option to provide a choice of fall velocity   
!  computations.          
!
!-----------------------------------------------------------------------
!

!
!-----------------------------------------------------------------------
!
!  INPUT:
!
!  irsg    Flag identifying hydrometeor field as rain, snow or hail
!          =  0 for rain;  qpt is rain
!          =  1 for snow;  qpt is snow
!          =  2 for hail;  qpt is hail
!  q       Hydrometeor field defined at the scalar point
!  rhobar  Air density defined at scalar point (g/cm**3)
!
!  OUTPUT:
!
!  vtr     Vertical velocity defined at the scalar point (cm/s)
!
!-----------------------------------------------------------------------
!
  INTEGER :: nx,ny,nz
  REAL :: q     (nx,ny,nz)
  REAL :: rhobar(nx,ny,nz)
  REAL :: vtr   (nx,ny,nz)
  REAL :: tema, temb, temc, temd, rho0cgs
  REAL :: temp, interp, f1, f2, rstep
  INTEGER :: INDEX
!
!-----------------------------------------------------------------------
!
!  Common variables. Defined in subroutine STCSTICE
!
!-----------------------------------------------------------------------
!
  COMMON/bterv/ zrc,zgc,zsc,vrc,vgc,vsc
  COMMON/b3cs/ ag,bg,as,bs,aww,bww,bgh,bgq,bsh,bsq,bwh,bwq

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

!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
!  Beginning of executable code...
!
!@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
!
  rho0cgs = rho0*0.001

  rstep = 1.0/50.0E-10

  temc = 0.0

  IF (irsg == 0) THEN    ! irsg = 0    for rain water (qr)

    DO k = 1,nz-1
      DO j = 1,ny-1
        DO i = 1,nx-1
          tema=SQRT(rho0cgs/rhobar(i,j,k))
          temb=rhobar(i,j,k)*q(i,j,k)

          IF (temb > 1.e-16) THEN

            IF (fallopt == 2) THEN  ! ferrier numbers.....

               temd=zrc/(sqrt(sqrt(temb)))
!              vtr(i,j,k)=vrc*tema*(temd**4)/((temd+1.95)**5)

              temp  = min( 50.0e-6, max(0.0,temb) ) * rstep
              index = int(temp)
              f1 = pwrlam195ratio(index)
              f2 = pwrlam195ratio(index+1)

              vtr(i,j,k)=vrc*tema*(f1 + (f2-f1)*(temp-index))

            ELSE   !  original Lin scheme

              temp  = min( 50.0e-6, max(0.0,temb) ) * rstep
              index = int(temp)
              f1 = pwr2(index)
              f2 = pwr2(index+1)

              vtr(i,j,k) = vrc * tema * (f1 + (f2-f1)*(temp-index))

            END IF

          ELSE
            vtr(i,j,k) = 0.0
          END IF
          temc = max(temc,vtr(i,j,k))
        END DO
      END DO
    END DO

!   print *,'max rain fall speed is=',temc,irsg

  ELSE IF( irsg == 1) THEN  !irsg = 1 for snow (qs)

    DO k = 1,nz-1
      DO j = 1,ny-1
        DO i = 1,nx-1
          tema=SQRT(rho0cgs/rhobar(i,j,k))
          temb=rhobar(i,j,k)*q(i,j,k)

          IF (temb > 1.e-16) THEN

            IF (fallopt == 2) THEN  ! Ferrier numbers.....
!              vtr(i,j,k) = vsc * tema * (temb**.105)
              vtr(i,j,k) = vsc * tema * (temb**.105)

              temp  = min( 50.0e-6, max(0.0,temb) ) * rstep
              index = int(temp)
              f1 = pwr105(index)
              f2 = pwr105(index+1)
              vtr(i,j,k) = vsc * tema * ( f1 + (f2-f1)*(temp-index) )

            ELSE          ! original Lin scheme.......
              temp  = min( 50.0e-6, max(0.0,temb) ) * rstep
              index = int(temp)
              f1 = pwr0625(index)
              f2 = pwr0625(index+1)
              vtr(i,j,k) = vsc * tema * ( f1 + (f2-f1)*(temp-index) )
            END IF

          ELSE
            vtr(i,j,k) = 0.0
          END IF

          temc = max(temc,vtr(i,j,k))
        END DO
      END DO
    END DO

!   print *,'max snow fall speed is=',temc,irsg

  ELSE IF( irsg == 2) THEN  ! irsg = 2   for graupel (qg)

    DO k = 1,nz-1
      DO j = 1,ny-1
        DO i = 1,nx-1
          tema=sqrt(rho0cgs/rhobar(i,j,k))
          temb=rhobar(i,j,k)*q(i,j,k)
          IF (temb > 1.e-16) THEN

            IF (fallopt == 2) THEN  ! new method

!              vtr(i,j,k) = vgc * tema * (temb**.1596)

              temp  = min( 50.0e-6, max(0.0,temb) ) * rstep
              index = int(temp)
              f1 = pwr1596(index)
              f2 = pwr1596(index+1)
              vtr(i,j,k) = vgc * tema * ( f1 + (f2-f1)*(temp-index) )

            ELSE    ! Original method
              temp  = min( 50.0e-6, max(0.0,temb) ) * rstep
              index = int(temp)
              f1 = pwr0625(index)
              f2 = pwr0625(index+1)
              interp = f1 + (f2 - f1) * (temp - index)
              vtr(i,j,k) = vgc / sqrt(rhobar(i,j,k)) * interp * interp
            END IF
            temc = max(temc,vtr(i,j,k))

          ELSE
            vtr(i,j,k) = 0.0
          END IF
            temc = max(temc,vtr(i,j,k))
        END DO
      END DO
    END DO

!   print *,'max hail fall speed is=',temc,irsg

  END IF

  RETURN
END SUBROUTINE terv