LCOV - code coverage report
Current view: | top level - physics/pumas-frozen - micro_mg1_0.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 0 | 1509 | 0.0 %
Date: | 2025-01-13 21:54:50 |  | Functions: | 0 | 3 | 0.0 %

              Line data    Source code
           1             : module micro_mg1_0
           2             :
           3             : !---------------------------------------------------------------------------------
           4             : ! Purpose:
           5             : !   MG microphysics
           6             : !
           7             : ! Author: Andrew Gettelman, Hugh Morrison.
           8             : ! Contributions from: Xiaohong Liu and Steve Ghan
           9             : ! December 2005-May 2010
          10             : ! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008)
          11             : !                 Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010)
          12             : ! for questions contact Hugh Morrison, Andrew Gettelman
          13             : ! e-mail: morrison@ucar.edu, andrew@ucar.edu
          14             : !
          15             : ! NOTE: Modified to allow other microphysics packages (e.g. CARMA) to do ice
          16             : ! microphysics in cooperation with the MG liquid microphysics. This is
          17             : ! controlled by the do_cldice variable.
          18             : !
          19             : ! NOTE: If do_cldice is false, then MG microphysics should not update CLDICE
          20             : ! or NUMICE; however, it is assumed that the other microphysics scheme will have
          21             : ! updated CLDICE and NUMICE. The other microphysics should handle the following
          22             : ! processes that would have been done by MG:
          23             : !   - Detrainment (liquid and ice)
          24             : !   - Homogeneous ice nucleation
          25             : !   - Heterogeneous ice nucleation
          26             : !   - Bergeron process
          27             : !   - Melting of ice
          28             : !   - Freezing of cloud drops
          29             : !   - Autoconversion (ice -> snow)
          30             : !   - Growth/Sublimation of ice
          31             : !   - Sedimentation of ice
          32             : !---------------------------------------------------------------------------------
          33             : ! modification for sub-columns, HM, (orig 8/11/10)
          34             : ! This is done using the logical 'microp_uniform' set to .true. = uniform for subcolumns
          35             : !---------------------------------------------------------------------------------
          36             :
          37             : ! Procedures required:
          38             : ! 1) An implementation of the gamma function (if not intrinsic).
          39             : ! 2) saturation vapor pressure to specific humidity formula
          40             : ! 3) svp over water
          41             : ! 4) svp over ice
          42             :
          43             : #ifndef HAVE_GAMMA_INTRINSICS
          44             : use shr_spfn_mod, only: gamma => shr_spfn_gamma
          45             : #endif
          46             :
          47             :   use wv_sat_methods, only: &
          48             :        svp_water => wv_sat_svp_water, &
          49             :        svp_ice => wv_sat_svp_ice, &
          50             :        svp_to_qsat => wv_sat_svp_to_qsat
          51             :
          52             :   use phys_control, only: phys_getopts
          53             :
          54             : implicit none
          55             : private
          56             : save
          57             :
          58             : ! Note: The liu_in option has been removed, as there was a serious bug with this
          59             : ! option being set to false. The code now behaves as if the default liu_in=.true.
          60             : ! is always on. Addition/reinstatement of ice nucleation options will likely be
          61             : ! done outside of this module.
          62             :
          63             : public :: &
          64             :      micro_mg_init, &
          65             :      micro_mg_get_cols, &
          66             :      micro_mg_tend
          67             :
          68             : integer, parameter :: r8 = selected_real_kind(12)      ! 8 byte real
          69             :
          70             : real(r8) :: g              !gravity
          71             : real(r8) :: r              !Dry air Gas constant
          72             : real(r8) :: rv             !water vapor gas contstant
          73             : real(r8) :: cpp            !specific heat of dry air
          74             : real(r8) :: rhow           !density of liquid water
          75             : real(r8) :: tmelt          ! Freezing point of water (K)
          76             : real(r8) :: xxlv           ! latent heat of vaporization
          77             : real(r8) :: xlf            !latent heat of freezing
          78             : real(r8) :: xxls           !latent heat of sublimation
          79             :
          80             : real(r8) :: rhosn  ! bulk density snow
          81             : real(r8) :: rhoi   ! bulk density ice
          82             :
          83             : real(r8) :: ac,bc,as,bs,ai,bi,ar,br  !fall speed parameters
          84             : real(r8) :: ci,di    !ice mass-diameter relation parameters
          85             : real(r8) :: cs,ds    !snow mass-diameter relation parameters
          86             : real(r8) :: cr,dr    !drop mass-diameter relation parameters
          87             : real(r8) :: f1s,f2s  !ventilation param for snow
          88             : real(r8) :: Eii      !collection efficiency aggregation of ice
          89             : real(r8) :: Ecr      !collection efficiency cloud droplets/rain
          90             : real(r8) :: f1r,f2r  !ventilation param for rain
          91             : real(r8) :: DCS      !autoconversion size threshold
          92             : real(r8) :: qsmall   !min mixing ratio
          93             : real(r8) :: bimm,aimm !immersion freezing
          94             : real(r8) :: rhosu     !typical 850mn air density
          95             : real(r8) :: mi0       ! new crystal mass
          96             : real(r8) :: rin       ! radius of contact nuclei
          97             : real(r8) :: pi       ! pi
          98             :
          99             : ! Additional constants to help speed up code
         100             :
         101             : real(r8) :: cons1
         102             : real(r8) :: cons4
         103             : real(r8) :: cons5
         104             : real(r8) :: cons6
         105             : real(r8) :: cons7
         106             : real(r8) :: cons8
         107             : real(r8) :: cons11
         108             : real(r8) :: cons13
         109             : real(r8) :: cons14
         110             : real(r8) :: cons16
         111             : real(r8) :: cons17
         112             : real(r8) :: cons22
         113             : real(r8) :: cons23
         114             : real(r8) :: cons24
         115             : real(r8) :: cons25
         116             : real(r8) :: cons27
         117             : real(r8) :: cons28
         118             :
         119             : real(r8) :: lammini
         120             : real(r8) :: lammaxi
         121             : real(r8) :: lamminr
         122             : real(r8) :: lammaxr
         123             : real(r8) :: lammins
         124             : real(r8) :: lammaxs
         125             :
         126             : ! parameters for snow/rain fraction for convective clouds
         127             : real(r8) :: tmax_fsnow ! max temperature for transition to convective snow
         128             : real(r8) :: tmin_fsnow ! min temperature for transition to convective snow
         129             :
         130             : !needed for findsp
         131             : real(r8) :: tt0       ! Freezing temperature
         132             :
         133             : real(r8) :: csmin,csmax,minrefl,mindbz
         134             :
         135             : real(r8) :: rhmini     ! Minimum rh for ice cloud fraction > 0.
         136             :
         137             : logical :: use_hetfrz_classnuc ! option to use heterogeneous freezing
         138             :
         139             : character(len=16)  :: micro_mg_precip_frac_method  ! type of precipitation fraction method
         140             : real(r8)           :: micro_mg_berg_eff_factor     ! berg efficiency factor
         141             :
         142             : ! Switches for specification rather than prediction of droplet and crystal number
         143             : ! note: number will be adjusted as needed to keep mean size within bounds,
         144             : ! even when specified droplet or ice number is used
         145             : !
         146             : ! If constant cloud ice number is set (nicons = .true.),
         147             : ! then all microphysical processes except mass transfer due to ice nucleation
         148             : ! (mnuccd) are based on the fixed cloud ice number. Calculation of
         149             : ! mnuccd follows from the prognosed ice crystal number ni.
         150             : logical :: nccons ! nccons=.true. to specify constant cloud droplet number
         151             : logical :: nicons ! nicons=.true. to specify constant cloud ice number
         152             :
         153             : ! parameters for specified ice and droplet number concentration
         154             : ! note: these are local in-cloud values, not grid-mean
         155             : real(r8) :: ncnst ! droplet num concentration when nccons=.true. (m-3)
         156             : real(r8) :: ninst ! ice num concentration when nicons=.true. (m-3)
         157             :
         158             : !===============================================================================
         159             : contains
         160             : !===============================================================================
         161             :
         162           0 : subroutine micro_mg_init( &
         163             :      kind, gravit, rair, rh2o, cpair,  &
         164             :      rhoh2o, tmelt_in, latvap, latice, &
         165             :      rhmini_in, micro_mg_dcs, use_hetfrz_classnuc_in, &
         166           0 :      micro_mg_precip_frac_method_in, micro_mg_berg_eff_factor_in, &
         167           0 :      nccons_in, nicons_in, ncnst_in, ninst_in, errstring)
         168             :
         169             : !-----------------------------------------------------------------------
         170             : !
         171             : ! Purpose:
         172             : ! initialize constants for the morrison microphysics
         173             : !
         174             : ! Author: Andrew Gettelman Dec 2005
         175             : !
         176             : !-----------------------------------------------------------------------
         177             :
         178             : integer,          intent(in)  :: kind            ! Kind used for reals
         179             : real(r8),         intent(in)  :: gravit
         180             : real(r8),         intent(in)  :: rair
         181             : real(r8),         intent(in)  :: rh2o
         182             : real(r8),         intent(in)  :: cpair
         183             : real(r8),         intent(in)  :: rhoh2o
         184             : real(r8),         intent(in)  :: tmelt_in        ! Freezing point of water (K)
         185             : real(r8),         intent(in)  :: latvap
         186             : real(r8),         intent(in)  :: latice
         187             : real(r8),         intent(in)  :: rhmini_in       ! Minimum rh for ice cloud fraction > 0.
         188             : real(r8),         intent(in)  :: micro_mg_dcs
         189             : logical,          intent(in)  :: use_hetfrz_classnuc_in
         190             : character(len=16),intent(in)  :: micro_mg_precip_frac_method_in  ! type of precipitation fraction method
         191             : real(r8),         intent(in)  :: micro_mg_berg_eff_factor_in     ! berg efficiency factor
         192             : logical,          intent(in)  :: nccons_in
         193             : logical,          intent(in)  :: nicons_in
         194             : real(r8),         intent(in)  :: ncnst_in
         195             : real(r8),         intent(in)  :: ninst_in
         196             :
         197             : character(128),   intent(out) :: errstring       ! Output status (non-blank for error return)
         198             :
         199             : integer k
         200             :
         201             : integer l,m, iaer
         202             : real(r8) surften       ! surface tension of water w/respect to air (N/m)
         203             : real(r8) arg
         204             : !-----------------------------------------------------------------------
         205             :
         206           0 : errstring = ' '
         207             :
         208           0 : if( kind .ne. r8 ) then
         209           0 :    errstring = 'micro_mg_init: KIND of reals does not match'
         210           0 :    return
         211             : end if
         212             :
         213             : !declarations for morrison codes (transforms variable names)
         214             :
         215           0 : g= gravit                  !gravity
         216           0 : r= rair                    !Dry air Gas constant: note units(phys_constants are in J/K/kmol)
         217           0 : rv= rh2o                   !water vapor gas contstant
         218           0 : cpp = cpair                !specific heat of dry air
         219           0 : rhow = rhoh2o              !density of liquid water
         220           0 : tmelt = tmelt_in
         221           0 : rhmini = rhmini_in
         222           0 : micro_mg_precip_frac_method = micro_mg_precip_frac_method_in
         223           0 : micro_mg_berg_eff_factor    = micro_mg_berg_eff_factor_in
         224             :
         225           0 : nccons = nccons_in
         226           0 : nicons = nicons_in
         227           0 : ncnst  = ncnst_in
         228           0 : ninst  = ninst_in
         229             :
         230             : ! latent heats
         231             :
         232           0 : xxlv = latvap         ! latent heat vaporization
         233           0 : xlf = latice          ! latent heat freezing
         234           0 : xxls = xxlv + xlf     ! latent heat of sublimation
         235             :
         236             : ! flags
         237           0 : use_hetfrz_classnuc = use_hetfrz_classnuc_in
         238             :
         239             : ! parameters for snow/rain fraction for convective clouds
         240             :
         241           0 : tmax_fsnow = tmelt
         242           0 : tmin_fsnow = tmelt-5._r8
         243             :
         244             : ! parameters below from Reisner et al. (1998)
         245             : ! density parameters (kg/m3)
         246             :
         247           0 : rhosn = 250._r8    ! bulk density snow  (++ ceh)
         248           0 : rhoi = 500._r8     ! bulk density ice
         249           0 : rhow = 1000._r8    ! bulk density liquid
         250             :
         251             :
         252             : ! fall speed parameters, V = aD^b
         253             : ! V is in m/s
         254             :
         255             : ! droplets
         256           0 : ac = 3.e7_r8
         257           0 : bc = 2._r8
         258             :
         259             : ! snow
         260           0 : as = 11.72_r8
         261           0 : bs = 0.41_r8
         262             :
         263             : ! cloud ice
         264           0 : ai = 700._r8
         265           0 : bi = 1._r8
         266             :
         267             : ! rain
         268           0 : ar = 841.99667_r8
         269           0 : br = 0.8_r8
         270             :
         271             : ! particle mass-diameter relationship
         272             : ! currently we assume spherical particles for cloud ice/snow
         273             : ! m = cD^d
         274             :
         275           0 : pi= 3.1415927_r8
         276             :
         277             : ! cloud ice mass-diameter relationship
         278             :
         279           0 : ci = rhoi*pi/6._r8
         280           0 : di = 3._r8
         281             :
         282             : ! snow mass-diameter relationship
         283             :
         284           0 : cs = rhosn*pi/6._r8
         285           0 : ds = 3._r8
         286             :
         287             : ! drop mass-diameter relationship
         288             :
         289           0 : cr = rhow*pi/6._r8
         290           0 : dr = 3._r8
         291             :
         292             : ! ventilation parameters for snow
         293             : ! hall and prupacher
         294             :
         295           0 : f1s = 0.86_r8
         296           0 : f2s = 0.28_r8
         297             :
         298             : ! collection efficiency, aggregation of cloud ice and snow
         299             :
         300           0 : Eii = 0.1_r8
         301             :
         302             : ! collection efficiency, accretion of cloud water by rain
         303             :
         304           0 : Ecr = 1.0_r8
         305             :
         306             : ! ventilation constants for rain
         307             :
         308           0 : f1r = 0.78_r8
         309           0 : f2r = 0.32_r8
         310             :
         311             : ! autoconversion size threshold for cloud ice to snow (m)
         312             :
         313           0 : Dcs = micro_mg_dcs
         314             :
         315             : ! smallest mixing ratio considered in microphysics
         316             :
         317           0 : qsmall = 1.e-18_r8
         318             :
         319             : ! immersion freezing parameters, bigg 1953
         320             :
         321           0 : bimm = 100._r8
         322           0 : aimm = 0.66_r8
         323             :
         324             : ! typical air density at 850 mb
         325             :
         326           0 : rhosu = 85000._r8/(rair * tmelt)
         327             :
         328             : ! mass of new crystal due to aerosol freezing and growth (kg)
         329             :
         330           0 : mi0 = 4._r8/3._r8*pi*rhoi*(10.e-6_r8)*(10.e-6_r8)*(10.e-6_r8)
         331             :
         332             : ! radius of contact nuclei aerosol (m)
         333             :
         334           0 : rin = 0.1e-6_r8
         335             :
         336             : ! freezing temperature
         337           0 : tt0=273.15_r8
         338             :
         339           0 : pi=4._r8*atan(1.0_r8)
         340             :
         341             : !Range of cloudsat reflectivities (dBz) for analytic simulator
         342           0 : csmin= -30._r8
         343           0 : csmax= 26._r8
         344           0 : mindbz = -99._r8
         345             : !      minrefl = 10._r8**(mindbz/10._r8)
         346           0 : minrefl = 1.26e-10_r8
         347             :
         348             : ! Define constants to help speed up code (limit calls to gamma function)
         349             :
         350           0 : cons1=gamma(1._r8+di)
         351           0 : cons4=gamma(1._r8+br)
         352           0 : cons5=gamma(4._r8+br)
         353           0 : cons6=gamma(1._r8+ds)
         354           0 : cons7=gamma(1._r8+bs)
         355           0 : cons8=gamma(4._r8+bs)
         356           0 : cons11=gamma(3._r8+bs)
         357           0 : cons13=gamma(5._r8/2._r8+br/2._r8)
         358           0 : cons14=gamma(5._r8/2._r8+bs/2._r8)
         359           0 : cons16=gamma(1._r8+bi)
         360           0 : cons17=gamma(4._r8+bi)
         361           0 : cons22=(4._r8/3._r8*pi*rhow*(25.e-6_r8)**3)
         362           0 : cons23=dcs**3
         363           0 : cons24=dcs**2
         364           0 : cons25=dcs**bs
         365           0 : cons27=xxlv**2
         366           0 : cons28=xxls**2
         367             :
         368           0 : lammaxi = 1._r8/10.e-6_r8
         369           0 : lammini = 1._r8/(2._r8*dcs)
         370           0 : lammaxr = 1._r8/20.e-6_r8
         371           0 : lamminr = 1._r8/500.e-6_r8
         372           0 : lammaxs = 1._r8/10.e-6_r8
         373           0 : lammins = 1._r8/2000.e-6_r8
         374             :
         375             : end subroutine micro_mg_init
         376             :
         377             : !===============================================================================
         378             : !microphysics routine for each timestep goes here...
         379             :
         380           0 : subroutine micro_mg_tend ( &
         381             :      microp_uniform, pcols, pver, ncol, top_lev, deltatin,&
         382           0 :      tn, qn, qc, qi, nc,                              &
         383           0 :      ni, p, pdel, cldn, liqcldf,                      &
         384           0 :      relvar, accre_enhan,                             &
         385           0 :      icecldf, rate1ord_cw2pr_st, naai, npccnin,       &
         386           0 :      rndst, nacon, tlat, qvlat, qctend,               &
         387           0 :      qitend, nctend, nitend, effc, effc_fn,           &
         388           0 :      effi, prect, preci, nevapr, evapsnow, am_evp_st, &
         389           0 :      prain, prodsnow, cmeout, deffi, pgamrad,         &
         390           0 :      lamcrad, qsout, dsout, rflx, sflx,               &
         391           0 :      qrout, reff_rain, reff_snow, qcsevap, qisevap,   &
         392           0 :      qvres, cmeiout, vtrmc, vtrmi, qcsedten,          &
         393           0 :      qisedten, prao, prco, mnuccco, mnuccto,          &
         394           0 :      msacwio, psacwso, bergso, bergo, melto,          &
         395           0 :      homoo, qcreso, prcio, praio, qireso,             &
         396           0 :      mnuccro, pracso, meltsdt, frzrdt, mnuccdo,       &
         397           0 :      nrout, nsout, refl, arefl, areflz,               &
         398           0 :      frefl, csrfl, acsrfl, fcsrfl, rercld,            &
         399           0 :      ncai, ncal, qrout2, qsout2, nrout2,              &
         400           0 :      nsout2, drout2, dsout2, freqs, freqr,            &
         401           0 :      nfice, prer_evap, do_cldice, errstring,          &
         402           0 :      tnd_qsnow, tnd_nsnow, re_ice,                    &
         403           0 :      frzimm, frzcnt, frzdep)
         404             :
         405             : ! input arguments
         406             : logical,  intent(in) :: microp_uniform  ! True = configure uniform for sub-columns  False = use w/o sub-columns (standard)
         407             : integer,  intent(in) :: pcols                ! size of column (first) index
         408             : integer,  intent(in) :: pver                 ! number of layers in columns
         409             : integer,  intent(in) :: ncol                 ! number of columns
         410             : integer,  intent(in) :: top_lev              ! top level microphys is applied
         411             : real(r8), intent(in) :: deltatin             ! time step (s)
         412             : real(r8), intent(in) :: tn(pcols,pver)       ! input temperature (K)
         413             : real(r8), intent(in) :: qn(pcols,pver)       ! input h20 vapor mixing ratio (kg/kg)
         414             : real(r8), intent(in) :: relvar(pcols,pver)   ! relative variance of cloud water (-)
         415             : real(r8), intent(in) :: accre_enhan(pcols,pver) ! optional accretion enhancement factor (-)
         416             :
         417             : ! note: all input cloud variables are grid-averaged
         418             : real(r8), intent(inout) :: qc(pcols,pver)    ! cloud water mixing ratio (kg/kg)
         419             : real(r8), intent(inout) :: qi(pcols,pver)    ! cloud ice mixing ratio (kg/kg)
         420             : real(r8), intent(inout) :: nc(pcols,pver)    ! cloud water number conc (1/kg)
         421             : real(r8), intent(inout) :: ni(pcols,pver)    ! cloud ice number conc (1/kg)
         422             : real(r8), intent(in) :: p(pcols,pver)        ! air pressure (pa)
         423             : real(r8), intent(in) :: pdel(pcols,pver)     ! pressure difference across level (pa)
         424             : real(r8), intent(in) :: cldn(pcols,pver)     ! cloud fraction
         425             : real(r8), intent(in) :: icecldf(pcols,pver)  ! ice cloud fraction
         426             : real(r8), intent(in) :: liqcldf(pcols,pver)  ! liquid cloud fraction
         427             :
         428             : real(r8), intent(out) :: rate1ord_cw2pr_st(pcols,pver) ! 1st order rate for direct cw to precip conversion
         429             : ! used for scavenging
         430             : ! Inputs for aerosol activation
         431             : real(r8), intent(in) :: naai(pcols,pver)      ! ice nulceation number (from microp_aero_ts)
         432             : real(r8), intent(in) :: npccnin(pcols,pver)   ! ccn activated number tendency (from microp_aero_ts)
         433             : real(r8), intent(in) :: rndst(pcols,pver,4)   ! radius of 4 dust bins for contact freezing (from microp_aero_ts)
         434             : real(r8), intent(in) :: nacon(pcols,pver,4)   ! number in 4 dust bins for contact freezing  (from microp_aero_ts)
         435             :
         436             : ! Used with CARMA cirrus microphysics
         437             : ! (or similar external microphysics model)
         438             : logical,  intent(in) :: do_cldice             ! Prognosing cldice
         439             :
         440             : ! output arguments
         441             :
         442             : real(r8), intent(out) :: tlat(pcols,pver)    ! latent heating rate       (W/kg)
         443             : real(r8), intent(out) :: qvlat(pcols,pver)   ! microphysical tendency qv (1/s)
         444             : real(r8), intent(out) :: qctend(pcols,pver)  ! microphysical tendency qc (1/s)
         445             : real(r8), intent(out) :: qitend(pcols,pver)  ! microphysical tendency qi (1/s)
         446             : real(r8), intent(out) :: nctend(pcols,pver)  ! microphysical tendency nc (1/(kg*s))
         447             : real(r8), intent(out) :: nitend(pcols,pver)  ! microphysical tendency ni (1/(kg*s))
         448             : real(r8), intent(out) :: effc(pcols,pver)    ! droplet effective radius (micron)
         449             : real(r8), intent(out) :: effc_fn(pcols,pver) ! droplet effective radius, assuming nc = 1.e8 kg-1
         450             : real(r8), intent(out) :: effi(pcols,pver)    ! cloud ice effective radius (micron)
         451             : real(r8), intent(out) :: prect(pcols)        ! surface precip rate (m/s)
         452             : real(r8), intent(out) :: preci(pcols)        ! cloud ice/snow precip rate (m/s)
         453             : real(r8), intent(out) :: nevapr(pcols,pver)  ! evaporation rate of rain + snow
         454             : real(r8), intent(out) :: evapsnow(pcols,pver)! sublimation rate of snow
         455             : real(r8), intent(out) :: am_evp_st(pcols,pver)! stratiform evaporation area
         456             : real(r8), intent(out) :: prain(pcols,pver)   ! production of rain + snow
         457             : real(r8), intent(out) :: prodsnow(pcols,pver)! production of snow
         458             : real(r8), intent(out) :: cmeout(pcols,pver)  ! evap/sub of cloud
         459             : real(r8), intent(out) :: deffi(pcols,pver)   ! ice effective diameter for optics (radiation)
         460             : real(r8), intent(out) :: pgamrad(pcols,pver) ! ice gamma parameter for optics (radiation)
         461             : real(r8), intent(out) :: lamcrad(pcols,pver) ! slope of droplet distribution for optics (radiation)
         462             : real(r8), intent(out) :: qsout(pcols,pver)   ! snow mixing ratio (kg/kg)
         463             : real(r8), intent(out) :: dsout(pcols,pver)   ! snow diameter (m)
         464             : real(r8), intent(out) :: rflx(pcols,pver+1)  ! grid-box average rain flux (kg m^-2 s^-1)
         465             : real(r8), intent(out) :: sflx(pcols,pver+1)  ! grid-box average snow flux (kg m^-2 s^-1)
         466             : real(r8), intent(out) :: qrout(pcols,pver)     ! grid-box average rain mixing ratio (kg/kg)
         467             : real(r8), intent(inout) :: reff_rain(pcols,pver) ! rain effective radius (micron)
         468             : real(r8), intent(inout) :: reff_snow(pcols,pver) ! snow effective radius (micron)
         469             : real(r8), intent(out) :: qcsevap(pcols,pver) ! cloud water evaporation due to sedimentation
         470             : real(r8), intent(out) :: qisevap(pcols,pver) ! cloud ice sublimation due to sublimation
         471             : real(r8), intent(out) :: qvres(pcols,pver) ! residual condensation term to ensure RH < 100%
         472             : real(r8), intent(out) :: cmeiout(pcols,pver) ! grid-mean cloud ice sub/dep
         473             : real(r8), intent(out) :: vtrmc(pcols,pver) ! mass-weighted cloud water fallspeed
         474             : real(r8), intent(out) :: vtrmi(pcols,pver) ! mass-weighted cloud ice fallspeed
         475             : real(r8), intent(out) :: qcsedten(pcols,pver) ! qc sedimentation tendency
         476             : real(r8), intent(out) :: qisedten(pcols,pver) ! qi sedimentation tendency
         477             : ! microphysical process rates for output (mixing ratio tendencies)
         478             : real(r8), intent(out) :: prao(pcols,pver) ! accretion of cloud by rain
         479             : real(r8), intent(out) :: prco(pcols,pver) ! autoconversion of cloud to rain
         480             : real(r8), intent(out) :: mnuccco(pcols,pver) ! mixing rat tend due to immersion freezing
         481             : real(r8), intent(out) :: mnuccto(pcols,pver) ! mixing ratio tend due to contact freezing
         482             : real(r8), intent(out) :: msacwio(pcols,pver) ! mixing ratio tend due to H-M splintering
         483             : real(r8), intent(out) :: psacwso(pcols,pver) ! collection of cloud water by snow
         484             : real(r8), intent(out) :: bergso(pcols,pver) ! bergeron process on snow
         485             : real(r8), intent(out) :: bergo(pcols,pver) ! bergeron process on cloud ice
         486             : real(r8), intent(out) :: melto(pcols,pver) ! melting of cloud ice
         487             : real(r8), intent(out) :: homoo(pcols,pver) ! homogeneos freezign cloud water
         488             : real(r8), intent(out) :: qcreso(pcols,pver) ! residual cloud condensation due to removal of excess supersat
         489             : real(r8), intent(out) :: prcio(pcols,pver) ! autoconversion of cloud ice to snow
         490             : real(r8), intent(out) :: praio(pcols,pver) ! accretion of cloud ice by snow
         491             : real(r8), intent(out) :: qireso(pcols,pver) ! residual ice deposition due to removal of excess supersat
         492             : real(r8), intent(out) :: mnuccro(pcols,pver) ! mixing ratio tendency due to heterogeneous freezing of rain to snow (1/s)
         493             : real(r8), intent(out) :: pracso (pcols,pver) ! mixing ratio tendency due to accretion of rain by snow (1/s)
         494             : real(r8), intent(out) :: meltsdt(pcols,pver) ! latent heating rate due to melting of snow  (W/kg)
         495             : real(r8), intent(out) :: frzrdt (pcols,pver) ! latent heating rate due to homogeneous freezing of rain (W/kg)
         496             : real(r8), intent(out) :: mnuccdo(pcols,pver) ! mass tendency from ice nucleation
         497             : real(r8), intent(out) :: nrout(pcols,pver) ! rain number concentration (1/m3)
         498             : real(r8), intent(out) :: nsout(pcols,pver) ! snow number concentration (1/m3)
         499             : real(r8), intent(out) :: refl(pcols,pver)    ! analytic radar reflectivity
         500             : real(r8), intent(out) :: arefl(pcols,pver)  !average reflectivity will zero points outside valid range
         501             : real(r8), intent(out) :: areflz(pcols,pver)  !average reflectivity in z.
         502             : real(r8), intent(out) :: frefl(pcols,pver)
         503             : real(r8), intent(out) :: csrfl(pcols,pver)   !cloudsat reflectivity
         504             : real(r8), intent(out) :: acsrfl(pcols,pver)  !cloudsat average
         505             : real(r8), intent(out) :: fcsrfl(pcols,pver)
         506             : real(r8), intent(out) :: rercld(pcols,pver) ! effective radius calculation for rain + cloud
         507             : real(r8), intent(out) :: ncai(pcols,pver) ! output number conc of ice nuclei available (1/m3)
         508             : real(r8), intent(out) :: ncal(pcols,pver) ! output number conc of CCN (1/m3)
         509             : real(r8), intent(out) :: qrout2(pcols,pver)
         510             : real(r8), intent(out) :: qsout2(pcols,pver)
         511             : real(r8), intent(out) :: nrout2(pcols,pver)
         512             : real(r8), intent(out) :: nsout2(pcols,pver)
         513             : real(r8), intent(out) :: drout2(pcols,pver) ! mean rain particle diameter (m)
         514             : real(r8), intent(out) :: dsout2(pcols,pver) ! mean snow particle diameter (m)
         515             : real(r8), intent(out) :: freqs(pcols,pver)
         516             : real(r8), intent(out) :: freqr(pcols,pver)
         517             : real(r8), intent(out) :: nfice(pcols,pver)
         518             : real(r8), intent(out) :: prer_evap(pcols,pver)
         519             :
         520           0 : real(r8) :: nevapr2(pcols,pver)
         521             :
         522             : character(128),   intent(out) :: errstring       ! Output status (non-blank for error return)
         523             :
         524             : ! Tendencies calculated by external schemes that can replace MG's native
         525             : ! process tendencies.
         526             :
         527             : ! Used with CARMA cirrus microphysics
         528             : ! (or similar external microphysics model)
         529             : real(r8), intent(in) :: tnd_qsnow(:,:) ! snow mass tendency (kg/kg/s)
         530             : real(r8), intent(in) :: tnd_nsnow(:,:) ! snow number tendency (#/kg/s)
         531             : real(r8), intent(in) :: re_ice(:,:)    ! ice effective radius (m)
         532             :
         533             : ! From external ice nucleation.
         534             : real(r8), intent(in) :: frzimm(:,:) ! Number tendency due to immersion freezing (1/cm3)
         535             : real(r8), intent(in) :: frzcnt(:,:) ! Number tendency due to contact freezing (1/cm3)
         536             : real(r8), intent(in) :: frzdep(:,:) ! Number tendency due to deposition nucleation (1/cm3)
         537             :
         538             : ! local workspace
         539             : ! all units mks unless otherwise stated
         540             :
         541             : ! Additional constants to help speed up code
         542             : real(r8) :: cons2
         543             : real(r8) :: cons3
         544             : real(r8) :: cons9
         545             : real(r8) :: cons10
         546             : real(r8) :: cons12
         547             : real(r8) :: cons15
         548             : real(r8) :: cons18
         549             : real(r8) :: cons19
         550             : real(r8) :: cons20
         551             :
         552             : ! temporary variables for sub-stepping
         553           0 : real(r8) :: t1(pcols,pver)
         554           0 : real(r8) :: q1(pcols,pver)
         555           0 : real(r8) :: qc1(pcols,pver)
         556           0 : real(r8) :: qi1(pcols,pver)
         557           0 : real(r8) :: nc1(pcols,pver)
         558           0 : real(r8) :: ni1(pcols,pver)
         559           0 : real(r8) :: tlat1(pcols,pver)
         560           0 : real(r8) :: qvlat1(pcols,pver)
         561           0 : real(r8) :: qctend1(pcols,pver)
         562           0 : real(r8) :: qitend1(pcols,pver)
         563           0 : real(r8) :: nctend1(pcols,pver)
         564           0 : real(r8) :: nitend1(pcols,pver)
         565           0 : real(r8) :: prect1(pcols)
         566           0 : real(r8) :: preci1(pcols)
         567             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         568             :
         569             : real(r8) :: deltat        ! sub-time step (s)
         570             : real(r8) :: omsm    ! number near unity for round-off issues
         571             : real(r8) :: dto2    ! dt/2 (s)
         572             : real(r8) :: mincld  ! minimum allowed cloud fraction
         573           0 : real(r8) :: q(pcols,pver) ! water vapor mixing ratio (kg/kg)
         574           0 : real(r8) :: t(pcols,pver) ! temperature (K)
         575           0 : real(r8) :: rho(pcols,pver) ! air density (kg m-3)
         576           0 : real(r8) :: dv(pcols,pver)  ! diffusivity of water vapor in air
         577           0 : real(r8) :: mu(pcols,pver)  ! viscocity of air
         578           0 : real(r8) :: sc(pcols,pver)  ! schmidt number
         579           0 : real(r8) :: kap(pcols,pver) ! thermal conductivity of air
         580           0 : real(r8) :: rhof(pcols,pver) ! air density correction factor for fallspeed
         581           0 : real(r8) :: cldmax(pcols,pver) ! precip fraction assuming maximum overlap
         582           0 : real(r8) :: cldm(pcols,pver)   ! cloud fraction
         583           0 : real(r8) :: icldm(pcols,pver)   ! ice cloud fraction
         584           0 : real(r8) :: lcldm(pcols,pver)   ! liq cloud fraction
         585             : real(r8) :: icwc(pcols)    ! in cloud water content (liquid+ice)
         586             : real(r8) :: calpha(pcols)  ! parameter for cond/evap (Zhang et al. 2003)
         587             : real(r8) :: cbeta(pcols) ! parameter for cond/evap (Zhang et al. 2003)
         588             : real(r8) :: cbetah(pcols) ! parameter for cond/evap (Zhang et al. 2003)
         589             : real(r8) :: cgamma(pcols) ! parameter for cond/evap (Zhang et al. 2003)
         590             : real(r8) :: cgamah(pcols) ! parameter for cond/evap (Zhang et al. 2003)
         591             : real(r8) :: rcgama(pcols) ! parameter for cond/evap (Zhang et al. 2003)
         592             : real(r8) :: cmec1(pcols) ! parameter for cond/evap (Zhang et al. 2003)
         593             : real(r8) :: cmec2(pcols) ! parameter for cond/evap (Zhang et al. 2003)
         594             : real(r8) :: cmec3(pcols) ! parameter for cond/evap (Zhang et al. 2003)
         595             : real(r8) :: cmec4(pcols) ! parameter for cond/evap (Zhang et al. 2003)
         596             : real(r8) :: qtmp ! dummy qv
         597             : real(r8) :: dum  ! temporary dummy variable
         598             :
         599           0 : real(r8) :: cme(pcols,pver)  ! total (liquid+ice) cond/evap rate of cloud
         600             :
         601           0 : real(r8) :: cmei(pcols,pver) ! dep/sublimation rate of cloud ice
         602           0 : real(r8) :: cwml(pcols,pver) ! cloud water mixing ratio
         603           0 : real(r8) :: cwmi(pcols,pver) ! cloud ice mixing ratio
         604           0 : real(r8) :: nnuccd(pver)   ! ice nucleation rate from deposition/cond.-freezing
         605           0 : real(r8) :: mnuccd(pver)   ! mass tendency from ice nucleation
         606             : real(r8) :: qcld              ! total cloud water
         607             : real(r8) :: lcldn(pcols,pver) ! fractional coverage of new liquid cloud
         608             : real(r8) :: lcldo(pcols,pver) ! fractional coverage of old liquid cloud
         609             : real(r8) :: nctend_mixnuc(pcols,pver)
         610             : real(r8) :: arg ! argument of erfc
         611             :
         612             : ! for calculation of rate1ord_cw2pr_st
         613           0 : real(r8) :: qcsinksum_rate1ord(pver)   ! sum over iterations of cw to precip sink
         614           0 : real(r8) :: qcsum_rate1ord(pver)    ! sum over iterations of cloud water
         615             :
         616             : real(r8) :: alpha
         617             :
         618             : real(r8) :: dum1,dum2   !general dummy variables
         619             :
         620           0 : real(r8) :: npccn(pver)     ! droplet activation rate
         621           0 : real(r8) :: qcic(pcols,pver) ! in-cloud cloud liquid mixing ratio
         622           0 : real(r8) :: qiic(pcols,pver) ! in-cloud cloud ice mixing ratio
         623           0 : real(r8) :: qniic(pcols,pver) ! in-precip snow mixing ratio
         624           0 : real(r8) :: qric(pcols,pver) ! in-precip rain mixing ratio
         625           0 : real(r8) :: ncic(pcols,pver) ! in-cloud droplet number conc
         626           0 : real(r8) :: niic(pcols,pver) ! in-cloud cloud ice number conc
         627           0 : real(r8) :: nsic(pcols,pver) ! in-precip snow number conc
         628           0 : real(r8) :: nric(pcols,pver) ! in-precip rain number conc
         629           0 : real(r8) :: lami(pver) ! slope of cloud ice size distr
         630           0 : real(r8) :: n0i(pver) ! intercept of cloud ice size distr
         631           0 : real(r8) :: lamc(pver) ! slope of cloud liquid size distr
         632             : real(r8) :: n0c(pver) ! intercept of cloud liquid size distr
         633           0 : real(r8) :: lams(pver) ! slope of snow size distr
         634           0 : real(r8) :: n0s(pver) ! intercept of snow size distr
         635           0 : real(r8) :: lamr(pver) ! slope of rain size distr
         636           0 : real(r8) :: n0r(pver) ! intercept of rain size distr
         637           0 : real(r8) :: cdist1(pver) ! size distr parameter to calculate droplet freezing
         638             : ! combined size of precip & cloud drops
         639           0 : real(r8) :: arcld(pcols,pver) ! averaging control flag
         640             : real(r8) :: Actmp  !area cross section of drops
         641             : real(r8) :: Artmp  !area cross section of rain
         642             :
         643           0 : real(r8) :: pgam(pver) ! spectral width parameter of droplet size distr
         644             : real(r8) :: lammax  ! maximum allowed slope of size distr
         645             : real(r8) :: lammin  ! minimum allowed slope of size distr
         646             : real(r8) :: nacnt   ! number conc of contact ice nuclei
         647           0 : real(r8) :: mnuccc(pver) ! mixing ratio tendency due to freezing of cloud water
         648           0 : real(r8) :: nnuccc(pver) ! number conc tendency due to freezing of cloud water
         649             :
         650           0 : real(r8) :: mnucct(pver) ! mixing ratio tendency due to contact freezing of cloud water
         651           0 : real(r8) :: nnucct(pver) ! number conc tendency due to contact freezing of cloud water
         652           0 : real(r8) :: msacwi(pver) ! mixing ratio tendency due to HM ice multiplication
         653           0 : real(r8) :: nsacwi(pver) ! number conc tendency due to HM ice multiplication
         654             :
         655           0 : real(r8) :: prc(pver) ! qc tendency due to autoconversion of cloud droplets
         656           0 : real(r8) :: nprc(pver) ! number conc tendency due to autoconversion of cloud droplets
         657           0 : real(r8) :: nprc1(pver) ! qr tendency due to autoconversion of cloud droplets
         658           0 : real(r8) :: nsagg(pver) ! ns tendency due to self-aggregation of snow
         659             : real(r8) :: dc0  ! mean size droplet size distr
         660             : real(r8) :: ds0  ! mean size snow size distr (area weighted)
         661             : real(r8) :: eci  ! collection efficiency for riming of snow by droplets
         662           0 : real(r8) :: psacws(pver) ! mixing rat tendency due to collection of droplets by snow
         663           0 : real(r8) :: npsacws(pver) ! number conc tendency due to collection of droplets by snow
         664             : real(r8) :: uni ! number-weighted cloud ice fallspeed
         665             : real(r8) :: umi ! mass-weighted cloud ice fallspeed
         666           0 : real(r8) :: uns(pver) ! number-weighted snow fallspeed
         667           0 : real(r8) :: ums(pver) ! mass-weighted snow fallspeed
         668           0 : real(r8) :: unr(pver) ! number-weighted rain fallspeed
         669           0 : real(r8) :: umr(pver) ! mass-weighted rain fallspeed
         670             : real(r8) :: unc ! number-weighted cloud droplet fallspeed
         671             : real(r8) :: umc ! mass-weighted cloud droplet fallspeed
         672           0 : real(r8) :: pracs(pver) ! mixing rat tendency due to collection of rain by snow
         673           0 : real(r8) :: npracs(pver) ! number conc tendency due to collection of rain by snow
         674           0 : real(r8) :: mnuccr(pver) ! mixing rat tendency due to freezing of rain
         675           0 : real(r8) :: nnuccr(pver) ! number conc tendency due to freezing of rain
         676           0 : real(r8) :: pra(pver) ! mixing rat tendnency due to accretion of droplets by rain
         677           0 : real(r8) :: npra(pver) ! nc tendnency due to accretion of droplets by rain
         678           0 : real(r8) :: nragg(pver) ! nr tendency due to self-collection of rain
         679           0 : real(r8) :: prci(pver) ! mixing rat tendency due to autoconversion of cloud ice to snow
         680           0 : real(r8) :: nprci(pver) ! number conc tendency due to autoconversion of cloud ice to snow
         681           0 : real(r8) :: prai(pver) ! mixing rat tendency due to accretion of cloud ice by snow
         682           0 : real(r8) :: nprai(pver) ! number conc tendency due to accretion of cloud ice by snow
         683             : real(r8) :: qvs ! liquid saturation vapor mixing ratio
         684             : real(r8) :: qvi ! ice saturation vapor mixing ratio
         685             : real(r8) :: dqsdt ! change of sat vapor mixing ratio with temperature
         686             : real(r8) :: dqsidt ! change of ice sat vapor mixing ratio with temperature
         687             : real(r8) :: ab ! correction factor for rain evap to account for latent heat
         688             : real(r8) :: qclr ! water vapor mixing ratio in clear air
         689             : real(r8) :: abi ! correction factor for snow sublimation to account for latent heat
         690             : real(r8) :: epss ! 1/ sat relaxation timescale for snow
         691             : real(r8) :: epsr ! 1/ sat relaxation timescale for rain
         692           0 : real(r8) :: pre(pver) ! rain mixing rat tendency due to evaporation
         693           0 : real(r8) :: prds(pver) ! snow mixing rat tendency due to sublimation
         694             : real(r8) :: qce ! dummy qc for conservation check
         695             : real(r8) :: qie ! dummy qi for conservation check
         696             : real(r8) :: nce ! dummy nc for conservation check
         697             : real(r8) :: nie ! dummy ni for conservation check
         698             : real(r8) :: ratio ! parameter for conservation check
         699           0 : real(r8) :: dumc(pcols,pver) ! dummy in-cloud qc
         700           0 : real(r8) :: dumnc(pcols,pver) ! dummy in-cloud nc
         701           0 : real(r8) :: dumi(pcols,pver) ! dummy in-cloud qi
         702           0 : real(r8) :: dumni(pcols,pver) ! dummy in-cloud ni
         703             : real(r8) :: dums(pcols,pver) ! dummy in-cloud snow mixing rat
         704             : real(r8) :: dumns(pcols,pver) ! dummy in-cloud snow number conc
         705             : real(r8) :: dumr(pcols,pver) ! dummy in-cloud rain mixing rat
         706             : real(r8) :: dumnr(pcols,pver) ! dummy in-cloud rain number conc
         707             : ! below are parameters for cloud water and cloud ice sedimentation calculations
         708             : real(r8) :: fr(pver)
         709             : real(r8) :: fnr(pver)
         710           0 : real(r8) :: fc(pver)
         711           0 : real(r8) :: fnc(pver)
         712           0 : real(r8) :: fi(pver)
         713           0 : real(r8) :: fni(pver)
         714             : real(r8) :: fs(pver)
         715             : real(r8) :: fns(pver)
         716             : real(r8) :: faloutr(pver)
         717             : real(r8) :: faloutnr(pver)
         718           0 : real(r8) :: faloutc(pver)
         719           0 : real(r8) :: faloutnc(pver)
         720           0 : real(r8) :: falouti(pver)
         721           0 : real(r8) :: faloutni(pver)
         722             : real(r8) :: falouts(pver)
         723             : real(r8) :: faloutns(pver)
         724             : real(r8) :: faltndr
         725             : real(r8) :: faltndnr
         726             : real(r8) :: faltndc
         727             : real(r8) :: faltndnc
         728             : real(r8) :: faltndi
         729             : real(r8) :: faltndni
         730             : real(r8) :: faltnds
         731             : real(r8) :: faltndns
         732             : real(r8) :: faltndqie
         733             : real(r8) :: faltndqce
         734             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
         735           0 : real(r8) :: relhum(pcols,pver) ! relative humidity
         736             : real(r8) :: csigma(pcols) ! parameter for cond/evap of cloud water/ice
         737             : real(r8) :: rgvm ! max fallspeed for all species
         738           0 : real(r8) :: arn(pcols,pver) ! air density corrected rain fallspeed parameter
         739           0 : real(r8) :: asn(pcols,pver) ! air density corrected snow fallspeed parameter
         740           0 : real(r8) :: acn(pcols,pver) ! air density corrected cloud droplet fallspeed parameter
         741           0 : real(r8) :: ain(pcols,pver) ! air density corrected cloud ice fallspeed parameter
         742           0 : real(r8) :: nsubi(pver) ! evaporation of cloud ice number
         743           0 : real(r8) :: nsubc(pver) ! evaporation of droplet number
         744           0 : real(r8) :: nsubs(pver) ! evaporation of snow number
         745           0 : real(r8) :: nsubr(pver) ! evaporation of rain number
         746             : real(r8) :: mtime ! factor to account for droplet activation timescale
         747           0 : real(r8) :: dz(pcols,pver) ! height difference across model vertical level
         748             :
         749             :
         750             : !! add precip flux variables for sub-stepping
         751           0 : real(r8) :: rflx1(pcols,pver+1)
         752           0 : real(r8) :: sflx1(pcols,pver+1)
         753             :
         754             : ! returns from function/subroutine calls
         755             : real(r8) :: tsp(pcols,pver)      ! saturation temp (K)
         756             : real(r8) :: qsp(pcols,pver)      ! saturation mixing ratio (kg/kg)
         757             : real(r8) :: qsphy(pcols,pver)      ! saturation mixing ratio (kg/kg): hybrid rh
         758           0 : real(r8) :: qs(pcols)            ! liquid-ice weighted sat mixing rat (kg/kg)
         759           0 : real(r8) :: es(pcols)            ! liquid-ice weighted sat vapor press (pa)
         760           0 : real(r8) :: esl(pcols,pver)      ! liquid sat vapor pressure (pa)
         761           0 : real(r8) :: esi(pcols,pver)      ! ice sat vapor pressure (pa)
         762             :
         763             : ! sum of source/sink terms for diagnostic precip
         764             :
         765           0 : real(r8) :: qnitend(pcols,pver) ! snow mixing ratio source/sink term
         766           0 : real(r8) :: nstend(pcols,pver)  ! snow number concentration source/sink term
         767           0 : real(r8) :: qrtend(pcols,pver) ! rain mixing ratio source/sink term
         768           0 : real(r8) :: nrtend(pcols,pver)  ! rain number concentration source/sink term
         769             : real(r8) :: qrtot ! vertically-integrated rain mixing rat source/sink term
         770             : real(r8) :: nrtot ! vertically-integrated rain number conc source/sink term
         771             : real(r8) :: qstot ! vertically-integrated snow mixing rat source/sink term
         772             : real(r8) :: nstot ! vertically-integrated snow number conc source/sink term
         773             :
         774             : ! new terms for Bergeron process
         775             :
         776             : real(r8) :: dumnnuc ! provisional ice nucleation rate (for calculating bergeron)
         777             : real(r8) :: ninew  ! provisional cloud ice number conc (for calculating bergeron)
         778             : real(r8) :: qinew ! provisional cloud ice mixing ratio (for calculating bergeron)
         779             : real(r8) :: qvl  ! liquid sat mixing ratio
         780             : real(r8) :: epsi ! 1/ sat relaxation timecale for cloud ice
         781             : real(r8) :: prd ! provisional deposition rate of cloud ice at water sat
         782           0 : real(r8) :: berg(pcols,pver) ! mixing rat tendency due to bergeron process for cloud ice
         783           0 : real(r8) :: bergs(pver) ! mixing rat tendency due to bergeron process for snow
         784             :
         785             : !bergeron terms
         786             : real(r8) :: bergtsf   !bergeron timescale to remove all liquid
         787             : real(r8) :: rhin      !modified RH for vapor deposition
         788             :
         789             : ! diagnostic rain/snow for output to history
         790             : ! values are in-precip (local) !!!!
         791             :
         792           0 : real(r8) :: drout(pcols,pver)     ! rain diameter (m)
         793             :
         794             : !averageed rain/snow for history
         795             : real(r8) :: dumfice
         796             :
         797             : !ice nucleation, droplet activation
         798           0 : real(r8) :: dum2i(pcols,pver) ! number conc of ice nuclei available (1/kg)
         799           0 : real(r8) :: dum2l(pcols,pver) ! number conc of CCN (1/kg)
         800             : real(r8) :: ncmax
         801             : real(r8) :: nimax
         802             :
         803             : real(r8) :: qcvar     ! 1/relative variance of sub-grid qc
         804             :
         805             : ! loop array variables
         806             : integer i,k,nstep,n, l
         807             : integer ii,kk, m
         808             :
         809             : ! loop variables for sub-step solution
         810           0 : integer iter,it,ltrue(pcols)
         811             :
         812             : ! used in contact freezing via dust particles
         813             : real(r8)  tcnt, viscosity, mfp
         814             : real(r8)  slip1, slip2, slip3, slip4
         815             : !        real(r8)  dfaer1, dfaer2, dfaer3, dfaer4
         816             : !        real(r8)  nacon1,nacon2,nacon3,nacon4
         817             : real(r8)  ndfaer1, ndfaer2, ndfaer3, ndfaer4
         818             : real(r8)  nslip1, nslip2, nslip3, nslip4
         819             :
         820             : ! used in ice effective radius
         821             : real(r8)  bbi, cci, ak, iciwc, rvi
         822             :
         823             : ! used in Bergeron processe and water vapor deposition
         824             : real(r8)  Tk, deles, Aprpr, Bprpr, Cice, qi0, Crate, qidep
         825             :
         826             : ! mean cloud fraction over the time step
         827           0 : real(r8)  cldmw(pcols,pver)
         828             :
         829             : ! used in secondary ice production
         830             : real(r8) ni_secp
         831             :
         832             : ! variabels to check for RH after rain evap
         833             :
         834             : real(r8) :: esn
         835             : real(r8) :: qsn
         836             : real(r8) :: ttmp
         837             :
         838             :
         839             :
         840           0 : real(r8) :: rainrt(pcols,pver)  ! rain rate for reflectivity calculation
         841           0 : real(r8) :: rainrt1(pcols,pver)
         842             : real(r8) :: tmp
         843             :
         844             : real(r8) dmc,ssmc,dstrn  ! variables for modal scheme.
         845             :
         846             : real(r8), parameter :: cdnl    = 0.e6_r8    ! cloud droplet number limiter
         847             :
         848             : ! heterogeneous freezing
         849           0 : real(r8) :: mnudep(pver) ! mixing ratio tendency due to deposition of water vapor
         850           0 : real(r8) :: nnudep(pver) ! number conc tendency due to deposition of water vapor
         851             : real(r8) :: con1 ! work cnstant
         852             : real(r8) :: r3lx ! Mean volume radius (m)
         853             : real(r8) :: mi0l
         854             : real(r8) :: frztmp
         855             :
         856             : logical  :: do_clubb_sgs
         857             :
         858             : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         859             :
         860             : ! Return error message
         861           0 : errstring = ' '
         862             :
         863           0 : call phys_getopts(do_clubb_sgs_out = do_clubb_sgs)
         864             :
         865             : ! initialize  output fields for number conc qand ice nucleation
         866           0 : ncai(1:ncol,1:pver)=0._r8
         867           0 : ncal(1:ncol,1:pver)=0._r8
         868             :
         869             : !Initialize rain size
         870           0 : rercld(1:ncol,1:pver)=0._r8
         871           0 : arcld(1:ncol,1:pver)=0._r8
         872             :
         873             : !initialize radiation output variables
         874           0 : pgamrad(1:ncol,1:pver)=0._r8 ! liquid gamma parameter for optics (radiation)
         875           0 : lamcrad(1:ncol,1:pver)=0._r8 ! slope of droplet distribution for optics (radiation)
         876           0 : deffi  (1:ncol,1:pver)=0._r8 ! slope of droplet distribution for optics (radiation)
         877             : !initialize radiation output variables
         878             : !initialize water vapor tendency term output
         879           0 : qcsevap(1:ncol,1:pver)=0._r8
         880           0 : qisevap(1:ncol,1:pver)=0._r8
         881           0 : qvres  (1:ncol,1:pver)=0._r8
         882           0 : cmeiout (1:ncol,1:pver)=0._r8
         883           0 : vtrmc (1:ncol,1:pver)=0._r8
         884           0 : vtrmi (1:ncol,1:pver)=0._r8
         885           0 : qcsedten (1:ncol,1:pver)=0._r8
         886           0 : qisedten (1:ncol,1:pver)=0._r8
         887             :
         888           0 : prao(1:ncol,1:pver)=0._r8
         889           0 : prco(1:ncol,1:pver)=0._r8
         890           0 : mnuccco(1:ncol,1:pver)=0._r8
         891           0 : mnuccto(1:ncol,1:pver)=0._r8
         892           0 : msacwio(1:ncol,1:pver)=0._r8
         893           0 : psacwso(1:ncol,1:pver)=0._r8
         894           0 : bergso(1:ncol,1:pver)=0._r8
         895           0 : bergo(1:ncol,1:pver)=0._r8
         896           0 : melto(1:ncol,1:pver)=0._r8
         897           0 : homoo(1:ncol,1:pver)=0._r8
         898           0 : qcreso(1:ncol,1:pver)=0._r8
         899           0 : prcio(1:ncol,1:pver)=0._r8
         900           0 : praio(1:ncol,1:pver)=0._r8
         901           0 : qireso(1:ncol,1:pver)=0._r8
         902           0 : mnuccro(1:ncol,1:pver)=0._r8
         903           0 : pracso (1:ncol,1:pver)=0._r8
         904           0 : meltsdt(1:ncol,1:pver)=0._r8
         905           0 : frzrdt (1:ncol,1:pver)=0._r8
         906           0 : mnuccdo(1:ncol,1:pver)=0._r8
         907             :
         908           0 : rflx(:,:)=0._r8
         909           0 : sflx(:,:)=0._r8
         910           0 : effc(:,:)=0._r8
         911           0 : effc_fn(:,:)=0._r8
         912           0 : effi(:,:)=0._r8
         913             :
         914             : ! assign variable deltat for sub-stepping...
         915           0 : deltat=deltatin
         916             :
         917             : ! parameters for scheme
         918             :
         919           0 : omsm=0.99999_r8
         920           0 : dto2=0.5_r8*deltat
         921           0 : mincld=0.0001_r8
         922             :
         923             : ! initialize multi-level fields
         924           0 : q(1:ncol,1:pver)=qn(1:ncol,1:pver)
         925           0 : t(1:ncol,1:pver)=tn(1:ncol,1:pver)
         926             :
         927             : ! initialize time-varying parameters
         928             :
         929           0 : do k=1,pver
         930           0 :    do i=1,ncol
         931           0 :       rho(i,k)=p(i,k)/(r*t(i,k))
         932           0 :       dv(i,k) = 8.794E-5_r8*t(i,k)**1.81_r8/p(i,k)
         933           0 :       mu(i,k) = 1.496E-6_r8*t(i,k)**1.5_r8/(t(i,k)+120._r8)
         934           0 :       sc(i,k) = mu(i,k)/(rho(i,k)*dv(i,k))
         935           0 :       kap(i,k) = 1.414e3_r8*1.496e-6_r8*t(i,k)**1.5_r8/(t(i,k)+120._r8)
         936             :
         937             :       ! air density adjustment for fallspeed parameters
         938             :       ! includes air density correction factor to the
         939             :       ! power of 0.54 following Heymsfield and Bansemer 2007
         940             :
         941           0 :       rhof(i,k)=(rhosu/rho(i,k))**0.54_r8
         942             :
         943           0 :       arn(i,k)=ar*rhof(i,k)
         944           0 :       asn(i,k)=as*rhof(i,k)
         945           0 :       acn(i,k)=ac*rhof(i,k)
         946           0 :       ain(i,k)=ai*rhof(i,k)
         947             :
         948             :       ! get dz from dp and hydrostatic approx
         949             :       ! keep dz positive (define as layer k-1 - layer k)
         950             :
         951           0 :       dz(i,k)= pdel(i,k)/(rho(i,k)*g)
         952             :
         953             :    end do
         954             : end do
         955             :
         956             : ! initialization
         957           0 : qc(1:ncol,1:top_lev-1) = 0._r8
         958           0 : qi(1:ncol,1:top_lev-1) = 0._r8
         959           0 : nc(1:ncol,1:top_lev-1) = 0._r8
         960           0 : ni(1:ncol,1:top_lev-1) = 0._r8
         961           0 : t1(1:ncol,1:pver) = t(1:ncol,1:pver)
         962           0 : q1(1:ncol,1:pver) = q(1:ncol,1:pver)
         963           0 : qc1(1:ncol,1:pver) = qc(1:ncol,1:pver)
         964           0 : qi1(1:ncol,1:pver) = qi(1:ncol,1:pver)
         965           0 : nc1(1:ncol,1:pver) = nc(1:ncol,1:pver)
         966           0 : ni1(1:ncol,1:pver) = ni(1:ncol,1:pver)
         967             :
         968             : ! initialize tendencies to zero
         969           0 : tlat1(1:ncol,1:pver)=0._r8
         970           0 : qvlat1(1:ncol,1:pver)=0._r8
         971           0 : qctend1(1:ncol,1:pver)=0._r8
         972           0 : qitend1(1:ncol,1:pver)=0._r8
         973           0 : nctend1(1:ncol,1:pver)=0._r8
         974           0 : nitend1(1:ncol,1:pver)=0._r8
         975             :
         976             : ! initialize precip output
         977           0 : qrout(1:ncol,1:pver)=0._r8
         978           0 : qsout(1:ncol,1:pver)=0._r8
         979           0 : nrout(1:ncol,1:pver)=0._r8
         980           0 : nsout(1:ncol,1:pver)=0._r8
         981           0 : dsout(1:ncol,1:pver)=0._r8
         982             :
         983           0 : drout(1:ncol,1:pver)=0._r8
         984             :
         985           0 : reff_rain(1:ncol,1:pver)=0._r8
         986           0 : reff_snow(1:ncol,1:pver)=0._r8
         987             :
         988             : ! initialize variables for trop_mozart
         989           0 : nevapr(1:ncol,1:pver) = 0._r8
         990           0 : nevapr2(1:ncol,1:pver) = 0._r8
         991           0 : evapsnow(1:ncol,1:pver) = 0._r8
         992           0 : prain(1:ncol,1:pver) = 0._r8
         993           0 : prodsnow(1:ncol,1:pver) = 0._r8
         994           0 : cmeout(1:ncol,1:pver) = 0._r8
         995             :
         996           0 : am_evp_st(1:ncol,1:pver) = 0._r8
         997             :
         998             : ! for refl calc
         999           0 : rainrt1(1:ncol,1:pver) = 0._r8
        1000             :
        1001             : ! initialize precip fraction and output tendencies
        1002           0 : cldmax(1:ncol,1:pver)=mincld
        1003             :
        1004             : !initialize aerosol number
        1005             : !        naer2(1:ncol,1:pver,:)=0._r8
        1006           0 : dum2l(1:ncol,1:pver)=0._r8
        1007           0 : dum2i(1:ncol,1:pver)=0._r8
        1008             :
        1009             : ! initialize avg precip rate
        1010           0 : prect1(1:ncol)=0._r8
        1011           0 : preci1(1:ncol)=0._r8
        1012             :
        1013             : !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1014             : !Get humidity and saturation vapor pressures
        1015             :
        1016           0 : do k=top_lev,pver
        1017             :
        1018           0 :    do i=1,ncol
        1019             :
        1020             :       ! find wet bulk temperature and saturation value for provisional t and q without
        1021             :       ! condensation
        1022             :
        1023           0 :       es(i) = svp_water(t(i,k))
        1024           0 :       qs(i) = svp_to_qsat(es(i), p(i,k))
        1025             :
        1026             :       ! Prevents negative values.
        1027           0 :       if (qs(i) < 0.0_r8) then
        1028           0 :          qs(i) = 1.0_r8
        1029           0 :          es(i) = p(i,k)
        1030             :       end if
        1031             :
        1032           0 :       esl(i,k)=svp_water(t(i,k))
        1033           0 :       esi(i,k)=svp_ice(t(i,k))
        1034             :
        1035             :       ! hm fix, make sure when above freezing that esi=esl, not active yet
        1036           0 :       if (t(i,k).gt.tmelt)esi(i,k)=esl(i,k)
        1037             :
        1038           0 :       relhum(i,k)=q(i,k)/qs(i)
        1039             :
        1040             :       ! get cloud fraction, check for minimum
        1041             :
        1042           0 :       cldm(i,k)=max(cldn(i,k),mincld)
        1043           0 :       cldmw(i,k)=max(cldn(i,k),mincld)
        1044             :
        1045           0 :       icldm(i,k)=max(icecldf(i,k),mincld)
        1046           0 :       lcldm(i,k)=max(liqcldf(i,k),mincld)
        1047             :
        1048             :       ! subcolumns, set cloud fraction variables to one
        1049             :       ! if cloud water or ice is present, if not present
        1050             :       ! set to mincld (mincld used instead of zero, to prevent
        1051             :       ! possible division by zero errors
        1052             :
        1053           0 :       if (microp_uniform) then
        1054             :
        1055           0 :          cldm(i,k)=mincld
        1056           0 :          cldmw(i,k)=mincld
        1057           0 :          icldm(i,k)=mincld
        1058           0 :          lcldm(i,k)=mincld
        1059             :
        1060           0 :          if (qc(i,k).ge.qsmall) then
        1061           0 :             lcldm(i,k)=1._r8
        1062           0 :             cldm(i,k)=1._r8
        1063           0 :             cldmw(i,k)=1._r8
        1064             :          end if
        1065             :
        1066           0 :          if (qi(i,k).ge.qsmall) then
        1067           0 :             cldm(i,k)=1._r8
        1068           0 :             icldm(i,k)=1._r8
        1069             :          end if
        1070             :
        1071             :       end if               ! sub-columns
        1072             :
        1073             :       ! calculate nfice based on liquid and ice mmr (no rain and snow mmr available yet)
        1074             :
        1075           0 :       nfice(i,k)=0._r8
        1076           0 :       dumfice=qc(i,k)+qi(i,k)
        1077           0 :       if (dumfice.gt.qsmall .and. qi(i,k).gt.qsmall) then
        1078           0 :          nfice(i,k)=qi(i,k)/dumfice
        1079             :       endif
        1080             :
        1081           0 :       if (do_cldice .and. (t(i,k).lt.tmelt - 5._r8)) then
        1082             :
        1083             :          ! if aerosols interact with ice set number of activated ice nuclei
        1084           0 :          dum2=naai(i,k)
        1085             :
        1086           0 :          dumnnuc=(dum2-ni(i,k)/icldm(i,k))/deltat*icldm(i,k)
        1087           0 :          dumnnuc=max(dumnnuc,0._r8)
        1088             :          ! get provisional ni and qi after nucleation in order to calculate
        1089             :          ! Bergeron process below
        1090           0 :          ninew=ni(i,k)+dumnnuc*deltat
        1091           0 :          qinew=qi(i,k)+dumnnuc*deltat*mi0
        1092             :
        1093             :          !T>268
        1094             :       else
        1095           0 :          ninew=ni(i,k)
        1096           0 :          qinew=qi(i,k)
        1097             :       end if
        1098             :
        1099             :       ! Initialize CME components
        1100             :
        1101           0 :       cme(i,k) = 0._r8
        1102           0 :       cmei(i,k)=0._r8
        1103             :
        1104             :
        1105             :       !-------------------------------------------------------------------
        1106             :       !Bergeron process
        1107             :
        1108             :       ! make sure to initialize bergeron process to zero
        1109           0 :       berg(i,k)=0._r8
        1110           0 :       prd = 0._r8
        1111             :
        1112             :       !condensation loop.
        1113             :
        1114             :       ! get in-cloud qi and ni after nucleation
        1115           0 :       if (icldm(i,k) .gt. 0._r8) then
        1116           0 :          qiic(i,k)=qinew/icldm(i,k)
        1117           0 :          niic(i,k)=ninew/icldm(i,k)
        1118             :       else
        1119           0 :          qiic(i,k)=0._r8
        1120           0 :          niic(i,k)=0._r8
        1121             :       endif
        1122             :
        1123           0 :       if (nicons) then
        1124           0 :         niic(i,k) = ninst/rho(i,k)
        1125             :       end if
        1126             :
        1127             :       !if T < 0 C then bergeron.
        1128           0 :       if (do_cldice .and. (t(i,k).lt.273.15_r8)) then
        1129             :
        1130             :          !if ice exists
        1131           0 :          if (qi(i,k).gt.qsmall) then
        1132             :
        1133           0 :             bergtsf = 0._r8 ! bergeron time scale (fraction of timestep)
        1134             :
        1135           0 :             qvi = svp_to_qsat(esi(i,k), p(i,k))
        1136           0 :             qvl = svp_to_qsat(esl(i,k), p(i,k))
        1137             :
        1138           0 :             dqsidt =  xxls*qvi/(rv*t(i,k)**2)
        1139           0 :             abi = 1._r8+dqsidt*xxls/cpp
        1140             :
        1141             :             ! get ice size distribution parameters
        1142             :
        1143           0 :             if (qiic(i,k).ge.qsmall) then
        1144             :                lami(k) = (cons1*ci* &
        1145           0 :                     niic(i,k)/qiic(i,k))**(1._r8/di)
        1146           0 :                n0i(k) = niic(i,k)*lami(k)
        1147             :
        1148             :                ! check for slope
        1149             :                ! adjust vars
        1150           0 :                if (lami(k).lt.lammini) then
        1151             :
        1152           0 :                   lami(k) = lammini
        1153           0 :                   n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
        1154           0 :                else if (lami(k).gt.lammaxi) then
        1155           0 :                   lami(k) = lammaxi
        1156           0 :                   n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
        1157             :                end if
        1158             :
        1159           0 :                epsi = 2._r8*pi*n0i(k)*rho(i,k)*Dv(i,k)/(lami(k)*lami(k))
        1160             :
        1161             :                !if liquid exists
        1162           0 :                if (qc(i,k).gt. qsmall) then
        1163             :
        1164             :                   !begin bergeron process
        1165             :                   !     do bergeron (vapor deposition with RHw=1)
        1166             :                   !     code to find berg (a rate) goes here
        1167             :
        1168             :                   ! calculate Bergeron process
        1169             :
        1170           0 :                   prd = epsi*(qvl-qvi)/abi
        1171             :
        1172             :                else
        1173             :                   prd = 0._r8
        1174             :                end if
        1175             :
        1176             :                ! multiply by cloud fraction
        1177             :
        1178           0 :                prd = prd*min(icldm(i,k),lcldm(i,k))
        1179             :
        1180             :                !     transfer of existing cloud liquid to ice
        1181             :
        1182           0 :                berg(i,k)=max(0._r8,prd)
        1183             :
        1184             :             end if  !end liquid exists bergeron
        1185             :
        1186           0 :             if (berg(i,k).gt.0._r8) then
        1187           0 :                bergtsf=max(0._r8,(qc(i,k)/berg(i,k))/deltat)
        1188             :
        1189           0 :                if(bergtsf.lt.1._r8) berg(i,k) = max(0._r8,qc(i,k)/deltat)
        1190             :
        1191             :             endif
        1192             :
        1193             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        1194             :
        1195           0 :             if (bergtsf.lt.1._r8.or.icldm(i,k).gt.lcldm(i,k)) then
        1196             :
        1197           0 :                if (qiic(i,k).ge.qsmall) then
        1198             :
        1199             :                   ! first case is for case when liquid water is present, but is completely depleted
        1200             :                   ! in time step, i.e., bergrsf > 0 but < 1
        1201             :
        1202           0 :                   if (qc(i,k).ge.qsmall) then
        1203           0 :                      rhin  = (1.0_r8 + relhum(i,k)) / 2._r8
        1204           0 :                      if ((rhin*esl(i,k)/esi(i,k)) > 1._r8) then
        1205           0 :                         prd = epsi*(rhin*qvl-qvi)/abi
        1206             :
        1207             :                         ! multiply by cloud fraction assuming liquid/ice maximum overlap
        1208           0 :                         prd = prd*min(icldm(i,k),lcldm(i,k))
        1209             :
        1210             :                         ! add to cmei
        1211           0 :                         cmei(i,k) = cmei(i,k) + (prd * (1._r8- bergtsf))
        1212             :
        1213             :                      end if ! rhin
        1214             :                   end if ! qc > qsmall
        1215             :
        1216             :                   ! second case is for pure ice cloud, either no liquid, or icldm > lcldm
        1217             :
        1218           0 :                   if (qc(i,k).lt.qsmall.or.icldm(i,k).gt.lcldm(i,k)) then
        1219             :
        1220             :                      ! note: for case of no liquid, need to set liquid cloud fraction to zero
        1221             :                      ! store liquid cloud fraction in 'dum'
        1222             :
        1223           0 :                      if (qc(i,k).lt.qsmall) then
        1224             :                         dum=0._r8
        1225             :                      else
        1226           0 :                         dum=lcldm(i,k)
        1227             :                      end if
        1228             :
        1229             :                      ! set RH to grid-mean value for pure ice cloud
        1230           0 :                      rhin = relhum(i,k)
        1231             :
        1232           0 :                      if ((rhin*esl(i,k)/esi(i,k)) > 1._r8) then
        1233             :
        1234           0 :                         prd = epsi*(rhin*qvl-qvi)/abi
        1235             :
        1236             :                         ! multiply by relevant cloud fraction for pure ice cloud
        1237             :                         ! assuming maximum overlap of liquid/ice
        1238           0 :                         prd = prd*max((icldm(i,k)-dum),0._r8)
        1239           0 :                         cmei(i,k) = cmei(i,k) + prd
        1240             :
        1241             :                      end if ! rhin
        1242             :                   end if ! qc or icldm > lcldm
        1243             :                end if ! qiic
        1244             :             end if ! bergtsf or icldm > lcldm
        1245             :
        1246             :             !     if deposition, it should not reduce grid mean rhi below 1.0
        1247           0 :             if(cmei(i,k) > 0.0_r8 .and. (relhum(i,k)*esl(i,k)/esi(i,k)) > 1._r8 ) &
        1248           0 :                  cmei(i,k)=min(cmei(i,k),(q(i,k)-qs(i)*esi(i,k)/esl(i,k))/abi/deltat)
        1249             :
        1250             :          end if            !end ice exists loop
        1251             :          !this ends temperature < 0. loop
        1252             :
        1253             :          !-------------------------------------------------------------------
        1254             :       end if  !
        1255             :       !..............................................................
        1256             :
        1257             :       ! evaporation should not exceed available water
        1258             :
        1259           0 :       if ((-berg(i,k)).lt.-qc(i,k)/deltat) berg(i,k) = max(qc(i,k)/deltat,0._r8)
        1260             :
        1261             :       !sublimation process...
        1262           0 :       if (do_cldice .and. ((relhum(i,k)*esl(i,k)/esi(i,k)).lt.1._r8 .and. qiic(i,k).ge.qsmall )) then
        1263             :
        1264           0 :          qvi = svp_to_qsat(esi(i,k), p(i,k))
        1265           0 :          qvl = svp_to_qsat(esl(i,k), p(i,k))
        1266           0 :          dqsidt =  xxls*qvi/(rv*t(i,k)**2)
        1267           0 :          abi = 1._r8+dqsidt*xxls/cpp
        1268             :
        1269             :          ! get ice size distribution parameters
        1270             :
        1271             :          lami(k) = (cons1*ci* &
        1272           0 :               niic(i,k)/qiic(i,k))**(1._r8/di)
        1273           0 :          n0i(k) = niic(i,k)*lami(k)
        1274             :
        1275             :          ! check for slope
        1276             :          ! adjust vars
        1277           0 :          if (lami(k).lt.lammini) then
        1278             :
        1279           0 :             lami(k) = lammini
        1280           0 :             n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
        1281           0 :          else if (lami(k).gt.lammaxi) then
        1282           0 :             lami(k) = lammaxi
        1283           0 :             n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
        1284             :          end if
        1285             :
        1286           0 :          epsi = 2._r8*pi*n0i(k)*rho(i,k)*Dv(i,k)/(lami(k)*lami(k))
        1287             :
        1288             :          ! modify for ice fraction below
        1289           0 :          prd = epsi*(relhum(i,k)*qvl-qvi)/abi * icldm(i,k)
        1290           0 :          cmei(i,k)=min(prd,0._r8)
        1291             :
        1292             :       endif
        1293             :
        1294             :       ! sublimation should not exceed available ice
        1295           0 :       if (cmei(i,k).lt.-qi(i,k)/deltat) cmei(i,k)=-qi(i,k)/deltat
        1296             :
        1297             :       ! sublimation should not increase grid mean rhi above 1.0
        1298           0 :       if(cmei(i,k) < 0.0_r8 .and. (relhum(i,k)*esl(i,k)/esi(i,k)) < 1._r8 ) &
        1299           0 :            cmei(i,k)=min(0._r8,max(cmei(i,k),(q(i,k)-qs(i)*esi(i,k)/esl(i,k))/abi/deltat))
        1300             :
        1301             :       ! limit cmei due for roundoff error
        1302             :
        1303           0 :       cmei(i,k)=cmei(i,k)*omsm
        1304             :
        1305             :       ! conditional for ice nucleation
        1306           0 :       if (do_cldice .and. (t(i,k).lt.(tmelt - 5._r8))) then
        1307             :
        1308             :          ! using Liu et al. (2007) ice nucleation with hooks into simulated aerosol
        1309             :          ! ice nucleation rate (dum2) has already been calculated and read in (naai)
        1310             :
        1311           0 :          dum2i(i,k)=naai(i,k)
        1312             :       else
        1313           0 :          dum2i(i,k)=0._r8
        1314             :       end if
        1315             :
        1316             :    end do ! i loop
        1317             : end do ! k loop
        1318             :
        1319             :
        1320             : !! initialize sub-step precip flux variables
        1321           0 : do i=1,ncol
        1322             :    !! flux is zero at top interface, so these should stay as 0.
        1323           0 :    rflx1(i,1)=0._r8
        1324           0 :    sflx1(i,1)=0._r8
        1325           0 :    do k=top_lev,pver
        1326             :
        1327             :       ! initialize normal and sub-step precip flux variables
        1328           0 :       rflx1(i,k+1)=0._r8
        1329           0 :       sflx1(i,k+1)=0._r8
        1330             :    end do ! i loop
        1331             : end do ! k loop
        1332             : !! initialize final precip flux variables.
        1333           0 : do i=1,ncol
        1334             :    !! flux is zero at top interface, so these should stay as 0.
        1335           0 :    rflx(i,1)=0._r8
        1336           0 :    sflx(i,1)=0._r8
        1337           0 :    do k=top_lev,pver
        1338             :       ! initialize normal and sub-step precip flux variables
        1339           0 :       rflx(i,k+1)=0._r8
        1340           0 :       sflx(i,k+1)=0._r8
        1341             :    end do ! i loop
        1342             : end do ! k loop
        1343             :
        1344           0 : do i=1,ncol
        1345           0 :    ltrue(i)=0
        1346           0 :    do k=top_lev,pver
        1347             :       ! skip microphysical calculations if no cloud water
        1348             :
        1349           0 :       if (qc(i,k).ge.qsmall.or.qi(i,k).ge.qsmall.or.cmei(i,k).ge.qsmall) ltrue(i)=1
        1350             :    end do
        1351             : end do
        1352             :
        1353             : ! assign number of sub-steps to iter
        1354             : ! use 2 sub-steps, following tests described in MG2008
        1355             : iter = 2
        1356             :
        1357             : ! get sub-step time step
        1358           0 : deltat=deltat/real(iter)
        1359             :
        1360             : ! since activation/nucleation processes are fast, need to take into account
        1361             : ! factor mtime = mixing timescale in cloud / model time step
        1362             : ! mixing time can be interpreted as cloud depth divided by sub-grid vertical velocity
        1363             : ! for now mixing timescale is assumed to be 1 timestep for modal aerosols, 20 min bulk
        1364             :
        1365             : !        note: mtime for bulk aerosols was set to: mtime=deltat/1200._r8
        1366             :
        1367             : mtime=1._r8
        1368           0 : rate1ord_cw2pr_st(:,:)=0._r8 ! rce 2010/05/01
        1369             :
        1370             : !!!! skip calculations if no cloud water
        1371           0 : do i=1,ncol
        1372           0 :    if (ltrue(i).eq.0) then
        1373           0 :       tlat(i,1:pver)=0._r8
        1374           0 :       qvlat(i,1:pver)=0._r8
        1375           0 :       qctend(i,1:pver)=0._r8
        1376           0 :       qitend(i,1:pver)=0._r8
        1377           0 :       qnitend(i,1:pver)=0._r8
        1378           0 :       qrtend(i,1:pver)=0._r8
        1379           0 :       nctend(i,1:pver)=0._r8
        1380           0 :       nitend(i,1:pver)=0._r8
        1381           0 :       nrtend(i,1:pver)=0._r8
        1382           0 :       nstend(i,1:pver)=0._r8
        1383           0 :       prect(i)=0._r8
        1384           0 :       preci(i)=0._r8
        1385           0 :       rflx(i,1:pver+1)=0._r8
        1386           0 :       sflx(i,1:pver+1)=0._r8
        1387           0 :       qniic(i,1:pver)=0._r8
        1388           0 :       qric(i,1:pver)=0._r8
        1389           0 :       nsic(i,1:pver)=0._r8
        1390           0 :       nric(i,1:pver)=0._r8
        1391           0 :       rainrt(i,1:pver)=0._r8
        1392             :       goto 300
        1393             :    end if
        1394             :
        1395           0 :    qcsinksum_rate1ord(1:pver)=0._r8
        1396           0 :    qcsum_rate1ord(1:pver)=0._r8
        1397             :
        1398             :
        1399             : !!!!!!!!! begin sub-step!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        1400             :    !.....................................................................................................
        1401           0 :    do it=1,iter
        1402             :
        1403             :       ! initialize sub-step microphysical tendencies
        1404             :
        1405           0 :       tlat(i,1:pver)=0._r8
        1406           0 :       qvlat(i,1:pver)=0._r8
        1407           0 :       qctend(i,1:pver)=0._r8
        1408           0 :       qitend(i,1:pver)=0._r8
        1409           0 :       qnitend(i,1:pver)=0._r8
        1410           0 :       qrtend(i,1:pver)=0._r8
        1411           0 :       nctend(i,1:pver)=0._r8
        1412           0 :       nitend(i,1:pver)=0._r8
        1413           0 :       nrtend(i,1:pver)=0._r8
        1414           0 :       nstend(i,1:pver)=0._r8
        1415             :
        1416             :       ! initialize diagnostic precipitation to zero
        1417             :
        1418           0 :       qniic(i,1:pver)=0._r8
        1419           0 :       qric(i,1:pver)=0._r8
        1420           0 :       nsic(i,1:pver)=0._r8
        1421           0 :       nric(i,1:pver)=0._r8
        1422             :
        1423           0 :       rainrt(i,1:pver)=0._r8
        1424             :
        1425             :
        1426             :       ! begin new i,k loop, calculate new cldmax after adjustment to cldm above
        1427             :
        1428             :       ! initialize vertically-integrated rain and snow tendencies
        1429             :
        1430           0 :       qrtot = 0._r8
        1431           0 :       nrtot = 0._r8
        1432           0 :       qstot = 0._r8
        1433           0 :       nstot = 0._r8
        1434             :
        1435             :       ! initialize precip at surface
        1436             :
        1437           0 :       prect(i)=0._r8
        1438           0 :       preci(i)=0._r8
        1439             :
        1440             :       ! initialize fluxes
        1441           0 :       rflx(i,1:pver+1)=0._r8
        1442           0 :       sflx(i,1:pver+1)=0._r8
        1443             :
        1444           0 :       do k=top_lev,pver
        1445             :
        1446           0 :          qcvar=relvar(i,k)
        1447           0 :          cons2=gamma(qcvar+2.47_r8)
        1448           0 :          cons3=gamma(qcvar)
        1449           0 :          cons9=gamma(qcvar+2._r8)
        1450           0 :          cons10=gamma(qcvar+1._r8)
        1451           0 :          cons12=gamma(qcvar+1.15_r8)
        1452           0 :          cons15=gamma(qcvar+bc/3._r8)
        1453           0 :          cons18=qcvar**2.47_r8
        1454           0 :          cons19=qcvar**2
        1455           0 :          cons20=qcvar**1.15_r8
        1456             :
        1457             :          ! set cwml and cwmi to current qc and qi
        1458             :
        1459           0 :          cwml(i,k)=qc(i,k)
        1460           0 :          cwmi(i,k)=qi(i,k)
        1461             :
        1462             :          ! initialize precip fallspeeds to zero
        1463             :
        1464           0 :          ums(k)=0._r8
        1465           0 :          uns(k)=0._r8
        1466           0 :          umr(k)=0._r8
        1467           0 :          unr(k)=0._r8
        1468             :
        1469             :          ! calculate precip fraction based on maximum overlap assumption
        1470             :
        1471             :          ! for sub-columns cldm has already been set to 1 if cloud
        1472             :          ! water or ice is present, so cldmax will be correctly set below
        1473             :          ! and nothing extra needs to be done here
        1474             :
        1475           0 :          if (k.eq.top_lev) then
        1476           0 :             cldmax(i,k)=cldm(i,k)
        1477             :          else
        1478             :             ! if rain or snow mix ratio is smaller than
        1479             :             ! threshold, then set cldmax to cloud fraction at current level
        1480             :
        1481           0 :             if (do_clubb_sgs) then
        1482           0 :                if (qc(i,k).ge.qsmall.or.qi(i,k).ge.qsmall) then
        1483           0 :                   cldmax(i,k)=cldm(i,k)
        1484             :                else
        1485           0 :                   cldmax(i,k)=cldmax(i,k-1)
        1486             :                end if
        1487             :             else
        1488             :
        1489           0 :                if (qric(i,k-1).ge.qsmall.or.qniic(i,k-1).ge.qsmall) then
        1490           0 :                   cldmax(i,k)=max(cldmax(i,k-1),cldm(i,k))
        1491             :                else
        1492           0 :                   cldmax(i,k)=cldm(i,k)
        1493             :                end if
        1494             :             endif
        1495             :          end if
        1496             :
        1497             :          ! decrease in number concentration due to sublimation/evap
        1498             :          ! divide by cloud fraction to get in-cloud decrease
        1499             :          ! don't reduce Nc due to bergeron process
        1500             :
        1501           0 :          if (cmei(i,k) < 0._r8 .and. qi(i,k) > qsmall .and. cldm(i,k) > mincld) then
        1502           0 :             nsubi(k)=cmei(i,k)/qi(i,k)*ni(i,k)/cldm(i,k)
        1503             :          else
        1504           0 :             nsubi(k)=0._r8
        1505             :          end if
        1506           0 :          nsubc(k)=0._r8
        1507             :
        1508             :
        1509             :          !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1510             :          !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1511             :
        1512             :          ! ice nucleation if activated nuclei exist at t<-5C AND rhmini + 5%
        1513             :
        1514           0 :          if (do_cldice .and. dum2i(i,k).gt.0._r8.and.t(i,k).lt.(tmelt - 5._r8).and. &
        1515             :               relhum(i,k)*esl(i,k)/esi(i,k).gt. rhmini+0.05_r8) then
        1516             :
        1517             :             !if NCAI > 0. then set numice = ncai (as before)
        1518             :             !note: this is gridbox averaged
        1519             :
        1520           0 :             nnuccd(k)=(dum2i(i,k)-ni(i,k)/icldm(i,k))/deltat*icldm(i,k)
        1521           0 :             nnuccd(k)=max(nnuccd(k),0._r8)
        1522           0 :             nimax = dum2i(i,k)*icldm(i,k)
        1523             :
        1524             :             !Calc mass of new particles using new crystal mass...
        1525             :             !also this will be multiplied by mtime as nnuccd is...
        1526             :
        1527           0 :             mnuccd(k) = nnuccd(k) * mi0
        1528             :
        1529             :             !  add mnuccd to cmei....
        1530           0 :             cmei(i,k)= cmei(i,k) + mnuccd(k) * mtime
        1531             :
        1532             :             !  limit cmei
        1533             :
        1534           0 :             qvi = svp_to_qsat(esi(i,k), p(i,k))
        1535           0 :             dqsidt =  xxls*qvi/(rv*t(i,k)**2)
        1536           0 :             abi = 1._r8+dqsidt*xxls/cpp
        1537           0 :             cmei(i,k)=min(cmei(i,k),(q(i,k)-qvi)/abi/deltat)
        1538             :
        1539             :             ! limit for roundoff error
        1540           0 :             cmei(i,k)=cmei(i,k)*omsm
        1541             :
        1542             :          else
        1543           0 :             nnuccd(k)=0._r8
        1544           0 :             nimax = 0._r8
        1545           0 :             mnuccd(k) = 0._r8
        1546             :          end if
        1547             :
        1548             :          !c............................................................................
        1549             :          !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1550             :          ! obtain in-cloud values of cloud water/ice mixing ratios and number concentrations
        1551             :          ! for microphysical process calculations
        1552             :          ! units are kg/kg for mixing ratio, 1/kg for number conc
        1553             :
        1554             :          ! limit in-cloud values to 0.005 kg/kg
        1555             :
        1556           0 :          qcic(i,k)=min(cwml(i,k)/lcldm(i,k),5.e-3_r8)
        1557           0 :          qiic(i,k)=min(cwmi(i,k)/icldm(i,k),5.e-3_r8)
        1558           0 :          ncic(i,k)=max(nc(i,k)/lcldm(i,k),0._r8)
        1559           0 :          niic(i,k)=max(ni(i,k)/icldm(i,k),0._r8)
        1560             :
        1561           0 :          if (nccons) then
        1562           0 :            ncic(i,k) = ncnst/rho(i,k)
        1563             :          end if
        1564           0 :          if (nicons) then
        1565           0 :            niic(i,k) = ninst/rho(i,k)
        1566             :          end if
        1567             :
        1568           0 :          if (qc(i,k) - berg(i,k)*deltat.lt.qsmall) then
        1569           0 :             qcic(i,k)=0._r8
        1570           0 :             ncic(i,k)=0._r8
        1571           0 :             if (qc(i,k)-berg(i,k)*deltat.lt.0._r8) then
        1572           0 :                berg(i,k)=qc(i,k)/deltat*omsm
        1573             :             end if
        1574             :          end if
        1575             :
        1576           0 :          if (do_cldice .and. qi(i,k)+(cmei(i,k)+berg(i,k))*deltat.lt.qsmall) then
        1577           0 :             qiic(i,k)=0._r8
        1578           0 :             niic(i,k)=0._r8
        1579           0 :             if (qi(i,k)+(cmei(i,k)+berg(i,k))*deltat.lt.0._r8) then
        1580           0 :                cmei(i,k)=(-qi(i,k)/deltat-berg(i,k))*omsm
        1581             :             end if
        1582             :          end if
        1583             :
        1584             :          ! add to cme output
        1585             :
        1586           0 :          cmeout(i,k) = cmeout(i,k)+cmei(i,k)
        1587             :
        1588             :          !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1589             :          ! droplet activation
        1590             :          ! calculate potential for droplet activation if cloud water is present
        1591             :          ! formulation from Abdul-Razzak and Ghan (2000) and Abdul-Razzak et al. (1998), AR98
        1592             :          ! number tendency (npccnin) is read in from companion routine
        1593             :
        1594             :          ! assume aerosols already activated are equal to number of existing droplets for simplicity
        1595             :          ! multiply by cloud fraction to obtain grid-average tendency
        1596             :
        1597           0 :          if (qcic(i,k).ge.qsmall) then
        1598           0 :             npccn(k) = max(0._r8,npccnin(i,k))
        1599           0 :             dum2l(i,k)=(nc(i,k)+npccn(k)*deltat)/lcldm(i,k)
        1600           0 :             dum2l(i,k)=max(dum2l(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3
        1601           0 :             ncmax = dum2l(i,k)*lcldm(i,k)
        1602             :          else
        1603           0 :             npccn(k)=0._r8
        1604           0 :             dum2l(i,k)=0._r8
        1605           0 :             ncmax = 0._r8
        1606             :          end if
        1607             :
        1608             :          !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1609             :          ! get size distribution parameters based on in-cloud cloud water/ice
        1610             :          ! these calculations also ensure consistency between number and mixing ratio
        1611             :          !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1612             :
        1613             :          !......................................................................
        1614             :          ! cloud ice
        1615             :
        1616           0 :          if (qiic(i,k).ge.qsmall) then
        1617             :
        1618             :             ! add upper limit to in-cloud number concentration to prevent numerical error
        1619           0 :             niic(i,k)=min(niic(i,k),qiic(i,k)*1.e20_r8)
        1620             :
        1621           0 :             lami(k) = (cons1*ci*niic(i,k)/qiic(i,k))**(1._r8/di)
        1622           0 :             n0i(k) = niic(i,k)*lami(k)
        1623             :
        1624             :             ! check for slope
        1625             :             ! adjust vars
        1626             :
        1627           0 :             if (lami(k).lt.lammini) then
        1628             :
        1629           0 :                lami(k) = lammini
        1630           0 :                n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
        1631           0 :                niic(i,k) = n0i(k)/lami(k)
        1632           0 :             else if (lami(k).gt.lammaxi) then
        1633           0 :                lami(k) = lammaxi
        1634           0 :                n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1)
        1635           0 :                niic(i,k) = n0i(k)/lami(k)
        1636             :             end if
        1637             :
        1638             :          else
        1639           0 :             lami(k) = 0._r8
        1640           0 :             n0i(k) = 0._r8
        1641             :          end if
        1642             :
        1643           0 :          if (qcic(i,k).ge.qsmall) then
        1644             :
        1645             :             ! add upper limit to in-cloud number concentration to prevent numerical error
        1646           0 :             ncic(i,k)=min(ncic(i,k),qcic(i,k)*1.e20_r8)
        1647             :
        1648           0 :             ncic(i,k)=max(ncic(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm
        1649             :
        1650             :             ! get pgam from fit to observations of martin et al. 1994
        1651             :
        1652           0 :             pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
        1653           0 :             pgam(k)=1._r8/(pgam(k)**2)-1._r8
        1654           0 :             pgam(k)=max(pgam(k),2._r8)
        1655           0 :             pgam(k)=min(pgam(k),15._r8)
        1656             :
        1657             :             ! calculate lamc
        1658             :
        1659             :             lamc(k) = (pi/6._r8*rhow*ncic(i,k)*gamma(pgam(k)+4._r8)/ &
        1660           0 :                  (qcic(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
        1661             :
        1662             :             ! lammin, 50 micron diameter max mean size
        1663             :
        1664           0 :             lammin = (pgam(k)+1._r8)/50.e-6_r8
        1665           0 :             lammax = (pgam(k)+1._r8)/2.e-6_r8
        1666             :
        1667           0 :             if (lamc(k).lt.lammin) then
        1668           0 :                lamc(k) = lammin
        1669             :                ncic(i,k) = 6._r8*lamc(k)**3*qcic(i,k)* &
        1670           0 :                     gamma(pgam(k)+1._r8)/(pi*rhow*gamma(pgam(k)+4._r8))
        1671           0 :             else if (lamc(k).gt.lammax) then
        1672           0 :                lamc(k) = lammax
        1673             :                ncic(i,k) = 6._r8*lamc(k)**3*qcic(i,k)* &
        1674           0 :                     gamma(pgam(k)+1._r8)/(pi*rhow*gamma(pgam(k)+4._r8))
        1675             :             end if
        1676             :
        1677             :             ! parameter to calculate droplet freezing
        1678             :
        1679           0 :             cdist1(k) = ncic(i,k)/gamma(pgam(k)+1._r8)
        1680             :
        1681             :          else
        1682           0 :             lamc(k) = 0._r8
        1683           0 :             cdist1(k) = 0._r8
        1684             :          end if
        1685             :
        1686             :          !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1687             :          !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1688             :          ! begin micropysical process calculations
        1689             :          !.................................................................
        1690             :          ! autoconversion of cloud liquid water to rain
        1691             :          ! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc
        1692             :          ! minimum qc of 1 x 10^-8 prevents floating point error
        1693             :
        1694           0 :          if (qcic(i,k).ge.1.e-8_r8) then
        1695             :
        1696             :             ! nprc is increase in rain number conc due to autoconversion
        1697             :             ! nprc1 is decrease in cloud droplet conc due to autoconversion
        1698             :
        1699             :             ! assume exponential sub-grid distribution of qc, resulting in additional
        1700             :             ! factor related to qcvar below
        1701             :
        1702             :             ! hm switch for sub-columns, don't include sub-grid qc
        1703           0 :             if (microp_uniform) then
        1704             :
        1705             :                prc(k) = 1350._r8*qcic(i,k)**2.47_r8* &
        1706           0 :                     (ncic(i,k)/1.e6_r8*rho(i,k))**(-1.79_r8)
        1707           0 :                nprc(k) = prc(k)/(4._r8/3._r8*pi*rhow*(25.e-6_r8)**3)
        1708           0 :                nprc1(k) = prc(k)/(qcic(i,k)/ncic(i,k))
        1709             :
        1710             :             else
        1711             :
        1712             :                prc(k) = cons2/(cons3*cons18)*1350._r8*qcic(i,k)**2.47_r8* &
        1713           0 :                     (ncic(i,k)/1.e6_r8*rho(i,k))**(-1.79_r8)
        1714           0 :                nprc(k) = prc(k)/cons22
        1715           0 :                nprc1(k) = prc(k)/(qcic(i,k)/ncic(i,k))
        1716             :
        1717             :             end if               ! sub-column switch
        1718             :
        1719             :          else
        1720           0 :             prc(k)=0._r8
        1721           0 :             nprc(k)=0._r8
        1722           0 :             nprc1(k)=0._r8
        1723             :          end if
        1724             :
        1725             :          ! add autoconversion to precip from above to get provisional rain mixing ratio
        1726             :          ! and number concentration (qric and nric)
        1727             :
        1728             :          ! 0.45 m/s is fallspeed of new rain drop (80 micron diameter)
        1729             :
        1730           0 :          dum=0.45_r8
        1731           0 :          dum1=0.45_r8
        1732             :
        1733           0 :          if (k.eq.top_lev) then
        1734           0 :             qric(i,k)=prc(k)*lcldm(i,k)*dz(i,k)/cldmax(i,k)/dum
        1735           0 :             nric(i,k)=nprc(k)*lcldm(i,k)*dz(i,k)/cldmax(i,k)/dum
        1736             :          else
        1737           0 :             if (qric(i,k-1).ge.qsmall) then
        1738           0 :                dum=umr(k-1)
        1739           0 :                dum1=unr(k-1)
        1740             :             end if
        1741             :
        1742             :             ! no autoconversion of rain number if rain/snow falling from above
        1743             :             ! this assumes that new drizzle drops formed by autoconversion are rapidly collected
        1744             :             ! by the existing rain/snow particles from above
        1745             :
        1746           0 :             if (qric(i,k-1).ge.1.e-9_r8.or.qniic(i,k-1).ge.1.e-9_r8) then
        1747           0 :                nprc(k)=0._r8
        1748             :             end if
        1749             :
        1750             :             qric(i,k) = (rho(i,k-1)*umr(k-1)*qric(i,k-1)*cldmax(i,k-1)+ &
        1751             :                  (rho(i,k)*dz(i,k)*((pra(k-1)+prc(k))*lcldm(i,k)+(pre(k-1)-pracs(k-1)-mnuccr(k-1))*cldmax(i,k))))&
        1752           0 :                  /(dum*rho(i,k)*cldmax(i,k))
        1753             :             nric(i,k) = (rho(i,k-1)*unr(k-1)*nric(i,k-1)*cldmax(i,k-1)+ &
        1754             :                  (rho(i,k)*dz(i,k)*(nprc(k)*lcldm(i,k)+(nsubr(k-1)-npracs(k-1)-nnuccr(k-1)+nragg(k-1))*cldmax(i,k))))&
        1755           0 :                  /(dum1*rho(i,k)*cldmax(i,k))
        1756             :
        1757             :          end if
        1758             :
        1759             :          !.......................................................................
        1760             :          ! Autoconversion of cloud ice to snow
        1761             :          ! similar to Ferrier (1994)
        1762             :
        1763           0 :          if (do_cldice) then
        1764           0 :             if (t(i,k).le.273.15_r8.and.qiic(i,k).ge.qsmall) then
        1765             :
        1766             :                ! note: assumes autoconversion timescale of 180 sec
        1767             :
        1768           0 :                nprci(k) = n0i(k)/(lami(k)*180._r8)*exp(-lami(k)*dcs)
        1769             :
        1770             :                prci(k) = pi*rhoi*n0i(k)/(6._r8*180._r8)* &
        1771             :                     (cons23/lami(k)+3._r8*cons24/lami(k)**2+ &
        1772           0 :                     6._r8*dcs/lami(k)**3+6._r8/lami(k)**4)*exp(-lami(k)*dcs)
        1773             :             else
        1774           0 :                prci(k)=0._r8
        1775           0 :                nprci(k)=0._r8
        1776             :             end if
        1777             :          else
        1778             :             ! Add in the particles that we have already converted to snow, and
        1779             :             ! don't do any further autoconversion of ice.
        1780           0 :             prci(k)  = tnd_qsnow(i, k) / cldm(i,k)
        1781           0 :             nprci(k) = tnd_nsnow(i, k) / cldm(i,k)
        1782             :          end if
        1783             :
        1784             :          ! add autoconversion to flux from level above to get provisional snow mixing ratio
        1785             :          ! and number concentration (qniic and nsic)
        1786             :
        1787           0 :          dum=(asn(i,k)*cons25)
        1788           0 :          dum1=(asn(i,k)*cons25)
        1789             :
        1790           0 :          if (k.eq.top_lev) then
        1791           0 :             qniic(i,k)=prci(k)*icldm(i,k)*dz(i,k)/cldmax(i,k)/dum
        1792           0 :             nsic(i,k)=nprci(k)*icldm(i,k)*dz(i,k)/cldmax(i,k)/dum
        1793             :          else
        1794           0 :             if (qniic(i,k-1).ge.qsmall) then
        1795           0 :                dum=ums(k-1)
        1796           0 :                dum1=uns(k-1)
        1797             :             end if
        1798             :
        1799             :             qniic(i,k) = (rho(i,k-1)*ums(k-1)*qniic(i,k-1)*cldmax(i,k-1)+ &
        1800             :                  (rho(i,k)*dz(i,k)*((prci(k)+prai(k-1)+psacws(k-1)+bergs(k-1))*icldm(i,k)+(prds(k-1)+ &
        1801             :                  pracs(k-1)+mnuccr(k-1))*cldmax(i,k))))&
        1802           0 :                  /(dum*rho(i,k)*cldmax(i,k))
        1803             :
        1804             :             nsic(i,k) = (rho(i,k-1)*uns(k-1)*nsic(i,k-1)*cldmax(i,k-1)+ &
        1805             :                  (rho(i,k)*dz(i,k)*(nprci(k)*icldm(i,k)+(nsubs(k-1)+nsagg(k-1)+nnuccr(k-1))*cldmax(i,k))))&
        1806           0 :                  /(dum1*rho(i,k)*cldmax(i,k))
        1807             :
        1808             :          end if
        1809             :
        1810             :          ! if precip mix ratio is zero so should number concentration
        1811             :
        1812           0 :          if (qniic(i,k).lt.qsmall) then
        1813           0 :             qniic(i,k)=0._r8
        1814           0 :             nsic(i,k)=0._r8
        1815             :          end if
        1816             :
        1817           0 :          if (qric(i,k).lt.qsmall) then
        1818           0 :             qric(i,k)=0._r8
        1819           0 :             nric(i,k)=0._r8
        1820             :          end if
        1821             :
        1822             :          ! make sure number concentration is a positive number to avoid
        1823             :          ! taking root of negative later
        1824             :
        1825           0 :          nric(i,k)=max(nric(i,k),0._r8)
        1826           0 :          nsic(i,k)=max(nsic(i,k),0._r8)
        1827             :
        1828             :          !.......................................................................
        1829             :          ! get size distribution parameters for precip
        1830             :          !......................................................................
        1831             :          ! rain
        1832             :
        1833           0 :          if (qric(i,k).ge.qsmall) then
        1834           0 :             lamr(k) = (pi*rhow*nric(i,k)/qric(i,k))**(1._r8/3._r8)
        1835           0 :             n0r(k) = nric(i,k)*lamr(k)
        1836             :
        1837             :             ! check for slope
        1838             :             ! adjust vars
        1839             :
        1840           0 :             if (lamr(k).lt.lamminr) then
        1841             :
        1842           0 :                lamr(k) = lamminr
        1843             :
        1844           0 :                n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
        1845           0 :                nric(i,k) = n0r(k)/lamr(k)
        1846           0 :             else if (lamr(k).gt.lammaxr) then
        1847           0 :                lamr(k) = lammaxr
        1848           0 :                n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
        1849           0 :                nric(i,k) = n0r(k)/lamr(k)
        1850             :             end if
        1851             :
        1852             :             ! provisional rain number and mass weighted mean fallspeed (m/s)
        1853             :
        1854           0 :             unr(k) = min(arn(i,k)*cons4/lamr(k)**br,9.1_r8*rhof(i,k))
        1855           0 :             umr(k) = min(arn(i,k)*cons5/(6._r8*lamr(k)**br),9.1_r8*rhof(i,k))
        1856             :
        1857             :          else
        1858           0 :             lamr(k) = 0._r8
        1859           0 :             n0r(k) = 0._r8
        1860           0 :             umr(k) = 0._r8
        1861           0 :             unr(k) = 0._r8
        1862             :          end if
        1863             :
        1864             :          !......................................................................
        1865             :          ! snow
        1866             :
        1867           0 :          if (qniic(i,k).ge.qsmall) then
        1868           0 :             lams(k) = (cons6*cs*nsic(i,k)/qniic(i,k))**(1._r8/ds)
        1869           0 :             n0s(k) = nsic(i,k)*lams(k)
        1870             :
        1871             :             ! check for slope
        1872             :             ! adjust vars
        1873             :
        1874           0 :             if (lams(k).lt.lammins) then
        1875           0 :                lams(k) = lammins
        1876           0 :                n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
        1877           0 :                nsic(i,k) = n0s(k)/lams(k)
        1878             :
        1879           0 :             else if (lams(k).gt.lammaxs) then
        1880           0 :                lams(k) = lammaxs
        1881           0 :                n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
        1882           0 :                nsic(i,k) = n0s(k)/lams(k)
        1883             :             end if
        1884             :
        1885             :             ! provisional snow number and mass weighted mean fallspeed (m/s)
        1886             :
        1887           0 :             ums(k) = min(asn(i,k)*cons8/(6._r8*lams(k)**bs),1.2_r8*rhof(i,k))
        1888           0 :             uns(k) = min(asn(i,k)*cons7/lams(k)**bs,1.2_r8*rhof(i,k))
        1889             :
        1890             :          else
        1891           0 :             lams(k) = 0._r8
        1892           0 :             n0s(k) = 0._r8
        1893           0 :             ums(k) = 0._r8
        1894           0 :             uns(k) = 0._r8
        1895             :          end if
        1896             :
        1897             :          !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        1898             :
        1899             :          ! heterogeneous freezing of cloud water
        1900             :
        1901           0 :          if (.not. use_hetfrz_classnuc) then
        1902             :
        1903           0 :             if (do_cldice .and. qcic(i,k).ge.qsmall .and. t(i,k).lt.269.15_r8) then
        1904             :
        1905             :                ! immersion freezing (Bigg, 1953)
        1906             :
        1907             :
        1908             :                ! subcolumns
        1909             :
        1910           0 :                if (microp_uniform) then
        1911             :
        1912             :                   mnuccc(k) = &
        1913             :                      pi*pi/36._r8*rhow* &
        1914             :                      cdist1(k)*gamma(7._r8+pgam(k))* &
        1915             :                      bimm*(exp(aimm*(273.15_r8-t(i,k)))-1._r8)/ &
        1916           0 :                      lamc(k)**3/lamc(k)**3
        1917             :
        1918             :                   nnuccc(k) = &
        1919             :                      pi/6._r8*cdist1(k)*gamma(pgam(k)+4._r8) &
        1920             :                      *bimm* &
        1921           0 :                      (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamc(k)**3
        1922             :
        1923             :                else
        1924             :
        1925             :                   mnuccc(k) = cons9/(cons3*cons19)* &
        1926             :                      pi*pi/36._r8*rhow* &
        1927             :                      cdist1(k)*gamma(7._r8+pgam(k))* &
        1928             :                      bimm*(exp(aimm*(273.15_r8-t(i,k)))-1._r8)/ &
        1929           0 :                      lamc(k)**3/lamc(k)**3
        1930             :
        1931             :                   nnuccc(k) = cons10/(cons3*qcvar)* &
        1932             :                      pi/6._r8*cdist1(k)*gamma(pgam(k)+4._r8) &
        1933             :                      *bimm* &
        1934           0 :                      (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamc(k)**3
        1935             :                end if           ! sub-columns
        1936             :
        1937             :
        1938             :                ! contact freezing (-40<T<-3 C) (Young, 1974) with hooks into simulated dust
        1939             :                ! dust size and number in 4 bins are read in from companion routine
        1940             :
        1941           0 :                tcnt=(270.16_r8-t(i,k))**1.3_r8
        1942           0 :                viscosity=1.8e-5_r8*(t(i,k)/298.0_r8)**0.85_r8    ! Viscosity (kg/m/s)
        1943             :                mfp=2.0_r8*viscosity/(p(i,k)  &                   ! Mean free path (m)
        1944           0 :                   *sqrt(8.0_r8*28.96e-3_r8/(pi*8.314409_r8*t(i,k))))
        1945             :
        1946           0 :                nslip1=1.0_r8+(mfp/rndst(i,k,1))*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rndst(i,k,1)/mfp))))! Slip correction factor
        1947           0 :                nslip2=1.0_r8+(mfp/rndst(i,k,2))*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rndst(i,k,2)/mfp))))
        1948           0 :                nslip3=1.0_r8+(mfp/rndst(i,k,3))*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rndst(i,k,3)/mfp))))
        1949           0 :                nslip4=1.0_r8+(mfp/rndst(i,k,4))*(1.257_r8+(0.4_r8*Exp(-(1.1_r8*rndst(i,k,4)/mfp))))
        1950             :
        1951           0 :                ndfaer1=1.381e-23_r8*t(i,k)*nslip1/(6._r8*pi*viscosity*rndst(i,k,1))  ! aerosol diffusivity (m2/s)
        1952           0 :                ndfaer2=1.381e-23_r8*t(i,k)*nslip2/(6._r8*pi*viscosity*rndst(i,k,2))
        1953           0 :                ndfaer3=1.381e-23_r8*t(i,k)*nslip3/(6._r8*pi*viscosity*rndst(i,k,3))
        1954           0 :                ndfaer4=1.381e-23_r8*t(i,k)*nslip4/(6._r8*pi*viscosity*rndst(i,k,4))
        1955             :
        1956             :
        1957           0 :                if (microp_uniform) then
        1958             :
        1959             :                   mnucct(k) = &
        1960             :                      (ndfaer1*(nacon(i,k,1)*tcnt)+ndfaer2*(nacon(i,k,2)*tcnt)+ &
        1961             :                      ndfaer3*(nacon(i,k,3)*tcnt)+ndfaer4*(nacon(i,k,4)*tcnt))*pi*pi/3._r8*rhow* &
        1962           0 :                      cdist1(k)*gamma(pgam(k)+5._r8)/lamc(k)**4
        1963             :
        1964             :                   nnucct(k) = (ndfaer1*(nacon(i,k,1)*tcnt)+ndfaer2*(nacon(i,k,2)*tcnt)+ &
        1965             :                      ndfaer3*(nacon(i,k,3)*tcnt)+ndfaer4*(nacon(i,k,4)*tcnt))*2._r8*pi*  &
        1966           0 :                      cdist1(k)*gamma(pgam(k)+2._r8)/lamc(k)
        1967             :
        1968             :                else
        1969             :
        1970             :                   mnucct(k) = gamma(qcvar+4._r8/3._r8)/(cons3*qcvar**(4._r8/3._r8))*  &
        1971             :                      (ndfaer1*(nacon(i,k,1)*tcnt)+ndfaer2*(nacon(i,k,2)*tcnt)+ &
        1972             :                      ndfaer3*(nacon(i,k,3)*tcnt)+ndfaer4*(nacon(i,k,4)*tcnt))*pi*pi/3._r8*rhow* &
        1973           0 :                      cdist1(k)*gamma(pgam(k)+5._r8)/lamc(k)**4
        1974             :
        1975             :                   nnucct(k) =  gamma(qcvar+1._r8/3._r8)/(cons3*qcvar**(1._r8/3._r8))*  &
        1976             :                      (ndfaer1*(nacon(i,k,1)*tcnt)+ndfaer2*(nacon(i,k,2)*tcnt)+ &
        1977             :                      ndfaer3*(nacon(i,k,3)*tcnt)+ndfaer4*(nacon(i,k,4)*tcnt))*2._r8*pi*  &
        1978           0 :                      cdist1(k)*gamma(pgam(k)+2._r8)/lamc(k)
        1979             :
        1980             :                end if      ! sub-column switch
        1981             :
        1982             :                ! make sure number of droplets frozen does not exceed available ice nuclei concentration
        1983             :                ! this prevents 'runaway' droplet freezing
        1984             :
        1985           0 :                if (nnuccc(k)*lcldm(i,k).gt.nnuccd(k)) then
        1986           0 :                   dum=(nnuccd(k)/(nnuccc(k)*lcldm(i,k)))
        1987             :                   ! scale mixing ratio of droplet freezing with limit
        1988           0 :                   mnuccc(k)=mnuccc(k)*dum
        1989           0 :                   nnuccc(k)=nnuccd(k)/lcldm(i,k)
        1990             :                end if
        1991             :
        1992             :             else
        1993           0 :                mnuccc(k)=0._r8
        1994           0 :                nnuccc(k)=0._r8
        1995           0 :                mnucct(k)=0._r8
        1996           0 :                nnucct(k)=0._r8
        1997             :             end if
        1998             :
        1999             :          else
        2000           0 :             if (do_cldice .and. qcic(i,k) >= qsmall) then
        2001           0 :                con1 = 1._r8/(1.333_r8*pi)**0.333_r8
        2002           0 :                r3lx = con1*(rho(i,k)*qcic(i,k)/(rhow*max(ncic(i,k)*rho(i,k), 1.0e6_r8)))**0.333_r8 ! in m
        2003           0 :                r3lx = max(4.e-6_r8, r3lx)
        2004           0 :                mi0l = 4._r8/3._r8*pi*rhow*r3lx**3_r8
        2005             :
        2006           0 :                nnuccc(k) = frzimm(i,k)*1.0e6_r8/rho(i,k)
        2007           0 :                mnuccc(k) = nnuccc(k)*mi0l
        2008             :
        2009           0 :                nnucct(k) = frzcnt(i,k)*1.0e6_r8/rho(i,k)
        2010           0 :                mnucct(k) = nnucct(k)*mi0l
        2011             :
        2012           0 :                nnudep(k) = frzdep(i,k)*1.0e6_r8/rho(i,k)
        2013           0 :                mnudep(k) = nnudep(k)*mi0
        2014             :             else
        2015           0 :                nnuccc(k) = 0._r8
        2016           0 :                mnuccc(k) = 0._r8
        2017             :
        2018           0 :                nnucct(k) = 0._r8
        2019           0 :                mnucct(k) = 0._r8
        2020             :
        2021           0 :                nnudep(k) = 0._r8
        2022           0 :                mnudep(k) = 0._r8
        2023             :             end if
        2024             :          endif
        2025             :
        2026             :
        2027             :          !.......................................................................
        2028             :          ! snow self-aggregation from passarelli, 1978, used by reisner, 1998
        2029             :          ! this is hard-wired for bs = 0.4 for now
        2030             :          ! ignore self-collection of cloud ice
        2031             :
        2032           0 :          if (qniic(i,k).ge.qsmall .and. t(i,k).le.273.15_r8) then
        2033             :             nsagg(k) = -1108._r8*asn(i,k)*Eii* &
        2034             :                  pi**((1._r8-bs)/3._r8)*rhosn**((-2._r8-bs)/3._r8)*rho(i,k)** &
        2035             :                  ((2._r8+bs)/3._r8)*qniic(i,k)**((2._r8+bs)/3._r8)* &
        2036             :                  (nsic(i,k)*rho(i,k))**((4._r8-bs)/3._r8)/ &
        2037           0 :                  (4._r8*720._r8*rho(i,k))
        2038             :          else
        2039           0 :             nsagg(k)=0._r8
        2040             :          end if
        2041             :
        2042             :          !.......................................................................
        2043             :          ! accretion of cloud droplets onto snow/graupel
        2044             :          ! here use continuous collection equation with
        2045             :          ! simple gravitational collection kernel
        2046             :          ! ignore collisions between droplets/cloud ice
        2047             :          ! since minimum size ice particle for accretion is 50 - 150 micron
        2048             :
        2049             :          ! ignore collision of snow with droplets above freezing
        2050             :
        2051           0 :          if (qniic(i,k).ge.qsmall .and. t(i,k).le.tmelt .and. &
        2052             :               qcic(i,k).ge.qsmall) then
        2053             :
        2054             :             ! put in size dependent collection efficiency
        2055             :             ! mean diameter of snow is area-weighted, since
        2056             :             ! accretion is function of crystal geometric area
        2057             :             ! collection efficiency is approximation based on stoke's law (Thompson et al. 2004)
        2058             :
        2059           0 :             dc0 = (pgam(k)+1._r8)/lamc(k)
        2060           0 :             ds0 = 1._r8/lams(k)
        2061           0 :             dum = dc0*dc0*uns(k)*rhow/(9._r8*mu(i,k)*ds0)
        2062           0 :             eci = dum*dum/((dum+0.4_r8)*(dum+0.4_r8))
        2063             :
        2064           0 :             eci = max(eci,0._r8)
        2065           0 :             eci = min(eci,1._r8)
        2066             :
        2067             :
        2068             :             ! no impact of sub-grid distribution of qc since psacws
        2069             :             ! is linear in qc
        2070             :
        2071             :             psacws(k) = pi/4._r8*asn(i,k)*qcic(i,k)*rho(i,k)* &
        2072             :                  n0s(k)*Eci*cons11/ &
        2073           0 :                  lams(k)**(bs+3._r8)
        2074             :             npsacws(k) = pi/4._r8*asn(i,k)*ncic(i,k)*rho(i,k)* &
        2075             :                  n0s(k)*Eci*cons11/ &
        2076           0 :                  lams(k)**(bs+3._r8)
        2077             :          else
        2078           0 :             psacws(k)=0._r8
        2079           0 :             npsacws(k)=0._r8
        2080             :          end if
        2081             :
        2082             :          ! add secondary ice production due to accretion of droplets by snow
        2083             :          ! (Hallet-Mossop process) (from Cotton et al., 1986)
        2084             :
        2085           0 :          if (.not. do_cldice) then
        2086           0 :             ni_secp   = 0.0_r8
        2087           0 :             nsacwi(k) = 0.0_r8
        2088           0 :             msacwi(k) = 0.0_r8
        2089           0 :          else if((t(i,k).lt.270.16_r8) .and. (t(i,k).ge.268.16_r8)) then
        2090           0 :             ni_secp   = 3.5e8_r8*(270.16_r8-t(i,k))/2.0_r8*psacws(k)
        2091           0 :             nsacwi(k) = ni_secp
        2092           0 :             msacwi(k) = min(ni_secp*mi0,psacws(k))
        2093           0 :          else if((t(i,k).lt.268.16_r8) .and. (t(i,k).ge.265.16_r8)) then
        2094           0 :             ni_secp   = 3.5e8_r8*(t(i,k)-265.16_r8)/3.0_r8*psacws(k)
        2095           0 :             nsacwi(k) = ni_secp
        2096           0 :             msacwi(k) = min(ni_secp*mi0,psacws(k))
        2097             :          else
        2098           0 :             ni_secp   = 0.0_r8
        2099           0 :             nsacwi(k) = 0.0_r8
        2100           0 :             msacwi(k) = 0.0_r8
        2101             :          endif
        2102           0 :          psacws(k) = max(0.0_r8,psacws(k)-ni_secp*mi0)
        2103             :
        2104             :          !.......................................................................
        2105             :          ! accretion of rain water by snow
        2106             :          ! formula from ikawa and saito, 1991, used by reisner et al., 1998
        2107             :
        2108           0 :          if (qric(i,k).ge.1.e-8_r8 .and. qniic(i,k).ge.1.e-8_r8 .and. &
        2109             :               t(i,k).le.273.15_r8) then
        2110             :
        2111             :             pracs(k) = pi*pi*ecr*(((1.2_r8*umr(k)-0.95_r8*ums(k))**2+ &
        2112             :                  0.08_r8*ums(k)*umr(k))**0.5_r8*rhow*rho(i,k)* &
        2113             :                  n0r(k)*n0s(k)* &
        2114             :                  (5._r8/(lamr(k)**6*lams(k))+ &
        2115             :                  2._r8/(lamr(k)**5*lams(k)**2)+ &
        2116           0 :                  0.5_r8/(lamr(k)**4*lams(k)**3)))
        2117             :
        2118             :             npracs(k) = pi/2._r8*rho(i,k)*ecr*(1.7_r8*(unr(k)-uns(k))**2+ &
        2119             :                  0.3_r8*unr(k)*uns(k))**0.5_r8*n0r(k)*n0s(k)* &
        2120             :                  (1._r8/(lamr(k)**3*lams(k))+ &
        2121             :                  1._r8/(lamr(k)**2*lams(k)**2)+ &
        2122           0 :                  1._r8/(lamr(k)*lams(k)**3))
        2123             :
        2124             :          else
        2125           0 :             pracs(k)=0._r8
        2126           0 :             npracs(k)=0._r8
        2127             :          end if
        2128             :
        2129             :          !.......................................................................
        2130             :          ! heterogeneous freezing of rain drops
        2131             :          ! follows from Bigg (1953)
        2132             :
        2133           0 :          if (t(i,k).lt.269.15_r8 .and. qric(i,k).ge.qsmall) then
        2134             :
        2135             :             mnuccr(k) = 20._r8*pi*pi*rhow*nric(i,k)*bimm* &
        2136             :                  (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamr(k)**3 &
        2137           0 :                  /lamr(k)**3
        2138             :
        2139             :             nnuccr(k) = pi*nric(i,k)*bimm* &
        2140           0 :                  (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamr(k)**3
        2141             :          else
        2142           0 :             mnuccr(k)=0._r8
        2143           0 :             nnuccr(k)=0._r8
        2144             :          end if
        2145             :
        2146             :          !.......................................................................
        2147             :          ! accretion of cloud liquid water by rain
        2148             :          ! formula from Khrouditnov and Kogan (2000)
        2149             :          ! gravitational collection kernel, droplet fall speed neglected
        2150             :
        2151           0 :          if (qric(i,k).ge.qsmall .and. qcic(i,k).ge.qsmall) then
        2152             :
        2153             :             ! include sub-grid distribution of cloud water
        2154             :
        2155             :             ! add sub-column switch
        2156             :
        2157           0 :             if (microp_uniform) then
        2158             :
        2159           0 :                pra(k) = 67._r8*(qcic(i,k)*qric(i,k))**1.15_r8
        2160           0 :                npra(k) = pra(k)/(qcic(i,k)/ncic(i,k))
        2161             :
        2162             :             else
        2163             :
        2164           0 :                pra(k) = accre_enhan(i,k)*(cons12/(cons3*cons20)*67._r8*(qcic(i,k)*qric(i,k))**1.15_r8)
        2165           0 :                npra(k) = pra(k)/(qcic(i,k)/ncic(i,k))
        2166             :
        2167             :             end if               ! sub-column switch
        2168             :
        2169             :          else
        2170           0 :             pra(k)=0._r8
        2171           0 :             npra(k)=0._r8
        2172             :          end if
        2173             :
        2174             :          !.......................................................................
        2175             :          ! Self-collection of rain drops
        2176             :          ! from Beheng(1994)
        2177             :
        2178           0 :          if (qric(i,k).ge.qsmall) then
        2179           0 :             nragg(k) = -8._r8*nric(i,k)*qric(i,k)*rho(i,k)
        2180             :          else
        2181           0 :             nragg(k)=0._r8
        2182             :          end if
        2183             :
        2184             :          !.......................................................................
        2185             :          ! Accretion of cloud ice by snow
        2186             :          ! For this calculation, it is assumed that the Vs >> Vi
        2187             :          ! and Ds >> Di for continuous collection
        2188             :
        2189             :          if (do_cldice .and. qniic(i,k).ge.qsmall.and.qiic(i,k).ge.qsmall &
        2190           0 :               .and.t(i,k).le.273.15_r8) then
        2191             :
        2192             :             prai(k) = pi/4._r8*asn(i,k)*qiic(i,k)*rho(i,k)* &
        2193             :                  n0s(k)*Eii*cons11/ &
        2194           0 :                  lams(k)**(bs+3._r8)
        2195             :             nprai(k) = pi/4._r8*asn(i,k)*niic(i,k)* &
        2196             :                  rho(i,k)*n0s(k)*Eii*cons11/ &
        2197           0 :                  lams(k)**(bs+3._r8)
        2198             :          else
        2199           0 :             prai(k)=0._r8
        2200           0 :             nprai(k)=0._r8
        2201             :          end if
        2202             :
        2203             :          !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        2204             :          ! calculate evaporation/sublimation of rain and snow
        2205             :          ! note: evaporation/sublimation occurs only in cloud-free portion of grid cell
        2206             :          ! in-cloud condensation/deposition of rain and snow is neglected
        2207             :          ! except for transfer of cloud water to snow through bergeron process
        2208             :
        2209             :          ! initialize evap/sub tendncies
        2210           0 :          pre(k)=0._r8
        2211           0 :          prds(k)=0._r8
        2212             :
        2213             :          ! evaporation of rain
        2214             :          ! only calculate if there is some precip fraction > cloud fraction
        2215             :
        2216           0 :          if (qcic(i,k)+qiic(i,k).lt.1.e-6_r8.or.cldmax(i,k).gt.lcldm(i,k)) then
        2217             :
        2218             :             ! set temporary cloud fraction to zero if cloud water + ice is very small
        2219             :             ! this will ensure that evaporation/sublimation of precip occurs over
        2220             :             ! entire grid cell, since min cloud fraction is specified otherwise
        2221           0 :             if (qcic(i,k)+qiic(i,k).lt.1.e-6_r8) then
        2222             :                dum=0._r8
        2223             :             else
        2224           0 :                dum=lcldm(i,k)
        2225             :             end if
        2226             :
        2227             :             ! saturation vapor pressure
        2228           0 :             esn=svp_water(t(i,k))
        2229           0 :             qsn=svp_to_qsat(esn, p(i,k))
        2230             :
        2231             :             ! recalculate saturation vapor pressure for liquid and ice
        2232           0 :             esl(i,k)=esn
        2233           0 :             esi(i,k)=svp_ice(t(i,k))
        2234             :             ! hm fix, make sure when above freezing that esi=esl, not active yet
        2235           0 :             if (t(i,k).gt.tmelt)esi(i,k)=esl(i,k)
        2236             :
        2237             :             ! calculate q for out-of-cloud region
        2238           0 :             qclr=(q(i,k)-dum*qsn)/(1._r8-dum)
        2239             :
        2240           0 :             if (qric(i,k).ge.qsmall) then
        2241             :
        2242           0 :                qvs=svp_to_qsat(esl(i,k), p(i,k))
        2243           0 :                dqsdt = xxlv*qvs/(rv*t(i,k)**2)
        2244           0 :                ab = 1._r8+dqsdt*xxlv/cpp
        2245             :                epsr = 2._r8*pi*n0r(k)*rho(i,k)*Dv(i,k)* &
        2246             :                     (f1r/(lamr(k)*lamr(k))+ &
        2247             :                     f2r*(arn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* &
        2248             :                     sc(i,k)**(1._r8/3._r8)*cons13/ &
        2249           0 :                     (lamr(k)**(5._r8/2._r8+br/2._r8)))
        2250             :
        2251           0 :                pre(k) = epsr*(qclr-qvs)/ab
        2252             :
        2253             :                ! only evaporate in out-of-cloud region
        2254             :                ! and distribute across cldmax
        2255           0 :                pre(k)=min(pre(k)*(cldmax(i,k)-dum),0._r8)
        2256           0 :                pre(k)=pre(k)/cldmax(i,k)
        2257           0 :                am_evp_st(i,k) = max(cldmax(i,k)-dum, 0._r8)
        2258             :             end if
        2259             :
        2260             :             ! sublimation of snow
        2261           0 :             if (qniic(i,k).ge.qsmall) then
        2262           0 :                qvi=svp_to_qsat(esi(i,k), p(i,k))
        2263           0 :                dqsidt =  xxls*qvi/(rv*t(i,k)**2)
        2264           0 :                abi = 1._r8+dqsidt*xxls/cpp
        2265             :                epss = 2._r8*pi*n0s(k)*rho(i,k)*Dv(i,k)* &
        2266             :                     (f1s/(lams(k)*lams(k))+ &
        2267             :                     f2s*(asn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* &
        2268             :                     sc(i,k)**(1._r8/3._r8)*cons14/ &
        2269           0 :                     (lams(k)**(5._r8/2._r8+bs/2._r8)))
        2270           0 :                prds(k) = epss*(qclr-qvi)/abi
        2271             :
        2272             :                ! only sublimate in out-of-cloud region and distribute over cldmax
        2273           0 :                prds(k)=min(prds(k)*(cldmax(i,k)-dum),0._r8)
        2274           0 :                prds(k)=prds(k)/cldmax(i,k)
        2275           0 :                am_evp_st(i,k) = max(cldmax(i,k)-dum, 0._r8)
        2276             :             end if
        2277             :
        2278             :             ! make sure RH not pushed above 100% due to rain evaporation/snow sublimation
        2279             :             ! get updated RH at end of time step based on cloud water/ice condensation/evap
        2280             :
        2281           0 :             qtmp=q(i,k)-(cmei(i,k)+(pre(k)+prds(k))*cldmax(i,k))*deltat
        2282             :             ttmp=t(i,k)+((pre(k)*cldmax(i,k))*xxlv+ &
        2283           0 :                  (cmei(i,k)+prds(k)*cldmax(i,k))*xxls)*deltat/cpp
        2284             :
        2285             :             !limit range of temperatures!
        2286           0 :             ttmp=max(180._r8,min(ttmp,323._r8))
        2287             :
        2288           0 :             esn=svp_water(ttmp)  ! use rhw to allow ice supersaturation
        2289           0 :             qsn=svp_to_qsat(esn, p(i,k))
        2290             :
        2291             :             ! modify precip evaporation rate if q > qsat
        2292           0 :             if (qtmp.gt.qsn) then
        2293           0 :                if (pre(k)+prds(k).lt.-1.e-20_r8) then
        2294           0 :                   dum1=pre(k)/(pre(k)+prds(k))
        2295             :                   ! recalculate q and t after cloud water cond but without precip evap
        2296           0 :                   qtmp=q(i,k)-(cmei(i,k))*deltat
        2297           0 :                   ttmp=t(i,k)+(cmei(i,k)*xxls)*deltat/cpp
        2298           0 :                   esn=svp_water(ttmp) ! use rhw to allow ice supersaturation
        2299           0 :                   qsn=svp_to_qsat(esn, p(i,k))
        2300           0 :                   dum=(qtmp-qsn)/(1._r8 + cons27*qsn/(cpp*rv*ttmp**2))
        2301           0 :                   dum=min(dum,0._r8)
        2302             :
        2303             :                   ! modify rates if needed, divide by cldmax to get local (in-precip) value
        2304           0 :                   pre(k)=dum*dum1/deltat/cldmax(i,k)
        2305             :
        2306             :                   ! do separately using RHI for prds....
        2307           0 :                   esn=svp_ice(ttmp) ! use rhi to allow ice supersaturation
        2308           0 :                   qsn=svp_to_qsat(esn, p(i,k))
        2309           0 :                   dum=(qtmp-qsn)/(1._r8 + cons28*qsn/(cpp*rv*ttmp**2))
        2310           0 :                   dum=min(dum,0._r8)
        2311             :
        2312             :                   ! modify rates if needed, divide by cldmax to get local (in-precip) value
        2313           0 :                   prds(k)=dum*(1._r8-dum1)/deltat/cldmax(i,k)
        2314             :                end if
        2315             :             end if
        2316             :          end if
        2317             :
        2318             :          ! bergeron process - evaporation of droplets and deposition onto snow
        2319             :
        2320           0 :          if (qniic(i,k).ge.qsmall.and.qcic(i,k).ge.qsmall.and.t(i,k).lt.tmelt) then
        2321           0 :             qvi=svp_to_qsat(esi(i,k), p(i,k))
        2322           0 :             qvs=svp_to_qsat(esl(i,k), p(i,k))
        2323           0 :             dqsidt =  xxls*qvi/(rv*t(i,k)**2)
        2324           0 :             abi = 1._r8+dqsidt*xxls/cpp
        2325             :             epss = 2._r8*pi*n0s(k)*rho(i,k)*Dv(i,k)* &
        2326             :                  (f1s/(lams(k)*lams(k))+ &
        2327             :                  f2s*(asn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* &
        2328             :                  sc(i,k)**(1._r8/3._r8)*cons14/ &
        2329           0 :                  (lams(k)**(5._r8/2._r8+bs/2._r8)))
        2330           0 :             bergs(k)=epss*(qvs-qvi)/abi
        2331             :          else
        2332           0 :             bergs(k)=0._r8
        2333             :          end if
        2334             :
        2335             :          !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        2336             :          ! conservation to ensure no negative values of cloud water/precipitation
        2337             :          ! in case microphysical process rates are large
        2338             :
        2339             :          ! make sure and use end-of-time step values for cloud water, ice, due
        2340             :          ! condensation/deposition
        2341             :
        2342             :          ! note: for check on conservation, processes are multiplied by omsm
        2343             :          ! to prevent problems due to round off error
        2344             :
        2345             :          ! include mixing timescale  (mtime)
        2346             :
        2347           0 :          qce=(qc(i,k) - berg(i,k)*deltat)
        2348           0 :          nce=(nc(i,k)+npccn(k)*deltat*mtime)
        2349           0 :          qie=(qi(i,k)+(cmei(i,k)+berg(i,k))*deltat)
        2350           0 :          nie=(ni(i,k)+nnuccd(k)*deltat*mtime)
        2351             :
        2352             :          ! conservation of qc
        2353             :
        2354             :          dum = (prc(k)+pra(k)+mnuccc(k)+mnucct(k)+msacwi(k)+ &
        2355           0 :               psacws(k)+bergs(k))*lcldm(i,k)*deltat
        2356             :
        2357           0 :          if (dum.gt.qce) then
        2358           0 :             ratio = qce/deltat/lcldm(i,k)/(prc(k)+pra(k)+mnuccc(k)+mnucct(k)+msacwi(k)+psacws(k)+bergs(k))*omsm
        2359             :
        2360           0 :             prc(k) = prc(k)*ratio
        2361           0 :             pra(k) = pra(k)*ratio
        2362           0 :             mnuccc(k) = mnuccc(k)*ratio
        2363           0 :             mnucct(k) = mnucct(k)*ratio
        2364           0 :             msacwi(k) = msacwi(k)*ratio
        2365           0 :             psacws(k) = psacws(k)*ratio
        2366           0 :             bergs(k) = bergs(k)*ratio
        2367             :          end if
        2368             :
        2369             :          ! conservation of nc
        2370             :
        2371             :          dum = (nprc1(k)+npra(k)+nnuccc(k)+nnucct(k)+ &
        2372           0 :               npsacws(k)-nsubc(k))*lcldm(i,k)*deltat
        2373             :
        2374           0 :          if (dum.gt.nce) then
        2375             :             ratio = nce/deltat/((nprc1(k)+npra(k)+nnuccc(k)+nnucct(k)+&
        2376           0 :                  npsacws(k)-nsubc(k))*lcldm(i,k))*omsm
        2377             :
        2378           0 :             nprc1(k) = nprc1(k)*ratio
        2379           0 :             npra(k) = npra(k)*ratio
        2380           0 :             nnuccc(k) = nnuccc(k)*ratio
        2381           0 :             nnucct(k) = nnucct(k)*ratio
        2382           0 :             npsacws(k) = npsacws(k)*ratio
        2383           0 :             nsubc(k)=nsubc(k)*ratio
        2384             :          end if
        2385             :
        2386             :          ! conservation of qi
        2387             :
        2388           0 :          if (do_cldice) then
        2389             :
        2390           0 :             frztmp = -mnuccc(k) - mnucct(k) - msacwi(k)
        2391           0 :             if (use_hetfrz_classnuc) frztmp = -mnuccc(k)-mnucct(k)-mnudep(k)-msacwi(k)
        2392           0 :             dum = ( frztmp*lcldm(i,k) + (prci(k)+prai(k))*icldm(i,k) )*deltat
        2393             :
        2394           0 :             if (dum.gt.qie) then
        2395             :
        2396           0 :                frztmp = mnuccc(k) + mnucct(k) + msacwi(k)
        2397           0 :                if (use_hetfrz_classnuc) frztmp = mnuccc(k) + mnucct(k) + mnudep(k) + msacwi(k)
        2398           0 :                ratio = (qie/deltat + frztmp*lcldm(i,k))/((prci(k)+prai(k))*icldm(i,k))*omsm
        2399           0 :                prci(k) = prci(k)*ratio
        2400           0 :                prai(k) = prai(k)*ratio
        2401             :             end if
        2402             :
        2403             :             ! conservation of ni
        2404           0 :             frztmp = -nnucct(k) - nsacwi(k)
        2405           0 :             if (use_hetfrz_classnuc) frztmp = -nnucct(k) - nnuccc(k) - nnudep(k) - nsacwi(k)
        2406           0 :             dum = ( frztmp*lcldm(i,k) + (nprci(k)+nprai(k)-nsubi(k))*icldm(i,k) )*deltat
        2407             :
        2408           0 :             if (dum.gt.nie) then
        2409             :
        2410           0 :                frztmp = nnucct(k) + nsacwi(k)
        2411           0 :                if (use_hetfrz_classnuc) frztmp = nnucct(k) + nnuccc(k) + nnudep(k) + nsacwi(k)
        2412             :                ratio = (nie/deltat + frztmp*lcldm(i,k))/ &
        2413           0 :                      ((nprci(k)+nprai(k)-nsubi(k))*icldm(i,k))*omsm
        2414           0 :                nprci(k) = nprci(k)*ratio
        2415           0 :                nprai(k) = nprai(k)*ratio
        2416           0 :                nsubi(k) = nsubi(k)*ratio
        2417             :             end if
        2418             :          end if
        2419             :
        2420             :          ! for precipitation conservation, use logic that vertical integral
        2421             :          ! of tendency from current level to top of model (i.e., qrtot) cannot be negative
        2422             :
        2423             :          ! conservation of rain mixing rat
        2424             :
        2425           0 :          if (((prc(k)+pra(k))*lcldm(i,k)+(-mnuccr(k)+pre(k)-pracs(k))*&
        2426             :               cldmax(i,k))*dz(i,k)*rho(i,k)+qrtot.lt.0._r8) then
        2427             :
        2428           0 :             if (-pre(k)+pracs(k)+mnuccr(k).ge.qsmall) then
        2429             :
        2430             :                ratio = (qrtot/(dz(i,k)*rho(i,k))+(prc(k)+pra(k))*lcldm(i,k))/&
        2431           0 :                     ((-pre(k)+pracs(k)+mnuccr(k))*cldmax(i,k))*omsm
        2432             :
        2433           0 :                pre(k) = pre(k)*ratio
        2434           0 :                pracs(k) = pracs(k)*ratio
        2435           0 :                mnuccr(k) = mnuccr(k)*ratio
        2436             :             end if
        2437             :          end if
        2438             :
        2439             :          ! conservation of nr
        2440             :          ! for now neglect evaporation of nr
        2441           0 :          nsubr(k)=0._r8
        2442             :
        2443           0 :          if ((nprc(k)*lcldm(i,k)+(-nnuccr(k)+nsubr(k)-npracs(k)&
        2444             :               +nragg(k))*cldmax(i,k))*dz(i,k)*rho(i,k)+nrtot.lt.0._r8) then
        2445             :
        2446           0 :             if (-nsubr(k)-nragg(k)+npracs(k)+nnuccr(k).ge.qsmall) then
        2447             :
        2448             :                ratio = (nrtot/(dz(i,k)*rho(i,k))+nprc(k)*lcldm(i,k))/&
        2449           0 :                     ((-nsubr(k)-nragg(k)+npracs(k)+nnuccr(k))*cldmax(i,k))*omsm
        2450             :
        2451           0 :                nsubr(k) = nsubr(k)*ratio
        2452           0 :                npracs(k) = npracs(k)*ratio
        2453           0 :                nnuccr(k) = nnuccr(k)*ratio
        2454           0 :                nragg(k) = nragg(k)*ratio
        2455             :             end if
        2456             :          end if
        2457             :
        2458             :          ! conservation of snow mix ratio
        2459             :
        2460           0 :          if (((bergs(k)+psacws(k))*lcldm(i,k)+(prai(k)+prci(k))*icldm(i,k)+(pracs(k)+&
        2461             :               mnuccr(k)+prds(k))*cldmax(i,k))*dz(i,k)*rho(i,k)+qstot.lt.0._r8) then
        2462             :
        2463           0 :             if (-prds(k).ge.qsmall) then
        2464             :
        2465             :                ratio = (qstot/(dz(i,k)*rho(i,k))+(bergs(k)+psacws(k))*lcldm(i,k)+(prai(k)+prci(k))*icldm(i,k)+&
        2466           0 :                     (pracs(k)+mnuccr(k))*cldmax(i,k))/(-prds(k)*cldmax(i,k))*omsm
        2467             :
        2468           0 :                prds(k) = prds(k)*ratio
        2469             :             end if
        2470             :          end if
        2471             :
        2472             :          ! conservation of ns
        2473             :
        2474             :          ! calculate loss of number due to sublimation
        2475             :          ! for now neglect sublimation of ns
        2476           0 :          nsubs(k)=0._r8
        2477             :
        2478           0 :          if ((nprci(k)*icldm(i,k)+(nnuccr(k)+nsubs(k)+nsagg(k))*cldmax(i,k))*&
        2479             :               dz(i,k)*rho(i,k)+nstot.lt.0._r8) then
        2480             :
        2481           0 :             if (-nsubs(k)-nsagg(k).ge.qsmall) then
        2482             :
        2483             :                ratio = (nstot/(dz(i,k)*rho(i,k))+nprci(k)*icldm(i,k)+&
        2484           0 :                     nnuccr(k)*cldmax(i,k))/((-nsubs(k)-nsagg(k))*cldmax(i,k))*omsm
        2485             :
        2486           0 :                nsubs(k) = nsubs(k)*ratio
        2487           0 :                nsagg(k) = nsagg(k)*ratio
        2488             :             end if
        2489             :          end if
        2490             :
        2491             :          ! get tendencies due to microphysical conversion processes
        2492             :          ! note: tendencies are multiplied by appropaiate cloud/precip
        2493             :          ! fraction to get grid-scale values
        2494             :          ! note: cmei is already grid-average values
        2495             :
        2496           0 :          qvlat(i,k) = qvlat(i,k)-(pre(k)+prds(k))*cldmax(i,k)-cmei(i,k)
        2497             :
        2498             :          tlat(i,k) = tlat(i,k)+((pre(k)*cldmax(i,k)) &
        2499             :               *xxlv+(prds(k)*cldmax(i,k)+cmei(i,k))*xxls+ &
        2500             :               ((bergs(k)+psacws(k)+mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k)+(mnuccr(k)+ &
        2501           0 :               pracs(k))*cldmax(i,k)+berg(i,k))*xlf)
        2502             :
        2503             :          qctend(i,k) = qctend(i,k)+ &
        2504             :               (-pra(k)-prc(k)-mnuccc(k)-mnucct(k)-msacwi(k)- &
        2505           0 :               psacws(k)-bergs(k))*lcldm(i,k)-berg(i,k)
        2506             :
        2507           0 :          if (do_cldice) then
        2508             :
        2509           0 :             frztmp = mnuccc(k) + mnucct(k) + msacwi(k)
        2510           0 :             if (use_hetfrz_classnuc) frztmp = mnuccc(k) + mnucct(k) + mnudep(k) + msacwi(k)
        2511             :             qitend(i,k) = qitend(i,k) + frztmp*lcldm(i,k) + &
        2512           0 :                (-prci(k)-prai(k))*icldm(i,k) + cmei(i,k) + berg(i,k)
        2513             :
        2514             :          end if
        2515             :
        2516             :          qrtend(i,k) = qrtend(i,k)+ &
        2517             :               (pra(k)+prc(k))*lcldm(i,k)+(pre(k)-pracs(k)- &
        2518           0 :               mnuccr(k))*cldmax(i,k)
        2519             :
        2520             :          qnitend(i,k) = qnitend(i,k)+ &
        2521             :               (prai(k)+prci(k))*icldm(i,k)+(psacws(k)+bergs(k))*lcldm(i,k)+(prds(k)+ &
        2522           0 :               pracs(k)+mnuccr(k))*cldmax(i,k)
        2523             :
        2524             :          ! add output for cmei (accumulate)
        2525           0 :          cmeiout(i,k) = cmeiout(i,k) + cmei(i,k)
        2526             :
        2527             :          ! assign variables for trop_mozart, these are grid-average
        2528             :          ! evaporation/sublimation is stored here as positive term
        2529             :
        2530           0 :          evapsnow(i,k) = evapsnow(i,k)-prds(k)*cldmax(i,k)
        2531           0 :          nevapr(i,k) = nevapr(i,k)-pre(k)*cldmax(i,k)
        2532           0 :          nevapr2(i,k) = nevapr2(i,k)-pre(k)*cldmax(i,k)
        2533             :
        2534             :          ! change to make sure prain is positive: do not remove snow from
        2535             :          ! prain used for wet deposition
        2536             :          prain(i,k) = prain(i,k)+(pra(k)+prc(k))*lcldm(i,k)+(-pracs(k)- &
        2537           0 :               mnuccr(k))*cldmax(i,k)
        2538             :          prodsnow(i,k) = prodsnow(i,k)+(prai(k)+prci(k))*icldm(i,k)+(psacws(k)+bergs(k))*lcldm(i,k)+(&
        2539           0 :               pracs(k)+mnuccr(k))*cldmax(i,k)
        2540             :
        2541             :          ! following are used to calculate 1st order conversion rate of cloud water
        2542             :          !    to rain and snow (1/s), for later use in aerosol wet removal routine
        2543             :          ! previously, wetdepa used (prain/qc) for this, and the qc in wetdepa may be smaller than the qc
        2544             :          !    used to calculate pra, prc, ... in this routine
        2545             :          ! qcsinksum_rate1ord = sum over iterations{ rate of direct transfer of cloud water to rain & snow }
        2546             :          !                      (no cloud ice or bergeron terms)
        2547             :          ! qcsum_rate1ord     = sum over iterations{ qc used in calculation of the transfer terms }
        2548             :
        2549           0 :          qcsinksum_rate1ord(k) = qcsinksum_rate1ord(k) + (pra(k)+prc(k)+psacws(k))*lcldm(i,k)
        2550           0 :          qcsum_rate1ord(k) = qcsum_rate1ord(k) + qc(i,k)
        2551             :
        2552             :          ! microphysics output, note this is grid-averaged
        2553           0 :          prao(i,k)=prao(i,k)+pra(k)*lcldm(i,k)
        2554           0 :          prco(i,k)=prco(i,k)+prc(k)*lcldm(i,k)
        2555           0 :          mnuccco(i,k)=mnuccco(i,k)+mnuccc(k)*lcldm(i,k)
        2556           0 :          mnuccto(i,k)=mnuccto(i,k)+mnucct(k)*lcldm(i,k)
        2557           0 :          mnuccdo(i,k)=mnuccdo(i,k)+mnuccd(k)*lcldm(i,k)
        2558           0 :          msacwio(i,k)=msacwio(i,k)+msacwi(k)*lcldm(i,k)
        2559           0 :          psacwso(i,k)=psacwso(i,k)+psacws(k)*lcldm(i,k)
        2560           0 :          bergso(i,k)=bergso(i,k)+bergs(k)*lcldm(i,k)
        2561           0 :          bergo(i,k)=bergo(i,k)+berg(i,k)
        2562           0 :          prcio(i,k)=prcio(i,k)+prci(k)*icldm(i,k)
        2563           0 :          praio(i,k)=praio(i,k)+prai(k)*icldm(i,k)
        2564           0 :          mnuccro(i,k)=mnuccro(i,k)+mnuccr(k)*cldmax(i,k)
        2565           0 :          pracso (i,k)=pracso (i,k)+pracs (k)*cldmax(i,k)
        2566             :
        2567             :          ! multiply activation/nucleation by mtime to account for fast timescale
        2568             :
        2569             :          nctend(i,k) = nctend(i,k)+ npccn(k)*mtime+&
        2570             :               (-nnuccc(k)-nnucct(k)-npsacws(k)+nsubc(k) &
        2571           0 :               -npra(k)-nprc1(k))*lcldm(i,k)
        2572             :
        2573           0 :          if (do_cldice) then
        2574             :
        2575           0 :             frztmp = nnucct(k) + nsacwi(k)
        2576           0 :             if (use_hetfrz_classnuc) frztmp = nnucct(k) + nnuccc(k) + nnudep(k) + nsacwi(k)
        2577             :             nitend(i,k) = nitend(i,k) + nnuccd(k)*mtime + &
        2578           0 :                   frztmp*lcldm(i,k) + (nsubi(k)-nprci(k)-nprai(k))*icldm(i,k)
        2579             :
        2580             :          end if
        2581             :
        2582             :          nstend(i,k) = nstend(i,k)+(nsubs(k)+ &
        2583           0 :               nsagg(k)+nnuccr(k))*cldmax(i,k)+nprci(k)*icldm(i,k)
        2584             :
        2585             :          nrtend(i,k) = nrtend(i,k)+ &
        2586             :               nprc(k)*lcldm(i,k)+(nsubr(k)-npracs(k)-nnuccr(k) &
        2587           0 :               +nragg(k))*cldmax(i,k)
        2588             :
        2589             :          ! make sure that nc and ni at advanced time step do not exceed
        2590             :          ! maximum (existing N + source terms*dt), which is possible due to
        2591             :          ! fast nucleation timescale
        2592             :
        2593           0 :          if (nctend(i,k).gt.0._r8.and.nc(i,k)+nctend(i,k)*deltat.gt.ncmax) then
        2594           0 :             nctend(i,k)=max(0._r8,(ncmax-nc(i,k))/deltat)
        2595             :          end if
        2596             :
        2597           0 :          if (do_cldice .and. nitend(i,k).gt.0._r8.and.ni(i,k)+nitend(i,k)*deltat.gt.nimax) then
        2598           0 :             nitend(i,k)=max(0._r8,(nimax-ni(i,k))/deltat)
        2599             :          end if
        2600             :
        2601             :          ! get final values for precipitation q and N, based on
        2602             :          ! flux of precip from above, source/sink term, and terminal fallspeed
        2603             :          ! see eq. 15-16 in MG2008
        2604             :
        2605             :          ! rain
        2606             :
        2607           0 :          if (qric(i,k).ge.qsmall) then
        2608           0 :             if (k.eq.top_lev) then
        2609           0 :                qric(i,k)=qrtend(i,k)*dz(i,k)/cldmax(i,k)/umr(k)
        2610           0 :                nric(i,k)=nrtend(i,k)*dz(i,k)/cldmax(i,k)/unr(k)
        2611             :             else
        2612           0 :                qric(i,k) = (rho(i,k-1)*umr(k-1)*qric(i,k-1)*cldmax(i,k-1)+ &
        2613           0 :                     (rho(i,k)*dz(i,k)*qrtend(i,k)))/(umr(k)*rho(i,k)*cldmax(i,k))
        2614             :                nric(i,k) = (rho(i,k-1)*unr(k-1)*nric(i,k-1)*cldmax(i,k-1)+ &
        2615           0 :                     (rho(i,k)*dz(i,k)*nrtend(i,k)))/(unr(k)*rho(i,k)*cldmax(i,k))
        2616             :
        2617             :             end if
        2618             :          else
        2619           0 :             qric(i,k)=0._r8
        2620           0 :             nric(i,k)=0._r8
        2621             :          end if
        2622             :
        2623             :          ! snow
        2624             :
        2625           0 :          if (qniic(i,k).ge.qsmall) then
        2626           0 :             if (k.eq.top_lev) then
        2627           0 :                qniic(i,k)=qnitend(i,k)*dz(i,k)/cldmax(i,k)/ums(k)
        2628           0 :                nsic(i,k)=nstend(i,k)*dz(i,k)/cldmax(i,k)/uns(k)
        2629             :             else
        2630           0 :                qniic(i,k) = (rho(i,k-1)*ums(k-1)*qniic(i,k-1)*cldmax(i,k-1)+ &
        2631           0 :                     (rho(i,k)*dz(i,k)*qnitend(i,k)))/(ums(k)*rho(i,k)*cldmax(i,k))
        2632             :                nsic(i,k) = (rho(i,k-1)*uns(k-1)*nsic(i,k-1)*cldmax(i,k-1)+ &
        2633           0 :                     (rho(i,k)*dz(i,k)*nstend(i,k)))/(uns(k)*rho(i,k)*cldmax(i,k))
        2634             :             end if
        2635             :          else
        2636           0 :             qniic(i,k)=0._r8
        2637           0 :             nsic(i,k)=0._r8
        2638             :          end if
        2639             :
        2640             :          ! calculate precipitation flux at surface
        2641             :          ! divide by density of water to get units of m/s
        2642             :
        2643             :          prect(i) = prect(i)+(qrtend(i,k)*dz(i,k)*rho(i,k)+&
        2644           0 :               qnitend(i,k)*dz(i,k)*rho(i,k))/rhow
        2645           0 :          preci(i) = preci(i)+qnitend(i,k)*dz(i,k)*rho(i,k)/rhow
        2646             :
        2647             :          ! convert rain rate from m/s to mm/hr
        2648             :
        2649           0 :          rainrt(i,k)=qric(i,k)*rho(i,k)*umr(k)/rhow*3600._r8*1000._r8
        2650             :
        2651             :          ! vertically-integrated precip source/sink terms (note: grid-averaged)
        2652             :
        2653           0 :          qrtot = max(qrtot+qrtend(i,k)*dz(i,k)*rho(i,k),0._r8)
        2654           0 :          qstot = max(qstot+qnitend(i,k)*dz(i,k)*rho(i,k),0._r8)
        2655           0 :          nrtot = max(nrtot+nrtend(i,k)*dz(i,k)*rho(i,k),0._r8)
        2656           0 :          nstot = max(nstot+nstend(i,k)*dz(i,k)*rho(i,k),0._r8)
        2657             :
        2658             :          ! calculate melting and freezing of precip
        2659             :
        2660             :          ! melt snow at +2 C
        2661             :
        2662           0 :          if (t(i,k)+tlat(i,k)/cpp*deltat > 275.15_r8) then
        2663           0 :             if (qstot > 0._r8) then
        2664             :
        2665             :                ! make sure melting snow doesn't reduce temperature below threshold
        2666           0 :                dum = -xlf/cpp*qstot/(dz(i,k)*rho(i,k))
        2667           0 :                if (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt.275.15_r8) then
        2668           0 :                   dum = (t(i,k)+tlat(i,k)/cpp*deltat-275.15_r8)*cpp/xlf
        2669           0 :                   dum = dum/(xlf/cpp*qstot/(dz(i,k)*rho(i,k)))
        2670           0 :                   dum = max(0._r8,dum)
        2671           0 :                   dum = min(1._r8,dum)
        2672             :                else
        2673             :                   dum = 1._r8
        2674             :                end if
        2675             :
        2676           0 :                qric(i,k)=qric(i,k)+dum*qniic(i,k)
        2677           0 :                nric(i,k)=nric(i,k)+dum*nsic(i,k)
        2678           0 :                qniic(i,k)=(1._r8-dum)*qniic(i,k)
        2679           0 :                nsic(i,k)=(1._r8-dum)*nsic(i,k)
        2680             :                ! heating tendency
        2681           0 :                tmp=-xlf*dum*qstot/(dz(i,k)*rho(i,k))
        2682           0 :                meltsdt(i,k)=meltsdt(i,k) + tmp
        2683             :
        2684           0 :                tlat(i,k)=tlat(i,k)+tmp
        2685           0 :                qrtot=qrtot+dum*qstot
        2686           0 :                nrtot=nrtot+dum*nstot
        2687           0 :                qstot=(1._r8-dum)*qstot
        2688           0 :                nstot=(1._r8-dum)*nstot
        2689           0 :                preci(i)=(1._r8-dum)*preci(i)
        2690             :             end if
        2691             :          end if
        2692             :
        2693             :          ! freeze all rain at -5C for Arctic
        2694             :
        2695           0 :          if (t(i,k)+tlat(i,k)/cpp*deltat < (tmelt - 5._r8)) then
        2696             :
        2697           0 :             if (qrtot > 0._r8) then
        2698             :
        2699             :                ! make sure freezing rain doesn't increase temperature above threshold
        2700           0 :                dum = xlf/cpp*qrtot/(dz(i,k)*rho(i,k))
        2701           0 :                if (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.(tmelt - 5._r8)) then
        2702           0 :                   dum = -(t(i,k)+tlat(i,k)/cpp*deltat-(tmelt-5._r8))*cpp/xlf
        2703           0 :                   dum = dum/(xlf/cpp*qrtot/(dz(i,k)*rho(i,k)))
        2704           0 :                   dum = max(0._r8,dum)
        2705           0 :                   dum = min(1._r8,dum)
        2706             :                else
        2707             :                   dum = 1._r8
        2708             :                end if
        2709             :
        2710           0 :                qniic(i,k)=qniic(i,k)+dum*qric(i,k)
        2711           0 :                nsic(i,k)=nsic(i,k)+dum*nric(i,k)
        2712           0 :                qric(i,k)=(1._r8-dum)*qric(i,k)
        2713           0 :                nric(i,k)=(1._r8-dum)*nric(i,k)
        2714             :                ! heating tendency
        2715           0 :                tmp = xlf*dum*qrtot/(dz(i,k)*rho(i,k))
        2716           0 :                frzrdt(i,k)=frzrdt(i,k) + tmp
        2717             :
        2718           0 :                tlat(i,k)=tlat(i,k)+tmp
        2719           0 :                qstot=qstot+dum*qrtot
        2720           0 :                qrtot=(1._r8-dum)*qrtot
        2721           0 :                nstot=nstot+dum*nrtot
        2722           0 :                nrtot=(1._r8-dum)*nrtot
        2723           0 :                preci(i)=preci(i)+dum*(prect(i)-preci(i))
        2724             :             end if
        2725             :          end if
        2726             :
        2727             :          ! Precip Flux Calculation (Diagnostic)
        2728           0 :          rflx(i,k+1)=(prect(i)-preci(i)) * rhow
        2729           0 :          sflx(i,k+1)=preci(i) * rhow
        2730             :
        2731             :          ! if rain/snow mix ratio is zero so should number concentration
        2732             :
        2733           0 :          if (qniic(i,k).lt.qsmall) then
        2734           0 :             qniic(i,k)=0._r8
        2735           0 :             nsic(i,k)=0._r8
        2736             :          end if
        2737             :
        2738           0 :          if (qric(i,k).lt.qsmall) then
        2739           0 :             qric(i,k)=0._r8
        2740           0 :             nric(i,k)=0._r8
        2741             :          end if
        2742             :
        2743             :          ! make sure number concentration is a positive number to avoid
        2744             :          ! taking root of negative
        2745             :
        2746           0 :          nric(i,k)=max(nric(i,k),0._r8)
        2747           0 :          nsic(i,k)=max(nsic(i,k),0._r8)
        2748             :
        2749             :          !.......................................................................
        2750             :          ! get size distribution parameters for fallspeed calculations
        2751             :          !......................................................................
        2752             :          ! rain
        2753             :
        2754           0 :          if (qric(i,k).ge.qsmall) then
        2755           0 :             lamr(k) = (pi*rhow*nric(i,k)/qric(i,k))**(1._r8/3._r8)
        2756           0 :             n0r(k) = nric(i,k)*lamr(k)
        2757             :
        2758             :             ! check for slope
        2759             :             ! change lammax and lammin for rain and snow
        2760             :             ! adjust vars
        2761             :
        2762           0 :             if (lamr(k).lt.lamminr) then
        2763             :
        2764           0 :                lamr(k) = lamminr
        2765             :
        2766           0 :                n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
        2767           0 :                nric(i,k) = n0r(k)/lamr(k)
        2768           0 :             else if (lamr(k).gt.lammaxr) then
        2769           0 :                lamr(k) = lammaxr
        2770           0 :                n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow)
        2771           0 :                nric(i,k) = n0r(k)/lamr(k)
        2772             :             end if
        2773             :
        2774             :
        2775             :             ! 'final' values of number and mass weighted mean fallspeed for rain (m/s)
        2776             :
        2777           0 :             unr(k) = min(arn(i,k)*cons4/lamr(k)**br,9.1_r8*rhof(i,k))
        2778           0 :             umr(k) = min(arn(i,k)*cons5/(6._r8*lamr(k)**br),9.1_r8*rhof(i,k))
        2779             :
        2780             :          else
        2781           0 :             lamr(k) = 0._r8
        2782           0 :             n0r(k) = 0._r8
        2783           0 :             umr(k)=0._r8
        2784           0 :             unr(k)=0._r8
        2785             :          end if
        2786             :
        2787             :          !calculate mean size of combined rain and snow
        2788             :
        2789           0 :          if (lamr(k).gt.0._r8) then
        2790           0 :             Artmp = n0r(k) * pi / (2._r8 * lamr(k)**3._r8)
        2791             :          else
        2792             :             Artmp = 0._r8
        2793             :          endif
        2794             :
        2795           0 :          if (lamc(k).gt.0._r8) then
        2796           0 :             Actmp = cdist1(k) * pi * gamma(pgam(k)+3._r8)/(4._r8 * lamc(k)**2._r8)
        2797             :          else
        2798             :             Actmp = 0._r8
        2799             :          endif
        2800             :
        2801           0 :          if (Actmp.gt.0_r8.or.Artmp.gt.0) then
        2802           0 :             rercld(i,k)=rercld(i,k) + 3._r8 *(qric(i,k) + qcic(i,k)) / (4._r8 * rhow * (Actmp + Artmp))
        2803           0 :             arcld(i,k)=arcld(i,k)+1._r8
        2804             :          endif
        2805             :
        2806             :          !......................................................................
        2807             :          ! snow
        2808             :
        2809           0 :          if (qniic(i,k).ge.qsmall) then
        2810             :             lams(k) = (cons6*cs*nsic(i,k)/ &
        2811           0 :                  qniic(i,k))**(1._r8/ds)
        2812           0 :             n0s(k) = nsic(i,k)*lams(k)
        2813             :
        2814             :             ! check for slope
        2815             :             ! adjust vars
        2816             :
        2817           0 :             if (lams(k).lt.lammins) then
        2818           0 :                lams(k) = lammins
        2819           0 :                n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
        2820           0 :                nsic(i,k) = n0s(k)/lams(k)
        2821             :
        2822           0 :             else if (lams(k).gt.lammaxs) then
        2823           0 :                lams(k) = lammaxs
        2824           0 :                n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6)
        2825           0 :                nsic(i,k) = n0s(k)/lams(k)
        2826             :             end if
        2827             :
        2828             :             ! 'final' values of number and mass weighted mean fallspeed for snow (m/s)
        2829             :
        2830           0 :             ums(k) = min(asn(i,k)*cons8/(6._r8*lams(k)**bs),1.2_r8*rhof(i,k))
        2831           0 :             uns(k) = min(asn(i,k)*cons7/lams(k)**bs,1.2_r8*rhof(i,k))
        2832             :
        2833             :          else
        2834           0 :             lams(k) = 0._r8
        2835           0 :             n0s(k) = 0._r8
        2836           0 :             ums(k) = 0._r8
        2837           0 :             uns(k) = 0._r8
        2838             :          end if
        2839             :
        2840             :          !c........................................................................
        2841             :          ! sum over sub-step for average process rates
        2842             :
        2843             :          ! convert rain/snow q and N for output to history, note,
        2844             :          ! output is for gridbox average
        2845             :
        2846           0 :          qrout(i,k)=qrout(i,k)+qric(i,k)*cldmax(i,k)
        2847           0 :          qsout(i,k)=qsout(i,k)+qniic(i,k)*cldmax(i,k)
        2848           0 :          nrout(i,k)=nrout(i,k)+nric(i,k)*rho(i,k)*cldmax(i,k)
        2849           0 :          nsout(i,k)=nsout(i,k)+nsic(i,k)*rho(i,k)*cldmax(i,k)
        2850             :
        2851           0 :          tlat1(i,k)=tlat1(i,k)+tlat(i,k)
        2852           0 :          qvlat1(i,k)=qvlat1(i,k)+qvlat(i,k)
        2853           0 :          qctend1(i,k)=qctend1(i,k)+qctend(i,k)
        2854           0 :          qitend1(i,k)=qitend1(i,k)+qitend(i,k)
        2855           0 :          nctend1(i,k)=nctend1(i,k)+nctend(i,k)
        2856           0 :          nitend1(i,k)=nitend1(i,k)+nitend(i,k)
        2857             :
        2858           0 :          t(i,k)=t(i,k)+tlat(i,k)*deltat/cpp
        2859           0 :          q(i,k)=q(i,k)+qvlat(i,k)*deltat
        2860           0 :          qc(i,k)=qc(i,k)+qctend(i,k)*deltat
        2861           0 :          qi(i,k)=qi(i,k)+qitend(i,k)*deltat
        2862           0 :          nc(i,k)=nc(i,k)+nctend(i,k)*deltat
        2863           0 :          ni(i,k)=ni(i,k)+nitend(i,k)*deltat
        2864             :
        2865           0 :          rainrt1(i,k)=rainrt1(i,k)+rainrt(i,k)
        2866             :
        2867             :          !divide rain radius over substeps for average
        2868           0 :          if (arcld(i,k) .gt. 0._r8) then
        2869           0 :             rercld(i,k)=rercld(i,k)/arcld(i,k)
        2870             :          end if
        2871             :
        2872             :          !! add to summing sub-stepping variable
        2873           0 :          rflx1(i,k+1)=rflx1(i,k+1)+rflx(i,k+1)
        2874           0 :          sflx1(i,k+1)=sflx1(i,k+1)+sflx(i,k+1)
        2875             :
        2876             :          !c........................................................................
        2877             :
        2878             :       end do ! k loop
        2879             :
        2880           0 :       prect1(i)=prect1(i)+prect(i)
        2881           0 :       preci1(i)=preci1(i)+preci(i)
        2882             :
        2883             :    end do ! it loop, sub-step
        2884             :
        2885           0 :    do k = top_lev, pver
        2886           0 :       rate1ord_cw2pr_st(i,k) = qcsinksum_rate1ord(k)/max(qcsum_rate1ord(k),1.0e-30_r8)
        2887             :    end do
        2888             :
        2889           0 : 300 continue  ! continue if no cloud water
        2890             : end do ! i loop
        2891             :
        2892             : ! convert dt from sub-step back to full time step
        2893           0 : deltat=deltat*real(iter)
        2894             :
        2895             : !c.............................................................................
        2896             :
        2897           0 : do i=1,ncol
        2898             :
        2899             :    ! skip all calculations if no cloud water
        2900           0 :    if (ltrue(i).eq.0) then
        2901             :
        2902           0 :       do k=1,top_lev-1
        2903             :          ! assign zero values for effective radius above 1 mbar
        2904           0 :          effc(i,k)=0._r8
        2905           0 :          effi(i,k)=0._r8
        2906           0 :          effc_fn(i,k)=0._r8
        2907           0 :          lamcrad(i,k)=0._r8
        2908           0 :          pgamrad(i,k)=0._r8
        2909           0 :          deffi(i,k)=0._r8
        2910             :       end do
        2911             :
        2912           0 :       do k=top_lev,pver
        2913             :          ! assign default values for effective radius
        2914           0 :          effc(i,k)=10._r8
        2915           0 :          effi(i,k)=25._r8
        2916           0 :          effc_fn(i,k)=10._r8
        2917           0 :          lamcrad(i,k)=0._r8
        2918           0 :          pgamrad(i,k)=0._r8
        2919           0 :          deffi(i,k)=0._r8
        2920             :       end do
        2921             :       goto 500
        2922             :    end if
        2923             :
        2924             :    ! initialize nstep for sedimentation sub-steps
        2925           0 :    nstep = 1
        2926             :
        2927             :    ! divide precip rate by number of sub-steps to get average over time step
        2928             :
        2929           0 :    prect(i)=prect1(i)/real(iter)
        2930           0 :    preci(i)=preci1(i)/real(iter)
        2931             :
        2932           0 :    do k=top_lev,pver
        2933             :
        2934             :       ! assign variables back to start-of-timestep values before updating after sub-steps
        2935             :
        2936           0 :       t(i,k)=t1(i,k)
        2937           0 :       q(i,k)=q1(i,k)
        2938           0 :       qc(i,k)=qc1(i,k)
        2939           0 :       qi(i,k)=qi1(i,k)
        2940           0 :       nc(i,k)=nc1(i,k)
        2941           0 :       ni(i,k)=ni1(i,k)
        2942             :
        2943             :       ! divide microphysical tendencies by number of sub-steps to get average over time step
        2944             :
        2945           0 :       tlat(i,k)=tlat1(i,k)/real(iter)
        2946           0 :       qvlat(i,k)=qvlat1(i,k)/real(iter)
        2947           0 :       qctend(i,k)=qctend1(i,k)/real(iter)
        2948           0 :       qitend(i,k)=qitend1(i,k)/real(iter)
        2949           0 :       nctend(i,k)=nctend1(i,k)/real(iter)
        2950           0 :       nitend(i,k)=nitend1(i,k)/real(iter)
        2951             :
        2952           0 :       rainrt(i,k)=rainrt1(i,k)/real(iter)
        2953             :
        2954             :       ! divide by number of sub-steps to find final values
        2955           0 :       rflx(i,k+1)=rflx1(i,k+1)/real(iter)
        2956           0 :       sflx(i,k+1)=sflx1(i,k+1)/real(iter)
        2957             :
        2958             :       ! divide output precip q and N by number of sub-steps to get average over time step
        2959             :
        2960           0 :       qrout(i,k)=qrout(i,k)/real(iter)
        2961           0 :       qsout(i,k)=qsout(i,k)/real(iter)
        2962           0 :       nrout(i,k)=nrout(i,k)/real(iter)
        2963           0 :       nsout(i,k)=nsout(i,k)/real(iter)
        2964             :
        2965             :       ! divide trop_mozart variables by number of sub-steps to get average over time step
        2966             :
        2967           0 :       nevapr(i,k) = nevapr(i,k)/real(iter)
        2968           0 :       nevapr2(i,k) = nevapr2(i,k)/real(iter)
        2969           0 :       evapsnow(i,k) = evapsnow(i,k)/real(iter)
        2970           0 :       prain(i,k) = prain(i,k)/real(iter)
        2971           0 :       prodsnow(i,k) = prodsnow(i,k)/real(iter)
        2972           0 :       cmeout(i,k) = cmeout(i,k)/real(iter)
        2973             :
        2974           0 :       cmeiout(i,k) = cmeiout(i,k)/real(iter)
        2975           0 :       meltsdt(i,k) = meltsdt(i,k)/real(iter)
        2976           0 :       frzrdt (i,k) = frzrdt (i,k)/real(iter)
        2977             :
        2978             :
        2979             :       ! microphysics output
        2980           0 :       prao(i,k)=prao(i,k)/real(iter)
        2981           0 :       prco(i,k)=prco(i,k)/real(iter)
        2982           0 :       mnuccco(i,k)=mnuccco(i,k)/real(iter)
        2983           0 :       mnuccto(i,k)=mnuccto(i,k)/real(iter)
        2984           0 :       msacwio(i,k)=msacwio(i,k)/real(iter)
        2985           0 :       psacwso(i,k)=psacwso(i,k)/real(iter)
        2986           0 :       bergso(i,k)=bergso(i,k)/real(iter)
        2987           0 :       bergo(i,k)=bergo(i,k)/real(iter)
        2988           0 :       prcio(i,k)=prcio(i,k)/real(iter)
        2989           0 :       praio(i,k)=praio(i,k)/real(iter)
        2990             :
        2991           0 :       mnuccro(i,k)=mnuccro(i,k)/real(iter)
        2992           0 :       pracso (i,k)=pracso (i,k)/real(iter)
        2993             :
        2994           0 :       mnuccdo(i,k)=mnuccdo(i,k)/real(iter)
        2995             :
        2996             :       ! modify to include snow. in prain & evap (diagnostic here: for wet dep)
        2997           0 :       nevapr(i,k) = nevapr(i,k) + evapsnow(i,k)
        2998           0 :       prer_evap(i,k) = nevapr2(i,k)
        2999           0 :       prain(i,k) = prain(i,k) + prodsnow(i,k)
        3000             :
        3001             :       !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        3002             :       ! calculate sedimentation for cloud water and ice
        3003             :       !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        3004             :
        3005             :       ! update in-cloud cloud mixing ratio and number concentration
        3006             :       ! with microphysical tendencies to calculate sedimentation, assign to dummy vars
        3007             :       ! note: these are in-cloud values***, hence we divide by cloud fraction
        3008             :
        3009           0 :       dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k)
        3010           0 :       dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)/icldm(i,k)
        3011           0 :       dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k),0._r8)
        3012           0 :       dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)/icldm(i,k),0._r8)
        3013             :
        3014           0 :       if (nccons) then
        3015           0 :         dumnc(i,k) = ncnst/rho(i,k)
        3016             :       end if
        3017           0 :       if (nicons) then
        3018           0 :         dumni(i,k) = ninst/rho(i,k)
        3019             :       end if
        3020             :
        3021             :       ! obtain new slope parameter to avoid possible singularity
        3022             :
        3023           0 :       if (dumi(i,k).ge.qsmall) then
        3024             :          ! add upper limit to in-cloud number concentration to prevent numerical error
        3025           0 :          dumni(i,k)=min(dumni(i,k),dumi(i,k)*1.e20_r8)
        3026             :
        3027             :          lami(k) = (cons1*ci* &
        3028           0 :               dumni(i,k)/dumi(i,k))**(1._r8/di)
        3029           0 :          lami(k)=max(lami(k),lammini)
        3030           0 :          lami(k)=min(lami(k),lammaxi)
        3031             :       else
        3032           0 :          lami(k)=0._r8
        3033             :       end if
        3034             :
        3035           0 :       if (dumc(i,k).ge.qsmall) then
        3036             :          ! add upper limit to in-cloud number concentration to prevent numerical error
        3037           0 :          dumnc(i,k)=min(dumnc(i,k),dumc(i,k)*1.e20_r8)
        3038             :          ! add lower limit to in-cloud number concentration
        3039           0 :          dumnc(i,k)=max(dumnc(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3
        3040           0 :          pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
        3041           0 :          pgam(k)=1._r8/(pgam(k)**2)-1._r8
        3042           0 :          pgam(k)=max(pgam(k),2._r8)
        3043           0 :          pgam(k)=min(pgam(k),15._r8)
        3044             :
        3045             :          lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ &
        3046           0 :               (dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
        3047           0 :          lammin = (pgam(k)+1._r8)/50.e-6_r8
        3048           0 :          lammax = (pgam(k)+1._r8)/2.e-6_r8
        3049           0 :          lamc(k)=max(lamc(k),lammin)
        3050           0 :          lamc(k)=min(lamc(k),lammax)
        3051             :       else
        3052           0 :          lamc(k)=0._r8
        3053             :       end if
        3054             :
        3055             :       ! calculate number and mass weighted fall velocity for droplets
        3056             :       ! include effects of sub-grid distribution of cloud water
        3057             :
        3058             :
        3059           0 :       if (dumc(i,k).ge.qsmall) then
        3060           0 :          unc = acn(i,k)*gamma(1._r8+bc+pgam(k))/(lamc(k)**bc*gamma(pgam(k)+1._r8))
        3061           0 :          umc = acn(i,k)*gamma(4._r8+bc+pgam(k))/(lamc(k)**bc*gamma(pgam(k)+4._r8))
        3062             :          ! fallspeed for output
        3063           0 :          vtrmc(i,k)=umc
        3064             :       else
        3065             :          umc = 0._r8
        3066             :          unc = 0._r8
        3067             :       end if
        3068             :
        3069             :       ! calculate number and mass weighted fall velocity for cloud ice
        3070             :
        3071           0 :       if (dumi(i,k).ge.qsmall) then
        3072           0 :          uni =  ain(i,k)*cons16/lami(k)**bi
        3073           0 :          umi = ain(i,k)*cons17/(6._r8*lami(k)**bi)
        3074           0 :          uni=min(uni,1.2_r8*rhof(i,k))
        3075           0 :          umi=min(umi,1.2_r8*rhof(i,k))
        3076             :
        3077             :          ! fallspeed
        3078           0 :          vtrmi(i,k)=umi
        3079             :       else
        3080             :          umi = 0._r8
        3081             :          uni = 0._r8
        3082             :       end if
        3083             :
        3084           0 :       fi(k) = g*rho(i,k)*umi
        3085           0 :       fni(k) = g*rho(i,k)*uni
        3086           0 :       fc(k) = g*rho(i,k)*umc
        3087           0 :       fnc(k) = g*rho(i,k)*unc
        3088             :
        3089             :       ! calculate number of split time steps to ensure courant stability criteria
        3090             :       ! for sedimentation calculations
        3091             :
        3092           0 :       rgvm = max(fi(k),fc(k),fni(k),fnc(k))
        3093           0 :       nstep = max(int(rgvm*deltat/pdel(i,k)+1._r8),nstep)
        3094             :
        3095             :       ! redefine dummy variables - sedimentation is calculated over grid-scale
        3096             :       ! quantities to ensure conservation
        3097             :
        3098           0 :       dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)
        3099           0 :       dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)
        3100           0 :       dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat),0._r8)
        3101           0 :       dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat),0._r8)
        3102             :
        3103           0 :       if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8
        3104           0 :       if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8
        3105             :
        3106             :    end do       !!! vertical loop
        3107           0 :    do n = 1,nstep  !! loop over sub-time step to ensure stability
        3108             :
        3109           0 :       do k = top_lev,pver
        3110           0 :          if (do_cldice) then
        3111           0 :             falouti(k) = fi(k)*dumi(i,k)
        3112           0 :             faloutni(k) = fni(k)*dumni(i,k)
        3113             :          else
        3114           0 :             falouti(k)  = 0._r8
        3115           0 :             faloutni(k) = 0._r8
        3116             :          end if
        3117             :
        3118           0 :          faloutc(k) = fc(k)*dumc(i,k)
        3119           0 :          faloutnc(k) = fnc(k)*dumnc(i,k)
        3120             :       end do
        3121             :
        3122             :       ! top of model
        3123             :
        3124           0 :       k = top_lev
        3125           0 :       faltndi = falouti(k)/pdel(i,k)
        3126           0 :       faltndni = faloutni(k)/pdel(i,k)
        3127           0 :       faltndc = faloutc(k)/pdel(i,k)
        3128           0 :       faltndnc = faloutnc(k)/pdel(i,k)
        3129             :
        3130             :       ! add fallout terms to microphysical tendencies
        3131             :
        3132           0 :       qitend(i,k) = qitend(i,k)-faltndi/nstep
        3133           0 :       nitend(i,k) = nitend(i,k)-faltndni/nstep
        3134           0 :       qctend(i,k) = qctend(i,k)-faltndc/nstep
        3135           0 :       nctend(i,k) = nctend(i,k)-faltndnc/nstep
        3136             :
        3137             :       ! sedimentation tendencies for output
        3138           0 :       qcsedten(i,k)=qcsedten(i,k)-faltndc/nstep
        3139           0 :       qisedten(i,k)=qisedten(i,k)-faltndi/nstep
        3140             :
        3141           0 :       dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep
        3142           0 :       dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep
        3143           0 :       dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep
        3144           0 :       dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep
        3145             :
        3146           0 :       do k = top_lev+1,pver
        3147             :
        3148             :          ! for cloud liquid and ice, if cloud fraction increases with height
        3149             :          ! then add flux from above to both vapor and cloud water of current level
        3150             :          ! this means that flux entering clear portion of cell from above evaporates
        3151             :          ! instantly
        3152             :
        3153           0 :          dum=lcldm(i,k)/lcldm(i,k-1)
        3154           0 :          dum=min(dum,1._r8)
        3155           0 :          dum1=icldm(i,k)/icldm(i,k-1)
        3156           0 :          dum1=min(dum1,1._r8)
        3157             :
        3158           0 :          faltndqie=(falouti(k)-falouti(k-1))/pdel(i,k)
        3159           0 :          faltndi=(falouti(k)-dum1*falouti(k-1))/pdel(i,k)
        3160           0 :          faltndni=(faloutni(k)-dum1*faloutni(k-1))/pdel(i,k)
        3161           0 :          faltndqce=(faloutc(k)-faloutc(k-1))/pdel(i,k)
        3162           0 :          faltndc=(faloutc(k)-dum*faloutc(k-1))/pdel(i,k)
        3163           0 :          faltndnc=(faloutnc(k)-dum*faloutnc(k-1))/pdel(i,k)
        3164             :
        3165             :          ! add fallout terms to eulerian tendencies
        3166             :
        3167           0 :          qitend(i,k) = qitend(i,k)-faltndi/nstep
        3168           0 :          nitend(i,k) = nitend(i,k)-faltndni/nstep
        3169           0 :          qctend(i,k) = qctend(i,k)-faltndc/nstep
        3170           0 :          nctend(i,k) = nctend(i,k)-faltndnc/nstep
        3171             :
        3172             :          ! sedimentation tendencies for output
        3173           0 :          qcsedten(i,k)=qcsedten(i,k)-faltndc/nstep
        3174           0 :          qisedten(i,k)=qisedten(i,k)-faltndi/nstep
        3175             :
        3176             :          ! add terms to to evap/sub of cloud water
        3177             :
        3178           0 :          qvlat(i,k)=qvlat(i,k)-(faltndqie-faltndi)/nstep
        3179             :          ! for output
        3180           0 :          qisevap(i,k)=qisevap(i,k)-(faltndqie-faltndi)/nstep
        3181           0 :          qvlat(i,k)=qvlat(i,k)-(faltndqce-faltndc)/nstep
        3182             :          ! for output
        3183           0 :          qcsevap(i,k)=qcsevap(i,k)-(faltndqce-faltndc)/nstep
        3184             :
        3185           0 :          tlat(i,k)=tlat(i,k)+(faltndqie-faltndi)*xxls/nstep
        3186           0 :          tlat(i,k)=tlat(i,k)+(faltndqce-faltndc)*xxlv/nstep
        3187             :
        3188           0 :          dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep
        3189           0 :          dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep
        3190           0 :          dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep
        3191           0 :          dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep
        3192             :
        3193           0 :          Fni(K)=MAX(Fni(K)/pdel(i,K),Fni(K-1)/pdel(i,K-1))*pdel(i,K)
        3194           0 :          FI(K)=MAX(FI(K)/pdel(i,K),FI(K-1)/pdel(i,K-1))*pdel(i,K)
        3195           0 :          fnc(k)=max(fnc(k)/pdel(i,k),fnc(k-1)/pdel(i,k-1))*pdel(i,k)
        3196           0 :          Fc(K)=MAX(Fc(K)/pdel(i,K),Fc(K-1)/pdel(i,K-1))*pdel(i,K)
        3197             :
        3198             :       end do   !! k loop
        3199             :
        3200             :       ! units below are m/s
        3201             :       ! cloud water/ice sedimentation flux at surface
        3202             :       ! is added to precip flux at surface to get total precip (cloud + precip water)
        3203             :       ! rate
        3204             :
        3205           0 :       prect(i) = prect(i)+(faloutc(pver)+falouti(pver))/g/nstep/1000._r8
        3206           0 :       preci(i) = preci(i)+(falouti(pver))/g/nstep/1000._r8
        3207             :
        3208             :       ! Add fallout to Precip Flux: note unit change m/s *kg/m3 = kg/m2
        3209           0 :       do k = top_lev,pver
        3210           0 :          rflx(i,k+1)=rflx(i,k+1)+(faloutc(k))/g/nstep/1000._r8 * rhow
        3211           0 :          sflx(i,k+1)=sflx(i,k+1)+(falouti(k))/g/nstep/1000._r8 * rhow
        3212             :       end do
        3213             :
        3214             :    end do   !! nstep loop
        3215             :
        3216             :    ! end sedimentation
        3217             :    !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
        3218             :
        3219             :    ! get new update for variables that includes sedimentation tendency
        3220             :    ! note : here dum variables are grid-average, NOT in-cloud
        3221             :
        3222           0 :    do k=top_lev,pver
        3223             :
        3224           0 :       dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)
        3225           0 :       dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)
        3226           0 :       dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)
        3227           0 :       dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)
        3228             :
        3229           0 :       if (nccons) then
        3230           0 :         dumnc(i,k) = ncnst/rho(i,k)*lcldm(i,k)
        3231             :       end if
        3232           0 :       if (nicons) then
        3233           0 :         dumni(i,k) = ninst/rho(i,k)*icldm(i,k)
        3234             :       end if
        3235             :
        3236           0 :       if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8
        3237           0 :       if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8
        3238             :
        3239             :       ! calculate instantaneous processes (melting, homogeneous freezing)
        3240           0 :       if (do_cldice) then
        3241             :
        3242           0 :          if (t(i,k)+tlat(i,k)/cpp*deltat > tmelt) then
        3243           0 :             if (dumi(i,k) > 0._r8) then
        3244             :
        3245             :                ! limit so that melting does not push temperature below freezing
        3246           0 :                dum = -dumi(i,k)*xlf/cpp
        3247           0 :                if (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt.tmelt) then
        3248           0 :                   dum = (t(i,k)+tlat(i,k)/cpp*deltat-tmelt)*cpp/xlf
        3249           0 :                   dum = dum/dumi(i,k)*xlf/cpp
        3250           0 :                   dum = max(0._r8,dum)
        3251           0 :                   dum = min(1._r8,dum)
        3252             :                else
        3253             :                   dum = 1._r8
        3254             :                end if
        3255             :
        3256           0 :                qctend(i,k)=qctend(i,k)+dum*dumi(i,k)/deltat
        3257             :
        3258             :                ! for output
        3259           0 :                melto(i,k)=dum*dumi(i,k)/deltat
        3260             :
        3261             :                ! assume melting ice produces droplet
        3262             :                ! mean volume radius of 8 micron
        3263             :
        3264             :                nctend(i,k)=nctend(i,k)+3._r8*dum*dumi(i,k)/deltat/ &
        3265           0 :                     (4._r8*pi*5.12e-16_r8*rhow)
        3266             :
        3267           0 :                qitend(i,k)=((1._r8-dum)*dumi(i,k)-qi(i,k))/deltat
        3268           0 :                nitend(i,k)=((1._r8-dum)*dumni(i,k)-ni(i,k))/deltat
        3269           0 :                tlat(i,k)=tlat(i,k)-xlf*dum*dumi(i,k)/deltat
        3270             :             end if
        3271             :          end if
        3272             :
        3273             :          ! homogeneously freeze droplets at -40 C
        3274             :
        3275           0 :          if (t(i,k)+tlat(i,k)/cpp*deltat < 233.15_r8) then
        3276           0 :             if (dumc(i,k) > 0._r8) then
        3277             :
        3278             :                ! limit so that freezing does not push temperature above threshold
        3279           0 :                dum = dumc(i,k)*xlf/cpp
        3280           0 :                if (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.233.15_r8) then
        3281           0 :                   dum = -(t(i,k)+tlat(i,k)/cpp*deltat-233.15_r8)*cpp/xlf
        3282           0 :                   dum = dum/dumc(i,k)*xlf/cpp
        3283           0 :                   dum = max(0._r8,dum)
        3284           0 :                   dum = min(1._r8,dum)
        3285             :                else
        3286             :                   dum = 1._r8
        3287             :                end if
        3288             :
        3289           0 :                qitend(i,k)=qitend(i,k)+dum*dumc(i,k)/deltat
        3290             :                ! for output
        3291           0 :                homoo(i,k)=dum*dumc(i,k)/deltat
        3292             :
        3293             :                ! assume 25 micron mean volume radius of homogeneously frozen droplets
        3294             :                ! consistent with size of detrained ice in stratiform.F90
        3295             :                nitend(i,k)=nitend(i,k)+dum*3._r8*dumc(i,k)/(4._r8*3.14_r8*1.563e-14_r8* &
        3296           0 :                     500._r8)/deltat
        3297           0 :                qctend(i,k)=((1._r8-dum)*dumc(i,k)-qc(i,k))/deltat
        3298           0 :                nctend(i,k)=((1._r8-dum)*dumnc(i,k)-nc(i,k))/deltat
        3299           0 :                tlat(i,k)=tlat(i,k)+xlf*dum*dumc(i,k)/deltat
        3300             :             end if
        3301             :          end if
        3302             :
        3303             :          ! remove any excess over-saturation, which is possible due to non-linearity when adding
        3304             :          ! together all microphysical processes
        3305             :          ! follow code similar to old CAM scheme
        3306             :
        3307           0 :          qtmp=q(i,k)+qvlat(i,k)*deltat
        3308           0 :          ttmp=t(i,k)+tlat(i,k)/cpp*deltat
        3309             :
        3310           0 :          esn = svp_water(ttmp)  ! use rhw to allow ice supersaturation
        3311           0 :          qsn = svp_to_qsat(esn, p(i,k))
        3312             :
        3313           0 :          if (qtmp > qsn .and. qsn > 0) then
        3314             :             ! expression below is approximate since there may be ice deposition
        3315           0 :             dum = (qtmp-qsn)/(1._r8+cons27*qsn/(cpp*rv*ttmp**2))/deltat
        3316             :             ! add to output cme
        3317           0 :             cmeout(i,k) = cmeout(i,k)+dum
        3318             :             ! now add to tendencies, partition between liquid and ice based on temperature
        3319           0 :             if (ttmp > 268.15_r8) then
        3320             :                dum1=0.0_r8
        3321             :                ! now add to tendencies, partition between liquid and ice based on te
        3322           0 :             else if (ttmp < 238.15_r8) then
        3323             :                dum1=1.0_r8
        3324             :             else
        3325           0 :                dum1=(268.15_r8-ttmp)/30._r8
        3326             :             end if
        3327             :
        3328             :             dum = (qtmp-qsn)/(1._r8+(xxls*dum1+xxlv*(1._r8-dum1))**2 &
        3329           0 :                  *qsn/(cpp*rv*ttmp**2))/deltat
        3330           0 :             qctend(i,k)=qctend(i,k)+dum*(1._r8-dum1)
        3331             :             ! for output
        3332           0 :             qcreso(i,k)=dum*(1._r8-dum1)
        3333           0 :             qitend(i,k)=qitend(i,k)+dum*dum1
        3334           0 :             qireso(i,k)=dum*dum1
        3335           0 :             qvlat(i,k)=qvlat(i,k)-dum
        3336             :             ! for output
        3337           0 :             qvres(i,k)=-dum
        3338           0 :             tlat(i,k)=tlat(i,k)+dum*(1._r8-dum1)*xxlv+dum*dum1*xxls
        3339             :          end if
        3340             :       end if
        3341             :
        3342             :       !...............................................................................
        3343             :       ! calculate effective radius for pass to radiation code
        3344             :       ! if no cloud water, default value is 10 micron for droplets,
        3345             :       ! 25 micron for cloud ice
        3346             :
        3347             :       ! update cloud variables after instantaneous processes to get effective radius
        3348             :       ! variables are in-cloud to calculate size dist parameters
        3349             :
        3350           0 :       dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)/lcldm(i,k)
        3351           0 :       dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)/icldm(i,k)
        3352           0 :       dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)/lcldm(i,k)
        3353           0 :       dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)/icldm(i,k)
        3354             :
        3355           0 :       if (nccons) then
        3356           0 :         dumnc(i,k) = ncnst/rho(i,k)
        3357             :       end if
        3358           0 :       if (nicons) then
        3359           0 :         dumni(i,k) = ninst/rho(i,k)
        3360             :       end if
        3361             :
        3362             :       ! limit in-cloud mixing ratio to reasonable value of 5 g kg-1
        3363             :
        3364           0 :       dumc(i,k)=min(dumc(i,k),5.e-3_r8)
        3365           0 :       dumi(i,k)=min(dumi(i,k),5.e-3_r8)
        3366             :
        3367             :       !...................
        3368             :       ! cloud ice effective radius
        3369             :
        3370           0 :       if (dumi(i,k).ge.qsmall) then
        3371             :
        3372           0 :          if (nicons) then
        3373             :            ! make sure ni is consistent with the constant N by adjusting
        3374             :            ! tendency, need to multiply by cloud fraction
        3375             :            ! note that nitend may be further adjusted below if mean crystal
        3376             :            ! size is out of bounds
        3377           0 :            nitend(i,k) = (ninst/rho(i,k)*icldm(i,k) - ni(i,k))/deltat
        3378             :          end if
        3379             :
        3380             :          ! add upper limit to in-cloud number concentration to prevent numerical error
        3381           0 :          dumni(i,k)=min(dumni(i,k),dumi(i,k)*1.e20_r8)
        3382           0 :          lami(k) = (cons1*ci*dumni(i,k)/dumi(i,k))**(1._r8/di)
        3383             :
        3384           0 :          if (lami(k).lt.lammini) then
        3385           0 :             lami(k) = lammini
        3386           0 :             n0i(k) = lami(k)**(di+1._r8)*dumi(i,k)/(ci*cons1)
        3387           0 :             niic(i,k) = n0i(k)/lami(k)
        3388             :             ! adjust number conc if needed to keep mean size in reasonable range
        3389           0 :             if (do_cldice) nitend(i,k)=(niic(i,k)*icldm(i,k)-ni(i,k))/deltat
        3390             :
        3391           0 :          else if (lami(k).gt.lammaxi) then
        3392           0 :             lami(k) = lammaxi
        3393           0 :             n0i(k) = lami(k)**(di+1._r8)*dumi(i,k)/(ci*cons1)
        3394           0 :             niic(i,k) = n0i(k)/lami(k)
        3395             :             ! adjust number conc if needed to keep mean size in reasonable range
        3396           0 :             if (do_cldice) nitend(i,k)=(niic(i,k)*icldm(i,k)-ni(i,k))/deltat
        3397             :          end if
        3398           0 :          effi(i,k) = 1.5_r8/lami(k)*1.e6_r8
        3399             :
        3400             :       else
        3401           0 :          effi(i,k) = 25._r8
        3402             :       end if
        3403             :
        3404             :       ! NOTE: If CARMA is doing the ice microphysics, then the ice effective
        3405             :       ! radius has already been determined from the size distribution.
        3406           0 :       if (.not. do_cldice) then
        3407           0 :          effi(i,k) = re_ice(i,k) * 1e6_r8      ! m -> um
        3408             :       end if
        3409             :
        3410             :       !...................
        3411             :       ! cloud droplet effective radius
        3412             :
        3413           0 :       if (dumc(i,k).ge.qsmall) then
        3414             :
        3415           0 :          if (nccons) then
        3416             :            ! make sure nc is consistent with the constant N by adjusting
        3417             :            ! tendency, need to multiply by cloud fraction
        3418             :            ! note that nctend may be further adjusted below if mean droplet
        3419             :            ! size is out of bounds
        3420           0 :            nctend(i,k) = (ncnst/rho(i,k)*lcldm(i,k) - nc(i,k))/deltat
        3421             :          end if
        3422             :
        3423             :          ! add upper limit to in-cloud number concentration to prevent numerical error
        3424           0 :          dumnc(i,k)=min(dumnc(i,k),dumc(i,k)*1.e20_r8)
        3425             :
        3426             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        3427             :          ! set tendency to ensure minimum droplet concentration
        3428             :          ! after update by microphysics, except when lambda exceeds bounds on mean drop
        3429             :          ! size or if there is no cloud water
        3430           0 :          if (dumnc(i,k).lt.cdnl/rho(i,k)) then
        3431           0 :             nctend(i,k)=(cdnl/rho(i,k)*lcldm(i,k)-nc(i,k))/deltat
        3432             :          end if
        3433           0 :          dumnc(i,k)=max(dumnc(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3
        3434             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
        3435           0 :          pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
        3436           0 :          pgam(k)=1._r8/(pgam(k)**2)-1._r8
        3437           0 :          pgam(k)=max(pgam(k),2._r8)
        3438           0 :          pgam(k)=min(pgam(k),15._r8)
        3439             :
        3440             :          lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ &
        3441           0 :               (dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
        3442           0 :          lammin = (pgam(k)+1._r8)/50.e-6_r8
        3443             :          ! Multiply by omsm to fit within RRTMG's table.
        3444           0 :          lammax = (pgam(k)+1._r8)*omsm/2.e-6_r8
        3445           0 :          if (lamc(k).lt.lammin) then
        3446           0 :             lamc(k) = lammin
        3447             :             ncic(i,k) = 6._r8*lamc(k)**3*dumc(i,k)* &
        3448             :                  gamma(pgam(k)+1._r8)/ &
        3449           0 :                  (pi*rhow*gamma(pgam(k)+4._r8))
        3450             :             ! adjust number conc if needed to keep mean size in reasonable range
        3451           0 :             nctend(i,k)=(ncic(i,k)*lcldm(i,k)-nc(i,k))/deltat
        3452             :
        3453           0 :          else if (lamc(k).gt.lammax) then
        3454           0 :             lamc(k) = lammax
        3455             :             ncic(i,k) = 6._r8*lamc(k)**3*dumc(i,k)* &
        3456             :                  gamma(pgam(k)+1._r8)/ &
        3457           0 :                  (pi*rhow*gamma(pgam(k)+4._r8))
        3458             :             ! adjust number conc if needed to keep mean size in reasonable range
        3459           0 :             nctend(i,k)=(ncic(i,k)*lcldm(i,k)-nc(i,k))/deltat
        3460             :          end if
        3461             :
        3462             :          effc(i,k) = &
        3463             :               gamma(pgam(k)+4._r8)/ &
        3464           0 :               gamma(pgam(k)+3._r8)/lamc(k)/2._r8*1.e6_r8
        3465             :          !assign output fields for shape here
        3466           0 :          lamcrad(i,k)=lamc(k)
        3467           0 :          pgamrad(i,k)=pgam(k)
        3468             :
        3469             :       else
        3470           0 :          effc(i,k) = 10._r8
        3471           0 :          lamcrad(i,k)=0._r8
        3472           0 :          pgamrad(i,k)=0._r8
        3473             :       end if
        3474             :
        3475             :       ! ice effective diameter for david mitchell's optics
        3476           0 :       if (do_cldice) then
        3477           0 :          deffi(i,k)=effi(i,k)*rhoi/917._r8*2._r8
        3478             :       else
        3479           0 :          deffi(i,k)=effi(i,k) * 2._r8
        3480             :       end if
        3481             :
        3482             :
        3483             : !!! recalculate effective radius for constant number, in order to separate
        3484             :       ! first and second indirect effects
        3485             :       ! assume constant number of 10^8 kg-1
        3486             :
        3487           0 :       dumnc(i,k)=1.e8_r8
        3488             :
        3489           0 :       if (dumc(i,k).ge.qsmall) then
        3490           0 :          pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8
        3491           0 :          pgam(k)=1._r8/(pgam(k)**2)-1._r8
        3492           0 :          pgam(k)=max(pgam(k),2._r8)
        3493           0 :          pgam(k)=min(pgam(k),15._r8)
        3494             :
        3495             :          lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ &
        3496           0 :               (dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8)
        3497           0 :          lammin = (pgam(k)+1._r8)/50.e-6_r8
        3498           0 :          lammax = (pgam(k)+1._r8)/2.e-6_r8
        3499           0 :          if (lamc(k).lt.lammin) then
        3500           0 :             lamc(k) = lammin
        3501           0 :          else if (lamc(k).gt.lammax) then
        3502           0 :             lamc(k) = lammax
        3503             :          end if
        3504             :          effc_fn(i,k) = &
        3505             :               gamma(pgam(k)+4._r8)/ &
        3506           0 :               gamma(pgam(k)+3._r8)/lamc(k)/2._r8*1.e6_r8
        3507             :
        3508             :       else
        3509           0 :          effc_fn(i,k) = 10._r8
        3510             :       end if
        3511             :
        3512             :
        3513             : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1!
        3514             :
        3515             :    end do ! vertical k loop
        3516             :
        3517             : 500 continue
        3518             :
        3519           0 :    do k=top_lev,pver
        3520             :       ! if updated q (after microphysics) is zero, then ensure updated n is also zero
        3521             :
        3522           0 :       if (qc(i,k)+qctend(i,k)*deltat.lt.qsmall) nctend(i,k)=-nc(i,k)/deltat
        3523           0 :       if (do_cldice .and. qi(i,k)+qitend(i,k)*deltat.lt.qsmall) nitend(i,k)=-ni(i,k)/deltat
        3524             :    end do
        3525             :
        3526             : end do ! i loop
        3527             :
        3528             : ! add snow ouptut
        3529           0 : do i = 1,ncol
        3530           0 :    do k=top_lev,pver
        3531           0 :       if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then
        3532           0 :          dsout(i,k)=3._r8*rhosn/917._r8*(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8)
        3533             :       endif
        3534             :    end do
        3535             : end do
        3536             :
        3537             : !calculate effective radius of rain and snow in microns for COSP using Eq. 9 of COSP v1.3 manual
        3538           0 : do i = 1,ncol
        3539           0 :    do k=top_lev,pver
        3540             :       !! RAIN
        3541           0 :       if (qrout(i,k).gt.1.e-7_r8.and.nrout(i,k).gt.0._r8) then
        3542           0 :          reff_rain(i,k)=1.5_r8*(pi * rhow * nrout(i,k)/qrout(i,k))**(-1._r8/3._r8)*1.e6_r8
        3543             :       endif
        3544             :       !! SNOW
        3545           0 :       if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then
        3546           0 :          reff_snow(i,k)=1.5_r8*(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8)*1.e6_r8
        3547             :       end if
        3548             :    end do
        3549             : end do
        3550             :
        3551             : ! analytic radar reflectivity
        3552             : ! formulas from Matthew Shupe, NOAA/CERES
        3553             : ! *****note: radar reflectivity is local (in-precip average)
        3554             : ! units of mm^6/m^3
        3555             :
        3556           0 : do i = 1,ncol
        3557           0 :    do k=top_lev,pver
        3558           0 :       if (qc(i,k)+qctend(i,k)*deltat.ge.qsmall .and. nc(i,k)+nctend(i,k)*deltat.gt.10._r8) then
        3559             :          dum=((qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)*1000._r8)**2 &
        3560           0 :               /(0.109_r8*(nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)/1.e6_r8)*lcldm(i,k)/cldmax(i,k)
        3561             :       else
        3562             :          dum=0._r8
        3563             :       end if
        3564           0 :       if (qi(i,k)+qitend(i,k)*deltat.ge.qsmall) then
        3565           0 :          dum1=((qi(i,k)+qitend(i,k)*deltat)*rho(i,k)/icldm(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)*icldm(i,k)/cldmax(i,k)
        3566             :       else
        3567             :          dum1=0._r8
        3568             :       end if
        3569             :
        3570           0 :       if (qsout(i,k).ge.qsmall) then
        3571           0 :          dum1=dum1+(qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)
        3572             :       end if
        3573             :
        3574           0 :       refl(i,k)=dum+dum1
        3575             :
        3576             :       ! add rain rate, but for 37 GHz formulation instead of 94 GHz
        3577             :       ! formula approximated from data of Matrasov (2007)
        3578             :       ! rainrt is the rain rate in mm/hr
        3579             :       ! reflectivity (dum) is in DBz
        3580             :
        3581           0 :       if (rainrt(i,k).ge.0.001_r8) then
        3582           0 :          dum=log10(rainrt(i,k)**6._r8)+16._r8
        3583             :
        3584             :          ! convert from DBz to mm^6/m^3
        3585             :
        3586           0 :          dum = 10._r8**(dum/10._r8)
        3587             :       else
        3588             :          ! don't include rain rate in R calculation for values less than 0.001 mm/hr
        3589             :          dum=0._r8
        3590             :       end if
        3591             :
        3592             :       ! add to refl
        3593             :
        3594           0 :       refl(i,k)=refl(i,k)+dum
        3595             :
        3596             :       !output reflectivity in Z.
        3597           0 :       areflz(i,k)=refl(i,k)
        3598             :
        3599             :       ! convert back to DBz
        3600             :
        3601           0 :       if (refl(i,k).gt.minrefl) then
        3602           0 :          refl(i,k)=10._r8*log10(refl(i,k))
        3603             :       else
        3604           0 :          refl(i,k)=-9999._r8
        3605             :       end if
        3606             :
        3607             :       !set averaging flag
        3608           0 :       if (refl(i,k).gt.mindbz) then
        3609           0 :          arefl(i,k)=refl(i,k)
        3610           0 :          frefl(i,k)=1.0_r8
        3611             :       else
        3612           0 :          arefl(i,k)=0._r8
        3613           0 :          areflz(i,k)=0._r8
        3614           0 :          frefl(i,k)=0._r8
        3615             :       end if
        3616             :
        3617             :       ! bound cloudsat reflectivity
        3618             :
        3619           0 :       csrfl(i,k)=min(csmax,refl(i,k))
        3620             :
        3621             :       !set averaging flag
        3622           0 :       if (csrfl(i,k).gt.csmin) then
        3623           0 :          acsrfl(i,k)=refl(i,k)
        3624           0 :          fcsrfl(i,k)=1.0_r8
        3625             :       else
        3626           0 :          acsrfl(i,k)=0._r8
        3627           0 :          fcsrfl(i,k)=0._r8
        3628             :       end if
        3629             :
        3630             :    end do
        3631             : end do
        3632             :
        3633             :
        3634             : ! averaging for snow and rain number and diameter
        3635             :
        3636           0 : qrout2(:,:)=0._r8
        3637           0 : qsout2(:,:)=0._r8
        3638           0 : nrout2(:,:)=0._r8
        3639           0 : nsout2(:,:)=0._r8
        3640           0 : drout2(:,:)=0._r8
        3641           0 : dsout2(:,:)=0._r8
        3642           0 : freqs(:,:)=0._r8
        3643           0 : freqr(:,:)=0._r8
        3644           0 : do i = 1,ncol
        3645           0 :    do k=top_lev,pver
        3646           0 :       if (qrout(i,k).gt.1.e-7_r8.and.nrout(i,k).gt.0._r8) then
        3647           0 :          qrout2(i,k)=qrout(i,k)
        3648           0 :          nrout2(i,k)=nrout(i,k)
        3649           0 :          drout2(i,k)=(pi * rhow * nrout(i,k)/qrout(i,k))**(-1._r8/3._r8)
        3650           0 :          freqr(i,k)=1._r8
        3651             :       endif
        3652           0 :       if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then
        3653           0 :          qsout2(i,k)=qsout(i,k)
        3654           0 :          nsout2(i,k)=nsout(i,k)
        3655           0 :          dsout2(i,k)=(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8)
        3656           0 :          freqs(i,k)=1._r8
        3657             :       endif
        3658             :    end do
        3659             : end do
        3660             :
        3661             : ! output activated liquid and ice (convert from #/kg -> #/m3)
        3662           0 : do i = 1,ncol
        3663           0 :    do k=top_lev,pver
        3664           0 :       ncai(i,k)=dum2i(i,k)*rho(i,k)
        3665           0 :       ncal(i,k)=dum2l(i,k)*rho(i,k)
        3666             :    end do
        3667             : end do
        3668             :
        3669             :
        3670             : !redefine fice here....
        3671           0 : nfice(:,:)=0._r8
        3672           0 : do k=top_lev,pver
        3673           0 :    do i=1,ncol
        3674           0 :       dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)
        3675           0 :       dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)
        3676           0 :       dumfice=qsout(i,k) + qrout(i,k) + dumc(i,k) + dumi(i,k)
        3677             :
        3678           0 :       if (dumfice.gt.qsmall.and.(qsout(i,k)+dumi(i,k).gt.qsmall)) then
        3679           0 :          nfice(i,k)=(qsout(i,k) + dumi(i,k))/dumfice
        3680             :       endif
        3681             :
        3682           0 :       if (nfice(i,k).gt.1._r8) then
        3683           0 :          nfice(i,k)=1._r8
        3684             :       endif
        3685             :
        3686             :    enddo
        3687             : enddo
        3688             :
        3689             :
        3690           0 : end subroutine micro_mg_tend
        3691             :
        3692             : !========================================================================
        3693             : !UTILITIES
        3694             : !========================================================================
        3695             :
        3696           0 : pure subroutine micro_mg_get_cols(ncol, nlev, top_lev, qcn, qin, &
        3697             :      mgncol, mgcols)
        3698             :
        3699             :   ! Determines which columns microphysics should operate over by
        3700             :   ! checking for non-zero cloud water/ice.
        3701             :
        3702             :   integer, intent(in) :: ncol      ! Number of columns with meaningful data
        3703             :   integer, intent(in) :: nlev      ! Number of levels to use
        3704             :   integer, intent(in) :: top_lev   ! Top level for microphysics
        3705             :
        3706             :   real(r8), intent(in) :: qcn(:,:) ! cloud water mixing ratio (kg/kg)
        3707             :   real(r8), intent(in) :: qin(:,:) ! cloud ice mixing ratio (kg/kg)
        3708             :
        3709             :   integer, intent(out) :: mgncol   ! Number of columns MG will use
        3710             :   integer, allocatable, intent(out) :: mgcols(:) ! column indices
        3711             :
        3712             :   integer :: lev_offset  ! top_lev - 1 (defined here for consistency)
        3713           0 :   logical :: ltrue(ncol) ! store tests for each column
        3714             :
        3715             :   integer :: i, ii ! column indices
        3716             :
        3717           0 :   if (allocated(mgcols)) deallocate(mgcols)
        3718             :
        3719           0 :   lev_offset = top_lev - 1
        3720             :
        3721             :   ! Using "any" along dimension 2 collapses across levels, but
        3722             :   ! not columns, so we know if water is present at any level
        3723             :   ! in each column.
        3724             :
        3725           0 :   ltrue = any(qcn(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2)
        3726           0 :   ltrue = ltrue .or. any(qin(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2)
        3727             :
        3728             :   ! Scan for true values to get a usable list of indices.
        3729             :
        3730           0 :   mgncol = count(ltrue)
        3731           0 :   allocate(mgcols(mgncol))
        3732           0 :   i = 0
        3733           0 :   do ii = 1,ncol
        3734           0 :      if (ltrue(ii)) then
        3735           0 :         i = i + 1
        3736           0 :         mgcols(i) = ii
        3737             :      end if
        3738             :   end do
        3739             :
        3740           0 : end subroutine micro_mg_get_cols
        3741             :
        3742             : end module micro_mg1_0

Generated by: LCOV version 1.14