
!###############################################################################
! ***** VERSION OF MICROPHYSICS DESIGNED FOR HIGHER RESOLUTION MESO ETA MODEL
!       (1) Represents sedimentation by preserving a portion of the precipitation
!           through top-down integration from cloud-top.  Modified procedure to
!           Zhao and Carr (1997).
!       (2) Microphysical equations are modified to be less sensitive to time
!           steps by use of Clausius-Clapeyron equation to account for changes in
!           saturation mixing ratios in response to latent heating/cooling.
!       (3) Prevent spurious temperature oscillations across 0C due to
!           microphysics.
!       (4) Uses lookup tables for: calculating two different ventilation
!           coefficients in condensation and deposition processes; accretion of
!           cloud water by precipitation; precipitation mass; precipitation rate
!           (and mass-weighted precipitation fall speeds).
!       (5) Assumes temperature-dependent variation in mean diameter of large ice
!           (Houze et al., 1979; Ryan et al., 1996).
!        -> 8/22/01: This relationship has been extended to colder temperatures
!           to parameterize smaller large-ice particles down to mean sizes of MDImin,
!           which is 50 microns reached at -55.9C.
!       (6) Attempts to differentiate growth of large and small ice, mainly for
!           improved transition from thin cirrus to thick, precipitating ice
!           anvils.
!        -> 8/22/01: This feature has been diminished by effectively adjusting to
!           ice saturation during depositional growth at temperatures colder than
!           -10C.  Ice sublimation is calculated more explicitly.  The logic is
!           that sources of are either poorly understood (e.g., nucleation for NWP)
!           or are not represented in the Eta model (e.g., detrainment of ice from
!           convection).  Otherwise the model is too wet compared to the radiosonde
!           observations based on 1 Feb - 18 March 2001 retrospective runs.
!       (7) Top-down integration also attempts to treat mixed-phase processes,
!           allowing a mixture of ice and water.  Based on numerous observational
!           studies, ice growth is based on nucleation at cloud top &
!           subsequent growth by vapor deposition and riming as the ice particles
!           fall through the cloud.  Effective nucleation rates are a function
!           of ice supersaturation following Meyers et al. (JAM, 1992).
!        -> 8/22/01: The simulated relative humidities were far too moist compared
!           to the rawinsonde observations.  This feature has been substantially
!           diminished, limited to a much narrower temperature range of 0 to -10C.
!       (8) Depositional growth of newly nucleated ice is calculated for large time
!           steps using Fig. 8 of Miller and Young (JAS, 1979), at 1 deg intervals
!           using their ice crystal masses calculated after 600 s of growth in water
!           saturated conditions.  The growth rates are normalized by time step
!           assuming 3D growth with time**1.5 following eq. (6.3) in Young (1993).
!        -> 8/22/01: This feature has been effectively limited to 0 to -10C.
!       (9) Ice precipitation rates can increase due to increase in response to
!           cloud water riming due to (a) increased density & mass of the rimed
!           ice, and (b) increased fall speeds of rimed ice.
!        -> 8/22/01: This feature has been effectively limited to 0 to -10C.
!###############################################################################
!###############################################################################

    SUBROUTINE GSMCOLUMN ( ARAIN, ASNOW, DTPH, I_index, J_index, LSFC, &
    P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col, &
    THICK_col, WC_col )

!###############################################################################
!###############################################################################

!-------------------------------------------------------------------------------
!----- NOTE:  Code is currently set up w/o threading!
!-------------------------------------------------------------------------------
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .
! SUBPROGRAM:  Grid-scale microphysical processes - condensation & precipitation
!   PRGRMMR: Ferrier         ORG: W/NP22     DATE: 08-2001
!-------------------------------------------------------------------------------
! ABSTRACT:
!   * Merges original GSCOND & PRECPD subroutines.
!   * Code has been substantially streamlined and restructured.
!   * Exchange between water vapor & small cloud condensate is calculated using
!     the original Asai (1965, J. Japan) algorithm.  See also references to
!     Yau and Austin (1979, JAS), Rutledge and Hobbs (1983, JAS), and Tao et al.
!     (1989, MWR).  This algorithm replaces the Sundqvist et al. (1989, MWR)
!     parameterization.
!-------------------------------------------------------------------------------

! USAGE:
!   * CALL GSMCOLUMN FROM SUBROUTINE GSMDRIVE
!   * SUBROUTINE GSMDRIVE CALLED FROM MAIN PROGRAM EBU

! INPUT ARGUMENT LIST:
!   DTPH       - physics time step (s)
!   I_index    - I index
!   J_index    - J index
!   LSFC       - Eta level of level above surface, ground
!   P_col      - vertical column of model pressure (Pa)
!   QI_col     - vertical column of model ice mixing ratio (kg/kg)
!   QR_col     - vertical column of model rain ratio (kg/kg)
!   QV_col     - vertical column of model water vapor specific humidity (kg/kg)
!   QW_col     - vertical column of model cloud water mixing ratio (kg/kg)
!   RimeF_col  - vertical column of rime factor for ice in model (ratio, defined below)
!   T_col      - vertical column of model temperature (deg K)
!   THICK_col  - vertical column of model mass thickness (density*height increment)
!   WC_col     - vertical column of model mixing ratio of total condensate (kg/kg)


! OUTPUT ARGUMENT LIST:
!   ARAIN      - accumulated rainfall at the surface (kg)
!   ASNOW      - accumulated snowfall at the surface (kg)
!   QV_col     - vertical column of model water vapor specific humidity (kg/kg)
!   WC_col     - vertical column of model mixing ratio of total condensate (kg/kg)
!   QW_col     - vertical column of model cloud water mixing ratio (kg/kg)
!   QI_col     - vertical column of model ice mixing ratio (kg/kg)
!   QR_col     - vertical column of model rain ratio (kg/kg)
!   RimeF_col  - vertical column of rime factor for ice in model (ratio, defined below)
!   T_col      - vertical column of model temperature (deg K)

! OUTPUT FILES:
!     NONE

! Subprograms & Functions called:
!   * Real Function CONDENSE  - cloud water condensation
!   * Real Function DEPOSIT   - ice deposition (not sublimation)

! UNIQUE: NONE

! LIBRARY: NONE

! COMMON BLOCKS:
!     CMICRO_CONS  - key constants initialized in GSMCONST
!     CMICRO_STATS - accumulated and maximum statistics
!     CMY_GROWTH   - lookup table for growth of ice crystals in
!                    water saturated conditions (Miller & Young, 1979)
!     IVENT_TABLES - lookup tables for ventilation effects of ice
!     IACCR_TABLES - lookup tables for accretion rates of ice
!     IMASS_TABLES - lookup tables for mass content of ice
!     IRATE_TABLES - lookup tables for precipitation rates of ice
!     IRIME_TABLES - lookup tables for increase in fall speed of rimed ice
!     RVENT_TABLES - lookup tables for ventilation effects of rain
!     RACCR_TABLES - lookup tables for accretion rates of rain
!     RMASS_TABLES - lookup tables for mass content of rain
!     RVELR_TABLES - lookup tables for fall speeds of rain
!     RRATE_TABLES - lookup tables for precipitation rates of rain

! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP


!-------------------------------------------------------------------------
!--------------- Arrays & constants in argument list ---------------------
!-------------------------------------------------------------------------

    INCLUDE "parmeta.f90"
    INCLUDE "mpp.h"
    REAL :: ARAIN, ASNOW, P_col(LM), QI_col(LM), QR_col(LM), QV_col(LM), &
    QW_col(LM), RimeF_col(LM), T_col(LM), THICK_col(LM), &
    WC_col(LM)

    INTEGER :: I_index, J_index, LSFC

!-------------------------------------------------------------------------
!-------------- Common blocks for microphysical statistics ---------------
!-------------------------------------------------------------------------

!-------------------------------------------------------------------------
!--------- Common blocks for constants initialized in GSMCONST ----------

    COMMON /CMICRO_CONS/ ABFR, CBFR, CIACW, CIACR, C_N0r0, &
    CN0r0, CN0r_DMRmin, CN0r_DMRmax, CRACW, CRAUT, ESW0, &
    QAUT0, RFmax, RHgrd, RQR_DR1, RQR_DR2, RQR_DR3, RQR_DRmin, &
    RQR_DRmax, RR_DRmin, RR_DR1, RR_DR2, RR_DR3, RR_DRmax

!--- The following variables are for microphysical statistics

    INTEGER, PARAMETER :: ITLO=-60, ITHI=40
    COMMON /CMICRO_STATS/ NSTATS(ITLO:ITHI,4), QMAX(ITLO:ITHI,5), &
    QTOT(ITLO:ITHI,22)
    INTEGER :: NSTATS
    REAL :: QMAX, QTOT

!-------------------------------------------------------------------------
!--------------- Common blocks for various lookup tables -----------------

!--- Discretized growth rates of small ice crystals after their nucleation
!     at 1 C intervals from -1 C to -35 C, based on calculations by Miller
!     and Young (1979, JAS) after 600 s of growth.  Resultant growth rates
!     are multiplied by physics time step in GSMCONST.

    INTEGER, PARAMETER :: MY_T1=1, MY_T2=35
    COMMON /CMY600/ MY_GROWTH(MY_T1:MY_T2)
    REAL :: MY_GROWTH

!-------------------------------------------------------------------------

!--- Mean ice-particle diameters varying from 50 microns to 1000 microns
!      (1 mm), assuming an exponential size distribution.

