module_bfmicrophysics.f

 
 
      MODULE module_microphysics
!
      USE machine , ONLY : kind_phys
      USE funcphys
      USE physcons, cp => con_cp, rd => con_rd, rv => con_rv            &
     &,             t0c => con_t0c, hvap => con_hvap, hfus => con_hfus  &
     &,             eps => con_eps, epsm1 => con_epsm1                  &
     &,             eps1 => con_fvirt, pi => con_pi, grav => con_g 
      implicit none
!
!--- Common block of constants used in column microphysics
!
      real,private ::  abfr, cbfr, ciacw, ciacr, c_n0r0,                &
     &CN0r0, CN0r_DMRmin, CN0r_DMRmax, CRACW, CRAUT, ESW0,              &
     &QAUTx, RFmax, RQR_DR1, RQR_DR2, RQR_DR3, RQR_DRmin,               &
     &RQR_DRmax, RR_DRmin, RR_DR1, RR_DR2, RR_DR3, RR_DRmax
!
      integer, private :: mic_step
!
!--- Common block for lookup table used in calculating growth rates of
!    nucleated ice crystals growing in water saturated conditions
!--- 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, PRIVATE,PARAMETER :: my_t1=1, my_t2=35
      REAL,PRIVATE,DIMENSION(MY_T1:MY_T2) :: my_growth
!
!--- Parameters for ice lookup tables, which establish the range of mean ice
!    particle diameters; from a minimum mean diameter of 0.05 mm (DMImin) to a
!    maximum mean diameter of 1.00 mm (DMImax).  The tables store solutions
!    at 1 micron intervals (DelDMI) of mean ice particle diameter.
!
      REAL, PRIVATE,PARAMETER :: dmimin=.05e-3,      dmimax=1.e-3,      &
     &                           xmimin=1.e6*dmimin, xmimax=1.e6*dmimax,&
     &                           deldmi=1.e-6
      INTEGER, PRIVATE,PARAMETER :: mdimin=xmimin, mdimax=xmimax
!
!--- Various ice lookup tables
!
      REAL, PRIVATE,DIMENSION(MDImin:MDImax) ::                         &
     &      ACCRI,MASSI,SDENS,VSNOWI,VENTI1,VENTI2
!
!--- Mean rain drop diameters varying from 50 microns (0.05 mm) to 450 microns
!      (0.45 mm), assuming an exponential size distribution.
!
      REAL, PRIVATE,PARAMETER :: dmrmin=.05e-3,      dmrmax=.45e-3,     &
     &                           xmrmin=1.e6*dmrmin, xmrmax=1.e6*dmrmax,&
     &                           deldmr=1.e-6,       nlimin=100.
!    &,                          NLImin=100., NLImax=20.E3
      INTEGER, PRIVATE,PARAMETER :: mdrmin=xmrmin, mdrmax=xmrmax
!
!--- Factor of 1.5 for RECImin, RESNOWmin, & RERAINmin accounts for
!    integrating exponential distributions for effective radius
!    (i.e., the r**3/r**2 moments).
!
!     INTEGER, PRIVATE, PARAMETER :: INDEXSmin=300
!!    INTEGER, PRIVATE, PARAMETER :: INDEXSmin=200
      INTEGER, PRIVATE, PARAMETER :: indexsmin=100
      REAL, PRIVATE, PARAMETER :: rerainmin=1.5*xmrmin                  &
!    &, RECImin=1.5*XMImin, RESNOWmin=1.5*INDEXSmin, RECWmin=8.0
!    &, RECImin=1.5*XMImin, RESNOWmin=1.5*INDEXSmin, RECWmin=7.5
     &, recimin=1.5*xmimin, resnowmin=1.5*indexsmin, recwmin=10.
!    &, RECImin=1.5*XMImin, RESNOWmin=1.5*INDEXSmin, RECWmin=15.
!    &, RECImin=1.5*XMImin, RESNOWmin=1.5*INDEXSmin, RECWmin=5.
 
!
!--- Various rain lookup tables
!--- Rain lookup tables for mean rain drop diameters from DMRmin to DMRmax,
!      assuming exponential size distributions for the rain drops
!
      REAL, PRIVATE,DIMENSION(MDRmin:MDRmax)::                          &
     &      ACCRR,MASSR,RRATE,VRAIN,VENTR1,VENTR2
!
!--- Common block for riming tables
!--- 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, PRIVATE,PARAMETER :: nrime=40
      REAL, DIMENSION(2:9,0:Nrime),PRIVATE :: vel_rf
!
!--- The following variables are for microphysical statistics
!
      INTEGER, PARAMETER :: itlo=-60, ithi=40
      INTEGER  nstats(itlo:ithi,4)
      REAL     qmax(itlo:ithi,5),  qtot(itlo:ithi,22)
!
      REAL, PRIVATE,  PARAMETER ::                                      &
!    &  T_ICE=-10., T_ICE_init=-5.      !- Ver1
!!!  &, T_ICE=-20.                      !- Ver2
     &  t_ice=-40., t_ice_init=-15.     !- Ver2
!    &  T_ICE=-30., T_ICE_init=-5.      !- Ver2
!
!     Some other miscellaneous parameters
!
      REAL, PRIVATE, PARAMETER :: thom=t_ice, tnw=50., toler=1.0e-20    &
!     REAL, PRIVATE, PARAMETER :: Thom=T_ICE, TNW=50., TOLER=5.E-7
!     REAL, PRIVATE, PARAMETER :: Thom=-35., TNW=50., TOLER=5.E-7
 
!    &, emisCU=.75/1.66       ! Used for convective cloud l/w emissivity
 
! Assume fixed cloud ice effective radius
     &, recice=recimin                                                  &
     &, epsq=1.0e-20                                                    &
!    &, EPSQ=1.E-12                                                     &
     &, flg0p1=0.1, flg0p2=0.2, flg1p0=1.0
!
!
      CONTAINS
!
      SUBROUTINE gsmconst (DTPG,mype,first)
!
      implicit none
!-------------------------------------------------------------------------------
!---  SUBPROGRAM DOCUMENTATION BLOCK
!   PRGRMMR: Ferrier         ORG: W/NP22     DATE: February 2001
!-------------------------------------------------------------------------------
! ABSTRACT:
!   * Reads various microphysical lookup tables used in COLUMN_MICRO
!   * Lookup tables were created "offline" and are read in during execution
!   * Creates lookup tables for saturation vapor pressure w/r/t water & ice
!-------------------------------------------------------------------------------
!
! USAGE: CALL GSMCONST FROM SUBROUTINE GSMDRIVE AT MODEL START TIME
!
!   INPUT ARGUMENT LIST:
!       DTPH - physics time step (s)
!
!   OUTPUT ARGUMENT LIST:
!     NONE
!
!   OUTPUT FILES:
!     NONE
!
!
!   SUBROUTINES:
!     MY_GROWTH_RATES - lookup table for growth of nucleated ice
!
!   UNIQUE: NONE
!
!   LIBRARY: NONE
!
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!
      integer mype
      real    dtpg
      logical first
!
!--- Parameters & data statement for local calculations
!
      REAL, PARAMETER :: C1=1./3., dmr1=.1e-3, dmr2=.2e-3, dmr3=.32e-3, &
     & n0r0=8.e6, n0s0=4.e6, rhol=1000., rhos=100.,                     &
     & xmr1=1.e6*dmr1, xmr2=1.e6*dmr2, xmr3=1.e6*dmr3
      INTEGER, PARAMETER :: MDR1=xmr1, mdr2=xmr2, mdr3=xmr3
!
      real dtph, bbfr
      integer i
!
!--- Added on 5/16/01 for Moorthi
!
      logical, parameter :: read_lookup=.false., write_lookup=.false.
!
!------------------------------------------------------------------------
!  *************  Parameters used in ETA model -- Not used in Global Model *****
!
!--- DPHD, DLMD are delta latitude and longitude at the model (NOT geodetic) equator
!    => "DX" is the hypotenuse of the model zonal & meridional grid increments.
!
!     DX=111.*(DPHD**2+DLMD**2)**.5         ! Resolution at MODEL equator (km)
!     DX=MIN(100., MAX(5., DX) )
!
!--- Assume the following functional relationship for key constants that
!    depend on grid resolution from DXmin (5 km) to DXmax (100 km) resolution:
!
!     DXmin=5.
!     DXmax=100.
!     DX=MIN(DXmax, MAX(DXmin, DX) )
!
!--- EXtune determines the degree to which the coefficients change with resolution.
!    The larger EXtune is, the more sensitive the parameter.
!
!     EXtune=1.
 
!
!--- FXtune ==> F(DX) is the grid-resolution tuning parameter (from 0 to 1)
!
!     FXtune=((DXmax-DX)/(DXmax-DXmin))**EXtune
!     FXtune=MAX(0., MIN(1., FXtune))
!
!--- Calculate grid-averaged RH for the onset of condensation (RHgrd) based on
!    simple ***ASSUMED*** (user-specified) values at DXmax and at DXmin.
!
!     RH_DXmax=.90              !-- 90% RH at DXmax=100 km
!     RH_DXmin=.98              !-- 98% RH at DXmin=5 km
!
!--- Note that RHgrd is right now fixed throughout the domain!!
!
!     RHgrd=RH_DXmax+(RH_DXmin-RH_DXmax)*FXtune
!   ********************************************************************************
!
!
      if (first) then
!
!--- Read in various lookup tables
!
      if ( read_lookup ) then
        OPEN (unit=1,file="eta_micro_lookup.dat",form="UNFORMATTED")
        READ(1) ventr1
        READ(1) ventr2
        READ(1) accrr
        READ(1) massr
        READ(1) vrain
        READ(1) rrate
        READ(1) venti1
        READ(1) venti2
        READ(1) accri
        READ(1) massi
        READ(1) vsnowi
        READ(1) vel_rf
!       read(1) my_growth    ! Applicable only for DTPH=180 s for offline testing
        CLOSE (1)
      else
        CALL ice_lookup                   ! Lookup tables for ice
        CALL rain_lookup                  ! Lookup tables for rain
        if (write_lookup) then
          open(unit=1,file='micro_lookup.dat',form='unformatted')
          write(1) ventr1
          write(1) ventr2
          write(1) accrr
          write(1) massr
          write(1) vrain
          write(1) rrate
          write(1) venti1
          write(1) venti2
          write(1) accri
          write(1) massi
          write(1) vsnowi
          write(1) vel_rf
!         write(1) my_growth    ! Applicable only for DTPH=180 s ????
          CLOSE (1)
        endif
      endif
!!
!--- Constants associated with Biggs (1953) freezing of rain, as parameterized
!    following Lin et al. (JCAM, 1983) & Reisner et al. (1998, QJRMS).
!
      abfr=-0.66
      bbfr=100.
      cbfr=20.*pi*pi*bbfr*rhol*1.e-21
!
!--- QAUT0 is the threshold cloud content for autoconversion to rain
!      needed for droplets to reach a diameter of 20 microns (following
!      Manton and Cotton, 1977; Banta and Hanson, 1987, JCAM).  It is
!      **STRONGLY** affected by the assumed droplet number concentrations
!     XNCW!  For example, QAUT0=1.2567, 0.8378, or 0.4189 g/m**3 for
!     droplet number concentrations of 300, 200, and 100 cm**-3, respectively.
!
!--- Calculate grid-averaged XNCW based on simple ***ASSUMED*** (user-specified)
!    values at DXmax and at DXmin.
!
!     XNCW_DXmax=50.E6          !--  50 /cm**3 at DXmax=100 km
!     XNCW_DXmin=200.E6         !-- 200 /cm**3 at DXmin=5 km
!
!--- Note that XNCW is right now fixed throughout the domain!!
!
!     XNCW=XNCW_DXmax+(XNCW_DXmin-XNCW_DXmax)*FXtune
!
!     QAUT0=PI*RHOL*XNCW*(20.E-6)**3/6.
      qautx=pi*rhol*1.0e6*(20.e-6)**3/6.
