module micro_mg1_0 !--------------------------------------------------------------------------------- ! Purpose: ! MG microphysics ! ! Author: Andrew Gettelman, Hugh Morrison. ! Contributions from: Xiaohong Liu and Steve Ghan ! December 2005-May 2010 ! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008) ! Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010) ! for questions contact Hugh Morrison, Andrew Gettelman ! e-mail: morrison@ucar.edu, andrew@ucar.edu ! ! NOTE: Modified to allow other microphysics packages (e.g. CARMA) to do ice ! microphysics in cooperation with the MG liquid microphysics. This is ! controlled by the do_cldice variable. ! ! NOTE: If do_cldice is false, then MG microphysics should not update CLDICE ! or NUMICE; however, it is assumed that the other microphysics scheme will have ! updated CLDICE and NUMICE. The other microphysics should handle the following ! processes that would have been done by MG: ! - Detrainment (liquid and ice) ! - Homogeneous ice nucleation ! - Heterogeneous ice nucleation ! - Bergeron process ! - Melting of ice ! - Freezing of cloud drops ! - Autoconversion (ice -> snow) ! - Growth/Sublimation of ice ! - Sedimentation of ice !--------------------------------------------------------------------------------- ! modification for sub-columns, HM, (orig 8/11/10) ! This is done using the logical 'microp_uniform' set to .true. = uniform for subcolumns !--------------------------------------------------------------------------------- ! Procedures required: ! 1) An implementation of the gamma function (if not intrinsic). ! 2) saturation vapor pressure to specific humidity formula ! 3) svp over water ! 4) svp over ice #ifndef HAVE_GAMMA_INTRINSICS use shr_spfn_mod, only: gamma => shr_spfn_gamma #endif use wv_sat_methods, only: & svp_water => wv_sat_svp_water, & svp_ice => wv_sat_svp_ice, & svp_to_qsat => wv_sat_svp_to_qsat use phys_control, only: phys_getopts implicit none private save ! Note: The liu_in option has been removed, as there was a serious bug with this ! option being set to false. The code now behaves as if the default liu_in=.true. ! is always on. Addition/reinstatement of ice nucleation options will likely be ! done outside of this module. public :: & micro_mg_init, & micro_mg_get_cols, & micro_mg_tend integer, parameter :: r8 = selected_real_kind(12) ! 8 byte real real(r8) :: g !gravity real(r8) :: r !Dry air Gas constant real(r8) :: rv !water vapor gas contstant real(r8) :: cpp !specific heat of dry air real(r8) :: rhow !density of liquid water real(r8) :: tmelt ! Freezing point of water (K) real(r8) :: xxlv ! latent heat of vaporization real(r8) :: xlf !latent heat of freezing real(r8) :: xxls !latent heat of sublimation real(r8) :: rhosn ! bulk density snow real(r8) :: rhoi ! bulk density ice real(r8) :: ac,bc,as,bs,ai,bi,ar,br !fall speed parameters real(r8) :: ci,di !ice mass-diameter relation parameters real(r8) :: cs,ds !snow mass-diameter relation parameters real(r8) :: cr,dr !drop mass-diameter relation parameters real(r8) :: f1s,f2s !ventilation param for snow real(r8) :: Eii !collection efficiency aggregation of ice real(r8) :: Ecr !collection efficiency cloud droplets/rain real(r8) :: f1r,f2r !ventilation param for rain real(r8) :: DCS !autoconversion size threshold real(r8) :: qsmall !min mixing ratio real(r8) :: bimm,aimm !immersion freezing real(r8) :: rhosu !typical 850mn air density real(r8) :: mi0 ! new crystal mass real(r8) :: rin ! radius of contact nuclei real(r8) :: pi ! pi ! Additional constants to help speed up code real(r8) :: cons1 real(r8) :: cons4 real(r8) :: cons5 real(r8) :: cons6 real(r8) :: cons7 real(r8) :: cons8 real(r8) :: cons11 real(r8) :: cons13 real(r8) :: cons14 real(r8) :: cons16 real(r8) :: cons17 real(r8) :: cons22 real(r8) :: cons23 real(r8) :: cons24 real(r8) :: cons25 real(r8) :: cons27 real(r8) :: cons28 real(r8) :: lammini real(r8) :: lammaxi real(r8) :: lamminr real(r8) :: lammaxr real(r8) :: lammins real(r8) :: lammaxs ! parameters for snow/rain fraction for convective clouds real(r8) :: tmax_fsnow ! max temperature for transition to convective snow real(r8) :: tmin_fsnow ! min temperature for transition to convective snow !needed for findsp real(r8) :: tt0 ! Freezing temperature real(r8) :: csmin,csmax,minrefl,mindbz real(r8) :: rhmini ! Minimum rh for ice cloud fraction > 0. logical :: use_hetfrz_classnuc ! option to use heterogeneous freezing character(len=16) :: micro_mg_precip_frac_method ! type of precipitation fraction method real(r8) :: micro_mg_berg_eff_factor ! berg efficiency factor !=============================================================================== contains !=============================================================================== subroutine micro_mg_init( & kind, gravit, rair, rh2o, cpair, & rhoh2o, tmelt_in, latvap, latice, & rhmini_in, micro_mg_dcs, use_hetfrz_classnuc_in, & micro_mg_precip_frac_method_in, micro_mg_berg_eff_factor_in, errstring) !----------------------------------------------------------------------- ! ! Purpose: ! initialize constants for the morrison microphysics ! ! Author: Andrew Gettelman Dec 2005 ! !----------------------------------------------------------------------- integer, intent(in) :: kind ! Kind used for reals real(r8), intent(in) :: gravit real(r8), intent(in) :: rair real(r8), intent(in) :: rh2o real(r8), intent(in) :: cpair real(r8), intent(in) :: rhoh2o real(r8), intent(in) :: tmelt_in ! Freezing point of water (K) real(r8), intent(in) :: latvap real(r8), intent(in) :: latice real(r8), intent(in) :: rhmini_in ! Minimum rh for ice cloud fraction > 0. real(r8), intent(in) :: micro_mg_dcs logical, intent(in) :: use_hetfrz_classnuc_in character(len=16),intent(in) :: micro_mg_precip_frac_method_in ! type of precipitation fraction method real(r8), intent(in) :: micro_mg_berg_eff_factor_in ! berg efficiency factor character(128), intent(out) :: errstring ! Output status (non-blank for error return) integer k integer l,m, iaer real(r8) surften ! surface tension of water w/respect to air (N/m) real(r8) arg !----------------------------------------------------------------------- errstring = ' ' if( kind .ne. r8 ) then errstring = 'micro_mg_init: KIND of reals does not match' return end if !declarations for morrison codes (transforms variable names) g= gravit !gravity r= rair !Dry air Gas constant: note units(phys_constants are in J/K/kmol) rv= rh2o !water vapor gas contstant cpp = cpair !specific heat of dry air rhow = rhoh2o !density of liquid water tmelt = tmelt_in rhmini = rhmini_in micro_mg_precip_frac_method = micro_mg_precip_frac_method_in micro_mg_berg_eff_factor = micro_mg_berg_eff_factor_in ! latent heats xxlv = latvap ! latent heat vaporization xlf = latice ! latent heat freezing xxls = xxlv + xlf ! latent heat of sublimation ! flags use_hetfrz_classnuc = use_hetfrz_classnuc_in ! parameters for snow/rain fraction for convective clouds tmax_fsnow = tmelt tmin_fsnow = tmelt-5._r8 ! parameters below from Reisner et al. (1998) ! density parameters (kg/m3) rhosn = 250._r8 ! bulk density snow (++ ceh) rhoi = 500._r8 ! bulk density ice rhow = 1000._r8 ! bulk density liquid ! fall speed parameters, V = aD^b ! V is in m/s ! droplets ac = 3.e7_r8 bc = 2._r8 ! snow as = 11.72_r8 bs = 0.41_r8 ! cloud ice ai = 700._r8 bi = 1._r8 ! rain ar = 841.99667_r8 br = 0.8_r8 ! particle mass-diameter relationship ! currently we assume spherical particles for cloud ice/snow ! m = cD^d pi= 3.1415927_r8 ! cloud ice mass-diameter relationship ci = rhoi*pi/6._r8 di = 3._r8 ! snow mass-diameter relationship cs = rhosn*pi/6._r8 ds = 3._r8 ! drop mass-diameter relationship cr = rhow*pi/6._r8 dr = 3._r8 ! ventilation parameters for snow ! hall and prupacher f1s = 0.86_r8 f2s = 0.28_r8 ! collection efficiency, aggregation of cloud ice and snow Eii = 0.1_r8 ! collection efficiency, accretion of cloud water by rain Ecr = 1.0_r8 ! ventilation constants for rain f1r = 0.78_r8 f2r = 0.32_r8 ! autoconversion size threshold for cloud ice to snow (m) Dcs = micro_mg_dcs ! smallest mixing ratio considered in microphysics qsmall = 1.e-18_r8 ! immersion freezing parameters, bigg 1953 bimm = 100._r8 aimm = 0.66_r8 ! typical air density at 850 mb rhosu = 85000._r8/(rair * tmelt) ! mass of new crystal due to aerosol freezing and growth (kg) mi0 = 4._r8/3._r8*pi*rhoi*(10.e-6_r8)*(10.e-6_r8)*(10.e-6_r8) ! radius of contact nuclei aerosol (m) rin = 0.1e-6_r8 ! freezing temperature tt0=273.15_r8 pi=4._r8*atan(1.0_r8) !Range of cloudsat reflectivities (dBz) for analytic simulator csmin= -30._r8 csmax= 26._r8 mindbz = -99._r8 ! minrefl = 10._r8**(mindbz/10._r8) minrefl = 1.26e-10_r8 ! Define constants to help speed up code (limit calls to gamma function) cons1=gamma(1._r8+di) cons4=gamma(1._r8+br) cons5=gamma(4._r8+br) cons6=gamma(1._r8+ds) cons7=gamma(1._r8+bs) cons8=gamma(4._r8+bs) cons11=gamma(3._r8+bs) cons13=gamma(5._r8/2._r8+br/2._r8) cons14=gamma(5._r8/2._r8+bs/2._r8) cons16=gamma(1._r8+bi) cons17=gamma(4._r8+bi) cons22=(4._r8/3._r8*pi*rhow*(25.e-6_r8)**3) cons23=dcs**3 cons24=dcs**2 cons25=dcs**bs cons27=xxlv**2 cons28=xxls**2 lammaxi = 1._r8/10.e-6_r8 lammini = 1._r8/(2._r8*dcs) lammaxr = 1._r8/20.e-6_r8 lamminr = 1._r8/500.e-6_r8 lammaxs = 1._r8/10.e-6_r8 lammins = 1._r8/2000.e-6_r8 end subroutine micro_mg_init !=============================================================================== !microphysics routine for each timestep goes here... subroutine micro_mg_tend ( & microp_uniform, pcols, pver, ncol, top_lev, deltatin,& tn, qn, qc, qi, nc, & ni, p, pdel, cldn, liqcldf, & relvar, accre_enhan, & icecldf, rate1ord_cw2pr_st, naai, npccnin, & rndst, nacon, tlat, qvlat, qctend, & qitend, nctend, nitend, effc, effc_fn, & effi, prect, preci, nevapr, evapsnow, am_evp_st, & prain, prodsnow, cmeout, deffi, pgamrad, & lamcrad, qsout, dsout, rflx, sflx, & qrout, reff_rain, reff_snow, qcsevap, qisevap, & qvres, cmeiout, vtrmc, vtrmi, qcsedten, & qisedten, prao, prco, mnuccco, mnuccto, & msacwio, psacwso, bergso, bergo, melto, & homoo, qcreso, prcio, praio, qireso, & mnuccro, pracso, meltsdt, frzrdt, mnuccdo, & nrout, nsout, refl, arefl, areflz, & frefl, csrfl, acsrfl, fcsrfl, rercld, & ncai, ncal, qrout2, qsout2, nrout2, & nsout2, drout2, dsout2, freqs, freqr, & nfice, prer_evap, do_cldice, errstring, & tnd_qsnow, tnd_nsnow, re_ice, & frzimm, frzcnt, frzdep) ! input arguments logical, intent(in) :: microp_uniform ! True = configure uniform for sub-columns False = use w/o sub-columns (standard) integer, intent(in) :: pcols ! size of column (first) index integer, intent(in) :: pver ! number of layers in columns integer, intent(in) :: ncol ! number of columns integer, intent(in) :: top_lev ! top level microphys is applied real(r8), intent(in) :: deltatin ! time step (s) real(r8), intent(in) :: tn(pcols,pver) ! input temperature (K) real(r8), intent(in) :: qn(pcols,pver) ! input h20 vapor mixing ratio (kg/kg) real(r8), intent(in) :: relvar(pcols,pver) ! relative variance of cloud water (-) real(r8), intent(in) :: accre_enhan(pcols,pver) ! optional accretion enhancement factor (-) ! note: all input cloud variables are grid-averaged real(r8), intent(inout) :: qc(pcols,pver) ! cloud water mixing ratio (kg/kg) real(r8), intent(inout) :: qi(pcols,pver) ! cloud ice mixing ratio (kg/kg) real(r8), intent(inout) :: nc(pcols,pver) ! cloud water number conc (1/kg) real(r8), intent(inout) :: ni(pcols,pver) ! cloud ice number conc (1/kg) real(r8), intent(in) :: p(pcols,pver) ! air pressure (pa) real(r8), intent(in) :: pdel(pcols,pver) ! pressure difference across level (pa) real(r8), intent(in) :: cldn(pcols,pver) ! cloud fraction real(r8), intent(in) :: icecldf(pcols,pver) ! ice cloud fraction real(r8), intent(in) :: liqcldf(pcols,pver) ! liquid cloud fraction real(r8), intent(out) :: rate1ord_cw2pr_st(pcols,pver) ! 1st order rate for direct cw to precip conversion ! used for scavenging ! Inputs for aerosol activation real(r8), intent(in) :: naai(pcols,pver) ! ice nulceation number (from microp_aero_ts) real(r8), intent(in) :: npccnin(pcols,pver) ! ccn activated number tendency (from microp_aero_ts) real(r8), intent(in) :: rndst(pcols,pver,4) ! radius of 4 dust bins for contact freezing (from microp_aero_ts) real(r8), intent(in) :: nacon(pcols,pver,4) ! number in 4 dust bins for contact freezing (from microp_aero_ts) ! Used with CARMA cirrus microphysics ! (or similar external microphysics model) logical, intent(in) :: do_cldice ! Prognosing cldice ! output arguments real(r8), intent(out) :: tlat(pcols,pver) ! latent heating rate (W/kg) real(r8), intent(out) :: qvlat(pcols,pver) ! microphysical tendency qv (1/s) real(r8), intent(out) :: qctend(pcols,pver) ! microphysical tendency qc (1/s) real(r8), intent(out) :: qitend(pcols,pver) ! microphysical tendency qi (1/s) real(r8), intent(out) :: nctend(pcols,pver) ! microphysical tendency nc (1/(kg*s)) real(r8), intent(out) :: nitend(pcols,pver) ! microphysical tendency ni (1/(kg*s)) real(r8), intent(out) :: effc(pcols,pver) ! droplet effective radius (micron) real(r8), intent(out) :: effc_fn(pcols,pver) ! droplet effective radius, assuming nc = 1.e8 kg-1 real(r8), intent(out) :: effi(pcols,pver) ! cloud ice effective radius (micron) real(r8), intent(out) :: prect(pcols) ! surface precip rate (m/s) real(r8), intent(out) :: preci(pcols) ! cloud ice/snow precip rate (m/s) real(r8), intent(out) :: nevapr(pcols,pver) ! evaporation rate of rain + snow real(r8), intent(out) :: evapsnow(pcols,pver)! sublimation rate of snow real(r8), intent(out) :: am_evp_st(pcols,pver)! stratiform evaporation area real(r8), intent(out) :: prain(pcols,pver) ! production of rain + snow real(r8), intent(out) :: prodsnow(pcols,pver)! production of snow real(r8), intent(out) :: cmeout(pcols,pver) ! evap/sub of cloud real(r8), intent(out) :: deffi(pcols,pver) ! ice effective diameter for optics (radiation) real(r8), intent(out) :: pgamrad(pcols,pver) ! ice gamma parameter for optics (radiation) real(r8), intent(out) :: lamcrad(pcols,pver) ! slope of droplet distribution for optics (radiation) real(r8), intent(out) :: qsout(pcols,pver) ! snow mixing ratio (kg/kg) real(r8), intent(out) :: dsout(pcols,pver) ! snow diameter (m) real(r8), intent(out) :: rflx(pcols,pver+1) ! grid-box average rain flux (kg m^-2 s^-1) real(r8), intent(out) :: sflx(pcols,pver+1) ! grid-box average snow flux (kg m^-2 s^-1) real(r8), intent(out) :: qrout(pcols,pver) ! grid-box average rain mixing ratio (kg/kg) real(r8), intent(inout) :: reff_rain(pcols,pver) ! rain effective radius (micron) real(r8), intent(inout) :: reff_snow(pcols,pver) ! snow effective radius (micron) real(r8), intent(out) :: qcsevap(pcols,pver) ! cloud water evaporation due to sedimentation real(r8), intent(out) :: qisevap(pcols,pver) ! cloud ice sublimation due to sublimation real(r8), intent(out) :: qvres(pcols,pver) ! residual condensation term to ensure RH < 100% real(r8), intent(out) :: cmeiout(pcols,pver) ! grid-mean cloud ice sub/dep real(r8), intent(out) :: vtrmc(pcols,pver) ! mass-weighted cloud water fallspeed real(r8), intent(out) :: vtrmi(pcols,pver) ! mass-weighted cloud ice fallspeed real(r8), intent(out) :: qcsedten(pcols,pver) ! qc sedimentation tendency real(r8), intent(out) :: qisedten(pcols,pver) ! qi sedimentation tendency ! microphysical process rates for output (mixing ratio tendencies) real(r8), intent(out) :: prao(pcols,pver) ! accretion of cloud by rain real(r8), intent(out) :: prco(pcols,pver) ! autoconversion of cloud to rain real(r8), intent(out) :: mnuccco(pcols,pver) ! mixing rat tend due to immersion freezing real(r8), intent(out) :: mnuccto(pcols,pver) ! mixing ratio tend due to contact freezing real(r8), intent(out) :: msacwio(pcols,pver) ! mixing ratio tend due to H-M splintering real(r8), intent(out) :: psacwso(pcols,pver) ! collection of cloud water by snow real(r8), intent(out) :: bergso(pcols,pver) ! bergeron process on snow real(r8), intent(out) :: bergo(pcols,pver) ! bergeron process on cloud ice real(r8), intent(out) :: melto(pcols,pver) ! melting of cloud ice real(r8), intent(out) :: homoo(pcols,pver) ! homogeneos freezign cloud water real(r8), intent(out) :: qcreso(pcols,pver) ! residual cloud condensation due to removal of excess supersat real(r8), intent(out) :: prcio(pcols,pver) ! autoconversion of cloud ice to snow real(r8), intent(out) :: praio(pcols,pver) ! accretion of cloud ice by snow real(r8), intent(out) :: qireso(pcols,pver) ! residual ice deposition due to removal of excess supersat real(r8), intent(out) :: mnuccro(pcols,pver) ! mixing ratio tendency due to heterogeneous freezing of rain to snow (1/s) real(r8), intent(out) :: pracso (pcols,pver) ! mixing ratio tendency due to accretion of rain by snow (1/s) real(r8), intent(out) :: meltsdt(pcols,pver) ! latent heating rate due to melting of snow (W/kg) real(r8), intent(out) :: frzrdt (pcols,pver) ! latent heating rate due to homogeneous freezing of rain (W/kg) real(r8), intent(out) :: mnuccdo(pcols,pver) ! mass tendency from ice nucleation real(r8), intent(out) :: nrout(pcols,pver) ! rain number concentration (1/m3) real(r8), intent(out) :: nsout(pcols,pver) ! snow number concentration (1/m3) real(r8), intent(out) :: refl(pcols,pver) ! analytic radar reflectivity real(r8), intent(out) :: arefl(pcols,pver) !average reflectivity will zero points outside valid range real(r8), intent(out) :: areflz(pcols,pver) !average reflectivity in z. real(r8), intent(out) :: frefl(pcols,pver) real(r8), intent(out) :: csrfl(pcols,pver) !cloudsat reflectivity real(r8), intent(out) :: acsrfl(pcols,pver) !cloudsat average real(r8), intent(out) :: fcsrfl(pcols,pver) real(r8), intent(out) :: rercld(pcols,pver) ! effective radius calculation for rain + cloud real(r8), intent(out) :: ncai(pcols,pver) ! output number conc of ice nuclei available (1/m3) real(r8), intent(out) :: ncal(pcols,pver) ! output number conc of CCN (1/m3) real(r8), intent(out) :: qrout2(pcols,pver) real(r8), intent(out) :: qsout2(pcols,pver) real(r8), intent(out) :: nrout2(pcols,pver) real(r8), intent(out) :: nsout2(pcols,pver) real(r8), intent(out) :: drout2(pcols,pver) ! mean rain particle diameter (m) real(r8), intent(out) :: dsout2(pcols,pver) ! mean snow particle diameter (m) real(r8), intent(out) :: freqs(pcols,pver) real(r8), intent(out) :: freqr(pcols,pver) real(r8), intent(out) :: nfice(pcols,pver) real(r8), intent(out) :: prer_evap(pcols,pver) real(r8) :: nevapr2(pcols,pver) character(128), intent(out) :: errstring ! Output status (non-blank for error return) ! Tendencies calculated by external schemes that can replace MG's native ! process tendencies. ! Used with CARMA cirrus microphysics ! (or similar external microphysics model) real(r8), intent(in) :: tnd_qsnow(:,:) ! snow mass tendency (kg/kg/s) real(r8), intent(in) :: tnd_nsnow(:,:) ! snow number tendency (#/kg/s) real(r8), intent(in) :: re_ice(:,:) ! ice effective radius (m) ! From external ice nucleation. real(r8), intent(in) :: frzimm(:,:) ! Number tendency due to immersion freezing (1/cm3) real(r8), intent(in) :: frzcnt(:,:) ! Number tendency due to contact freezing (1/cm3) real(r8), intent(in) :: frzdep(:,:) ! Number tendency due to deposition nucleation (1/cm3) ! local workspace ! all units mks unless otherwise stated ! Additional constants to help speed up code real(r8) :: cons2 real(r8) :: cons3 real(r8) :: cons9 real(r8) :: cons10 real(r8) :: cons12 real(r8) :: cons15 real(r8) :: cons18 real(r8) :: cons19 real(r8) :: cons20 ! temporary variables for sub-stepping real(r8) :: t1(pcols,pver) real(r8) :: q1(pcols,pver) real(r8) :: qc1(pcols,pver) real(r8) :: qi1(pcols,pver) real(r8) :: nc1(pcols,pver) real(r8) :: ni1(pcols,pver) real(r8) :: tlat1(pcols,pver) real(r8) :: qvlat1(pcols,pver) real(r8) :: qctend1(pcols,pver) real(r8) :: qitend1(pcols,pver) real(r8) :: nctend1(pcols,pver) real(r8) :: nitend1(pcols,pver) real(r8) :: prect1(pcols) real(r8) :: preci1(pcols) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real(r8) :: deltat ! sub-time step (s) real(r8) :: omsm ! number near unity for round-off issues real(r8) :: dto2 ! dt/2 (s) real(r8) :: mincld ! minimum allowed cloud fraction real(r8) :: q(pcols,pver) ! water vapor mixing ratio (kg/kg) real(r8) :: t(pcols,pver) ! temperature (K) real(r8) :: rho(pcols,pver) ! air density (kg m-3) real(r8) :: dv(pcols,pver) ! diffusivity of water vapor in air real(r8) :: mu(pcols,pver) ! viscocity of air real(r8) :: sc(pcols,pver) ! schmidt number real(r8) :: kap(pcols,pver) ! thermal conductivity of air real(r8) :: rhof(pcols,pver) ! air density correction factor for fallspeed real(r8) :: cldmax(pcols,pver) ! precip fraction assuming maximum overlap real(r8) :: cldm(pcols,pver) ! cloud fraction real(r8) :: icldm(pcols,pver) ! ice cloud fraction real(r8) :: lcldm(pcols,pver) ! liq cloud fraction real(r8) :: icwc(pcols) ! in cloud water content (liquid+ice) real(r8) :: calpha(pcols) ! parameter for cond/evap (Zhang et al. 2003) real(r8) :: cbeta(pcols) ! parameter for cond/evap (Zhang et al. 2003) real(r8) :: cbetah(pcols) ! parameter for cond/evap (Zhang et al. 2003) real(r8) :: cgamma(pcols) ! parameter for cond/evap (Zhang et al. 2003) real(r8) :: cgamah(pcols) ! parameter for cond/evap (Zhang et al. 2003) real(r8) :: rcgama(pcols) ! parameter for cond/evap (Zhang et al. 2003) real(r8) :: cmec1(pcols) ! parameter for cond/evap (Zhang et al. 2003) real(r8) :: cmec2(pcols) ! parameter for cond/evap (Zhang et al. 2003) real(r8) :: cmec3(pcols) ! parameter for cond/evap (Zhang et al. 2003) real(r8) :: cmec4(pcols) ! parameter for cond/evap (Zhang et al. 2003) real(r8) :: qtmp ! dummy qv real(r8) :: dum ! temporary dummy variable real(r8) :: cme(pcols,pver) ! total (liquid+ice) cond/evap rate of cloud real(r8) :: cmei(pcols,pver) ! dep/sublimation rate of cloud ice real(r8) :: cwml(pcols,pver) ! cloud water mixing ratio real(r8) :: cwmi(pcols,pver) ! cloud ice mixing ratio real(r8) :: nnuccd(pver) ! ice nucleation rate from deposition/cond.-freezing real(r8) :: mnuccd(pver) ! mass tendency from ice nucleation real(r8) :: qcld ! total cloud water real(r8) :: lcldn(pcols,pver) ! fractional coverage of new liquid cloud real(r8) :: lcldo(pcols,pver) ! fractional coverage of old liquid cloud real(r8) :: nctend_mixnuc(pcols,pver) real(r8) :: arg ! argument of erfc ! for calculation of rate1ord_cw2pr_st real(r8) :: qcsinksum_rate1ord(pver) ! sum over iterations of cw to precip sink real(r8) :: qcsum_rate1ord(pver) ! sum over iterations of cloud water real(r8) :: alpha real(r8) :: dum1,dum2 !general dummy variables real(r8) :: npccn(pver) ! droplet activation rate real(r8) :: qcic(pcols,pver) ! in-cloud cloud liquid mixing ratio real(r8) :: qiic(pcols,pver) ! in-cloud cloud ice mixing ratio real(r8) :: qniic(pcols,pver) ! in-precip snow mixing ratio real(r8) :: qric(pcols,pver) ! in-precip rain mixing ratio real(r8) :: ncic(pcols,pver) ! in-cloud droplet number conc real(r8) :: niic(pcols,pver) ! in-cloud cloud ice number conc real(r8) :: nsic(pcols,pver) ! in-precip snow number conc real(r8) :: nric(pcols,pver) ! in-precip rain number conc real(r8) :: lami(pver) ! slope of cloud ice size distr real(r8) :: n0i(pver) ! intercept of cloud ice size distr real(r8) :: lamc(pver) ! slope of cloud liquid size distr real(r8) :: n0c(pver) ! intercept of cloud liquid size distr real(r8) :: lams(pver) ! slope of snow size distr real(r8) :: n0s(pver) ! intercept of snow size distr real(r8) :: lamr(pver) ! slope of rain size distr real(r8) :: n0r(pver) ! intercept of rain size distr real(r8) :: cdist1(pver) ! size distr parameter to calculate droplet freezing ! combined size of precip & cloud drops real(r8) :: arcld(pcols,pver) ! averaging control flag real(r8) :: Actmp !area cross section of drops real(r8) :: Artmp !area cross section of rain real(r8) :: pgam(pver) ! spectral width parameter of droplet size distr real(r8) :: lammax ! maximum allowed slope of size distr real(r8) :: lammin ! minimum allowed slope of size distr real(r8) :: nacnt ! number conc of contact ice nuclei real(r8) :: mnuccc(pver) ! mixing ratio tendency due to freezing of cloud water real(r8) :: nnuccc(pver) ! number conc tendency due to freezing of cloud water real(r8) :: mnucct(pver) ! mixing ratio tendency due to contact freezing of cloud water real(r8) :: nnucct(pver) ! number conc tendency due to contact freezing of cloud water real(r8) :: msacwi(pver) ! mixing ratio tendency due to HM ice multiplication real(r8) :: nsacwi(pver) ! number conc tendency due to HM ice multiplication real(r8) :: prc(pver) ! qc tendency due to autoconversion of cloud droplets real(r8) :: nprc(pver) ! number conc tendency due to autoconversion of cloud droplets real(r8) :: nprc1(pver) ! qr tendency due to autoconversion of cloud droplets real(r8) :: nsagg(pver) ! ns tendency due to self-aggregation of snow real(r8) :: dc0 ! mean size droplet size distr real(r8) :: ds0 ! mean size snow size distr (area weighted) real(r8) :: eci ! collection efficiency for riming of snow by droplets real(r8) :: psacws(pver) ! mixing rat tendency due to collection of droplets by snow real(r8) :: npsacws(pver) ! number conc tendency due to collection of droplets by snow real(r8) :: uni ! number-weighted cloud ice fallspeed real(r8) :: umi ! mass-weighted cloud ice fallspeed real(r8) :: uns(pver) ! number-weighted snow fallspeed real(r8) :: ums(pver) ! mass-weighted snow fallspeed real(r8) :: unr(pver) ! number-weighted rain fallspeed real(r8) :: umr(pver) ! mass-weighted rain fallspeed real(r8) :: unc ! number-weighted cloud droplet fallspeed real(r8) :: umc ! mass-weighted cloud droplet fallspeed real(r8) :: pracs(pver) ! mixing rat tendency due to collection of rain by snow real(r8) :: npracs(pver) ! number conc tendency due to collection of rain by snow real(r8) :: mnuccr(pver) ! mixing rat tendency due to freezing of rain real(r8) :: nnuccr(pver) ! number conc tendency due to freezing of rain real(r8) :: pra(pver) ! mixing rat tendnency due to accretion of droplets by rain real(r8) :: npra(pver) ! nc tendnency due to accretion of droplets by rain real(r8) :: nragg(pver) ! nr tendency due to self-collection of rain real(r8) :: prci(pver) ! mixing rat tendency due to autoconversion of cloud ice to snow real(r8) :: nprci(pver) ! number conc tendency due to autoconversion of cloud ice to snow real(r8) :: prai(pver) ! mixing rat tendency due to accretion of cloud ice by snow real(r8) :: nprai(pver) ! number conc tendency due to accretion of cloud ice by snow real(r8) :: qvs ! liquid saturation vapor mixing ratio real(r8) :: qvi ! ice saturation vapor mixing ratio real(r8) :: dqsdt ! change of sat vapor mixing ratio with temperature real(r8) :: dqsidt ! change of ice sat vapor mixing ratio with temperature real(r8) :: ab ! correction factor for rain evap to account for latent heat real(r8) :: qclr ! water vapor mixing ratio in clear air real(r8) :: abi ! correction factor for snow sublimation to account for latent heat real(r8) :: epss ! 1/ sat relaxation timescale for snow real(r8) :: epsr ! 1/ sat relaxation timescale for rain real(r8) :: pre(pver) ! rain mixing rat tendency due to evaporation real(r8) :: prds(pver) ! snow mixing rat tendency due to sublimation real(r8) :: qce ! dummy qc for conservation check real(r8) :: qie ! dummy qi for conservation check real(r8) :: nce ! dummy nc for conservation check real(r8) :: nie ! dummy ni for conservation check real(r8) :: ratio ! parameter for conservation check real(r8) :: dumc(pcols,pver) ! dummy in-cloud qc real(r8) :: dumnc(pcols,pver) ! dummy in-cloud nc real(r8) :: dumi(pcols,pver) ! dummy in-cloud qi real(r8) :: dumni(pcols,pver) ! dummy in-cloud ni real(r8) :: dums(pcols,pver) ! dummy in-cloud snow mixing rat real(r8) :: dumns(pcols,pver) ! dummy in-cloud snow number conc real(r8) :: dumr(pcols,pver) ! dummy in-cloud rain mixing rat real(r8) :: dumnr(pcols,pver) ! dummy in-cloud rain number conc ! below are parameters for cloud water and cloud ice sedimentation calculations real(r8) :: fr(pver) real(r8) :: fnr(pver) real(r8) :: fc(pver) real(r8) :: fnc(pver) real(r8) :: fi(pver) real(r8) :: fni(pver) real(r8) :: fs(pver) real(r8) :: fns(pver) real(r8) :: faloutr(pver) real(r8) :: faloutnr(pver) real(r8) :: faloutc(pver) real(r8) :: faloutnc(pver) real(r8) :: falouti(pver) real(r8) :: faloutni(pver) real(r8) :: falouts(pver) real(r8) :: faloutns(pver) real(r8) :: faltndr real(r8) :: faltndnr real(r8) :: faltndc real(r8) :: faltndnc real(r8) :: faltndi real(r8) :: faltndni real(r8) :: faltnds real(r8) :: faltndns real(r8) :: faltndqie real(r8) :: faltndqce !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! real(r8) :: relhum(pcols,pver) ! relative humidity real(r8) :: csigma(pcols) ! parameter for cond/evap of cloud water/ice real(r8) :: rgvm ! max fallspeed for all species real(r8) :: arn(pcols,pver) ! air density corrected rain fallspeed parameter real(r8) :: asn(pcols,pver) ! air density corrected snow fallspeed parameter real(r8) :: acn(pcols,pver) ! air density corrected cloud droplet fallspeed parameter real(r8) :: ain(pcols,pver) ! air density corrected cloud ice fallspeed parameter real(r8) :: nsubi(pver) ! evaporation of cloud ice number real(r8) :: nsubc(pver) ! evaporation of droplet number real(r8) :: nsubs(pver) ! evaporation of snow number real(r8) :: nsubr(pver) ! evaporation of rain number real(r8) :: mtime ! factor to account for droplet activation timescale real(r8) :: dz(pcols,pver) ! height difference across model vertical level !! add precip flux variables for sub-stepping real(r8) :: rflx1(pcols,pver+1) real(r8) :: sflx1(pcols,pver+1) ! returns from function/subroutine calls real(r8) :: tsp(pcols,pver) ! saturation temp (K) real(r8) :: qsp(pcols,pver) ! saturation mixing ratio (kg/kg) real(r8) :: qsphy(pcols,pver) ! saturation mixing ratio (kg/kg): hybrid rh real(r8) :: qs(pcols) ! liquid-ice weighted sat mixing rat (kg/kg) real(r8) :: es(pcols) ! liquid-ice weighted sat vapor press (pa) real(r8) :: esl(pcols,pver) ! liquid sat vapor pressure (pa) real(r8) :: esi(pcols,pver) ! ice sat vapor pressure (pa) ! sum of source/sink terms for diagnostic precip real(r8) :: qnitend(pcols,pver) ! snow mixing ratio source/sink term real(r8) :: nstend(pcols,pver) ! snow number concentration source/sink term real(r8) :: qrtend(pcols,pver) ! rain mixing ratio source/sink term real(r8) :: nrtend(pcols,pver) ! rain number concentration source/sink term real(r8) :: qrtot ! vertically-integrated rain mixing rat source/sink term real(r8) :: nrtot ! vertically-integrated rain number conc source/sink term real(r8) :: qstot ! vertically-integrated snow mixing rat source/sink term real(r8) :: nstot ! vertically-integrated snow number conc source/sink term ! new terms for Bergeron process real(r8) :: dumnnuc ! provisional ice nucleation rate (for calculating bergeron) real(r8) :: ninew ! provisional cloud ice number conc (for calculating bergeron) real(r8) :: qinew ! provisional cloud ice mixing ratio (for calculating bergeron) real(r8) :: qvl ! liquid sat mixing ratio real(r8) :: epsi ! 1/ sat relaxation timecale for cloud ice real(r8) :: prd ! provisional deposition rate of cloud ice at water sat real(r8) :: berg(pcols,pver) ! mixing rat tendency due to bergeron process for cloud ice real(r8) :: bergs(pver) ! mixing rat tendency due to bergeron process for snow !bergeron terms real(r8) :: bergtsf !bergeron timescale to remove all liquid real(r8) :: rhin !modified RH for vapor deposition ! diagnostic rain/snow for output to history ! values are in-precip (local) !!!! real(r8) :: drout(pcols,pver) ! rain diameter (m) !averageed rain/snow for history real(r8) :: dumfice !ice nucleation, droplet activation real(r8) :: dum2i(pcols,pver) ! number conc of ice nuclei available (1/kg) real(r8) :: dum2l(pcols,pver) ! number conc of CCN (1/kg) real(r8) :: ncmax real(r8) :: nimax real(r8) :: qcvar ! 1/relative variance of sub-grid qc ! loop array variables integer i,k,nstep,n, l integer ii,kk, m ! loop variables for sub-step solution integer iter,it,ltrue(pcols) ! used in contact freezing via dust particles real(r8) tcnt, viscosity, mfp real(r8) slip1, slip2, slip3, slip4 ! real(r8) dfaer1, dfaer2, dfaer3, dfaer4 ! real(r8) nacon1,nacon2,nacon3,nacon4 real(r8) ndfaer1, ndfaer2, ndfaer3, ndfaer4 real(r8) nslip1, nslip2, nslip3, nslip4 ! used in ice effective radius real(r8) bbi, cci, ak, iciwc, rvi ! used in Bergeron processe and water vapor deposition real(r8) Tk, deles, Aprpr, Bprpr, Cice, qi0, Crate, qidep ! mean cloud fraction over the time step real(r8) cldmw(pcols,pver) ! used in secondary ice production real(r8) ni_secp ! variabels to check for RH after rain evap real(r8) :: esn real(r8) :: qsn real(r8) :: ttmp real(r8) :: rainrt(pcols,pver) ! rain rate for reflectivity calculation real(r8) :: rainrt1(pcols,pver) real(r8) :: tmp real(r8) dmc,ssmc,dstrn ! variables for modal scheme. real(r8), parameter :: cdnl = 0.e6_r8 ! cloud droplet number limiter ! heterogeneous freezing real(r8) :: mnudep(pver) ! mixing ratio tendency due to deposition of water vapor real(r8) :: nnudep(pver) ! number conc tendency due to deposition of water vapor real(r8) :: con1 ! work cnstant real(r8) :: r3lx ! Mean volume radius (m) real(r8) :: mi0l real(r8) :: frztmp logical :: do_clubb_sgs !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! Return error message errstring = ' ' call phys_getopts(do_clubb_sgs_out = do_clubb_sgs) ! initialize output fields for number conc qand ice nucleation ncai(1:ncol,1:pver)=0._r8 ncal(1:ncol,1:pver)=0._r8 !Initialize rain size rercld(1:ncol,1:pver)=0._r8 arcld(1:ncol,1:pver)=0._r8 !initialize radiation output variables pgamrad(1:ncol,1:pver)=0._r8 ! liquid gamma parameter for optics (radiation) lamcrad(1:ncol,1:pver)=0._r8 ! slope of droplet distribution for optics (radiation) deffi (1:ncol,1:pver)=0._r8 ! slope of droplet distribution for optics (radiation) !initialize radiation output variables !initialize water vapor tendency term output qcsevap(1:ncol,1:pver)=0._r8 qisevap(1:ncol,1:pver)=0._r8 qvres (1:ncol,1:pver)=0._r8 cmeiout (1:ncol,1:pver)=0._r8 vtrmc (1:ncol,1:pver)=0._r8 vtrmi (1:ncol,1:pver)=0._r8 qcsedten (1:ncol,1:pver)=0._r8 qisedten (1:ncol,1:pver)=0._r8 prao(1:ncol,1:pver)=0._r8 prco(1:ncol,1:pver)=0._r8 mnuccco(1:ncol,1:pver)=0._r8 mnuccto(1:ncol,1:pver)=0._r8 msacwio(1:ncol,1:pver)=0._r8 psacwso(1:ncol,1:pver)=0._r8 bergso(1:ncol,1:pver)=0._r8 bergo(1:ncol,1:pver)=0._r8 melto(1:ncol,1:pver)=0._r8 homoo(1:ncol,1:pver)=0._r8 qcreso(1:ncol,1:pver)=0._r8 prcio(1:ncol,1:pver)=0._r8 praio(1:ncol,1:pver)=0._r8 qireso(1:ncol,1:pver)=0._r8 mnuccro(1:ncol,1:pver)=0._r8 pracso (1:ncol,1:pver)=0._r8 meltsdt(1:ncol,1:pver)=0._r8 frzrdt (1:ncol,1:pver)=0._r8 mnuccdo(1:ncol,1:pver)=0._r8 rflx(:,:)=0._r8 sflx(:,:)=0._r8 effc(:,:)=0._r8 effc_fn(:,:)=0._r8 effi(:,:)=0._r8 ! assign variable deltat for sub-stepping... deltat=deltatin ! parameters for scheme omsm=0.99999_r8 dto2=0.5_r8*deltat mincld=0.0001_r8 ! initialize multi-level fields q(1:ncol,1:pver)=qn(1:ncol,1:pver) t(1:ncol,1:pver)=tn(1:ncol,1:pver) ! initialize time-varying parameters do k=1,pver do i=1,ncol rho(i,k)=p(i,k)/(r*t(i,k)) dv(i,k) = 8.794E-5_r8*t(i,k)**1.81_r8/p(i,k) mu(i,k) = 1.496E-6_r8*t(i,k)**1.5_r8/(t(i,k)+120._r8) sc(i,k) = mu(i,k)/(rho(i,k)*dv(i,k)) kap(i,k) = 1.414e3_r8*1.496e-6_r8*t(i,k)**1.5_r8/(t(i,k)+120._r8) ! air density adjustment for fallspeed parameters ! includes air density correction factor to the ! power of 0.54 following Heymsfield and Bansemer 2007 rhof(i,k)=(rhosu/rho(i,k))**0.54_r8 arn(i,k)=ar*rhof(i,k) asn(i,k)=as*rhof(i,k) acn(i,k)=ac*rhof(i,k) ain(i,k)=ai*rhof(i,k) ! get dz from dp and hydrostatic approx ! keep dz positive (define as layer k-1 - layer k) dz(i,k)= pdel(i,k)/(rho(i,k)*g) end do end do ! initialization qc(1:ncol,1:top_lev-1) = 0._r8 qi(1:ncol,1:top_lev-1) = 0._r8 nc(1:ncol,1:top_lev-1) = 0._r8 ni(1:ncol,1:top_lev-1) = 0._r8 t1(1:ncol,1:pver) = t(1:ncol,1:pver) q1(1:ncol,1:pver) = q(1:ncol,1:pver) qc1(1:ncol,1:pver) = qc(1:ncol,1:pver) qi1(1:ncol,1:pver) = qi(1:ncol,1:pver) nc1(1:ncol,1:pver) = nc(1:ncol,1:pver) ni1(1:ncol,1:pver) = ni(1:ncol,1:pver) ! initialize tendencies to zero tlat1(1:ncol,1:pver)=0._r8 qvlat1(1:ncol,1:pver)=0._r8 qctend1(1:ncol,1:pver)=0._r8 qitend1(1:ncol,1:pver)=0._r8 nctend1(1:ncol,1:pver)=0._r8 nitend1(1:ncol,1:pver)=0._r8 ! initialize precip output qrout(1:ncol,1:pver)=0._r8 qsout(1:ncol,1:pver)=0._r8 nrout(1:ncol,1:pver)=0._r8 nsout(1:ncol,1:pver)=0._r8 dsout(1:ncol,1:pver)=0._r8 drout(1:ncol,1:pver)=0._r8 reff_rain(1:ncol,1:pver)=0._r8 reff_snow(1:ncol,1:pver)=0._r8 ! initialize variables for trop_mozart nevapr(1:ncol,1:pver) = 0._r8 nevapr2(1:ncol,1:pver) = 0._r8 evapsnow(1:ncol,1:pver) = 0._r8 prain(1:ncol,1:pver) = 0._r8 prodsnow(1:ncol,1:pver) = 0._r8 cmeout(1:ncol,1:pver) = 0._r8 am_evp_st(1:ncol,1:pver) = 0._r8 ! for refl calc rainrt1(1:ncol,1:pver) = 0._r8 ! initialize precip fraction and output tendencies cldmax(1:ncol,1:pver)=mincld !initialize aerosol number ! naer2(1:ncol,1:pver,:)=0._r8 dum2l(1:ncol,1:pver)=0._r8 dum2i(1:ncol,1:pver)=0._r8 ! initialize avg precip rate prect1(1:ncol)=0._r8 preci1(1:ncol)=0._r8 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !Get humidity and saturation vapor pressures do k=top_lev,pver do i=1,ncol ! find wet bulk temperature and saturation value for provisional t and q without ! condensation es(i) = svp_water(t(i,k)) qs(i) = svp_to_qsat(es(i), p(i,k)) ! Prevents negative values. if (qs(i) < 0.0_r8) then qs(i) = 1.0_r8 es(i) = p(i,k) end if esl(i,k)=svp_water(t(i,k)) esi(i,k)=svp_ice(t(i,k)) ! hm fix, make sure when above freezing that esi=esl, not active yet if (t(i,k).gt.tmelt)esi(i,k)=esl(i,k) relhum(i,k)=q(i,k)/qs(i) ! get cloud fraction, check for minimum cldm(i,k)=max(cldn(i,k),mincld) cldmw(i,k)=max(cldn(i,k),mincld) icldm(i,k)=max(icecldf(i,k),mincld) lcldm(i,k)=max(liqcldf(i,k),mincld) ! subcolumns, set cloud fraction variables to one ! if cloud water or ice is present, if not present ! set to mincld (mincld used instead of zero, to prevent ! possible division by zero errors if (microp_uniform) then cldm(i,k)=mincld cldmw(i,k)=mincld icldm(i,k)=mincld lcldm(i,k)=mincld if (qc(i,k).ge.qsmall) then lcldm(i,k)=1._r8 cldm(i,k)=1._r8 cldmw(i,k)=1._r8 end if if (qi(i,k).ge.qsmall) then cldm(i,k)=1._r8 icldm(i,k)=1._r8 end if end if ! sub-columns ! calculate nfice based on liquid and ice mmr (no rain and snow mmr available yet) nfice(i,k)=0._r8 dumfice=qc(i,k)+qi(i,k) if (dumfice.gt.qsmall .and. qi(i,k).gt.qsmall) then nfice(i,k)=qi(i,k)/dumfice endif if (do_cldice .and. (t(i,k).lt.tmelt - 5._r8)) then ! if aerosols interact with ice set number of activated ice nuclei dum2=naai(i,k) dumnnuc=(dum2-ni(i,k)/icldm(i,k))/deltat*icldm(i,k) dumnnuc=max(dumnnuc,0._r8) ! get provisional ni and qi after nucleation in order to calculate ! Bergeron process below ninew=ni(i,k)+dumnnuc*deltat qinew=qi(i,k)+dumnnuc*deltat*mi0 !T>268 else ninew=ni(i,k) qinew=qi(i,k) end if ! Initialize CME components cme(i,k) = 0._r8 cmei(i,k)=0._r8 !------------------------------------------------------------------- !Bergeron process ! make sure to initialize bergeron process to zero berg(i,k)=0._r8 prd = 0._r8 !condensation loop. ! get in-cloud qi and ni after nucleation if (icldm(i,k) .gt. 0._r8) then qiic(i,k)=qinew/icldm(i,k) niic(i,k)=ninew/icldm(i,k) else qiic(i,k)=0._r8 niic(i,k)=0._r8 endif !if T < 0 C then bergeron. if (do_cldice .and. (t(i,k).lt.273.15_r8)) then !if ice exists if (qi(i,k).gt.qsmall) then bergtsf = 0._r8 ! bergeron time scale (fraction of timestep) qvi = svp_to_qsat(esi(i,k), p(i,k)) qvl = svp_to_qsat(esl(i,k), p(i,k)) dqsidt = xxls*qvi/(rv*t(i,k)**2) abi = 1._r8+dqsidt*xxls/cpp ! get ice size distribution parameters if (qiic(i,k).ge.qsmall) then lami(k) = (cons1*ci* & niic(i,k)/qiic(i,k))**(1._r8/di) n0i(k) = niic(i,k)*lami(k) ! check for slope ! adjust vars if (lami(k).lt.lammini) then lami(k) = lammini n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1) else if (lami(k).gt.lammaxi) then lami(k) = lammaxi n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1) end if epsi = 2._r8*pi*n0i(k)*rho(i,k)*Dv(i,k)/(lami(k)*lami(k)) !if liquid exists if (qc(i,k).gt. qsmall) then !begin bergeron process ! do bergeron (vapor deposition with RHw=1) ! code to find berg (a rate) goes here ! calculate Bergeron process prd = epsi*(qvl-qvi)/abi else prd = 0._r8 end if ! multiply by cloud fraction prd = prd*min(icldm(i,k),lcldm(i,k)) ! transfer of existing cloud liquid to ice berg(i,k)=max(0._r8,prd) end if !end liquid exists bergeron if (berg(i,k).gt.0._r8) then bergtsf=max(0._r8,(qc(i,k)/berg(i,k))/deltat) if(bergtsf.lt.1._r8) berg(i,k) = max(0._r8,qc(i,k)/deltat) endif !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! if (bergtsf.lt.1._r8.or.icldm(i,k).gt.lcldm(i,k)) then if (qiic(i,k).ge.qsmall) then ! first case is for case when liquid water is present, but is completely depleted ! in time step, i.e., bergrsf > 0 but < 1 if (qc(i,k).ge.qsmall) then rhin = (1.0_r8 + relhum(i,k)) / 2._r8 if ((rhin*esl(i,k)/esi(i,k)) > 1._r8) then prd = epsi*(rhin*qvl-qvi)/abi ! multiply by cloud fraction assuming liquid/ice maximum overlap prd = prd*min(icldm(i,k),lcldm(i,k)) ! add to cmei cmei(i,k) = cmei(i,k) + (prd * (1._r8- bergtsf)) end if ! rhin end if ! qc > qsmall ! second case is for pure ice cloud, either no liquid, or icldm > lcldm if (qc(i,k).lt.qsmall.or.icldm(i,k).gt.lcldm(i,k)) then ! note: for case of no liquid, need to set liquid cloud fraction to zero ! store liquid cloud fraction in 'dum' if (qc(i,k).lt.qsmall) then dum=0._r8 else dum=lcldm(i,k) end if ! set RH to grid-mean value for pure ice cloud rhin = relhum(i,k) if ((rhin*esl(i,k)/esi(i,k)) > 1._r8) then prd = epsi*(rhin*qvl-qvi)/abi ! multiply by relevant cloud fraction for pure ice cloud ! assuming maximum overlap of liquid/ice prd = prd*max((icldm(i,k)-dum),0._r8) cmei(i,k) = cmei(i,k) + prd end if ! rhin end if ! qc or icldm > lcldm end if ! qiic end if ! bergtsf or icldm > lcldm ! if deposition, it should not reduce grid mean rhi below 1.0 if(cmei(i,k) > 0.0_r8 .and. (relhum(i,k)*esl(i,k)/esi(i,k)) > 1._r8 ) & cmei(i,k)=min(cmei(i,k),(q(i,k)-qs(i)*esi(i,k)/esl(i,k))/abi/deltat) end if !end ice exists loop !this ends temperature < 0. loop !------------------------------------------------------------------- end if ! !.............................................................. ! evaporation should not exceed available water if ((-berg(i,k)).lt.-qc(i,k)/deltat) berg(i,k) = max(qc(i,k)/deltat,0._r8) !sublimation process... if (do_cldice .and. ((relhum(i,k)*esl(i,k)/esi(i,k)).lt.1._r8 .and. qiic(i,k).ge.qsmall )) then qvi = svp_to_qsat(esi(i,k), p(i,k)) qvl = svp_to_qsat(esl(i,k), p(i,k)) dqsidt = xxls*qvi/(rv*t(i,k)**2) abi = 1._r8+dqsidt*xxls/cpp ! get ice size distribution parameters lami(k) = (cons1*ci* & niic(i,k)/qiic(i,k))**(1._r8/di) n0i(k) = niic(i,k)*lami(k) ! check for slope ! adjust vars if (lami(k).lt.lammini) then lami(k) = lammini n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1) else if (lami(k).gt.lammaxi) then lami(k) = lammaxi n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1) end if epsi = 2._r8*pi*n0i(k)*rho(i,k)*Dv(i,k)/(lami(k)*lami(k)) ! modify for ice fraction below prd = epsi*(relhum(i,k)*qvl-qvi)/abi * icldm(i,k) cmei(i,k)=min(prd,0._r8) endif ! sublimation should not exceed available ice if (cmei(i,k).lt.-qi(i,k)/deltat) cmei(i,k)=-qi(i,k)/deltat ! sublimation should not increase grid mean rhi above 1.0 if(cmei(i,k) < 0.0_r8 .and. (relhum(i,k)*esl(i,k)/esi(i,k)) < 1._r8 ) & cmei(i,k)=min(0._r8,max(cmei(i,k),(q(i,k)-qs(i)*esi(i,k)/esl(i,k))/abi/deltat)) ! limit cmei due for roundoff error cmei(i,k)=cmei(i,k)*omsm ! conditional for ice nucleation if (do_cldice .and. (t(i,k).lt.(tmelt - 5._r8))) then ! using Liu et al. (2007) ice nucleation with hooks into simulated aerosol ! ice nucleation rate (dum2) has already been calculated and read in (naai) dum2i(i,k)=naai(i,k) else dum2i(i,k)=0._r8 end if end do ! i loop end do ! k loop !! initialize sub-step precip flux variables do i=1,ncol !! flux is zero at top interface, so these should stay as 0. rflx1(i,1)=0._r8 sflx1(i,1)=0._r8 do k=top_lev,pver ! initialize normal and sub-step precip flux variables rflx1(i,k+1)=0._r8 sflx1(i,k+1)=0._r8 end do ! i loop end do ! k loop !! initialize final precip flux variables. do i=1,ncol !! flux is zero at top interface, so these should stay as 0. rflx(i,1)=0._r8 sflx(i,1)=0._r8 do k=top_lev,pver ! initialize normal and sub-step precip flux variables rflx(i,k+1)=0._r8 sflx(i,k+1)=0._r8 end do ! i loop end do ! k loop do i=1,ncol ltrue(i)=0 do k=top_lev,pver ! skip microphysical calculations if no cloud water if (qc(i,k).ge.qsmall.or.qi(i,k).ge.qsmall.or.cmei(i,k).ge.qsmall) ltrue(i)=1 end do end do ! assign number of sub-steps to iter ! use 2 sub-steps, following tests described in MG2008 iter = 2 ! get sub-step time step deltat=deltat/real(iter) ! since activation/nucleation processes are fast, need to take into account ! factor mtime = mixing timescale in cloud / model time step ! mixing time can be interpreted as cloud depth divided by sub-grid vertical velocity ! for now mixing timescale is assumed to be 1 timestep for modal aerosols, 20 min bulk ! note: mtime for bulk aerosols was set to: mtime=deltat/1200._r8 mtime=1._r8 rate1ord_cw2pr_st(:,:)=0._r8 ! rce 2010/05/01 !!!! skip calculations if no cloud water do i=1,ncol if (ltrue(i).eq.0) then tlat(i,1:pver)=0._r8 qvlat(i,1:pver)=0._r8 qctend(i,1:pver)=0._r8 qitend(i,1:pver)=0._r8 qnitend(i,1:pver)=0._r8 qrtend(i,1:pver)=0._r8 nctend(i,1:pver)=0._r8 nitend(i,1:pver)=0._r8 nrtend(i,1:pver)=0._r8 nstend(i,1:pver)=0._r8 prect(i)=0._r8 preci(i)=0._r8 rflx(i,1:pver+1)=0._r8 sflx(i,1:pver+1)=0._r8 qniic(i,1:pver)=0._r8 qric(i,1:pver)=0._r8 nsic(i,1:pver)=0._r8 nric(i,1:pver)=0._r8 rainrt(i,1:pver)=0._r8 goto 300 end if qcsinksum_rate1ord(1:pver)=0._r8 qcsum_rate1ord(1:pver)=0._r8 !!!!!!!!! begin sub-step!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !..................................................................................................... do it=1,iter ! initialize sub-step microphysical tendencies tlat(i,1:pver)=0._r8 qvlat(i,1:pver)=0._r8 qctend(i,1:pver)=0._r8 qitend(i,1:pver)=0._r8 qnitend(i,1:pver)=0._r8 qrtend(i,1:pver)=0._r8 nctend(i,1:pver)=0._r8 nitend(i,1:pver)=0._r8 nrtend(i,1:pver)=0._r8 nstend(i,1:pver)=0._r8 ! initialize diagnostic precipitation to zero qniic(i,1:pver)=0._r8 qric(i,1:pver)=0._r8 nsic(i,1:pver)=0._r8 nric(i,1:pver)=0._r8 rainrt(i,1:pver)=0._r8 ! begin new i,k loop, calculate new cldmax after adjustment to cldm above ! initialize vertically-integrated rain and snow tendencies qrtot = 0._r8 nrtot = 0._r8 qstot = 0._r8 nstot = 0._r8 ! initialize precip at surface prect(i)=0._r8 preci(i)=0._r8 ! initialize fluxes rflx(i,1:pver+1)=0._r8 sflx(i,1:pver+1)=0._r8 do k=top_lev,pver qcvar=relvar(i,k) cons2=gamma(qcvar+2.47_r8) cons3=gamma(qcvar) cons9=gamma(qcvar+2._r8) cons10=gamma(qcvar+1._r8) cons12=gamma(qcvar+1.15_r8) cons15=gamma(qcvar+bc/3._r8) cons18=qcvar**2.47_r8 cons19=qcvar**2 cons20=qcvar**1.15_r8 ! set cwml and cwmi to current qc and qi cwml(i,k)=qc(i,k) cwmi(i,k)=qi(i,k) ! initialize precip fallspeeds to zero ums(k)=0._r8 uns(k)=0._r8 umr(k)=0._r8 unr(k)=0._r8 ! calculate precip fraction based on maximum overlap assumption ! for sub-columns cldm has already been set to 1 if cloud ! water or ice is present, so cldmax will be correctly set below ! and nothing extra needs to be done here if (k.eq.top_lev) then cldmax(i,k)=cldm(i,k) else ! if rain or snow mix ratio is smaller than ! threshold, then set cldmax to cloud fraction at current level if (do_clubb_sgs) then if (qc(i,k).ge.qsmall.or.qi(i,k).ge.qsmall) then cldmax(i,k)=cldm(i,k) else cldmax(i,k)=cldmax(i,k-1) end if else if (qric(i,k-1).ge.qsmall.or.qniic(i,k-1).ge.qsmall) then cldmax(i,k)=max(cldmax(i,k-1),cldm(i,k)) else cldmax(i,k)=cldm(i,k) end if endif end if ! decrease in number concentration due to sublimation/evap ! divide by cloud fraction to get in-cloud decrease ! don't reduce Nc due to bergeron process if (cmei(i,k) < 0._r8 .and. qi(i,k) > qsmall .and. cldm(i,k) > mincld) then nsubi(k)=cmei(i,k)/qi(i,k)*ni(i,k)/cldm(i,k) else nsubi(k)=0._r8 end if nsubc(k)=0._r8 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! ice nucleation if activated nuclei exist at t<-5C AND rhmini + 5% if (do_cldice .and. dum2i(i,k).gt.0._r8.and.t(i,k).lt.(tmelt - 5._r8).and. & relhum(i,k)*esl(i,k)/esi(i,k).gt. rhmini+0.05_r8) then !if NCAI > 0. then set numice = ncai (as before) !note: this is gridbox averaged nnuccd(k)=(dum2i(i,k)-ni(i,k)/icldm(i,k))/deltat*icldm(i,k) nnuccd(k)=max(nnuccd(k),0._r8) nimax = dum2i(i,k)*icldm(i,k) !Calc mass of new particles using new crystal mass... !also this will be multiplied by mtime as nnuccd is... mnuccd(k) = nnuccd(k) * mi0 ! add mnuccd to cmei.... cmei(i,k)= cmei(i,k) + mnuccd(k) * mtime ! limit cmei qvi = svp_to_qsat(esi(i,k), p(i,k)) dqsidt = xxls*qvi/(rv*t(i,k)**2) abi = 1._r8+dqsidt*xxls/cpp cmei(i,k)=min(cmei(i,k),(q(i,k)-qvi)/abi/deltat) ! limit for roundoff error cmei(i,k)=cmei(i,k)*omsm else nnuccd(k)=0._r8 nimax = 0._r8 mnuccd(k) = 0._r8 end if !c............................................................................ !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! obtain in-cloud values of cloud water/ice mixing ratios and number concentrations ! for microphysical process calculations ! units are kg/kg for mixing ratio, 1/kg for number conc ! limit in-cloud values to 0.005 kg/kg qcic(i,k)=min(cwml(i,k)/lcldm(i,k),5.e-3_r8) qiic(i,k)=min(cwmi(i,k)/icldm(i,k),5.e-3_r8) ncic(i,k)=max(nc(i,k)/lcldm(i,k),0._r8) niic(i,k)=max(ni(i,k)/icldm(i,k),0._r8) if (qc(i,k) - berg(i,k)*deltat.lt.qsmall) then qcic(i,k)=0._r8 ncic(i,k)=0._r8 if (qc(i,k)-berg(i,k)*deltat.lt.0._r8) then berg(i,k)=qc(i,k)/deltat*omsm end if end if if (do_cldice .and. qi(i,k)+(cmei(i,k)+berg(i,k))*deltat.lt.qsmall) then qiic(i,k)=0._r8 niic(i,k)=0._r8 if (qi(i,k)+(cmei(i,k)+berg(i,k))*deltat.lt.0._r8) then cmei(i,k)=(-qi(i,k)/deltat-berg(i,k))*omsm end if end if ! add to cme output cmeout(i,k) = cmeout(i,k)+cmei(i,k) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! droplet activation ! calculate potential for droplet activation if cloud water is present ! formulation from Abdul-Razzak and Ghan (2000) and Abdul-Razzak et al. (1998), AR98 ! number tendency (npccnin) is read in from companion routine ! assume aerosols already activated are equal to number of existing droplets for simplicity ! multiply by cloud fraction to obtain grid-average tendency if (qcic(i,k).ge.qsmall) then npccn(k) = max(0._r8,npccnin(i,k)) dum2l(i,k)=(nc(i,k)+npccn(k)*deltat)/lcldm(i,k) dum2l(i,k)=max(dum2l(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3 ncmax = dum2l(i,k)*lcldm(i,k) else npccn(k)=0._r8 dum2l(i,k)=0._r8 ncmax = 0._r8 end if !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! get size distribution parameters based on in-cloud cloud water/ice ! these calculations also ensure consistency between number and mixing ratio !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !...................................................................... ! cloud ice if (qiic(i,k).ge.qsmall) then ! add upper limit to in-cloud number concentration to prevent numerical error niic(i,k)=min(niic(i,k),qiic(i,k)*1.e20_r8) lami(k) = (cons1*ci*niic(i,k)/qiic(i,k))**(1._r8/di) n0i(k) = niic(i,k)*lami(k) ! check for slope ! adjust vars if (lami(k).lt.lammini) then lami(k) = lammini n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1) niic(i,k) = n0i(k)/lami(k) else if (lami(k).gt.lammaxi) then lami(k) = lammaxi n0i(k) = lami(k)**(di+1._r8)*qiic(i,k)/(ci*cons1) niic(i,k) = n0i(k)/lami(k) end if else lami(k) = 0._r8 n0i(k) = 0._r8 end if if (qcic(i,k).ge.qsmall) then ! add upper limit to in-cloud number concentration to prevent numerical error ncic(i,k)=min(ncic(i,k),qcic(i,k)*1.e20_r8) ncic(i,k)=max(ncic(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm ! get pgam from fit to observations of martin et al. 1994 pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8 pgam(k)=1._r8/(pgam(k)**2)-1._r8 pgam(k)=max(pgam(k),2._r8) pgam(k)=min(pgam(k),15._r8) ! calculate lamc lamc(k) = (pi/6._r8*rhow*ncic(i,k)*gamma(pgam(k)+4._r8)/ & (qcic(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8) ! lammin, 50 micron diameter max mean size lammin = (pgam(k)+1._r8)/50.e-6_r8 lammax = (pgam(k)+1._r8)/2.e-6_r8 if (lamc(k).lt.lammin) then lamc(k) = lammin ncic(i,k) = 6._r8*lamc(k)**3*qcic(i,k)* & gamma(pgam(k)+1._r8)/(pi*rhow*gamma(pgam(k)+4._r8)) else if (lamc(k).gt.lammax) then lamc(k) = lammax ncic(i,k) = 6._r8*lamc(k)**3*qcic(i,k)* & gamma(pgam(k)+1._r8)/(pi*rhow*gamma(pgam(k)+4._r8)) end if ! parameter to calculate droplet freezing cdist1(k) = ncic(i,k)/gamma(pgam(k)+1._r8) else lamc(k) = 0._r8 cdist1(k) = 0._r8 end if !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! begin micropysical process calculations !................................................................. ! autoconversion of cloud liquid water to rain ! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc ! minimum qc of 1 x 10^-8 prevents floating point error if (qcic(i,k).ge.1.e-8_r8) then ! nprc is increase in rain number conc due to autoconversion ! nprc1 is decrease in cloud droplet conc due to autoconversion ! assume exponential sub-grid distribution of qc, resulting in additional ! factor related to qcvar below ! hm switch for sub-columns, don't include sub-grid qc if (microp_uniform) then prc(k) = 1350._r8*qcic(i,k)**2.47_r8* & (ncic(i,k)/1.e6_r8*rho(i,k))**(-1.79_r8) nprc(k) = prc(k)/(4._r8/3._r8*pi*rhow*(25.e-6_r8)**3) nprc1(k) = prc(k)/(qcic(i,k)/ncic(i,k)) else prc(k) = cons2/(cons3*cons18)*1350._r8*qcic(i,k)**2.47_r8* & (ncic(i,k)/1.e6_r8*rho(i,k))**(-1.79_r8) nprc(k) = prc(k)/cons22 nprc1(k) = prc(k)/(qcic(i,k)/ncic(i,k)) end if ! sub-column switch else prc(k)=0._r8 nprc(k)=0._r8 nprc1(k)=0._r8 end if ! add autoconversion to precip from above to get provisional rain mixing ratio ! and number concentration (qric and nric) ! 0.45 m/s is fallspeed of new rain drop (80 micron diameter) dum=0.45_r8 dum1=0.45_r8 if (k.eq.top_lev) then qric(i,k)=prc(k)*lcldm(i,k)*dz(i,k)/cldmax(i,k)/dum nric(i,k)=nprc(k)*lcldm(i,k)*dz(i,k)/cldmax(i,k)/dum else if (qric(i,k-1).ge.qsmall) then dum=umr(k-1) dum1=unr(k-1) end if ! no autoconversion of rain number if rain/snow falling from above ! this assumes that new drizzle drops formed by autoconversion are rapidly collected ! by the existing rain/snow particles from above if (qric(i,k-1).ge.1.e-9_r8.or.qniic(i,k-1).ge.1.e-9_r8) then nprc(k)=0._r8 end if qric(i,k) = (rho(i,k-1)*umr(k-1)*qric(i,k-1)*cldmax(i,k-1)+ & (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))))& /(dum*rho(i,k)*cldmax(i,k)) nric(i,k) = (rho(i,k-1)*unr(k-1)*nric(i,k-1)*cldmax(i,k-1)+ & (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))))& /(dum1*rho(i,k)*cldmax(i,k)) end if !....................................................................... ! Autoconversion of cloud ice to snow ! similar to Ferrier (1994) if (do_cldice) then if (t(i,k).le.273.15_r8.and.qiic(i,k).ge.qsmall) then ! note: assumes autoconversion timescale of 180 sec nprci(k) = n0i(k)/(lami(k)*180._r8)*exp(-lami(k)*dcs) prci(k) = pi*rhoi*n0i(k)/(6._r8*180._r8)* & (cons23/lami(k)+3._r8*cons24/lami(k)**2+ & 6._r8*dcs/lami(k)**3+6._r8/lami(k)**4)*exp(-lami(k)*dcs) else prci(k)=0._r8 nprci(k)=0._r8 end if else ! Add in the particles that we have already converted to snow, and ! don't do any further autoconversion of ice. prci(k) = tnd_qsnow(i, k) / cldm(i,k) nprci(k) = tnd_nsnow(i, k) / cldm(i,k) end if ! add autoconversion to flux from level above to get provisional snow mixing ratio ! and number concentration (qniic and nsic) dum=(asn(i,k)*cons25) dum1=(asn(i,k)*cons25) if (k.eq.top_lev) then qniic(i,k)=prci(k)*icldm(i,k)*dz(i,k)/cldmax(i,k)/dum nsic(i,k)=nprci(k)*icldm(i,k)*dz(i,k)/cldmax(i,k)/dum else if (qniic(i,k-1).ge.qsmall) then dum=ums(k-1) dum1=uns(k-1) end if qniic(i,k) = (rho(i,k-1)*ums(k-1)*qniic(i,k-1)*cldmax(i,k-1)+ & (rho(i,k)*dz(i,k)*((prci(k)+prai(k-1)+psacws(k-1)+bergs(k-1))*icldm(i,k)+(prds(k-1)+ & pracs(k-1)+mnuccr(k-1))*cldmax(i,k))))& /(dum*rho(i,k)*cldmax(i,k)) nsic(i,k) = (rho(i,k-1)*uns(k-1)*nsic(i,k-1)*cldmax(i,k-1)+ & (rho(i,k)*dz(i,k)*(nprci(k)*icldm(i,k)+(nsubs(k-1)+nsagg(k-1)+nnuccr(k-1))*cldmax(i,k))))& /(dum1*rho(i,k)*cldmax(i,k)) end if ! if precip mix ratio is zero so should number concentration if (qniic(i,k).lt.qsmall) then qniic(i,k)=0._r8 nsic(i,k)=0._r8 end if if (qric(i,k).lt.qsmall) then qric(i,k)=0._r8 nric(i,k)=0._r8 end if ! make sure number concentration is a positive number to avoid ! taking root of negative later nric(i,k)=max(nric(i,k),0._r8) nsic(i,k)=max(nsic(i,k),0._r8) !....................................................................... ! get size distribution parameters for precip !...................................................................... ! rain if (qric(i,k).ge.qsmall) then lamr(k) = (pi*rhow*nric(i,k)/qric(i,k))**(1._r8/3._r8) n0r(k) = nric(i,k)*lamr(k) ! check for slope ! adjust vars if (lamr(k).lt.lamminr) then lamr(k) = lamminr n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow) nric(i,k) = n0r(k)/lamr(k) else if (lamr(k).gt.lammaxr) then lamr(k) = lammaxr n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow) nric(i,k) = n0r(k)/lamr(k) end if ! provisional rain number and mass weighted mean fallspeed (m/s) unr(k) = min(arn(i,k)*cons4/lamr(k)**br,9.1_r8*rhof(i,k)) umr(k) = min(arn(i,k)*cons5/(6._r8*lamr(k)**br),9.1_r8*rhof(i,k)) else lamr(k) = 0._r8 n0r(k) = 0._r8 umr(k) = 0._r8 unr(k) = 0._r8 end if !...................................................................... ! snow if (qniic(i,k).ge.qsmall) then lams(k) = (cons6*cs*nsic(i,k)/qniic(i,k))**(1._r8/ds) n0s(k) = nsic(i,k)*lams(k) ! check for slope ! adjust vars if (lams(k).lt.lammins) then lams(k) = lammins n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6) nsic(i,k) = n0s(k)/lams(k) else if (lams(k).gt.lammaxs) then lams(k) = lammaxs n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6) nsic(i,k) = n0s(k)/lams(k) end if ! provisional snow number and mass weighted mean fallspeed (m/s) ums(k) = min(asn(i,k)*cons8/(6._r8*lams(k)**bs),1.2_r8*rhof(i,k)) uns(k) = min(asn(i,k)*cons7/lams(k)**bs,1.2_r8*rhof(i,k)) else lams(k) = 0._r8 n0s(k) = 0._r8 ums(k) = 0._r8 uns(k) = 0._r8 end if !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! heterogeneous freezing of cloud water if (.not. use_hetfrz_classnuc) then if (do_cldice .and. qcic(i,k).ge.qsmall .and. t(i,k).lt.269.15_r8) then ! immersion freezing (Bigg, 1953) ! subcolumns if (microp_uniform) then mnuccc(k) = & pi*pi/36._r8*rhow* & cdist1(k)*gamma(7._r8+pgam(k))* & bimm*(exp(aimm*(273.15_r8-t(i,k)))-1._r8)/ & lamc(k)**3/lamc(k)**3 nnuccc(k) = & pi/6._r8*cdist1(k)*gamma(pgam(k)+4._r8) & *bimm* & (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamc(k)**3 else mnuccc(k) = cons9/(cons3*cons19)* & pi*pi/36._r8*rhow* & cdist1(k)*gamma(7._r8+pgam(k))* & bimm*(exp(aimm*(273.15_r8-t(i,k)))-1._r8)/ & lamc(k)**3/lamc(k)**3 nnuccc(k) = cons10/(cons3*qcvar)* & pi/6._r8*cdist1(k)*gamma(pgam(k)+4._r8) & *bimm* & (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamc(k)**3 end if ! sub-columns ! contact freezing (-40= qsmall) then con1 = 1._r8/(1.333_r8*pi)**0.333_r8 r3lx = con1*(rho(i,k)*qcic(i,k)/(rhow*max(ncic(i,k)*rho(i,k), 1.0e6_r8)))**0.333_r8 ! in m r3lx = max(4.e-6_r8, r3lx) mi0l = 4._r8/3._r8*pi*rhow*r3lx**3_r8 nnuccc(k) = frzimm(i,k)*1.0e6_r8/rho(i,k) mnuccc(k) = nnuccc(k)*mi0l nnucct(k) = frzcnt(i,k)*1.0e6_r8/rho(i,k) mnucct(k) = nnucct(k)*mi0l nnudep(k) = frzdep(i,k)*1.0e6_r8/rho(i,k) mnudep(k) = nnudep(k)*mi0 else nnuccc(k) = 0._r8 mnuccc(k) = 0._r8 nnucct(k) = 0._r8 mnucct(k) = 0._r8 nnudep(k) = 0._r8 mnudep(k) = 0._r8 end if endif !....................................................................... ! snow self-aggregation from passarelli, 1978, used by reisner, 1998 ! this is hard-wired for bs = 0.4 for now ! ignore self-collection of cloud ice if (qniic(i,k).ge.qsmall .and. t(i,k).le.273.15_r8) then nsagg(k) = -1108._r8*asn(i,k)*Eii* & pi**((1._r8-bs)/3._r8)*rhosn**((-2._r8-bs)/3._r8)*rho(i,k)** & ((2._r8+bs)/3._r8)*qniic(i,k)**((2._r8+bs)/3._r8)* & (nsic(i,k)*rho(i,k))**((4._r8-bs)/3._r8)/ & (4._r8*720._r8*rho(i,k)) else nsagg(k)=0._r8 end if !....................................................................... ! accretion of cloud droplets onto snow/graupel ! here use continuous collection equation with ! simple gravitational collection kernel ! ignore collisions between droplets/cloud ice ! since minimum size ice particle for accretion is 50 - 150 micron ! ignore collision of snow with droplets above freezing if (qniic(i,k).ge.qsmall .and. t(i,k).le.tmelt .and. & qcic(i,k).ge.qsmall) then ! put in size dependent collection efficiency ! mean diameter of snow is area-weighted, since ! accretion is function of crystal geometric area ! collection efficiency is approximation based on stoke's law (Thompson et al. 2004) dc0 = (pgam(k)+1._r8)/lamc(k) ds0 = 1._r8/lams(k) dum = dc0*dc0*uns(k)*rhow/(9._r8*mu(i,k)*ds0) eci = dum*dum/((dum+0.4_r8)*(dum+0.4_r8)) eci = max(eci,0._r8) eci = min(eci,1._r8) ! no impact of sub-grid distribution of qc since psacws ! is linear in qc psacws(k) = pi/4._r8*asn(i,k)*qcic(i,k)*rho(i,k)* & n0s(k)*Eci*cons11/ & lams(k)**(bs+3._r8) npsacws(k) = pi/4._r8*asn(i,k)*ncic(i,k)*rho(i,k)* & n0s(k)*Eci*cons11/ & lams(k)**(bs+3._r8) else psacws(k)=0._r8 npsacws(k)=0._r8 end if ! add secondary ice production due to accretion of droplets by snow ! (Hallet-Mossop process) (from Cotton et al., 1986) if (.not. do_cldice) then ni_secp = 0.0_r8 nsacwi(k) = 0.0_r8 msacwi(k) = 0.0_r8 else if((t(i,k).lt.270.16_r8) .and. (t(i,k).ge.268.16_r8)) then ni_secp = 3.5e8_r8*(270.16_r8-t(i,k))/2.0_r8*psacws(k) nsacwi(k) = ni_secp msacwi(k) = min(ni_secp*mi0,psacws(k)) else if((t(i,k).lt.268.16_r8) .and. (t(i,k).ge.265.16_r8)) then ni_secp = 3.5e8_r8*(t(i,k)-265.16_r8)/3.0_r8*psacws(k) nsacwi(k) = ni_secp msacwi(k) = min(ni_secp*mi0,psacws(k)) else ni_secp = 0.0_r8 nsacwi(k) = 0.0_r8 msacwi(k) = 0.0_r8 endif psacws(k) = max(0.0_r8,psacws(k)-ni_secp*mi0) !....................................................................... ! accretion of rain water by snow ! formula from ikawa and saito, 1991, used by reisner et al., 1998 if (qric(i,k).ge.1.e-8_r8 .and. qniic(i,k).ge.1.e-8_r8 .and. & t(i,k).le.273.15_r8) then pracs(k) = pi*pi*ecr*(((1.2_r8*umr(k)-0.95_r8*ums(k))**2+ & 0.08_r8*ums(k)*umr(k))**0.5_r8*rhow*rho(i,k)* & n0r(k)*n0s(k)* & (5._r8/(lamr(k)**6*lams(k))+ & 2._r8/(lamr(k)**5*lams(k)**2)+ & 0.5_r8/(lamr(k)**4*lams(k)**3))) npracs(k) = pi/2._r8*rho(i,k)*ecr*(1.7_r8*(unr(k)-uns(k))**2+ & 0.3_r8*unr(k)*uns(k))**0.5_r8*n0r(k)*n0s(k)* & (1._r8/(lamr(k)**3*lams(k))+ & 1._r8/(lamr(k)**2*lams(k)**2)+ & 1._r8/(lamr(k)*lams(k)**3)) else pracs(k)=0._r8 npracs(k)=0._r8 end if !....................................................................... ! heterogeneous freezing of rain drops ! follows from Bigg (1953) if (t(i,k).lt.269.15_r8 .and. qric(i,k).ge.qsmall) then mnuccr(k) = 20._r8*pi*pi*rhow*nric(i,k)*bimm* & (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamr(k)**3 & /lamr(k)**3 nnuccr(k) = pi*nric(i,k)*bimm* & (exp(aimm*(273.15_r8-t(i,k)))-1._r8)/lamr(k)**3 else mnuccr(k)=0._r8 nnuccr(k)=0._r8 end if !....................................................................... ! accretion of cloud liquid water by rain ! formula from Khrouditnov and Kogan (2000) ! gravitational collection kernel, droplet fall speed neglected if (qric(i,k).ge.qsmall .and. qcic(i,k).ge.qsmall) then ! include sub-grid distribution of cloud water ! add sub-column switch if (microp_uniform) then pra(k) = 67._r8*(qcic(i,k)*qric(i,k))**1.15_r8 npra(k) = pra(k)/(qcic(i,k)/ncic(i,k)) else pra(k) = accre_enhan(i,k)*(cons12/(cons3*cons20)*67._r8*(qcic(i,k)*qric(i,k))**1.15_r8) npra(k) = pra(k)/(qcic(i,k)/ncic(i,k)) end if ! sub-column switch else pra(k)=0._r8 npra(k)=0._r8 end if !....................................................................... ! Self-collection of rain drops ! from Beheng(1994) if (qric(i,k).ge.qsmall) then nragg(k) = -8._r8*nric(i,k)*qric(i,k)*rho(i,k) else nragg(k)=0._r8 end if !....................................................................... ! Accretion of cloud ice by snow ! For this calculation, it is assumed that the Vs >> Vi ! and Ds >> Di for continuous collection if (do_cldice .and. qniic(i,k).ge.qsmall.and.qiic(i,k).ge.qsmall & .and.t(i,k).le.273.15_r8) then prai(k) = pi/4._r8*asn(i,k)*qiic(i,k)*rho(i,k)* & n0s(k)*Eii*cons11/ & lams(k)**(bs+3._r8) nprai(k) = pi/4._r8*asn(i,k)*niic(i,k)* & rho(i,k)*n0s(k)*Eii*cons11/ & lams(k)**(bs+3._r8) else prai(k)=0._r8 nprai(k)=0._r8 end if !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! calculate evaporation/sublimation of rain and snow ! note: evaporation/sublimation occurs only in cloud-free portion of grid cell ! in-cloud condensation/deposition of rain and snow is neglected ! except for transfer of cloud water to snow through bergeron process ! initialize evap/sub tendncies pre(k)=0._r8 prds(k)=0._r8 ! evaporation of rain ! only calculate if there is some precip fraction > cloud fraction if (qcic(i,k)+qiic(i,k).lt.1.e-6_r8.or.cldmax(i,k).gt.lcldm(i,k)) then ! set temporary cloud fraction to zero if cloud water + ice is very small ! this will ensure that evaporation/sublimation of precip occurs over ! entire grid cell, since min cloud fraction is specified otherwise if (qcic(i,k)+qiic(i,k).lt.1.e-6_r8) then dum=0._r8 else dum=lcldm(i,k) end if ! saturation vapor pressure esn=svp_water(t(i,k)) qsn=svp_to_qsat(esn, p(i,k)) ! recalculate saturation vapor pressure for liquid and ice esl(i,k)=esn esi(i,k)=svp_ice(t(i,k)) ! hm fix, make sure when above freezing that esi=esl, not active yet if (t(i,k).gt.tmelt)esi(i,k)=esl(i,k) ! calculate q for out-of-cloud region qclr=(q(i,k)-dum*qsn)/(1._r8-dum) if (qric(i,k).ge.qsmall) then qvs=svp_to_qsat(esl(i,k), p(i,k)) dqsdt = xxlv*qvs/(rv*t(i,k)**2) ab = 1._r8+dqsdt*xxlv/cpp epsr = 2._r8*pi*n0r(k)*rho(i,k)*Dv(i,k)* & (f1r/(lamr(k)*lamr(k))+ & f2r*(arn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* & sc(i,k)**(1._r8/3._r8)*cons13/ & (lamr(k)**(5._r8/2._r8+br/2._r8))) pre(k) = epsr*(qclr-qvs)/ab ! only evaporate in out-of-cloud region ! and distribute across cldmax pre(k)=min(pre(k)*(cldmax(i,k)-dum),0._r8) pre(k)=pre(k)/cldmax(i,k) am_evp_st(i,k) = max(cldmax(i,k)-dum, 0._r8) end if ! sublimation of snow if (qniic(i,k).ge.qsmall) then qvi=svp_to_qsat(esi(i,k), p(i,k)) dqsidt = xxls*qvi/(rv*t(i,k)**2) abi = 1._r8+dqsidt*xxls/cpp epss = 2._r8*pi*n0s(k)*rho(i,k)*Dv(i,k)* & (f1s/(lams(k)*lams(k))+ & f2s*(asn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* & sc(i,k)**(1._r8/3._r8)*cons14/ & (lams(k)**(5._r8/2._r8+bs/2._r8))) prds(k) = epss*(qclr-qvi)/abi ! only sublimate in out-of-cloud region and distribute over cldmax prds(k)=min(prds(k)*(cldmax(i,k)-dum),0._r8) prds(k)=prds(k)/cldmax(i,k) am_evp_st(i,k) = max(cldmax(i,k)-dum, 0._r8) end if ! make sure RH not pushed above 100% due to rain evaporation/snow sublimation ! get updated RH at end of time step based on cloud water/ice condensation/evap qtmp=q(i,k)-(cmei(i,k)+(pre(k)+prds(k))*cldmax(i,k))*deltat ttmp=t(i,k)+((pre(k)*cldmax(i,k))*xxlv+ & (cmei(i,k)+prds(k)*cldmax(i,k))*xxls)*deltat/cpp !limit range of temperatures! ttmp=max(180._r8,min(ttmp,323._r8)) esn=svp_water(ttmp) ! use rhw to allow ice supersaturation qsn=svp_to_qsat(esn, p(i,k)) ! modify precip evaporation rate if q > qsat if (qtmp.gt.qsn) then if (pre(k)+prds(k).lt.-1.e-20_r8) then dum1=pre(k)/(pre(k)+prds(k)) ! recalculate q and t after cloud water cond but without precip evap qtmp=q(i,k)-(cmei(i,k))*deltat ttmp=t(i,k)+(cmei(i,k)*xxls)*deltat/cpp esn=svp_water(ttmp) ! use rhw to allow ice supersaturation qsn=svp_to_qsat(esn, p(i,k)) dum=(qtmp-qsn)/(1._r8 + cons27*qsn/(cpp*rv*ttmp**2)) dum=min(dum,0._r8) ! modify rates if needed, divide by cldmax to get local (in-precip) value pre(k)=dum*dum1/deltat/cldmax(i,k) ! do separately using RHI for prds.... esn=svp_ice(ttmp) ! use rhi to allow ice supersaturation qsn=svp_to_qsat(esn, p(i,k)) dum=(qtmp-qsn)/(1._r8 + cons28*qsn/(cpp*rv*ttmp**2)) dum=min(dum,0._r8) ! modify rates if needed, divide by cldmax to get local (in-precip) value prds(k)=dum*(1._r8-dum1)/deltat/cldmax(i,k) end if end if end if ! bergeron process - evaporation of droplets and deposition onto snow if (qniic(i,k).ge.qsmall.and.qcic(i,k).ge.qsmall.and.t(i,k).lt.tmelt) then qvi=svp_to_qsat(esi(i,k), p(i,k)) qvs=svp_to_qsat(esl(i,k), p(i,k)) dqsidt = xxls*qvi/(rv*t(i,k)**2) abi = 1._r8+dqsidt*xxls/cpp epss = 2._r8*pi*n0s(k)*rho(i,k)*Dv(i,k)* & (f1s/(lams(k)*lams(k))+ & f2s*(asn(i,k)*rho(i,k)/mu(i,k))**0.5_r8* & sc(i,k)**(1._r8/3._r8)*cons14/ & (lams(k)**(5._r8/2._r8+bs/2._r8))) bergs(k)=epss*(qvs-qvi)/abi else bergs(k)=0._r8 end if !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! conservation to ensure no negative values of cloud water/precipitation ! in case microphysical process rates are large ! make sure and use end-of-time step values for cloud water, ice, due ! condensation/deposition ! note: for check on conservation, processes are multiplied by omsm ! to prevent problems due to round off error ! include mixing timescale (mtime) qce=(qc(i,k) - berg(i,k)*deltat) nce=(nc(i,k)+npccn(k)*deltat*mtime) qie=(qi(i,k)+(cmei(i,k)+berg(i,k))*deltat) nie=(ni(i,k)+nnuccd(k)*deltat*mtime) ! conservation of qc dum = (prc(k)+pra(k)+mnuccc(k)+mnucct(k)+msacwi(k)+ & psacws(k)+bergs(k))*lcldm(i,k)*deltat if (dum.gt.qce) then ratio = qce/deltat/lcldm(i,k)/(prc(k)+pra(k)+mnuccc(k)+mnucct(k)+msacwi(k)+psacws(k)+bergs(k))*omsm prc(k) = prc(k)*ratio pra(k) = pra(k)*ratio mnuccc(k) = mnuccc(k)*ratio mnucct(k) = mnucct(k)*ratio msacwi(k) = msacwi(k)*ratio psacws(k) = psacws(k)*ratio bergs(k) = bergs(k)*ratio end if ! conservation of nc dum = (nprc1(k)+npra(k)+nnuccc(k)+nnucct(k)+ & npsacws(k)-nsubc(k))*lcldm(i,k)*deltat if (dum.gt.nce) then ratio = nce/deltat/((nprc1(k)+npra(k)+nnuccc(k)+nnucct(k)+& npsacws(k)-nsubc(k))*lcldm(i,k))*omsm nprc1(k) = nprc1(k)*ratio npra(k) = npra(k)*ratio nnuccc(k) = nnuccc(k)*ratio nnucct(k) = nnucct(k)*ratio npsacws(k) = npsacws(k)*ratio nsubc(k)=nsubc(k)*ratio end if ! conservation of qi if (do_cldice) then frztmp = -mnuccc(k) - mnucct(k) - msacwi(k) if (use_hetfrz_classnuc) frztmp = -mnuccc(k)-mnucct(k)-mnudep(k)-msacwi(k) dum = ( frztmp*lcldm(i,k) + (prci(k)+prai(k))*icldm(i,k) )*deltat if (dum.gt.qie) then frztmp = mnuccc(k) + mnucct(k) + msacwi(k) if (use_hetfrz_classnuc) frztmp = mnuccc(k) + mnucct(k) + mnudep(k) + msacwi(k) ratio = (qie/deltat + frztmp*lcldm(i,k))/((prci(k)+prai(k))*icldm(i,k))*omsm prci(k) = prci(k)*ratio prai(k) = prai(k)*ratio end if ! conservation of ni frztmp = -nnucct(k) - nsacwi(k) if (use_hetfrz_classnuc) frztmp = -nnucct(k) - nnuccc(k) - nnudep(k) - nsacwi(k) dum = ( frztmp*lcldm(i,k) + (nprci(k)+nprai(k)-nsubi(k))*icldm(i,k) )*deltat if (dum.gt.nie) then frztmp = nnucct(k) + nsacwi(k) if (use_hetfrz_classnuc) frztmp = nnucct(k) + nnuccc(k) + nnudep(k) + nsacwi(k) ratio = (nie/deltat + frztmp*lcldm(i,k))/ & ((nprci(k)+nprai(k)-nsubi(k))*icldm(i,k))*omsm nprci(k) = nprci(k)*ratio nprai(k) = nprai(k)*ratio nsubi(k) = nsubi(k)*ratio end if end if ! for precipitation conservation, use logic that vertical integral ! of tendency from current level to top of model (i.e., qrtot) cannot be negative ! conservation of rain mixing rat if (((prc(k)+pra(k))*lcldm(i,k)+(-mnuccr(k)+pre(k)-pracs(k))*& cldmax(i,k))*dz(i,k)*rho(i,k)+qrtot.lt.0._r8) then if (-pre(k)+pracs(k)+mnuccr(k).ge.qsmall) then ratio = (qrtot/(dz(i,k)*rho(i,k))+(prc(k)+pra(k))*lcldm(i,k))/& ((-pre(k)+pracs(k)+mnuccr(k))*cldmax(i,k))*omsm pre(k) = pre(k)*ratio pracs(k) = pracs(k)*ratio mnuccr(k) = mnuccr(k)*ratio end if end if ! conservation of nr ! for now neglect evaporation of nr nsubr(k)=0._r8 if ((nprc(k)*lcldm(i,k)+(-nnuccr(k)+nsubr(k)-npracs(k)& +nragg(k))*cldmax(i,k))*dz(i,k)*rho(i,k)+nrtot.lt.0._r8) then if (-nsubr(k)-nragg(k)+npracs(k)+nnuccr(k).ge.qsmall) then ratio = (nrtot/(dz(i,k)*rho(i,k))+nprc(k)*lcldm(i,k))/& ((-nsubr(k)-nragg(k)+npracs(k)+nnuccr(k))*cldmax(i,k))*omsm nsubr(k) = nsubr(k)*ratio npracs(k) = npracs(k)*ratio nnuccr(k) = nnuccr(k)*ratio nragg(k) = nragg(k)*ratio end if end if ! conservation of snow mix ratio if (((bergs(k)+psacws(k))*lcldm(i,k)+(prai(k)+prci(k))*icldm(i,k)+(pracs(k)+& mnuccr(k)+prds(k))*cldmax(i,k))*dz(i,k)*rho(i,k)+qstot.lt.0._r8) then if (-prds(k).ge.qsmall) then ratio = (qstot/(dz(i,k)*rho(i,k))+(bergs(k)+psacws(k))*lcldm(i,k)+(prai(k)+prci(k))*icldm(i,k)+& (pracs(k)+mnuccr(k))*cldmax(i,k))/(-prds(k)*cldmax(i,k))*omsm prds(k) = prds(k)*ratio end if end if ! conservation of ns ! calculate loss of number due to sublimation ! for now neglect sublimation of ns nsubs(k)=0._r8 if ((nprci(k)*icldm(i,k)+(nnuccr(k)+nsubs(k)+nsagg(k))*cldmax(i,k))*& dz(i,k)*rho(i,k)+nstot.lt.0._r8) then if (-nsubs(k)-nsagg(k).ge.qsmall) then ratio = (nstot/(dz(i,k)*rho(i,k))+nprci(k)*icldm(i,k)+& nnuccr(k)*cldmax(i,k))/((-nsubs(k)-nsagg(k))*cldmax(i,k))*omsm nsubs(k) = nsubs(k)*ratio nsagg(k) = nsagg(k)*ratio end if end if ! get tendencies due to microphysical conversion processes ! note: tendencies are multiplied by appropaiate cloud/precip ! fraction to get grid-scale values ! note: cmei is already grid-average values qvlat(i,k) = qvlat(i,k)-(pre(k)+prds(k))*cldmax(i,k)-cmei(i,k) tlat(i,k) = tlat(i,k)+((pre(k)*cldmax(i,k)) & *xxlv+(prds(k)*cldmax(i,k)+cmei(i,k))*xxls+ & ((bergs(k)+psacws(k)+mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k)+(mnuccr(k)+ & pracs(k))*cldmax(i,k)+berg(i,k))*xlf) qctend(i,k) = qctend(i,k)+ & (-pra(k)-prc(k)-mnuccc(k)-mnucct(k)-msacwi(k)- & psacws(k)-bergs(k))*lcldm(i,k)-berg(i,k) if (do_cldice) then frztmp = mnuccc(k) + mnucct(k) + msacwi(k) if (use_hetfrz_classnuc) frztmp = mnuccc(k) + mnucct(k) + mnudep(k) + msacwi(k) qitend(i,k) = qitend(i,k) + frztmp*lcldm(i,k) + & (-prci(k)-prai(k))*icldm(i,k) + cmei(i,k) + berg(i,k) end if qrtend(i,k) = qrtend(i,k)+ & (pra(k)+prc(k))*lcldm(i,k)+(pre(k)-pracs(k)- & mnuccr(k))*cldmax(i,k) qnitend(i,k) = qnitend(i,k)+ & (prai(k)+prci(k))*icldm(i,k)+(psacws(k)+bergs(k))*lcldm(i,k)+(prds(k)+ & pracs(k)+mnuccr(k))*cldmax(i,k) ! add output for cmei (accumulate) cmeiout(i,k) = cmeiout(i,k) + cmei(i,k) ! assign variables for trop_mozart, these are grid-average ! evaporation/sublimation is stored here as positive term evapsnow(i,k) = evapsnow(i,k)-prds(k)*cldmax(i,k) nevapr(i,k) = nevapr(i,k)-pre(k)*cldmax(i,k) nevapr2(i,k) = nevapr2(i,k)-pre(k)*cldmax(i,k) ! change to make sure prain is positive: do not remove snow from ! prain used for wet deposition prain(i,k) = prain(i,k)+(pra(k)+prc(k))*lcldm(i,k)+(-pracs(k)- & mnuccr(k))*cldmax(i,k) prodsnow(i,k) = prodsnow(i,k)+(prai(k)+prci(k))*icldm(i,k)+(psacws(k)+bergs(k))*lcldm(i,k)+(& pracs(k)+mnuccr(k))*cldmax(i,k) ! following are used to calculate 1st order conversion rate of cloud water ! to rain and snow (1/s), for later use in aerosol wet removal routine ! previously, wetdepa used (prain/qc) for this, and the qc in wetdepa may be smaller than the qc ! used to calculate pra, prc, ... in this routine ! qcsinksum_rate1ord = sum over iterations{ rate of direct transfer of cloud water to rain & snow } ! (no cloud ice or bergeron terms) ! qcsum_rate1ord = sum over iterations{ qc used in calculation of the transfer terms } qcsinksum_rate1ord(k) = qcsinksum_rate1ord(k) + (pra(k)+prc(k)+psacws(k))*lcldm(i,k) qcsum_rate1ord(k) = qcsum_rate1ord(k) + qc(i,k) ! microphysics output, note this is grid-averaged prao(i,k)=prao(i,k)+pra(k)*lcldm(i,k) prco(i,k)=prco(i,k)+prc(k)*lcldm(i,k) mnuccco(i,k)=mnuccco(i,k)+mnuccc(k)*lcldm(i,k) mnuccto(i,k)=mnuccto(i,k)+mnucct(k)*lcldm(i,k) mnuccdo(i,k)=mnuccdo(i,k)+mnuccd(k)*lcldm(i,k) msacwio(i,k)=msacwio(i,k)+msacwi(k)*lcldm(i,k) psacwso(i,k)=psacwso(i,k)+psacws(k)*lcldm(i,k) bergso(i,k)=bergso(i,k)+bergs(k)*lcldm(i,k) bergo(i,k)=bergo(i,k)+berg(i,k) prcio(i,k)=prcio(i,k)+prci(k)*icldm(i,k) praio(i,k)=praio(i,k)+prai(k)*icldm(i,k) mnuccro(i,k)=mnuccro(i,k)+mnuccr(k)*cldmax(i,k) pracso (i,k)=pracso (i,k)+pracs (k)*cldmax(i,k) ! multiply activation/nucleation by mtime to account for fast timescale nctend(i,k) = nctend(i,k)+ npccn(k)*mtime+& (-nnuccc(k)-nnucct(k)-npsacws(k)+nsubc(k) & -npra(k)-nprc1(k))*lcldm(i,k) if (do_cldice) then frztmp = nnucct(k) + nsacwi(k) if (use_hetfrz_classnuc) frztmp = nnucct(k) + nnuccc(k) + nnudep(k) + nsacwi(k) nitend(i,k) = nitend(i,k) + nnuccd(k)*mtime + & frztmp*lcldm(i,k) + (nsubi(k)-nprci(k)-nprai(k))*icldm(i,k) end if nstend(i,k) = nstend(i,k)+(nsubs(k)+ & nsagg(k)+nnuccr(k))*cldmax(i,k)+nprci(k)*icldm(i,k) nrtend(i,k) = nrtend(i,k)+ & nprc(k)*lcldm(i,k)+(nsubr(k)-npracs(k)-nnuccr(k) & +nragg(k))*cldmax(i,k) ! make sure that nc and ni at advanced time step do not exceed ! maximum (existing N + source terms*dt), which is possible due to ! fast nucleation timescale if (nctend(i,k).gt.0._r8.and.nc(i,k)+nctend(i,k)*deltat.gt.ncmax) then nctend(i,k)=max(0._r8,(ncmax-nc(i,k))/deltat) end if if (do_cldice .and. nitend(i,k).gt.0._r8.and.ni(i,k)+nitend(i,k)*deltat.gt.nimax) then nitend(i,k)=max(0._r8,(nimax-ni(i,k))/deltat) end if ! get final values for precipitation q and N, based on ! flux of precip from above, source/sink term, and terminal fallspeed ! see eq. 15-16 in MG2008 ! rain if (qric(i,k).ge.qsmall) then if (k.eq.top_lev) then qric(i,k)=qrtend(i,k)*dz(i,k)/cldmax(i,k)/umr(k) nric(i,k)=nrtend(i,k)*dz(i,k)/cldmax(i,k)/unr(k) else qric(i,k) = (rho(i,k-1)*umr(k-1)*qric(i,k-1)*cldmax(i,k-1)+ & (rho(i,k)*dz(i,k)*qrtend(i,k)))/(umr(k)*rho(i,k)*cldmax(i,k)) nric(i,k) = (rho(i,k-1)*unr(k-1)*nric(i,k-1)*cldmax(i,k-1)+ & (rho(i,k)*dz(i,k)*nrtend(i,k)))/(unr(k)*rho(i,k)*cldmax(i,k)) end if else qric(i,k)=0._r8 nric(i,k)=0._r8 end if ! snow if (qniic(i,k).ge.qsmall) then if (k.eq.top_lev) then qniic(i,k)=qnitend(i,k)*dz(i,k)/cldmax(i,k)/ums(k) nsic(i,k)=nstend(i,k)*dz(i,k)/cldmax(i,k)/uns(k) else qniic(i,k) = (rho(i,k-1)*ums(k-1)*qniic(i,k-1)*cldmax(i,k-1)+ & (rho(i,k)*dz(i,k)*qnitend(i,k)))/(ums(k)*rho(i,k)*cldmax(i,k)) nsic(i,k) = (rho(i,k-1)*uns(k-1)*nsic(i,k-1)*cldmax(i,k-1)+ & (rho(i,k)*dz(i,k)*nstend(i,k)))/(uns(k)*rho(i,k)*cldmax(i,k)) end if else qniic(i,k)=0._r8 nsic(i,k)=0._r8 end if ! calculate precipitation flux at surface ! divide by density of water to get units of m/s prect(i) = prect(i)+(qrtend(i,k)*dz(i,k)*rho(i,k)+& qnitend(i,k)*dz(i,k)*rho(i,k))/rhow preci(i) = preci(i)+qnitend(i,k)*dz(i,k)*rho(i,k)/rhow ! convert rain rate from m/s to mm/hr rainrt(i,k)=qric(i,k)*rho(i,k)*umr(k)/rhow*3600._r8*1000._r8 ! vertically-integrated precip source/sink terms (note: grid-averaged) qrtot = max(qrtot+qrtend(i,k)*dz(i,k)*rho(i,k),0._r8) qstot = max(qstot+qnitend(i,k)*dz(i,k)*rho(i,k),0._r8) nrtot = max(nrtot+nrtend(i,k)*dz(i,k)*rho(i,k),0._r8) nstot = max(nstot+nstend(i,k)*dz(i,k)*rho(i,k),0._r8) ! calculate melting and freezing of precip ! melt snow at +2 C if (t(i,k)+tlat(i,k)/cpp*deltat > 275.15_r8) then if (qstot > 0._r8) then ! make sure melting snow doesn't reduce temperature below threshold dum = -xlf/cpp*qstot/(dz(i,k)*rho(i,k)) if (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt.275.15_r8) then dum = (t(i,k)+tlat(i,k)/cpp*deltat-275.15_r8)*cpp/xlf dum = dum/(xlf/cpp*qstot/(dz(i,k)*rho(i,k))) dum = max(0._r8,dum) dum = min(1._r8,dum) else dum = 1._r8 end if qric(i,k)=qric(i,k)+dum*qniic(i,k) nric(i,k)=nric(i,k)+dum*nsic(i,k) qniic(i,k)=(1._r8-dum)*qniic(i,k) nsic(i,k)=(1._r8-dum)*nsic(i,k) ! heating tendency tmp=-xlf*dum*qstot/(dz(i,k)*rho(i,k)) meltsdt(i,k)=meltsdt(i,k) + tmp tlat(i,k)=tlat(i,k)+tmp qrtot=qrtot+dum*qstot nrtot=nrtot+dum*nstot qstot=(1._r8-dum)*qstot nstot=(1._r8-dum)*nstot preci(i)=(1._r8-dum)*preci(i) end if end if ! freeze all rain at -5C for Arctic if (t(i,k)+tlat(i,k)/cpp*deltat < (tmelt - 5._r8)) then if (qrtot > 0._r8) then ! make sure freezing rain doesn't increase temperature above threshold dum = xlf/cpp*qrtot/(dz(i,k)*rho(i,k)) if (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.(tmelt - 5._r8)) then dum = -(t(i,k)+tlat(i,k)/cpp*deltat-(tmelt-5._r8))*cpp/xlf dum = dum/(xlf/cpp*qrtot/(dz(i,k)*rho(i,k))) dum = max(0._r8,dum) dum = min(1._r8,dum) else dum = 1._r8 end if qniic(i,k)=qniic(i,k)+dum*qric(i,k) nsic(i,k)=nsic(i,k)+dum*nric(i,k) qric(i,k)=(1._r8-dum)*qric(i,k) nric(i,k)=(1._r8-dum)*nric(i,k) ! heating tendency tmp = xlf*dum*qrtot/(dz(i,k)*rho(i,k)) frzrdt(i,k)=frzrdt(i,k) + tmp tlat(i,k)=tlat(i,k)+tmp qstot=qstot+dum*qrtot qrtot=(1._r8-dum)*qrtot nstot=nstot+dum*nrtot nrtot=(1._r8-dum)*nrtot preci(i)=preci(i)+dum*(prect(i)-preci(i)) end if end if ! Precip Flux Calculation (Diagnostic) rflx(i,k+1)=(prect(i)-preci(i)) * rhow sflx(i,k+1)=preci(i) * rhow ! if rain/snow mix ratio is zero so should number concentration if (qniic(i,k).lt.qsmall) then qniic(i,k)=0._r8 nsic(i,k)=0._r8 end if if (qric(i,k).lt.qsmall) then qric(i,k)=0._r8 nric(i,k)=0._r8 end if ! make sure number concentration is a positive number to avoid ! taking root of negative nric(i,k)=max(nric(i,k),0._r8) nsic(i,k)=max(nsic(i,k),0._r8) !....................................................................... ! get size distribution parameters for fallspeed calculations !...................................................................... ! rain if (qric(i,k).ge.qsmall) then lamr(k) = (pi*rhow*nric(i,k)/qric(i,k))**(1._r8/3._r8) n0r(k) = nric(i,k)*lamr(k) ! check for slope ! change lammax and lammin for rain and snow ! adjust vars if (lamr(k).lt.lamminr) then lamr(k) = lamminr n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow) nric(i,k) = n0r(k)/lamr(k) else if (lamr(k).gt.lammaxr) then lamr(k) = lammaxr n0r(k) = lamr(k)**4*qric(i,k)/(pi*rhow) nric(i,k) = n0r(k)/lamr(k) end if ! 'final' values of number and mass weighted mean fallspeed for rain (m/s) unr(k) = min(arn(i,k)*cons4/lamr(k)**br,9.1_r8*rhof(i,k)) umr(k) = min(arn(i,k)*cons5/(6._r8*lamr(k)**br),9.1_r8*rhof(i,k)) else lamr(k) = 0._r8 n0r(k) = 0._r8 umr(k)=0._r8 unr(k)=0._r8 end if !calculate mean size of combined rain and snow if (lamr(k).gt.0._r8) then Artmp = n0r(k) * pi / (2._r8 * lamr(k)**3._r8) else Artmp = 0._r8 endif if (lamc(k).gt.0._r8) then Actmp = cdist1(k) * pi * gamma(pgam(k)+3._r8)/(4._r8 * lamc(k)**2._r8) else Actmp = 0._r8 endif if (Actmp.gt.0_r8.or.Artmp.gt.0) then rercld(i,k)=rercld(i,k) + 3._r8 *(qric(i,k) + qcic(i,k)) / (4._r8 * rhow * (Actmp + Artmp)) arcld(i,k)=arcld(i,k)+1._r8 endif !...................................................................... ! snow if (qniic(i,k).ge.qsmall) then lams(k) = (cons6*cs*nsic(i,k)/ & qniic(i,k))**(1._r8/ds) n0s(k) = nsic(i,k)*lams(k) ! check for slope ! adjust vars if (lams(k).lt.lammins) then lams(k) = lammins n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6) nsic(i,k) = n0s(k)/lams(k) else if (lams(k).gt.lammaxs) then lams(k) = lammaxs n0s(k) = lams(k)**(ds+1._r8)*qniic(i,k)/(cs*cons6) nsic(i,k) = n0s(k)/lams(k) end if ! 'final' values of number and mass weighted mean fallspeed for snow (m/s) ums(k) = min(asn(i,k)*cons8/(6._r8*lams(k)**bs),1.2_r8*rhof(i,k)) uns(k) = min(asn(i,k)*cons7/lams(k)**bs,1.2_r8*rhof(i,k)) else lams(k) = 0._r8 n0s(k) = 0._r8 ums(k) = 0._r8 uns(k) = 0._r8 end if !c........................................................................ ! sum over sub-step for average process rates ! convert rain/snow q and N for output to history, note, ! output is for gridbox average qrout(i,k)=qrout(i,k)+qric(i,k)*cldmax(i,k) qsout(i,k)=qsout(i,k)+qniic(i,k)*cldmax(i,k) nrout(i,k)=nrout(i,k)+nric(i,k)*rho(i,k)*cldmax(i,k) nsout(i,k)=nsout(i,k)+nsic(i,k)*rho(i,k)*cldmax(i,k) tlat1(i,k)=tlat1(i,k)+tlat(i,k) qvlat1(i,k)=qvlat1(i,k)+qvlat(i,k) qctend1(i,k)=qctend1(i,k)+qctend(i,k) qitend1(i,k)=qitend1(i,k)+qitend(i,k) nctend1(i,k)=nctend1(i,k)+nctend(i,k) nitend1(i,k)=nitend1(i,k)+nitend(i,k) t(i,k)=t(i,k)+tlat(i,k)*deltat/cpp q(i,k)=q(i,k)+qvlat(i,k)*deltat qc(i,k)=qc(i,k)+qctend(i,k)*deltat qi(i,k)=qi(i,k)+qitend(i,k)*deltat nc(i,k)=nc(i,k)+nctend(i,k)*deltat ni(i,k)=ni(i,k)+nitend(i,k)*deltat rainrt1(i,k)=rainrt1(i,k)+rainrt(i,k) !divide rain radius over substeps for average if (arcld(i,k) .gt. 0._r8) then rercld(i,k)=rercld(i,k)/arcld(i,k) end if !! add to summing sub-stepping variable rflx1(i,k+1)=rflx1(i,k+1)+rflx(i,k+1) sflx1(i,k+1)=sflx1(i,k+1)+sflx(i,k+1) !c........................................................................ end do ! k loop prect1(i)=prect1(i)+prect(i) preci1(i)=preci1(i)+preci(i) end do ! it loop, sub-step do k = top_lev, pver rate1ord_cw2pr_st(i,k) = qcsinksum_rate1ord(k)/max(qcsum_rate1ord(k),1.0e-30_r8) end do 300 continue ! continue if no cloud water end do ! i loop ! convert dt from sub-step back to full time step deltat=deltat*real(iter) !c............................................................................. do i=1,ncol ! skip all calculations if no cloud water if (ltrue(i).eq.0) then do k=1,top_lev-1 ! assign zero values for effective radius above 1 mbar effc(i,k)=0._r8 effi(i,k)=0._r8 effc_fn(i,k)=0._r8 lamcrad(i,k)=0._r8 pgamrad(i,k)=0._r8 deffi(i,k)=0._r8 end do do k=top_lev,pver ! assign default values for effective radius effc(i,k)=10._r8 effi(i,k)=25._r8 effc_fn(i,k)=10._r8 lamcrad(i,k)=0._r8 pgamrad(i,k)=0._r8 deffi(i,k)=0._r8 end do goto 500 end if ! initialize nstep for sedimentation sub-steps nstep = 1 ! divide precip rate by number of sub-steps to get average over time step prect(i)=prect1(i)/real(iter) preci(i)=preci1(i)/real(iter) do k=top_lev,pver ! assign variables back to start-of-timestep values before updating after sub-steps t(i,k)=t1(i,k) q(i,k)=q1(i,k) qc(i,k)=qc1(i,k) qi(i,k)=qi1(i,k) nc(i,k)=nc1(i,k) ni(i,k)=ni1(i,k) ! divide microphysical tendencies by number of sub-steps to get average over time step tlat(i,k)=tlat1(i,k)/real(iter) qvlat(i,k)=qvlat1(i,k)/real(iter) qctend(i,k)=qctend1(i,k)/real(iter) qitend(i,k)=qitend1(i,k)/real(iter) nctend(i,k)=nctend1(i,k)/real(iter) nitend(i,k)=nitend1(i,k)/real(iter) rainrt(i,k)=rainrt1(i,k)/real(iter) ! divide by number of sub-steps to find final values rflx(i,k+1)=rflx1(i,k+1)/real(iter) sflx(i,k+1)=sflx1(i,k+1)/real(iter) ! divide output precip q and N by number of sub-steps to get average over time step qrout(i,k)=qrout(i,k)/real(iter) qsout(i,k)=qsout(i,k)/real(iter) nrout(i,k)=nrout(i,k)/real(iter) nsout(i,k)=nsout(i,k)/real(iter) ! divide trop_mozart variables by number of sub-steps to get average over time step nevapr(i,k) = nevapr(i,k)/real(iter) nevapr2(i,k) = nevapr2(i,k)/real(iter) evapsnow(i,k) = evapsnow(i,k)/real(iter) prain(i,k) = prain(i,k)/real(iter) prodsnow(i,k) = prodsnow(i,k)/real(iter) cmeout(i,k) = cmeout(i,k)/real(iter) cmeiout(i,k) = cmeiout(i,k)/real(iter) meltsdt(i,k) = meltsdt(i,k)/real(iter) frzrdt (i,k) = frzrdt (i,k)/real(iter) ! microphysics output prao(i,k)=prao(i,k)/real(iter) prco(i,k)=prco(i,k)/real(iter) mnuccco(i,k)=mnuccco(i,k)/real(iter) mnuccto(i,k)=mnuccto(i,k)/real(iter) msacwio(i,k)=msacwio(i,k)/real(iter) psacwso(i,k)=psacwso(i,k)/real(iter) bergso(i,k)=bergso(i,k)/real(iter) bergo(i,k)=bergo(i,k)/real(iter) prcio(i,k)=prcio(i,k)/real(iter) praio(i,k)=praio(i,k)/real(iter) mnuccro(i,k)=mnuccro(i,k)/real(iter) pracso (i,k)=pracso (i,k)/real(iter) mnuccdo(i,k)=mnuccdo(i,k)/real(iter) ! modify to include snow. in prain & evap (diagnostic here: for wet dep) nevapr(i,k) = nevapr(i,k) + evapsnow(i,k) prer_evap(i,k) = nevapr2(i,k) prain(i,k) = prain(i,k) + prodsnow(i,k) !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! calculate sedimentation for cloud water and ice !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! update in-cloud cloud mixing ratio and number concentration ! with microphysical tendencies to calculate sedimentation, assign to dummy vars ! note: these are in-cloud values***, hence we divide by cloud fraction dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k) dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)/icldm(i,k) dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k),0._r8) dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)/icldm(i,k),0._r8) ! obtain new slope parameter to avoid possible singularity if (dumi(i,k).ge.qsmall) then ! add upper limit to in-cloud number concentration to prevent numerical error dumni(i,k)=min(dumni(i,k),dumi(i,k)*1.e20_r8) lami(k) = (cons1*ci* & dumni(i,k)/dumi(i,k))**(1._r8/di) lami(k)=max(lami(k),lammini) lami(k)=min(lami(k),lammaxi) else lami(k)=0._r8 end if if (dumc(i,k).ge.qsmall) then ! add upper limit to in-cloud number concentration to prevent numerical error dumnc(i,k)=min(dumnc(i,k),dumc(i,k)*1.e20_r8) ! add lower limit to in-cloud number concentration dumnc(i,k)=max(dumnc(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3 pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8 pgam(k)=1._r8/(pgam(k)**2)-1._r8 pgam(k)=max(pgam(k),2._r8) pgam(k)=min(pgam(k),15._r8) lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ & (dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8) lammin = (pgam(k)+1._r8)/50.e-6_r8 lammax = (pgam(k)+1._r8)/2.e-6_r8 lamc(k)=max(lamc(k),lammin) lamc(k)=min(lamc(k),lammax) else lamc(k)=0._r8 end if ! calculate number and mass weighted fall velocity for droplets ! include effects of sub-grid distribution of cloud water if (dumc(i,k).ge.qsmall) then unc = acn(i,k)*gamma(1._r8+bc+pgam(k))/(lamc(k)**bc*gamma(pgam(k)+1._r8)) umc = acn(i,k)*gamma(4._r8+bc+pgam(k))/(lamc(k)**bc*gamma(pgam(k)+4._r8)) ! fallspeed for output vtrmc(i,k)=umc else umc = 0._r8 unc = 0._r8 end if ! calculate number and mass weighted fall velocity for cloud ice if (dumi(i,k).ge.qsmall) then uni = ain(i,k)*cons16/lami(k)**bi umi = ain(i,k)*cons17/(6._r8*lami(k)**bi) uni=min(uni,1.2_r8*rhof(i,k)) umi=min(umi,1.2_r8*rhof(i,k)) ! fallspeed vtrmi(i,k)=umi else umi = 0._r8 uni = 0._r8 end if fi(k) = g*rho(i,k)*umi fni(k) = g*rho(i,k)*uni fc(k) = g*rho(i,k)*umc fnc(k) = g*rho(i,k)*unc ! calculate number of split time steps to ensure courant stability criteria ! for sedimentation calculations rgvm = max(fi(k),fc(k),fni(k),fnc(k)) nstep = max(int(rgvm*deltat/pdel(i,k)+1._r8),nstep) ! redefine dummy variables - sedimentation is calculated over grid-scale ! quantities to ensure conservation dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat) dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat) dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat),0._r8) dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat),0._r8) if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8 if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8 end do !!! vertical loop do n = 1,nstep !! loop over sub-time step to ensure stability do k = top_lev,pver if (do_cldice) then falouti(k) = fi(k)*dumi(i,k) faloutni(k) = fni(k)*dumni(i,k) else falouti(k) = 0._r8 faloutni(k) = 0._r8 end if faloutc(k) = fc(k)*dumc(i,k) faloutnc(k) = fnc(k)*dumnc(i,k) end do ! top of model k = top_lev faltndi = falouti(k)/pdel(i,k) faltndni = faloutni(k)/pdel(i,k) faltndc = faloutc(k)/pdel(i,k) faltndnc = faloutnc(k)/pdel(i,k) ! add fallout terms to microphysical tendencies qitend(i,k) = qitend(i,k)-faltndi/nstep nitend(i,k) = nitend(i,k)-faltndni/nstep qctend(i,k) = qctend(i,k)-faltndc/nstep nctend(i,k) = nctend(i,k)-faltndnc/nstep ! sedimentation tendencies for output qcsedten(i,k)=qcsedten(i,k)-faltndc/nstep qisedten(i,k)=qisedten(i,k)-faltndi/nstep dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep do k = top_lev+1,pver ! for cloud liquid and ice, if cloud fraction increases with height ! then add flux from above to both vapor and cloud water of current level ! this means that flux entering clear portion of cell from above evaporates ! instantly dum=lcldm(i,k)/lcldm(i,k-1) dum=min(dum,1._r8) dum1=icldm(i,k)/icldm(i,k-1) dum1=min(dum1,1._r8) faltndqie=(falouti(k)-falouti(k-1))/pdel(i,k) faltndi=(falouti(k)-dum1*falouti(k-1))/pdel(i,k) faltndni=(faloutni(k)-dum1*faloutni(k-1))/pdel(i,k) faltndqce=(faloutc(k)-faloutc(k-1))/pdel(i,k) faltndc=(faloutc(k)-dum*faloutc(k-1))/pdel(i,k) faltndnc=(faloutnc(k)-dum*faloutnc(k-1))/pdel(i,k) ! add fallout terms to eulerian tendencies qitend(i,k) = qitend(i,k)-faltndi/nstep nitend(i,k) = nitend(i,k)-faltndni/nstep qctend(i,k) = qctend(i,k)-faltndc/nstep nctend(i,k) = nctend(i,k)-faltndnc/nstep ! sedimentation tendencies for output qcsedten(i,k)=qcsedten(i,k)-faltndc/nstep qisedten(i,k)=qisedten(i,k)-faltndi/nstep ! add terms to to evap/sub of cloud water qvlat(i,k)=qvlat(i,k)-(faltndqie-faltndi)/nstep ! for output qisevap(i,k)=qisevap(i,k)-(faltndqie-faltndi)/nstep qvlat(i,k)=qvlat(i,k)-(faltndqce-faltndc)/nstep ! for output qcsevap(i,k)=qcsevap(i,k)-(faltndqce-faltndc)/nstep tlat(i,k)=tlat(i,k)+(faltndqie-faltndi)*xxls/nstep tlat(i,k)=tlat(i,k)+(faltndqce-faltndc)*xxlv/nstep dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep Fni(K)=MAX(Fni(K)/pdel(i,K),Fni(K-1)/pdel(i,K-1))*pdel(i,K) FI(K)=MAX(FI(K)/pdel(i,K),FI(K-1)/pdel(i,K-1))*pdel(i,K) fnc(k)=max(fnc(k)/pdel(i,k),fnc(k-1)/pdel(i,k-1))*pdel(i,k) Fc(K)=MAX(Fc(K)/pdel(i,K),Fc(K-1)/pdel(i,K-1))*pdel(i,K) end do !! k loop ! units below are m/s ! cloud water/ice sedimentation flux at surface ! is added to precip flux at surface to get total precip (cloud + precip water) ! rate prect(i) = prect(i)+(faloutc(pver)+falouti(pver))/g/nstep/1000._r8 preci(i) = preci(i)+(falouti(pver))/g/nstep/1000._r8 ! Add fallout to Precip Flux: note unit change m/s *kg/m3 = kg/m2 do k = top_lev,pver rflx(i,k+1)=rflx(i,k+1)+(faloutc(k))/g/nstep/1000._r8 * rhow sflx(i,k+1)=sflx(i,k+1)+(falouti(k))/g/nstep/1000._r8 * rhow end do end do !! nstep loop ! end sedimentation !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! get new update for variables that includes sedimentation tendency ! note : here dum variables are grid-average, NOT in-cloud do k=top_lev,pver dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8) dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8) dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8) dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8) if (dumc(i,k).lt.qsmall) dumnc(i,k)=0._r8 if (dumi(i,k).lt.qsmall) dumni(i,k)=0._r8 ! calculate instantaneous processes (melting, homogeneous freezing) if (do_cldice) then if (t(i,k)+tlat(i,k)/cpp*deltat > tmelt) then if (dumi(i,k) > 0._r8) then ! limit so that melting does not push temperature below freezing dum = -dumi(i,k)*xlf/cpp if (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt.tmelt) then dum = (t(i,k)+tlat(i,k)/cpp*deltat-tmelt)*cpp/xlf dum = dum/dumi(i,k)*xlf/cpp dum = max(0._r8,dum) dum = min(1._r8,dum) else dum = 1._r8 end if qctend(i,k)=qctend(i,k)+dum*dumi(i,k)/deltat ! for output melto(i,k)=dum*dumi(i,k)/deltat ! assume melting ice produces droplet ! mean volume radius of 8 micron nctend(i,k)=nctend(i,k)+3._r8*dum*dumi(i,k)/deltat/ & (4._r8*pi*5.12e-16_r8*rhow) qitend(i,k)=((1._r8-dum)*dumi(i,k)-qi(i,k))/deltat nitend(i,k)=((1._r8-dum)*dumni(i,k)-ni(i,k))/deltat tlat(i,k)=tlat(i,k)-xlf*dum*dumi(i,k)/deltat end if end if ! homogeneously freeze droplets at -40 C if (t(i,k)+tlat(i,k)/cpp*deltat < 233.15_r8) then if (dumc(i,k) > 0._r8) then ! limit so that freezing does not push temperature above threshold dum = dumc(i,k)*xlf/cpp if (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.233.15_r8) then dum = -(t(i,k)+tlat(i,k)/cpp*deltat-233.15_r8)*cpp/xlf dum = dum/dumc(i,k)*xlf/cpp dum = max(0._r8,dum) dum = min(1._r8,dum) else dum = 1._r8 end if qitend(i,k)=qitend(i,k)+dum*dumc(i,k)/deltat ! for output homoo(i,k)=dum*dumc(i,k)/deltat ! assume 25 micron mean volume radius of homogeneously frozen droplets ! consistent with size of detrained ice in stratiform.F90 nitend(i,k)=nitend(i,k)+dum*3._r8*dumc(i,k)/(4._r8*3.14_r8*1.563e-14_r8* & 500._r8)/deltat qctend(i,k)=((1._r8-dum)*dumc(i,k)-qc(i,k))/deltat nctend(i,k)=((1._r8-dum)*dumnc(i,k)-nc(i,k))/deltat tlat(i,k)=tlat(i,k)+xlf*dum*dumc(i,k)/deltat end if end if ! remove any excess over-saturation, which is possible due to non-linearity when adding ! together all microphysical processes ! follow code similar to old CAM scheme qtmp=q(i,k)+qvlat(i,k)*deltat ttmp=t(i,k)+tlat(i,k)/cpp*deltat esn = svp_water(ttmp) ! use rhw to allow ice supersaturation qsn = svp_to_qsat(esn, p(i,k)) if (qtmp > qsn .and. qsn > 0) then ! expression below is approximate since there may be ice deposition dum = (qtmp-qsn)/(1._r8+cons27*qsn/(cpp*rv*ttmp**2))/deltat ! add to output cme cmeout(i,k) = cmeout(i,k)+dum ! now add to tendencies, partition between liquid and ice based on temperature if (ttmp > 268.15_r8) then dum1=0.0_r8 ! now add to tendencies, partition between liquid and ice based on te else if (ttmp < 238.15_r8) then dum1=1.0_r8 else dum1=(268.15_r8-ttmp)/30._r8 end if dum = (qtmp-qsn)/(1._r8+(xxls*dum1+xxlv*(1._r8-dum1))**2 & *qsn/(cpp*rv*ttmp**2))/deltat qctend(i,k)=qctend(i,k)+dum*(1._r8-dum1) ! for output qcreso(i,k)=dum*(1._r8-dum1) qitend(i,k)=qitend(i,k)+dum*dum1 qireso(i,k)=dum*dum1 qvlat(i,k)=qvlat(i,k)-dum ! for output qvres(i,k)=-dum tlat(i,k)=tlat(i,k)+dum*(1._r8-dum1)*xxlv+dum*dum1*xxls end if end if !............................................................................... ! calculate effective radius for pass to radiation code ! if no cloud water, default value is 10 micron for droplets, ! 25 micron for cloud ice ! update cloud variables after instantaneous processes to get effective radius ! variables are in-cloud to calculate size dist parameters dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)/lcldm(i,k) dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)/icldm(i,k) dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)/lcldm(i,k) dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)/icldm(i,k) ! limit in-cloud mixing ratio to reasonable value of 5 g kg-1 dumc(i,k)=min(dumc(i,k),5.e-3_r8) dumi(i,k)=min(dumi(i,k),5.e-3_r8) !................... ! cloud ice effective radius if (dumi(i,k).ge.qsmall) then ! add upper limit to in-cloud number concentration to prevent numerical error dumni(i,k)=min(dumni(i,k),dumi(i,k)*1.e20_r8) lami(k) = (cons1*ci*dumni(i,k)/dumi(i,k))**(1._r8/di) if (lami(k).lt.lammini) then lami(k) = lammini n0i(k) = lami(k)**(di+1._r8)*dumi(i,k)/(ci*cons1) niic(i,k) = n0i(k)/lami(k) ! adjust number conc if needed to keep mean size in reasonable range if (do_cldice) nitend(i,k)=(niic(i,k)*icldm(i,k)-ni(i,k))/deltat else if (lami(k).gt.lammaxi) then lami(k) = lammaxi n0i(k) = lami(k)**(di+1._r8)*dumi(i,k)/(ci*cons1) niic(i,k) = n0i(k)/lami(k) ! adjust number conc if needed to keep mean size in reasonable range if (do_cldice) nitend(i,k)=(niic(i,k)*icldm(i,k)-ni(i,k))/deltat end if effi(i,k) = 1.5_r8/lami(k)*1.e6_r8 else effi(i,k) = 25._r8 end if ! NOTE: If CARMA is doing the ice microphysics, then the ice effective ! radius has already been determined from the size distribution. if (.not. do_cldice) then effi(i,k) = re_ice(i,k) * 1e6_r8 ! m -> um end if !................... ! cloud droplet effective radius if (dumc(i,k).ge.qsmall) then ! add upper limit to in-cloud number concentration to prevent numerical error dumnc(i,k)=min(dumnc(i,k),dumc(i,k)*1.e20_r8) !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ! set tendency to ensure minimum droplet concentration ! after update by microphysics, except when lambda exceeds bounds on mean drop ! size or if there is no cloud water if (dumnc(i,k).lt.cdnl/rho(i,k)) then nctend(i,k)=(cdnl/rho(i,k)*lcldm(i,k)-nc(i,k))/deltat end if dumnc(i,k)=max(dumnc(i,k),cdnl/rho(i,k)) ! sghan minimum in #/cm3 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8 pgam(k)=1._r8/(pgam(k)**2)-1._r8 pgam(k)=max(pgam(k),2._r8) pgam(k)=min(pgam(k),15._r8) lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ & (dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8) lammin = (pgam(k)+1._r8)/50.e-6_r8 ! Multiply by omsm to fit within RRTMG's table. lammax = (pgam(k)+1._r8)*omsm/2.e-6_r8 if (lamc(k).lt.lammin) then lamc(k) = lammin ncic(i,k) = 6._r8*lamc(k)**3*dumc(i,k)* & gamma(pgam(k)+1._r8)/ & (pi*rhow*gamma(pgam(k)+4._r8)) ! adjust number conc if needed to keep mean size in reasonable range nctend(i,k)=(ncic(i,k)*lcldm(i,k)-nc(i,k))/deltat else if (lamc(k).gt.lammax) then lamc(k) = lammax ncic(i,k) = 6._r8*lamc(k)**3*dumc(i,k)* & gamma(pgam(k)+1._r8)/ & (pi*rhow*gamma(pgam(k)+4._r8)) ! adjust number conc if needed to keep mean size in reasonable range nctend(i,k)=(ncic(i,k)*lcldm(i,k)-nc(i,k))/deltat end if effc(i,k) = & gamma(pgam(k)+4._r8)/ & gamma(pgam(k)+3._r8)/lamc(k)/2._r8*1.e6_r8 !assign output fields for shape here lamcrad(i,k)=lamc(k) pgamrad(i,k)=pgam(k) else effc(i,k) = 10._r8 lamcrad(i,k)=0._r8 pgamrad(i,k)=0._r8 end if ! ice effective diameter for david mitchell's optics if (do_cldice) then deffi(i,k)=effi(i,k)*rhoi/917._r8*2._r8 else deffi(i,k)=effi(i,k) * 2._r8 end if !!! recalculate effective radius for constant number, in order to separate ! first and second indirect effects ! assume constant number of 10^8 kg-1 dumnc(i,k)=1.e8_r8 if (dumc(i,k).ge.qsmall) then pgam(k)=0.0005714_r8*(ncic(i,k)/1.e6_r8*rho(i,k))+0.2714_r8 pgam(k)=1._r8/(pgam(k)**2)-1._r8 pgam(k)=max(pgam(k),2._r8) pgam(k)=min(pgam(k),15._r8) lamc(k) = (pi/6._r8*rhow*dumnc(i,k)*gamma(pgam(k)+4._r8)/ & (dumc(i,k)*gamma(pgam(k)+1._r8)))**(1._r8/3._r8) lammin = (pgam(k)+1._r8)/50.e-6_r8 lammax = (pgam(k)+1._r8)/2.e-6_r8 if (lamc(k).lt.lammin) then lamc(k) = lammin else if (lamc(k).gt.lammax) then lamc(k) = lammax end if effc_fn(i,k) = & gamma(pgam(k)+4._r8)/ & gamma(pgam(k)+3._r8)/lamc(k)/2._r8*1.e6_r8 else effc_fn(i,k) = 10._r8 end if !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!1! end do ! vertical k loop 500 continue do k=top_lev,pver ! if updated q (after microphysics) is zero, then ensure updated n is also zero if (qc(i,k)+qctend(i,k)*deltat.lt.qsmall) nctend(i,k)=-nc(i,k)/deltat if (do_cldice .and. qi(i,k)+qitend(i,k)*deltat.lt.qsmall) nitend(i,k)=-ni(i,k)/deltat end do end do ! i loop ! add snow ouptut do i = 1,ncol do k=top_lev,pver if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then dsout(i,k)=3._r8*rhosn/917._r8*(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8) endif end do end do !calculate effective radius of rain and snow in microns for COSP using Eq. 9 of COSP v1.3 manual do i = 1,ncol do k=top_lev,pver !! RAIN if (qrout(i,k).gt.1.e-7_r8.and.nrout(i,k).gt.0._r8) then reff_rain(i,k)=1.5_r8*(pi * rhow * nrout(i,k)/qrout(i,k))**(-1._r8/3._r8)*1.e6_r8 endif !! SNOW if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then reff_snow(i,k)=1.5_r8*(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8)*1.e6_r8 end if end do end do ! analytic radar reflectivity ! formulas from Matthew Shupe, NOAA/CERES ! *****note: radar reflectivity is local (in-precip average) ! units of mm^6/m^3 do i = 1,ncol do k=top_lev,pver if (qc(i,k)+qctend(i,k)*deltat.ge.qsmall .and. nc(i,k)+nctend(i,k)*deltat.gt.10._r8) then dum=((qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k)*rho(i,k)*1000._r8)**2 & /(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) else dum=0._r8 end if if (qi(i,k)+qitend(i,k)*deltat.ge.qsmall) then 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) else dum1=0._r8 end if if (qsout(i,k).ge.qsmall) then dum1=dum1+(qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8) end if refl(i,k)=dum+dum1 ! add rain rate, but for 37 GHz formulation instead of 94 GHz ! formula approximated from data of Matrasov (2007) ! rainrt is the rain rate in mm/hr ! reflectivity (dum) is in DBz if (rainrt(i,k).ge.0.001_r8) then dum=log10(rainrt(i,k)**6._r8)+16._r8 ! convert from DBz to mm^6/m^3 dum = 10._r8**(dum/10._r8) else ! don't include rain rate in R calculation for values less than 0.001 mm/hr dum=0._r8 end if ! add to refl refl(i,k)=refl(i,k)+dum !output reflectivity in Z. areflz(i,k)=refl(i,k) ! convert back to DBz if (refl(i,k).gt.minrefl) then refl(i,k)=10._r8*log10(refl(i,k)) else refl(i,k)=-9999._r8 end if !set averaging flag if (refl(i,k).gt.mindbz) then arefl(i,k)=refl(i,k) frefl(i,k)=1.0_r8 else arefl(i,k)=0._r8 areflz(i,k)=0._r8 frefl(i,k)=0._r8 end if ! bound cloudsat reflectivity csrfl(i,k)=min(csmax,refl(i,k)) !set averaging flag if (csrfl(i,k).gt.csmin) then acsrfl(i,k)=refl(i,k) fcsrfl(i,k)=1.0_r8 else acsrfl(i,k)=0._r8 fcsrfl(i,k)=0._r8 end if end do end do ! averaging for snow and rain number and diameter qrout2(:,:)=0._r8 qsout2(:,:)=0._r8 nrout2(:,:)=0._r8 nsout2(:,:)=0._r8 drout2(:,:)=0._r8 dsout2(:,:)=0._r8 freqs(:,:)=0._r8 freqr(:,:)=0._r8 do i = 1,ncol do k=top_lev,pver if (qrout(i,k).gt.1.e-7_r8.and.nrout(i,k).gt.0._r8) then qrout2(i,k)=qrout(i,k) nrout2(i,k)=nrout(i,k) drout2(i,k)=(pi * rhow * nrout(i,k)/qrout(i,k))**(-1._r8/3._r8) freqr(i,k)=1._r8 endif if (qsout(i,k).gt.1.e-7_r8.and.nsout(i,k).gt.0._r8) then qsout2(i,k)=qsout(i,k) nsout2(i,k)=nsout(i,k) dsout2(i,k)=(pi * rhosn * nsout(i,k)/qsout(i,k))**(-1._r8/3._r8) freqs(i,k)=1._r8 endif end do end do ! output activated liquid and ice (convert from #/kg -> #/m3) do i = 1,ncol do k=top_lev,pver ncai(i,k)=dum2i(i,k)*rho(i,k) ncal(i,k)=dum2l(i,k)*rho(i,k) end do end do !redefine fice here.... nfice(:,:)=0._r8 do k=top_lev,pver do i=1,ncol dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat) dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat) dumfice=qsout(i,k) + qrout(i,k) + dumc(i,k) + dumi(i,k) if (dumfice.gt.qsmall.and.(qsout(i,k)+dumi(i,k).gt.qsmall)) then nfice(i,k)=(qsout(i,k) + dumi(i,k))/dumfice endif if (nfice(i,k).gt.1._r8) then nfice(i,k)=1._r8 endif enddo enddo end subroutine micro_mg_tend !======================================================================== !UTILITIES !======================================================================== pure subroutine micro_mg_get_cols(ncol, nlev, top_lev, qcn, qin, & mgncol, mgcols) ! Determines which columns microphysics should operate over by ! checking for non-zero cloud water/ice. integer, intent(in) :: ncol ! Number of columns with meaningful data integer, intent(in) :: nlev ! Number of levels to use integer, intent(in) :: top_lev ! Top level for microphysics real(r8), intent(in) :: qcn(:,:) ! cloud water mixing ratio (kg/kg) real(r8), intent(in) :: qin(:,:) ! cloud ice mixing ratio (kg/kg) integer, intent(out) :: mgncol ! Number of columns MG will use integer, allocatable, intent(out) :: mgcols(:) ! column indices integer :: lev_offset ! top_lev - 1 (defined here for consistency) logical :: ltrue(ncol) ! store tests for each column integer :: i, ii ! column indices if (allocated(mgcols)) deallocate(mgcols) lev_offset = top_lev - 1 ! Using "any" along dimension 2 collapses across levels, but ! not columns, so we know if water is present at any level ! in each column. ltrue = any(qcn(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2) ltrue = ltrue .or. any(qin(:ncol,top_lev:(nlev+lev_offset)) >= qsmall, 2) ! Scan for true values to get a usable list of indices. mgncol = count(ltrue) allocate(mgcols(mgncol)) i = 0 do ii = 1,ncol if (ltrue(ii)) then i = i + 1 mgcols(i) = ii end if end do end subroutine micro_mg_get_cols end module micro_mg1_0