!---- Meaning of the following arrays:
!        - mdiam - mean diameter (m)
!        - VENTI1 - integrated quantity associated w/ ventilation effects
!                   (capacitance only) for calculating vapor deposition onto ice
!        - VENTI2 - integrated quantity associated w/ ventilation effects
!                   (with fall speed) for calculating vapor deposition onto ice
!        - ACCRI  - integrated quantity associated w/ cloud water collection by ice
!        - MASSI  - integrated quantity associated w/ ice mass
!        - VSNOWI - mass-weighted fall speed of snow (large ice), used to calculate
!                   precipitation rates

    REAL, PARAMETER :: DMImin=.05e-3, DMImax=1.e-3, DelDMI=1.e-6, &
    XMImin=1.e6*DMImin, XMImax=1.e6*DMImax
    INTEGER, PARAMETER :: MDImin=XMImin, MDImax=XMImax

    COMMON /IACCR_TABLES/ ACCRI(MDImin:MDImax)
    COMMON /IMASS_TABLES/ MASSI(MDImin:MDImax)
    REAL :: MASSI
    COMMON /IRATE_TABLES/ VSNOWI(MDImin:MDImax)
    COMMON /IVENT_TABLES/ VENTI1(MDImin:MDImax), VENTI2(MDImin:MDImax)

!-------------------------------------------------------------------------

!--- VEL_RF - velocity increase of rimed particles as functions of crude
!      particle size categories (at 0.1 mm intervals of mean ice particle
!      sizes) and rime factor (different values of Rime Factor of 1.1**N,
!      where N=0 to Nrime).

    INTEGER, PARAMETER :: Nrime=40
    COMMON /IRIME_TABLES/ VEL_RF(2:9,0:Nrime)

!-------------------------------------------------------------------------

!--- Mean rain drop diameters varying from 50 microns (0.05 mm) to 450 microns
!      (0.45 mm), assuming an exponential size distribution.

    REAL, PARAMETER :: DMRmin=.05e-3, DMRmax=.45e-3, DelDMR=1.e-6, &
    XMRmin=1.e6*DMRmin, XMRmax=1.e6*DMRmax
    INTEGER, PARAMETER :: MDRmin=XMRmin, MDRmax=XMRmax

    COMMON /RACCR_TABLES/ ACCRR(MDRmin:MDRmax)
    COMMON /RMASS_TABLES/ MASSR(MDRmin:MDRmax)
    REAL :: MASSR
    COMMON /RRATE_TABLES/ RRATE(MDRmin:MDRmax)
    COMMON /RVELR_TABLES/ VRAIN(MDRmin:MDRmax)
    COMMON /RVENT_TABLES/ VENTR1(MDRmin:MDRmax), VENTR2(MDRmin:MDRmax)

!-------------------------------------------------------------------------
!------- Key parameters, local variables, & important comments ---------
!-----------------------------------------------------------------------

!--- KEY Parameters:

!---- Comments on 2 August 2000
!    * EPSQ=1.E-20 is the lower limit for specific humidity and cloud
!      condensate.  The value of EPSQ will need to be changed in the other
!      subroutines in order to make it consistent throughout the Eta code.
!- NLImax - maximum number concentration of large ice crystals (20,000 /m**3, 20 per liter)
!- NLImin - minimum number concentration of large ice crystals (100 /m**3, 0.1 per liter)

    REAL, PARAMETER :: CP=1004.6, EPSQ=1.E-20, EPSW=1.E-12, RD=287.04, &
    RHOL=1000., RV=461.5, T0C=273.15, XLS=2.834E6, EPS=RD/RV, &
    NLImax=20.E3, NLImin=100., T_ICE=-10., T_ICE_init=-5., &
    TOLER=5.E-7, &
    CLIMIT=10.*EPSQ, CLIMIT1=-CLIMIT, EPS1=RV/RD-1., RCP=1./CP, &
    RCPRV=RCP/RV, RRHOL=1./RHOL, XLS1=XLS*RCP, XLS2=XLS*XLS*RCPRV, &
    XLS3=XLS*XLS/RV, &
    C1=1./3., C2=1./6., C3=3.31/6., &
    DMR1=.1E-3, DMR2=.2E-3, DMR3=.32E-3, N0r0=8.E6, N0rmin=1.e4, &
    N0s0=4.E6, RHO0=1.194, XMR1=1.e6*DMR1, XMR2=1.e6*DMR2, &
    XMR3=1.e6*DMR3, Xratio=.025
    INTEGER, PARAMETER :: MDR1=XMR1, MDR2=XMR2, MDR3=XMR3

!--- If BLEND=1:
!      precipitation (large) ice amounts are estimated at each level as a
!      blend of ice falling from the grid point above and the precip ice
!      present at the start of the time step (see TOT_ICE below).
!--- If BLEND=0:
!      precipitation (large) ice amounts are estimated to be the precip
!      ice present at the start of the time step.

!--- Extended to include sedimentation of rain on 2/5/01

    REAL, PARAMETER :: BLEND=1.

!--- This variable is for debugging purposes (if .true.)

    LOGICAL, PARAMETER :: PRINT_diag= .TRUE. 

!--- Local variables

    REAL :: EMAIRI, N0r, NLICE, NSmICE
    LOGICAL :: CLEAR, ICE_logical, DBG_logical, DIAG_print, RAIN_logical


!#######################################################################
!########################## Begin Execution ############################
!#######################################################################


    ARAIN=0.                ! Accumulated rainfall into grid box from above (kg/m**2)
    ASNOW=0.                ! Accumulated snowfall into grid box from above (kg/m**2)
    DIAG_print= .FALSE. 