!
!--- Based on rain lookup tables for mean diameters from 0.05 to 0.45 mm
!    * Four different functional relationships of mean drop diameter as
!      a function of rain rate (RR), derived based on simple fits to
!      mass-weighted fall speeds of rain as functions of mean diameter
!      from the lookup tables.
!
      rr_drmin=n0r0*rrate(mdrmin)     
      rr_dr1=n0r0*rrate(mdr1)         
      rr_dr2=n0r0*rrate(mdr2)         
      rr_dr3=n0r0*rrate(mdr3)         
      rr_drmax=n0r0*rrate(mdrmax)     
!
      rqr_drmin=n0r0*massr(mdrmin)    
      rqr_dr1=n0r0*massr(mdr1)        
      rqr_dr2=n0r0*massr(mdr2)        
      rqr_dr3=n0r0*massr(mdr3)        
      rqr_drmax=n0r0*massr(mdrmax)    
      c_n0r0=pi*rhol*n0r0
      cn0r0=1.e6/c_n0r0**.25
      cn0r_dmrmin=1./(pi*rhol*dmrmin**4)
      cn0r_dmrmax=1./(pi*rhol*dmrmax**4)
!
      endif                     !  If (first) then loop ends here
!
!     Find out what microphysics time step should be
!
      mic_step = max(1, int(dtpg/600.0+0.5))
!     mic_step = max(1, int(dtpg/300.0+0.5))
      dtph     = dtpg / mic_step
      if (mype == 0) print *,' DTPG=',dtpg,' mic_step=',mic_step        &
     &,                ' dtph=',dtph
!
!--- Calculates coefficients for growth rates of ice nucleated in water
!    saturated conditions, scaled by physics time step (lookup table)
!
      CALL my_growth_rates (dtph)
!
!--- CIACW is used in calculating riming rates
!      The assumed effective collection efficiency of cloud water rimed onto
!      ice is =0.5 below:
!
!Moor CIACW=DTPH*0.25*PI*0.5*(1.E5)**C1   ! commented on 20050422
!      ice is =0.1 below:
      ciacw=dtph*0.25*pi*0.1*(1.e5)**c1
!     CIACW = 0.0      ! Brad's suggestion 20040614
!
!--- CIACR is used in calculating freezing of rain colliding with large ice
!      The assumed collection efficiency is 1.0
!
      ciacr=pi*dtph
!
!--- CRACW is used in calculating collection of cloud water by rain (an
!      assumed collection efficiency of 1.0)
!
!Moor CRACW=DTPH*0.25*PI*1.0                 ! commented on 20050422
!
!      assumed collection efficiency of 0.1)
      cracw=dtph*0.25*pi*0.1
!     CRACW = 0.0      ! Brad's suggestion 20040614
!
      esw0=fpvsl(t0c)           
      rfmax=1.1**nrime          
!
!------------------------------------------------------------------------
!--------------- Constants passed through argument list -----------------
!------------------------------------------------------------------------
!
!--- Important parameters for self collection (autoconversion) of
!    cloud water to rain.
!
!--- CRAUT is proportional to the rate that cloud water is converted by
!      self collection to rain (autoconversion rate)
!
      craut=1.-exp(-1.e-3*dtph)
!
!     IF (MYPE == 0)
!    & WRITE(6,"(A, A,F6.2,A, A,F5.4, A,F7.3,A, A,F6.2,A, A,F6.3,A)")
!    &   'KEY MICROPHYSICAL PARAMETERS FOR '
!    &  ,'DX=',DX,' KM:'
!    &  ,'   FXtune=',FXtune
!    &  ,'   RHgrd=',100.*RHgrd,' %'
!    &  ,'   NCW=',1.E-6*XNCW,' /cm**3'
!    &  ,'   QAUT0=',1.E3*QAUT0,' g/kg'
!
!--- For calculating snow optical depths by considering bulk density of
!      snow based on emails from Q. Fu (6/27-28/01), where optical
!      depth (T) = 1.5*SWP/(Reff*DENS), SWP is snow water path, Reff
!      is effective radius, and DENS is the bulk density of snow.
!
!    SWP (kg/m**2)=(1.E-3 kg/g)*SWPrad, SWPrad in g/m**2 used in radiation
!    T = 1.5*1.E3*SWPrad/(Reff*DENS)
!
!    See derivation for MASSI(INDEXS), note equal to RHO*QSNOW/NSNOW
!
!    SDENS=1.5e3/DENS, DENS=MASSI(INDEXS)/[PI*(1.E-6*INDEXS)**3]
!
      DO i=mdimin,mdimax
!MoorthiSDENS(I)=PI*1.5E-15*FLOAT(I*I*I)/MASSI(I)
        sdens(i)=pi*1.0e-15*float(i*i*i)/massi(i)
      ENDDO
!
!-----------------------------------------------------------------------
!
      END subroutine gsmconst
 
!
      SUBROUTINE my_growth_rates (DTPH)
!
      implicit none
!
!--- Below are tabulated values for the predicted mass of ice crystals
!    after 600 s of growth in water saturated conditions, based on 
!    calculations from Miller and Young (JAS, 1979).  These values are
!    crudely estimated from tabulated curves at 600 s from Fig. 6.9 of
!    Young (1993).  Values at temperatures colder than -27C were 
!    assumed to be invariant with temperature.  
!
!--- Used to normalize Miller & Young (1979) calculations of ice growth
!    over large time steps using their tabulated values at 600 s.
!    Assumes 3D growth with time**1.5 following eq. (6.3) in Young (1993).
!
      real dtph, dt_ice
      REAL MY_600(MY_T1:MY_T2)
!
!-- 20090714: These values are in g and need to be converted to kg below
      DATA my_600 /                                                     &
     & 5.5e-8,  1.4e-7,  2.8e-7, 6.e-7,   3.3e-6,                       & !  -1 to  -5 deg C
     & 2.e-6,   9.e-7,   8.8e-7, 8.2e-7,  9.4e-7,                       & !  -6 to -10 deg C
     & 1.2e-6,  1.85e-6, 5.5e-6, 1.5e-5,  1.7e-5,                       & ! -11 to -15 deg C
     & 1.5e-5,  1.e-5,   3.4e-6, 1.85e-6, 1.35e-6,                      & ! -16 to -20 deg C
     & 1.05e-6, 1.e-6,   9.5e-7, 9.0e-7 , 9.5e-7,                       &  ! -21 to -25 deg C
     & 9.5e-7,  9.e-7,   9.e-7,  9.e-7,   9.e-7,                        &  ! -26 to -30 deg C
     & 9.e-7,   9.e-7,   9.e-7,  9.e-7,   9.e-7 /                         ! -31 to -35 deg C
!
!-----------------------------------------------------------------------
!
      dt_ice=(dtph/600.)**1.5
!     MY_GROWTH=DT_ICE*MY_600          ! original version
      my_growth=dt_ice*my_600*1.e-3    !-- 20090714: Convert from g to kg
!
!-----------------------------------------------------------------------
!
      END subroutine my_growth_rates
!
      subroutine ice_lookup
!
      implicit none
!-----------------------------------------------------------------------------------
!
!---- Key diameter values in mm
!
!-----------------------------------------------------------------------------------
!
!---- Key concepts:
!       - Actual physical diameter of particles (D)
!       - Ratio of actual particle diameters to mean diameter (x=D/MD)
!       - Mean diameter of exponentially distributed particles, which is the
!         same as 1./LAMDA of the distribution (MD)
!       - All quantitative relationships relating ice particle characteristics as
!         functions of their diameter (e.g., ventilation coefficients, normalized
!         accretion rates, ice content, and mass-weighted fall speeds) are a result
!         of using composite relationships for ice crystals smaller than 1.5 mm
!         diameter merged with analogous relationships for larger sized aggregates.
!         Relationships are derived as functions of mean ice particle sizes assuming
!         exponential size spectra and assuming the properties of ice crystals at
!         sizes smaller than 1.5 mm and aggregates at larger sizes.  
!
!-----------------------------------------------------------------------------------
!
!---- Actual ice particle diameters for which integrated distributions are derived
!       - DminI - minimum diameter for integration (.02 mm, 20 microns)
!       - DmaxI - maximum diameter for integration (2 cm)
!       - DdelI - interval for integration (1 micron)
!
      real, parameter :: DminI=.02e-3, dmaxi=20.e-3, ddeli=1.e-6,       &
     &  ximin=1.e6*dmini, ximax=1.e6*dmaxi
      integer, parameter :: IDImin=ximin, idimax=ximax
!
!---- Meaning of the following arrays:
!        - diam - ice particle diameter (m)
!        - mass - ice particle mass (kg)
!        - vel  - ice particle fall speeds (m/s)
!        - vent1 - 1st term in ice particle ventilation factor
!        - vent2 - 2nd term in ice particle ventilation factor
!
      real diam(IDImin:IDImax),mass(IDImin:IDImax),vel(IDImin:IDImax),  &
     & vent1(IDImin:IDImax),vent2(IDImin:IDImax)
!
!-----------------------------------------------------------------------------------
!
!---- Found from trial & error that the m(D) & V(D) mass & velocity relationships
!       between the ice crystals and aggregates overlapped & merged near a particle
!       diameter sizes of 1.5 mm.  Thus, ice crystal relationships are used for
!       sizes smaller than 1.5 mm and aggregate relationships for larger sizes.
!
      real, parameter :: d_crystal_max=1.5
!
!---- The quantity xmax represents the maximum value of "x" in which the
!       integrated values are calculated.  For xmax=20., this means that
!       integrated ventilation, accretion, mass, and precipitation rates are
!       calculated for ice particle sizes less than 20.*mdiam, the mean particle diameter.
!
      real, parameter :: xmax=20.
!
!-----------------------------------------------------------------------------------
!
!---- 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, used to calculate precip rates
!
!--- Mean ice-particle diameters varying from 50 microns to 1000 microns (1 mm), 
!      assuming an exponential size distribution.  
!
      real mdiam
!
!-----------------------------------------------------------------------------------
!------------- Constants & parameters for ventilation factors of ice ---------------
!-----------------------------------------------------------------------------------
!
!---- These parameters are used for calculating the ventilation factors for ice
!       crystals between 0.2 and 1.5 mm diameter (Hall and Pruppacher, JAS, 1976).  
!       From trial & error calculations, it was determined that the ventilation
!       factors of smaller ice crystals could be approximated by a simple linear
!       increase in the ventilation coefficient from 1.0 at 50 microns (.05 mm) to 
!       1.1 at 200 microns (0.2 mm), rather than using the more complex function of
!       1.0 + .14*(Sc**.33*Re**.5)**2 recommended by Hall & Pruppacher.
!
      real, parameter :: cvent1i=.86, cvent2i=.28
!
!---- These parameters are used for calculating the ventilation factors for larger
!       aggregates, where D>=1.5 mm (see Rutledge and Hobbs, JAS, 1983; 
!       Thorpe and Mason, 1966).
!
      real, parameter :: cvent1a=.65, cvent2a=.44
!
      real m_agg,m_bullet,m_column,m_ice,m_plate
!
!---- Various constants
!
      real, parameter :: c1=2./3., cexp=1./3.
!
      logical :: iprint
      logical, parameter :: print_diag=.false.
!
!-----------------------------------------------------------------------------------
!- Constants & parameters for calculating the increase in fall speed of rimed ice --
!-----------------------------------------------------------------------------------
!
!---- Constants & arrays for estimating increasing fall speeds of rimed ice.
!     Based largely on theory and results from Bohm (JAS, 1989, 2419-2427).
!
!-------------------- Standard atmosphere conditions at 1000 mb --------------------
!
      real, parameter :: t_std=288., dens_std=1000.e2/(287.04*288.)
!
!---- These "bulk densities" are the actual ice densities in the ice portion of the 
!     lattice.  They are based on text associated w/ (12) on p. 2425 of Bohm (JAS, 
!     1989).  Columns, plates, & bullets are assumed to have an average bulk density 
!     of 850 kg/m**3.  Aggregates were assumed to have a slightly larger bulk density 
!     of 600 kg/m**3 compared with dendrites (i.e., the least dense, most "lacy" & 
!     tenous ice crystal, which was assumed to be ~500 kg/m**3 in Bohm).  
!
      real, parameter :: dens_crystal=850., dens_agg=600.
