#define MODAL_AERO #define WRF_PORT ! Updated to CESM1.0.3 (CAM5.1.01) by Balwinder.Singh@pnnl.gov module cldwat2m_micro !--------------------------------------------------------------------------------- ! Purpose: ! CAM Interface for 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 !--------------------------------------------------------------------------------- ! modification for sub-columns, HM, (orig 8/11/10) ! This is done using the logical 'sub_column' set to .true. = subcolumns !--------------------------------------------------------------------------------- use shr_kind_mod, only: r8=>shr_kind_r8 #ifndef WRF_PORT use spmd_utils, only: masterproc use ppgrid, only: pcols, pver, pverp #else use module_cam_support, only: masterproc, pcols, pver, pverp #endif use physconst, only: gravit, rair, tmelt, cpair, rh2o, rhoh2o use physconst, only: latvap, latice use wv_saturation, only: polysvp, epsqs #ifndef WRF_PORT use cam_history, only: addfld, add_default, phys_decomp, outfld use cam_history_support, only : fillvalue use cam_logfile, only: iulog use phys_control, only: phys_getopts use cldwat2m_macro, only: rhmini #else use module_cam_support, only: addfld, add_default, phys_decomp, outfld, & fillvalue, iulog #endif implicit none private save #ifdef WRF_PORT !Macrophysics is not ported yet therefore rhmini is defined here (obtained from Cam5's cldwat2_macro)-Balwinder.Singh@pnnl.gov real(r8), public, parameter :: rhmini = 0.80_r8 ! Minimum rh for ice cloud fraction > 0. #endif logical, public :: liu_in = .true. ! True = Liu et al 2007 Ice nucleation ! False = cooper fixed ice nucleation (MG2008) public :: ini_micro, mmicro_pcond,gamma !constants remaped real(r8), private:: g !gravity real(r8), private:: r !Dry air Gas constant real(r8), private:: rv !water vapor gas contstant real(r8), private:: cpp !specific heat of dry air real(r8), private:: rhow !density of liquid water real(r8), private:: xxlv ! latent heat of vaporization real(r8), private:: xlf !latent heat of freezing real(r8), private:: xxls !latent heat of sublimation real(r8), private:: rhosn ! bulk density snow real(r8), private:: rhoi ! bulk density ice real(r8), private:: ac,bc,as,bs,ai,bi,ar,br !fall speed parameters real(r8), private:: ci,di !ice mass-diameter relation parameters real(r8), private:: cs,ds !snow mass-diameter relation parameters real(r8), private:: cr,dr !drop mass-diameter relation parameters real(r8), private:: f1s,f2s !ventilation param for snow real(r8), private:: Eii !collection efficiency aggregation of ice real(r8), private:: Ecc !collection efficiency real(r8), private:: Ecr !collection efficiency cloud droplets/rain real(r8), private:: f1r,f2r !ventilation param for rain real(r8), private:: DCS !autoconversion size threshold real(r8), private:: qsmall !min mixing ratio real(r8), private:: bimm,aimm !immersion freezing real(r8), private:: rhosu !typical 850mn air density real(r8), private:: mi0 ! new crystal mass real(r8), private:: rin ! radius of contact nuclei real(r8), private:: qcvar ! 1/relative variance of sub-grid qc real(r8), private:: pi ! pi ! Additional constants to help speed up code real(r8), private:: cons1 real(r8), private:: cons2 real(r8), private:: cons3 real(r8), private:: cons4 real(r8), private:: cons5 real(r8), private:: cons6 real(r8), private:: cons7 real(r8), private:: cons8 real(r8), private:: cons9 real(r8), private:: cons10 real(r8), private:: cons11 real(r8), private:: cons12 real(r8), private:: cons13 real(r8), private:: cons14 real(r8), private:: cons15 real(r8), private:: cons16 real(r8), private:: cons17 real(r8), private:: cons18 real(r8), private:: cons19 real(r8), private:: cons20 real(r8), private:: cons21 real(r8), private:: cons22 real(r8), private:: cons23 real(r8), private:: cons24 real(r8), private:: cons25 real(r8), private:: cons27 real(r8), private:: cons28 real(r8), private:: lammini real(r8), private:: lammaxi real(r8), private:: lamminr real(r8), private:: lammaxr real(r8), private:: lammins real(r8), private:: lammaxs ! parameters for snow/rain fraction for convective clouds real(r8), private:: tmax_fsnow ! max temperature for transition to convective snow real(r8), private:: tmin_fsnow ! min temperature for transition to convective snow !needed for findsp real(r8), private:: tt0 ! Freezing temperature real(r8), private:: csmin,csmax,minrefl,mindbz contains !=============================================================================== subroutine ini_micro !----------------------------------------------------------------------- ! ! Purpose: ! initialize constants for the morrison microphysics ! called from stratiform.F90 ! ! Author: Andrew Gettelman Dec 2005 ! !----------------------------------------------------------------------- integer k integer l,m, iaer real(r8) surften ! surface tension of water w/respect to air (N/m) real(r8) arg character(len=16) :: eddy_scheme = ' ' logical :: history_microphysics ! output variables for microphysics diagnostics package #ifndef WRF_PORT ! Query the PBL eddy scheme call phys_getopts(eddy_scheme_out = eddy_scheme, & history_microphysics_out = history_microphysics ) #else history_microphysics = .FALSE. #endif ! diagnostic precip call addfld ('QRAIN ','kg/kg ',pver, 'A','Diagnostic grid-mean rain mixing ratio' ,phys_decomp) call addfld ('QSNOW ','kg/kg ',pver, 'A','Diagnostic grid-mean snow mixing ratio' ,phys_decomp) call addfld ('NRAIN ','m-3 ',pver, 'A','Diagnostic grid-mean rain number conc' ,phys_decomp) call addfld ('NSNOW ','m-3 ',pver, 'A','Diagnostic grid-mean snow number conc' ,phys_decomp) ! size of precip call addfld ('RERCLD ','m ',pver, 'A','Diagnostic effective radius of Liquid Cloud and Rain' ,phys_decomp) call addfld ('DSNOW ','m ',pver, 'A','Diagnostic grid-mean snow diameter' ,phys_decomp) ! diagnostic radar reflectivity, cloud-averaged call addfld ('REFL ','DBz ',pver, 'A','94 GHz radar reflectivity' ,phys_decomp) call addfld ('AREFL ','DBz ',pver, 'A','Average 94 GHz radar reflectivity' ,phys_decomp) call addfld ('FREFL ','fraction ',pver, 'A','Fractional occurance of radar reflectivity' ,phys_decomp) call addfld ('CSRFL ','DBz ',pver, 'A','94 GHz radar reflectivity (CloudSat thresholds)' ,phys_decomp) call addfld ('ACSRFL ','DBz ',pver, 'A','Average 94 GHz radar reflectivity (CloudSat thresholds)' ,phys_decomp) call addfld ('FCSRFL ','fraction ',pver, 'A','Fractional occurance of radar reflectivity (CloudSat thresholds)' & ,phys_decomp) call addfld ('AREFLZ ','mm^6/m^3 ',pver, 'A','Average 94 GHz radar reflectivity' ,phys_decomp) ! Aerosol information call addfld ('NCAL ','#/m3 ',pver, 'A','Number Concentation Activated for Liquid',phys_decomp) call addfld ('NCAI ','#/m3 ',pver, 'A','Number Concentation Activated for Ice',phys_decomp) ! Average rain and snow mixing ratio (Q), number (N) and diameter (D), with frequency call addfld ('AQRAIN ','kg/kg ',pver, 'A','Average rain mixing ratio' ,phys_decomp) call addfld ('AQSNOW ','kg/kg ',pver, 'A','Average snow mixing ratio' ,phys_decomp) call addfld ('ANRAIN ','m-3 ',pver, 'A','Average rain number conc' ,phys_decomp) call addfld ('ANSNOW ','m-3 ',pver, 'A','Average snow number conc' ,phys_decomp) call addfld ('ADRAIN ','Micron ',pver, 'A','Average rain effective Diameter' ,phys_decomp) call addfld ('ADSNOW ','Micron ',pver, 'A','Average snow effective Diameter' ,phys_decomp) call addfld ('FREQR ','fraction ',pver, 'A','Fractional occurance of rain' ,phys_decomp) call addfld ('FREQS ','fraction ',pver, 'A','Fractional occurance of snow' ,phys_decomp) if ( history_microphysics) then call add_default ('AQSNOW ', 1, ' ') call add_default ('FREQR ', 1, ' ') call add_default ('FREQS ', 1, ' ') call add_default ('AQRAIN ', 1, ' ') call add_default ('AQSNOW ', 1, ' ') call add_default ('ANRAIN ', 1, ' ') call add_default ('ANSNOW ', 1, ' ') 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 ! latent heats xxlv = latvap ! latent heat vaporization xlf = latice ! latent heat freezing xxls = xxlv + xlf ! latent heat of sublimation ! 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 = 400.e-6_r8 ! 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 ! 1 / relative variance of sub-grid cloud water distribution ! see morrison and gettelman, 2007, J. Climate for details qcvar = 2._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) cons2=gamma(qcvar+2.47_r8) cons3=gamma(qcvar) 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) cons9=gamma(qcvar+2._r8) cons10=gamma(qcvar+1._r8) cons11=gamma(3._r8+bs) cons12=gamma(qcvar+1.15_r8) cons13=gamma(5._r8/2._r8+br/2._r8) cons14=gamma(5._r8/2._r8+bs/2._r8) cons15=gamma(qcvar+bc/3._r8) cons16=gamma(1._r8+bi) cons17=gamma(4._r8+bi) cons18=qcvar**2.47_r8 cons19=qcvar**2 cons20=qcvar**1.15_r8 cons21=qcvar**(bc/3._r8) 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 return end subroutine ini_micro !=============================================================================== !microphysics routine for each timestep goes here... subroutine mmicro_pcond ( sub_column, & lchnk, ncol, deltatin, tn, & qn, qc, qi, & nc, ni, p, pdel, cldn, & liqcldf, icecldf, & cldo, & rate1ord_cw2pr_st, & naai, npccnin, rndst,nacon, & tlat, qvlat, & qctend, qitend, nctend, nitend, effc, & effc_fn, effi, prect, preci, & nevapr, evapsnow, & 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 & #ifdef WRF_PORT , nsout, nrout & #endif ) !Author: Hugh Morrison, Andrew Gettelman, NCAR ! e-mail: morrison@ucar.edu, andrew@ucar.edu use wv_saturation, only: vqsatd_water ! input arguments logical, intent(in) :: sub_column ! True = configure for sub-columns False = use w/o sub-columns (standard) integer, intent(in) :: lchnk integer, intent(in) :: ncol 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) ! 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(inout) :: cldo(pcols,pver) ! old 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) ! 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) :: 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) :: 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) :: 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 ! local workspace ! all units mks unless otherwise stated ! 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) :: rercld(pcols,pver) ! effective radius calculation for rain + cloud 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 !fice variable real(r8) :: nfice(pcols,pver) !! 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) real(r8) :: gammas(pcols) ! parameter for cond/evap of cloud water ! 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), intent(out) :: qrout(pcols,pver) ! grid-box average rain mixing ratio (kg/kg) real(r8), intent(out) :: reff_rain(pcols,pver) ! rain effective radius (micron) real(r8), intent(out) :: reff_snow(pcols,pver) ! snow effective radius (micron) real(r8) :: drout(pcols,pver) ! rain diameter (m) #ifndef WRF_PORT real(r8) :: nrout(pcols,pver) ! rain number concentration (1/m3) real(r8) :: nsout(pcols,pver) ! snow number concentration (1/m3) #else !WRF needs nrout and nsout to be the ouput variables real(r8), intent(out) :: nrout(pcols,pver) ! rain number concentration (1/m3) real(r8), intent(out) :: nsout(pcols,pver) ! snow number concentration (1/m3) #endif real(r8), intent(out) :: dsout(pcols,pver) ! snow diameter (m) real(r8), intent(out) :: qsout(pcols,pver) ! snow mixing ratio (kg/kg) !averageed rain/snow for history real(r8) :: qrout2(pcols,pver) real(r8) :: qsout2(pcols,pver) real(r8) :: nrout2(pcols,pver) real(r8) :: nsout2(pcols,pver) real(r8) :: freqs(pcols,pver) real(r8) :: freqr(pcols,pver) real(r8) :: dumfice real(r8) :: drout2(pcols,pver) ! mean rain particle diameter (m) real(r8) :: dsout2(pcols,pver) ! mean snow particle diameter (m) !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 !output fields for number conc real(r8) :: ncai(pcols,pver) ! output number conc of ice nuclei available (1/m3) real(r8) :: ncal(pcols,pver) ! output number conc of CCN (1/m3) ! 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) :: refl(pcols,pver) ! analytic radar reflectivity real(r8) :: rainrt(pcols,pver) ! rain rate for reflectivity calculation real(r8) :: rainrt1(pcols,pver) real(r8) :: csrfl(pcols,pver) !cloudsat reflectivity real(r8) :: arefl(pcols,pver) !average reflectivity will zero points outside valid range real(r8) :: acsrfl(pcols,pver) !cloudsat average real(r8) :: frefl(pcols,pver) real(r8) :: fcsrfl(pcols,pver) real(r8) :: areflz(pcols,pver) !average reflectivity in z. real(r8) :: tmp real(r8) dmc,ssmc,dstrn ! variables for modal scheme. real(r8), parameter :: cdnl = 0.e6_r8 ! cloud droplet number limiter !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc ! 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 ! 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 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 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 !! initialize as fillvalue to avoid Floating Exceptions reff_rain(1:ncol,1:pver)=fillvalue reff_snow(1:ncol,1:pver)=fillvalue ! initialize variables for trop_mozart nevapr(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 ! 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=1,pver ! find wet bulk temperature and saturation value for provisional t and q without ! condensation call vqsatd_water(t(1,k),p(1,k),es,qs,gammas,ncol) ! use rhw do i=1,ncol esl(i,k)=polysvp(t(i,k),0) esi(i,k)=polysvp(t(i,k),1) ! 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 (sub_column) 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. cldm(i,k)=1. cldmw(i,k)=1. end if if (qi(i,k).ge.qsmall) then cldm(i,k)=1. icldm(i,k)=1. 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 (t(i,k).lt.tmelt - 5._r8) then if (liu_in) then ! if aerosols interact with ice set number of activated ice nuclei dum2=naai(i,k) else ! cooper curve (factor of 1000 is to convert from L-1 to m-3) dum2=0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k)))*1000._r8 ! put limit on number of nucleated crystals, set to number at T=-30 C ! cooper (limit to value at -35 C) dum2=min(dum,208.9e3_r8)/rho(i,k) ! convert from m-3 to kg-1 endif 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 (t(i,k).lt.273.15) then !if ice exists if (qi(i,k).gt.qsmall) then bergtsf = 0._r8 ! bergeron time scale (fraction of timestep) qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k)) qvl=epsqs*esl(i,k)/(p(i,k)-(1._r8-epsqs)*esl(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 ((relhum(i,k)*esl(i,k)/esi(i,k)).lt.1._r8 .and. qiic(i,k).ge.qsmall ) then qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k)) qvl=epsqs*esl(i,k)/(p(i,k)-(1._r8-epsqs)*esl(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 (t(i,k).lt.(tmelt - 5._r8)) then if ( liu_in ) 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 ! cooper curve (factor of 1000 is to convert from L-1 to m-3) dum2i(i,k)=0.005_r8*exp(0.304_r8*(273.15_r8-t(i,k)))*1000._r8 ! put limit on number of nucleated crystals, set to number at T=-30 C ! cooper (limit to value at -35 C) dum2i(i,k)=min(dum2i(i,k),208.9e3_r8)/rho(i,k) ! convert from m-3 to kg-1 endif else dum2i(i,k)=0._r8 end if end do ! i loop end do ! k loop cldo(:ncol,:)=cldn(:ncol,:) !! 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=1,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=1,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=1,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 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 do k=1,pver ! 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.1) 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 (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 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 (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=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(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 (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 (sub_column) 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.1) 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 (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 ! 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.1) 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 (qcic(i,k).ge.qsmall .and. t(i,k).lt.269.15_r8) then ! immersion freezing (Bigg, 1953) ! subcolumns if (sub_column) 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> Vi ! and Ds >> Di for continuous collection if (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=polysvp(t(i,k),0) qsn=min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8) ! recalculate saturation vapor pressure for liquid and ice esl(i,k)=esn esi(i,k)=polysvp(t(i,k),1) ! 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=epsqs*esl(i,k)/(p(i,k)-(1._r8-epsqs)*esl(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) end if ! sublimation of snow if (qniic(i,k).ge.qsmall) then qvi=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(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) 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=polysvp(ttmp,0) ! use rhw to allow ice supersaturation qsn=min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8) ! modify precip evaporation rate if q > qsat if (qtmp.gt.qsn) then if (pre(k)+prds(k).lt.-1.e-20) 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=polysvp(ttmp,0) ! use rhw to allow ice supersaturation qsn=min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8) 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=polysvp(ttmp,1) ! use rhi to allow ice supersaturation qsn=min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8) 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=epsqs*esi(i,k)/(p(i,k)-(1._r8-epsqs)*esi(i,k)) qvs=epsqs*esl(i,k)/(p(i,k)-(1._r8-epsqs)*esl(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 !PMA if (abs(prc(k)+pra(k)+mnuccc(k)+mnucct(k)+msacwi(k)+psacws(k)+bergs(k)) > 1.e-20_r8) then ratio = qce/deltat/lcldm(i,k)/(prc(k)+pra(k)+mnuccc(k)+mnucct(k)+msacwi(k)+psacws(k)+bergs(k))*omsm else ratio = 0.5_r8 end if 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 !PMA if (abs(nprc1(k)+npra(k)+nnuccc(k)+nnucct(k)+npsacws(k)-nsubc(k)) > 1.e-20_r8) then ratio = nce/deltat/((nprc1(k)+npra(k)+nnuccc(k)+nnucct(k)+& npsacws(k)-nsubc(k))*lcldm(i,k))*omsm else ratio = 0.5_r8 end if 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 dum = ((-mnuccc(k)-mnucct(k)-msacwi(k))*lcldm(i,k)+(prci(k)+ & prai(k))*icldm(i,k))*deltat if (dum.gt.qie) then !jdf if(abs(prci(k)+prai(k)).gt.1.0e-20) then ratio = (qie/deltat+(mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k))/((prci(k)+prai(k))*icldm(i,k))*omsm else ratio=0.5 endif !jdf ratio = (qie/deltat+(mnuccc(k)+mnucct(k)+msacwi(k))*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 dum = ((-nnucct(k)-nsacwi(k))*lcldm(i,k)+(nprci(k)+ & nprai(k)-nsubi(k))*icldm(i,k))*deltat if (dum.gt.nie) then if(abs( ((nprci(k)+nprai(k)-nsubi(k))*icldm(i,k))*omsm ) .gt. 1.0e-20_r8) then ratio = (nie/deltat+(nnucct(k)+nsacwi(k))*lcldm(i,k))/ & ((nprci(k)+nprai(k)-nsubi(k))*icldm(i,k))*omsm else ratio = 0.5 endif nprci(k) = nprci(k)*ratio nprai(k) = nprai(k)*ratio nsubi(k) = nsubi(k)*ratio end if ! for preciptiation 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) qitend(i,k) = qitend(i,k)+ & (mnuccc(k)+mnucct(k)+msacwi(k))*lcldm(i,k)+(-prci(k)- & prai(k))*icldm(i,k)+cmei(i,k)+berg(i,k) 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) ! 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) nitend(i,k) = nitend(i,k)+ nnuccd(k)*mtime+ & (nnucct(k)+nsacwi(k))*lcldm(i,k)+(nsubi(k)-nprci(k)- & nprai(k))*icldm(i,k) 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 (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.1) 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.1) 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 ! 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 * 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 !calculate precip fluxes and adding them to summing sub-stepping variables !! flux is zero at top interface rflx(i,1)=0.0_r8 sflx(i,1)=0.0_r8 !! calculating the precip flux (kg/m2/s) as mixingratio(kg/kg)*airdensity(kg/m3)*massweightedfallspeed(m/s) rflx(i,k+1)=qrout(i,k)*rho(i,k)*umr(k) sflx(i,k+1)=qsout(i,k)*rho(i,k)*ums(k) !! 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 = 1, 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,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=1,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) 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) 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 = 1,pver falouti(k) = fi(k)*dumi(i,k) faloutni(k) = fni(k)*dumni(i,k) faloutc(k) = fc(k)*dumc(i,k) faloutnc(k) = fnc(k)*dumnc(i,k) end do ! top of model k = 1 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 = 2,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 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=1,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 (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 = polysvp(ttmp,0) ! use rhw to allow ice supersaturation qsn = min(epsqs*esn/(p(i,k)-(1._r8-epsqs)*esn),1._r8) 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 !............................................................................... ! 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 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 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 !................... ! 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 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*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 deffi(i,k)=effi(i,k)*rhoi/917._r8*2._r8 !!! 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 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=1,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 (qi(i,k)+qitend(i,k)*deltat.lt.qsmall) nitend(i,k)=-ni(i,k)/deltat end do end do ! i loop ! hm add rain/snow mixing ratio and number concentration as diagnostic call outfld('QRAIN',qrout, pcols, lchnk) call outfld('QSNOW',qsout, pcols, lchnk) call outfld('NRAIN',nrout, pcols, lchnk) call outfld('NSNOW',nsout, pcols, lchnk) ! add snow ouptut do i = 1,ncol do k=1,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 call outfld('DSNOW',dsout, pcols, lchnk) !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=1,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 endif 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=1,pver if (qc(i,k)+qctend(i,k)*deltat.ge.qsmall) 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) 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. 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 call outfld('REFL',refl, pcols, lchnk) call outfld('AREFL',arefl, pcols, lchnk) call outfld('AREFLZ',areflz, pcols, lchnk) call outfld('FREFL',frefl, pcols, lchnk) call outfld('CSRFL',csrfl, pcols, lchnk) call outfld('ACSRFL',acsrfl, pcols, lchnk) call outfld('FCSRFL',fcsrfl, pcols, lchnk) call outfld('RERCLD',rercld, pcols, lchnk) ! 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=1,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=1,pver ncai(i,k)=dum2i(i,k)*rho(i,k) ncal(i,k)=dum2l(i,k)*rho(i,k) end do end do call outfld('NCAL',ncal, pcols,lchnk) call outfld('NCAI',ncai, pcols,lchnk) !add averaged output fields. call outfld('AQRAIN',qrout2, pcols,lchnk) call outfld('AQSNOW',qsout2, pcols,lchnk) call outfld('ANRAIN',nrout2, pcols,lchnk) call outfld('ANSNOW',nsout2, pcols,lchnk) call outfld('ADRAIN',drout2, pcols,lchnk) call outfld('ADSNOW',dsout2, pcols,lchnk) call outfld('FREQR',freqr, pcols,lchnk) call outfld('FREQS',freqs, pcols,lchnk) !redefine fice here.... nfice(:,:)=0._r8 do k=1,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 call outfld('FICE',nfice, pcols, lchnk) return end subroutine mmicro_pcond ! !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc FUNCTION GAMMA(X) !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc !D DOUBLE PRECISION FUNCTION DGAMMA(X) !---------------------------------------------------------------------- ! ! THIS ROUTINE CALCULATES THE GAMMA FUNCTION FOR A REAL ARGUMENT X. ! COMPUTATION IS BASED ON AN ALGORITHM OUTLINED IN REFERENCE 1. ! THE PROGRAM USES RATIONAL FUNCTIONS THAT APPROXIMATE THE GAMMA ! FUNCTION TO AT LEAST 20 SIGNIFICANT DECIMAL DIGITS. COEFFICIENTS ! FOR THE APPROXIMATION OVER THE INTERVAL (1,2) ARE UNPUBLISHED. ! THOSE FOR THE APPROXIMATION FOR X .GE. 12 ARE FROM REFERENCE 2. ! THE ACCURACY ACHIEVED DEPENDS ON THE ARITHMETIC SYSTEM, THE ! COMPILER, THE INTRINSIC FUNCTIONS, AND PROPER SELECTION OF THE ! MACHINE-DEPENDENT CONSTANTS. ! ! !******************************************************************* !******************************************************************* ! ! EXPLANATION OF MACHINE-DEPENDENT CONSTANTS ! ! BETA - RADIX FOR THE FLOATING-POINT REPRESENTATION ! MAXEXP - THE SMALLEST POSITIVE POWER OF BETA THAT OVERFLOWS ! XBIG - THE LARGEST ARGUMENT FOR WHICH GAMMA(X) IS REPRESENTABLE ! IN THE MACHINE, I.E., THE SOLUTION TO THE EQUATION ! GAMMA(XBIG) = BETA**MAXEXP ! XINF - THE LARGEST MACHINE REPRESENTABLE FLOATING-POINT NUMBER; ! APPROXIMATELY BETA**MAXEXP ! EPS - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT ! 1.0+EPS .GT. 1.0 ! XMININ - THE SMALLEST POSITIVE FLOATING-POINT NUMBER SUCH THAT ! 1/XMININ IS MACHINE REPRESENTABLE ! ! APPROXIMATE VALUES FOR SOME IMPORTANT MACHINES ARE: ! ! BETA MAXEXP XBIG ! ! CRAY-1 (S.P.) 2 8191 966.961 ! CYBER 180/855 ! UNDER NOS (S.P.) 2 1070 177.803 ! IEEE (IBM/XT, ! SUN, ETC.) (S.P.) 2 128 35.040 ! IEEE (IBM/XT, ! SUN, ETC.) (D.P.) 2 1024 171.624 ! IBM 3033 (D.P.) 16 63 57.574 ! VAX D-FORMAT (D.P.) 2 127 34.844 ! VAX G-FORMAT (D.P.) 2 1023 171.489 ! ! XINF EPS XMININ ! ! CRAY-1 (S.P.) 5.45E+2465 7.11E-15 1.84E-2466 ! CYBER 180/855 ! UNDER NOS (S.P.) 1.26E+322 3.55E-15 3.14E-294 ! IEEE (IBM/XT, ! SUN, ETC.) (S.P.) 3.40E+38 1.19E-7 1.18E-38 ! IEEE (IBM/XT, ! SUN, ETC.) (D.P.) 1.79D+308 2.22D-16 2.23D-308 ! IBM 3033 (D.P.) 7.23D+75 2.22D-16 1.39D-76 ! VAX D-FORMAT (D.P.) 1.70D+38 1.39D-17 5.88D-39 ! VAX G-FORMAT (D.P.) 8.98D+307 1.11D-16 1.12D-308 ! !******************************************************************* !******************************************************************* ! ! ERROR RETURNS ! ! THE PROGRAM RETURNS THE VALUE XINF FOR SINGULARITIES OR ! WHEN OVERFLOW WOULD OCCUR. THE COMPUTATION IS BELIEVED ! TO BE FREE OF UNDERFLOW AND OVERFLOW. ! ! ! INTRINSIC FUNCTIONS REQUIRED ARE: ! ! INT, DBLE, EXP, LOG, REAL, SIN ! ! ! REFERENCES: AN OVERVIEW OF SOFTWARE DEVELOPMENT FOR SPECIAL ! FUNCTIONS W. J. CODY, LECTURE NOTES IN MATHEMATICS, ! 506, NUMERICAL ANALYSIS DUNDEE, 1975, G. A. WATSON ! (ED.), SPRINGER VERLAG, BERLIN, 1976. ! ! COMPUTER APPROXIMATIONS, HART, ET. AL., WILEY AND ! SONS, NEW YORK, 1968. ! ! LATEST MODIFICATION: OCTOBER 12, 1989 ! ! AUTHORS: W. J. CODY AND L. STOLTZ ! APPLIED MATHEMATICS DIVISION ! ARGONNE NATIONAL LABORATORY ! ARGONNE, IL 60439 ! !---------------------------------------------------------------------- INTEGER I,N LOGICAL PARITY real(r8) gamma REAL(r8) & !D DOUBLE PRECISION C,CONV,EPS,FACT,HALF,ONE,P,PI,Q,RES,SQRTPI,SUM,TWELVE, & TWO,X,XBIG,XDEN,XINF,XMININ,XNUM,Y,Y1,YSQ,Z,ZERO DIMENSION C(7),P(8),Q(8) !---------------------------------------------------------------------- ! MATHEMATICAL CONSTANTS !---------------------------------------------------------------------- DATA ONE,HALF,TWELVE,TWO,ZERO/1.0E0_r8,0.5E0_r8,12.0E0_r8,2.0E0_r8,0.0E0_r8/, & SQRTPI/0.9189385332046727417803297E0_r8/, & PI/3.1415926535897932384626434E0_r8/ !D DATA ONE,HALF,TWELVE,TWO,ZERO/1.0D0,0.5D0,12.0D0,2.0D0,0.0D0/, !D 1 SQRTPI/0.9189385332046727417803297D0/, !D 2 PI/3.1415926535897932384626434D0/ !---------------------------------------------------------------------- ! MACHINE DEPENDENT PARAMETERS !---------------------------------------------------------------------- DATA XBIG,XMININ,EPS/35.040E0_r8,1.18E-38_r8,1.19E-7_r8/, & XINF/3.4E38_r8/ !D DATA XBIG,XMININ,EPS/171.624D0,2.23D-308,2.22D-16/, !D 1 XINF/1.79D308/ !---------------------------------------------------------------------- ! NUMERATOR AND DENOMINATOR COEFFICIENTS FOR RATIONAL MINIMAX ! APPROXIMATION OVER (1,2). !---------------------------------------------------------------------- DATA P/-1.71618513886549492533811E+0_r8,2.47656508055759199108314E+1_r8,& -3.79804256470945635097577E+2_r8,6.29331155312818442661052E+2_r8,& 8.66966202790413211295064E+2_r8,-3.14512729688483675254357E+4_r8,& -3.61444134186911729807069E+4_r8,6.64561438202405440627855E+4_r8/ DATA Q/-3.08402300119738975254353E+1_r8,3.15350626979604161529144E+2_r8,& -1.01515636749021914166146E+3_r8,-3.10777167157231109440444E+3_r8,& 2.25381184209801510330112E+4_r8,4.75584627752788110767815E+3_r8,& -1.34659959864969306392456E+5_r8,-1.15132259675553483497211E+5_r8/ !D DATA P/-1.71618513886549492533811D+0,2.47656508055759199108314D+1, !D 1 -3.79804256470945635097577D+2,6.29331155312818442661052D+2, !D 2 8.66966202790413211295064D+2,-3.14512729688483675254357D+4, !D 3 -3.61444134186911729807069D+4,6.64561438202405440627855D+4/ !D DATA Q/-3.08402300119738975254353D+1,3.15350626979604161529144D+2, !D 1 -1.01515636749021914166146D+3,-3.10777167157231109440444D+3, !D 2 2.25381184209801510330112D+4,4.75584627752788110767815D+3, !D 3 -1.34659959864969306392456D+5,-1.15132259675553483497211D+5/ !---------------------------------------------------------------------- ! COEFFICIENTS FOR MINIMAX APPROXIMATION OVER (12, INF). !---------------------------------------------------------------------- DATA C/-1.910444077728E-03_r8,8.4171387781295E-04_r8, & -5.952379913043012E-04_r8,7.93650793500350248E-04_r8,& -2.777777777777681622553E-03_r8,8.333333333333333331554247E-02_r8,& 5.7083835261E-03_r8/ !D DATA C/-1.910444077728D-03,8.4171387781295D-04, !D 1 -5.952379913043012D-04,7.93650793500350248D-04, !D 2 -2.777777777777681622553D-03,8.333333333333333331554247D-02, !D 3 5.7083835261D-03/ !---------------------------------------------------------------------- ! STATEMENT FUNCTIONS FOR CONVERSION BETWEEN INTEGER AND FLOAT !---------------------------------------------------------------------- CONV(I) = REAL(I,r8) !D CONV(I) = DBLE(I) PARITY=.FALSE. FACT=ONE N=0 Y=X IF(Y.LE.ZERO)THEN !---------------------------------------------------------------------- ! ARGUMENT IS NEGATIVE !---------------------------------------------------------------------- Y=-X Y1=AINT(Y) RES=Y-Y1 IF(RES.NE.ZERO)THEN IF(Y1.NE.AINT(Y1*HALF)*TWO)PARITY=.TRUE. FACT=-PI/SIN(PI*RES) Y=Y+ONE ELSE RES=XINF GOTO 900 ENDIF ENDIF !---------------------------------------------------------------------- ! ARGUMENT IS POSITIVE !---------------------------------------------------------------------- IF(Y.LT.EPS)THEN !---------------------------------------------------------------------- ! ARGUMENT .LT. EPS !---------------------------------------------------------------------- IF(Y.GE.XMININ)THEN RES=ONE/Y ELSE RES=XINF GOTO 900 ENDIF ELSEIF(Y.LT.TWELVE)THEN Y1=Y IF(Y.LT.ONE)THEN !---------------------------------------------------------------------- ! 0.0 .LT. ARGUMENT .LT. 1.0 !---------------------------------------------------------------------- Z=Y Y=Y+ONE ELSE !---------------------------------------------------------------------- ! 1.0 .LT. ARGUMENT .LT. 12.0, REDUCE ARGUMENT IF NECESSARY !---------------------------------------------------------------------- N=INT(Y)-1 Y=Y-CONV(N) Z=Y-ONE ENDIF !---------------------------------------------------------------------- ! EVALUATE APPROXIMATION FOR 1.0 .LT. ARGUMENT .LT. 2.0 !---------------------------------------------------------------------- XNUM=ZERO XDEN=ONE DO 260 I=1,8 XNUM=(XNUM+P(I))*Z XDEN=XDEN*Z+Q(I) 260 CONTINUE RES=XNUM/XDEN+ONE IF(Y1.LT.Y)THEN !---------------------------------------------------------------------- ! ADJUST RESULT FOR CASE 0.0 .LT. ARGUMENT .LT. 1.0 !---------------------------------------------------------------------- RES=RES/Y1 ELSEIF(Y1.GT.Y)THEN !---------------------------------------------------------------------- ! ADJUST RESULT FOR CASE 2.0 .LT. ARGUMENT .LT. 12.0 !---------------------------------------------------------------------- DO 290 I=1,N RES=RES*Y Y=Y+ONE 290 CONTINUE ENDIF ELSE !---------------------------------------------------------------------- ! EVALUATE FOR ARGUMENT .GE. 12.0, !---------------------------------------------------------------------- IF(Y.LE.XBIG)THEN YSQ=Y*Y SUM=C(7) DO 350 I=1,6 SUM=SUM/YSQ+C(I) 350 CONTINUE SUM=SUM/Y-Y+SQRTPI SUM=SUM+(Y-HALF)*LOG(Y) RES=EXP(SUM) ELSE RES=XINF GOTO 900 ENDIF ENDIF !---------------------------------------------------------------------- ! FINAL ADJUSTMENTS AND RETURN !---------------------------------------------------------------------- IF(PARITY)RES=-RES IF(FACT.NE.ONE)RES=FACT/RES 900 GAMMA=RES !D900 DGAMMA = RES RETURN ! ---------- LAST LINE OF GAMMA ---------- END function gamma end module cldwat2m_micro