!-----------------------------------------------------------------------
!------------ Loop from top (L=1) to surface (L=LSFC) ------------------
!-----------------------------------------------------------------------

    DO 10 L=1,LSFC

    !--- Skip this level and go to the next lower level if no condensate
    !      and very low specific humidities
    
        IF (QV_col(L) <= EPSQ .AND. WC_col(L) <= EPSQ) GO TO 10
    
    !-----------------------------------------------------------------------
    !------------ Proceed with cloud microphysics calculations -------------
    !-----------------------------------------------------------------------
    
        TK=T_col(L)         ! Temperature (deg K)
        TC=TK-T0C           ! Temperature (deg C)
        PP=P_col(L)         ! Pressure (Pa)
        QV=QV_col(L)        ! Specific humidity of water vapor (kg/kg)
        WV=QV/(1.-QV)       ! Water vapor mixing ratio (kg/kg)
        WC=WC_col(L)        ! Grid-scale mixing ratio of total condensate (water or ice; kg/kg)
    
    !-----------------------------------------------------------------------
    !--- Moisture variables below are mixing ratios & not specifc humidities
    !-----------------------------------------------------------------------
    
        CLEAR= .TRUE. 
    
    !--- This check is to determine grid-scale saturation when no condensate is present
    
        ESW=1000.*FPVS0(TK)              ! Saturation vapor pressure w/r/t water
        QSW=EPS*ESW/(PP-ESW)             ! Saturation mixing ratio w/r/t water
        WS=QSW                           ! General saturation mixing ratio (water/ice)
        IF (TC < 0.) THEN
            ESI=1000.*FPVS(TK)             ! Saturation vapor pressure w/r/t ice
            QSI=EPS*ESI/(PP-ESI)           ! Saturation mixing ratio w/r/t water
            WS=QSI                         ! General saturation mixing ratio (water/ice)
        ENDIF
    
    !--- Effective grid-scale Saturation mixing ratios
    
        QSWgrd=RHgrd*QSW
        QSIgrd=RHgrd*QSI
        WSgrd=RHgrd*WS
    
    !--- Check if air is subsaturated and w/o condensate
    
        IF (WV > WSgrd .OR. WC > EPSW) CLEAR= .FALSE. 
    
    !--- Check if any rain is falling into layer from above
    
        IF (ARAIN > CLIMIT) THEN
            CLEAR= .FALSE. 
        ELSE
            ARAIN=0.
        ENDIF
    
    !--- Check if any ice is falling into layer from above
    
    !--- NOTE that "SNOW" in variable names is synonomous with
    !    large, precipitation ice particles
    
        IF (ASNOW > CLIMIT) THEN
            CLEAR= .FALSE. 
        ELSE
            ASNOW=0.
        ENDIF
    
    !-----------------------------------------------------------------------
    !-- Loop to the end if in clear, subsaturated air free of condensate ---
    !-----------------------------------------------------------------------
    
        IF (CLEAR) GO TO 10
    
    !-----------------------------------------------------------------------
    !--------- Initialize RHO, THICK & microphysical processes -------------
    !-----------------------------------------------------------------------
    
    
    !--- Virtual temperature, TV=T*(1./EPS-1)*Q, Q is specific humidity;
    !    (see pp. 63-65 in Fleagle & Businger, 1963)
    
        RHO=PP/(RD*TK*(1.+EPS1*QV))   ! Air density (kg/m**3)
        RRHO=1./RHO                ! Reciprocal of air density
        DTRHO=DTPH*RHO             ! Time step * air density
        BLDTRH=BLEND*DTRHO         ! Blend parameter * time step * air density
        THICK=THICK_col(L)         ! Layer thickness = RHO*DZ = -DP/G = (Psfc-Ptop)*D_ETA/(G*ETA_sfc)
    
        ARAINnew=0.                ! Updated accumulated rainfall
        ASNOWnew=0.                ! Updated accumulated snowfall
        QI=QI_col(L)               ! Ice mixing ratio
        QInew=0.                   ! Updated ice mixing ratio
        QR=QR_col(L)               ! Rain mixing ratio
        QRnew=0.                   ! Updated rain ratio
        QW=QW_col(L)               ! Cloud water mixing ratio
        QWnew=0.                   ! Updated cloud water ratio
    
        PCOND=0.                   ! Condensation (>0) or evaporation (<0) of cloud water (kg/kg)
        PIDEP=0.                   ! Deposition (>0) or sublimation (<0) of ice crystals (kg/kg)
        PIACW=0.                   ! Cloud water collection (riming) by precipitation ice (kg/kg; >0)
        PIACWI=0.                  ! Growth of precip ice by riming (kg/kg; >0)
        PIACWR=0.                  ! Shedding of accreted cloud water to form rain (kg/kg; >0)
        PIACR=0.                   ! Freezing of rain onto large ice at supercooled temps (kg/kg; >0)
        PICND=0.                   ! Condensation (>0) onto wet, melting ice (kg/kg)
        PIEVP=0.                   ! Evaporation (<0) from wet, melting ice (kg/kg)
        PIMLT=0.                   ! Melting ice (kg/kg; >0)
        PRAUT=0.                   ! Cloud water autoconversion to rain (kg/kg; >0)
        PRACW=0.                   ! Cloud water collection (accretion) by rain (kg/kg; >0)
        PREVP=0.                   ! Rain evaporation (kg/kg; <0)
    
    !--- Double check input hydrometeor mixing ratios
    
    !          DUM=WC-(QI+QW+QR)
    !          DUM1=ABS(DUM)
    !          DUM2=TOLER*MIN(WC, QI+QW+QR)
    !          IF (DUM1 .GT. DUM2) THEN
    !            WRITE(6,"(/2(a,i4),a,i2)") '{@ i=',I_index,' j=',J_index,
    !     &                                 ' L=',L
    !            WRITE(6,"(4(a12,g11.4,1x))")
    !     & '{@ TCold=',TC,'P=',.01*PP,'DIFF=',DUM,'WCold=',WC,
    !     & '{@ QIold=',QI,'QWold=',QW,'QRold=',QR
    !          ENDIF
    
    !***********************************************************************
    !*********** MAIN MICROPHYSICS CALCULATIONS NOW FOLLOW! ****************
    !***********************************************************************
    
    !--- Calculate a few variables, which are used more than once below
    
    !--- Latent heat of vaporization as a function of temperature from
    !      Bolton (1980, JAS)
    
        XLV=3.148E6-2370*TK        ! Latent heat of vaporization (Lv)
        XLF=XLS-XLV                ! Latent heat of fusion (Lf)
        XLV1=XLV*RCP               ! Lv/Cp
        XLF1=XLF*RCP               ! Lf/Cp
        TK2=1./(TK*TK)             ! 1./TK**2
        XLV2=XLV*XLV*QSW*TK2/RV    ! Lv**2*Qsw/(Rv*TK**2)
        DENOMW=1.+XLV2*RCP         ! Denominator term, Clausius-Clapeyron correction
    
    !--- Basic thermodynamic quantities
    !      * DYNVIS - dynamic viscosity  [ kg/(m*s) ]
    !      * THERM_COND - thermal conductivity  [ J/(m*s*K) ]
    !      * DIFFUS - diffusivity of water vapor  [ m**2/s ]
    
        TFACTOR=TK**1.5/(TK+120.)
        DYNVIS=1.496E-6*TFACTOR
        THERM_COND=2.116E-3*TFACTOR
        DIFFUS=8.794E-5*TK**1.81/PP
    
    !--- Air resistance term for the fall speed of ice following the
    !      basic research by Heymsfield, Kajikawa, others
    
        GAMMAS=(1.E5/PP)**C1
    
    !--- Air resistance for rain fall speed (Beard, 1985, JAS, p.470)
    
        GAMMAR=(RHO0/RHO)**.4
    
    !----------------------------------------------------------------------
    !-------------  IMPORTANT MICROPHYSICS DECISION TREE  -----------------
    !----------------------------------------------------------------------
    
    !--- Determine if conditions supporting ice are present
    
        IF (TC < 0. .OR. QI > CLIMIT .OR. ASNOW > CLIMIT) THEN
            ICE_logical= .TRUE. 
        ELSE
            ICE_logical= .FALSE. 
            QLICE=0.
            QTICE=0.
        ENDIF
    
    !--- Determine if rain is present
    
        RAIN_logical= .FALSE. 
        IF (ARAIN > CLIMIT .OR. QR > CLIMIT) RAIN_logical= .TRUE. 
    
        IF (ICE_logical) THEN
        
        !--- IMPORTANT:  Estimate time-averaged properties.
        
        !---
        !  * FLARGE  - ratio of number of large ice to total (large & small) ice
        !  * FSMALL  - ratio of number of small ice crystals to large ice particles
        !  ->  Small ice particles are assumed to have a mean diameter of 50 microns.
        !  * XSIMASS - used for calculating small ice mixing ratio
        !---
        !  * TOT_ICE - total mass (small & large) ice before microphysics,
        !              which is the sum of the total mass of large ice in the
        !              current layer and the input flux of ice from above
        !  * PILOSS  - greatest loss (<0) of total (small & large) ice by
        !              sublimation, removing all of the ice falling from above
        !              and the ice within the layer
        !  * RimeF1  - Rime Factor, which is the mass ratio of total (unrimed & rimed)
        !              ice mass to the unrimed ice mass (>=1)
        !  * VrimeF  - the velocity increase due to rime factor or melting (ratio, >=1)
        !  * VSNOW   - Fall speed of rimed snow w/ air resistance correction
        !  * EMAIRI  - equivalent mass of air associated layer and with fall of snow into layer
        !  * XLIMASS - used for calculating large ice mixing ratio
        !  * FLIMASS - mass fraction of large ice
        !  * QTICE   - time-averaged mixing ratio of total ice
        !  * QLICE   - time-averaged mixing ratio of large ice
        !  * NLICE   - time-averaged number concentration of large ice
        !  * NSmICE  - number concentration of small ice crystals at current level
        !---
        !--- Assumed number fraction of large ice particles to total (large & small)
        !    ice particles, which is based on a general impression of the literature.
        
            WVQW=WV+QW                ! Water vapor & cloud water
        
            IF (TC >= 0. .OR. WVQW < QSIgrd) THEN
            
            !--- Eliminate small ice particle contributions for melting & sublimation
            
                FLARGE=1.
            ELSE
            
            !--- Enhanced number of small ice particles during depositional growth
            !    (effective only when 0C > T >= T_ice [-10C] )
            
                FLARGE=.2
            
            !--- Larger number of small ice particles due to rime splintering
            
                IF (TC >= -8. .AND. TC <= -3.) FLARGE=.5*FLARGE
            
            ENDIF            ! End IF (TC >= 0. .OR. WVQW < QSIgrd)
            FSMALL=(1.-FLARGE)/FLARGE
            XSIMASS=RRHO*MASSI(MDImin)*FSMALL
            IF (QI <= CLIMIT .AND. ASNOW <= CLIMIT) THEN
                INDEXS=MDImin
                TOT_ICE=0.
                PILOSS=0.
                RimeF1=1.
                VrimeF=1.
                VEL_INC=GAMMAS
                VSNOW=0.
                EMAIRI=THICK
                XLIMASS=RRHO*RimeF1*MASSI(INDEXS)
                FLIMASS=XLIMASS/(XLIMASS+XSIMASS)
                QLICE=0.
                QTICE=0.
                NLICE=0.
                NSmICE=0.
            ELSE
            
            !--- For T<0C mean particle size follows Houze et al. (JAS, 1979, p. 160),
            !    converted from Fig. 5 plot of LAMDAs.  Similar set of relationships
            !    also shown in Fig. 8 of Ryan (BAMS, 1996, p. 66).
            
                DUM=XMImax*EXP(.0536*TC)
                INDEXS=MIN(MDImax, MAX(MDImin, INT(DUM) ) )
                TOT_ICE=THICK*QI+BLEND*ASNOW
                PILOSS=-TOT_ICE/THICK
                LBEF=MAX(1,L-1)
                DUM1=RimeF_col(LBEF)
                DUM2=RimeF_col(L)
                RimeF1=(DUM2*THICK*QI+DUM1*BLEND*ASNOW)/TOT_ICE
                RimeF1=MIN(RimeF1, RFmax)
                DO IPASS=0,1
                    IF (RimeF1 <= 1.) THEN
                        RimeF1=1.
                        VrimeF=1.
                    ELSE
                        IXS=MAX(2, MIN(INDEXS/100, 9))
                        XRF=10.492*ALOG(RimeF1)
                        IXRF=MAX(0, MIN(INT(XRF), Nrime))
                        IF (IXRF >= Nrime) THEN
                            VrimeF=VEL_RF(IXS,Nrime)
                        ELSE
                            VrimeF=VEL_RF(IXS,IXRF)+(XRF-FLOAT(IXRF))* &
                            (VEL_RF(IXS,IXRF+1)-VEL_RF(IXS,IXRF))
                        ENDIF
                    ENDIF            ! End IF (RimeF1 <= 1.)
                    VEL_INC=GAMMAS*VrimeF
                    VSNOW=VEL_INC*VSNOWI(INDEXS)
                    EMAIRI=THICK+BLDTRH*VSNOW
                    XLIMASS=RRHO*RimeF1*MASSI(INDEXS)
                    FLIMASS=XLIMASS/(XLIMASS+XSIMASS)
                    QTICE=TOT_ICE/EMAIRI
                    QLICE=FLIMASS*QTICE
                    NLICE=QLICE/XLIMASS
                    NSmICE=Fsmall*NLICE
                
                    IF ( (NLICE >= NLImin .AND. NLICE <= NLImax) &
                     .OR. IPASS == 1) THEN
                        EXIT
                    ELSE
                    
                    !--- Reduce excessive accumulation of ice at upper levels
                    !    associated with strong grid-resolved ascent
                    
                    !--- Force NLICE to be between NLImin and NLImax
                    
                        DUM=MAX(NLImin, MIN(NLImax, NLICE) )
                        XLI=RHO*(QTICE/DUM-XSIMASS)/RimeF1
                        IF (XLI <= MASSI(MDImin) ) THEN
                            INDEXS=MDImin
                        ELSE IF (XLI <= MASSI(450) ) THEN
                            DLI=9.5885E5*XLI**.42066         ! DLI in microns
                            INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
                        ELSE IF (XLI <= MASSI(MDImax) ) THEN
                            DLI=3.9751E6*XLI**.49870         ! DLI in microns
                            INDEXS=MIN(MDImax, MAX(MDImin, INT(DLI) ) )
                        ELSE
                            INDEXS=MDImax
                        
                        !--- 8/22/01: Increase density of large ice if maximum limits
                        !    are reached for number concentration (NLImax) and mean size
                        !    (MDImax).  Done to increase fall out of ice.
                        
                            IF (DUM >= NLImax) &
                            RimeF1=RHO*(QTICE/NLImax-XSIMASS)/MASSI(INDEXS)
                        ENDIF             ! End IF (XLI <= MASSI(MDImin) )
                    !            WRITE(6,"(4(a12,g11.4,1x))")
                    !     & '{$ TC=',TC,'P=',.01*PP,'NLICE=',NLICE,'DUM=',DUM,
                    !     & '{$ XLI=',XLI,'INDEXS=',FLOAT(INDEXS),'RHO=',RHO,'QTICE=',QTICE,
                    !     & '{$ XSIMASS=',XSIMASS,'RimeF1=',RimeF1
                    ENDIF                  ! End IF ( (NLICE >= NLImin .AND. NLICE <= NLImax) ...
                ENDDO                    ! End DO IPASS=0,1
            ENDIF                      ! End IF (QI <= CLIMIT .AND. ASNOW <= CLIMIT)
        ENDIF                        ! End IF (ICE_logical)
    
    !----------------------------------------------------------------------
    !--------------- Calculate individual processes -----------------------
    !----------------------------------------------------------------------
    
    !--- Cloud water autoconversion to rain and collection by rain
    
        IF (QW > CLIMIT .AND. TC >= T_ICE) THEN
        
        !--- QW0 could be modified based on land/sea properties,
        !      presence of convection, etc.  This is why QAUT0 and CRAUT
        !      are passed into the subroutine as externally determined
        !      parameters.  Can be changed in the future if desired.
        
            QW0=QAUT0*RRHO
            PRAUT=MAX(0., QW-QW0)*CRAUT
            IF (QLICE > CLIMIT) THEN
            
            !--- Collection of cloud water by large ice particles ("snow")
            !    PIACWI=PIACW for riming, PIACWI=0 for shedding
            
                FWS=MIN(1., CIACW*VEL_INC*NLICE*ACCRI(INDEXS)/PP**C1)
                PIACW=FWS*QW
                IF (TC < 0.) PIACWI=PIACW    ! Large ice riming
            ENDIF           ! End IF (QLICE > CLIMIT)
        ENDIF             ! End IF (QW > CLIMIT .AND. TC >= T_ICE)
    
    !----------------------------------------------------------------------
    !--- Loop around some of the ice-phase processes if no ice should be present
    !----------------------------------------------------------------------
    
        IF (ICE_logical .EQV. .FALSE. ) GO TO 20
    
    !--- Now the pretzel logic of calculating ice deposition
    
        IF (TC < T_ICE .AND. (WV > QSIgrd .OR. QW > CLIMIT)) THEN
        
        !--- Adjust to ice saturation at T<T_ICE (-10C) if supersaturated.
        !    Sources of ice due to nucleation and convective detrainment are
        !    either poorly understood, poorly resolved at typical NWP
        !    resolutions, or are not represented (e.g., no detrained
        !    condensate in BMJ Cu scheme).
        
            PCOND=-QW
            DUM1=TK+XLV1*PCOND                 ! Updated (dummy) temperature (deg K)
            DUM2=WV+QW                         ! Updated (dummy) water vapor mixing ratio
            DUM=1000.*FPVS(DUM1)               ! Updated (dummy) saturation vapor pressure w/r/t ice
            DUM=RHgrd*EPS*DUM/(PP-DUM)         ! Updated (dummy) saturation mixing ratio w/r/t ice
            IF (DUM2 > DUM) PIDEP=DEPOSIT (PP, RHgrd, DUM1, DUM2)
            DWVi=0.    ! Used only for debugging
        
        ELSE IF (TC < 0.) THEN
        
        !--- These quantities are handy for ice deposition/sublimation
        !    PIDEP_max - max deposition or minimum sublimation to ice saturation
        
            DENOMI=1.+XLS2*QSI*TK2
            DWVi=MIN(WVQW,QSW)-QSI
            PIDEP_max=MAX(PILOSS, DWVi/DENOMI)
            IF (QTICE > 0.) THEN
            
            !--- Calculate ice deposition/sublimation
            !      * SFACTOR - [VEL_INC**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5],
            !        where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
            !      * Units: SFACTOR - s**.5/m ;  ABI - m**2/s ;  NLICE - m**-3 ;
            !               VENTIL, VENTIS - m**-2 ;  VENTI1 - m ;
            !               VENTI2 - m**2/s**.5 ; DIDEP - unitless
            
                SFACTOR=VEL_INC**.5*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2
                ABI=1./(RHO*XLS3*QSI*TK2/THERM_COND+1./DIFFUS)
            
            !--- VENTIL - Number concentration * ventilation factors for large ice
            !--- VENTIS - Number concentration * ventilation factors for small ice
            
            !--- Variation in the number concentration of ice with time is not
            !      accounted for in these calculations (could be in the future).
            
                VENTIL=(VENTI1(INDEXS)+SFACTOR*VENTI2(INDEXS))*NLICE
                VENTIS=(VENTI1(MDImin)+SFACTOR*VENTI2(MDImin))*NSmICE
                DIDEP=ABI*(VENTIL+VENTIS)*DTPH
            
            !--- Account for change in water vapor supply w/ time
            
                IF (DIDEP >= Xratio) &
                DIDEP=(1.-EXP(-DIDEP*DENOMI))/DENOMI
                IF (DWVi > 0.) THEN
                    PIDEP=MIN(DWVi*DIDEP, PIDEP_max)
                ELSE IF (DWVi < 0.) THEN
                    PIDEP=MAX(DWVi*DIDEP, PIDEP_max)
                ENDIF
            
            ELSE IF (WVQW > QSI .AND. TC <= T_ICE_init) THEN
            
            !--- Ice nucleation in near water-saturated conditions.  Ice crystal
            !    growth during time step calculated using Miller & Young (1979, JAS).
            !--- These deposition rates could drive conditions below water saturation,
            !    which is the basis of these calculations.  Intended to approximate
            !    more complex & computationally intensive calculations.
            
                INDEX_MY=MAX(MY_T1, MIN( INT(.5-TC), MY_T2 ) )
            
            !--- DUM1 is the supersaturation w/r/t ice at water-saturated conditions
            
            !--- DUM2 is the number of ice crystals nucleated at water-saturated
            !    conditions based on Meyers et al. (JAM, 1992).
            
            !--- Prevent unrealistically large ice initiation (limited by PIDEP_max)
            !      if DUM2 values are increased in future experiments
            
                DUM1=QSW/QSI-1.
                DUM2=1.E3*EXP(12.96*DUM1-.039)
                PIDEP=MIN(PIDEP_max, DUM2*MY_GROWTH(INDEX_MY)*RRHO)
            
            ENDIF       ! End IF (QTICE > 0.)
        
        ENDIF         ! End IF (TC < T_ICE .AND. (WV > QSIgrd .OR. QW > CLIMIT))
    
    !------------------------------------------------------------------------
    
        20 CONTINUE     ! Jump here if conditions for ice are not present
    
    !------------------------------------------------------------------------
    
    !--- Cloud water condensation
    
        IF (TC >= T_ICE .AND. (QW > CLIMIT .OR. WV > QSWgrd)) THEN
            IF (PIACWI == 0. .AND. PIDEP == 0.) THEN
                PCOND=CONDENSE (PP, QW, RHgrd, TK, WV)
            ELSE
            
            !--- Modify cloud condensation in response to ice processes
            
                DUM=XLV*QSWgrd*RCPRV*TK2
                DENOMWI=1.+XLS*DUM
                DENOMF=XLF*DUM
                DUM=MAX(0., PIDEP)
                PCOND=(WV-QSWgrd-DENOMWI*DUM-DENOMF*PIACWI)/DENOMW
                DUM1=-QW
                DUM2=PCOND-PIACW
                IF (DUM2 < DUM1) THEN
                
                !--- Limit cloud water sinks
                
                    DUM=DUM1/DUM2
                    PCOND=DUM*PCOND
                    PIACW=DUM*PIACW
                    PIACWI=DUM*PIACWI
                ENDIF        ! End IF (DUM2 < DUM1)
            ENDIF          ! End IF (PIACWI == 0. .AND. PIDEP == 0.)
        ENDIF            ! End IF (TC >= T_ICE .AND. (QW > CLIMIT .OR. WV > QSWgrd))
    
    !--- Limit freezing of accreted rime to prevent temperature oscillations,
    !    a crude Schumann-Ludlam limit (p. 209 of Young, 1993).
    
        TCC=TC+XLV1*PCOND+XLS1*PIDEP+XLF1*PIACWI
        IF (TCC > 0.) THEN
            PIACWI=0.
            TCC=TC+XLV1*PCOND+XLS1*PIDEP
        ENDIF
    
        IF (TC > 0. .AND. TCC > 0. .AND. ICE_logical) THEN
        
        !--- Calculate melting and evaporation/condensation
        !      * Units: SFACTOR - s**.5/m ;  ABI - m**2/s ;  NLICE - m**-3 ;
        !               VENTIL - m**-2 ;  VENTI1 - m ;
        !               VENTI2 - m**2/s**.5 ; CIEVP - /s
        
            SFACTOR=VEL_INC**.5*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2
            VENTIL=NLICE*(VENTI1(INDEXS)+SFACTOR*VENTI2(INDEXS))
            AIEVP=VENTIL*DIFFUS*DTPH
            IF (AIEVP < Xratio) THEN
                DIEVP=AIEVP
            ELSE
                DIEVP=1.-EXP(-AIEVP)
            ENDIF
            QSW0=EPS*ESW0/(PP-ESW0)
            DWV0=MIN(WV,QSW)-QSW0
            DUM=QW+PCOND
            IF (WV < QSW .AND. DUM <= CLIMIT) THEN
            
            !--- Evaporation from melting snow (sink of snow) or shedding
            !    of water condensed onto melting snow (source of rain)
            
                DUM=DWV0*DIEVP
                PIEVP=MAX( MIN(0., DUM), PILOSS)
                PICND=MAX(0., DUM)
            ENDIF            ! End IF (WV < QSW .AND. DUM <= CLIMIT)
            PIMLT=THERM_COND*TCC*VENTIL*RRHO*DTPH/XLF
        
        !--- Limit melting to prevent temperature oscillations across 0C
        
            DUM1=MAX( 0., (TCC+XLV1*PIEVP)/XLF1 )
            PIMLT=MIN(PIMLT, DUM1)
        
        !--- Limit loss of snow by melting (>0) and evaporation
        
            DUM=PIEVP-PIMLT
            IF (DUM < PILOSS) THEN
                DUM1=PILOSS/DUM
                PIMLT=PIMLT*DUM1
                PIEVP=PIEVP*DUM1
            ENDIF           ! End IF (DUM > QTICE)
        ENDIF             ! End IF (TC > 0. .AND. TCC > 0. .AND. ICE_logical)
    
    !--- IMPORTANT:  Estimate time-averaged properties.
    
    !  * TOT_RAIN - total mass of rain before microphysics, which is the sum of
    !               the total mass of rain in the current layer and the input
    !               flux of rain from above
    !  * VRAIN1   - fall speed of rain into grid from above (with air resistance correction)
    !  * QTRAIN   - time-averaged mixing ratio of rain (kg/kg)
    !  * PRLOSS   - greatest loss (<0) of rain, removing all rain falling from
    !               above and the rain within the layer
    !  * RQR      - rain content (kg/m**3)
    !  * INDEXR   - mean size of rain drops to the nearest 1 micron in size
    !  * N0r      - intercept of rain size distribution (typically 10**6 m**-4)
    
        TOT_RAIN=0.
        VRAIN1=0.
        QTRAIN=0.
        PRLOSS=0.
        RQR=0.
        N0r=0.
        INDEXR1=INDEXR    ! For debugging only
        INDEXR=MDRmin
        IF (RAIN_logical) THEN
            IF (ARAIN <= 0.) THEN
                INDEXR=MDRmin
                VRAIN1=0.
            ELSE
            
            !--- INDEXR (related to mean diameter) & N0r could be modified
            !      by land/sea properties, presence of convection, etc.
            
            !--- Rain rate normalized to a density of 1.194 kg/m**3
            
                RR=ARAIN/(DTPH*GAMMAR)
            
                IF (RR <= RR_DRmin) THEN
                
                !--- Assume fixed mean diameter of rain (0.2 mm) for low rain rates,
                !      instead vary N0r with rain rate
                
                    INDEXR=MDRmin
                ELSE IF (RR <= RR_DR1) THEN
                
                !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
                !      for mean diameters (Dr) between 0.05 and 0.10 mm:
                !      V(Dr)=5.6023e4*Dr**1.136, V in m/s and Dr in m
                !      RR = PI*1000.*N0r0*5.6023e4*Dr**(4+1.136) = 1.408e15*Dr**5.136,
                !        RR in kg/(m**2*s)
                !      Dr (m) = 1.123e-3*RR**.1947 -> Dr (microns) = 1.123e3*RR**.1947
                
                    INDEXR=INT( 1.123E3*RR**.1947 + .5 )
                    INDEXR=MAX( MDRmin, MIN(INDEXR, MDR1) )
                ELSE IF (RR <= RR_DR2) THEN
                
                !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
                !      for mean diameters (Dr) between 0.10 and 0.20 mm:
                !      V(Dr)=1.0867e4*Dr**.958, V in m/s and Dr in m
                !      RR = PI*1000.*N0r0*1.0867e4*Dr**(4+.958) = 2.731e14*Dr**4.958,
                !        RR in kg/(m**2*s)
                !      Dr (m) = 1.225e-3*RR**.2017 -> Dr (microns) = 1.225e3*RR**.2017
                
                    INDEXR=INT( 1.225E3*RR**.2017 + .5 )
                    INDEXR=MAX( MDR1, MIN(INDEXR, MDR2) )
                ELSE IF (RR <= RR_DR3) THEN
                
                !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
                !      for mean diameters (Dr) between 0.20 and 0.32 mm:
                !      V(Dr)=2831.*Dr**.80, V in m/s and Dr in m
                !      RR = PI*1000.*N0r0*2831.*Dr**(4+.80) = 7.115e13*Dr**4.80,
                !        RR in kg/(m**2*s)
                !      Dr (m) = 1.3006e-3*RR**.2083 -> Dr (microns) = 1.3006e3*RR**.2083
                
                    INDEXR=INT( 1.3006E3*RR**.2083 + .5 )
                    INDEXR=MAX( MDR2, MIN(INDEXR, MDR3) )
                ELSE IF (RR <= RR_DRmax) THEN
                
                !--- Best fit to mass-weighted fall speeds (V) from rain lookup tables
                !      for mean diameters (Dr) between 0.32 and 0.45 mm:
                !      V(Dr)=944.8*Dr**.6636, V in m/s and Dr in m
                !      RR = PI*1000.*N0r0*944.8*Dr**(4+.6636) = 2.3745e13*Dr**4.6636,
                !        RR in kg/(m**2*s)
                !      Dr (m) = 1.355e-3*RR**.2144 -> Dr (microns) = 1.355e3*RR**.2144
                
                    INDEXR=INT( 1.355E3*RR**.2144 + .5 )
                    INDEXR=MAX( MDR3, MIN(INDEXR, MDRmax) )
                ELSE
                
                !--- Assume fixed mean diameter of rain (0.45 mm) for high rain rates,
                !      instead vary N0r with rain rate
                
                    INDEXR=MDRmax
                ENDIF              ! End IF (RR <= RR_DRmin) etc.
                VRAIN1=GAMMAR*VRAIN(INDEXR)
            ENDIF              ! End IF (ARAIN <= 0.)
            INDEXR1=INDEXR     ! For debugging only
            TOT_RAIN=THICK*QR+BLEND*ARAIN
            QTRAIN=TOT_RAIN/(THICK+BLDTRH*VRAIN1)
            PRLOSS=-TOT_RAIN/THICK
            RQR=RHO*QTRAIN
        
        !--- RQR - time-averaged rain content (kg/m**3)
        
            IF (RQR <= RQR_DRmin) THEN
                N0r=MAX(N0rmin, CN0r_DMRmin*RQR)
                INDEXR=MDRmin
            ELSE IF (RQR >= RQR_DRmax) THEN
                N0r=CN0r_DMRmax*RQR
                INDEXR=MDRmax
            ELSE
                N0r=N0r0
                INDEXR=MAX( XMRmin, MIN(CN0r0*RQR**.25, XMRmax) )
            ENDIF
        
            IF (TC < T_ICE) THEN
                PIACR=-PRLOSS
            ELSE
                DWVr=WV-PCOND-QSW
                DUM=QW+PCOND
                IF (DWVr < 0. .AND. DUM <= CLIMIT) THEN
                
                !--- Rain evaporation
                
                !    * RFACTOR - [GAMMAR**.5]*[Schmidt**(1./3.)]*[(RHO/DYNVIS)**.5],
                !        where Schmidt (Schmidt Number) =DYNVIS/(RHO*DIFFUS)
                
                !    * Units: RFACTOR - s**.5/m ;  ABW - m**2/s ;  VENTR - m**-2 ;
                !             N0r - m**-4 ;  VENTR1 - m**2 ;  VENTR2 - m**3/s**.5 ;
                !             CREVP - unitless
                
                    RFACTOR=GAMMAR**.5*(RHO/(DIFFUS*DIFFUS*DYNVIS))**C2
                    ABW=1./(RHO*XLV2/THERM_COND+1./DIFFUS)
                
                !--- Note that VENTR1, VENTR2 lookup tables do not include the
                !      1/Davg multiplier as in the ice tables
                
                    VENTR=N0r*(VENTR1(INDEXR)+RFACTOR*VENTR2(INDEXR))
                    CREVP=ABW*VENTR*DTPH
                    IF (CREVP < Xratio) THEN
                        DUM=DWVr*CREVP
                    ELSE
                        DUM=DWVr*(1.-EXP(-CREVP*DENOMW))/DENOMW
                    ENDIF
                    PREVP=MAX(DUM, PRLOSS)
                ELSE IF (QW > CLIMIT) THEN
                    FWR=CRACW*GAMMAR*N0r*ACCRR(INDEXR)
                    PRACW=MIN(1.,FWR)*QW
                ENDIF           ! End IF (DWVr < 0. .AND. DUM <= CLIMIT)
            
                IF (TC < 0. .AND. TCC < 0.) THEN
                
                !--- Biggs (1953) heteorogeneous freezing (e.g., Lin et al., 1983)
                
                    PIACR=CBFR*N0r*RRHO*(EXP(ABFR*TC)-1.)*(FLOAT(INDEXR))**7
                    IF (QLICE > CLIMIT) THEN
                    
                    !--- Freezing of rain by collisions w/ large ice
                    
                        DUM=GAMMAR*VRAIN(INDEXR)
                        DUM1=DUM-VSNOW
                    
                    !--- DUM2 - Difference in spectral fall speeds of rain and
                    !      large ice, parameterized following eq. (48) on p. 112 of
                    !      Murakami (J. Meteor. Soc. Japan, 1990)
                    
                        DUM2=(DUM1*DUM1+.04*DUM*VSNOW)**.5
                        DUM1=5.E-12*INDEXR*INDEXR+2.E-12*INDEXR*INDEXS &
                        +.5E-12*INDEXS*INDEXS
                        FIR=MIN(1., CIACR*NLICE*DUM1*DUM2)
                    
                    !--- Future?  Should COLLECTION BY SMALL ICE SHOULD BE INCLUDED???
                    
                        PIACR=MIN(PIACR+FIR*QTRAIN, QTRAIN)
                    ENDIF        ! End IF (QLICE > CLIMIT)
                    DUM=PREVP-PIACR
                    If (DUM < PRLOSS) THEN
                        DUM1=PRLOSS/DUM
                        PREVP=DUM1*PREVP
                        PIACR=DUM1*PIACR
                    ENDIF        ! End If (DUM < PRLOSS)
                ENDIF          ! End IF (TC < 0. .AND. TCC < 0.)
            ENDIF            ! End IF (TC < T_ICE)
        ENDIF              ! End IF (RAIN_logical)
    
    !----------------------------------------------------------------------
    !---------------------- Main Budget Equations -------------------------
    !----------------------------------------------------------------------
    
    
    !-----------------------------------------------------------------------
    !--- Update fields, determine characteristics for next lower layer ----
    !-----------------------------------------------------------------------
    
    !--- Carefully limit sinks of cloud water
    
        DUM1=PIACW+PRAUT+PRACW-MIN(0.,PCOND)
        IF (DUM1 > QW) THEN
            DUM=QW/DUM1
            PIACW=DUM*PIACW
            PIACWI=DUM*PIACWI
            PRAUT=DUM*PRAUT
            PRACW=DUM*PRACW
            IF (PCOND < 0.) PCOND=DUM*PCOND
        ENDIF
        PIACWR=PIACW-PIACWI          ! TC >= 0C
    
    !--- QWnew - updated cloud water mixing ratio
    
        DELW=PCOND-PIACW-PRAUT-PRACW
        QWnew=QW+DELW
        IF (QWnew <= CLIMIT) QWnew=0.
        IF (QW > 0. .AND. QWnew /= 0.) THEN
            DUM=QWnew/QW
            IF (DUM < TOLER) QWnew=0.
        ENDIF
    
    !--- Update temperature and water vapor mixing ratios
    
        DELT= XLV1*(PCOND+PIEVP+PICND+PREVP) &
        +XLS1*PIDEP+XLF1*(PIACWI+PIACR-PIMLT)
        Tnew=TK+DELT
    
        DELV=-PCOND-PIDEP-PIEVP-PICND-PREVP
        WVnew=WV+DELV
    
    !--- Update ice mixing ratios
    
    !---
    !  * TOT_ICEnew - total mass (small & large) ice after microphysics,
    !                 which is the sum of the total mass of large ice in the
    !                 current layer and the flux of ice out of the grid box below
    !  * RimeF      - Rime Factor, which is the mass ratio of total (unrimed &
    !                 rimed) ice mass to the unrimed ice mass (>=1)
    !  * QInew      - updated mixing ratio of total (large & small) ice in layer
    !      -> TOT_ICEnew=QInew*THICK+BLDTRH*QLICEnew*VSNOW
    !        -> But QLICEnew=QInew*FLIMASS, so
    !      -> TOT_ICEnew=QInew*(THICK+BLDTRH*FLIMASS*VSNOW)
    !  * ASNOWnew   - updated accumulation of snow at bottom of grid cell
    !---
    
        DELI=0.
        RimeF=1.
        IF (ICE_logical) THEN
            DELI=PIDEP+PIEVP+PIACWI+PIACR-PIMLT
            TOT_ICEnew=TOT_ICE+THICK*DELI
            IF (TOT_ICE > 0. .AND. TOT_ICEnew /= 0.) THEN
                DUM=TOT_ICEnew/TOT_ICE
                IF (DUM < TOLER) TOT_ICEnew=0.
            ENDIF
            IF (TOT_ICEnew <= CLIMIT) THEN
                TOT_ICEnew=0.
                RimeF=1.
                QInew=0.
                ASNOWnew=0.
            ELSE
            
            !--- Update rime factor if appropriate
            
                DUM=PIACWI+PIACR
                IF (DUM <= CLIMIT .AND. PIDEP <= CLIMIT) THEN
                    RimeF=RimeF1
                ELSE
                
                !--- Rime Factor, RimeF = (Total ice mass)/(Total unrimed ice mass)
                !      DUM1 - Total ice mass, rimed & unrimed
                !      DUM2 - Estimated mass of *unrimed* ice
                
                    DUM1=TOT_ICE+THICK*(PIDEP+DUM)
                    DUM2=TOT_ICE/RimeF1+THICK*PIDEP
                    IF (DUM2 <= 0.) THEN
                        RimeF=RFmax
                    ELSE
                        RimeF=MIN(RFmax, MAX(1., DUM1/DUM2) )
                    ENDIF
                ENDIF       ! End IF (DUM <= CLIMIT .AND. PIDEP <= CLIMIT)
                QInew=TOT_ICEnew/(THICK+BLDTRH*FLIMASS*VSNOW)
                IF (QInew <= CLIMIT) QInew=0.
                IF (QI > 0. .AND. QInew /= 0.) THEN
                    DUM=QInew/QI
                    IF (DUM < TOLER) QInew=0.
                ENDIF
                ASNOWnew=BLDTRH*FLIMASS*VSNOW*QInew
                IF (ASNOW > 0. .AND. ASNOWnew /= 0.) THEN
                    DUM=ASNOWnew/ASNOW
                    IF (DUM < TOLER) ASNOWnew=0.
                ENDIF
            ENDIF         ! End IF (TOT_ICEnew <= CLIMIT)
        ENDIF           ! End IF (ICE_logical)
    
    !--- Update rain mixing ratios
    
    !---
    ! * TOT_RAINnew - total mass of rain after microphysics
    !                 current layer and the input flux of ice from above
    ! * VRAIN2      - time-averaged fall speed of rain in grid and below
    !                 (with air resistance correction)
    ! * QRnew       - updated rain mixing ratio in layer
    !      -> TOT_RAINnew=QRnew*(THICK+BLDTRH*VRAIN2)
    !  * ARAINnew  - updated accumulation of rain at bottom of grid cell
    !---
    
        DELR=PRAUT+PRACW+PIACWR-PIACR+PIMLT+PREVP+PICND
        TOT_RAINnew=TOT_RAIN+THICK*DELR
        IF (TOT_RAIN > 0. .AND. TOT_RAINnew /= 0.) THEN
            DUM=TOT_RAINnew/TOT_RAIN
            IF (DUM < TOLER) TOT_RAINnew=0.
        ENDIF
        IF (TOT_RAINnew <= CLIMIT) THEN
            TOT_RAINnew=0.
            VRAIN2=0.
            QRnew=0.
            ARAINnew=0.
        ELSE
        
        !--- 1st guess time-averaged rain rate at bottom of grid box
        
            RR=TOT_RAINnew/(DTPH*GAMMAR)
        
        !--- Use same algorithm as above for calculating mean drop diameter
        !      (IDR, in microns), which is used to estimate the time-averaged
        !      fall speed of rain drops at the bottom of the grid layer.  This
        !      isnt perfect, but the alternative is solving a transcendental
        !      equation that is numerically inefficient and nasty to program
        !      (coded in earlier versions of GSMCOLUMN prior to 8-22-01).
        
            IF (RR <= RR_DRmin) THEN
                IDR=MDRmin
            ELSE IF (RR <= RR_DR1) THEN
                IDR=INT( 1.123E3*RR**.1947 + .5 )
                IDR=MAX( MDRmin, MIN(IDR, MDR1) )
            ELSE IF (RR <= RR_DR2) THEN
                IDR=INT( 1.225E3*RR**.2017 + .5 )
                IDR=MAX( MDR1, MIN(IDR, MDR2) )
            ELSE IF (RR <= RR_DR3) THEN
                IDR=INT( 1.3006E3*RR**.2083 + .5 )
                IDR=MAX( MDR2, MIN(IDR, MDR3) )
            ELSE IF (RR <= RR_DRmax) THEN
                IDR=INT( 1.355E3*RR**.2144 + .5 )
                IDR=MAX( MDR3, MIN(IDR, MDRmax) )
            ELSE
                IDR=MDRmax
            ENDIF              ! End IF (RR <= RR_DRmin)
            VRAIN2=GAMMAR*VRAIN(IDR)
            QRnew=TOT_RAINnew/(THICK+BLDTRH*VRAIN2)
            IF (QRnew <= CLIMIT) QRnew=0.
            IF (QR > 0. .AND. QRnew /= 0.) THEN
                DUM=QRnew/QR
                IF (DUM < TOLER) QRnew=0.
            ENDIF
            ARAINnew=BLDTRH*VRAIN2*QRnew
            IF (ARAIN > 0. .AND. ARAINnew /= 0.) THEN
                DUM=ARAINnew/ARAIN
                IF (DUM < TOLER) ARAINnew=0.
            ENDIF
        ENDIF
    
        WCnew=QWnew+QRnew+QInew
    
    !----------------------------------------------------------------------
    !-------------- Begin debugging & verification ------------------------
    !----------------------------------------------------------------------
    
    !--- QT, QTnew - total water (vapor & condensate) before & after microphysics, resp.
    
        QT=THICK*(WV+WC)+ARAIN+ASNOW
        QTnew=THICK*(WVnew+WCnew)+ARAINnew+ASNOWnew
        BUDGET=QT-QTnew
    
    !--- Additional check on budget preservation, accounting for truncation effects
    
        DBG_logical= .FALSE. 
    !          DUM=ABS(BUDGET)
    !          IF (DUM .GT. TOLER) THEN
    !            DUM=DUM/MIN(QT, QTnew)
    !            IF (DUM .GT. TOLER) DBG_logical=.TRUE.
    !          ENDIF
    !!
    !          DUM=(RHgrd+.001)*QSInew
    !          IF ( (QWnew.GT.CLIMIT .OR. QRnew.GT.CLIMIT .OR. WVnew.GT.DUM)
    !     &        .AND. TC.LT.T_ICE )  DBG_logical=.TRUE.
    
    !          IF (TC.GT.5. .AND. QInew.GT.CLIMIT) DBG_logical=.TRUE.
    
        IF ((WVnew < CLIMIT .OR. DBG_logical) .AND. PRINT_diag) THEN
        
            WRITE(6,"(/2(a,i4),2(a,i2))") '{} i=',I_index,' j=',J_index, &
            ' L=',L,' LSFC=',LSFC
        
            ESW=1000.*FPVS0(Tnew)
            QSWnew=EPS*ESW/(PP-ESW)
            IF (TC < 0. .OR. Tnew < 0.) THEN
                ESI=1000.*FPVS(Tnew)
                QSInew=EPS*ESI/(PP-ESI)
            ELSE
                QSI=QSW
                QSInew=QSWnew
            ENDIF
            WSnew=QSInew
            WRITE(6,"(4(a12,g11.4,1x))") &
            '{} TCold=',TC,'TCnew=',Tnew-T0C,'P=',.01*PP,'RHO=',RHO, &
            '{} THICK=',THICK,'RHold=',WV/WS,'RHnew=',WVnew/WSnew, &
            'RHgrd=',RHgrd, &
            '{} RHWold=',WV/QSW,'RHWnew=',WVnew/QSWnew,'RHIold=',WV/QSI, &
            'RHInew=',WVnew/QSInew, &
            '{} QSWold=',QSW,'QSWnew=',QSWnew,'QSIold=',QSI,'QSInew=',QSInew, &
            '{} WSold=',WS,'WSnew=',WSnew,'WVold=',WV,'WVnew=',WVnew, &
            '{} WCold=',WC,'WCnew=',WCnew,'QWold=',QW,'QWnew=',QWnew, &
            '{} QIold=',QI,'QInew=',QInew,'QRold=',QR,'QRnew=',QRnew, &
            '{} ARAINold=',ARAIN,'ARAINnew=',ARAINnew,'ASNOWold=',ASNOW, &
            'ASNOWnew=',ASNOWnew, &
            '{} TOT_RAIN=',TOT_RAIN,'TOT_RAINnew=',TOT_RAINnew, &
            'TOT_ICE=',TOT_ICE,'TOT_ICEnew=',TOT_ICEnew, &
            '{} BUDGET=',BUDGET,'QTold=',QT,'QTnew=',QTnew
        
            WRITE(6,"(4(a12,g11.4,1x))") &
            '{} DELT=',DELT,'DELV=',DELV,'DELW=',DELW,'DELI=',DELI, &
            '{} DELR=',DELR,'PCOND=',PCOND,'PIDEP=',PIDEP,'PIEVP=',PIEVP, &
            '{} PICND=',PICND,'PREVP=',PREVP,'PRAUT=',PRAUT,'PRACW=',PRACW, &
            '{} PIACW=',PIACW,'PIACWI=',PIACWI,'PIACWR=',PIACWR,'PIMLT=', &
            PIMLT, &
            '{} PIACR=',PIACR
        
            IF (ICE_logical) WRITE(6,"(4(a12,g11.4,1x))") &
            '{} RimeF1=',RimeF1,'GAMMAS=',GAMMAS,'VrimeF=',VrimeF, &
            'VSNOW=',VSNOW, &
            '{} INDEXS=',FLOAT(INDEXS),'FLARGE=',FLARGE,'FSMALL=',FSMALL, &
            'FLIMASS=',FLIMASS, &
            '{} XSIMASS=',XSIMASS,'XLIMASS=',XLIMASS,'QLICE=',QLICE, &
            'QTICE=',QTICE, &
            '{} NLICE=',NLICE,'NSmICE=',NSmICE,'PILOSS=',PILOSS, &
            'EMAIRI=',EMAIRI, &
            '{} RimeF=',RimeF
        
            IF (TOT_RAIN > 0. .OR. TOT_RAINnew > 0.) &
            WRITE(6,"(4(a12,g11.4,1x))") &
            '{} INDEXR1=',FLOAT(INDEXR1),'INDEXR=',FLOAT(INDEXR), &
            'GAMMAR=',GAMMAR,'N0r=',N0r, &
            '{} VRAIN1=',VRAIN1,'VRAIN2=',VRAIN2,'QTRAIN=',QTRAIN,'RQR=',RQR, &
            '{} PRLOSS=',PRLOSS,'VOLR1=',THICK+BLDTRH*VRAIN1, &
            'VOLR2=',THICK+BLDTRH*VRAIN2
        
            IF (PRAUT > 0.) WRITE(6,"(a12,g11.4,1x)") '{} QW0=',QW0
        
            IF (PRACW > 0.) WRITE(6,"(a12,g11.4,1x)") '{} FWR=',FWR
        
            IF (PIACR > 0.) WRITE(6,"(a12,g11.4,1x)") '{} FIR=',FIR
        
            DUM=PIMLT+PICND-PREVP-PIEVP
            IF (DUM > 0. .OR. DWVi /= 0.) &
            WRITE(6,"(4(a12,g11.4,1x))") &
            '{} TFACTOR=',TFACTOR,'DYNVIS=',DYNVIS, &
            'THERM_CON=',THERM_COND,'DIFFUS=',DIFFUS
        
            IF (PREVP < 0.) WRITE(6,"(4(a12,g11.4,1x))") &
            '{} RFACTOR=',RFACTOR,'ABW=',ABW,'VENTR=',VENTR,'CREVP=',CREVP, &
            '{} DWVr=',DWVr,'DENOMW=',DENOMW
        
            IF (PIDEP /= 0. .AND. DWVi /= 0.) &
            WRITE(6,"(4(a12,g11.4,1x))") &
            '{} DWVi=',DWVi,'DENOMI=',DENOMI,'PIDEP_max=',PIDEP_max, &
            'SFACTOR=',SFACTOR, &
            '{} ABI=',ABI,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
            'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
            '{} VENTIS=',VENTIS,'DIDEP=',DIDEP
        
            IF (PIDEP > 0. .AND. PCOND /= 0.) &
            WRITE(6,"(4(a12,g11.4,1x))") &
            '{} DENOMW=',DENOMW,'DENOMWI=',DENOMWI,'DENOMF=',DENOMF, &
            'DUM2=',PCOND-PIACW
        
            IF (FWS > 0.) WRITE(6,"(4(a12,g11.4,1x))") &
            '{} FWS=',FWS
        
            DUM=PIMLT+PICND-PIEVP
            IF (DUM > 0.) WRITE(6,"(4(a12,g11.4,1x))") &
            '{} SFACTOR=',SFACTOR,'VENTIL=',VENTIL,'VENTIL1=',VENTI1(INDEXS), &
            'VENTIL2=',SFACTOR*VENTI2(INDEXS), &
            '{} AIEVP=',AIEVP,'DIEVP=',DIEVP,'QSW0=',QSW0,'DWV0=',DWV0
        
        ENDIF
    
    !----------------------------------------------------------------------
    !-------------- Water budget statistics & maximum values --------------
    !----------------------------------------------------------------------
    
        IF (PRINT_diag) THEN
            ITdx=MAX( ITLO, MIN( INT(Tnew-T0C), ITHI ) )
            IF (QInew > CLIMIT) NSTATS(ITdx,1)=NSTATS(ITdx,1)+1
            IF (QInew > CLIMIT .AND.  QRnew+QWnew > CLIMIT) &
            NSTATS(ITdx,2)=NSTATS(ITdx,2)+1
            IF (QWnew > CLIMIT) NSTATS(ITdx,3)=NSTATS(ITdx,3)+1
            IF (QRnew > CLIMIT) NSTATS(ITdx,4)=NSTATS(ITdx,4)+1
        
            QMAX(ITdx,1)=MAX(QMAX(ITdx,1), QInew)
            QMAX(ITdx,2)=MAX(QMAX(ITdx,2), QWnew)
            QMAX(ITdx,3)=MAX(QMAX(ITdx,3), QRnew)
            QMAX(ITdx,4)=MAX(QMAX(ITdx,4), ASNOWnew)
            QMAX(ITdx,5)=MAX(QMAX(ITdx,5), ARAINnew)
            QTOT(ITdx,1)=QTOT(ITdx,1)+QInew*THICK
            QTOT(ITdx,2)=QTOT(ITdx,2)+QWnew*THICK
            QTOT(ITdx,3)=QTOT(ITdx,3)+QRnew*THICK
        
            QTOT(ITdx,4)=QTOT(ITdx,4)+PCOND*THICK
            QTOT(ITdx,5)=QTOT(ITdx,5)+PICND*THICK
            QTOT(ITdx,6)=QTOT(ITdx,6)+PIEVP*THICK
            QTOT(ITdx,7)=QTOT(ITdx,7)+PIDEP*THICK
            QTOT(ITdx,8)=QTOT(ITdx,8)+PREVP*THICK
            QTOT(ITdx,9)=QTOT(ITdx,9)+PRAUT*THICK
            QTOT(ITdx,10)=QTOT(ITdx,10)+PRACW*THICK
            QTOT(ITdx,11)=QTOT(ITdx,11)+PIMLT*THICK
            QTOT(ITdx,12)=QTOT(ITdx,12)+PIACW*THICK
            QTOT(ITdx,13)=QTOT(ITdx,13)+PIACWI*THICK
            QTOT(ITdx,14)=QTOT(ITdx,14)+PIACWR*THICK
            QTOT(ITdx,15)=QTOT(ITdx,15)+PIACR*THICK
        
            QTOT(ITdx,16)=QTOT(ITdx,16)+(WVnew-WV)*THICK
            QTOT(ITdx,17)=QTOT(ITdx,17)+(QWnew-QW)*THICK
            QTOT(ITdx,18)=QTOT(ITdx,18)+(QInew-QI)*THICK
            QTOT(ITdx,19)=QTOT(ITdx,19)+(QRnew-QR)*THICK
            QTOT(ITdx,20)=QTOT(ITdx,20)+(ARAINnew-ARAIN)
            QTOT(ITdx,21)=QTOT(ITdx,21)+(ASNOWnew-ASNOW)
            IF (QInew > 0.) &
            QTOT(ITdx,22)=QTOT(ITdx,22)+QInew*THICK/RimeF
        
        ENDIF
    
    !----------------------------------------------------------------------
    !------------------------- Update arrays ------------------------------
    !----------------------------------------------------------------------
    
        T_col(L)=Tnew                           ! Updated temperature
    
        QV_col(L)=max(EPSQ, WVnew/(1.+WVnew))   ! Updated specific humidity
        WC_col(L)=max(EPSQ, WCnew)              ! Updated total condensate mixing ratio
        QI_col(L)=max(EPSQ, QInew)              ! Updated ice mixing ratio
        QR_col(L)=max(EPSQ, QRnew)              ! Updated rain mixing ratio
        QW_col(L)=max(EPSQ, QWnew)              ! Updated cloud water mixing ratio
        RimeF_col(L)=RimeF                      ! Updated rime factor
        ASNOW=ASNOWnew                          ! Updated accumulated snow
        ARAIN=ARAINnew                          ! Updated accumulated rain
    
    !#######################################################################
    
        CONTINUE         ! ##### End "L" loop through model levels #####
    10 END DO

!#######################################################################

!-----------------------------------------------------------------------
!--------------------------- Return to GSMDRIVE -----------------------
!-----------------------------------------------------------------------

    RETURN
    END SUBROUTINE GSMCOLUMN

!#######################################################################
!--------- Produces accurate calculation of cloud condensation ---------
!#######################################################################

    REAL FUNCTION CONDENSE (PP, QW, RHgrd, TK, WV)

!---------------------------------------------------------------------------------
!------   The Asai (1965) algorithm takes into consideration the release of ------
!------   latent heat in increasing the temperature & in increasing the     ------
!------   saturation mixing ratio (following the Clausius-Clapeyron eqn.).  ------
!---------------------------------------------------------------------------------

    INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
    REAL (KIND=HIGH_PRES), PARAMETER :: CLIMIT=1.E-20, RHLIMIT=.001, RHLIMIT1=-RHLIMIT
    REAL, PARAMETER :: CP=1004.6, RD=287.04, RV=461.5, EPS=RD/RV, RCP=1./CP, RCPRV=RCP/RV
    REAL (KIND=HIGH_PRES) :: COND, SSAT, WCdum

!-----------------------------------------------------------------------

!--- LV (T) is from Bolton (JAS, 1980)

    XLV=3.148E6-2370.*TK
    XLV1=XLV*RCP
    XLV2=XLV*XLV*RCPRV
    Tdum=TK
    WVdum=WV
    WCdum=QW
    ESW=1000.*FPVS0(Tdum)                     ! Saturation vapor press w/r/t water
    WS=RHgrd*EPS*ESW/(PP-ESW)                 ! Saturation mixing ratio w/r/t water
    DWV=WVdum-WS                              ! Deficit grid-scale water vapor mixing ratio
    SSAT=DWV/WS                               ! Supersaturation ratio
    CONDENSE=0.
    DO WHILE ((SSAT < RHLIMIT1 .AND. WCdum > CLIMIT) &
         .OR. SSAT > RHLIMIT)
        COND=DWV/(1.+XLV2*WS/(Tdum*Tdum))       ! Asai (1965, J. Japan)
        COND=MAX(COND, -WCdum)                  ! Limit cloud water evaporation
        Tdum=Tdum+XLV1*COND                     ! Updated temperature
        WVdum=WVdum-COND                        ! Updated water vapor mixing ratio
        WCdum=WCdum+COND                        ! Updated cloud water mixing ratio
        CONDENSE=CONDENSE+COND                  ! Total cloud water condensation
        ESW=1000.*FPVS0(Tdum)                   ! Updated saturation vapor press w/r/t water
        WS=RHgrd*EPS*ESW/(PP-ESW)               ! Updated saturation mixing ratio w/r/t water
        DWV=WVdum-WS                            ! Deficit grid-scale water vapor mixing ratio
        SSAT=DWV/WS                             ! Grid-scale supersaturation ratio
    ENDDO
    RETURN
    END FUNCTION CONDENSE

!#######################################################################
!---------------- Calculate ice deposition at T<T_ICE ------------------
!#######################################################################

    REAL FUNCTION DEPOSIT (PP, RHgrd, Tdum, WVdum)

!--- Also uses the Asai (1965) algorithm, but uses a different target
!      vapor pressure for the adjustment

    INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
    REAL (KIND=HIGH_PRES), PARAMETER :: RHLIMIT=.001, RHLIMIT1=-RHLIMIT
    REAL, PARAMETER :: CP=1004.6, RD=287.04, RV=461.5, XLS=2.834E6, EPS=RD/RV, RCP=1./CP, RCPRV=RCP/RV, XLS1=XLS*RCP, &
    XLS2=XLS*XLS*RCPRV
    REAL (KIND=HIGH_PRES) :: DEP, SSAT

!-----------------------------------------------------------------------

    ESI=1000.*FPVS(Tdum)                      ! Saturation vapor press w/r/t ice
    WS=RHgrd*EPS*ESI/(PP-ESI)                 ! Saturation mixing ratio
    DWV=WVdum-WS                              ! Deficit grid-scale water vapor mixing ratio
    SSAT=DWV/WS                               ! Supersaturation ratio
    DEPOSIT=0.
    DO WHILE (SSAT > RHLIMIT .OR. SSAT < RHLIMIT1)
    
    !--- Note that XLVS2=LS*LV/(CP*RV)=LV*WS/(RV*T*T)*(LS/CP*DEP1),
    !     where WS is the saturation mixing ratio following Clausius-
    !     Clapeyron (see Asai,1965; Young,1993,p.405)
    
        DEP=DWV/(1.+XLS2*WS/(Tdum*Tdum))        ! Asai (1965, J. Japan)
        Tdum=Tdum+XLS1*DEP                      ! Updated temperature
        WVdum=WVdum-DEP                         ! Updated ice mixing ratio
        DEPOSIT=DEPOSIT+DEP                     ! Total ice deposition
        ESI=1000.*FPVS(Tdum)                    ! Updated saturation vapor press w/r/t ice
        WS=RHgrd*EPS*ESI/(PP-ESI)               ! Updated saturation mixing ratio w/r/t ice
        DWV=WVdum-WS                            ! Deficit grid-scale water vapor mixing ratio
        SSAT=DWV/WS                             ! Grid-scale supersaturation ratio
    ENDDO
    RETURN
    END FUNCTION DEPOSIT