!
!--- A value of Nrime=40 for a logarithmic ratio of 1.1 yields a maximum rime factor
!      of 1.1**40 = 45.26 that is resolved in these tables.  This allows the largest
!      ice particles with a mean diameter of MDImax=1000 microns to achieve bulk 
!      densities of 900 kg/m**3 for rimed ice.  
!
!     integer, parameter :: Nrime=40
      real m_rime,                                                      &
     &     rime_factor(0:Nrime), rime_vel(0:Nrime),                     &
     &     vel_rime(IDImin:IDImax,Nrime), ivel_rime(MDImin:MDImax,Nrime)
!
      integer i, j, jj, k, icount
      real c2,      cbulk, cbulk_ice, px, dynvis_std, crime1            &
     &,    crime2,  crime3, crime4, crime5, d, c_avg, c_agg             &
     &,    c_bullet, c_column, c_plate, cl_agg, cl_bullet               &
     &,    cl_column, cl_plate, v_agg, v_bullet, v_column               &
     &,    v_plate,   wd,       ecc_column                              &
     &,    cvent1,    cvent2, crime_best, rime_m1, rime_m2              &
     &,    x_rime,    re_rime, smom3, pratei, expf                      &
     &,    bulk_dens, xmass,  xmdiam, ecc_plate, dx
!
!-----------------------------------------------------------------------------------
!----------------------------- BEGIN EXECUTION -------------------------------------
!-----------------------------------------------------------------------------------
!
!
      c2=1./sqrt(3.)
!     pi=acos(-1.)
      cbulk=6./pi
      cbulk_ice=900.*pi/6.    ! Maximum bulk ice density allowed of 900 kg/m**3
      px=.4**cexp             ! Convert fall speeds from 400 mb (Starr & Cox) to 1000 mb
!
!--------------------- Dynamic viscosity (1000 mb, 288 K) --------------------------
!
      dynvis_std=1.496e-6*t_std**1.5/(t_std+120.)
      crime1=pi/24.
      crime2=8.*9.81*dens_std/(pi*dynvis_std**2)
      crime3=crime1*dens_crystal
      crime4=crime1*dens_agg
      crime5=dynvis_std/dens_std
      do i=0,nrime
        rime_factor(i)=1.1**i
      enddo
!
!#######################################################################
!      Characteristics as functions of actual ice particle diameter 
!#######################################################################
!
!----   M(D) & V(D) for 3 categories of ice crystals described by Starr 
!----   & Cox (1985). 
!
!----   Capacitance & characteristic lengths for Reynolds Number calculations
!----   are based on Young (1993; p. 144 & p. 150).  c-axis & a-axis 
!----   relationships are from Heymsfield (JAS, 1972; Table 1, p. 1351).
!
      icount=60
!
      if (print_diag)                                                   & 
     &  write(7,"(2a)") '---- Increase in fall speeds of rimed ice',    &
     &    ' particles as function of ice particle diameter ----'
      do i=idimin,idimax
        if (icount == 60 .and. print_diag) then
          write(6,"(/2a/3a)") 'Particle masses (mg), fall speeds ',     &
     &      '(m/s), and ventilation factors',                           &
     &      '  D(mm)  CR_mass   Mass_bull   Mass_col  Mass_plat ',      &
     &      '  Mass_agg   CR_vel  V_bul CR_col CR_pla Aggreg',          &
     &      '    Vent1      Vent2 '                               
          write(7,"(3a)") '        <----------------------------------',&
     &      '---------------  Rime Factor  --------------------------', &
     &      '--------------------------->'
          write(7,"(a,23f5.2)") '  D(mm)',(rime_factor(k), k=1,5),      &
     &       (rime_factor(k), k=6,40,2)
          icount=0
        endif
        d=(float(i)+.5)*1.e-3    ! in mm
        c_avg=0.
        c_agg=0.
        c_bullet=0.
        c_column=0.
        c_plate=0.
        cl_agg=0.
        cl_bullet=0.
        cl_column=0.
        cl_plate=0.
        m_agg=0.
        m_bullet=0.
        m_column=0.
        m_plate=0.
        v_agg=0.
        v_bullet=0.
        v_column=0.
        v_plate=0.
        if (d < d_crystal_max) then
!
!---- This block of code calculates bulk characteristics based on average
!     characteristics of bullets, plates, & column ice crystals <1.5 mm size
!
!---- Mass-diameter relationships from Heymsfield (1972) & used
!       in Starr & Cox (1985), units in mg
!---- "d" is maximum dimension size of crystal in mm, 
!
! Mass of pure ice for spherical particles, used as an upper limit for the
!   mass of small columns (<~ 80 microns) & plates (<~ 35 microns)
!
          m_ice=.48*d**3   ! Mass of pure ice for spherical particle
!
          m_bullet=min(.044*d**3, m_ice)
          m_column=min(.017*d**1.7, m_ice)
          m_plate=min(.026*d**2.5, m_ice)
!
          mass(i)=m_bullet+m_column+m_plate
!
!---- These relationships are from Starr & Cox (1985), applicable at 400 mb
!---- "d" is maximum dimension size of crystal in mm, dx in microns
!
          dx=1000.*d            ! Convert from mm to microns
          if (dx <= 200.) then
            v_column=8.114e-5*dx**1.585
            v_bullet=5.666e-5*dx**1.663
            v_plate=1.e-3*dx
          else if (dx <= 400.) then
            v_column=4.995e-3*dx**.807
            v_bullet=3.197e-3*dx**.902
            v_plate=1.48e-3*dx**.926
          else if (dx <= 600.) then
            v_column=2.223e-2*dx**.558
            v_bullet=2.977e-2*dx**.529
            v_plate=9.5e-4*dx
          else if (dx <= 800.) then
            v_column=4.352e-2*dx**.453
            v_bullet=2.144e-2*dx**.581
            v_plate=3.161e-3*dx**.812
          else 
            v_column=3.833e-2*dx**.472
            v_bullet=3.948e-2*dx**.489
            v_plate=7.109e-3*dx**.691
          endif
!
!---- Reduce fall speeds from 400 mb to 1000 mb
!
          v_column=px*v_column
          v_bullet=px*v_bullet
          v_plate=px*v_plate
!
!---- DIFFERENT VERSION!  CALCULATES MASS-WEIGHTED CRYSTAL FALL SPEEDS
!
          vel(i)=(m_bullet*v_bullet+m_column*v_column+m_plate*v_plate)/ &
     &           mass(i)
          mass(i)=mass(i)/3.
!
!---- Shape factor and characteristic length of various ice habits,
!     capacitance is equal to 4*PI*(Shape factor)
!       See Young (1993, pp. 143-152 for guidance)
!
!---- Bullets:
!
!---- Shape factor for bullets (Heymsfield, 1975)
          c_bullet=.5*d
!---- Length-width functions for bullets from Heymsfield (JAS, 1972)
          if (d > 0.3) then
            wd=.25*d**.7856     ! Width (mm); a-axis
          else
            wd=.185*d**.552
          endif
!---- Characteristic length for bullets (see first multiplicative term on right
!       side of eq. 7 multiplied by crystal width on p. 821 of Heymsfield, 1975)
          cl_bullet=.5*pi*wd*(.25*wd+d)/(d+wd)
!
!---- Plates:
!
!---- Length-width function for plates from Heymsfield (JAS, 1972)
          wd=.0449*d**.449      ! Width or thickness (mm); c-axis
!---- Eccentricity & shape factor for thick plates following Young (1993, p. 144)
          ecc_plate=sqrt(1.-wd*wd/(d*d))         ! Eccentricity
          c_plate=d*ecc_plate/asin(ecc_plate)    ! Shape factor
!---- Characteristic length for plates following Young (1993, p. 150, eq. 6.6)
          cl_plate=d+2.*wd      ! Characteristic lengths for plates
!
!---- Columns:
!
!---- Length-width function for columns from Heymsfield (JAS, 1972)
          if (d > 0.2) then
            wd=.1973*d**.414    ! Width (mm); a-axis
          else
            wd=.5*d             ! Width (mm); a-axis
          endif
!---- Eccentricity & shape factor for columns following Young (1993, p. 144)
          ecc_column=sqrt(1.-wd*wd/(d*d))                     ! Eccentricity
          c_column=ecc_column*d/log((1.+ecc_column)*d/wd)     ! Shape factor
!---- Characteristic length for columns following Young (1993, p. 150, eq. 6.7)
          cl_column=(wd+2.*d)/(c1+c2*d/wd)       ! Characteristic lengths for columns
!
!---- Convert shape factor & characteristic lengths from mm to m for 
!       ventilation calculations
!
          c_bullet=.001*c_bullet
          c_plate=.001*c_plate
          c_column=.001*c_column
          cl_bullet=.001*cl_bullet
          cl_plate=.001*cl_plate
          cl_column=.001*cl_column
!
!---- Make a smooth transition between a ventilation coefficient of 1.0 at 50 microns
!       to 1.1 at 200 microns
!
          if (d > 0.2) then
            cvent1=cvent1i
            cvent2=cvent2i/3.
          else
            cvent1=1.0+.1*max(0., d-.05)/.15
            cvent2=0.
          endif
!
!---- Ventilation factors for ice crystals:
!
          vent1(i)=cvent1*(c_bullet+c_plate+c_column)/3.
          vent2(i)=cvent2*(c_bullet*sqrt(v_bullet*cl_bullet)            &
     &                    +c_plate*sqrt(v_plate*cl_plate)               &
     &                    +c_column*sqrt(v_column*cl_column) )
          crime_best=crime3     ! For calculating Best No. of rimed ice crystals
        else
!
!---- This block of code calculates bulk characteristics based on average
!     characteristics of unrimed aggregates >= 1.5 mm using Locatelli & 
!     Hobbs (JGR, 1974, 2185-2197) data.
!
!----- This category is a composite of aggregates of unrimed radiating 
!-----   assemblages of dendrites or dendrites; aggregates of unrimed
!-----   radiating assemblages of plates, side planes, bullets, & columns;
!-----   aggregates of unrimed side planes (mass in mg, velocity in m/s)
!
          m_agg=(.073*d**1.4+.037*d**1.9+.04*d**1.4)/3.
          v_agg=(.8*d**.16+.69*d**.41+.82*d**.12)/3.
          mass(i)=m_agg
          vel(i)=v_agg
!
!---- Assume spherical aggregates
!
!---- Shape factor is the same as for bullets, = D/2
          c_agg=.001*.5*d         ! Units of m
!---- Characteristic length is surface area divided by perimeter
!       (.25*PI*D**2)/(PI*D**2) = D/4
          cl_agg=.5*c_agg         ! Units of m
!
!---- Ventilation factors for aggregates:
!
          vent1(i)=cvent1a*c_agg
          vent2(i)=cvent2a*c_agg*sqrt(v_agg*cl_agg)
          crime_best=crime4     ! For calculating Best No. of rimed aggregates
        endif
!
!---- Convert from shape factor to capacitance for ventilation factors
!
        vent1(i)=4.*pi*vent1(i)
        vent2(i)=4.*pi*vent2(i)
        diam(i)=1.e-3*d             ! Convert from mm to m
        mass(i)=1.e-6*mass(i)       ! Convert from mg to kg
!
!---- Calculate increase in fall speeds of individual rimed ice particles
!
        do k=0,nrime
!---- Mass of rimed ice particle associated with rime_factor(k)
          rime_m1=rime_factor(k)*mass(i)
          rime_m2=cbulk_ice*diam(i)**3
          m_rime=min(rime_m1, rime_m2)
!---- Best Number (X) of rimed ice particle combining eqs. (8) & (12) in Bohm
          x_rime=crime2*m_rime*(crime_best/m_rime)**.25
