1 2! KGEN-generated Fortran source file 3! 4! Filename : micro_mg2_0.F90 5! Generated at: 2015-03-31 09:44:40 6! KGEN version: 0.4.5 7 8 9 10 MODULE micro_mg2_0 11 !--------------------------------------------------------------------------------- 12 ! Purpose: 13 ! MG microphysics version 2.0 - Update of MG microphysics with 14 ! prognostic precipitation. 15 ! 16 ! Author: Andrew Gettelman, Hugh Morrison. 17 ! Contributions from: Peter Caldwell, Xiaohong Liu and Steve Ghan 18 ! Version 2 history: Sep 2011: Development begun. 19 ! Feb 2013: Added of prognostic precipitation. 20 ! invoked in 1 by specifying -microphys=mg2.0 21 ! 22 ! for questions contact Hugh Morrison, Andrew Gettelman 23 ! e-mail: morrison@ucar.edu, andrew@ucar.edu 24 !--------------------------------------------------------------------------------- 25 ! 26 ! NOTE: Modified to allow other microphysics packages (e.g. CARMA) to do ice 27 ! microphysics in cooperation with the MG liquid microphysics. This is 28 ! controlled by the do_cldice variable. 29 ! 30 ! If do_cldice is false, then MG microphysics should not update CLDICE or 31 ! NUMICE; it is assumed that the other microphysics scheme will have updated 32 ! CLDICE and NUMICE. The other microphysics should handle the following 33 ! processes that would have been done by MG: 34 ! - Detrainment (liquid and ice) 35 ! - Homogeneous ice nucleation 36 ! - Heterogeneous ice nucleation 37 ! - Bergeron process 38 ! - Melting of ice 39 ! - Freezing of cloud drops 40 ! - Autoconversion (ice -> snow) 41 ! - Growth/Sublimation of ice 42 ! - Sedimentation of ice 43 ! 44 ! This option has not been updated since the introduction of prognostic 45 ! precipitation, and probably should be adjusted to cover snow as well. 46 ! 47 !--------------------------------------------------------------------------------- 48 ! Based on micro_mg (restructuring of former cldwat2m_micro) 49 ! Author: Andrew Gettelman, Hugh Morrison. 50 ! Contributions from: Xiaohong Liu and Steve Ghan 51 ! December 2005-May 2010 52 ! Description in: Morrison and Gettelman, 2008. J. Climate (MG2008) 53 ! Gettelman et al., 2010 J. Geophys. Res. - Atmospheres (G2010) 54 ! for questions contact Hugh Morrison, Andrew Gettelman 55 ! e-mail: morrison@ucar.edu, andrew@ucar.edu 56 !--------------------------------------------------------------------------------- 57 ! Code comments added by HM, 093011 58 ! General code structure: 59 ! 60 ! Code is divided into two main subroutines: 61 ! subroutine micro_mg_init --> initializes microphysics routine, should be called 62 ! once at start of simulation 63 ! subroutine micro_mg_tend --> main microphysics routine to be called each time step 64 ! this also calls several smaller subroutines to calculate 65 ! microphysical processes and other utilities 66 ! 67 ! List of external functions: 68 ! qsat_water --> for calculating saturation vapor pressure with respect to liquid water 69 ! qsat_ice --> for calculating saturation vapor pressure with respect to ice 70 ! gamma --> standard mathematical gamma function 71 ! ......................................................................... 72 ! List of inputs through use statement in fortran90: 73 ! Variable Name Description Units 74 ! ......................................................................... 75 ! gravit acceleration due to gravity m s-2 76 ! rair dry air gas constant for air J kg-1 K-1 77 ! tmelt temperature of melting point for water K 78 ! cpair specific heat at constant pressure for dry air J kg-1 K-1 79 ! rh2o gas constant for water vapor J kg-1 K-1 80 ! latvap latent heat of vaporization J kg-1 81 ! latice latent heat of fusion J kg-1 82 ! qsat_water external function for calculating liquid water 83 ! saturation vapor pressure/humidity - 84 ! qsat_ice external function for calculating ice 85 ! saturation vapor pressure/humidity pa 86 ! rhmini relative humidity threshold parameter for 87 ! nucleating ice - 88 ! ......................................................................... 89 ! NOTE: List of all inputs/outputs passed through the call/subroutine statement 90 ! for micro_mg_tend is given below at the start of subroutine micro_mg_tend. 91 !--------------------------------------------------------------------------------- 92 ! Procedures required: 93 ! 1) An implementation of the gamma function (if not intrinsic). 94 ! 2) saturation vapor pressure and specific humidity over water 95 ! 3) svp over ice 96 USE shr_spfn_mod, ONLY: gamma => shr_spfn_gamma 97 USE wv_sat_methods, ONLY: qsat_water => wv_sat_qsat_water 98 USE wv_sat_methods, ONLY: qsat_ice => wv_sat_qsat_ice 99 ! Parameters from the utilities module. 100 USE micro_mg_utils, ONLY: r8 101 USE micro_mg_utils, ONLY: qsmall 102 USE micro_mg_utils, ONLY: mincld 103 USE micro_mg_utils, ONLY: ar 104 USE micro_mg_utils, ONLY: as 105 USE micro_mg_utils, ONLY: rhow 106 USE micro_mg_utils, ONLY: ai 107 USE micro_mg_utils, ONLY: mi0 108 USE micro_mg_utils, ONLY: br 109 USE micro_mg_utils, ONLY: bs 110 USE micro_mg_utils, ONLY: pi 111 USE micro_mg_utils, ONLY: rhosn 112 USE micro_mg_utils, ONLY: omsm 113 USE micro_mg_utils, ONLY: rising_factorial 114 USE micro_mg_utils, ONLY: bc 115 USE micro_mg_utils, ONLY: bi 116 USE micro_mg_utils, ONLY: rhows 117 USE micro_mg_utils, ONLY: rhoi 118 IMPLICIT NONE 119 PRIVATE 120 PUBLIC micro_mg_tend 121 ! switch for specification rather than prediction of droplet and crystal number 122 ! note: number will be adjusted as needed to keep mean size within bounds, 123 ! even when specified droplet or ice number is used 124 ! If constant cloud ice number is set (nicons = .true.), 125 ! then all microphysical processes except mass transfer due to ice nucleation 126 ! (mnuccd) are based on the fixed cloud ice number. Calculation of 127 ! mnuccd follows from the prognosed ice crystal number ni. 128 ! nccons = .true. to specify constant cloud droplet number 129 ! nicons = .true. to specify constant cloud ice number 130 LOGICAL, parameter, public :: nccons = .false. 131 LOGICAL, parameter, public :: nicons = .false. 132 !========================================================= 133 ! Private module parameters 134 !========================================================= 135 ! parameters for specified ice and droplet number concentration 136 ! note: these are local in-cloud values, not grid-mean 137 REAL(KIND=r8), parameter :: ncnst = 100.e6_r8 ! droplet num concentration when nccons=.true. (m-3) 138 REAL(KIND=r8), parameter :: ninst = 0.1e6_r8 ! ice num concentration when nicons=.true. (m-3) 139 !Range of cloudsat reflectivities (dBz) for analytic simulator 140 REAL(KIND=r8), parameter :: csmin = -30._r8 141 REAL(KIND=r8), parameter :: csmax = 26._r8 142 REAL(KIND=r8), parameter :: mindbz = -99._r8 143 REAL(KIND=r8), parameter :: minrefl = 1.26e-10_r8 ! minrefl = 10._r8**(mindbz/10._r8) 144 ! autoconversion size threshold for cloud ice to snow (m) 145 REAL(KIND=r8) :: dcs 146 ! minimum mass of new crystal due to freezing of cloud droplets done 147 ! externally (kg) 148 REAL(KIND=r8), parameter :: mi0l_min = 4._r8/3._r8*pi*rhow*(4.e-6_r8)**3 149 !========================================================= 150 ! Constants set in initialization 151 !========================================================= 152 ! Set using arguments to micro_mg_init 153 REAL(KIND=r8) :: g ! gravity 154 REAL(KIND=r8) :: r ! dry air gas constant 155 REAL(KIND=r8) :: rv ! water vapor gas constant 156 REAL(KIND=r8) :: cpp ! specific heat of dry air 157 REAL(KIND=r8) :: tmelt ! freezing point of water (K) 158 ! latent heats of: 159 REAL(KIND=r8) :: xxlv ! vaporization 160 REAL(KIND=r8) :: xlf ! freezing 161 REAL(KIND=r8) :: xxls ! sublimation 162 REAL(KIND=r8) :: rhmini ! Minimum rh for ice cloud fraction > 0. 163 ! flags 164 LOGICAL :: microp_uniform 165 LOGICAL :: do_cldice 166 LOGICAL :: use_hetfrz_classnuc 167 REAL(KIND=r8) :: rhosu ! typical 850mn air density 168 REAL(KIND=r8) :: icenuct ! ice nucleation temperature: currently -5 degrees C 169 REAL(KIND=r8) :: snowmelt ! what temp to melt all snow: currently 2 degrees C 170 REAL(KIND=r8) :: rainfrze ! what temp to freeze all rain: currently -5 degrees C 171 ! additional constants to help speed up code 172 REAL(KIND=r8) :: gamma_br_plus1 173 REAL(KIND=r8) :: gamma_br_plus4 174 REAL(KIND=r8) :: gamma_bs_plus1 175 REAL(KIND=r8) :: gamma_bs_plus4 176 REAL(KIND=r8) :: gamma_bi_plus1 177 REAL(KIND=r8) :: gamma_bi_plus4 178 REAL(KIND=r8) :: xxlv_squared 179 REAL(KIND=r8) :: xxls_squared 180 CHARACTER(LEN=16) :: micro_mg_precip_frac_method ! type of precipitation fraction method 181 REAL(KIND=r8) :: micro_mg_berg_eff_factor ! berg efficiency factor 182 !=============================================================================== 183 PUBLIC kgen_read_externs_micro_mg2_0 184 CONTAINS 185 186 ! write subroutines 187 ! No subroutines 188 189 ! module extern variables 190 191 SUBROUTINE kgen_read_externs_micro_mg2_0(kgen_unit) 192 INTEGER, INTENT(IN) :: kgen_unit 193 READ(UNIT=kgen_unit) dcs 194 READ(UNIT=kgen_unit) g 195 READ(UNIT=kgen_unit) r 196 READ(UNIT=kgen_unit) rv 197 READ(UNIT=kgen_unit) cpp 198 READ(UNIT=kgen_unit) tmelt 199 READ(UNIT=kgen_unit) xxlv 200 READ(UNIT=kgen_unit) xlf 201 READ(UNIT=kgen_unit) xxls 202 READ(UNIT=kgen_unit) rhmini 203 READ(UNIT=kgen_unit) microp_uniform 204 READ(UNIT=kgen_unit) do_cldice 205 READ(UNIT=kgen_unit) use_hetfrz_classnuc 206 READ(UNIT=kgen_unit) rhosu 207 READ(UNIT=kgen_unit) icenuct 208 READ(UNIT=kgen_unit) snowmelt 209 READ(UNIT=kgen_unit) rainfrze 210 READ(UNIT=kgen_unit) gamma_br_plus1 211 READ(UNIT=kgen_unit) gamma_br_plus4 212 READ(UNIT=kgen_unit) gamma_bs_plus1 213 READ(UNIT=kgen_unit) gamma_bs_plus4 214 READ(UNIT=kgen_unit) gamma_bi_plus1 215 READ(UNIT=kgen_unit) gamma_bi_plus4 216 READ(UNIT=kgen_unit) xxlv_squared 217 READ(UNIT=kgen_unit) xxls_squared 218 READ(UNIT=kgen_unit) micro_mg_precip_frac_method 219 READ(UNIT=kgen_unit) micro_mg_berg_eff_factor 220 END SUBROUTINE kgen_read_externs_micro_mg2_0 221 222 !=============================================================================== 223 224 !=============================================================================== 225 !microphysics routine for each timestep goes here... 226 227 SUBROUTINE micro_mg_tend(mgncol, nlev, deltatin, t, q, qcn, qin, ncn, nin, qrn, qsn, nrn, nsn, relvar, accre_enhan, p, & 228 pdel, cldn, liqcldf, icecldf, qcsinksum_rate1ord, naai, npccn, rndst, nacon, tlat, qvlat, qctend, qitend, nctend, nitend, & 229 qrtend, qstend, nrtend, nstend, effc, effc_fn, effi, prect, preci, nevapr, evapsnow, prain, prodsnow, cmeout, deffi, & 230 pgamrad, lamcrad, qsout, dsout, rflx, sflx, qrout, reff_rain, reff_snow, qcsevap, qisevap, qvres, cmeitot, vtrmc, vtrmi, & 231 umr, ums, qcsedten, qisedten, qrsedten, qssedten, pratot, prctot, mnuccctot, mnuccttot, msacwitot, psacwstot, bergstot, & 232 bergtot, melttot, homotot, qcrestot, prcitot, praitot, qirestot, mnuccrtot, pracstot, meltsdttot, frzrdttot, mnuccdtot, & 233 nrout, nsout, refl, arefl, areflz, frefl, csrfl, acsrfl, fcsrfl, rercld, ncai, ncal, qrout2, qsout2, nrout2, nsout2, & 234 drout2, dsout2, freqs, freqr, nfice, qcrat, errstring, tnd_qsnow, tnd_nsnow, re_ice, prer_evap, frzimm, frzcnt, frzdep) 235 ! Below arguments are "optional" (pass null pointers to omit). 236 ! Constituent properties. 237 USE micro_mg_utils, ONLY: mg_liq_props 238 USE micro_mg_utils, ONLY: mg_ice_props 239 USE micro_mg_utils, ONLY: mg_rain_props 240 USE micro_mg_utils, ONLY: mg_snow_props 241 ! Size calculation functions. 242 USE micro_mg_utils, ONLY: size_dist_param_liq 243 USE micro_mg_utils, ONLY: size_dist_param_basic 244 USE micro_mg_utils, ONLY: avg_diameter 245 ! Microphysical processes. 246 USE micro_mg_utils, ONLY: kk2000_liq_autoconversion 247 USE micro_mg_utils, ONLY: ice_autoconversion 248 USE micro_mg_utils, ONLY: immersion_freezing 249 USE micro_mg_utils, ONLY: contact_freezing 250 USE micro_mg_utils, ONLY: snow_self_aggregation 251 USE micro_mg_utils, ONLY: accrete_cloud_water_snow 252 USE micro_mg_utils, ONLY: secondary_ice_production 253 USE micro_mg_utils, ONLY: accrete_rain_snow 254 USE micro_mg_utils, ONLY: heterogeneous_rain_freezing 255 USE micro_mg_utils, ONLY: accrete_cloud_water_rain 256 USE micro_mg_utils, ONLY: self_collection_rain 257 USE micro_mg_utils, ONLY: accrete_cloud_ice_snow 258 USE micro_mg_utils, ONLY: evaporate_sublimate_precip 259 USE micro_mg_utils, ONLY: bergeron_process_snow 260 USE micro_mg_utils, ONLY: ice_deposition_sublimation 261 !Authors: Hugh Morrison, Andrew Gettelman, NCAR, Peter Caldwell, LLNL 262 ! e-mail: morrison@ucar.edu, andrew@ucar.edu 263 ! input arguments 264 INTEGER, intent(in) :: mgncol ! number of microphysics columns 265 INTEGER, intent(in) :: nlev ! number of layers 266 REAL(KIND=r8), intent(in) :: deltatin ! time step (s) 267 REAL(KIND=r8), intent(in) :: t(:,:) ! input temperature (K) 268 REAL(KIND=r8), intent(in) :: q(:,:) ! input h20 vapor mixing ratio (kg/kg) 269 ! note: all input cloud variables are grid-averaged 270 REAL(KIND=r8), intent(in) :: qcn(:,:) ! cloud water mixing ratio (kg/kg) 271 REAL(KIND=r8), intent(in) :: qin(:,:) ! cloud ice mixing ratio (kg/kg) 272 REAL(KIND=r8), intent(in) :: ncn(:,:) ! cloud water number conc (1/kg) 273 REAL(KIND=r8), intent(in) :: nin(:,:) ! cloud ice number conc (1/kg) 274 REAL(KIND=r8), intent(in) :: qrn(:,:) ! rain mixing ratio (kg/kg) 275 REAL(KIND=r8), intent(in) :: qsn(:,:) ! snow mixing ratio (kg/kg) 276 REAL(KIND=r8), intent(in) :: nrn(:,:) ! rain number conc (1/kg) 277 REAL(KIND=r8), intent(in) :: nsn(:,:) ! snow number conc (1/kg) 278 REAL(KIND=r8), intent(in) :: relvar(:,:) ! cloud water relative variance (-) 279 REAL(KIND=r8), intent(in) :: accre_enhan(:,:) ! optional accretion 280 ! enhancement factor (-) 281 REAL(KIND=r8), intent(in) :: p(:,:) ! air pressure (pa) 282 REAL(KIND=r8), intent(in) :: pdel(:,:) ! pressure difference across level (pa) 283 REAL(KIND=r8), intent(in) :: cldn(:,:) ! cloud fraction (no units) 284 REAL(KIND=r8), intent(in) :: liqcldf(:,:) ! liquid cloud fraction (no units) 285 REAL(KIND=r8), intent(in) :: icecldf(:,:) ! ice cloud fraction (no units) 286 ! used for scavenging 287 ! Inputs for aerosol activation 288 REAL(KIND=r8), intent(in) :: naai(:,:) ! ice nucleation number (from microp_aero_ts) (1/kg) 289 REAL(KIND=r8), intent(in) :: npccn(:,:) ! ccn activated number tendency (from microp_aero_ts) (1/kg*s) 290 ! Note that for these variables, the dust bin is assumed to be the last index. 291 ! (For example, in 1, the last dimension is always size 4.) 292 REAL(KIND=r8), intent(in) :: rndst(:,:,:) ! radius of each dust bin, for contact freezing (from microp_aero_ts) (m) 293 REAL(KIND=r8), intent(in) :: nacon(:,:,:) ! number in each dust bin, for contact freezing (from microp_aero_ts) (1/m^3) 294 ! output arguments 295 REAL(KIND=r8), intent(out) :: qcsinksum_rate1ord(:,:) ! 1st order rate for 296 ! direct cw to precip conversion 297 REAL(KIND=r8), intent(out) :: tlat(:,:) ! latent heating rate (W/kg) 298 REAL(KIND=r8), intent(out) :: qvlat(:,:) ! microphysical tendency qv (1/s) 299 REAL(KIND=r8), intent(out) :: qctend(:,:) ! microphysical tendency qc (1/s) 300 REAL(KIND=r8), intent(out) :: qitend(:,:) ! microphysical tendency qi (1/s) 301 REAL(KIND=r8), intent(out) :: nctend(:,:) ! microphysical tendency nc (1/(kg*s)) 302 REAL(KIND=r8), intent(out) :: nitend(:,:) ! microphysical tendency ni (1/(kg*s)) 303 REAL(KIND=r8), intent(out) :: qrtend(:,:) ! microphysical tendency qr (1/s) 304 REAL(KIND=r8), intent(out) :: qstend(:,:) ! microphysical tendency qs (1/s) 305 REAL(KIND=r8), intent(out) :: nrtend(:,:) ! microphysical tendency nr (1/(kg*s)) 306 REAL(KIND=r8), intent(out) :: nstend(:,:) ! microphysical tendency ns (1/(kg*s)) 307 REAL(KIND=r8), intent(out) :: effc(:,:) ! droplet effective radius (micron) 308 REAL(KIND=r8), intent(out) :: effc_fn(:,:) ! droplet effective radius, assuming nc = 1.e8 kg-1 309 REAL(KIND=r8), intent(out) :: effi(:,:) ! cloud ice effective radius (micron) 310 REAL(KIND=r8), intent(out) :: prect(:) ! surface precip rate (m/s) 311 REAL(KIND=r8), intent(out) :: preci(:) ! cloud ice/snow precip rate (m/s) 312 REAL(KIND=r8), intent(out) :: nevapr(:,:) ! evaporation rate of rain + snow (1/s) 313 REAL(KIND=r8), intent(out) :: evapsnow(:,:) ! sublimation rate of snow (1/s) 314 REAL(KIND=r8), intent(out) :: prain(:,:) ! production of rain + snow (1/s) 315 REAL(KIND=r8), intent(out) :: prodsnow(:,:) ! production of snow (1/s) 316 REAL(KIND=r8), intent(out) :: cmeout(:,:) ! evap/sub of cloud (1/s) 317 REAL(KIND=r8), intent(out) :: deffi(:,:) ! ice effective diameter for optics (radiation) (micron) 318 REAL(KIND=r8), intent(out) :: pgamrad(:,:) ! ice gamma parameter for optics (radiation) (no units) 319 REAL(KIND=r8), intent(out) :: lamcrad(:,:) ! slope of droplet distribution for optics (radiation) (1/m) 320 REAL(KIND=r8), intent(out) :: qsout(:,:) ! snow mixing ratio (kg/kg) 321 REAL(KIND=r8), intent(out) :: dsout(:,:) ! snow diameter (m) 322 REAL(KIND=r8), intent(out) :: rflx(:,:) ! grid-box average rain flux (kg m^-2 s^-1) 323 REAL(KIND=r8), intent(out) :: sflx(:,:) ! grid-box average snow flux (kg m^-2 s^-1) 324 REAL(KIND=r8), intent(out) :: qrout(:,:) ! grid-box average rain mixing ratio (kg/kg) 325 REAL(KIND=r8), intent(out) :: reff_rain(:,:) ! rain effective radius (micron) 326 REAL(KIND=r8), intent(out) :: reff_snow(:,:) ! snow effective radius (micron) 327 REAL(KIND=r8), intent(out) :: qcsevap(:,:) ! cloud water evaporation due to sedimentation (1/s) 328 REAL(KIND=r8), intent(out) :: qisevap(:,:) ! cloud ice sublimation due to sublimation (1/s) 329 REAL(KIND=r8), intent(out) :: qvres(:,:) ! residual condensation term to ensure RH < 100% (1/s) 330 REAL(KIND=r8), intent(out) :: cmeitot(:,:) ! grid-mean cloud ice sub/dep (1/s) 331 REAL(KIND=r8), intent(out) :: vtrmc(:,:) ! mass-weighted cloud water fallspeed (m/s) 332 REAL(KIND=r8), intent(out) :: vtrmi(:,:) ! mass-weighted cloud ice fallspeed (m/s) 333 REAL(KIND=r8), intent(out) :: umr(:,:) ! mass weighted rain fallspeed (m/s) 334 REAL(KIND=r8), intent(out) :: ums(:,:) ! mass weighted snow fallspeed (m/s) 335 REAL(KIND=r8), intent(out) :: qcsedten(:,:) ! qc sedimentation tendency (1/s) 336 REAL(KIND=r8), intent(out) :: qisedten(:,:) ! qi sedimentation tendency (1/s) 337 REAL(KIND=r8), intent(out) :: qrsedten(:,:) ! qr sedimentation tendency (1/s) 338 REAL(KIND=r8), intent(out) :: qssedten(:,:) ! qs sedimentation tendency (1/s) 339 ! microphysical process rates for output (mixing ratio tendencies) (all have units of 1/s) 340 REAL(KIND=r8), intent(out) :: pratot(:,:) ! accretion of cloud by rain 341 REAL(KIND=r8), intent(out) :: prctot(:,:) ! autoconversion of cloud to rain 342 REAL(KIND=r8), intent(out) :: mnuccctot(:,:) ! mixing ratio tend due to immersion freezing 343 REAL(KIND=r8), intent(out) :: mnuccttot(:,:) ! mixing ratio tend due to contact freezing 344 REAL(KIND=r8), intent(out) :: msacwitot(:,:) ! mixing ratio tend due to H-M splintering 345 REAL(KIND=r8), intent(out) :: psacwstot(:,:) ! collection of cloud water by snow 346 REAL(KIND=r8), intent(out) :: bergstot(:,:) ! bergeron process on snow 347 REAL(KIND=r8), intent(out) :: bergtot(:,:) ! bergeron process on cloud ice 348 REAL(KIND=r8), intent(out) :: melttot(:,:) ! melting of cloud ice 349 REAL(KIND=r8), intent(out) :: homotot(:,:) ! homogeneous freezing cloud water 350 REAL(KIND=r8), intent(out) :: qcrestot(:,:) ! residual cloud condensation due to removal of excess supersat 351 REAL(KIND=r8), intent(out) :: prcitot(:,:) ! autoconversion of cloud ice to snow 352 REAL(KIND=r8), intent(out) :: praitot(:,:) ! accretion of cloud ice by snow 353 REAL(KIND=r8), intent(out) :: qirestot(:,:) ! residual ice deposition due to removal of excess supersat 354 REAL(KIND=r8), intent(out) :: mnuccrtot(:,:) ! mixing ratio tendency due to heterogeneous freezing of rain to snow (1/s) 355 REAL(KIND=r8), intent(out) :: pracstot(:,:) ! mixing ratio tendency due to accretion of rain by snow (1/s) 356 REAL(KIND=r8), intent(out) :: meltsdttot(:,:) ! latent heating rate due to melting of snow (W/kg) 357 REAL(KIND=r8), intent(out) :: frzrdttot(:,:) ! latent heating rate due to homogeneous freezing of rain (W/kg) 358 REAL(KIND=r8), intent(out) :: mnuccdtot(:,:) ! mass tendency from ice nucleation 359 REAL(KIND=r8), intent(out) :: nrout(:,:) ! rain number concentration (1/m3) 360 REAL(KIND=r8), intent(out) :: nsout(:,:) ! snow number concentration (1/m3) 361 REAL(KIND=r8), intent(out) :: refl(:,:) ! analytic radar reflectivity 362 REAL(KIND=r8), intent(out) :: arefl(:,:) ! average reflectivity will zero points outside valid range 363 REAL(KIND=r8), intent(out) :: areflz(:,:) ! average reflectivity in z. 364 REAL(KIND=r8), intent(out) :: frefl(:,:) ! fractional occurrence of radar reflectivity 365 REAL(KIND=r8), intent(out) :: csrfl(:,:) ! cloudsat reflectivity 366 REAL(KIND=r8), intent(out) :: acsrfl(:,:) ! cloudsat average 367 REAL(KIND=r8), intent(out) :: fcsrfl(:,:) ! cloudsat fractional occurrence of radar reflectivity 368 REAL(KIND=r8), intent(out) :: rercld(:,:) ! effective radius calculation for rain + cloud 369 REAL(KIND=r8), intent(out) :: ncai(:,:) ! output number conc of ice nuclei available (1/m3) 370 REAL(KIND=r8), intent(out) :: ncal(:,:) ! output number conc of CCN (1/m3) 371 REAL(KIND=r8), intent(out) :: qrout2(:,:) ! copy of qrout as used to compute drout2 372 REAL(KIND=r8), intent(out) :: qsout2(:,:) ! copy of qsout as used to compute dsout2 373 REAL(KIND=r8), intent(out) :: nrout2(:,:) ! copy of nrout as used to compute drout2 374 REAL(KIND=r8), intent(out) :: nsout2(:,:) ! copy of nsout as used to compute dsout2 375 REAL(KIND=r8), intent(out) :: drout2(:,:) ! mean rain particle diameter (m) 376 REAL(KIND=r8), intent(out) :: dsout2(:,:) ! mean snow particle diameter (m) 377 REAL(KIND=r8), intent(out) :: freqs(:,:) ! fractional occurrence of snow 378 REAL(KIND=r8), intent(out) :: freqr(:,:) ! fractional occurrence of rain 379 REAL(KIND=r8), intent(out) :: nfice(:,:) ! fractional occurrence of ice 380 REAL(KIND=r8), intent(out) :: qcrat(:,:) ! limiter for qc process rates (1=no limit --> 0. no qc) 381 REAL(KIND=r8), intent(out) :: prer_evap(:,:) 382 CHARACTER(LEN=128), intent(out) :: errstring ! output status (non-blank for error return) 383 ! Tendencies calculated by external schemes that can replace MG's native 384 ! process tendencies. 385 ! Used with CARMA cirrus microphysics 386 ! (or similar external microphysics model) 387 REAL(KIND=r8), intent(in), pointer :: tnd_qsnow(:,:) ! snow mass tendency (kg/kg/s) 388 REAL(KIND=r8), intent(in), pointer :: tnd_nsnow(:,:) ! snow number tendency (#/kg/s) 389 REAL(KIND=r8), intent(in), pointer :: re_ice(:,:) ! ice effective radius (m) 390 ! From external ice nucleation. 391 REAL(KIND=r8), intent(in), pointer :: frzimm(:,:) ! Number tendency due to immersion freezing (1/cm3) 392 REAL(KIND=r8), intent(in), pointer :: frzcnt(:,:) ! Number tendency due to contact freezing (1/cm3) 393 REAL(KIND=r8), intent(in), pointer :: frzdep(:,:) ! Number tendency due to deposition nucleation (1/cm3) 394 ! local workspace 395 ! all units mks unless otherwise stated 396 ! local copies of input variables 397 REAL(KIND=r8) :: qc(mgncol,nlev) ! cloud liquid mixing ratio (kg/kg) 398 REAL(KIND=r8) :: qi(mgncol,nlev) ! cloud ice mixing ratio (kg/kg) 399 REAL(KIND=r8) :: nc(mgncol,nlev) ! cloud liquid number concentration (1/kg) 400 REAL(KIND=r8) :: ni(mgncol,nlev) ! cloud liquid number concentration (1/kg) 401 REAL(KIND=r8) :: qr(mgncol,nlev) ! rain mixing ratio (kg/kg) 402 REAL(KIND=r8) :: qs(mgncol,nlev) ! snow mixing ratio (kg/kg) 403 REAL(KIND=r8) :: nr(mgncol,nlev) ! rain number concentration (1/kg) 404 REAL(KIND=r8) :: ns(mgncol,nlev) ! snow number concentration (1/kg) 405 ! general purpose variables 406 REAL(KIND=r8) :: deltat ! sub-time step (s) 407 REAL(KIND=r8) :: mtime ! the assumed ice nucleation timescale 408 ! physical properties of the air at a given point 409 REAL(KIND=r8) :: rho(mgncol,nlev) ! density (kg m-3) 410 REAL(KIND=r8) :: dv(mgncol,nlev) ! diffusivity of water vapor 411 REAL(KIND=r8) :: mu(mgncol,nlev) ! viscosity 412 REAL(KIND=r8) :: sc(mgncol,nlev) ! schmidt number 413 REAL(KIND=r8) :: rhof(mgncol,nlev) ! density correction factor for fallspeed 414 ! cloud fractions 415 REAL(KIND=r8) :: precip_frac(mgncol,nlev) ! precip fraction assuming maximum overlap 416 REAL(KIND=r8) :: cldm(mgncol,nlev) ! cloud fraction 417 REAL(KIND=r8) :: icldm(mgncol,nlev) ! ice cloud fraction 418 REAL(KIND=r8) :: lcldm(mgncol,nlev) ! liq cloud fraction 419 ! mass mixing ratios 420 REAL(KIND=r8) :: qcic(mgncol,nlev) ! in-cloud cloud liquid 421 REAL(KIND=r8) :: qiic(mgncol,nlev) ! in-cloud cloud ice 422 REAL(KIND=r8) :: qsic(mgncol,nlev) ! in-precip snow 423 REAL(KIND=r8) :: qric(mgncol,nlev) ! in-precip rain 424 ! number concentrations 425 REAL(KIND=r8) :: ncic(mgncol,nlev) ! in-cloud droplet 426 REAL(KIND=r8) :: niic(mgncol,nlev) ! in-cloud cloud ice 427 REAL(KIND=r8) :: nsic(mgncol,nlev) ! in-precip snow 428 REAL(KIND=r8) :: nric(mgncol,nlev) ! in-precip rain 429 ! maximum allowed ni value 430 REAL(KIND=r8) :: nimax(mgncol,nlev) 431 ! Size distribution parameters for: 432 ! cloud ice 433 REAL(KIND=r8) :: lami(mgncol,nlev) ! slope 434 REAL(KIND=r8) :: n0i(mgncol,nlev) ! intercept 435 ! cloud liquid 436 REAL(KIND=r8) :: lamc(mgncol,nlev) ! slope 437 REAL(KIND=r8) :: pgam(mgncol,nlev) ! spectral width parameter 438 ! snow 439 REAL(KIND=r8) :: lams(mgncol,nlev) ! slope 440 REAL(KIND=r8) :: n0s(mgncol,nlev) ! intercept 441 ! rain 442 REAL(KIND=r8) :: lamr(mgncol,nlev) ! slope 443 REAL(KIND=r8) :: n0r(mgncol,nlev) ! intercept 444 ! Rates/tendencies due to: 445 ! Instantaneous snow melting 446 REAL(KIND=r8) :: minstsm(mgncol,nlev) ! mass mixing ratio 447 REAL(KIND=r8) :: ninstsm(mgncol,nlev) ! number concentration 448 ! Instantaneous rain freezing 449 REAL(KIND=r8) :: minstrf(mgncol,nlev) ! mass mixing ratio 450 REAL(KIND=r8) :: ninstrf(mgncol,nlev) ! number concentration 451 ! deposition of cloud ice 452 REAL(KIND=r8) :: vap_dep(mgncol,nlev) ! deposition from vapor to ice PMC 12/3/12 453 ! sublimation of cloud ice 454 REAL(KIND=r8) :: ice_sublim(mgncol,nlev) ! sublimation from ice to vapor PMC 12/3/12 455 ! ice nucleation 456 REAL(KIND=r8) :: nnuccd(mgncol,nlev) ! number rate from deposition/cond.-freezing 457 REAL(KIND=r8) :: mnuccd(mgncol,nlev) ! mass mixing ratio 458 ! freezing of cloud water 459 REAL(KIND=r8) :: mnuccc(mgncol,nlev) ! mass mixing ratio 460 REAL(KIND=r8) :: nnuccc(mgncol,nlev) ! number concentration 461 ! contact freezing of cloud water 462 REAL(KIND=r8) :: mnucct(mgncol,nlev) ! mass mixing ratio 463 REAL(KIND=r8) :: nnucct(mgncol,nlev) ! number concentration 464 ! deposition nucleation in mixed-phase clouds (from external scheme) 465 REAL(KIND=r8) :: mnudep(mgncol,nlev) ! mass mixing ratio 466 REAL(KIND=r8) :: nnudep(mgncol,nlev) ! number concentration 467 ! ice multiplication 468 REAL(KIND=r8) :: msacwi(mgncol,nlev) ! mass mixing ratio 469 REAL(KIND=r8) :: nsacwi(mgncol,nlev) ! number concentration 470 ! autoconversion of cloud droplets 471 REAL(KIND=r8) :: prc(mgncol,nlev) ! mass mixing ratio 472 REAL(KIND=r8) :: nprc(mgncol,nlev) ! number concentration (rain) 473 REAL(KIND=r8) :: nprc1(mgncol,nlev) ! number concentration (cloud droplets) 474 ! self-aggregation of snow 475 REAL(KIND=r8) :: nsagg(mgncol,nlev) ! number concentration 476 ! self-collection of rain 477 REAL(KIND=r8) :: nragg(mgncol,nlev) ! number concentration 478 ! collection of droplets by snow 479 REAL(KIND=r8) :: psacws(mgncol,nlev) ! mass mixing ratio 480 REAL(KIND=r8) :: npsacws(mgncol,nlev) ! number concentration 481 ! collection of rain by snow 482 REAL(KIND=r8) :: pracs(mgncol,nlev) ! mass mixing ratio 483 REAL(KIND=r8) :: npracs(mgncol,nlev) ! number concentration 484 ! freezing of rain 485 REAL(KIND=r8) :: mnuccr(mgncol,nlev) ! mass mixing ratio 486 REAL(KIND=r8) :: nnuccr(mgncol,nlev) ! number concentration 487 ! freezing of rain to form ice (mg add 4/26/13) 488 REAL(KIND=r8) :: mnuccri(mgncol,nlev) ! mass mixing ratio 489 REAL(KIND=r8) :: nnuccri(mgncol,nlev) ! number concentration 490 ! accretion of droplets by rain 491 REAL(KIND=r8) :: pra(mgncol,nlev) ! mass mixing ratio 492 REAL(KIND=r8) :: npra(mgncol,nlev) ! number concentration 493 ! autoconversion of cloud ice to snow 494 REAL(KIND=r8) :: prci(mgncol,nlev) ! mass mixing ratio 495 REAL(KIND=r8) :: nprci(mgncol,nlev) ! number concentration 496 ! accretion of cloud ice by snow 497 REAL(KIND=r8) :: prai(mgncol,nlev) ! mass mixing ratio 498 REAL(KIND=r8) :: nprai(mgncol,nlev) ! number concentration 499 ! evaporation of rain 500 REAL(KIND=r8) :: pre(mgncol,nlev) ! mass mixing ratio 501 ! sublimation of snow 502 REAL(KIND=r8) :: prds(mgncol,nlev) ! mass mixing ratio 503 ! number evaporation 504 REAL(KIND=r8) :: nsubi(mgncol,nlev) ! cloud ice 505 REAL(KIND=r8) :: nsubc(mgncol,nlev) ! droplet 506 REAL(KIND=r8) :: nsubs(mgncol,nlev) ! snow 507 REAL(KIND=r8) :: nsubr(mgncol,nlev) ! rain 508 ! bergeron process 509 REAL(KIND=r8) :: berg(mgncol,nlev) ! mass mixing ratio (cloud ice) 510 REAL(KIND=r8) :: bergs(mgncol,nlev) ! mass mixing ratio (snow) 511 ! fallspeeds 512 ! number-weighted 513 REAL(KIND=r8) :: uns(mgncol,nlev) ! snow 514 REAL(KIND=r8) :: unr(mgncol,nlev) ! rain 515 ! air density corrected fallspeed parameters 516 REAL(KIND=r8) :: arn(mgncol,nlev) ! rain 517 REAL(KIND=r8) :: asn(mgncol,nlev) ! snow 518 REAL(KIND=r8) :: acn(mgncol,nlev) ! cloud droplet 519 REAL(KIND=r8) :: ain(mgncol,nlev) ! cloud ice 520 ! Mass of liquid droplets used with external heterogeneous freezing. 521 REAL(KIND=r8) :: mi0l(mgncol) 522 ! saturation vapor pressures 523 REAL(KIND=r8) :: esl(mgncol,nlev) ! liquid 524 REAL(KIND=r8) :: esi(mgncol,nlev) ! ice 525 REAL(KIND=r8) :: esn ! checking for RH after rain evap 526 ! saturation vapor mixing ratios 527 REAL(KIND=r8) :: qvl(mgncol,nlev) ! liquid 528 REAL(KIND=r8) :: qvi(mgncol,nlev) ! ice 529 REAL(KIND=r8) :: qvn ! checking for RH after rain evap 530 ! relative humidity 531 REAL(KIND=r8) :: relhum(mgncol,nlev) 532 ! parameters for cloud water and cloud ice sedimentation calculations 533 REAL(KIND=r8) :: fc(nlev) 534 REAL(KIND=r8) :: fnc(nlev) 535 REAL(KIND=r8) :: fi(nlev) 536 REAL(KIND=r8) :: fni(nlev) 537 REAL(KIND=r8) :: fr(nlev) 538 REAL(KIND=r8) :: fnr(nlev) 539 REAL(KIND=r8) :: fs(nlev) 540 REAL(KIND=r8) :: fns(nlev) 541 REAL(KIND=r8) :: faloutc(nlev) 542 REAL(KIND=r8) :: faloutnc(nlev) 543 REAL(KIND=r8) :: falouti(nlev) 544 REAL(KIND=r8) :: faloutni(nlev) 545 REAL(KIND=r8) :: faloutr(nlev) 546 REAL(KIND=r8) :: faloutnr(nlev) 547 REAL(KIND=r8) :: falouts(nlev) 548 REAL(KIND=r8) :: faloutns(nlev) 549 REAL(KIND=r8) :: faltndc 550 REAL(KIND=r8) :: faltndnc 551 REAL(KIND=r8) :: faltndi 552 REAL(KIND=r8) :: faltndni 553 REAL(KIND=r8) :: faltndqie 554 REAL(KIND=r8) :: faltndqce 555 REAL(KIND=r8) :: faltndr 556 REAL(KIND=r8) :: faltndnr 557 REAL(KIND=r8) :: faltnds 558 REAL(KIND=r8) :: faltndns 559 REAL(KIND=r8) :: rainrt(mgncol,nlev) ! rain rate for reflectivity calculation 560 ! dummy variables 561 REAL(KIND=r8) :: dum 562 REAL(KIND=r8) :: dum1 563 REAL(KIND=r8) :: dum2 564 ! dummies for checking RH 565 REAL(KIND=r8) :: qtmp 566 REAL(KIND=r8) :: ttmp 567 ! dummies for conservation check 568 REAL(KIND=r8) :: ratio 569 REAL(KIND=r8) :: tmpfrz 570 ! dummies for in-cloud variables 571 REAL(KIND=r8) :: dumc(mgncol,nlev) ! qc 572 REAL(KIND=r8) :: dumnc(mgncol,nlev) ! nc 573 REAL(KIND=r8) :: dumi(mgncol,nlev) ! qi 574 REAL(KIND=r8) :: dumni(mgncol,nlev) ! ni 575 REAL(KIND=r8) :: dumr(mgncol,nlev) ! rain mixing ratio 576 REAL(KIND=r8) :: dumnr(mgncol,nlev) ! rain number concentration 577 REAL(KIND=r8) :: dums(mgncol,nlev) ! snow mixing ratio 578 REAL(KIND=r8) :: dumns(mgncol,nlev) ! snow number concentration 579 ! Array dummy variable 580 REAL(KIND=r8) :: dum_2d(mgncol,nlev) 581 ! loop array variables 582 ! "i" and "k" are column/level iterators for internal (MG) variables 583 ! "n" is used for other looping (currently just sedimentation) 584 INTEGER :: k 585 INTEGER :: i 586 INTEGER :: n 587 ! number of sub-steps for loops over "n" (for sedimentation) 588 INTEGER :: nstep 589 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 590 ! default return error message 591 errstring = ' ' 592 IF (.not. (do_cldice .or. (associated(tnd_qsnow) .and. associated(tnd_nsnow) .and. associated(re_ice)))) THEN 593 errstring = "MG's native cloud ice processes are disabled, but no replacement values were passed in." 594 END IF 595 IF (use_hetfrz_classnuc .and. (.not. (associated(frzimm) .and. associated(frzcnt) .and. associated(frzdep)))) THEN 596 errstring = "External heterogeneous freezing is enabled, but the required tendencies were not all passed in." 597 END IF 598 ! Process inputs 599 ! assign variable deltat to deltatin 600 deltat = deltatin 601 ! Copies of input concentrations that may be changed internally. 602 qc = qcn 603 nc = ncn 604 qi = qin 605 ni = nin 606 qr = qrn 607 nr = nrn 608 qs = qsn 609 ns = nsn 610 ! cldn: used to set cldm, unused for subcolumns 611 ! liqcldf: used to set lcldm, unused for subcolumns 612 ! icecldf: used to set icldm, unused for subcolumns 613 IF (microp_uniform) THEN 614 ! subcolumns, set cloud fraction variables to one 615 ! if cloud water or ice is present, if not present 616 ! set to mincld (mincld used instead of zero, to prevent 617 ! possible division by zero errors). 618 WHERE ( qc >= qsmall ) 619 lcldm = 1._r8 620 ELSEWHERE 621 lcldm = mincld 622 END WHERE 623 WHERE ( qi >= qsmall ) 624 icldm = 1._r8 625 ELSEWHERE 626 icldm = mincld 627 END WHERE 628 cldm = max(icldm, lcldm) 629 ELSE 630 ! get cloud fraction, check for minimum 631 cldm = max(cldn,mincld) 632 lcldm = max(liqcldf,mincld) 633 icldm = max(icecldf,mincld) 634 END IF 635 ! Initialize local variables 636 ! local physical properties 637 rho = p/(r*t) 638 dv = 8.794e-5_r8 * t**1.81_r8 / p 639 mu = 1.496e-6_r8 * t**1.5_r8 / (t + 120._r8) 640 sc = mu/(rho*dv) 641 ! air density adjustment for fallspeed parameters 642 ! includes air density correction factor to the 643 ! power of 0.54 following Heymsfield and Bansemer 2007 644 rhof = (rhosu/rho)**0.54_r8 645 arn = ar*rhof 646 asn = as*rhof 647 acn = g*rhow/(18._r8*mu) 648 ain = ai*(rhosu/rho)**0.35_r8 649 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 650 ! Get humidity and saturation vapor pressures 651 DO k=1,nlev 652 DO i=1,mgncol 653 CALL qsat_water(t(i,k), p(i,k), esl(i,k), qvl(i,k)) 654 ! make sure when above freezing that esi=esl, not active yet 655 IF (t(i,k) >= tmelt) THEN 656 esi(i,k) = esl(i,k) 657 qvi(i,k) = qvl(i,k) 658 ELSE 659 CALL qsat_ice(t(i,k), p(i,k), esi(i,k), qvi(i,k)) 660 END IF 661 END DO 662 END DO 663 relhum = q / max(qvl, qsmall) 664 !=============================================== 665 ! set mtime here to avoid answer-changing 666 mtime = deltat 667 ! initialize microphysics output 668 qcsevap = 0._r8 669 qisevap = 0._r8 670 qvres = 0._r8 671 cmeitot = 0._r8 672 vtrmc = 0._r8 673 vtrmi = 0._r8 674 qcsedten = 0._r8 675 qisedten = 0._r8 676 qrsedten = 0._r8 677 qssedten = 0._r8 678 pratot = 0._r8 679 prctot = 0._r8 680 mnuccctot = 0._r8 681 mnuccttot = 0._r8 682 msacwitot = 0._r8 683 psacwstot = 0._r8 684 bergstot = 0._r8 685 bergtot = 0._r8 686 melttot = 0._r8 687 homotot = 0._r8 688 qcrestot = 0._r8 689 prcitot = 0._r8 690 praitot = 0._r8 691 qirestot = 0._r8 692 mnuccrtot = 0._r8 693 pracstot = 0._r8 694 meltsdttot = 0._r8 695 frzrdttot = 0._r8 696 mnuccdtot = 0._r8 697 rflx = 0._r8 698 sflx = 0._r8 699 ! initialize precip output 700 qrout = 0._r8 701 qsout = 0._r8 702 nrout = 0._r8 703 nsout = 0._r8 704 ! for refl calc 705 rainrt = 0._r8 706 ! initialize rain size 707 rercld = 0._r8 708 qcsinksum_rate1ord = 0._r8 709 ! initialize variables for trop_mozart 710 nevapr = 0._r8 711 prer_evap = 0._r8 712 evapsnow = 0._r8 713 prain = 0._r8 714 prodsnow = 0._r8 715 cmeout = 0._r8 716 precip_frac = mincld 717 lamc = 0._r8 718 ! initialize microphysical tendencies 719 tlat = 0._r8 720 qvlat = 0._r8 721 qctend = 0._r8 722 qitend = 0._r8 723 qstend = 0._r8 724 qrtend = 0._r8 725 nctend = 0._r8 726 nitend = 0._r8 727 nrtend = 0._r8 728 nstend = 0._r8 729 ! initialize in-cloud and in-precip quantities to zero 730 qcic = 0._r8 731 qiic = 0._r8 732 qsic = 0._r8 733 qric = 0._r8 734 ncic = 0._r8 735 niic = 0._r8 736 nsic = 0._r8 737 nric = 0._r8 738 ! initialize precip at surface 739 prect = 0._r8 740 preci = 0._r8 741 ! initialize precip fallspeeds to zero 742 ums = 0._r8 743 uns = 0._r8 744 umr = 0._r8 745 unr = 0._r8 746 ! initialize limiter for output 747 qcrat = 1._r8 748 ! Many outputs have to be initialized here at the top to work around 749 ! ifort problems, even if they are always overwritten later. 750 effc = 10._r8 751 lamcrad = 0._r8 752 pgamrad = 0._r8 753 effc_fn = 10._r8 754 effi = 25._r8 755 deffi = 50._r8 756 qrout2 = 0._r8 757 nrout2 = 0._r8 758 drout2 = 0._r8 759 qsout2 = 0._r8 760 nsout2 = 0._r8 761 dsout = 0._r8 762 dsout2 = 0._r8 763 freqr = 0._r8 764 freqs = 0._r8 765 reff_rain = 0._r8 766 reff_snow = 0._r8 767 refl = -9999._r8 768 arefl = 0._r8 769 areflz = 0._r8 770 frefl = 0._r8 771 csrfl = 0._r8 772 acsrfl = 0._r8 773 fcsrfl = 0._r8 774 ncal = 0._r8 775 ncai = 0._r8 776 nfice = 0._r8 777 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 778 ! droplet activation 779 ! get provisional droplet number after activation. This is used for 780 ! all microphysical process calculations, for consistency with update of 781 ! droplet mass before microphysics 782 ! calculate potential for droplet activation if cloud water is present 783 ! tendency from activation (npccn) is read in from companion routine 784 ! output activated liquid and ice (convert from #/kg -> #/m3) 785 !-------------------------------------------------- 786 WHERE ( qc >= qsmall ) 787 nc = max(nc + npccn*deltat, 0._r8) 788 ncal = nc*rho/lcldm ! sghan minimum in #/cm3 789 ELSEWHERE 790 ncal = 0._r8 791 END WHERE 792 WHERE ( t < icenuct ) 793 ncai = naai*rho 794 ELSEWHERE 795 ncai = 0._r8 796 END WHERE 797 !=============================================== 798 ! ice nucleation if activated nuclei exist at t<-5C AND rhmini + 5% 799 !------------------------------------------------------- 800 IF (do_cldice) THEN 801 WHERE ( naai > 0._r8 .and. t < icenuct .and. relhum*esl/esi > rhmini+0.05_r8 ) 802 !if NAAI > 0. then set numice = naai (as before) 803 !note: this is gridbox averaged 804 nnuccd = (naai-ni/icldm)/mtime*icldm 805 nnuccd = max(nnuccd,0._r8) 806 nimax = naai*icldm 807 !Calc mass of new particles using new crystal mass... 808 !also this will be multiplied by mtime as nnuccd is... 809 mnuccd = nnuccd * mi0 810 ELSEWHERE 811 nnuccd = 0._r8 812 nimax = 0._r8 813 mnuccd = 0._r8 814 END WHERE 815 END IF 816 !============================================================================= 817 pre_vert_loop: DO k=1,nlev 818 pre_col_loop: DO i=1,mgncol 819 ! calculate instantaneous precip processes (melting and homogeneous freezing) 820 ! melting of snow at +2 C 821 IF (t(i,k) > snowmelt) THEN 822 IF (qs(i,k) > 0._r8) THEN 823 ! make sure melting snow doesn't reduce temperature below threshold 824 dum = -xlf/cpp*qs(i,k) 825 IF (t(i,k)+dum < snowmelt) THEN 826 dum = (t(i,k)-snowmelt)*cpp/xlf 827 dum = dum/qs(i,k) 828 dum = max(0._r8,dum) 829 dum = min(1._r8,dum) 830 ELSE 831 dum = 1._r8 832 END IF 833 minstsm(i,k) = dum*qs(i,k) 834 ninstsm(i,k) = dum*ns(i,k) 835 dum1 = -xlf*minstsm(i,k)/deltat 836 tlat(i,k) = tlat(i,k)+dum1 837 meltsdttot(i,k) = meltsdttot(i,k) + dum1 838 qs(i,k) = max(qs(i,k) - minstsm(i,k), 0._r8) 839 ns(i,k) = max(ns(i,k) - ninstsm(i,k), 0._r8) 840 qr(i,k) = max(qr(i,k) + minstsm(i,k), 0._r8) 841 nr(i,k) = max(nr(i,k) + ninstsm(i,k), 0._r8) 842 END IF 843 END IF 844 ! freezing of rain at -5 C 845 IF (t(i,k) < rainfrze) THEN 846 IF (qr(i,k) > 0._r8) THEN 847 ! make sure freezing rain doesn't increase temperature above threshold 848 dum = xlf/cpp*qr(i,k) 849 IF (t(i,k)+dum > rainfrze) THEN 850 dum = -(t(i,k)-rainfrze)*cpp/xlf 851 dum = dum/qr(i,k) 852 dum = max(0._r8,dum) 853 dum = min(1._r8,dum) 854 ELSE 855 dum = 1._r8 856 END IF 857 minstrf(i,k) = dum*qr(i,k) 858 ninstrf(i,k) = dum*nr(i,k) 859 ! heating tendency 860 dum1 = xlf*minstrf(i,k)/deltat 861 tlat(i,k) = tlat(i,k)+dum1 862 frzrdttot(i,k) = frzrdttot(i,k) + dum1 863 qr(i,k) = max(qr(i,k) - minstrf(i,k), 0._r8) 864 nr(i,k) = max(nr(i,k) - ninstrf(i,k), 0._r8) 865 qs(i,k) = max(qs(i,k) + minstrf(i,k), 0._r8) 866 ns(i,k) = max(ns(i,k) + ninstrf(i,k), 0._r8) 867 END IF 868 END IF 869 ! obtain in-cloud values of cloud water/ice mixing ratios and number concentrations 870 !------------------------------------------------------- 871 ! for microphysical process calculations 872 ! units are kg/kg for mixing ratio, 1/kg for number conc 873 IF (qc(i,k).ge.qsmall) THEN 874 ! limit in-cloud values to 0.005 kg/kg 875 qcic(i,k) = min(qc(i,k)/lcldm(i,k),5.e-3_r8) 876 ncic(i,k) = max(nc(i,k)/lcldm(i,k),0._r8) 877 ! specify droplet concentration 878 IF (nccons) THEN 879 ncic(i,k) = ncnst/rho(i,k) 880 END IF 881 ELSE 882 qcic(i,k) = 0._r8 883 ncic(i,k) = 0._r8 884 END IF 885 IF (qi(i,k).ge.qsmall) THEN 886 ! limit in-cloud values to 0.005 kg/kg 887 qiic(i,k) = min(qi(i,k)/icldm(i,k),5.e-3_r8) 888 niic(i,k) = max(ni(i,k)/icldm(i,k),0._r8) 889 ! switch for specification of cloud ice number 890 IF (nicons) THEN 891 niic(i,k) = ninst/rho(i,k) 892 END IF 893 ELSE 894 qiic(i,k) = 0._r8 895 niic(i,k) = 0._r8 896 END IF 897 END DO pre_col_loop 898 END DO pre_vert_loop 899 !======================================================================== 900 ! for sub-columns cldm has already been set to 1 if cloud 901 ! water or ice is present, so precip_frac will be correctly set below 902 ! and nothing extra needs to be done here 903 precip_frac = cldm 904 micro_vert_loop: DO k=1,nlev 905 IF (trim(micro_mg_precip_frac_method) == 'in_cloud') THEN 906 IF (k /= 1) THEN 907 WHERE ( qc(:,k) < qsmall .and. qi(:,k) < qsmall ) 908 precip_frac(:,k) = precip_frac(:,k-1) 909 END WHERE 910 END IF 911 ELSE IF (trim(micro_mg_precip_frac_method) == 'max_overlap') THEN 912 ! calculate precip fraction based on maximum overlap assumption 913 ! if rain or snow mix ratios are smaller than threshold, 914 ! then leave precip_frac as cloud fraction at current level 915 IF (k /= 1) THEN 916 WHERE ( qr(:,k-1) >= qsmall .or. qs(:,k-1) >= qsmall ) 917 precip_frac(:,k) = max(precip_frac(:,k-1),precip_frac(:,k)) 918 END WHERE 919 END IF 920 END IF 921 DO i = 1, mgncol 922 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 923 ! get size distribution parameters based on in-cloud cloud water 924 ! these calculations also ensure consistency between number and mixing ratio 925 !cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 926 ! cloud liquid 927 !------------------------------------------- 928 CALL size_dist_param_liq(mg_liq_props, qcic(i,k), ncic(i,k), rho(i,k), pgam(i,k), lamc(i,k)) 929 END DO 930 !======================================================================== 931 ! autoconversion of cloud liquid water to rain 932 ! formula from Khrouditnov and Kogan (2000), modified for sub-grid distribution of qc 933 ! minimum qc of 1 x 10^-8 prevents floating point error 934 CALL kk2000_liq_autoconversion(microp_uniform, qcic(:,k), ncic(:,k), rho(:,k), relvar(:,k), prc(:,k), nprc(:,k), & 935 nprc1(:,k)) 936 ! assign qric based on prognostic qr, using assumed precip fraction 937 ! note: this could be moved above for consistency with qcic and qiic calculations 938 qric(:,k) = qr(:,k)/precip_frac(:,k) 939 nric(:,k) = nr(:,k)/precip_frac(:,k) 940 ! limit in-precip mixing ratios to 10 g/kg 941 qric(:,k) = min(qric(:,k),0.01_r8) 942 ! add autoconversion to precip from above to get provisional rain mixing ratio 943 ! and number concentration (qric and nric) 944 WHERE ( qric(:,k).lt.qsmall ) 945 qric(:,k) = 0._r8 946 nric(:,k) = 0._r8 947 END WHERE 948 ! make sure number concentration is a positive number to avoid 949 ! taking root of negative later 950 nric(:,k) = max(nric(:,k),0._r8) 951 ! Get size distribution parameters for cloud ice 952 CALL size_dist_param_basic(mg_ice_props, qiic(:,k), niic(:,k), lami(:,k), n0i(:,k)) 953 !....................................................................... 954 ! Autoconversion of cloud ice to snow 955 ! similar to Ferrier (1994) 956 IF (do_cldice) THEN 957 CALL ice_autoconversion(t(:,k), qiic(:,k), lami(:,k), n0i(:,k), dcs, prci(:,k), nprci(:,k)) 958 ELSE 959 ! Add in the particles that we have already converted to snow, and 960 ! don't do any further autoconversion of ice. 961 prci(:,k) = tnd_qsnow(:,k) / cldm(:,k) 962 nprci(:,k) = tnd_nsnow(:,k) / cldm(:,k) 963 END IF 964 ! note, currently we don't have this 965 ! inside the do_cldice block, should be changed later 966 ! assign qsic based on prognostic qs, using assumed precip fraction 967 qsic(:,k) = qs(:,k)/precip_frac(:,k) 968 nsic(:,k) = ns(:,k)/precip_frac(:,k) 969 ! limit in-precip mixing ratios to 10 g/kg 970 qsic(:,k) = min(qsic(:,k),0.01_r8) 971 ! if precip mix ratio is zero so should number concentration 972 WHERE ( qsic(:,k) < qsmall ) 973 qsic(:,k) = 0._r8 974 nsic(:,k) = 0._r8 975 END WHERE 976 ! make sure number concentration is a positive number to avoid 977 ! taking root of negative later 978 nsic(:,k) = max(nsic(:,k),0._r8) 979 !....................................................................... 980 ! get size distribution parameters for precip 981 !...................................................................... 982 ! rain 983 CALL size_dist_param_basic(mg_rain_props, qric(:,k), nric(:,k), lamr(:,k), n0r(:,k)) 984 WHERE ( lamr(:,k) >= qsmall ) 985 ! provisional rain number and mass weighted mean fallspeed (m/s) 986 unr(:,k) = min(arn(:,k)*gamma_br_plus1/lamr(:,k)**br,9.1_r8*rhof(:,k)) 987 umr(:,k) = min(arn(:,k)*gamma_br_plus4/(6._r8*lamr(:,k)**br),9.1_r8*rhof(:,k)) 988 ELSEWHERE 989 umr(:,k) = 0._r8 990 unr(:,k) = 0._r8 991 END WHERE 992 !...................................................................... 993 ! snow 994 CALL size_dist_param_basic(mg_snow_props, qsic(:,k), nsic(:,k), lams(:,k), n0s(:,k)) 995 WHERE ( lams(:,k) > 0._r8 ) 996 ! provisional snow number and mass weighted mean fallspeed (m/s) 997 ums(:,k) = min(asn(:,k)*gamma_bs_plus4/(6._r8*lams(:,k)**bs),1.2_r8*rhof(:,k)) 998 uns(:,k) = min(asn(:,k)*gamma_bs_plus1/lams(:,k)**bs,1.2_r8*rhof(:,k)) 999 ELSEWHERE 1000 ums(:,k) = 0._r8 1001 uns(:,k) = 0._r8 1002 END WHERE 1003 IF (do_cldice) THEN 1004 IF (.not. use_hetfrz_classnuc) THEN 1005 ! heterogeneous freezing of cloud water 1006 !---------------------------------------------- 1007 CALL immersion_freezing(microp_uniform, t(:,k), pgam(:,k), lamc(:,k), qcic(:,k), ncic(:,k), relvar(:,k), & 1008 mnuccc(:,k), nnuccc(:,k)) 1009 ! make sure number of droplets frozen does not exceed available ice nuclei concentration 1010 ! this prevents 'runaway' droplet freezing 1011 WHERE ( qcic(:,k).ge.qsmall .and. t(:,k).lt.269.15_r8 ) 1012 WHERE ( nnuccc(:,k)*lcldm(:,k).gt.nnuccd(:,k) ) 1013 ! scale mixing ratio of droplet freezing with limit 1014 mnuccc(:,k) = mnuccc(:,k)*(nnuccd(:,k)/(nnuccc(:,k)*lcldm(:,k))) 1015 nnuccc(:,k) = nnuccd(:,k)/lcldm(:,k) 1016 END WHERE 1017 END WHERE 1018 CALL contact_freezing(microp_uniform, t(:,k), p(:,k), rndst(:,k,:), nacon(:,k,:), pgam(:,k), lamc(:,k), & 1019 qcic(:,k), ncic(:,k), relvar(:,k), mnucct(:,k), nnucct(:,k)) 1020 mnudep(:,k) = 0._r8 1021 nnudep(:,k) = 0._r8 1022 ELSE 1023 ! Mass of droplets frozen is the average droplet mass, except 1024 ! with two limiters: concentration must be at least 1/cm^3, and 1025 ! mass must be at least the minimum defined above. 1026 mi0l = qcic(:,k)/max(ncic(:,k), 1.0e6_r8/rho(:,k)) 1027 mi0l = max(mi0l_min, mi0l) 1028 WHERE ( qcic(:,k) >= qsmall ) 1029 nnuccc(:,k) = frzimm(:,k)*1.0e6_r8/rho(:,k) 1030 mnuccc(:,k) = nnuccc(:,k)*mi0l 1031 nnucct(:,k) = frzcnt(:,k)*1.0e6_r8/rho(:,k) 1032 mnucct(:,k) = nnucct(:,k)*mi0l 1033 nnudep(:,k) = frzdep(:,k)*1.0e6_r8/rho(:,k) 1034 mnudep(:,k) = nnudep(:,k)*mi0 1035 ELSEWHERE 1036 nnuccc(:,k) = 0._r8 1037 mnuccc(:,k) = 0._r8 1038 nnucct(:,k) = 0._r8 1039 mnucct(:,k) = 0._r8 1040 nnudep(:,k) = 0._r8 1041 mnudep(:,k) = 0._r8 1042 END WHERE 1043 END IF 1044 ELSE 1045 mnuccc(:,k) = 0._r8 1046 nnuccc(:,k) = 0._r8 1047 mnucct(:,k) = 0._r8 1048 nnucct(:,k) = 0._r8 1049 mnudep(:,k) = 0._r8 1050 nnudep(:,k) = 0._r8 1051 END IF 1052 CALL snow_self_aggregation(t(:,k), rho(:,k), asn(:,k), rhosn, qsic(:,k), nsic(:,k), nsagg(:,k)) 1053 CALL accrete_cloud_water_snow(t(:,k), rho(:,k), asn(:,k), uns(:,k), mu(:,k), qcic(:,k), ncic(:,k), qsic(:,k), & 1054 pgam(:,k), lamc(:,k), lams(:,k), n0s(:,k), psacws(:,k), npsacws(:,k)) 1055 IF (do_cldice) THEN 1056 CALL secondary_ice_production(t(:,k), psacws(:,k), msacwi(:,k), nsacwi(:,k)) 1057 ELSE 1058 nsacwi(:,k) = 0.0_r8 1059 msacwi(:,k) = 0.0_r8 1060 END IF 1061 CALL accrete_rain_snow(t(:,k), rho(:,k), umr(:,k), ums(:,k), unr(:,k), uns(:,k), qric(:,k), qsic(:,k), lamr(:,k), & 1062 n0r(:,k), lams(:,k), n0s(:,k), pracs(:,k), npracs(:,k)) 1063 CALL heterogeneous_rain_freezing(t(:,k), qric(:,k), nric(:,k), lamr(:,k), mnuccr(:,k), nnuccr(:,k)) 1064 CALL accrete_cloud_water_rain(microp_uniform, qric(:,k), qcic(:,k), ncic(:,k), relvar(:,k), accre_enhan(:,k), pra(& 1065 :,k), npra(:,k)) 1066 CALL self_collection_rain(rho(:,k), qric(:,k), nric(:,k), nragg(:,k)) 1067 IF (do_cldice) THEN 1068 CALL accrete_cloud_ice_snow(t(:,k), rho(:,k), asn(:,k), qiic(:,k), niic(:,k), qsic(:,k), lams(:,k), n0s(:,k), & 1069 prai(:,k), nprai(:,k)) 1070 ELSE 1071 prai(:,k) = 0._r8 1072 nprai(:,k) = 0._r8 1073 END IF 1074 CALL evaporate_sublimate_precip(t(:,k), rho(:,k), dv(:,k), mu(:,k), sc(:,k), q(:,k), qvl(:,k), qvi(:,k), lcldm(:,& 1075 k), precip_frac(:,k), arn(:,k), asn(:,k), qcic(:,k), qiic(:,k), qric(:,k), qsic(:,k), lamr(:,k), n0r(:,k), lams(:,& 1076 k), n0s(:,k), pre(:,k), prds(:,k)) 1077 CALL bergeron_process_snow(t(:,k), rho(:,k), dv(:,k), mu(:,k), sc(:,k), qvl(:,k), qvi(:,k), asn(:,k), qcic(:,k), & 1078 qsic(:,k), lams(:,k), n0s(:,k), bergs(:,k)) 1079 bergs(:,k) = bergs(:,k)*micro_mg_berg_eff_factor 1080 !+++PMC 12/3/12 - NEW VAPOR DEP/SUBLIMATION GOES HERE!!! 1081 IF (do_cldice) THEN 1082 CALL ice_deposition_sublimation(t(:,k), q(:,k), qi(:,k), ni(:,k), icldm(:,k), rho(:,k), dv(:,k), qvl(:,k), & 1083 qvi(:,k), berg(:,k), vap_dep(:,k), ice_sublim(:,k)) 1084 berg(:,k) = berg(:,k)*micro_mg_berg_eff_factor 1085 WHERE ( vap_dep(:,k) < 0._r8 .and. qi(:,k) > qsmall .and. icldm(:,k) > mincld ) 1086 nsubi(:,k) = vap_dep(:,k) / qi(:,k) * ni(:,k) / icldm(:,k) 1087 ELSEWHERE 1088 nsubi(:,k) = 0._r8 1089 END WHERE 1090 ! bergeron process should not reduce nc unless 1091 ! all ql is removed (which is handled elsewhere) 1092 !in fact, nothing in this entire file makes nsubc nonzero. 1093 nsubc(:,k) = 0._r8 1094 END IF !do_cldice 1095 !---PMC 12/3/12 1096 DO i=1,mgncol 1097 ! conservation to ensure no negative values of cloud water/precipitation 1098 ! in case microphysical process rates are large 1099 !=================================================================== 1100 ! note: for check on conservation, processes are multiplied by omsm 1101 ! to prevent problems due to round off error 1102 ! conservation of qc 1103 !------------------------------------------------------------------- 1104 dum = ((prc(i,k)+pra(i,k)+mnuccc(i,k)+mnucct(i,k)+msacwi(i,k)+ psacws(i,k)+bergs(i,k))*lcldm(i,k)& 1105 +berg(i,k))*deltat 1106 IF (dum.gt.qc(i,k)) THEN 1107 ratio = qc(i,k)/deltat/((prc(i,k)+pra(i,k)+mnuccc(i,k)+mnucct(i,k)+ msacwi(i,k)+psacws(i,& 1108 k)+bergs(i,k))*lcldm(i,k)+berg(i,k))*omsm 1109 prc(i,k) = prc(i,k)*ratio 1110 pra(i,k) = pra(i,k)*ratio 1111 mnuccc(i,k) = mnuccc(i,k)*ratio 1112 mnucct(i,k) = mnucct(i,k)*ratio 1113 msacwi(i,k) = msacwi(i,k)*ratio 1114 psacws(i,k) = psacws(i,k)*ratio 1115 bergs(i,k) = bergs(i,k)*ratio 1116 berg(i,k) = berg(i,k)*ratio 1117 qcrat(i,k) = ratio 1118 ELSE 1119 qcrat(i,k) = 1._r8 1120 END IF 1121 !PMC 12/3/12: ratio is also frac of step w/ liquid. 1122 !thus we apply berg for "ratio" of timestep and vapor 1123 !deposition for the remaining frac of the timestep. 1124 IF (qc(i,k) >= qsmall) THEN 1125 vap_dep(i,k) = vap_dep(i,k)*(1._r8-qcrat(i,k)) 1126 END IF 1127 END DO 1128 DO i=1,mgncol 1129 !================================================================= 1130 ! apply limiter to ensure that ice/snow sublimation and rain evap 1131 ! don't push conditions into supersaturation, and ice deposition/nucleation don't 1132 ! push conditions into sub-saturation 1133 ! note this is done after qc conservation since we don't know how large 1134 ! vap_dep is before then 1135 ! estimates are only approximate since other process terms haven't been limited 1136 ! for conservation yet 1137 ! first limit ice deposition/nucleation vap_dep + mnuccd 1138 dum1 = vap_dep(i,k) + mnuccd(i,k) 1139 IF (dum1 > 1.e-20_r8) THEN 1140 dum = (q(i,k)-qvi(i,k))/(1._r8 + xxls_squared*qvi(i,k)/(cpp*rv*t(i,k)**2))/deltat 1141 dum = max(dum,0._r8) 1142 IF (dum1 > dum) THEN 1143 ! Allocate the limited "dum" tendency to mnuccd and vap_dep 1144 ! processes. Don't divide by cloud fraction; these are grid- 1145 ! mean rates. 1146 dum1 = mnuccd(i,k) / (vap_dep(i,k)+mnuccd(i,k)) 1147 mnuccd(i,k) = dum*dum1 1148 vap_dep(i,k) = dum - mnuccd(i,k) 1149 END IF 1150 END IF 1151 END DO 1152 DO i=1,mgncol 1153 !=================================================================== 1154 ! conservation of nc 1155 !------------------------------------------------------------------- 1156 dum = (nprc1(i,k)+npra(i,k)+nnuccc(i,k)+nnucct(i,k)+ npsacws(i,k)-nsubc(i,k))*lcldm(i,k)*deltat 1157 IF (dum.gt.nc(i,k)) THEN 1158 ratio = nc(i,k)/deltat/((nprc1(i,k)+npra(i,k)+nnuccc(i,k)+nnucct(i,k)+ npsacws(i,k)-nsubc(& 1159 i,k))*lcldm(i,k))*omsm 1160 nprc1(i,k) = nprc1(i,k)*ratio 1161 npra(i,k) = npra(i,k)*ratio 1162 nnuccc(i,k) = nnuccc(i,k)*ratio 1163 nnucct(i,k) = nnucct(i,k)*ratio 1164 npsacws(i,k) = npsacws(i,k)*ratio 1165 nsubc(i,k) = nsubc(i,k)*ratio 1166 END IF 1167 mnuccri(i,k) = 0._r8 1168 nnuccri(i,k) = 0._r8 1169 IF (do_cldice) THEN 1170 ! freezing of rain to produce ice if mean rain size is smaller than Dcs 1171 IF (lamr(i,k) > qsmall .and. 1._r8/lamr(i,k) < dcs) THEN 1172 mnuccri(i,k) = mnuccr(i,k) 1173 nnuccri(i,k) = nnuccr(i,k) 1174 mnuccr(i,k) = 0._r8 1175 nnuccr(i,k) = 0._r8 1176 END IF 1177 END IF 1178 END DO 1179 DO i=1,mgncol 1180 ! conservation of rain mixing ratio 1181 !------------------------------------------------------------------- 1182 dum = ((-pre(i,k)+pracs(i,k)+mnuccr(i,k)+mnuccri(i,k))*precip_frac(i,k)- (pra(i,k)+prc(i,k))& 1183 *lcldm(i,k))*deltat 1184 ! note that qrtend is included below because of instantaneous freezing/melt 1185 IF (dum.gt.qr(i,k).and. (-pre(i,k)+pracs(i,k)+mnuccr(i,k)+mnuccri(i,k)).ge.qsmall) THEN 1186 ratio = (qr(i,k)/deltat+(pra(i,k)+prc(i,k))*lcldm(i,k))/ precip_frac(i,k)/(-pre(i,k)& 1187 +pracs(i,k)+mnuccr(i,k)+mnuccri(i,k))*omsm 1188 pre(i,k) = pre(i,k)*ratio 1189 pracs(i,k) = pracs(i,k)*ratio 1190 mnuccr(i,k) = mnuccr(i,k)*ratio 1191 mnuccri(i,k) = mnuccri(i,k)*ratio 1192 END IF 1193 END DO 1194 DO i=1,mgncol 1195 ! conservation of rain number 1196 !------------------------------------------------------------------- 1197 ! Add evaporation of rain number. 1198 IF (pre(i,k) < 0._r8) THEN 1199 dum = pre(i,k)*deltat/qr(i,k) 1200 dum = max(-1._r8,dum) 1201 nsubr(i,k) = dum*nr(i,k)/deltat 1202 ELSE 1203 nsubr(i,k) = 0._r8 1204 END IF 1205 END DO 1206 DO i=1,mgncol 1207 dum = ((-nsubr(i,k)+npracs(i,k)+nnuccr(i,k)+nnuccri(i,k)-nragg(i,k))*precip_frac(i,k)- nprc(i,k)& 1208 *lcldm(i,k))*deltat 1209 IF (dum.gt.nr(i,k)) THEN 1210 ratio = (nr(i,k)/deltat+nprc(i,k)*lcldm(i,k)/precip_frac(i,k))/ (-nsubr(i,k)+npracs(i,k)& 1211 +nnuccr(i,k)+nnuccri(i,k)-nragg(i,k))*omsm 1212 nragg(i,k) = nragg(i,k)*ratio 1213 npracs(i,k) = npracs(i,k)*ratio 1214 nnuccr(i,k) = nnuccr(i,k)*ratio 1215 nsubr(i,k) = nsubr(i,k)*ratio 1216 nnuccri(i,k) = nnuccri(i,k)*ratio 1217 END IF 1218 END DO 1219 IF (do_cldice) THEN 1220 DO i=1,mgncol 1221 ! conservation of qi 1222 !------------------------------------------------------------------- 1223 dum = ((-mnuccc(i,k)-mnucct(i,k)-mnudep(i,k)-msacwi(i,k))*lcldm(i,k)+(prci(i,k)+ prai(i,k)& 1224 )*icldm(i,k)-mnuccri(i,k)*precip_frac(i,k) -ice_sublim(i,k)-vap_dep(i,k)-berg(i,k)-mnuccd(& 1225 i,k))*deltat 1226 IF (dum.gt.qi(i,k)) THEN 1227 ratio = (qi(i,k)/deltat+vap_dep(i,k)+berg(i,k)+mnuccd(i,k)+ (mnuccc(i,k)+mnucct(i,& 1228 k)+mnudep(i,k)+msacwi(i,k))*lcldm(i,k)+ mnuccri(i,k)*precip_frac(i,k))/ & 1229 ((prci(i,k)+prai(i,k))*icldm(i,k)-ice_sublim(i,k))*omsm 1230 prci(i,k) = prci(i,k)*ratio 1231 prai(i,k) = prai(i,k)*ratio 1232 ice_sublim(i,k) = ice_sublim(i,k)*ratio 1233 END IF 1234 END DO 1235 END IF 1236 IF (do_cldice) THEN 1237 DO i=1,mgncol 1238 ! conservation of ni 1239 !------------------------------------------------------------------- 1240 IF (use_hetfrz_classnuc) THEN 1241 tmpfrz = nnuccc(i,k) 1242 ELSE 1243 tmpfrz = 0._r8 1244 END IF 1245 dum = ((-nnucct(i,k)-tmpfrz-nnudep(i,k)-nsacwi(i,k))*lcldm(i,k)+(nprci(i,k)+ nprai(i,k)& 1246 -nsubi(i,k))*icldm(i,k)-nnuccri(i,k)*precip_frac(i,k)- nnuccd(i,k))*deltat 1247 IF (dum.gt.ni(i,k)) THEN 1248 ratio = (ni(i,k)/deltat+nnuccd(i,k)+ (nnucct(i,k)+tmpfrz+nnudep(i,k)+nsacwi(i,k))& 1249 *lcldm(i,k)+ nnuccri(i,k)*precip_frac(i,k))/ ((nprci(i,k)+nprai(& 1250 i,k)-nsubi(i,k))*icldm(i,k))*omsm 1251 nprci(i,k) = nprci(i,k)*ratio 1252 nprai(i,k) = nprai(i,k)*ratio 1253 nsubi(i,k) = nsubi(i,k)*ratio 1254 END IF 1255 END DO 1256 END IF 1257 DO i=1,mgncol 1258 ! conservation of snow mixing ratio 1259 !------------------------------------------------------------------- 1260 dum = (-(prds(i,k)+pracs(i,k)+mnuccr(i,k))*precip_frac(i,k)-(prai(i,k)+prci(i,k))*icldm(i,k) -(& 1261 bergs(i,k)+psacws(i,k))*lcldm(i,k))*deltat 1262 IF (dum.gt.qs(i,k).and.-prds(i,k).ge.qsmall) THEN 1263 ratio = (qs(i,k)/deltat+(prai(i,k)+prci(i,k))*icldm(i,k)+ (bergs(i,k)+psacws(i,k))*lcldm(& 1264 i,k)+(pracs(i,k)+mnuccr(i,k))*precip_frac(i,k))/ precip_frac(i,k)/(-prds(i,k))*omsm 1265 prds(i,k) = prds(i,k)*ratio 1266 END IF 1267 END DO 1268 DO i=1,mgncol 1269 ! conservation of snow number 1270 !------------------------------------------------------------------- 1271 ! calculate loss of number due to sublimation 1272 ! for now neglect sublimation of ns 1273 nsubs(i,k) = 0._r8 1274 dum = ((-nsagg(i,k)-nsubs(i,k)-nnuccr(i,k))*precip_frac(i,k)-nprci(i,k)*icldm(i,k))*deltat 1275 IF (dum.gt.ns(i,k)) THEN 1276 ratio = (ns(i,k)/deltat+nnuccr(i,k)* precip_frac(i,k)+nprci(i,k)*icldm(i,k))/precip_frac(& 1277 i,k)/ (-nsubs(i,k)-nsagg(i,k))*omsm 1278 nsubs(i,k) = nsubs(i,k)*ratio 1279 nsagg(i,k) = nsagg(i,k)*ratio 1280 END IF 1281 END DO 1282 DO i=1,mgncol 1283 ! next limit ice and snow sublimation and rain evaporation 1284 ! get estimate of q and t at end of time step 1285 ! don't include other microphysical processes since they haven't 1286 ! been limited via conservation checks yet 1287 IF ((pre(i,k)+prds(i,k))*precip_frac(i,k)+ice_sublim(i,k) < -1.e-20_r8) THEN 1288 qtmp = q(i,k)-(ice_sublim(i,k)+vap_dep(i,k)+mnuccd(i,k)+ (pre(i,k)+prds(i,k))*precip_frac(& 1289 i,k))*deltat 1290 ttmp = t(i,k)+((pre(i,k)*precip_frac(i,k))*xxlv+ (prds(i,k)*precip_frac(i,k)+vap_dep(i,k)& 1291 +ice_sublim(i,k)+mnuccd(i,k))*xxls)*deltat/cpp 1292 ! use rhw to allow ice supersaturation 1293 CALL qsat_water(ttmp, p(i,k), esn, qvn) 1294 ! modify ice/precip evaporation rate if q > qsat 1295 IF (qtmp > qvn) THEN 1296 dum1 = pre(i,k)*precip_frac(i,k)/((pre(i,k)+prds(i,k))*precip_frac(i,k)+ice_sublim(i,k)) 1297 dum2 = prds(i,k)*precip_frac(i,k)/((pre(i,k)+prds(i,k))*precip_frac(i,k)+ice_sublim(i,k)) 1298 ! recalculate q and t after vap_dep and mnuccd but without evap or sublim 1299 qtmp = q(i,k)-(vap_dep(i,k)+mnuccd(i,k))*deltat 1300 ttmp = t(i,k)+((vap_dep(i,k)+mnuccd(i,k))*xxls)*deltat/cpp 1301 ! use rhw to allow ice supersaturation 1302 CALL qsat_water(ttmp, p(i,k), esn, qvn) 1303 dum = (qtmp-qvn)/(1._r8 + xxlv_squared*qvn/(cpp*rv*ttmp**2)) 1304 dum = min(dum,0._r8) 1305 ! modify rates if needed, divide by precip_frac to get local (in-precip) value 1306 pre(i,k) = dum*dum1/deltat/precip_frac(i,k) 1307 ! do separately using RHI for prds and ice_sublim 1308 CALL qsat_ice(ttmp, p(i,k), esn, qvn) 1309 dum = (qtmp-qvn)/(1._r8 + xxls_squared*qvn/(cpp*rv*ttmp**2)) 1310 dum = min(dum,0._r8) 1311 ! modify rates if needed, divide by precip_frac to get local (in-precip) value 1312 prds(i,k) = dum*dum2/deltat/precip_frac(i,k) 1313 ! don't divide ice_sublim by cloud fraction since it is grid-averaged 1314 dum1 = (1._r8-dum1-dum2) 1315 ice_sublim(i,k) = dum*dum1/deltat 1316 END IF 1317 END IF 1318 END DO 1319 ! Big "administration" loop enforces conservation, updates variables 1320 ! that accumulate over substeps, and sets output variables. 1321 DO i=1,mgncol 1322 ! get tendencies due to microphysical conversion processes 1323 !========================================================== 1324 ! note: tendencies are multiplied by appropriate cloud/precip 1325 ! fraction to get grid-scale values 1326 ! note: vap_dep is already grid-average values 1327 ! The net tendencies need to be added to rather than overwritten, 1328 ! because they may have a value already set for instantaneous 1329 ! melting/freezing. 1330 qvlat(i,k) = qvlat(i,k)-(pre(i,k)+prds(i,k))*precip_frac(i,k)- vap_dep(i,k)-ice_sublim(i,k)& 1331 -mnuccd(i,k)-mnudep(i,k)*lcldm(i,k) 1332 tlat(i,k) = tlat(i,k)+((pre(i,k)*precip_frac(i,k)) *xxlv+(prds(i,k)*precip_frac(i,k)+vap_dep(i,k)& 1333 +ice_sublim(i,k)+mnuccd(i,k)+mnudep(i,k)*lcldm(i,k))*xxls+ ((bergs(i,k)+psacws(i,k)+mnuccc(i,k)& 1334 +mnucct(i,k)+msacwi(i,k))*lcldm(i,k)+(mnuccr(i,k)+ pracs(i,k)+mnuccri(i,k))*precip_frac(i,k)& 1335 +berg(i,k))*xlf) 1336 qctend(i,k) = qctend(i,k)+ (-pra(i,k)-prc(i,k)-mnuccc(i,k)-mnucct(i,k)-msacwi(i,k)- & 1337 psacws(i,k)-bergs(i,k))*lcldm(i,k)-berg(i,k) 1338 IF (do_cldice) THEN 1339 qitend(i,k) = qitend(i,k)+ (mnuccc(i,k)+mnucct(i,k)+mnudep(i,k)+msacwi(i,k))*lcldm(i,k)+(& 1340 -prci(i,k)- prai(i,k))*icldm(i,k)+vap_dep(i,k)+berg(i,k)+ice_sublim(i,k)+ & 1341 mnuccd(i,k)+mnuccri(i,k)*precip_frac(i,k) 1342 END IF 1343 qrtend(i,k) = qrtend(i,k)+ (pra(i,k)+prc(i,k))*lcldm(i,k)+(pre(i,k)-pracs(i,k)- & 1344 mnuccr(i,k)-mnuccri(i,k))*precip_frac(i,k) 1345 qstend(i,k) = qstend(i,k)+ (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k)+(& 1346 prds(i,k)+ pracs(i,k)+mnuccr(i,k))*precip_frac(i,k) 1347 cmeout(i,k) = vap_dep(i,k) + ice_sublim(i,k) + mnuccd(i,k) 1348 ! add output for cmei (accumulate) 1349 cmeitot(i,k) = vap_dep(i,k) + ice_sublim(i,k) + mnuccd(i,k) 1350 ! assign variables for trop_mozart, these are grid-average 1351 !------------------------------------------------------------------- 1352 ! evaporation/sublimation is stored here as positive term 1353 evapsnow(i,k) = -prds(i,k)*precip_frac(i,k) 1354 nevapr(i,k) = -pre(i,k)*precip_frac(i,k) 1355 prer_evap(i,k) = -pre(i,k)*precip_frac(i,k) 1356 ! change to make sure prain is positive: do not remove snow from 1357 ! prain used for wet deposition 1358 prain(i,k) = (pra(i,k)+prc(i,k))*lcldm(i,k)+(-pracs(i,k)- mnuccr(i,k)-mnuccri(i,k))*precip_frac(& 1359 i,k) 1360 prodsnow(i,k) = (prai(i,k)+prci(i,k))*icldm(i,k)+(psacws(i,k)+bergs(i,k))*lcldm(i,k)+(pracs(i,k)+mnuccr(i,k))& 1361 *precip_frac(i,k) 1362 ! following are used to calculate 1st order conversion rate of cloud water 1363 ! to rain and snow (1/s), for later use in aerosol wet removal routine 1364 ! previously, wetdepa used (prain/qc) for this, and the qc in wetdepa may be smaller than the qc 1365 ! used to calculate pra, prc, ... in this routine 1366 ! qcsinksum_rate1ord = { rate of direct transfer of cloud water to rain & snow } 1367 ! (no cloud ice or bergeron terms) 1368 qcsinksum_rate1ord(i,k) = (pra(i,k)+prc(i,k)+psacws(i,k))*lcldm(i,k) 1369 ! Avoid zero/near-zero division. 1370 qcsinksum_rate1ord(i,k) = qcsinksum_rate1ord(i,k) / max(qc(i,k),1.0e-30_r8) 1371 ! microphysics output, note this is grid-averaged 1372 pratot(i,k) = pra(i,k)*lcldm(i,k) 1373 prctot(i,k) = prc(i,k)*lcldm(i,k) 1374 mnuccctot(i,k) = mnuccc(i,k)*lcldm(i,k) 1375 mnuccttot(i,k) = mnucct(i,k)*lcldm(i,k) 1376 msacwitot(i,k) = msacwi(i,k)*lcldm(i,k) 1377 psacwstot(i,k) = psacws(i,k)*lcldm(i,k) 1378 bergstot(i,k) = bergs(i,k)*lcldm(i,k) 1379 bergtot(i,k) = berg(i,k) 1380 prcitot(i,k) = prci(i,k)*icldm(i,k) 1381 praitot(i,k) = prai(i,k)*icldm(i,k) 1382 mnuccdtot(i,k) = mnuccd(i,k)*icldm(i,k) 1383 pracstot(i,k) = pracs(i,k)*precip_frac(i,k) 1384 mnuccrtot(i,k) = mnuccr(i,k)*precip_frac(i,k) 1385 nctend(i,k) = nctend(i,k)+ (-nnuccc(i,k)-nnucct(i,k)-npsacws(i,k)+nsubc(i,k) -npra(i,& 1386 k)-nprc1(i,k))*lcldm(i,k) 1387 IF (do_cldice) THEN 1388 IF (use_hetfrz_classnuc) THEN 1389 tmpfrz = nnuccc(i,k) 1390 ELSE 1391 tmpfrz = 0._r8 1392 END IF 1393 nitend(i,k) = nitend(i,k)+ nnuccd(i,k)+ (nnucct(i,k)+tmpfrz+nnudep(i,k)+nsacwi(i,k))& 1394 *lcldm(i,k)+(nsubi(i,k)-nprci(i,k)- nprai(i,k))*icldm(i,k)+nnuccri(i,k)*precip_frac(i,k) 1395 END IF 1396 nstend(i,k) = nstend(i,k)+(nsubs(i,k)+ nsagg(i,k)+nnuccr(i,k))*precip_frac(i,k)+nprci(i,k)*icldm(& 1397 i,k) 1398 nrtend(i,k) = nrtend(i,k)+ nprc(i,k)*lcldm(i,k)+(nsubr(i,k)-npracs(i,k)-nnuccr(i,k) & 1399 -nnuccri(i,k)+nragg(i,k))*precip_frac(i,k) 1400 ! make sure that ni at advanced time step does not exceed 1401 ! maximum (existing N + source terms*dt), which is possible if mtime < deltat 1402 ! note that currently mtime = deltat 1403 !================================================================ 1404 IF (do_cldice .and. nitend(i,k).gt.0._r8.and.ni(i,k)+nitend(i,k)*deltat.gt.nimax(i,k)) THEN 1405 nitend(i,k) = max(0._r8,(nimax(i,k)-ni(i,k))/deltat) 1406 END IF 1407 END DO 1408 ! End of "administration" loop 1409 END DO micro_vert_loop ! end k loop 1410 !----------------------------------------------------- 1411 ! convert rain/snow q and N for output to history, note, 1412 ! output is for gridbox average 1413 qrout = qr 1414 nrout = nr * rho 1415 qsout = qs 1416 nsout = ns * rho 1417 ! calculate precip fluxes 1418 ! calculate the precip flux (kg/m2/s) as mixingratio(kg/kg)*airdensity(kg/m3)*massweightedfallspeed(m/s) 1419 ! --------------------------------------------------------------------- 1420 rflx(:,2:) = rflx(:,2:) + (qric*rho*umr*precip_frac) 1421 sflx(:,2:) = sflx(:,2:) + (qsic*rho*ums*precip_frac) 1422 ! calculate n0r and lamr from rain mass and number 1423 ! divide by precip fraction to get in-precip (local) values of 1424 ! rain mass and number, divide by rhow to get rain number in kg^-1 1425 CALL size_dist_param_basic(mg_rain_props, qric, nric, lamr, n0r) 1426 ! Calculate rercld 1427 ! calculate mean size of combined rain and cloud water 1428 CALL calc_rercld(lamr, n0r, lamc, pgam, qric, qcic, ncic, rercld) 1429 ! Assign variables back to start-of-timestep values 1430 ! Some state variables are changed before the main microphysics loop 1431 ! to make "instantaneous" adjustments. Afterward, we must move those changes 1432 ! back into the tendencies. 1433 ! These processes: 1434 ! - Droplet activation (npccn, impacts nc) 1435 ! - Instantaneous snow melting (minstsm/ninstsm, impacts qr/qs/nr/ns) 1436 ! - Instantaneous rain freezing (minstfr/ninstrf, impacts qr/qs/nr/ns) 1437 !================================================================================ 1438 ! Re-apply droplet activation tendency 1439 nc = ncn 1440 nctend = nctend + npccn 1441 ! Re-apply rain freezing and snow melting. 1442 dum_2d = qs 1443 qs = qsn 1444 qstend = qstend + (dum_2d-qs)/deltat 1445 dum_2d = ns 1446 ns = nsn 1447 nstend = nstend + (dum_2d-ns)/deltat 1448 dum_2d = qr 1449 qr = qrn 1450 qrtend = qrtend + (dum_2d-qr)/deltat 1451 dum_2d = nr 1452 nr = nrn 1453 nrtend = nrtend + (dum_2d-nr)/deltat 1454 !............................................................................. 1455 !================================================================================ 1456 ! modify to include snow. in prain & evap (diagnostic here: for wet dep) 1457 nevapr = nevapr + evapsnow 1458 prain = prain + prodsnow 1459 sed_col_loop: DO i=1,mgncol 1460 DO k=1,nlev 1461 ! calculate sedimentation for cloud water and ice 1462 !================================================================================ 1463 ! update in-cloud cloud mixing ratio and number concentration 1464 ! with microphysical tendencies to calculate sedimentation, assign to dummy vars 1465 ! note: these are in-cloud values***, hence we divide by cloud fraction 1466 dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat)/lcldm(i,k) 1467 dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat)/icldm(i,k) 1468 dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat)/lcldm(i,k),0._r8) 1469 dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat)/icldm(i,k),0._r8) 1470 dumr(i,k) = (qr(i,k)+qrtend(i,k)*deltat)/precip_frac(i,k) 1471 dumnr(i,k) = max((nr(i,k)+nrtend(i,k)*deltat)/precip_frac(i,k),0._r8) 1472 dums(i,k) = (qs(i,k)+qstend(i,k)*deltat)/precip_frac(i,k) 1473 dumns(i,k) = max((ns(i,k)+nstend(i,k)*deltat)/precip_frac(i,k),0._r8) 1474 ! switch for specification of droplet and crystal number 1475 IF (nccons) THEN 1476 dumnc(i,k) = ncnst/rho(i,k) 1477 END IF 1478 ! switch for specification of cloud ice number 1479 IF (nicons) THEN 1480 dumni(i,k) = ninst/rho(i,k) 1481 END IF 1482 ! obtain new slope parameter to avoid possible singularity 1483 CALL size_dist_param_basic(mg_ice_props, dumi(i,k), dumni(i,k), lami(i,k)) 1484 CALL size_dist_param_liq(mg_liq_props, dumc(i,k), dumnc(i,k), rho(i,k), pgam(i,k), lamc(i,k)) 1485 ! calculate number and mass weighted fall velocity for droplets and cloud ice 1486 !------------------------------------------------------------------- 1487 IF (dumc(i,k).ge.qsmall) THEN 1488 vtrmc(i,k) = acn(i,k)*gamma(4._r8+bc+pgam(i,k))/ (lamc(i,k)**bc*gamma(pgam(i,k)+4._r8)) 1489 fc(k) = g*rho(i,k)*vtrmc(i,k) 1490 fnc(k) = g*rho(i,k)* acn(i,k)*gamma(1._r8+bc+pgam(i,k))/ (lamc(i,k)& 1491 **bc*gamma(pgam(i,k)+1._r8)) 1492 ELSE 1493 fc(k) = 0._r8 1494 fnc(k) = 0._r8 1495 END IF 1496 ! calculate number and mass weighted fall velocity for cloud ice 1497 IF (dumi(i,k).ge.qsmall) THEN 1498 vtrmi(i,k) = min(ain(i,k)*gamma_bi_plus4/(6._r8*lami(i,k)**bi), 1.2_r8*rhof(i,k)) 1499 fi(k) = g*rho(i,k)*vtrmi(i,k) 1500 fni(k) = g*rho(i,k)* min(ain(i,k)*gamma_bi_plus1/lami(i,k)**bi,1.2_r8*rhof(i,k)) 1501 ELSE 1502 fi(k) = 0._r8 1503 fni(k) = 0._r8 1504 END IF 1505 ! fallspeed for rain 1506 CALL size_dist_param_basic(mg_rain_props, dumr(i,k), dumnr(i,k), lamr(i,k)) 1507 IF (lamr(i,k).ge.qsmall) THEN 1508 ! 'final' values of number and mass weighted mean fallspeed for rain (m/s) 1509 unr(i,k) = min(arn(i,k)*gamma_br_plus1/lamr(i,k)**br,9.1_r8*rhof(i,k)) 1510 umr(i,k) = min(arn(i,k)*gamma_br_plus4/(6._r8*lamr(i,k)**br),9.1_r8*rhof(i,k)) 1511 fr(k) = g*rho(i,k)*umr(i,k) 1512 fnr(k) = g*rho(i,k)*unr(i,k) 1513 ELSE 1514 fr(k) = 0._r8 1515 fnr(k) = 0._r8 1516 END IF 1517 ! fallspeed for snow 1518 CALL size_dist_param_basic(mg_snow_props, dums(i,k), dumns(i,k), lams(i,k)) 1519 IF (lams(i,k).ge.qsmall) THEN 1520 ! 'final' values of number and mass weighted mean fallspeed for snow (m/s) 1521 ums(i,k) = min(asn(i,k)*gamma_bs_plus4/(6._r8*lams(i,k)**bs),1.2_r8*rhof(i,k)) 1522 uns(i,k) = min(asn(i,k)*gamma_bs_plus1/lams(i,k)**bs,1.2_r8*rhof(i,k)) 1523 fs(k) = g*rho(i,k)*ums(i,k) 1524 fns(k) = g*rho(i,k)*uns(i,k) 1525 ELSE 1526 fs(k) = 0._r8 1527 fns(k) = 0._r8 1528 END IF 1529 ! redefine dummy variables - sedimentation is calculated over grid-scale 1530 ! quantities to ensure conservation 1531 dumc(i,k) = (qc(i,k)+qctend(i,k)*deltat) 1532 dumnc(i,k) = max((nc(i,k)+nctend(i,k)*deltat),0._r8) 1533 dumi(i,k) = (qi(i,k)+qitend(i,k)*deltat) 1534 dumni(i,k) = max((ni(i,k)+nitend(i,k)*deltat),0._r8) 1535 dumr(i,k) = (qr(i,k)+qrtend(i,k)*deltat) 1536 dumnr(i,k) = max((nr(i,k)+nrtend(i,k)*deltat),0._r8) 1537 dums(i,k) = (qs(i,k)+qstend(i,k)*deltat) 1538 dumns(i,k) = max((ns(i,k)+nstend(i,k)*deltat),0._r8) 1539 IF (dumc(i,k).lt.qsmall) dumnc(i,k) = 0._r8 1540 IF (dumi(i,k).lt.qsmall) dumni(i,k) = 0._r8 1541 IF (dumr(i,k).lt.qsmall) dumnr(i,k) = 0._r8 1542 IF (dums(i,k).lt.qsmall) dumns(i,k) = 0._r8 1543 END DO !!! vertical loop 1544 ! initialize nstep for sedimentation sub-steps 1545 ! calculate number of split time steps to ensure courant stability criteria 1546 ! for sedimentation calculations 1547 !------------------------------------------------------------------- 1548 nstep = 1 + int(max( maxval( fi/pdel(i,:)), maxval(fni/pdel(i,:))) * deltat) 1549 ! loop over sedimentation sub-time step to ensure stability 1550 !============================================================== 1551 DO n = 1,nstep 1552 IF (do_cldice) THEN 1553 falouti = fi * dumi(i,:) 1554 faloutni = fni * dumni(i,:) 1555 ELSE 1556 falouti = 0._r8 1557 faloutni = 0._r8 1558 END IF 1559 ! top of model 1560 k = 1 1561 ! add fallout terms to microphysical tendencies 1562 faltndi = falouti(k)/pdel(i,k) 1563 faltndni = faloutni(k)/pdel(i,k) 1564 qitend(i,k) = qitend(i,k)-faltndi/nstep 1565 nitend(i,k) = nitend(i,k)-faltndni/nstep 1566 ! sedimentation tendency for output 1567 qisedten(i,k) = qisedten(i,k)-faltndi/nstep 1568 dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep 1569 dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep 1570 DO k = 2,nlev 1571 ! for cloud liquid and ice, if cloud fraction increases with height 1572 ! then add flux from above to both vapor and cloud water of current level 1573 ! this means that flux entering clear portion of cell from above evaporates 1574 ! instantly 1575 ! note: this is not an issue with precip, since we assume max overlap 1576 dum1 = icldm(i,k)/icldm(i,k-1) 1577 dum1 = min(dum1,1._r8) 1578 faltndqie = (falouti(k)-falouti(k-1))/pdel(i,k) 1579 faltndi = (falouti(k)-dum1*falouti(k-1))/pdel(i,k) 1580 faltndni = (faloutni(k)-dum1*faloutni(k-1))/pdel(i,k) 1581 ! add fallout terms to eulerian tendencies 1582 qitend(i,k) = qitend(i,k)-faltndi/nstep 1583 nitend(i,k) = nitend(i,k)-faltndni/nstep 1584 ! sedimentation tendency for output 1585 qisedten(i,k) = qisedten(i,k)-faltndi/nstep 1586 ! add terms to to evap/sub of cloud water 1587 qvlat(i,k) = qvlat(i,k)-(faltndqie-faltndi)/nstep 1588 ! for output 1589 qisevap(i,k) = qisevap(i,k)-(faltndqie-faltndi)/nstep 1590 tlat(i,k) = tlat(i,k)+(faltndqie-faltndi)*xxls/nstep 1591 dumi(i,k) = dumi(i,k)-faltndi*deltat/nstep 1592 dumni(i,k) = dumni(i,k)-faltndni*deltat/nstep 1593 END DO 1594 ! units below are m/s 1595 ! sedimentation flux at surface is added to precip flux at surface 1596 ! to get total precip (cloud + precip water) rate 1597 prect(i) = prect(i)+falouti(nlev)/g/real(nstep)/1000._r8 1598 preci(i) = preci(i)+falouti(nlev)/g/real(nstep)/1000._r8 1599 END DO 1600 ! calculate number of split time steps to ensure courant stability criteria 1601 ! for sedimentation calculations 1602 !------------------------------------------------------------------- 1603 nstep = 1 + int(max( maxval( fc/pdel(i,:)), maxval(fnc/pdel(i,:))) * deltat) 1604 ! loop over sedimentation sub-time step to ensure stability 1605 !============================================================== 1606 DO n = 1,nstep 1607 faloutc = fc * dumc(i,:) 1608 faloutnc = fnc * dumnc(i,:) 1609 ! top of model 1610 k = 1 1611 ! add fallout terms to microphysical tendencies 1612 faltndc = faloutc(k)/pdel(i,k) 1613 faltndnc = faloutnc(k)/pdel(i,k) 1614 qctend(i,k) = qctend(i,k)-faltndc/nstep 1615 nctend(i,k) = nctend(i,k)-faltndnc/nstep 1616 ! sedimentation tendency for output 1617 qcsedten(i,k) = qcsedten(i,k)-faltndc/nstep 1618 dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep 1619 dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep 1620 DO k = 2,nlev 1621 dum = lcldm(i,k)/lcldm(i,k-1) 1622 dum = min(dum,1._r8) 1623 faltndqce = (faloutc(k)-faloutc(k-1))/pdel(i,k) 1624 faltndc = (faloutc(k)-dum*faloutc(k-1))/pdel(i,k) 1625 faltndnc = (faloutnc(k)-dum*faloutnc(k-1))/pdel(i,k) 1626 ! add fallout terms to eulerian tendencies 1627 qctend(i,k) = qctend(i,k)-faltndc/nstep 1628 nctend(i,k) = nctend(i,k)-faltndnc/nstep 1629 ! sedimentation tendency for output 1630 qcsedten(i,k) = qcsedten(i,k)-faltndc/nstep 1631 ! add terms to to evap/sub of cloud water 1632 qvlat(i,k) = qvlat(i,k)-(faltndqce-faltndc)/nstep 1633 ! for output 1634 qcsevap(i,k) = qcsevap(i,k)-(faltndqce-faltndc)/nstep 1635 tlat(i,k) = tlat(i,k)+(faltndqce-faltndc)*xxlv/nstep 1636 dumc(i,k) = dumc(i,k)-faltndc*deltat/nstep 1637 dumnc(i,k) = dumnc(i,k)-faltndnc*deltat/nstep 1638 END DO 1639 prect(i) = prect(i)+faloutc(nlev)/g/real(nstep)/1000._r8 1640 END DO 1641 ! calculate number of split time steps to ensure courant stability criteria 1642 ! for sedimentation calculations 1643 !------------------------------------------------------------------- 1644 nstep = 1 + int(max( maxval( fr/pdel(i,:)), maxval(fnr/pdel(i,:))) * deltat) 1645 ! loop over sedimentation sub-time step to ensure stability 1646 !============================================================== 1647 DO n = 1,nstep 1648 faloutr = fr * dumr(i,:) 1649 faloutnr = fnr * dumnr(i,:) 1650 ! top of model 1651 k = 1 1652 ! add fallout terms to microphysical tendencies 1653 faltndr = faloutr(k)/pdel(i,k) 1654 faltndnr = faloutnr(k)/pdel(i,k) 1655 qrtend(i,k) = qrtend(i,k)-faltndr/nstep 1656 nrtend(i,k) = nrtend(i,k)-faltndnr/nstep 1657 ! sedimentation tendency for output 1658 qrsedten(i,k) = qrsedten(i,k)-faltndr/nstep 1659 dumr(i,k) = dumr(i,k)-faltndr*deltat/real(nstep) 1660 dumnr(i,k) = dumnr(i,k)-faltndnr*deltat/real(nstep) 1661 DO k = 2,nlev 1662 faltndr = (faloutr(k)-faloutr(k-1))/pdel(i,k) 1663 faltndnr = (faloutnr(k)-faloutnr(k-1))/pdel(i,k) 1664 ! add fallout terms to eulerian tendencies 1665 qrtend(i,k) = qrtend(i,k)-faltndr/nstep 1666 nrtend(i,k) = nrtend(i,k)-faltndnr/nstep 1667 ! sedimentation tendency for output 1668 qrsedten(i,k) = qrsedten(i,k)-faltndr/nstep 1669 dumr(i,k) = dumr(i,k)-faltndr*deltat/real(nstep) 1670 dumnr(i,k) = dumnr(i,k)-faltndnr*deltat/real(nstep) 1671 END DO 1672 prect(i) = prect(i)+faloutr(nlev)/g/real(nstep)/1000._r8 1673 END DO 1674 ! calculate number of split time steps to ensure courant stability criteria 1675 ! for sedimentation calculations 1676 !------------------------------------------------------------------- 1677 nstep = 1 + int(max( maxval( fs/pdel(i,:)), maxval(fns/pdel(i,:))) * deltat) 1678 ! loop over sedimentation sub-time step to ensure stability 1679 !============================================================== 1680 DO n = 1,nstep 1681 falouts = fs * dums(i,:) 1682 faloutns = fns * dumns(i,:) 1683 ! top of model 1684 k = 1 1685 ! add fallout terms to microphysical tendencies 1686 faltnds = falouts(k)/pdel(i,k) 1687 faltndns = faloutns(k)/pdel(i,k) 1688 qstend(i,k) = qstend(i,k)-faltnds/nstep 1689 nstend(i,k) = nstend(i,k)-faltndns/nstep 1690 ! sedimentation tendency for output 1691 qssedten(i,k) = qssedten(i,k)-faltnds/nstep 1692 dums(i,k) = dums(i,k)-faltnds*deltat/real(nstep) 1693 dumns(i,k) = dumns(i,k)-faltndns*deltat/real(nstep) 1694 DO k = 2,nlev 1695 faltnds = (falouts(k)-falouts(k-1))/pdel(i,k) 1696 faltndns = (faloutns(k)-faloutns(k-1))/pdel(i,k) 1697 ! add fallout terms to eulerian tendencies 1698 qstend(i,k) = qstend(i,k)-faltnds/nstep 1699 nstend(i,k) = nstend(i,k)-faltndns/nstep 1700 ! sedimentation tendency for output 1701 qssedten(i,k) = qssedten(i,k)-faltnds/nstep 1702 dums(i,k) = dums(i,k)-faltnds*deltat/real(nstep) 1703 dumns(i,k) = dumns(i,k)-faltndns*deltat/real(nstep) 1704 END DO !! k loop 1705 prect(i) = prect(i)+falouts(nlev)/g/real(nstep)/1000._r8 1706 preci(i) = preci(i)+falouts(nlev)/g/real(nstep)/1000._r8 1707 END DO !! nstep loop 1708 ! end sedimentation 1709 !ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 1710 ! get new update for variables that includes sedimentation tendency 1711 ! note : here dum variables are grid-average, NOT in-cloud 1712 DO k=1,nlev 1713 dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8) 1714 dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8) 1715 dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8) 1716 dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8) 1717 dumr(i,k) = max(qr(i,k)+qrtend(i,k)*deltat,0._r8) 1718 dumnr(i,k) = max(nr(i,k)+nrtend(i,k)*deltat,0._r8) 1719 dums(i,k) = max(qs(i,k)+qstend(i,k)*deltat,0._r8) 1720 dumns(i,k) = max(ns(i,k)+nstend(i,k)*deltat,0._r8) 1721 ! switch for specification of droplet and crystal number 1722 IF (nccons) THEN 1723 dumnc(i,k) = ncnst/rho(i,k)*lcldm(i,k) 1724 END IF 1725 ! switch for specification of cloud ice number 1726 IF (nicons) THEN 1727 dumni(i,k) = ninst/rho(i,k)*icldm(i,k) 1728 END IF 1729 IF (dumc(i,k).lt.qsmall) dumnc(i,k) = 0._r8 1730 IF (dumi(i,k).lt.qsmall) dumni(i,k) = 0._r8 1731 IF (dumr(i,k).lt.qsmall) dumnr(i,k) = 0._r8 1732 IF (dums(i,k).lt.qsmall) dumns(i,k) = 0._r8 1733 ! calculate instantaneous processes (melting, homogeneous freezing) 1734 !==================================================================== 1735 ! melting of snow at +2 C 1736 IF (t(i,k)+tlat(i,k)/cpp*deltat > snowmelt) THEN 1737 IF (dums(i,k) > 0._r8) THEN 1738 ! make sure melting snow doesn't reduce temperature below threshold 1739 dum = -xlf/cpp*dums(i,k) 1740 IF (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt. snowmelt) THEN 1741 dum = (t(i,k)+tlat(i,k)/cpp*deltat-snowmelt)*cpp/xlf 1742 dum = dum/dums(i,k) 1743 dum = max(0._r8,dum) 1744 dum = min(1._r8,dum) 1745 ELSE 1746 dum = 1._r8 1747 END IF 1748 qstend(i,k) = qstend(i,k)-dum*dums(i,k)/deltat 1749 nstend(i,k) = nstend(i,k)-dum*dumns(i,k)/deltat 1750 qrtend(i,k) = qrtend(i,k)+dum*dums(i,k)/deltat 1751 nrtend(i,k) = nrtend(i,k)+dum*dumns(i,k)/deltat 1752 dum1 = -xlf*dum*dums(i,k)/deltat 1753 tlat(i,k) = tlat(i,k)+dum1 1754 meltsdttot(i,k) = meltsdttot(i,k) + dum1 1755 END IF 1756 END IF 1757 ! freezing of rain at -5 C 1758 IF (t(i,k)+tlat(i,k)/cpp*deltat < rainfrze) THEN 1759 IF (dumr(i,k) > 0._r8) THEN 1760 ! make sure freezing rain doesn't increase temperature above threshold 1761 dum = xlf/cpp*dumr(i,k) 1762 IF (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.rainfrze) THEN 1763 dum = -(t(i,k)+tlat(i,k)/cpp*deltat-rainfrze)*cpp/xlf 1764 dum = dum/dumr(i,k) 1765 dum = max(0._r8,dum) 1766 dum = min(1._r8,dum) 1767 ELSE 1768 dum = 1._r8 1769 END IF 1770 qrtend(i,k) = qrtend(i,k)-dum*dumr(i,k)/deltat 1771 nrtend(i,k) = nrtend(i,k)-dum*dumnr(i,k)/deltat 1772 ! get mean size of rain = 1/lamr, add frozen rain to either snow or cloud ice 1773 ! depending on mean rain size 1774 CALL size_dist_param_basic(mg_rain_props, dumr(i,k), dumnr(i,k), lamr(i,k)) 1775 IF (lamr(i,k) < 1._r8/dcs) THEN 1776 qstend(i,k) = qstend(i,k)+dum*dumr(i,k)/deltat 1777 nstend(i,k) = nstend(i,k)+dum*dumnr(i,k)/deltat 1778 ELSE 1779 qitend(i,k) = qitend(i,k)+dum*dumr(i,k)/deltat 1780 nitend(i,k) = nitend(i,k)+dum*dumnr(i,k)/deltat 1781 END IF 1782 ! heating tendency 1783 dum1 = xlf*dum*dumr(i,k)/deltat 1784 frzrdttot(i,k) = frzrdttot(i,k) + dum1 1785 tlat(i,k) = tlat(i,k)+dum1 1786 END IF 1787 END IF 1788 IF (do_cldice) THEN 1789 IF (t(i,k)+tlat(i,k)/cpp*deltat > tmelt) THEN 1790 IF (dumi(i,k) > 0._r8) THEN 1791 ! limit so that melting does not push temperature below freezing 1792 !----------------------------------------------------------------- 1793 dum = -dumi(i,k)*xlf/cpp 1794 IF (t(i,k)+tlat(i,k)/cpp*deltat+dum.lt.tmelt) THEN 1795 dum = (t(i,k)+tlat(i,k)/cpp*deltat-tmelt)*cpp/xlf 1796 dum = dum/dumi(i,k) 1797 dum = max(0._r8,dum) 1798 dum = min(1._r8,dum) 1799 ELSE 1800 dum = 1._r8 1801 END IF 1802 qctend(i,k) = qctend(i,k)+dum*dumi(i,k)/deltat 1803 ! for output 1804 melttot(i,k) = dum*dumi(i,k)/deltat 1805 ! assume melting ice produces droplet 1806 ! mean volume radius of 8 micron 1807 nctend(i,k) = nctend(i,k)+3._r8*dum*dumi(i,k)/deltat/ (& 1808 4._r8*pi*5.12e-16_r8*rhow) 1809 qitend(i,k) = ((1._r8-dum)*dumi(i,k)-qi(i,k))/deltat 1810 nitend(i,k) = ((1._r8-dum)*dumni(i,k)-ni(i,k))/deltat 1811 tlat(i,k) = tlat(i,k)-xlf*dum*dumi(i,k)/deltat 1812 END IF 1813 END IF 1814 ! homogeneously freeze droplets at -40 C 1815 !----------------------------------------------------------------- 1816 IF (t(i,k)+tlat(i,k)/cpp*deltat < 233.15_r8) THEN 1817 IF (dumc(i,k) > 0._r8) THEN 1818 ! limit so that freezing does not push temperature above threshold 1819 dum = dumc(i,k)*xlf/cpp 1820 IF (t(i,k)+tlat(i,k)/cpp*deltat+dum.gt.233.15_r8) THEN 1821 dum = -(t(i,k)+tlat(i,k)/cpp*deltat-233.15_r8)*cpp/xlf 1822 dum = dum/dumc(i,k) 1823 dum = max(0._r8,dum) 1824 dum = min(1._r8,dum) 1825 ELSE 1826 dum = 1._r8 1827 END IF 1828 qitend(i,k) = qitend(i,k)+dum*dumc(i,k)/deltat 1829 ! for output 1830 homotot(i,k) = dum*dumc(i,k)/deltat 1831 ! assume 25 micron mean volume radius of homogeneously frozen droplets 1832 ! consistent with size of detrained ice in stratiform.F90 1833 nitend(i,k) = nitend(i,k)+dum*3._r8*dumc(i,k)/(4._r8*3.14_r8*1.563e-14_r8* & 1834 500._r8)/deltat 1835 qctend(i,k) = ((1._r8-dum)*dumc(i,k)-qc(i,k))/deltat 1836 nctend(i,k) = ((1._r8-dum)*dumnc(i,k)-nc(i,k))/deltat 1837 tlat(i,k) = tlat(i,k)+xlf*dum*dumc(i,k)/deltat 1838 END IF 1839 END IF 1840 ! remove any excess over-saturation, which is possible due to non-linearity when adding 1841 ! together all microphysical processes 1842 !----------------------------------------------------------------- 1843 ! follow code similar to old 1 scheme 1844 qtmp = q(i,k)+qvlat(i,k)*deltat 1845 ttmp = t(i,k)+tlat(i,k)/cpp*deltat 1846 ! use rhw to allow ice supersaturation 1847 CALL qsat_water(ttmp, p(i,k), esn, qvn) 1848 IF (qtmp > qvn .and. qvn > 0) THEN 1849 ! expression below is approximate since there may be ice deposition 1850 dum = (qtmp-qvn)/(1._r8+xxlv_squared*qvn/(cpp*rv*ttmp**2))/deltat 1851 ! add to output cme 1852 cmeout(i,k) = cmeout(i,k)+dum 1853 ! now add to tendencies, partition between liquid and ice based on temperature 1854 IF (ttmp > 268.15_r8) THEN 1855 dum1 = 0.0_r8 1856 ! now add to tendencies, partition between liquid and ice based on te 1857 !------------------------------------------------------- 1858 ELSE IF (ttmp < 238.15_r8) THEN 1859 dum1 = 1.0_r8 1860 ELSE 1861 dum1 = (268.15_r8-ttmp)/30._r8 1862 END IF 1863 dum = (qtmp-qvn)/(1._r8+(xxls*dum1+xxlv*(1._r8-dum1))**2 *qvn/(cpp*rv*ttmp**2))/deltat 1864 qctend(i,k) = qctend(i,k)+dum*(1._r8-dum1) 1865 ! for output 1866 qcrestot(i,k) = dum*(1._r8-dum1) 1867 qitend(i,k) = qitend(i,k)+dum*dum1 1868 qirestot(i,k) = dum*dum1 1869 qvlat(i,k) = qvlat(i,k)-dum 1870 ! for output 1871 qvres(i,k) = -dum 1872 tlat(i,k) = tlat(i,k)+dum*(1._r8-dum1)*xxlv+dum*dum1*xxls 1873 END IF 1874 END IF 1875 ! calculate effective radius for pass to radiation code 1876 !========================================================= 1877 ! if no cloud water, default value is 10 micron for droplets, 1878 ! 25 micron for cloud ice 1879 ! update cloud variables after instantaneous processes to get effective radius 1880 ! variables are in-cloud to calculate size dist parameters 1881 dumc(i,k) = max(qc(i,k)+qctend(i,k)*deltat,0._r8)/lcldm(i,k) 1882 dumi(i,k) = max(qi(i,k)+qitend(i,k)*deltat,0._r8)/icldm(i,k) 1883 dumnc(i,k) = max(nc(i,k)+nctend(i,k)*deltat,0._r8)/lcldm(i,k) 1884 dumni(i,k) = max(ni(i,k)+nitend(i,k)*deltat,0._r8)/icldm(i,k) 1885 dumr(i,k) = max(qr(i,k)+qrtend(i,k)*deltat,0._r8)/precip_frac(i,k) 1886 dumnr(i,k) = max(nr(i,k)+nrtend(i,k)*deltat,0._r8)/precip_frac(i,k) 1887 dums(i,k) = max(qs(i,k)+qstend(i,k)*deltat,0._r8)/precip_frac(i,k) 1888 dumns(i,k) = max(ns(i,k)+nstend(i,k)*deltat,0._r8)/precip_frac(i,k) 1889 ! switch for specification of droplet and crystal number 1890 IF (nccons) THEN 1891 dumnc(i,k) = ncnst/rho(i,k) 1892 END IF 1893 ! switch for specification of cloud ice number 1894 IF (nicons) THEN 1895 dumni(i,k) = ninst/rho(i,k) 1896 END IF 1897 ! limit in-cloud mixing ratio to reasonable value of 5 g kg-1 1898 dumc(i,k) = min(dumc(i,k),5.e-3_r8) 1899 dumi(i,k) = min(dumi(i,k),5.e-3_r8) 1900 ! limit in-precip mixing ratios 1901 dumr(i,k) = min(dumr(i,k),10.e-3_r8) 1902 dums(i,k) = min(dums(i,k),10.e-3_r8) 1903 ! cloud ice effective radius 1904 !----------------------------------------------------------------- 1905 IF (do_cldice) THEN 1906 IF (dumi(i,k).ge.qsmall) THEN 1907 dum_2d(i,k) = dumni(i,k) 1908 CALL size_dist_param_basic(mg_ice_props, dumi(i,k), dumni(i,k), lami(i,k)) 1909 IF (dumni(i,k) /=dum_2d(i,k)) THEN 1910 ! adjust number conc if needed to keep mean size in reasonable range 1911 nitend(i,k) = (dumni(i,k)*icldm(i,k)-ni(i,k))/deltat 1912 END IF 1913 effi(i,k) = 1.5_r8/lami(i,k)*1.e6_r8 1914 ELSE 1915 effi(i,k) = 25._r8 1916 END IF 1917 ! ice effective diameter for david mitchell's optics 1918 deffi(i,k) = effi(i,k)*rhoi/rhows*2._r8 1919 ELSE 1920 ! NOTE: If CARMA is doing the ice microphysics, then the ice effective 1921 ! radius has already been determined from the size distribution. 1922 effi(i,k) = re_ice(i,k) * 1.e6_r8 ! m -> um 1923 deffi(i,k) = effi(i,k) * 2._r8 1924 END IF 1925 ! cloud droplet effective radius 1926 !----------------------------------------------------------------- 1927 IF (dumc(i,k).ge.qsmall) THEN 1928 ! switch for specification of droplet and crystal number 1929 IF (nccons) THEN 1930 ! make sure nc is consistence with the constant N by adjusting tendency, need 1931 ! to multiply by cloud fraction 1932 ! note that nctend may be further adjusted below if mean droplet size is 1933 ! out of bounds 1934 nctend(i,k) = (ncnst/rho(i,k)*lcldm(i,k)-nc(i,k))/deltat 1935 END IF 1936 dum = dumnc(i,k) 1937 CALL size_dist_param_liq(mg_liq_props, dumc(i,k), dumnc(i,k), rho(i,k), pgam(i,k), lamc(i,k)) 1938 IF (dum /= dumnc(i,k)) THEN 1939 ! adjust number conc if needed to keep mean size in reasonable range 1940 nctend(i,k) = (dumnc(i,k)*lcldm(i,k)-nc(i,k))/deltat 1941 END IF 1942 effc(i,k) = (pgam(i,k)+3._r8)/lamc(i,k)/2._r8*1.e6_r8 1943 !assign output fields for shape here 1944 lamcrad(i,k) = lamc(i,k) 1945 pgamrad(i,k) = pgam(i,k) 1946 ! recalculate effective radius for constant number, in order to separate 1947 ! first and second indirect effects 1948 !====================================== 1949 ! assume constant number of 10^8 kg-1 1950 dumnc(i,k) = 1.e8_r8 1951 ! Pass in "false" adjust flag to prevent number from being changed within 1952 ! size distribution subroutine. 1953 CALL size_dist_param_liq(mg_liq_props, dumc(i,k), dumnc(i,k), rho(i,k), pgam(i,k), lamc(i,k)) 1954 effc_fn(i,k) = (pgam(i,k)+3._r8)/lamc(i,k)/2._r8*1.e6_r8 1955 ELSE 1956 effc(i,k) = 10._r8 1957 lamcrad(i,k) = 0._r8 1958 pgamrad(i,k) = 0._r8 1959 effc_fn(i,k) = 10._r8 1960 END IF 1961 ! recalculate 'final' rain size distribution parameters 1962 ! to ensure that rain size is in bounds, adjust rain number if needed 1963 IF (dumr(i,k).ge.qsmall) THEN 1964 dum = dumnr(i,k) 1965 CALL size_dist_param_basic(mg_rain_props, dumr(i,k), dumnr(i,k), lamr(i,k)) 1966 IF (dum /= dumnr(i,k)) THEN 1967 ! adjust number conc if needed to keep mean size in reasonable range 1968 nrtend(i,k) = (dumnr(i,k)*precip_frac(i,k)-nr(i,k))/deltat 1969 END IF 1970 END IF 1971 ! recalculate 'final' snow size distribution parameters 1972 ! to ensure that snow size is in bounds, adjust snow number if needed 1973 IF (dums(i,k).ge.qsmall) THEN 1974 dum = dumns(i,k) 1975 CALL size_dist_param_basic(mg_snow_props, dums(i,k), dumns(i,k), lams(i,k)) 1976 IF (dum /= dumns(i,k)) THEN 1977 ! adjust number conc if needed to keep mean size in reasonable range 1978 nstend(i,k) = (dumns(i,k)*precip_frac(i,k)-ns(i,k))/deltat 1979 END IF 1980 END IF 1981 END DO ! vertical k loop 1982 DO k=1,nlev 1983 ! if updated q (after microphysics) is zero, then ensure updated n is also zero 1984 !================================================================================= 1985 IF (qc(i,k)+qctend(i,k)*deltat.lt.qsmall) nctend(i,k) = -nc(i,k)/deltat 1986 IF (do_cldice .and. qi(i,k)+qitend(i,k)*deltat.lt.qsmall) nitend(i,k) = -ni(i,k)/deltat 1987 IF (qr(i,k)+qrtend(i,k)*deltat.lt.qsmall) nrtend(i,k) = -nr(i,k)/deltat 1988 IF (qs(i,k)+qstend(i,k)*deltat.lt.qsmall) nstend(i,k) = -ns(i,k)/deltat 1989 END DO 1990 END DO sed_col_loop ! i loop 1991 ! DO STUFF FOR OUTPUT: 1992 !================================================== 1993 ! qc and qi are only used for output calculations past here, 1994 ! so add qctend and qitend back in one more time 1995 qc = qc + qctend*deltat 1996 qi = qi + qitend*deltat 1997 ! averaging for snow and rain number and diameter 1998 !-------------------------------------------------- 1999 ! drout2/dsout2: 2000 ! diameter of rain and snow 2001 ! dsout: 2002 ! scaled diameter of snow (passed to radiation in 1) 2003 ! reff_rain/reff_snow: 2004 ! calculate effective radius of rain and snow in microns for COSP using Eq. 9 of COSP v1.3 manual 2005 WHERE ( qrout .gt. 1.e-7_r8 .and. nrout.gt.0._r8 ) 2006 qrout2 = qrout * precip_frac 2007 nrout2 = nrout * precip_frac 2008 ! The avg_diameter call does the actual calculation; other diameter 2009 ! outputs are just drout2 times constants. 2010 drout2 = avg_diameter(qrout, nrout, rho, rhow) 2011 freqr = precip_frac 2012 reff_rain = 1.5_r8*drout2*1.e6_r8 2013 ELSEWHERE 2014 qrout2 = 0._r8 2015 nrout2 = 0._r8 2016 drout2 = 0._r8 2017 freqr = 0._r8 2018 reff_rain = 0._r8 2019 END WHERE 2020 WHERE ( qsout .gt. 1.e-7_r8 .and. nsout.gt.0._r8 ) 2021 qsout2 = qsout * precip_frac 2022 nsout2 = nsout * precip_frac 2023 ! The avg_diameter call does the actual calculation; other diameter 2024 ! outputs are just dsout2 times constants. 2025 dsout2 = avg_diameter(qsout, nsout, rho, rhosn) 2026 freqs = precip_frac 2027 dsout = 3._r8*rhosn/rhows*dsout2 2028 reff_snow = 1.5_r8*dsout2*1.e6_r8 2029 ELSEWHERE 2030 dsout = 0._r8 2031 qsout2 = 0._r8 2032 nsout2 = 0._r8 2033 dsout2 = 0._r8 2034 freqs = 0._r8 2035 reff_snow = 0._r8 2036 END WHERE 2037 ! analytic radar reflectivity 2038 !-------------------------------------------------- 2039 ! formulas from Matthew Shupe, NOAA/CERES 2040 ! *****note: radar reflectivity is local (in-precip average) 2041 ! units of mm^6/m^3 2042 DO i = 1,mgncol 2043 DO k=1,nlev 2044 IF (qc(i,k).ge.qsmall) THEN 2045 dum = (qc(i,k)/lcldm(i,k)*rho(i,k)*1000._r8)**2 /(0.109_r8*(nc(i,k)+nctend(i,k)*deltat)& 2046 /lcldm(i,k)*rho(i,k)/1.e6_r8)*lcldm(i,k)/precip_frac(i,k) 2047 ELSE 2048 dum = 0._r8 2049 END IF 2050 IF (qi(i,k).ge.qsmall) THEN 2051 dum1 = (qi(i,k)*rho(i,k)/icldm(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8)*icldm(i,k)/precip_frac(i,k) 2052 ELSE 2053 dum1 = 0._r8 2054 END IF 2055 IF (qsout(i,k).ge.qsmall) THEN 2056 dum1 = dum1+(qsout(i,k)*rho(i,k)*1000._r8/0.1_r8)**(1._r8/0.63_r8) 2057 END IF 2058 refl(i,k) = dum+dum1 2059 ! add rain rate, but for 37 GHz formulation instead of 94 GHz 2060 ! formula approximated from data of Matrasov (2007) 2061 ! rainrt is the rain rate in mm/hr 2062 ! reflectivity (dum) is in DBz 2063 IF (rainrt(i,k).ge.0.001_r8) THEN 2064 dum = log10(rainrt(i,k)**6._r8)+16._r8 2065 ! convert from DBz to mm^6/m^3 2066 dum = 10._r8**(dum/10._r8) 2067 ELSE 2068 ! don't include rain rate in R calculation for values less than 0.001 mm/hr 2069 dum = 0._r8 2070 END IF 2071 ! add to refl 2072 refl(i,k) = refl(i,k)+dum 2073 !output reflectivity in Z. 2074 areflz(i,k) = refl(i,k) * precip_frac(i,k) 2075 ! convert back to DBz 2076 IF (refl(i,k).gt.minrefl) THEN 2077 refl(i,k) = 10._r8*log10(refl(i,k)) 2078 ELSE 2079 refl(i,k) = -9999._r8 2080 END IF 2081 !set averaging flag 2082 IF (refl(i,k).gt.mindbz) THEN 2083 arefl(i,k) = refl(i,k) * precip_frac(i,k) 2084 frefl(i,k) = precip_frac(i,k) 2085 ELSE 2086 arefl(i,k) = 0._r8 2087 areflz(i,k) = 0._r8 2088 frefl(i,k) = 0._r8 2089 END IF 2090 ! bound cloudsat reflectivity 2091 csrfl(i,k) = min(csmax,refl(i,k)) 2092 !set averaging flag 2093 IF (csrfl(i,k).gt.csmin) THEN 2094 acsrfl(i,k) = refl(i,k) * precip_frac(i,k) 2095 fcsrfl(i,k) = precip_frac(i,k) 2096 ELSE 2097 acsrfl(i,k) = 0._r8 2098 fcsrfl(i,k) = 0._r8 2099 END IF 2100 END DO 2101 END DO 2102 !redefine fice here.... 2103 dum_2d = qsout + qrout + qc + qi 2104 dumi = qsout + qi 2105 WHERE ( dumi .gt. qsmall .and. dum_2d .gt. qsmall ) 2106 nfice = min(dumi/dum_2d,1._r8) 2107 ELSEWHERE 2108 nfice = 0._r8 2109 END WHERE 2110 END SUBROUTINE micro_mg_tend 2111 !======================================================================== 2112 !OUTPUT CALCULATIONS 2113 !======================================================================== 2114 2115 elemental SUBROUTINE calc_rercld(lamr, n0r, lamc, pgam, qric, qcic, ncic, rercld) 2116 REAL(KIND=r8), intent(in) :: lamr ! rain size parameter (slope) 2117 REAL(KIND=r8), intent(in) :: n0r ! rain size parameter (intercept) 2118 REAL(KIND=r8), intent(in) :: lamc ! size distribution parameter (slope) 2119 REAL(KIND=r8), intent(in) :: pgam ! droplet size parameter 2120 REAL(KIND=r8), intent(in) :: qric ! in-cloud rain mass mixing ratio 2121 REAL(KIND=r8), intent(in) :: qcic ! in-cloud cloud liquid 2122 REAL(KIND=r8), intent(in) :: ncic ! in-cloud droplet number concentration 2123 REAL(KIND=r8), intent(inout) :: rercld ! effective radius calculation for rain + cloud 2124 ! combined size of precip & cloud drops 2125 REAL(KIND=r8) :: atmp 2126 ! Rain drops 2127 IF (lamr > 0._r8) THEN 2128 atmp = n0r * pi / (2._r8 * lamr**3._r8) 2129 ELSE 2130 atmp = 0._r8 2131 END IF 2132 ! Add cloud drops 2133 IF (lamc > 0._r8) THEN 2134 atmp = atmp + ncic * pi * rising_factorial(pgam+1._r8, 2)/(4._r8 * lamc**2._r8) 2135 END IF 2136 IF (atmp > 0._r8) THEN 2137 rercld = rercld + 3._r8 *(qric + qcic) / (4._r8 * rhow * atmp) 2138 END IF 2139 END SUBROUTINE calc_rercld 2140 !======================================================================== 2141 !UTILITIES 2142 !======================================================================== 2143 2144 END MODULE micro_mg2_0 2145