!---- Reynolds Number for rimed ice particle using eq. (11) in Bohm
          re_rime=8.5*(sqrt(1.+.1519*sqrt(x_rime))-1.)**2
          rime_vel(k)=crime5*re_rime/diam(i)
        enddo
        do k=1,nrime
          vel_rime(i,k)=rime_vel(k)/rime_vel(0)
        enddo
        if (print_diag) then
   !
   !---- Determine if statistics should be printed out.
   !
          iprint=.false.
          if (d <= 1.) then
            if (mod(i,10) == 0) iprint=.true.
          else
            if (mod(i,100) == 0) iprint=.true.
          endif
          if (iprint) then
            write(6,"(f7.4,5e11.4,1x,5f7.4,1x,2e11.4)")                 & 
     &        d,1.e6*mass(i),m_bullet,m_column,m_plate,m_agg,           &
     &        vel(i),v_bullet,v_column,v_plate,v_agg,                   &
     &        vent1(i),vent2(i)
            write(7,"(f7.4,23f5.2)") d,(vel_rime(i,k), k=1,5),          &
     &        (vel_rime(i,k), k=6,40,2)
            icount=icount+1
          endif
        endif
      enddo
!
!#######################################################################
!      Characteristics as functions of mean particle diameter
!#######################################################################
!
      venti1=0.
      venti2=0.
      accri=0.
      massi=0.
      vsnowi=0.
      vel_rf=0.
      ivel_rime=0.
      icount=0
      if (print_diag) then
        icount=60
        write(6,"(/2a)") '------------- Statistics as functions of ',   &
     &    'mean particle diameter -------------'
        write(7,"(/2a)") '------ Increase in fall speeds of rimed ice', &
     &    ' particles as functions of mean particle diameter -----'
      endif
      do j=mdimin,mdimax
        if (icount == 60 .AND. print_diag) then
          write(6,"(/2a)") 'D(mm)    Vent1      Vent2    ',             &
     &       'Accrete       Mass     Vel  Dens  '
          write(7,"(3a)") '      <----------------------------------',  &
     &      '---------------  Rime Factor  --------------------------', &
     &      '--------------------------->'
          write(7,"(a,23f5.2)") 'D(mm)',(rime_factor(k), k=1,5),        &
     &       (rime_factor(k), k=6,40,2)
          icount=0
        endif
        mdiam=deldmi*float(j)       ! in m
        smom3=0.
        pratei=0.
        rime_vel=0.                 ! Note that this array is being reused!
        do i=idimin,idimax
          dx=diam(i)/mdiam
          if (dx <= xmax) then      ! To prevent arithmetic underflows
            expf=exp(-dx)*ddeli
            venti1(j)=venti1(j)+vent1(i)*expf
            venti2(j)=venti2(j)+vent2(i)*expf
            accri(j)=accri(j)+diam(i)*diam(i)*vel(i)*expf
            xmass=mass(i)*expf
            do k=1,nrime
              rime_vel(k)=rime_vel(k)+xmass*vel_rime(i,k)
            enddo
            massi(j)=massi(j)+xmass
            pratei=pratei+xmass*vel(i)
            smom3=smom3+diam(i)**3*expf
          else
            exit
          endif
        enddo
   !
   !--- Increased fall velocities functions of mean diameter (j),
   !      normalized by ice content, and rime factor (k) 
   !
        do k=1,nrime
          ivel_rime(j,k)=rime_vel(k)/massi(j)
        enddo
   !
   !--- Increased fall velocities functions of ice content at 0.1 mm
   !      intervals (j_100) and rime factor (k); accumulations here
   !
        jj=j/100
        if (jj >= 2 .AND. jj <= 9) then
          do k=1,nrime
            vel_rf(jj,k)=vel_rf(jj,k)+ivel_rime(j,k)
          enddo
        endif
        bulk_dens=cbulk*massi(j)/smom3
        venti1(j)=venti1(j)/mdiam
        venti2(j)=venti2(j)/mdiam
        accri(j)=accri(j)/mdiam
        vsnowi(j)=pratei/massi(j)
        massi(j)=massi(j)/mdiam
        if (mod(j,10) == 0 .AND. print_diag) then
          xmdiam=1.e3*mdiam
          write(6,"(f6.3,4e11.4,f6.3,f8.3)") xmdiam,venti1(j),venti2(j),&
     &      accri(j),massi(j),vsnowi(j),bulk_dens
          write(7,"(f6.3,23f5.2)") xmdiam,(ivel_rime(j,k), k=1,5),      &
     &       (ivel_rime(j,k), k=6,40,2)
          icount=icount+1
        endif
      enddo
!
!--- Average increase in fall velocities rimed ice as functions of mean
!      particle diameter (j, only need 0.1 mm intervals) and rime factor (k)
!
      if (print_diag) then
        write(7,"(/2a)") ' ------- Increase in fall speeds of rimed ',  &
     &    'ice particles at reduced, 0.1-mm intervals  --------'
        write(7,"(3a)") '        <----------------------------------',  &
     &    '---------------  Rime Factor  --------------------------',   &
     &    '--------------------------->'
        write(7,"(a,23f5.2)") 'D(mm)',(rime_factor(k), k=1,5),          &
     &    (rime_factor(k), k=6,40,2)
      endif
      do j=2,9
        vel_rf(j,0)=1.
        do k=1,nrime
          vel_rf(j,k)=.01*vel_rf(j,k)
        enddo
        if (print_diag) write(7,"(f4.1,2x,23f5.2)") 0.1*j,              &
     &    (vel_rf(j,k), k=1,5),(vel_rf(j,k), k=6,40,2)
      enddo
!
!-----------------------------------------------------------------------------------
!
      end subroutine ice_lookup
!
      subroutine rain_lookup
      implicit none
!
!--- Parameters & arrays for fall speeds of rain as a function of rain drop
!      diameter.  These quantities are integrated over exponential size
!      distributions of rain drops at 1 micron intervals (DdelR) from minimum 
!      drop sizes of .05 mm (50 microns, DminR) to maximum drop sizes of 10 mm 
!      (DmaxR). 
!
      real, parameter :: DminR=.05e-3, dmaxr=10.e-3, ddelr=1.e-6,       &
     & xrmin=1.e6*dminr, xrmax=1.e6*dmaxr
      integer, parameter :: IDRmin=xrmin, idrmax=xrmax
      real diam(IDRmin:IDRmax), vel(IDRmin:IDRmax)
!
!--- Parameters rain lookup tables, which establish the range of mean drop
!      diameters; from a minimum mean diameter of 0.05 mm (DMRmin) to a 
!      maximum mean diameter of 0.45 mm (DMRmax).  The tables store solutions
!      at 1 micron intervals (DelDMR) of mean drop diameter.  
!
      real mdiam, mass
!
      logical, parameter :: print_diag=.false.
!
      real d, cmass, pi2, expf
      integer i, j, i1, i2
!
!-----------------------------------------------------------------------
!------- Fall speeds of rain as function of rain drop diameter ---------
!-----------------------------------------------------------------------
!
      do i=idrmin,idrmax
        diam(i)=float(i)*ddelr
        d=100.*diam(i)         ! Diameter in cm
        if (d <= .42) then
   !
   !--- Rutledge & Hobbs (1983); vel (m/s), d (cm)
   !
          vel(i)=max(0., -.267+51.5*d-102.25*d*d+75.7*d**3)
        else if (d > 0.42 .and. d <= .58) then
   !
   !--- Linear interpolation of Gunn & Kinzer (1949) data
   !
          vel(i)=8.92+.25/(.58-.42)*(d-.42)
        else
          vel(i)=9.17
        endif
      enddo
      do i=1,100
        i1=(i-1)*100+idrmin
        i2=i1+90
   !
   !--- Print out rain fall speeds only for D<=5.8 mm (.58 cm)
   !
        if (diam(i1) > .58e-2) exit
        if (print_diag) then
          write(6,"(1x)")
          write(6,"('D(mm)->  ',10f7.3)") (1000.*diam(j), j=i1,i2,10)
          write(6,"('V(m/s)-> ',10f7.3)") (vel(j), j=i1,i2,10)
        endif
      enddo
!
!-----------------------------------------------------------------------
!------------------- Lookup tables for rain processes ------------------
!-----------------------------------------------------------------------
!
!     pi=acos(-1.)
      pi2=2.*pi
      cmass=1000.*pi/6.
      if (print_diag) then
        write(6,"(/'Diam - Mean diameter (mm)'                          &
     &          /'VENTR1 - 1st ventilation coefficient (m**2)'          &
     &          /'VENTR2 - 2nd ventilation coefficient (m**3/s**.5)'    &
     &          /'ACCRR - accretion moment (m**4/s)'                    &
     &          /'RHO*QR - mass content (kg/m**3) for N0r=8e6'          &
     &          /'RRATE - rain rate moment (m**5/s)'                    &
     &          /'VR - mass-weighted rain fall speed (m/s)'             &
     &    /' Diam      VENTR1      VENTR2       ACCRR      ',           &
     &    'RHO*QR       RRATE    VR')")
      endif
      do j=mdrmin,mdrmax
        mdiam=float(j)*deldmr
        ventr2(j)=0.
        accrr(j)=0.
        massr(j)=0.
        rrate(j)=0.
        do i=idrmin,idrmax
          expf=exp(-diam(i)/mdiam)*ddelr
          ventr2(j)=ventr2(j)+diam(i)**1.5*vel(i)**.5*expf
          accrr(j)=accrr(j)+diam(i)*diam(i)*vel(i)*expf
          massr(j)=massr(j)+diam(i)**3*expf
          rrate(j)=rrate(j)+diam(i)**3*vel(i)*expf
        enddo
   !
   !--- Derived based on ventilation, F(D)=0.78+.31*Schmidt**(1/3)*Reynold**.5,
   !      where Reynold=(V*D*rho/dyn_vis), V is velocity, D is particle diameter,
   !      rho is air density, & dyn_vis is dynamic viscosity.  Only terms 
   !      containing velocity & diameter are retained in these tables.  
   !
        ventr1(j)=.78*pi2*mdiam**2
        ventr2(j)=.31*pi2*ventr2(j)
   !
        massr(j)=cmass*massr(j)
        rrate(j)=cmass*rrate(j)
        vrain(j)=rrate(j)/massr(j)
        if (print_diag) write(6,"(f6.3,5g12.5,f6.3)") 1000.*mdiam,      &
     &    ventr1(j),ventr2(j),accrr(j),8.e6*massr(j),rrate(j),vrain(j)
      enddo
!
!-----------------------------------------------------------------------
!
      end subroutine rain_lookup
!
!###############################################################################
! ***** 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 ( ARAING, ASNOWG, DTPG, I_index, J_index,    &
     & LSFC, P_col, QI_col, QR_col, QV_col, QW_col, RimeF_col, T_col,   &
     & THICK_col, WC_col, LM, RHC_col, XNCW, FLGmin, PRINT_diag, psfc)
!
      implicit none
!
!###############################################################################
!###############################################################################
!
!-------------------------------------------------------------------------------
!----- NOTE:  In this version of the Code threading should be done outside!  
!-------------------------------------------------------------------------------
!$$$  SUBPROGRAM DOCUMENTATION BLOCK
!                .      .    .     
! SUBPROGRAM:  Grid-scale microphysical processes - condensation & precipitation
!   PRGRMMR: Ferrier         ORG: W/NP22     DATE: 08-2001
!   Updated: Moorthi for GFS application
!-------------------------------------------------------------------------------
! 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
!  
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!   MACHINE : IBM SP
!
!------------------------------------------------------------------------- 
!--------------- Arrays & constants in argument list --------------------- 
!------------------------------------------------------------------------- 
!
      integer lm
      REAL ARAING, ASNOWG, 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), RHC_col(LM), XNCW(LM), ARAIN, ASNOW, dtpg, psfc
      real flgmin
!
      INTEGER I_index, J_index, LSFC
!
!
!------------------------------------------------------------------------- 
!
!--- 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
!
!------------------------------------------------------------------------- 
!------- Key parameters, local variables, & important comments ---------
!-----------------------------------------------------------------------
!
!--- KEY Parameters:
!
!---- Comments on 14 March 2002
!    * Set EPSQ to the universal value of 1.e-12 throughout the code
!      condensate.  The value of EPSQ will need to be changed in the other 
!      subroutines in order to make it consistent throughout the Eta code.  
!    * Set CLIMIT=10.*EPSQ as the lower limit for the total mass of 
!      condensate in the current layer and the input flux of condensate
!      from above (TOT_ICE, TOT_ICEnew, TOT_RAIN, and TOT_RAINnew).
!
!-- 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 ::   RHOL=1000.,  xls=hvap+hfus                   &
 
!    &, T_ICE=-10.          !- Ver1
!    &, T_ICE_init=-5.      !- Ver1
!!!  &, T_ICE=-20.          !- Ver2
!    &, T_ICE=-40.          !- Ver2
!    &, T_ICE_init=-15.,    !- Ver2
!
!    & CLIMIT=10.*EPSQ, EPS1=RV/RD-1., RCP=1./CP,
 
     &,climit=10.*epsq, 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  PRINT_diag
!
!--- Local variables
!
      REAL    EMAIRI, N0r,         NLICE,       NSmICE, NLImax, pfac
      LOGICAL CLEAR,  ICE_logical, DBG_logical, RAIN_logical
 
      integer lbef, ipass, ixrf, ixs, itdx, idr                         &
     &,       index_my, indexr, indexr1, indexs                         &
     &,       i, j, k, l, ntimes, item
!    &,       i, j, k, my_600, i1, i2, l, ntimes
 
      real flimass,  xlimass, vsnow,   qi_min, dum,    piloss           &
     &,    tot_ice,  xsimass, vel_inc, vrimef, rimef1, dum1             &
     &,    dum2,     fws,     denomi,  dwv                              &
     &,    xrf,      qw0,     dli,     xli,    fsmall                   &
     &,    prevp,    tk2,     dtph                                      &
     &,    pievp,    picnd,   piacr,   pracw                            &
     &,    praut,    pimlt,   qtice,   qlice                            &
     &,    gammar,   flarge,  wvqw,    dynvis                           &
     &,    tfactor,  denom,   gammas,  diffus, therm_cond               &
     &,    wvnew,    delv,    tnew,    tot_icenew, rimef                &
     &,    deli,     fwr,     crevp,   ventr,      delt                 &
     &,    delw,     fir,     delr,    qsinew,     qswnew               &
     &,    budget,   wsnew,   vrain2,  tot_rainnew                      &
     &,    qtnew,    qt,      wcnew,   abw                              &
     &,    aievp,    tcc,     denomf,  abi                              &
     &,    sfactor,  pidep_max,        didep,      ventis, ventil       &
     &,    dievp,    rqr,     rfactor, dwvr,       rr,     tot_rain     &
     &,    dwv0,     qsw0,    prloss,  qtrain,     vrain1               &
     &,    qsw,      ws,      esi,     esw, wv, wc, rhgrd, rho          &
     &,    rrho,     dtrho,   wsgrd,   qsi, qswgrd, qsigrd              &
     &,    tk,       tc,      pp,      bldtrh                           &
     &,    xlv,      xlv1,    xlf,     xlf1,  xlv2, denomw, denomwi     &
     &,    qwnew,    pcond,   pidep,   qrnew, qi,   qr,     qw          &
     &,    piacw,    piacwi,  piacwr,  qv,    dwvi                      &
     &,    arainnew, thick,   asnownew                                  &
     &,    qinew,    qi_min_0c, QSW_l, QSI_l, QSW0_l, SCHMIT_FAC
    
!
!
!#######################################################################
!########################## Begin Execution ############################
!#######################################################################
!
      dtph   = dtpg / mic_step
      araing = 0.    ! Total Accumulated rainfall at surface (kg/m**2)
      asnowg = 0.    ! Total Accumulated snowfall at surface (kg/m**2)
!
      do ntimes =1,mic_step
!
        qi_min_0c = 10.e3*massi(mdimin)   !- Ver5
        arain = 0.   ! Accumulated rainfall at surface for this step (kg/m**2)
        asnow = 0.   ! Accumulated snowfall at surface for this step (kg/m**2)
 
        indexr   = mdrmin
!
!-----------------------------------------------------------------------
!
        DO l=1,lsfc      !      Loop from top (L=1) to surface (L=LSFC)
 
!---      Skip this level and go to the next lower level if no condensate 
!         and very low specific humidities
!
          IF (qv_col(l) > epsq .OR. wc_col(l) > epsq) THEN
!
!-----------------------------------------------------------------------
!------------ 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)
            wv = qv               ! Water vapor specific humidity (kg/kg)
            wc = wc_col(l)        ! Grid-scale mixing ratio of total condensate
                                ! (water or ice; kg/kg)
!           WC = WC/(1.-WC)
            rhgrd = rhc_col(l)
 
!
!   Pressure dependen scaling factor for flgmin (tunable)
!
!!!         pfac = max(0.5, (min(1.0, pp*0.00002))**2)   ! commented on 02182011
!go         pfac = max(0.5, (sqrt(min(1.0, pp*0.00004))))
            pfac = 1.0
!
            clear = .true.
!    
!--- Check grid-scale saturation when no condensate is present
!    
            esw = min(pp, fpvsl(tk))     ! Saturation vapor pressure w/r/t water
!           QSW = EPS*ESW/(PP-ESW)       ! Saturation mixing ratio w/r/t water
            qsw = eps*esw/(pp+epsm1*esw) ! Saturation specific humidity  w/r/t water
            ws  = qsw                    ! General saturation mixing ratio (water/ice)
            qsi = qsw
            IF (tc < 0.) THEN
              esi = min(pp,fpvsi(tk))      ! Saturation vapor pressure w/r/t ice
!             QSI = EPS*ESI/(PP-ESI)       ! Saturation mixing ratio w/r/t water
              qsi = eps*esi/(pp+epsm1*esi) ! Saturation specific humidity w/r/t water
              ws  = qsi                    ! General saturation mixing ratio (water/ice)
              if (pp <= esi) ws = wv / rhgrd
            ENDIF
!
            dum  = min(pp, esw0)
            qsw0 = eps*dum/(pp+epsm1*dum)  ! Saturation specific Humidity at 0C
!
!--- Effective grid-scale Saturation mixing ratios
!
            qswgrd = rhgrd*qsw
            qsigrd = rhgrd*qsi
            wsgrd  = rhgrd*ws
            qsw_l  = qswgrd
            qsi_l  = qsigrd
            qsw0_l = qsw0*rhgrd
!
!--- Check if air is subsaturated and w/o condensate
!
            IF (wv > wsgrd .OR. wc > epsq) clear = .false.  ! Cloudy case
            IF (arain > climit) THEN ! If any rain is falling into layer from above
              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 (.not. clear) THEN
!
!-----------------------------------------------------------------------
!--------- 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
!
              arainnew = 0.           ! Updated accumulated rainfall at surface
              asnownew = 0.           ! Updated accumulated snowfall at surface
              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 >  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_l*tk2/rv  ! Lv**2*Qsw_l/(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.)
              tfactor    = tk*sqrt(tk)/(tk+120.)
              dynvis     = 1.496e-6*tfactor
              therm_cond = 2.116e-3*tfactor
              diffus     = 8.794e-5*tk**1.81/pp
              schmit_fac = (rho/(diffus*diffus*dynvis))**c2
!
!--- 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, JAOT, p.470)
!
              gammar = (rho0/rho)**0.4
!
!----------------------------------------------------------------------
!-------------  IMPORTANT MICROPHYSICS DECISION TREE  -----------------
!----------------------------------------------------------------------
!
!--- Determine if conditions supporting ice are present
!
              IF (tc < 0. .OR. qi > epsq .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 > epsq) 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
!
!--- 6/19/03 - Deleted some code here ....
!
!  *********************************************************
 
!               IF (TC >= 0. .OR. WVQW < QSIgrd) THEN
!  !
!  !--- Eliminate small ice particle contributions for melting & sublimation
!  !
!                 FLARGE = FLARGE1
!               ELSE
!  !
!  !--- Enhanced number of small ice particles during depositional growth
!  !    (effective only when 0C > T >= T_ice [-10C] )
!  !
!                 FLARGE = FLARGE2
!  !
!  !--- 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 <= epsq .AND. asnow <= climit) THEN
                  indexs  = mdimin
                  flarge  = 0.                   !--- Begin 6/19/03 changes
                  fsmall  = 1.
                  xsimass = rrho*massi(mdimin)   !--- End 6/19/03 changes
                  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).
   !
   !--- Begin 6/19/03 changes => allow NLImax to increase & FLARGE to 
   !    decrease at COLDER temperatures; set FLARGE to zero (i.e., only small
   !    ice) if the ice mixing ratio is below QI_min
 
!                 DUM    = MAX(0.05, MIN(1., EXP(.0536*TC)) )
                  dum    = max(0.05, min(1., exp(.0564*tc)) )
                  indexs = min(mdimax, max(mdimin, int(xmimax*dum) ) )
!
                  dum    = max(flgmin*pfac, dum)
 
                  qi_min = qi_min_0c * dum  !- Ver5    ----Not used ----
!!                QI_min = QI_min_0C        !- Ver5
!!!               QI_min = QI_min_0C/DUM    !- Ver5
 
                  nlimax = 10.e3/sqrt(dum)  !- Ver3
                  IF (tc < 0.) THEN
                    flarge = dum            !- Ver4
                  ELSE
                    flarge = 1.
                  ENDIF
                  fsmall  = (1.-flarge)/flarge
                  xsimass = rrho*massi(mdimin)*fsmall
                  tot_ice = thick*qi + blend*asnow
                  piloss  = -tot_ice/thick
                  lbef    = max(1,l-1)
                  rimef1  = (rimef_col(l)*thick*qi                      &
     &                    + rimef_col(lbef)*blend*asnow)/tot_ice
                  rimef1  = min(rimef1, rfmax)
                  vsnow   = 0.0
                  DO ipass=0,1
                    IF (rimef1 .LE. 1.) THEN
                      rimef1 = 1.
                      vrimef = 1.
                    ELSE
                      ixs  = max(2, min(indexs/100, 9))
                      xrf  = 10.492*log(rimef1)
                      ixrf = max(0, min(int(xrf), nrime))
                      IF (ixrf .GE. 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
                      IF(tc < 0) THEN
                        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
                        ENDIF             ! End IF (XLI <= MASSI(MDImin) ) 
                      ENDIF               ! End IF (TC < 0)
!
!--- Reduce excessive accumulation of ice at upper levels
!    associated with strong grid-resolved ascent
!
!--- Force NLICE to be between NLImin and NLImax
!
!--- 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.
!
!
 
                      dum = max(nlimin, min(nlimax, nlice) )
                      IF (dum >= nlimax .AND. indexs >= mdimax)         &
     &                 rimef1 = rho*(qtice/nlimax-xsimass)/massi(indexs)
!
!                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 <= EPSQ .AND. ASNOW <= CLIMIT)
              ENDIF          ! End IF (ICE_logical)
!
!----------------------------------------------------------------------
!--------------- Calculate individual processes -----------------------
!----------------------------------------------------------------------
!
!--- Cloud water autoconversion to rain and collection by rain
!
              IF (qw > epsq .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
                qw0   = qautx*rrho*xncw(l)
                praut = max(0., qw-qw0)*craut
                IF (qlice  >  epsq) THEN
      !
      !--- Collection of cloud water by large ice particles ("snow")
      !    PIACWI=PIACW for riming, PIACWI=0 for shedding
      !
!Moor              FWS   = MIN(1., CIACW*VEL_INC*NLICE*ACCRI(INDEXS)/PP**C1) ! 20050422
                   fws   = min(0.1, ciacw*vel_inc*nlice*accri(indexs)   &
     &                                                  /pp**c1)
                   piacw = fws*qw
                   IF (tc  < 0.) piacwi = piacw    ! Large ice riming
 
                ENDIF             ! End IF (QLICE > EPSQ)
              ENDIF               ! End IF (QW > EPSQ .AND. TC >= T_ICE)
!
!----------------------------------------------------------------------
!--- Loop around some of the ice-phase processes if no ice should be present
!----------------------------------------------------------------------
!
              IF (ice_logical) THEN
!
!--- Now the pretzel logic of calculating ice deposition
!
                IF (tc < t_ice .AND. (wv > qsigrd .OR. qw > epsq)) 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   = min(pp,fpvsi(dum1))          ! Updated (dummy) saturation vapor pressure w/r/t ice
                  dum   = rhgrd*eps*dum/(pp+epsm1*dum) ! Updated (dummy) saturation specific humidity 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_l*tk2
                  dwvi      = min(wvqw,qsw_l)-qsi_l
                  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
                    sfactor = sqrt(vel_inc)*schmit_fac
                    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_l .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-0.639)
                    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 > EPSQ))
!
!------------------------------------------------------------------------
!
              ENDIF     ! End of IF (ICE_logical)then loop
!
!------------------------------------------------------------------------
!
!--- Cloud water condensation
!
              IF (tc >= t_ice .AND. (qw > epsq .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 ! EndIF (DUM2 <  DUM1)
                ENDIF   ! EndIF (PIACWI == 0. .AND. PIDEP == 0.)
              ENDIF     ! EndIF (TC >= T_ICE .AND. (QW > EPSQ .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
                sfactor = sqrt(vel_inc)*schmit_fac
                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)
!               QSW0 = EPS*ESW0/(PP+epsm1*ESW0)
!!              dum  = min(PP, ESW0)
!!              QSW0 = EPS*dum/(PP+epsm1*dum)
!               DWV0 = MIN(WV,QSW)-QSW0
                dwv0 = min(wv,qsw_l)-qsw0_l
                dum  = qw + pcond
                IF (wv < qsw_l .AND. dum <= epsq) 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_l .AND. DUM <= EPSQ)
                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.
              indexr   = mdrmin
              indexr1  = indexr    ! For debugging only
              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) )
                  item   = cn0r0*sqrt(sqrt(rqr))               ! Moorthi 07/31/08
                  indexr = max( mdrmin, min(item, mdrmax) )    ! Moorthi 07/31/08
                ENDIF
   !
                IF (tc < t_ice) THEN
                  piacr = -prloss
                ELSE
                  dwvr = wv - pcond - qsw_l
                  dum  = qw + pcond
                  IF (dwvr < 0. .AND. dum <= epsq) 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
                    rfactor = sqrt(gammar)*schmit_fac
                    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 > epsq) THEN
                    fwr   = cracw*gammar*n0r*accrr(indexr)
!Moor               PRACW = MIN(1.,FWR)*QW               ! 20050422
                    pracw = min(0.1,fwr)*qw
                  ENDIF           ! End IF (DWVr < 0. .AND. DUM <= EPSQ)
!
                  IF (tc < 0. .AND. tcc < 0.) THEN
!
!--- Biggs (1953) heteorogeneous freezing (e.g., Lin et al., 1983)
!   - Rescaled mean drop diameter from microns (INDEXR) to mm (DUM) to prevent underflow
!
                    dum   = .001*float(indexr)
                    dum1  = dum * dum
                    dum   = (exp(abfr*tc)-1.)*dum1*dum1*dum1*dum
                    piacr = min(cbfr*n0r*rrho*dum, qtrain)
                    IF (qlice > epsq) 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 >  EPSQ)
                    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 <=  epsq) 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 <= epsq .AND. pidep <= epsq) 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 <= EPSQ .AND. PIDEP <= EPSQ)
                  qinew = tot_icenew/(thick+bldtrh*flimass*vsnow)
                  IF (qinew  <=  epsq) 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
   !      isn't 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 <= epsq) 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                ! End IF (TOT_RAINnew < CLIMIT)
!
              wcnew = qwnew + qrnew + qinew
!
!----------------------------------------------------------------------
!-------------- Begin debugging & verification ------------------------
!----------------------------------------------------------------------
!
!--- QT, QTnew - total water (vapor & condensate) before & after microphysics, resp.
!
!             QT=THICK*(QV+WC_col(l))+ARAIN+ASNOW
!             QTnew  = THICK*(WVnew/(1.+WVnew)+WCnew/(1.+wcnew))
!    &               + ARAINnew + ASNOWnew
 
              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 > TOLER) THEN
!               DUM=DUM/MIN(QT, QTnew)
!               IF (DUM > TOLER) DBG_logical=.TRUE.
!             ENDIF
!
!             DUM=(RHgrd+.001)*QSInew
!             IF ( (QWnew > EPSQ .OR. QRnew > EPSQ .OR. WVnew > DUM)
!     &           .AND. TC < T_ICE )  DBG_logical=.TRUE.
!
!             IF (TC > 5. .AND. QInewr > EPSQ) DBG_logical=.TRUE.
!
              IF ((wvnew < epsq .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    = min(pp, fpvsl(tnew))
!               QSWnew = EPS*ESW/(PP-ESW)
                qswnew = eps*esw/(pp+epsm1*esw)
                IF (tc < 0. .OR. tnew  <  0.) THEN
                  esi    = min(pp, fpvsi(tnew))
!                 QSInew = EPS*ESI/(PP-ESI)
                  qsinew = eps*esi/(pp+epsm1*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 > epsq) nstats(itdx,1) = nstats(itdx,1)+1
                IF (qinew > epsq .AND.  qrnew+qwnew > epsq)             &
     &            nstats(itdx,2) = nstats(itdx,2)+1
                IF (qwnew > epsq) nstats(itdx,3) = nstats(itdx,3)+1 
                IF (qrnew > epsq) 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                        ! temperature
!
!             QV_col(L)    = max(EPSQ, WVnew/(1.+WVnew)) ! specific humidity
              qv_col(l)    = max(0.0, wvnew           )  ! specific humidity
              wc_col(l)    = max(0.0, wcnew)             ! total condensate mixing ratio
              qi_col(l)    = max(0.0, qinew)             ! ice mixing ratio
              qr_col(l)    = max(0.0, qrnew)             ! rain mixing ratio
              qw_col(l)    = max(0.0, qwnew)             ! cloud water mixing ratio
              rimef_col(l) = rimef                       ! rime factor
              asnow        = asnownew                    ! accumulated snow
              arain        = arainnew                    ! accumulated rain
!
!#######################################################################
!
            ENDIF   ! End of IF (.NOT. CLEAR) THEN
          ENDIF     ! End of IF (QV_col(L) <= EPSQ .AND. WC_col(L) <= EPSQ) THEN
!
        ENDDO       ! ##### End "L" loop through model levels #####
!
        araing = araing + arain
        asnowg = asnowg + asnow
      enddo              ! do for ntimes=1,mic_step
!
!#######################################################################
!
!-----------------------------------------------------------------------
!--------------------------- Return to GSMDRIVE -----------------------
!-----------------------------------------------------------------------
!
      CONTAINS
!     END  SUBROUTINE GSMCOLUMN
!
!#######################################################################
!--------- Produces accurate calculation of cloud condensation ---------
!#######################################################################
!
      REAL FUNCTION CONDENSE (PP, QW, RHgrd, TK, WV)
!
      implicit none
!
!---------------------------------------------------------------------------------
!------ 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.).  ------
!---------------------------------------------------------------------------------
!
      real pp, qw, rhgrd, tk, wv
      INTEGER, PARAMETER :: HIGH_PRES=kind_phys
!     INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
      REAL (KIND=high_pres), PARAMETER ::                               &
     & rhlimit=.001, rhlimit1=-rhlimit
      REAL, PARAMETER :: RCP=1./cp, rcprv=rcp/rv
      REAL (KIND=high_pres) :: cond, ssat, wcdum, tsq
      real wvdum, tdum, xlv, xlv1, xlv2, ws, dwv, esw, rfac
!
!-----------------------------------------------------------------------
!
!--- 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      = min(pp, fpvsl(tdum))          ! Saturation vapor press w/r/t water
!     WS       = RHgrd*EPS*ESW/(PP-ESW)        ! Saturation mixing ratio w/r/t water
      ws       = rhgrd*eps*esw/(pp+epsm1*esw)  ! Saturation specific hum w/r/t water
      dwv      = wvdum - ws                    ! Deficit grid-scale specific humidity
      ssat     = dwv / ws                      ! Supersaturation ratio
      condense = 0.
      rfac     = 0.5                           ! converges faster with 0.5
      DO WHILE ((ssat < rhlimit1 .AND. wcdum > epsq)                    &
     &           .OR. ssat > rhlimit)
!
        xlv   = 3.148e6-2370.*tdum
        xlv1  = xlv*rcp
        xlv2  = xlv*xlv*rcprv
!
!       COND  = DWV/(1.+XLV2*WS/(Tdum*Tdum))  ! Asai (1965, J. Japan)
        tsq   = tdum*tdum
        cond  = rfac*dwv*tsq/(tsq+xlv2*ws)    ! Asai (1965, J. Japan)
!       COND  =      DWV*tsq/(tsq+XLV2*WS)    ! 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  = min(pp, fpvsl(tdum))           ! Updated saturation vapor press w/r/t water
!       WS   = RHgrd*EPS*ESW/(PP-ESW)         ! Updated saturation mixing ratio w/r/t water
        ws   = rhgrd*eps*esw/(pp+epsm1*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
        rfac = 1.0
      ENDDO
 
      END FUNCTION condense
!
!#######################################################################
!---------------- Calculate ice deposition at T<T_ICE ------------------
!#######################################################################
!
      REAL FUNCTION DEPOSIT (PP, RHgrd, Tdum, WVdum)
!
      implicit none
!
!--- Also uses the Asai (1965) algorithm, but uses a different target
!      vapor pressure for the adjustment
!
      REAL PP, RHgrd, Tdum, WVdum
      INTEGER, PARAMETER :: HIGH_PRES=kind_phys
!     INTEGER, PARAMETER :: HIGH_PRES=Selected_Real_Kind(15)
      REAL (KIND=high_pres), PARAMETER :: rhlimit=.001,                 & 
     & rhlimit1=-rhlimit
      REAL, PARAMETER :: RCP=1./cp, rcprv=rcp/rv, xls=hvap+hfus         &
     &,                  xls1=xls*rcp, xls2=xls*xls*rcprv
      REAL (KIND=high_pres) :: dep, ssat
      real esi, ws, dwv
!
!-----------------------------------------------------------------------
!
      esi=min(pp, fpvsi(tdum))                  ! Saturation vapor press w/r/t ice
!     WS=RHgrd*EPS*ESI/(PP-ESI)                 ! Saturation mixing ratio
      ws=rhgrd*eps*esi/(pp+epsm1*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=min(pp, fpvsi(tdum))                ! Updated saturation vapor press w/r/t ice
!       WS=RHgrd*EPS*ESI/(PP-ESI)               ! Updated saturation mixing ratio w/r/t ice
        ws=rhgrd*eps*esi/(pp+epsm1*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
      END FUNCTION deposit
!
      END  SUBROUTINE gsmcolumn
 
 
      SUBROUTINE rsipath(im, ix, ix2, levs, prsl, prsi, t, q, clw       &
     &,                  f_ice, f_rain, f_rime, flgmin                  &
     &,                  cwatp, cicep, rainp, snowp                     &
     &,                  recwat, rerain, resnow, lprnt, ipr)
!
      implicit none
!
!--------------------CLOUD----------------------------------------------
      integer im, ix, ix2, levs, ipr
      real    prsl(ix,levs), prsi(ix,levs+1), t(ix,levs), q(ix,levs)    &
     &,       clw(ix2,levs), f_ice(ix2,levs), f_rain(ix2,levs)          &
     &,       f_rime(ix2,levs)                                          &
     &,       cwatp(ix,levs), rainp(ix,levs),  cicep(ix,levs)           &
     &,       snowp(ix,levs), recwat(ix,levs), resnow(ix,levs)          &
     &,       rerain(ix,levs)
      real    flgmin
      real    frice, frrain, qcice, qcwat,  qrain, qsnow,   qtot, sden  &
     &,       cpath, rho,    dsnow, flarge, rimef, xsimass, nlice       &
     &,       tc,    recw1,  drain, xli,    dum,   NLImax, pfac, pp     &
     &,       snofac, tem
!
      real, parameter :: cexp=1./3.
      integer i, l, indexs
      logical lprnt
!
 
      recw1 = 620.3505 / tnw**cexp         ! cloud droplet effective radius
 
      do l=1,levs
        do i=1,im
                                           !--- HYDROMETEOR'S OPTICAL PATH
           cwatp(i,l) = 0.
           cicep(i,l) = 0.
           rainp(i,l) = 0.
           snowp(i,l) = 0.
                                           !--- HYDROMETEOR'S EFFECTIVE RADIUS
           recwat(i,l) = recwmin
           rerain(i,l) = rerainmin
           resnow(i,l) = resnowmin
        ENDDO
      ENDDO
 
      do l=1,levs
        DO i=1,im
 
   !        Assume parameterized condensate is 
   !         all water for T>=-10C,
   !         all ice for T<=-30C,
   !         and a linear mixture at -10C > T > -30C
   !
   !    * Determine hydrometeor composition of total condensate (QTOT)
   !
!         pp   = prsl(i,l) * 1000.0
          pp   = prsl(i,l) / prsi(i,levs+1)
!         pfac = max(0.25, sqrt(sqrt(min(1.0, pp*0.000025))))
!         pfac = max(0.5, sqrt(sqrt(min(1.0, pp*0.000025))))
!         pfac = max(0.5, sqrt(sqrt(min(1.0, pp*0.00002))))
!         pfac = max(0.25, sqrt(sqrt(min(1.0, pp*0.00001))))
!         pfac = max(0.25, sqrt(sqrt(min(1.0, pp))))
!         pfac = max(0.1, sqrt(min(1.0, pp*0.00001)))
!         pfac = max(0.5, sqrt(sqrt(min(1.0, pp*0.000033))))
!         pfac = max(0.5, sqrt(sqrt(min(1.0, pp*0.00004))))
!go       pfac = max(0.5, (sqrt(min(1.0, pp*0.000025))))
          pfac = 1.0
          tc   = t(i,l) - t0c
          qtot = clw(i,l)
          IF (qtot > epsq) THEN
             qcwat=0.
             qcice=0.
             qrain=0.
             qsnow=0.
             frice  = max(0.0, min(1.0, f_ice(i,l)))
             frrain = max(0.0, min(1.0, f_rain(i,l)))
             IF(tc <= thom) then
                qcice = qtot
             ELSE
                qcice = frice  * qtot
                qcwat = qtot   - qcice
                qrain = frrain * qcwat
                qcwat = qcwat  - qrain
             ENDIF
    !
    !--- Air density (RHO), model mass thickness (CPATH)
    !
             rho   = prsl(i,l)/(rd*t(i,l)*(1.+eps1*q(i,l)))
             cpath = (prsi(i,l+1)-prsi(i,l))*(1000000.0/grav)
 
    !! CLOUD WATER
    !
    !--- Effective radius (RECWAT) & total water path (CWATP)
    !    Assume monodisperse distribution of droplets (no factor of 1.5)
    !
             IF(qcwat > 0.) THEN
                recwat(i,l) = max(recwmin, recw1*(rho*qcwat)**cexp)
                cwatp(i,l)  = cpath*qcwat         ! cloud water path
!               tem         = 5.0*(1+max(0.0,min(1.0,-0.05*tc)))
!               RECWAT(I,L) = max(RECWAT(I,L), tem)
             ENDIF
 
    !! RAIN
    !
    !--- Effective radius (RERAIN) & total water path (RAINP)
    !--- Factor of 1.5 accounts for r**3/r**2 moments for exponentially
    !    distributed drops in effective radius calculations
    !    (from M.D. Chou's code provided to Y.-T. Hou)
    !
             IF(qrain > 0.) THEN
                drain       = cn0r0*sqrt(sqrt((rho*qrain)))
                rerain(i,l) = 1.5*max(xmrmin, min(drain, xmrmax))
                rainp(i,l)  = cpath*qrain         ! rain water path
             ENDIF
 
    !! SNOW (large ice) & CLOUD ICE
    !
    !--- Effective radius (RESNOW) & total ice path (SNOWP)
    !--- Total ice path (CICEP) for cloud ice 
    !--- Factor of 1.5 accounts for r**3/r**2 moments for exponentially
    !    distributed ice particles in effective radius calculations 
    !
    !--- Separation of cloud ice & "snow" uses algorithm from
    !    subroutine GSMCOLUMN
    !
             IF(qcice > 0.) THEN
    !
    !--- Mean particle size following Houze et al. (JAS, 1979, p. 160), 
    !    converted from Fig. 5 plot of LAMDAs.  An analogous set of
    !    relationships also shown by Fig. 8 of Ryan (BAMS, 1996, p. 66),
    !    but with a variety of different relationships that parallel the
    !    Houze curves.
    !
!               DUM=MAX(0.05, MIN(1., EXP(.0536*TC)) )
                dum=max(0.05, min(1., exp(.0564*tc)) )
                indexs=min(mdimax, max(mdimin, int(xmimax*dum) ) )
!               indexs=max(INDEXSmin, indexs)
!               NLImax=5.E3/sqrt(DUM)        !- Ver3
                dum=max(flgmin*pfac, dum)
!               DUM=MAX(FLGmin, DUM)
!               NLImax=20.E3            !- Ver3
!               NLImax=50.E3            !- Ver3 => comment this line out
                nlimax=10.e3/sqrt(dum)        !- Ver3
!               NLImax=5.E3/sqrt(DUM)        !- Ver3
!               NLImax=6.E3/sqrt(DUM)        !- Ver3
!               NLImax=7.5E3/sqrt(DUM)       !- Ver3
!               NLImax=20.E3/DUM        !- Ver3
!               NLImax=20.E3/max(0.2,DUM)        !- Ver3
!               NLImax=2.0E3/max(0.1,DUM)        !- Ver3
!               NLImax=2.5E3/max(0.1,DUM)        !- Ver3
!               NLImax=10.E3/max(0.2,DUM)        !- Ver3
!               NLImax=4.E3/max(0.2,DUM)        !- Ver3
!Moorthi        DSNOW  = XMImax*EXP(.0536*TC)
!Moorthi        INDEXS = MAX(INDEXSmin, MIN(MDImax, INT(DSNOW)))
    !
    !--- Assumed number fraction of large ice to total (large & small)
    !    ice particles, which is based on a general impression of the
    !    literature.
    !
    !    Small ice are assumed to have a mean diameter of 50 microns.
    !
                IF(tc >= 0.) THEN
                  flarge=flg1p0
                ELSE
                  flarge = dum
                ENDIF
!------------------------Commented by Moorthi -----------------------------
!               ELSEIF (TC >= -25.) THEN
!
!--- Note that absence of cloud water (QCWAT) is used as a quick
!    substitute for calculating water subsaturation as in GSMCOLUMN
!
!                  IF(QCWAT <= 0. .OR. TC < -8.
!    &                            .OR. TC > -3.)THEN
!                     FLARGE=FLG0P2
!                  ELSE
      
!--- Parameterize effects of rime splintering by increasing
!    number of small ice particles
!
!                     FLARGE=FLG0P1
!                  ENDIF
!               ELSEIF (TC <= -50.) THEN
!                  FLARGE=.01
!               ELSE
!                  FLARGE=.2*EXP(.1198*(TC+25))
!               ENDIF
!____________________________________________________________________________
 
                rimef=max(1., f_rime(i,l))
                xsimass=massi(mdimin)*(1.-flarge)/flarge
!     if (lprnt) print *,' rimef=',rimef,' xsimass=',xsimass
!    &,' indexs=',indexs,' massi=',massi(indexs),' flarge=',flarge
                nlice=rho*qcice/(xsimass+rimef*massi(indexs))
    !
    !--- From subroutine GSMCOLUMN:
    !--- Minimum number concentration for large ice of NLImin=10/m**3 
    !    at T>=0C.  Done in order to prevent unrealistically small 
    !    melting rates and tiny amounts of snow from falling to 
    !    unrealistically warm temperatures.
    !
                IF(tc >= 0.) THEN
                   nlice=max(nlimin, nlice)
                ELSEIF (nlice > nlimax) THEN
      !
      !--- Ferrier 6/13/01:  Prevent excess accumulation of ice
      !
                   xli=(rho*qcice/nlimax-xsimass)/rimef
 
                   IF(xli <= massi(450) ) THEN
                      dsnow=9.5885e5*xli**.42066
                   ELSE
                      dsnow=3.9751e6*xli**.49870
                   ENDIF
 
                   indexs=min(mdimax, max(indexs, int(dsnow)))
                   nlice=rho*qcice/(xsimass+rimef*massi(indexs))
                ENDIF
 
!               if (tc > -20.0 .and. indexs >= indexsmin) then
!                 snofac = max(0.0, min(1.0, exp(1.0*(tc+20.0))))
!               if (indexs >= indexsmin) then
!               if (tc > -20.0 .or. indexs >= indexsmin) then
!               if (tc > -40.0) then
!               if (tc >= -40.0 .or. prsl(i,l) > 50.0) then
!!              if (tc >= -20.0) then
!               if (tc >= -20.0 .or. prsl(i,l) > 50.0) then
!               if ((tc >= -20.0 .or.
!    &              prsi(i,levs+1)-prsi(i,l) < 30.0)
                if (prsi(i,levs+1)-prsi(i,l) < 40.0                     &
!               if (prsi(i,levs+1)-prsi(i,l) < 70.0
     &              .and. indexs >= indexsmin) then
!    &              prsi(i,levs)-prsl(i,l) < 20.0) then
!    &              prsi(i,levs)-prsl(i,l) < 30.0) then
!    &              prsi(i,levs)-prsl(i,l) < 40.0) then
!                 snofac = max(0.0, min(1.0, 0.05*(tc+40.0)))
!                 snofac = max(0.0, min(1.0, 0.1*(tc+25.0)))
!                 snofac = max(0.0, min(1.0, 0.0667*(tc+25.0)))
!               if (indexs > indexsmin) then
                  qsnow     = min(qcice, nlice*rimef*massi(indexs)/rho)
!    &                      * snofac
                endif
!               qsnow       = qcice
                qcice       = max(0., qcice-qsnow)
!               qsnow       = 0.0
                cicep(i,l) = cpath*qcice          ! cloud ice path
                resnow(i,l) = 1.5*float(indexs)
                sden        = sdens(indexs)/rimef  ! 1/snow density
                snowp(i,l) = cpath*qsnow*sden     ! snow path / snow density
!               SNOWP (I,L) = CPATH*QSNOW          ! snow path / snow density
!               if (lprnt .and. i == ipr) then
!                 print *,' L=',L,' snowp=',snowp(i,l),' cpath=',cpath
!    &,' qsnow=',qsnow,' sden=',sden,' rimef=',rimef,' indexs=',indexs
!    &,' sdens=',sdens(indexs),' resnow=',resnow(i,l)
!    &,' qcice=',qcice,' cicep=',cicep(i,l)
!               endif
 
                
             ENDIF                                 ! END QCICE BLOCK
           ENDIF                                   ! QTOT IF BLOCK
 
        ENDDO
      ENDDO
!
      END SUBROUTINE rsipath
 
 
 
!-----------------------------------
!\section gen_rsipath rsipath2 General Algorithm
      subroutine rsipath2                                               &
     &     ( plyr, plvl, tlyr, qlyr, qcwat, qcice, qrain, rrime,        &     ! inputs
     &       im, levs, iflip, flgmin,                                   &
     &       cwatp, cicep, rainp, snowp, recwat, rerain, resnow, snden  &     ! outputs
     &     )
 
! =================   subprogram documentation block   ================ !
!                                                                       !
! abstract:  this program is a modified version of ferrier's original   !
!   "rsipath" subprogram.  it computes layer's cloud liquid, ice, rain, !
!   and snow water condensate path and the partical effective radius    !
!   for liquid droplet, rain drop, and snow flake.                      !
!                                                                       !
!  ====================  defination of variables  ====================  !
!                                                                       !
! input variables:                                                      !
!   plyr  (IM,LEVS) : model layer mean pressure in mb (100Pa)           !
!   plvl  (IM,LEVS+1):model level pressure in mb (100Pa)                !
!   tlyr  (IM,LEVS) : model layer mean temperature in k                 !
!   qlyr  (IM,LEVS) : layer specific humidity in gm/gm                  !
!   qcwat (IM,LEVS) : layer cloud liquid water condensate amount        !
!   qcice (IM,LEVS) : layer cloud ice water condensate amount           !
!   qrain (IM,LEVS) : layer rain drop water amount                      !
!   rrime (IM,LEVS) : mass ratio of total to unrimed ice ( >= 1 )       !
!   IM              : horizontal dimention                              !
!   LEVS            : vertical layer dimensions                         !
!   iflip           : control flag for in/out vertical indexing         !
!                     =0: index from toa to surface                     !
!                     =1: index from surface to toa                     !
!   flgmin          : Minimum large ice fraction                        !
!   lprnt           : logical check print control flag                  !
!                                                                       !
! output variables:                                                     !
!   cwatp (IM,LEVS) : layer cloud liquid water path                     !
!   cicep (IM,LEVS) : layer cloud ice water path                        !
!   rainp (IM,LEVS) : layer rain water path                             !
!   snowp (IM,LEVS) : layer snow water path                             !
!   recwat(IM,LEVS) : layer cloud eff radius for liqid water (micron)   !
!   rerain(IM,LEVS) : layer rain water effective radius      (micron)   !
!   resnow(IM,LEVS) : layer snow flake effective radius      (micron)   !
!   snden (IM,LEVS) : 1/snow density                                    !
!                                                                       !
!                                                                       !
! usage:     call rsipath2                                              !
!                                                                       !
! subroutines called:  none                                             !
!                                                                       !
! program history log:                                                  !
!      xx-xx-2001   b. ferrier     - original program                   !
!      xx-xx-2004   s. moorthi     - modified for use in gfs model      !
!      05-20-2004   y. hou         - modified, added vertical index flag!
!                     to reduce data flipping, and rearrange code to    !
!                     be comformable with radiation part programs.      !
!                                                                       !
!  ====================    end of description    =====================  !
!
 
      implicit none
 
!  ---  constant parameter:
      real, parameter :: CEXP= 1.0/3.0
 
!  ---  inputs:
      real, dimension(:,:), intent(in) ::                               &
     &       plyr, plvl, tlyr, qlyr, qcwat, qcice, qrain, rrime
 
      integer, intent(in) :: IM, LEVS, iflip
      real, dimension(:),   intent(in) :: flgmin
!     logical, intent(in) :: lprnt
 
!  ---  output:
      real, dimension(:,:), intent(out) ::                              &
     &       cwatp, cicep, rainp, snowp, recwat, rerain, resnow, snden
 
!  ---  locals:
!     real,    dimension(IM,LEVS) :: delp, pp1, pp2
 
      real    :: recw1, dsnow, qsnow, qqcice, flarge, xsimass, pfac,    &
     &           nlice, xli, nlimax, dum, tem,                          &
     &           rho, cpath, rc, totcnd, tc
 
      integer :: i, k, indexs, ksfc, k1
!
!===>  ...  begin here
!
      recw1 = 620.3505 / tnw**cexp         ! cloud droplet effective radius
 
      do k = 1, levs
        do i = 1, im
                                           !--- hydrometeor's optical path
           cwatp(i,k) = 0.0
           cicep(i,k) = 0.0
           rainp(i,k) = 0.0
           snowp(i,k) = 0.0
           snden(i,k) = 0.0
                                           !--- hydrometeor's effective radius
           recwat(i,k) = recwmin
           rerain(i,k) = rerainmin
           resnow(i,k) = resnowmin
        enddo
      enddo
 
!  ---  set up pressure related arrays, convert unit from mb to cb (10Pa)
!       cause the rest part uses cb in computation
 
      if (iflip == 0) then        ! data from toa to sfc
        ksfc = levs + 1
        k1   = 0
      else                        ! data from sfc to top
        ksfc = 1
        k1   = 1
      endif                       ! end_if_iflip
!
      do k = 1, levs
        do i = 1, im
          totcnd = qcwat(i,k) + qcice(i,k) + qrain(i,k)
          qsnow = 0.0
          if(totcnd > epsq) then
 
!  ---  air density (rho), model mass thickness (cpath), temperature in c (tc)
 
            rho   = 0.1 * plyr(i,k)                                     &
     &            / (rd* tlyr(i,k) * (1.0 + eps1*qlyr(i,k)))
            cpath = abs(plvl(i,k+1) - plvl(i,k)) * (100000.0 / grav)
            tc    = tlyr(i,k) - t0c
 
!! cloud water
!
!  ---  effective radius (recwat) & total water path (cwatp):
!       assume monodisperse distribution of droplets (no factor of 1.5)
 
            if (qcwat(i,k) > 0.0) then
              recwat(i,k) = max(recwmin,recw1*(rho*qcwat(i,k))**cexp)
              cwatp(i,k) = cpath * qcwat(i,k)           ! cloud water path
!             tem         = 5.0*(1.0 + max(0.0, min(1.0,-0.05*tc)))
!             recwat(i,k) = max(recwat(i,k), tem)
            endif
 
!! rain
!
!  ---  effective radius (rerain) & total water path (rainp):
!       factor of 1.5 accounts for r**3/r**2 moments for exponentially
!       distributed drops in effective radius calculations
!       (from m.d. chou's code provided to y.-t. hou)
 
            if (qrain(i,k) > 0.0) then
              tem         = cn0r0 * sqrt(sqrt(rho*qrain(i,k)))
              rerain(i,k) = 1.5 * max(xmrmin, min(xmrmax, tem))
              rainp(i,k) = cpath * qrain(i,k)           ! rain water path
            endif
 
!! snow (large ice) & cloud ice
!
!  ---  effective radius (resnow) & total ice path (snowp) for snow, and
!       total ice path (cicep) for cloud ice:
!       factor of 1.5 accounts for r**3/r**2 moments for exponentially
!       distributed ice particles in effective radius calculations
!       separation of cloud ice & "snow" uses algorithm from subroutine gsmcolumn
 
!           pfac = max(0.5, sqrt(sqrt(min(1.0, pp1(i,k)*0.00004))))
!go         pfac = max(0.5, (sqrt(min(1.0, pp1(i,k)*0.000025))))
            pfac = 1.0
 
            if (qcice(i,k) > 0.0) then
 
!  ---  mean particle size following houze et al. (jas, 1979, p. 160),
!       converted from fig. 5 plot of lamdas.  an analogous set of
!       relationships also shown by fig. 8 of ryan (bams, 1996, p. 66),
!       but with a variety of different relationships that parallel
!       the houze curves.
 
!             dum = max(0.05, min(1.0, exp(0.0536*tc) ))
              dum = max(0.05, min(1.0, exp(0.0564*tc) ))
              indexs = min(mdimax, max(mdimin, int(xmimax*dum) ))
              dum=max(flgmin(i)*pfac, dum)
 
!  ---  assumed number fraction of large ice to total (large & small) ice
!       particles, which is based on a general impression of the literature.
!       small ice are assumed to have a mean diameter of 50 microns.
 
              if (tc >= 0.0) then
                flarge = flg1p0
              else
                flarge = dum
!               flarge = max(FLGmin*pfac, dum)
              endif
!------------------------commented by moorthi -----------------------------
!             elseif (tc >= -25.0) then
!
!  ---  note that absence of cloud water (qcwat) is used as a quick
!       substitute for calculating water subsaturation as in gsmcolumn
!
!               if (qcwat(i,k) <= 0.0 .or. tc < -8.0                 &
!    &                                .or. tc > -3.0) then
!                 flarge = FLG0P2
!               else
!
!  ---  parameterize effects of rime splintering by increasing
!       number of small ice particles
!
!                 flarge = FLG0P1
!               endif
!             elseif (tc <= -50.0) then
!               flarge = 0.01
!             else
!               flarge = 0.2 * exp(0.1198*(tc+25.0))
!             endif
!____________________________________________________________________________
 
              xsimass = massi(mdimin) * (1.0 - flarge) / flarge
!             nlimax = 20.0e3                                      !- ver3
!             NLImax=50.E3                 !- Ver3 => comment this line out
              nlimax=10.e3/sqrt(dum)       !- Ver3
!             NLImax=5.E3/sqrt(DUM)        !- Ver3
!             NLImax=6.E3/sqrt(DUM)        !- Ver3
!             NLImax=7.5E3/sqrt(DUM)       !- Ver3
 
!             indexs = min(MDImax, max(MDImin, int(XMImax*dum) ))
!moorthi      dsnow  = XMImax * exp(0.0536*tc)
!moorthi      indexs = max(INDEXSmin, min(MDImax, int(dsnow)))
 
!             if (lprnt) print *,' rrime=',rrime,' xsimass=',xsimass,   &
!    &       ' indexs=',indexs,' massi=',massi(indexs),' flarge=',flarge
 
              tem = rho * qcice(i,k)
              nlice = tem / (xsimass +rrime(i,k)*massi(indexs))
 
!  ---  from subroutine gsmcolumn:
!       minimum number concentration for large ice of NLImin=10/m**3
!       at t>=0c.  done in order to prevent unrealistically small
!       melting rates and tiny amounts of snow from falling to
!       unrealistically warm temperatures.
 
              if (tc >= 0.0) then
 
                nlice = max(nlimin, nlice)
 
              elseif (nlice > nlimax) then
 
!  ---  ferrier 6/13/01:  prevent excess accumulation of ice
 
                xli = (tem/nlimax - xsimass) / rrime(i,k)
 
                if (xli <= massi(450) ) then
                  dsnow = 9.5885e5 * xli**0.42066
                else
                  dsnow = 3.9751e6 * xli** 0.49870
                endif
 
                indexs = min(mdimax, max(indexs, int(dsnow)))
                nlice = tem / (xsimass + rrime(i,k)*massi(indexs))
 
              endif                               ! end if_tc block
 
!             if (abs(plvl(i,ksfc)-plvl(i,k+k1)) < 300.0                &
!             if (abs(plvl(i,ksfc)-plvl(i,k+k1)) < 400.0                &
!             if (plvl(i,k+k1) > 600.0                                  &
!    &                            .and. indexs >= INDEXSmin) then
!             if (tc > -20.0 .and. indexs >= indexsmin) then
              if (plvl(i,ksfc) > 850.0 .and.                            &
!    &            plvl(i,k+k1) > 600.0 .and. indexs >= indexsmin) then
     &            plvl(i,k+k1) > 700.0 .and. indexs >= indexsmin) then ! 20060516
!!            if (plvl(i,ksfc) > 800.0 .and.                            &
!!   &            plvl(i,k+k1) > 700.0 .and. indexs >= indexsmin) then
!             if (plvl(i,ksfc) > 700.0 .and.                            &
!    &            plvl(i,k+k1) > 600.0 .and. indexs >= indexsmin) then
                qsnow = min( qcice(i,k),                                &
     &                       nlice*rrime(i,k)*massi(indexs)/rho )
              endif
 
              qqcice      = max(0.0, qcice(i,k)-qsnow)
              cicep(i,k) = cpath * qqcice          ! cloud ice path
              resnow(i,k) = 1.5 * float(indexs)
              snden(i,k) = sdens(indexs) / rrime(i,k)   ! 1/snow density
              snowp(i,k) = cpath*qsnow             ! snow path
!             snowp (i,k) = cpath*qsnow*snden(i,k)  ! snow path / snow density
 
!             if (lprnt .and. i == ipr) then
!             if (i == 2) then
!               print *,' L=',k,' snowp=',snowp(i,k),' cpath=',cpath,   &
!    &         ' qsnow=',qsnow,' sden=',snden(i,k),' rrime=',rrime(i,k),&
!    &         ' indexs=',indexs,' sdens=',sdens(indexs),' resnow=',    &
!    &           resnow(i,k),' qcice=',qqcice,' cicep=',cicep(i,k)
!           endif
 
            endif                                 ! end if_qcice block
          endif                                   ! end if_totcnd block
 
        enddo
      enddo
!
!...................................
      end subroutine rsipath2
!-----------------------------------
 
      end MODULE module_microphysics