From 06b4dc54f2fcaae76a91658260b36882e45e6c9b Mon Sep 17 00:00:00 2001 From: Ted Mansell Date: Fri, 21 Feb 2025 17:32:46 -0600 Subject: [PATCH 1/8] NSSL-MP scheme update: - More accurate saturation mixing ratio calculation (iqvsopt=1) - Changed default droplet renucleation to irenuc=5, which allows extra nucleation at high supersaturation - A default explicit rain breakup for 3-moment rain (irainbreak=2) has been added to reduce excessive median drop diameters in rain cores, which in turn increases evaporation and cool pool temperature deficits (were too warm). - Imposed reflectivity conservation in graupel->hail conversion (ihlcnh=2) and Bigg freezing (both 2- and 3-moment) - Optional diagnostic outputs for supersaturation wrt to liquid (ssat) and ice (ssati); namelist option nssl_ssat_output --- Registry/Registry.EM_COMMON | 6 + doc/README.NSSLmp | 20 +- dyn_em/solve_em.F | 3 + phys/module_microphysics_driver.F | 9 +- phys/module_mp_nssl_2mom.F | 1913 ++++++++++++++++++++--------- phys/module_physics_init.F | 11 +- 6 files changed, 1406 insertions(+), 556 deletions(-) diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index 89874dc1a9..aacb32b5a6 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -1612,6 +1612,9 @@ state real itype ikj misc 1 - hdu "i state real itype_2 ikj misc 1 - hdu "ice_type2" "Diagnostic ice type cat 2 ISHMAEL microphysics" "" state real itype_3 ikj misc 1 - hdu "ice_type3" "Diagnostic ice type cat 3 ISHMAEL microphysics" "" +state real ssat ikj dyn_em 1 - hdu "ssat" "Supersaturation wrt liquid" "%" +state real ssati ikj dyn_em 1 - hdu "ssati" "Supersaturation wrt ice" "%" + # LIGHTNING NUDGING #state real ltg_dat ij misc 1 - r "ltg_dat" "gridded lightning data" "Flash per xkm x xkm per LAD_INT sec" # END LIGHTNING NUDGING @@ -2422,6 +2425,7 @@ rconfig integer nssl_ccn_is_ccna namelist,physics 1 0 rconfig integer nssl_2moment_on namelist,physics 1 -1 rh "NSSL 2-moment flag" "" "" rconfig integer nssl_3moment namelist,physics 1 0 rh "NSSL 3-moment flag" "" "" rconfig integer nssl_density_on namelist,physics 1 -1 rh "NSSL graupel/hail density flag" "" "" +rconfig integer nssl_ssat_output namelist,physics 1 0 rh "NSSL ssat output flag" "" "" @@ -3055,6 +3059,8 @@ package nssl_hail1m nssl_hail_on==2 - moist:qh; package nssl_ccn_opt nssl_ccn_on==1 - scalar:qnn package nssl_graupelvol nssl_density_on==1 - scalar:qvolg package nssl_hailvol nssl_density_on==2 - scalar:qvolg,qvolh +package nssl_ssat_out nssl_ssat_output==1 - state:ssat +package nssl_ssati_out nssl_ssat_output==2 - state:ssat,ssati package radar_refl compute_radar_ref==1 - state:refl_10cm,refd_max diff --git a/doc/README.NSSLmp b/doc/README.NSSLmp index e9b673653e..31ef9409aa 100644 --- a/doc/README.NSSLmp +++ b/doc/README.NSSLmp @@ -98,14 +98,23 @@ Cloud concentration nuclei (CCN) concentration is predicted as in Mansell et al. Droplet activation option method is controlled by the 'irenuc' option (internal to NSSL module). The default option (2) depletes CCN from the unactivated CCN field. A new option (7) instead counts the number of activated CCN (nucleated droplets) with the assumption of an initial constant CCN number mixing ratio. Option 7 better handles supersaturation at low CCN (e.g., maritime) concentrations by allowing extra droplet activation at high SS. irenuc : (nssl_mp_params namelist) - 2 = ccn field is UNactivated aerosol (default; old droplet activation) + 2 = ccn field is UNactivated aerosol (old default; old droplet activation) Can switch to counting activated CCN with nssl_ccn_is_ccna=1 + 5 = ccn field must be ACTVIATED aerosol (new default as of Feb. 2025) + Must have nssl_ccn_on=1 for irenuc=5 + Allows activation beyond limit of nssl_cccn at higher supersaturation + as an approximation of nucleation mode aerosol being activated. (Mainly + an issue for low CCN concentration with deep updrafts.) + If more strict limitation of activation is desired, use option 7. 7 = ccn field must be ACTVIATED aerosol (new droplet activation) Must have nssl_ccn_on=1 for irenuc=7 Excessive size sorting (common in 2-moment schemes) is effectively controlled by an adaptive breakup method that prevents reflectivity growth by sedimentation (Mansell 2010). For 2-moment, infall=4 (default; nssl_mp_params namelist) is recommended. For 3-moment, infall only really applies to droplets, cloud ice, and snow. +3-moment active rain breakup (WRF 4.7.x, 2025): The 3-moment rain without explicit breakup can result in cold pools that are too warm and rain median diameters that are too large in rain cores. A bin-model-based breakup parameterization for rain was implemented to address these issues. Very low rain rates (sparse drops) are largely unaffected (e.g., maintains Zdr arc feature). The breakup coefficient (rainbreakfac) has a default value of 1.0e6 and can reasonably be increased up to around 2.5e6 if desired (nssl_mp_params). Active breakup is automatically turned on for 3-moment (irainbreak=2) but not for 2-moment. Option irainbreak=2 is not recommended for 2-moment, but a user may experiment with irainbreak=11, which breaks up large drops in the tail of the spectrum starting at D=draintail (default 10.e-3 m). + Graupel -> hail conversion: The parameter ihlcnh selects the method of converting graupel (hail embryos) to the hail category. The default value is -1 for automatic setting. The original option (ihlcnh=1) is replaced by a new option (ihlcnh=3) as of May 2023. ihlcnh=3 converts from the graupel spectrum itself based on the wet growth diameter, which generally results in fewer initiated hailstones with larger diameters (and larger mean diameter at the ground). If hail size seems excessive, try setting ihlcnh=1, which tends to generate higher hail number concentrations and thus smaller diameters. +UPDATE (4.7.x/2025): The conversion has been updated to conserve reflectivity of the new hailstones compared to the graupel. This results in new hail that is smaller than previously but prevents spurious increases in reflectivity. (Active for both 2- and 3-moment) The June 2023 (WRF 4.6) update introduces changes in the default options for graupel/hail fall speeds and collection efficiencies. The original fall speed options (icdx=3; icdxhl=3) from Mansell et al. (2010) are switched to the Milbrandt and Morrison (2013) fall speed curves (icdx=6; icdxhl=6). Because the fall speeds are generally a bit lower, a partially compensating increase in maximum collection efficiency is set by default: ehw0/ehlw0 increased to 0.9. One effect is somewhat reduced total precipitation and cold pool intensity for supercell storms. @@ -138,6 +147,15 @@ Snow self-collection (aggregation) has been curbed in the 4.6 version by reducin Snow reflectivity formerly had a default setting that turned on a crude bright band enhancement (iusewetsnow=1). This is now turned off by default (iusewetsnow=0) These snow parameters can be accessed through the nssl_mp_params namelist. +Saturation mixing ratio (WRF 4.7.x, 2025): New formulation (iqvsopt=1) is more consistent with other microphysics schemes. Previously (iqvsopt=0), the quantity e/(p-e) was approximated as e/p, but the new default restores the full equation and uses slightly more accurate (Bolton) coefficients for the saturation (wrt liquid) tables. + +New options (Feb. 2025) (not enabled by default): + - Option (nsplinter=1001) for ice crystal production by drop freezing/shattering (Sullivan et al. 2018) + - Option (incwet = 1) to treat wet growth only for D > Dwet rather than all or nothing; results in greater hail production due to maintaining dry growth at D < Dwet + +New diagnostic output option: + nssl_ssat_output : (default 0); 1 = Supersaturation wrt liquid; 2 = also supersat. wrt ice + References: Mansell, E. R., C. L. Ziegler, and E. C. Bruning, 2010: Simulated electrification diff --git a/dyn_em/solve_em.F b/dyn_em/solve_em.F index 8ce0089897..cafc08b37a 100644 --- a/dyn_em/solve_em.F +++ b/dyn_em/solve_em.F @@ -3729,6 +3729,9 @@ END SUBROUTINE CMAQ_DRIVER & ,P8W=p8w ,P=p_phy ,PI_PHY=pi_phy & & ,RHO=grid%rho ,SPEC_ZONE=grid%spec_zone & & ,SR=grid%sr ,TH=th_phy & + & ,ssat=grid%ssat, ssati=grid%ssati & + & ,f_ssat=( SIZE(grid%ssat) .GT. 1 ) & + & ,f_ssati=( SIZE(grid%ssati) .GT. 1 ) & & ,refl_10cm=grid%refl_10cm & ! hm, 9/22/09 for refl & ,vmi3d=grid%vmi3d & ! for P3 & ,di3d=grid%di3d & ! for P3 diff --git a/phys/module_microphysics_driver.F b/phys/module_microphysics_driver.F index fec0628131..ab0345b8c4 100644 --- a/phys/module_microphysics_driver.F +++ b/phys/module_microphysics_driver.F @@ -111,6 +111,8 @@ SUBROUTINE microphysics_driver( & #endif ,qnwfa2d, qnifa2d, qnbca2d & ! for water/ice-friendly/black carbon aerosols ,qnocbb2d, qnbcbb2d & ! for biomass burning aerosols + ,ssat,ssati & + ,f_ssat,f_ssati & ,refl_10cm & ! HM, 9/22/09, add for refl ,vmi3d & ! for P3 ,di3d & ! for P3 @@ -609,7 +611,7 @@ SUBROUTINE microphysics_driver( & ! ! Optional ! - REAL, OPTIONAL, DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(OUT) :: refl_10cm + REAL, OPTIONAL, DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(OUT) :: refl_10cm,ssat,ssati REAL, OPTIONAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(INOUT) :: & ! for ntu3m qdcn_curr,qtcn_curr,qccn_curr,qrcn_curr,qnin_curr, & ! for ntu3m fi_curr,fs_curr,vi_curr,vs_curr,vg_curr,ai_curr, & ! for ntu3m @@ -728,7 +730,8 @@ SUBROUTINE microphysics_driver( & ,f_qvoli,f_qaoli & ! for Jensen ISHMAEL ,f_qvoli2,f_qaoli2 & ! for Jensen ISHMAEL ,f_qi3,f_qni3,f_qvoli3,f_qaoli3 & ! for Jensen ISHMAEL - ,f_qnwfa, f_qnifa, f_qnbca ! Added by G. Thompson + ,f_qnwfa, f_qnifa, f_qnbca & ! Added by G. Thompson + ,f_ssat,f_ssati LOGICAL, OPTIONAL, INTENT(IN) :: diagflag @@ -2076,6 +2079,8 @@ SUBROUTINE microphysics_driver( & GRPLNCV = GRAUPELNCV, & SR=SR, & dbz = refl_10cm, & + ssat3d = ssat, f_ssat=f_ssat, & + ssati = ssati, f_ssati=f_ssati, & #if ( WRF_CHEM == 1 ) WETSCAV_ON = config_flags%wetscav_onoff == 1, & EVAPPROD=evapprod,RAINPROD=rainprod, & diff --git a/phys/module_mp_nssl_2mom.F b/phys/module_mp_nssl_2mom.F index d89baf3571..980106d236 100644 --- a/phys/module_mp_nssl_2mom.F +++ b/phys/module_mp_nssl_2mom.F @@ -1,6 +1,6 @@ !WRF:MODEL_LAYER:PHYSICS -! prepocessed on "Aug 14 2023" at "16:15:23" +! prepocessed on "Feb 19 2025" at "16:07:57" @@ -71,6 +71,19 @@ ! ! !--------------------------------------------------------------------- +! Feb. 2025 +! - More accurate saturation mixing ratio calculation (iqvsopt=1) +! - Changed default droplet renucleation to irenuc=5, which allows extra nucleation at high supersaturation +! - Default explicit rain breakup for 3-moment (irainbreak=2) +! - Imposed reflectivity conservation in graupel->hail conversion (ihlcnh=3) and Bigg +! freezing (both 2- and 3-moment) +! - Option (nsplinter=1001) for ice crystal production by drop freezing/shattering (Sullivan et al. 2018) +! - Option (incwet = 1) to treat wet growth only for D > Dwet rather than all or nothing; results in +! slightly greater hail production due to maintaining dry growth at D < Dwet +! - Improved logic for sedimentation +! - Separated flushing of small masses into its own subroutine (smallvalues) +! - Some syntax fixes for issues with old versions of gfortran +!--------------------------------------------------------------------- ! Apr. 2023 (WRF-4.6) ! - Update to 3-moment for rain, graupel, and hail ! - Change default graupel/hail fall speeds to icdx/icdxhl=6 (Milbrandt & Morrison 2013) @@ -177,6 +190,7 @@ MODULE module_mp_nssl_2mom public nssl_2mom_driver public nssl_2mom_init + private gamma_sp,gamxinf,GAML02, GAML02d300, GAML02d500, fqvs, fqis private gamma_dp, gamxinfdp, gamma_dpr private delbk, delabk @@ -265,6 +279,7 @@ MODULE module_mp_nssl_2mom integer, private :: itfall = 0 integer, private :: iscfall = 1 integer, private :: irfall = -1 + integer, private :: iifall = 0 integer, private :: isfall = 2 ! default limit with method II (more restrictive) logical, private :: do_accurate_sedimentation = .false. ! if true, recalculate fall speeds on sub time steps; (more expensive) ! if false, reuse fall speeds on multiple steps (can have a noticeable speedup) @@ -338,8 +353,9 @@ MODULE module_mp_nssl_2mom ! (first nucleation is done with a KW sat. adj. step) integer, private :: issfilt = 0 ! flag to turn on filtering of supersaturation field integer, private :: icnuclimit = 0 ! limit droplet nucleation based on Konwar et al. (2012) and Chandrakar et al. (2016) - integer, private :: irenuc = 2 ! =1 to always allow renucleation of droplets within the cloud (do no use, obsolete) + integer, private :: irenuc = 5 ! =1 to always allow renucleation of droplets within the cloud (do no use, obsolete) ! =2 renucleation following Twomey/Cohard&Pinty + ! =5 Similar to 7 but can produce extra activated nuclei from the 'smaller' CCN at higher SS ! =7 New renucleation that requires prediction of the number of activated nuclei ! i.e., not only at cloud base integer, private :: irenuc3d = 0 ! =1 to include horizontal gradient in renucleation of droplets within the cloud @@ -400,6 +416,8 @@ MODULE module_mp_nssl_2mom integer, private :: ierw = 1 ! for single-moment rain (LFO/Z) integer, private :: iehr0c = 0 ! 0 -> no collection for T > 0C; 1 -> turn on collection/shedding for T > 0C integer, private :: iehlr0c = 0 ! 0 -> no collection for T > 0C; 1 -> turn on collection/shedding for T > 0C + real , private :: eiw0 = 0.5 ! constant or max assumed ice-crystal-droplet collection efficiency + real , private :: esw0 = 0.5 ! constant or max assumed snow-droplet collection efficiency real , private :: ehw0 = 0.9 ! 0.5 ! constant or max assumed graupel-droplet collection efficiency real , private :: erw0 = 1.0 ! constant assumed rain-droplet collection efficiency real , private :: ehlw0 = 0.9 ! 0.75 ! constant or max assumed hail-droplet collection efficiency @@ -512,6 +530,8 @@ MODULE module_mp_nssl_2mom real :: sheddiamlg = 10.0e-03 ! diameter of hail to use fwmlarge real :: sheddiam0 = 20.0e-03 ! diameter of hail at which all water is shed + real :: fwmhtmptem = -15. ! temperature at which fwmhtmp fully switches to liquid water only being on large particles + integer :: ifwmhtmptemopt = 1 ! option to use fwmhtmptem (1) or dwet (2) for max liquid at T < 0. integer :: ifwmhopt = 2 ! option for calculating maximum liquid fraction when fwmh and/or fwmhl is set to -1 ! 1 = maximum based on size of maximum mass diameter ! 2 = integrate over spectrum for maximum liquid (experimental) @@ -549,6 +569,8 @@ MODULE module_mp_nssl_2mom real , private :: dwtempmin = 242. ! lowest temperature to allow wet growth conversion to hail real , private :: dwehwmin = 0. ! Minimum ehw to use to find wet growth diameter (if > ehw0, then wet growth diam becomes smaller) real , private :: dg0thresh = 0.15 ! graupel wet growth diameter above which we say do not bother + real , private :: wetgrthtoffset = -1. ! maximum temperature (Celcius) for wet growth (shedding) + real , private :: hailcnvtoffset = -2. ! maximum temperature (Celcius) for hail conversion integer :: ifddenfac = 0 ! = 1 to use density threshold to count FD as GR when converting to HL real :: fddenthresh = 500. ! if ifddenfac > 0, then hail from FD with lower density are considered to come from graupel integer :: icvhl2h = 0 ! allow conversion of hail back to graupel when hail density gets close to minimum allowed @@ -567,8 +589,18 @@ MODULE module_mp_nssl_2mom ! = 1 use mean diameter for breakup ! = 2 use maximum mass diameter for breakup ! = 3 use mass-weighted diameter for breakup - integer :: iraintailbreak = 0 ! 1 = on - real :: draintail = 8.e-3 ! starting size for rain breakup + integer :: irainbreak = -1 ! -1 : auto sets off for 2-moment and on (=2) for 3-moment + ! 0 = off + ! 1 = on (no diameter dependence) (recommend using option 2) + ! 2 = (recommended) as for 1, but apply factor of 1-ec0 to turn off a smaller diameter (ec0 is rain self-coll factor) + ! 10 = as for 1, but sets ec0=1 for rain self-collection (i.e., no passive breakup); set higher rainbreakfac for this option + ! 11 = breakup for DSD tail only; uses draintail etc. + integer :: ibincracr = 0 + real :: rainbreakfac = 1.0e6 ! 1.e6 for irainbreak=2 (reduce double counting); 2.0e6 for lower hand fit for irainbreak=10; 2.542e6 for 'best' fit + real :: draintail = 10.e-3 ! starting size for rain breakup (irainbreak = 11) + real :: drsmall = 1.e-3 ! size of small drops from breakup (irainbreak = 11) + real :: qrbrthresh1 = 0.1e-3 ! lower threshold rain content (kg/m^3) for large drop breakup (irainbreak=11) + real :: qrbrthresh2 = 1.0e-3 ! upper threshold rain content (kg/m^3) for large drop breakup (irainbreak=11) integer, private :: dmrauto = 0 ! = -1 no limiter on crcnw ! = 0 limit crcnw when qr > 1.2*L (Cohard-Pinty 2002) @@ -634,7 +666,7 @@ MODULE module_mp_nssl_2mom real, private :: alphasmlr0 = 14.0 ! shape parameter for drops formed from melting/shedding snow real, private :: delta_alphamlr = 0.5 ! offset from alphamax at which melting does not further collapse the shape parameter - integer :: iqvsopt = 0 ! =0 use old default for tabqvs; =1 use Bolton formulation (Rogers and Yau) + integer, private :: iqvsopt = 1 ! =0 use old default for tabqvs with e/p approx; =1 use Bolton formulation (Rogers and Yau) with e/(p-e) integer :: imaxsupopt = 4 ! how to treat saturation adjustment in two-moment droplets ! 1 = add droplets with same mean mass as current droplets @@ -662,6 +694,8 @@ MODULE module_mp_nssl_2mom integer, private :: lccn = 9 ! 0 or 9, other indices adjusted accordingly integer, private :: lccnuf = 0 integer, private :: lccna = 0 + integer, private :: lccnaco = 0 + integer, private :: lccnanu = 0 integer, private :: lcina = 0 integer, private :: lcin = 0 integer, private :: lnc = 9 @@ -791,6 +825,9 @@ MODULE module_mp_nssl_2mom real, private :: delqxw = 1.0e-10! 1.0e-12 ! real :: tindmn = 233, tindmx = 298.0 ! min and max temperatures where inductive charging is allowed + integer, private :: imorrgdnglimit = 0 ! flag to impose limit on graupel slope parameter + real, private :: morrdnglimit = 2000.E-6 + ! ! gamma function lookup table ! @@ -834,6 +871,7 @@ MODULE module_mp_nssl_2mom integer lvol(lc:lqmx) integer lz(lc:lqmx) integer lliq(li:lqmx) + integer linfall(lc:lqmx) integer denscale(lc:lqmx) ! flag for density scaling (mixing ratio conversion) integer ido(lc:lqmx) @@ -888,11 +926,11 @@ MODULE module_mp_nssl_2mom real xvfmn, xvfmx ! min, max frozen drop volumes real xvgmn, xvgmx ! min, max graupel volumes real xvhmn, xvhmn0, xvhmx, xvhmx0 ! min, max hail volumes - real xvhlmn, xvhlmx ! min, max lg hail volumes + real xvhlmn, xvhlmx, xvhlmx0 ! min, max lg hail volumes - real, parameter :: dhlmn = 0.3e-3, dhlmx = 40.e-3 + real, parameter :: dhlmn = 0.3e-3 real, parameter :: dhmn0 = 0.3e-3 - real, private :: dhmn = dhmn0, dhmx = -1. + real, private :: dhmn = dhmn0, dhmx = -1., dhlmx = -1. ! 40.e-3 real, parameter :: cwradn = 2.0e-6, xcradmn = cwradn ! minimum radius real, parameter :: cwradx = 60.e-6, xcradmx = cwradx ! maximum radius @@ -915,7 +953,7 @@ MODULE module_mp_nssl_2mom parameter( xvfmn=0.523599*(0.1e-3)**3, xvfmx=0.523599*(10.e-3)**3 ) ! mks xvfmx = (pi/6)*(10mm)**3 parameter( xvgmn=0.523599*(0.1e-3)**3, xvgmx=0.523599*(10.e-3)**3 ) ! mks xvfmx = (pi/6)*(10mm)**3 parameter( xvhmn0=0.523599*(0.3e-3)**3, xvhmx0=0.523599*(20.e-3)**3 ) ! mks xvfmx = (pi/6)*(20mm)**3 - parameter( xvhlmn=0.523599*(dhlmn)**3, xvhlmx=0.523599*(dhlmx)**3 ) ! mks xvfmx = (pi/6)*(40mm)**3 + parameter( xvhlmn=0.523599*(dhlmn)**3, xvhlmx0=0.523599*(40.e-3)**3 ) ! mks xvfmx = (pi/6)*(40mm)**3 ! ! electrical permitivity of air C / (N m**2) - check the units @@ -941,6 +979,7 @@ MODULE module_mp_nssl_2mom real, parameter :: cbwbolton = 29.65 ! constants for Bolton formulation real, parameter :: cawbolton = 17.67 + real, parameter :: esbolton = 6.112e2 real, parameter :: tfrh = 233.15 real, parameter :: tfr = 273.15 @@ -992,14 +1031,50 @@ MODULE module_mp_nssl_2mom logical, parameter :: do_satadj_for_wrfchem = .true. - integer, parameter :: ac_opt = 0 ! option flag for alternate aerosol (for NUWRF only) + integer, private :: lcn_nu = 0 ! 27 ! need to check no conflict with other variables + integer, private :: lcn_ac = 0 ! 28 + integer, private :: lcn_co = 0 ! 29 + integer, private :: lcinp = 0 ! 30 + integer, private :: ac_opt = 0 ! option flag for: (1 and 2 currently for NUWRF only) + ! 0 : normal NSSL CCN physics + ! 1 : accumulation mode CN following Fridland et al. (2012, 2017), + ! where CCN number is sum of unactivated CCN and droplet concentrations + ! 2 : As for 1 but have three modes (but does not partition activated CCN) + ! 11: As for 1 but track activated CCN as a separate category (CN category advects only) + ! 22: As for 11 but 3 modes, each with its own activation tracer + real, private :: ac_wthresh = 10.0 ! for W < ac_wthresh, use max of sswater and diagnosed SS; otherwise use sswater logical, private :: nuaccoinp = .false. +! T.Iguchi Y2021 Update +! logical :: ac_only = .true. ! flag for considering ac_mode of CN only, or all nu,ac,co modes (still under construction) + + logical, private :: arg_para = .true. ! flag for Abdul-Razzak_and_Ghan parameterization works similarly to flag_qndrop, and neglects irenuc, ccna(mgs), and cnuc(mgs) + real, private :: nu_pmr = 7.5 * 1.e-3 * 1.e-6 ! aerosol radius (meter); these parameter values follow Cheng et al. (2007QJ) + real, private :: nu_pgw = 0.53 ! Unlike original Abdul-Razzak_and_Ghan, this value is used without log (Cheng et al. 2007QJ) + real, private :: nu_kappa = 0.07 ! ammonium sulfate as CCN (Petters and Kreidenweis, 2007ACP) + real, private :: ac_pmr = 3.8 * 1.e-2 * 1.e-6 ! aerosol radius (meter) + real, private :: ac_pgw = 0.69 + real, private :: ac_kappa = 0.61 ! ammonium sulfate as CCN (Petters and Kreidenweis, 2007ACP) + real, private :: co_pmr = 0.51 * 1.e-6 ! aerosol radius (meter) + real, private :: co_pgw = 0.77 + real, private :: co_kappa = 0.61 ! ammonium sulfate as CCN (Petters and Kreidenweis, 2007ACP) + + real, parameter :: cn_minlimit = 1.e3 ! 1.e3 m-3 = 0.001 cm-3 + + logical :: dm15_para = .false. ! flag for DeMott et al. (2015) parameterization for heterogenous freezing, regardless of "ibfc" + + ! Note to users: Many of these options are for development and not guaranteed to perform well. ! Some may not be functional depending on the version of the code. ! Some may be useful for ensemble physics diversity. Feel free to contact Ted Mansell if you have questions ! in that regard. NAMELIST /nssl_mp_params/ & +! nuwrf 3-mode params + ac_opt,arg_para, & + ac_kappa, ac_pmr, ac_pgw, & + nu_kappa, nu_pmr, nu_pgw, & + co_kappa, co_pmr, co_pgw, & +! --- ndebug, ncdebug,& iusewetgraupel, & iusewethail, & @@ -1007,16 +1082,17 @@ MODULE module_mp_nssl_2mom idbzci, & vtmaxsed, & itfall,iscfall, & - infall,irfall,isfall, & + infall,irfall,isfall,iifall, & rssflg, & sssflg, & hssflg, & hlssflg, & + irainbreak, rainbreakfac, & irimdenopt,rimdenvwgt, & rimc1, rimc2, rimc3, rimc4, & idiagnosecnu, & icnuclimit, & - irenuc, & + irenuc, ccn, & restoreccn, ccntimeconst, cck, & decayufccn, ufccntimeconst, & switchccn, old_cccn, & @@ -1078,6 +1154,7 @@ MODULE module_mp_nssl_2mom ehimax, & ehsmax, & ecollmx, & + eiw0, esw0, & ehw0, ehlw0, & ehr0, ehlr0, & erw0, & @@ -1087,7 +1164,7 @@ MODULE module_mp_nssl_2mom iqcinit, & ssmxinit, & xvdmx, & - dhmn, dhmx, & + dhmn, dhmx, dhlmx, & fwms,fwmh,fwmhl, & ifwmhopt, & ihxw2rain, & @@ -1103,7 +1180,8 @@ MODULE module_mp_nssl_2mom rescale_low_alphah, & rescale_low_alphahl, & rescale_high_alpha, & - ihlcnh, hldia1,iusedw, dwehwmin, dwmin, dwmax, dwtempmin, dg0thresh, & + ihlcnh, hldia1,iusedw, dwehwmin, dwmin, dwmax, dwtempmin, dg0thresh, incwet, & + wetgrthtoffset, hailcnvtoffset, & icvhl2h, hldnmn,hdnmn, & hlcnhdia, hlcnhqmin, & isedonly, & @@ -1159,6 +1237,8 @@ REAL FUNCTION fqis(t) END FUNCTION fqis +!==========================================================================================! + @@ -1167,6 +1247,7 @@ END FUNCTION fqis ! ##################################################################### SUBROUTINE nssl_2mom_init( & & ims,ime, jms,jme, kms,kme, nssl_params, ipctmp, mixphase,ihvol,idoniconlytmp, & + & namelist_filename, & & nssl_graupelfallfac, & & nssl_hailfallfac, & & nssl_ehw0, & @@ -1199,14 +1280,19 @@ SUBROUTINE nssl_2mom_init( & & nssl_alphahl, & & nssl_alphar integer, intent(in), optional :: & - & nssl_icdx, & + & nssl_icdx, & & nssl_icdxhl, myrank, mpiroot, & - & nssl_ufccn - logical, intent(in), optional :: nssl_density_on, nssl_hail_on, nssl_ccn_on, nssl_icecrystals_on + & nssl_ufccn, & + & nssl_ccn_on + logical, intent(in), optional :: nssl_density_on, nssl_hail_on, nssl_icecrystals_on integer, intent(inout), optional :: ccn_is_ccna integer, intent(in),optional :: infileunit + integer,parameter::strsize=512 + character(len=strsize), intent(in), optional :: namelist_filename + character(len=strsize) :: namelist_inputfile + integer, intent(in), optional :: ims,ime, jms,jme, kms,kme real, intent(in), dimension(20), optional :: nssl_params @@ -1223,7 +1309,7 @@ SUBROUTINE nssl_2mom_init( & integer :: hail_on = -1, density_on = -1, icecrystals_on = 1 integer :: ccn_on = -1 - double precision :: arg + double precision :: arg,cwch real :: temq integer :: igam integer :: i,il,j,l @@ -1294,7 +1380,7 @@ SUBROUTINE nssl_2mom_init( & ! hack to switch CCN field to CCNA (activated ccn) ! invertccn = .true. turn_on_ccna = .true. - irenuc = 7 + irenuc = 5 ENDIF ccnuf = Abs( nssl_params(14) ) IF ( present(nssl_ufccn) ) iufccn = nssl_ufccn @@ -1344,19 +1430,34 @@ SUBROUTINE nssl_2mom_init( & ENDIF ENDIF - + ! turn on active rain breakup by default for 3-moment rain since it has no implicit breakup from sedimentation + ! Check this after namelist read so that user can set irainbreak=0 to turn off + IF ( irainbreak == -1 ) THEN + IF ( ipconc >= 6 ) THEN + irainbreak = 2 + ELSE + irainbreak = 0 + ENDIF + ENDIF + + namelist_inputfile = 'namelist.input' ! default for WRF/cm1 + IF ( present( namelist_filename ) ) THEN + namelist_inputfile = namelist_filename + ELSE + ENDIF + IF ( .true. ) THEN ! set to true to enable internal namelist read - open(15,file='namelist.input',status='old',form='formatted',action='read') + open(15,file=namelist_inputfile,status='old',form='formatted',action='read') rewind(15) read(15,NML=nssl_mp_params,iostat=istat) close(15) IF ( istat /= 0 ) THEN #ifdef WRF_ELEC IF ( wrf_dm_on_monitor() ) THEN - write(0,*) 'NSSL_2MOM_INIT: PROBLEM WITH NSSL_MP_PARAMS namelist: not found or bad token' + write(0,*) 'NSSL_2MOM_INIT: NSSL_MP_PARAMS namelist: not found or bad token' ENDIF #else ! write(0,*) 'NSSL_2MOM_INIT: PROBLEM WITH NSSL_MP_PARAMS namelist: not found or bad token' @@ -1371,7 +1472,6 @@ SUBROUTINE nssl_2mom_init( & ENDIF - IF ( iufccn > 0 ) THEN ! make sure to use option that uses UF ccn irenuc = 7 IF ( ccnuf <= 0.0 ) decayufccn = .true. ! assume surface emission and need decay @@ -1382,7 +1482,7 @@ SUBROUTINE nssl_2mom_init( & ENDIF IF ( present( nssl_ccn_on ) ) THEN - IF ( nssl_ccn_on ) THEN + IF ( nssl_ccn_on > 0 ) THEN ccn_on = 1 ELSE ccn_on = 0 @@ -1393,7 +1493,7 @@ SUBROUTINE nssl_2mom_init( & IF ( irenuc >= 5 ) THEN turn_on_ccna = .true. IF ( present( nssl_ccn_on ) ) THEN - IF ( .not. nssl_ccn_on ) THEN + IF ( nssl_ccn_on <= 0 ) THEN write(0,*) 'NSSL_MP Error: Must have nssl_ccn_on=1 for irenuc >= 5!' STOP ENDIF @@ -1470,13 +1570,14 @@ SUBROUTINE nssl_2mom_init( & dtabqvs(l) = ((-caw*(-273.15 + temq))/(temq - cbw)**2 + & & caw/(temq - cbw))*tabqvs(l) ELSE - tabqvs(l) = exp(caw*(temq-273.15)/(temq-cbw)) + tabqvs(l) = exp(cawbolton*(temq-273.15)/(temq-cbwbolton)) dtabqvs(l) = ((-cawbolton*(-273.15 + temq))/(temq - cbwbolton)**2 + & & cawbolton/(temq - cbwbolton))*tabqvs(l) ENDIF tabqis(l) = exp(cai*(temq-273.15)/(temq-cbi)) dtabqis(l) = ((-cai*(-273.15 + temq))/(temq - cbi)**2 + & & cai/(temq - cbi))*tabqis(l) + end do bx(lr) = 0.85 @@ -1528,8 +1629,8 @@ SUBROUTINE nssl_2mom_init( & gmoi(igam) = gamma_dp(arg) end do - ! build lookup table to compute the number and mass fractions of rain drops - ! (imurain=1) greater than a given diameter. Used for qiacr and ciacr + ! build lookup table to compute the number and mass fractions of particles + ! (mu=1) greater than a given diameter. Used for qiacr and ciacr ! Uses incomplete gamma functions ! The terms with bxh or bxhl will be off if the actual bxh or bxhl is different from the base value (icdx=6 option) @@ -1617,6 +1718,8 @@ SUBROUTINE nssl_2mom_init( & lccn = 0 lccnuf = 0 lccna = 0 + lccnaco = 0 + lccnanu = 0 lnc = 0 lnr = 0 lni = 0 @@ -1780,7 +1883,8 @@ SUBROUTINE nssl_2mom_init( & !debug write(0,*) 'Setting lcin to ',lcin ENDIF na = ltmp - + + ln(:) = 0 ln(lc) = lnc ln(lr) = lnr ln(li) = lni @@ -1788,6 +1892,7 @@ SUBROUTINE nssl_2mom_init( & ln(lh) = lnh IF ( lhl .gt. 1 ) ln(lhl) = lnhl + ipc(:) = 0 ipc(lc) = 2 ipc(lr) = 3 ipc(li) = 1 @@ -1990,9 +2095,19 @@ SUBROUTINE nssl_2mom_init( & ido(lh) = idohw IF ( lhl .gt. 1 ) ido(lhl) = idohl + linfall(:) = infall + linfall(lc) = 0 IF ( irfall .lt. 0 ) irfall = infall IF ( isfall .lt. 0 ) isfall = infall + IF ( iifall .lt. 0 ) iifall = infall IF ( lzr > 0 ) irfall = 0 + IF ( lzs > 0 ) isfall = 0 + IF ( lzh > 0 ) linfall(lh) = 0 + IF ( lzhl > 0 .and. lhl > 0 ) linfall(lhl) = 0 + IF ( lzr > 0 .and. lf > 0 ) linfall(lf) = 0 + linfall(lr) = irfall + linfall(ls) = isfall + linfall(li) = iifall qccn = ccn/rho00 qccnuf = ccnuf/rho00 @@ -2023,6 +2138,19 @@ SUBROUTINE nssl_2mom_init( & ELSE xvhmx = 0.523599*(dhmx)**3 ENDIF + + IF ( dhlmx <= 0.0 ) THEN + xvhlmx = xvhlmx0 + ELSE + xvhlmx = 0.523599*(dhlmx)**3 + ENDIF + + IF ( ipconc == 5 .and. imorrgdnglimit >= 1 ) THEN + ! convert morrdnglimit to xvhmx equivalent + cwch = ((3. + alphah)*(2. + alphah)*(1.0 + alphah))**(-1./3.) + xvhmx = pi/6.0*(morrdnglimit/cwch)**3 + dhmx = morrdnglimit/cwch + ENDIF IF ( qhdpvdn < 0. ) qhdpvdn = xdnmn(lh) IF ( qhacidn < 0. ) qhacidn = xdnmn(lh) @@ -2224,7 +2352,9 @@ END SUBROUTINE nssl_2mom_init SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw, chl, & cn, vhw, vhl, cna, cni, f_cn, f_cna, f_cina, & f_qc, f_qr, f_qi, f_qs, f_qh, f_qhl, & - cnuf, f_cnuf, & + cn_nu, cn_co, cinp, f_cnnu, f_cnco, f_cinp, & + cna_co, cna_nu, f_cnaco, f_cnanu, & + cnuf, f_cnuf, cn_ac, f_cnac, & zrw, zhw, zhl, f_zrw, f_zhw, f_zhl, f_vhw, f_vhl, & qsw, qhw, qhlw, & tt, th, pii, p, w, dn, dz, dtp, itimestep, & @@ -2246,8 +2376,9 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw rscghis_2d,rscghis_2dp,rscghis_2dn, & scr,scw,sci,scs,sch,schl,sctot, & elec_physics, & - induc,elecz,scion,sciona, & + induc,elecz,scion,sciona,f_scion,f_sciona, & noninduc,noninducp,noninducn, & + ssat3d,ssati,f_ssat,f_ssati, & pcc2, pre2, depsubr, & mnucf2, melr2, ctr2, & rim1_2, rim2_2,rim3_2, & @@ -2261,6 +2392,7 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw ! re_liquid, re_graupel, re_hail, re_icesnow, & ! vtcloud, vtrain, vtsnow, vtgraupel, vthail, & ipelectmp, & + isedonly_in, & diagflag,ke_diag, & nssl_progn, & ! wrf-chem ! 20130903 acd_mb_washout start @@ -2294,7 +2426,11 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw qi,qhl,ccw,crw,cci,csw,chw,chl,vhw,vhl integer, optional, intent(in) :: is_theta_or_temp logical, optional, intent(in) :: f_zrw, f_zhw, f_zhl, f_vhw, f_vhl ! not used yet + logical, optional, intent(in) :: f_ssat,f_ssati real, dimension(ims:ime, kms:kme, jms:jme), optional, intent(inout):: dbz, vzf, cn, cna, cni, cnuf + real, dimension(ims:ime, kms:kme, jms:jme), optional, intent(inout):: cn_nu, cn_ac, cn_co, cinp, cna_co, cna_nu + logical, optional, intent(in) :: f_cnnu, f_cnac, f_cnco, f_cinp, f_cnaco, f_cnanu + real, dimension(ims:ime, jms:jme), optional, intent(inout):: compdbz real, dimension(ims:ime, jms:jme), optional, intent(inout):: rscghis_2d, & ! 2D accumulation arrays for vertically-integrated charging rate rscghis_2dp, & ! 2D accumulation arrays for vertically-integrated charging rate (positive only) @@ -2306,11 +2442,12 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw real, dimension(ims:ime, kms:kme, jms:jme), optional, intent(inout):: & induc,noninduc,noninducp,noninducn ! charging rates: inductive, noninductive (all, positive, negative to graupel) real, dimension(ims:ime, kms:kme, jms:jme), optional, intent(in) :: elecz ! elecsave = Ez - real, dimension(ims:ime, kms:kme, jms:jme,2),optional, intent(inout) :: scion + real, dimension(ims:ime, kms:kme, jms:jme, 2),optional, intent(inout) :: scion real, dimension(ims:ime, kms:kme, jms:jme), intent(in):: p,w,dz,dn real, dimension(ims:ime, kms:kme, jms:jme), intent(in):: pii real, dimension(ims:ime, kms:kme, jms:jme), optional, intent(inout):: & + ssat3d, ssati, & pcc2, pre2, depsubr, & mnucf2, melr2, ctr2, & rim1_2, rim2_2,rim3_2, & @@ -2345,13 +2482,14 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw real, dimension(ims:ime, jms:jme), intent(out), optional :: & rainncw2, rainnci2 ! liquid rain, ice, accumulation rates real, optional, intent(in) :: dx,dy - real, intent(in):: dtp - integer, intent(in):: itimestep !, ccntype + real, intent(in) :: dtp + integer, intent(in) :: itimestep !, ccntype integer, intent(in), optional :: ntmul, ntcnt logical, optional, intent(in) :: lastloop logical, optional, intent(in) :: diagflag, f_cna, f_cn, f_cina, f_cnuf logical, optional, intent(in) :: f_qc, f_qr, f_qi, f_qs, f_qh, f_qhl - integer, optional, intent(in) :: ipelectmp, ke_diag + logical, optional, intent(in) :: f_scion,f_sciona + integer, optional, intent(in) :: ipelectmp, ke_diag, isedonly_in LOGICAL, INTENT(IN), OPTIONAL :: nssl_progn ! flags for wrf-chem @@ -2423,16 +2561,18 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw double precision :: wvol5,wvol10 real :: tmp,dv,dv1,tmpchg real :: rdt + real :: temp1, c1 double precision :: dt1,dt2 double precision :: timesed,timesed1,timesed2,timesed3, timegs, timenucond, timedbz,zmaxsed double precision :: timevtcalc,timesetvt - logical :: f_cnatmp, f_cinatmp + logical :: f_cnatmp, f_cinatmp, f_cnacotmp, f_cnanutmp logical :: has_wetscav integer :: kediagloc integer :: iunit + integer :: isedonly_local real :: ycent, y, emissrate, emissrate0, emissrate1, z, fac, factot real :: fach(kts:kte) @@ -2520,6 +2660,13 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw ELSE ipelec = 0 ENDIF + + IF ( present( isedonly_in ) ) THEN + isedonly_local = isedonly_in + ELSE + isedonly_local = 0 + ENDIF + ! IF ( present( dbz ) ) THEN ! DO jy = jts,jte ! DO kz = kts,kte @@ -2643,10 +2790,9 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw infdo = 0 ENDIF - IF ( infall .ge. 3 .or. ipconc .ge. 6 ) THEN + IF ( Any(linfall(:) .ge. 3 ) .or. ipconc .ge. 6 ) THEN infdo = 2 ENDIF - IF ( present( HAILNCV ) .and. lhl < 1 ) THEN ! for WRF 3.1 compatibility HAILNCV(its:ite,jts:jte) = 0. @@ -2689,7 +2835,7 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw ! write(0,*) 'N2M: load an, jy,lccn = ',jy,lccn,qccn - IF ( present( pcc2 ) .and. makediag ) THEN + IF ( ( present( pcc2 ) .or. present( axtra ) ) .and. makediag ) THEN axtra2d(its:ite,1,kts:kte,:) = 0.0 ENDIF @@ -2749,7 +2895,7 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw IF ( lccna > 1 ) THEN IF ( present( cna ) .and. f_cnatmp ) THEN - an(ix,1,kz,lccna) = cna(ix,kz,jy) + an(ix,1,kz,lccna) = Max(0.0, cna(ix,kz,jy) ) ENDIF ENDIF @@ -2835,7 +2981,11 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw ! ! saturation mixing ratio ! - t8s = t00(ix,1,kz)*tabqvs(ltemq) !saturation mixing ratio wrt water + IF ( iqvsopt == 0 ) THEN + t8s = t00(ix,1,kz)*tabqvs(ltemq) !saturation mixing ratio wrt water + ELSE + t8s = rdorv*esbolton*tabqvs(ltemq)/(pn(ix,1,kz) - esbolton*tabqvs(ltemq)) + ENDIF t9s = t00(ix,1,kz)*tabqis(ltemq) !saturation mixing ratio wrt ice ! @@ -2952,12 +3102,10 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw ! IF ( .true. ) THEN -! #ifndef CM1 ! for real cases when hydrometeor mixing ratios have been initialized without concentrations IF ( itimestep == 1 .and. ipconc > 0 .and. loopcnt == 1 ) THEN call calcnfromq(nx,ny,nz,an,na,nor,nor,dn1) ENDIF -! #endif IF ( present(cu_used) .and. & ( present( qrcuten ) .or. present( qscuten ) .or. & @@ -2977,24 +3125,13 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw call calcnfromcuten(nx,ny,nz,ancuten,an,na,nor,nor,dn1) - DO kz = kts,kte - DO ix = its,ite - - - IF ( ipconc >= 6 ) THEN -! IF ( lzr > 0 ) an(ix,1,kz,lzr) = an(ix,1,kz,lzr) + ancuten(ix,1,kz,lzr) - ENDIF - - ENDDO - ENDDO - ENDIF !} ENDIF !} - - + IF ( isedonly_local == 0 ) THEN + call sediment1d(dtp,nx,ny,nz,an,na,nor,nor,xfall,dn1,dz2d,dz2dinv, & & t0,t7,infdo,jy,its,jts & & ,timesed1,timesed2,timesed3,zmaxsed,timesetvt) @@ -3015,11 +3152,11 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw & xfall(ix,1,lh)*1000./xdn0(lr) ) ENDIF IF ( present ( rainncw2 ) ) THEN ! rain only - rainncw2(ix,jy) = rainncw2(ix,jy) + dtp*dn1(ix,1,1)*xfall(ix,1,lr) + rainncw2(ix,jy) = rainncw2(ix,jy) + dtp*dn1(ix,1,1)*xfall(ix,1,lr) ENDIF IF ( present ( rainnci2 ) ) THEN ! ice only IF ( lhl > 1 ) THEN - rainnci2(ix,jy) =rainnci2(ix,jy) + dtp*dn1(ix,1,1)*(xfall(ix,1,ls)*1000./xdn0(lr) + & + rainnci2(ix,jy) = rainnci2(ix,jy) + dtp*dn1(ix,1,1)*(xfall(ix,1,ls)*1000./xdn0(lr) + & & xfall(ix,1,lh)*1000./xdn0(lr) + xfall(ix,1,lhl)*1000./xdn0(lr) ) ELSE rainnci2(ix,jy) = rainnci2(ix,jy) + dtp*dn1(ix,1,1)*(xfall(ix,1,ls)*1000./xdn0(lr) + & @@ -3040,11 +3177,7 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw SNOWNC(ix,jy) = SNOWNC(ix,jy) + SNOWNCV(ix,jy) ENDIF IF ( lhl > 1 ) THEN -!#ifdef CM1 -! IF ( .true. ) THEN -!#else IF ( present( HAILNC ) ) THEN -!#endif HAILNCV(ix,jy) = dtp*dn1(ix,1,1)*xfall(ix,1,lhl)*1000./xdn0(lr) IF ( loopcnt == loopmax ) HAILNC(ix,jy) = HAILNC(ix,jy) + HAILNCV(ix,jy) ! ELSEIF ( present( GRPLNCV ) ) THEN ! if no separate hail accum, then add to graupel @@ -3063,6 +3196,8 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw ENDIF ENDDO + ENDIF ! isedonly_local + ! ENDIF ! .false. IF ( isedonly /= 1 ) THEN @@ -3110,6 +3245,14 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw & ,axtra2d, makediag & & ,ssat,t00,t77,flag_qndrop) +! Clean up tiny values of mixing ratio and final checks on max/min sizes + CALL smallvalues & + & (nx,ny,nz,na,jy & + & ,nor,nor,dtp,nx & + & ,t0 & + & ,an,dn1,wn & + & ,t77,flag_qndrop) + ENDIF @@ -3128,6 +3271,45 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw ENDDO ENDIF + IF ( ( ( present( ssat3d ) .and. present( f_ssat ) ) .or. & + ( present( ssati ) .and. present( f_ssati ) ) ) .and. makediag ) THEN + DO kz = kts,kte + DO ix = its,ite + + ! updated temperature and qv + temp1 = t0(ix,1,kz) ! an(ix,1,kz,lt)*t77(ix,1,kz) + ltemq = Int( (temp1-163.15)/fqsat+1.5 ) + ltemq = Min( nqsat, Max(1,ltemq) ) + + IF ( present( ssat3d ) .and. present( f_ssat ) ) THEN + IF ( f_ssat ) THEN + +! c1 = t00(ix,1,kz)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + c1 = (380.0/pn(ix,1,kz))*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + c1 = rdorv*esbolton*tabqvs(ltemq)/(pn(ix,1,kz) - esbolton*tabqvs(ltemq)) + ENDIF + + IF ( c1 > 0. ) THEN + ssat3d(ix,kz,jy) = 100.*(an(ix,1,kz,lv)/c1 - 1.0) ! from "new" values + ENDIF + + ENDIF + ENDIF + + IF ( present( ssati ) .and. present( f_ssati ) ) THEN + IF ( f_ssati ) THEN + t9s = (380.0/pn(ix,1,kz))*tabqis(ltemq) !saturation mixing ratio wrt ice + ssati(ix,kz,jy) = 100.*(an(ix,1,kz,lv)/t9s - 1.0) ! Min(t8s,max(an(ix,1,kz,lv),0.0))/t9s ! qv/qvi + ENDIF + ENDIF + + ENDDO + ENDDO + ENDIF + + ! compute diagnostic S-band reflectivity if needed IF ( present( dbz ) .and. makediag .and. lastlooptmp ) THEN @@ -3391,17 +3573,15 @@ REAL FUNCTION GAMMA_SP(xx) real xx integer j -! Double precision ser,stp,tmp,x,y,cof(6) - - real*8 ser,stp,tmp,x,y,cof(6) + double precision :: ser,stp,tmp,x,y,cof(6) SAVE cof,stp - DATA cof,stp/76.18009172947146d+0, & + DATA cof /76.18009172947146d+0, & & -86.50532032941677d0, & & 24.01409824083091d0, & & -1.231739572450155d0, & & 0.1208650973866179d-2,& - & -0.5395239384953d-5, & - & 2.5066282746310005d0/ + & -0.5395239384953d-5/ + DATA stp/2.5066282746310005d0/ IF ( xx <= 0.0 ) THEN write(0,*) 'Argument to gamma must be > 0!! xx = ',xx @@ -3774,20 +3954,18 @@ END function BETA DOUBLE PRECISION FUNCTION GAMMA_DP(xx) implicit none - double precision xx + double precision :: xx integer j -! Double precision ser,stp,tmp,x,y,cof(6) - - real*8 ser,stp,tmp,x,y,cof(6) + double precision ser,stp,tmp,x,y,cof(6) SAVE cof,stp - DATA cof,stp/76.18009172947146d+0, & + DATA cof /76.18009172947146d+0, & & -86.50532032941677d0, & & 24.01409824083091d0, & & -1.231739572450155d0, & & 0.1208650973866179d-2,& - & -0.5395239384953d-5, & - & 2.5066282746310005d0/ + & -0.5395239384953d-5/ + DATA stp/2.5066282746310005d0/ x = xx y = x @@ -4063,7 +4241,7 @@ subroutine hailmaxd(dtp,nx,ny,nz,an,na,nor,norz,alpha2d,dn, & integer :: plo, phi integer :: ialp, i, j - logical :: debug_mpi = .TRUE. + logical :: debug_mpi = .false. ! ################################################################### @@ -4268,7 +4446,6 @@ subroutine sediment1d(dtp,nx,ny,nz,an,na,nor,norz,xfall,dn,dz3d,dz3dinv, & ! real gz(-nor+ng1:nz+nor),z1d(-nor+ng1:nz+nor,4) real dtp real xfall(nx,ny,na) ! array for stuff landing on the ground -! real xfall0(nx,ny) ! dummy array integer infdo integer jslab ! which line of xfall to use @@ -4276,14 +4453,6 @@ subroutine sediment1d(dtp,nx,ny,nz,an,na,nor,norz,xfall,dn,dz3d,dz3dinv, & real tmp, vtmax, dtptmp, dtfrac real, parameter :: dz = 200. -! real :: xvt(nz+1,nx,3,lc:lhab) ! (nx,nz,2,lc:lhab) ! 1=mass-weighted, 2=number-weighted -! real :: tmpn(-nor+ng1:nx+nor,-nor+ng1:ny+nor,-norz+ng1:nz+norz) -! real :: tmpn2(-nor+ng1:nx+nor,-nor+ng1:ny+nor,-norz+ng1:nz+norz) -! real :: z(-nor+ng1:nx+nor,-norz+ng1:nz+norz,lr:lhab) -! real :: db1(nx,nz+1),dtz1(nz+1,nx,0:1),dz2dinv(nz+1,nx),db1inv(nx,nz+1) - -! real :: rhovtzx(nz,nx) - real, allocatable :: db1(:,:), dtz1(:,:,:),dz2dinv(:,:),db1inv(:,:) ! db1(nx,nz+1),dtz1(nz+1,nx,0:1),dz2dinv(nz+1,nx),db1inv(nx,nz+1) real, allocatable :: rhovtzx(:,:) real, allocatable :: xfall0(:,:), xvt(:,:,:,:),tmpn(:,:,:),tmpn2(:,:,:),z(:,:,:) @@ -4294,33 +4463,6 @@ subroutine sediment1d(dtp,nx,ny,nz,an,na,nor,norz,xfall,dn,dz3d,dz3dinv, & integer :: ngs ! = 512 integer :: ngscnt,mgs,ipconc0 -! real :: qx(ngs,lv:lhab) -! real :: qxw(ngs,ls:lhab) -! real :: cx(ngs,lc:lhab) -! real :: xv(ngs,lc:lhab) -! real :: vtxbar(ngs,lc:lhab,3) -! real :: xmas(ngs,lc:lhab) -! real :: xdn(ngs,lc:lhab) -! real :: xdia(ngs,lc:lhab,3) -! real :: vx(ngs,li:lhab) -! real :: alpha(ngs,lc:lhab) -! real :: zx(ngs,lr:lhab) -! logical :: hasmass(nx,lc+1:lhab) -! -! integer igs(ngs),kgs(ngs) -! -! real rho0(ngs),temcg(ngs) -! -! real temg(ngs) -! -! real rhovt(ngs) -! -! real cwnc(ngs),cinc(ngs) -! real fadvisc(ngs),cwdia(ngs),cipmas(ngs) -! -! real cimasn,cimasx,cnina(ngs),cimas(ngs) -! -! real cnostmp(ngs) real, allocatable :: qx(:,:) real, allocatable :: qxw(:,:) @@ -4359,8 +4501,6 @@ subroutine sediment1d(dtp,nx,ny,nz,an,na,nor,norz,xfall,dn,dz3d,dz3dinv, & integer :: ixe, jye, kze integer :: plo, phi - logical :: debug_mpi = .TRUE. - ! ################################################################### @@ -4558,12 +4698,9 @@ subroutine sediment1d(dtp,nx,ny,nz,an,na,nor,norz,xfall,dn,dz3d,dz3dinv, & ENDIF ! (n .ge. 2) - IF ( il >= lr .and. ( infall .eq. 3 .or. infall .eq. 4 ) .and. ln(il) > 0 ) THEN - IF ( (il .eq. lr .and. irfall .eq. infall .and. lzr < 1) .or. & - (il .ge. lh .and. lz(il) .lt. 1 ) .or. (il == ls .and. isfall == infall ) ) THEN + IF ( il >= lr .and. ( linfall(il) .eq. 3 .or. linfall(il) .eq. 4 ) .and. ln(il) > 0 ) THEN call calczgr1d(nx,ny,nz,nor,na,an,ixe,kze, & & z,db1,jgs,ipconc, dnu(il), il, ln(il), qxmin(il), xvmn(il), xvmx(il), lvol(il), xdn0(il), ix ) - ENDIF ENDIF if (ndebug .gt. 0 ) write(0,*) 'dbg = 1b' @@ -4600,34 +4737,25 @@ subroutine sediment1d(dtp,nx,ny,nz,an,na,nor,norz,xfall,dn,dz3d,dz3dinv, & IF ( ipconc .gt. 0 ) THEN !{ IF ( ipconc .ge. ipc(il) ) THEN - IF ( ( infall .ge. 2 .or. (infall .eq. 0 .and. il .lt. lh) ) .and. lz(il) .lt. 1) THEN !{ + IF ( ( linfall(il) .ge. 2 ) .and. lz(il) .lt. 1) THEN !{ ! ! load number conc. into tmpn to do fallout by mass-weighted mean fall speed ! to put a lower bound on number conc. ! - IF ( ( infall .eq. 3 .or. infall .eq. 4 ) .and. ( (il == ls .and. isfall .eq. infall ) & - & .or. il .eq. lh .or. il .eq. lhl .or. il == lf .or. & - & ( il .eq. lr .and. irfall .eq. infall) ) ) THEN - - ! set up for method I+II + IF ( linfall(il) == 3 .or. linfall(il) == 4 ) THEN + ! set up for method I or I+II DO kz = kzb,kze -! DO ix = ixb,ixe tmpn2(ix,jy,kz) = z(ix,kz,il) -! ENDDO ENDDO DO kz = kzb,kze -! DO ix = ixb,ixe tmpn(ix,jy,kz) = an(ix,jy,kz,ln(il)) -! ENDDO ENDDO ELSE ! set up for method II only DO kz = kzb,kze -! DO ix = ixb,ixe tmpn(ix,jy,kz) = an(ix,jy,kz,ln(il)) -! ENDDO ENDDO ENDIF @@ -4638,17 +4766,14 @@ subroutine sediment1d(dtp,nx,ny,nz,an,na,nor,norz,xfall,dn,dz3d,dz3dinv, & if (ndebug .gt. 0 ) write(0,*) 'dbg = 3f' in = 2 - IF ( infall .eq. 1 ) in = 1 + IF ( linfall(il) .eq. 1 ) in = 1 call fallout1d(nx,ny,nz,nor,na,dtptmp,dtfrac,jgs,xvt(1,1,in,il), & & an,db1,ln(il),0,xfall,dtz1,ix) - IF ( lz(il) .lt. 1 ) THEN ! if not 3-moment, run one of the correction schemes - IF ( (infall .ge. 2 .or. infall .eq. 3) .and. .not. (infall .eq. 0 .and. il .ge. lh) & - & .and. ( il .eq. lr .or. (il .ge. li .and. il .le. lhab) )) THEN -! : .or. il .eq. lhl )) THEN - + IF ( lz(il) .lt. 1 ) THEN ! { if not 3-moment, run one of the correction schemes + IF ( linfall(il) >= 2 ) THEN xfall0(:,jgs) = 0.0 IF ( ( infall .eq. 3 .or. infall .eq. 4 ) .and. & @@ -4663,42 +4788,37 @@ subroutine sediment1d(dtp,nx,ny,nz,an,na,nor,norz,xfall,dn,dz3d,dz3dinv, & & tmpn,db1,1,0,xfall0,dtz1,ix) ENDIF - IF ( ( infall .eq. 3 .or. infall .eq. 4 ) .and. ( (il .eq. lr .and. irfall .eq. infall) & - & .or. il .ge. lh .or. (il == ls .and. isfall .eq. infall ) ) ) THEN -! "Method I" - dbz correction - + IF ( linfall(il) == 3 .or. linfall(il) == 4 ) THEN + ! "Method I" - dbz correction + ! Uses input tmpn2 (temp. Z-moment) to determine if new N and q values in an(:,:,:,ln(il)) + ! cause an increase in reflectivity moment. If so, either use N from mass-wgt Vt (tmpn) to replace + ! new N (infall=3; I) or use smaller N from tmpn or calculated from q and temporary Z (infall=4; I+II) + ! Uses 'z' array to check if new reflectivity is greater than pre-sedimentation reflectivity call calcnfromz1d(nx,ny,nz,nor,na,an,tmpn2,ixe,kze, & & z,db1,jgs,ipconc, dnu(il), il, ln(il), qxmin(il), xvmn(il), xvmx(il),tmpn, & & lvol(il), xdn0(il), infall, ix) - ELSEIF ( infall .eq. 5 .and. il .ge. lh .or. ( il == lr .and. irfall == 5 ) ) THEN + ELSEIF ( linfall(il) .eq. 5 .and. il .ge. lh .or. ( il == lr .and. irfall == 5 ) ) THEN DO kz = kzb,kze -! DO ix = ixb,ixe an(ix,jgs,kz,ln(il)) = Max( an(ix,jgs,kz,ln(il)), 0.5* ( an(ix,jgs,kz,ln(il)) + tmpn(ix,jy,kz) )) - -! ENDDO ENDDO ELSEIF ( .not. (il .eq. lr .and. irfall .eq. 0) .and. .not. (il .eq. ls .and. isfall .eq. 0) ) THEN ! "Method II" M-wgt N-fallout correction DO kz = kzb,kze -! DO ix = ixb,ixe - an(ix,jgs,kz,ln(il)) = Max( an(ix,jgs,kz,ln(il)), tmpn(ix,jy,kz) ) - -! ENDDO ENDDO - ENDIF - ENDIF ! lz(il) .lt. 1 + ENDIF !} + ENDIF - ENDIF - ENDIF + ENDIF !} lz(il) .lt. 1 + ENDIF ! ipconc > ipc - ENDIF !} + ENDIF !} (ipconc > 0) ENDDO ! n=1,ndfall @@ -4784,8 +4904,6 @@ subroutine fallout1d(nx,ny,nz,nor,na,dtp,dtfrac,jgs,vt, & integer :: ixb, jyb, kzb integer :: ixe, jye, kze - logical :: debug_mpi = .TRUE. - ! ################################################################### jy = 1 @@ -5205,7 +5323,7 @@ subroutine calcnfromq(nx,ny,nz,an,na,nor,norz,dn, & integer ix,jy,kz double precision vr,q,nrx,nrx2,rd,g1h,g1hl,g1r,g1s,zx,z,znew,zt,zxt,n1,laminv1 double precision :: zr, zs, zh, dninv - real, parameter :: xn0s = 3.0e6, xn0r = 8.0e6, xn0h = 2.0e5, xn0hl = 4.0e4 + real, parameter :: xn0s = 3.0e8, xn0r = 8.0e6, xn0h = 2.0e5, xn0hl = 4.0e4 real, parameter :: xdnr = 1000., xdns = 100. ,xdnh = 700.0, xdnhl = 900.0 real, parameter :: zhlfac = 1./(pi*xdnhl*xn0hl) real, parameter :: zhfac = 1./(pi*xdnh*xn0h) @@ -6083,7 +6201,7 @@ SUBROUTINE QVEXCESS(ngs,mgs,qwvp0,qv0,qcw1,pres,thetap0,theta0, & integer itertd integer ltemq - real gamss + real gamss, tmp real theta(ngs), qvap(ngs), pqs(ngs), qcw(ngs), qwv(ngs) real qcwtmp(ngs), qss(ngs), qvs(ngs), qwvp(ngs) real dqcw(ngs), dqwv(ngs), dqvcnd(ngs) @@ -6124,7 +6242,13 @@ SUBROUTINE QVEXCESS(ngs,mgs,qwvp0,qv0,qcw1,pres,thetap0,theta0, & ltemq = (temg(mgs)-163.15)/fqsat+1.5 ltemq = Min( nqsat, Max(1,ltemq) ) - qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + pqs(mgs) = (380.0)/(pres(mgs)) + qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + qvs(mgs) = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF + qss(mgs) = (0.01*ss1 + 1.0)*qvs(mgs) ! ! iterate adjustment @@ -6184,7 +6308,11 @@ SUBROUTINE QVEXCESS(ngs,mgs,qwvp0,qv0,qcw1,pres,thetap0,theta0, & ! tqvcon = temg(mgs)-cbw ltemq = (temg(mgs)-163.15)/fqsat+1.5 ltemq = Min( nqsat, Max(1,ltemq) ) - qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + qvs(mgs) = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF qcw(mgs) = max( 0.0, qcw(mgs) ) qwv(mgs) = max( 0.0, qvap(mgs)) qss(mgs) = (0.01*ss1 + 1.0)*qvs(mgs) @@ -6210,9 +6338,9 @@ END SUBROUTINE QVEXCESS ! SUBROUTINE setvtz(ngscnt,qx,qxmin,qxw,cx,rho0,rhovt,xdia,cno,cnostmp, & & xmas,vtxbar,xdn,xvmn0,xvmx0,xv,cdx,cdxgs, & - & ipconc1,ndebug1,ngs,nz,kgs,fadvisc, & + & ipconc1,ndebug1,ngs,nz,igs,kgs,fadvisc, & & cwmasn,cwmasx,cwradn,cnina,cimna,cimxa, & - & itype1a,itype2a,temcg,infdo,alpha,ildo,axx,bxx) + & itype1a,itype2a,temcg,infdo,alpha,axx,bxx,ildo) ! & itype1a,itype2a,temcg,infdo,alpha,ildo,axh,bxh,axhl,bxhl) @@ -6241,7 +6369,7 @@ SUBROUTINE setvtz(ngscnt,qx,qxmin,qxw,cx,rho0,rhovt,xdia,cno,cnostmp, & real cwc1, cimna, cimxa real cnina(ngs) - integer kgs(ngs) + integer igs(ngs),kgs(ngs) real fadvisc(ngs) real fsw @@ -7475,8 +7603,8 @@ SUBROUTINE setvtz(ngscnt,qx,qxmin,qxw,cx,rho0,rhovt,xdia,cno,cnostmp, & & (aax*(xdia(mgs,il,1) )**bbx * & & x)/y ! & Gamma(7.0 + alpha(mgs,il) + bbx)/Gamma(7. + alpha(mgs,il)) - IF ( .not. (vtxbar(mgs,il,1) > -1. .and. vtxbar(mgs,il,1) < 200. ) .or. & - .not. (vtxbar(mgs,il,3) > -1. .and. vtxbar(mgs,il,3) < 200. ) ) THEN + IF ( .not. (vtxbar(mgs,il,1) > -1. .and. vtxbar(mgs,il,1) < 250. ) .or. & + .not. (vtxbar(mgs,il,3) > -1. .and. vtxbar(mgs,il,3) < 250. ) ) THEN write(0,*) 'Setvtz: problem with vtxbar1/3: ',il,vtxbar(mgs,il,1),vtxbar(mgs,il,3),aax,bbx,x,y write(0,*) 'q, number, diam1,3(mm) = ', qx(mgs,il),cx(mgs,il),1000.*xdia(mgs,il,1),1000.*xdia(mgs,il,3) ! call commasmpi_abort() @@ -7712,7 +7840,7 @@ subroutine ziegfall1d(nx,ny,nz,nor,norz,na,dtp,jgs,ixcol, & logical ldoliq - real chw, qr, z, rd, alp, z1, g1, vr, nrx, tmp + real chw, qr, z, rd, alp, z1, g1, vr, nrx, tmp, tmpc, tmpz real vtmax real xvbarmax @@ -7735,7 +7863,7 @@ subroutine ziegfall1d(nx,ny,nz,nor,norz,na,dtp,jgs,ixcol, & logical :: debug_mpi = .false. - if (ndebugzf .gt. 0 ) write(0,*) "ZIEGFALL: ENTERED SUBROUTINE" + if (ndebugzf .gt. 0 ) write(0,*) "ZIEGFALL1D: ENTERED SUBROUTINE" ! ##################################################################### ! BEGIN EXECUTABLE @@ -7804,13 +7932,12 @@ subroutine ziegfall1d(nx,ny,nz,nor,norz,na,dtp,jgs,ixcol, & ngscnt = 0 - do kz = nzmpb,nz + do kz = 1,nz do ix = ixcol,ixcol flag = .false. - DO il = l1,l2 - flag = flag .or. ( an(ix,jy,kz,il) .gt. qxmin(il) ) + flag = flag .or. ( an(ix,jy,kz,il) > 0.0 ) ENDDO if ( flag ) then @@ -7819,7 +7946,7 @@ subroutine ziegfall1d(nx,ny,nz,nor,norz,na,dtp,jgs,ixcol, & ngscnt = ngscnt + 1 igs(ngscnt) = ix kgs(ngscnt) = kz - if ( ngscnt .eq. ngs ) goto 1100 + if ( ngscnt .eq. nz ) goto 1100 end if end do !!ix nxmpb = 1 @@ -7845,11 +7972,11 @@ subroutine ziegfall1d(nx,ny,nz,nor,norz,na,dtp,jgs,ixcol, & temg(mgs) = t0(igs(mgs),jy,kgs(mgs)) temcg(mgs) = temg(mgs) - tfr - + ! end do ! -! only need fadvisc for +! only need fadvisc for droplets IF ( lc .gt. 1 .and. (ildo == 0 .or. ildo == lc ) ) then do mgs = 1,ngscnt fadvisc(mgs) = advisc0*(416.16/(temg(mgs)+120.0))* & @@ -7898,58 +8025,52 @@ subroutine ziegfall1d(nx,ny,nz,nor,norz,na,dtp,jgs,ixcol, & cx(mgs,li) = Max(an(igs(mgs),jy,kgs(mgs),lni), 0.0) end do end if + if ( ipconc .ge. 2 .and. lc .gt. 1 .and. (ildo == 0 .or. ildo == lc ) ) then do mgs = 1,ngscnt cx(mgs,lc) = Max(an(igs(mgs),jy,kgs(mgs),lnc), 0.0) -! cx(mgs,lc) = Min( ccwmx, cx(mgs,lc) ) end do end if + if ( ipconc .ge. 3 .and. lr .gt. 1 .and. (ildo == 0 .or. ildo == lr ) ) then do mgs = 1,ngscnt cx(mgs,lr) = Max(an(igs(mgs),jy,kgs(mgs),lnr), 0.0) -! IF ( qx(mgs,lr) .le. qxmin(lr) ) THEN -! ELSE -! cx(mgs,lr) = Max( 0.0, cx(mgs,lr) ) -! ENDIF end do end if + if ( ipconc .ge. 4 .and. ls .gt. 1 .and. (ildo == 0 .or. ildo == ls ) ) then do mgs = 1,ngscnt cx(mgs,ls) = Max(an(igs(mgs),jy,kgs(mgs),lns), 0.0) -! IF ( qx(mgs,ls) .le. qxmin(ls) ) THEN -! ELSE -! cx(mgs,ls) = Max( 0.0, cx(mgs,ls) ) -! ENDIF end do end if if ( ipconc .ge. 5 .and. lh .gt. 1 .and. (ildo == 0 .or. ildo == lh ) ) then do mgs = 1,ngscnt - cx(mgs,lh) = Max(an(igs(mgs),jy,kgs(mgs),lnh), 0.0) -! IF ( qx(mgs,lh) .le. qxmin(lh) ) THEN -! ELSE -! cx(mgs,lh) = Max( 0.0, cx(mgs,lh) ) -! ENDIF - end do ENDIF if ( ipconc .ge. 5 .and. lhl .gt. 1 .and. (ildo == 0 .or. ildo == lhl ) ) then do mgs = 1,ngscnt - cx(mgs,lhl) = Max(an(igs(mgs),jy,kgs(mgs),lnhl), 0.0) -! IF ( qx(mgs,lhl) .le. qxmin(lhl) ) THEN -! cx(mgs,lhl) = 0.0 -! ELSEIF ( cx(mgs,lhl) .eq. 0.0 .and. qx(mgs,lhl) .lt. 3.0*qxmin(lhl) ) THEN -! qx(mgs,lhl) = 0.0 -! ELSE -! cx(mgs,lhl) = Max( 0.0, cx(mgs,lhl) ) -! ENDIF - end do end if - + + ! Vaporize tiny values + DO il = l1,l2 + IF ( lz(il) < 1 ) THEN + do mgs = 1,ngscnt + IF ( cx(mgs,il) <= cxmin .or. qx(mgs,il) < qxmin(il) ) THEN + cx(mgs,il) = 0.0 + an(igs(mgs),jgs,kgs(mgs),lv) = an(igs(mgs),jgs,kgs(mgs),lv) + an(igs(mgs),jgs,kgs(mgs),il) + qx(mgs,il) = 0.0 + an(igs(mgs),jgs,kgs(mgs),il) = qx(mgs,il) + an(igs(mgs),jgs,kgs(mgs),ln(il)) = cx(mgs,il) + ENDIF + end do + ENDIF + ENDDO + do mgs = 1,ngscnt xdn(mgs,lc) = xdn0(lc) xdn(mgs,lr) = xdn0(lr) @@ -8377,7 +8498,7 @@ subroutine ziegfall1d(nx,ny,nz,nor,norz,na,dtp,jgs,ixcol, & ! check for artificial breakup (graupel/hail larger than allowed max size) - IF ( imaxdiaopt == 1 ) THEN + IF ( imaxdiaopt == 1 .or. il /= lr ) THEN xvbarmax = xvmx(il) ELSEIF ( imaxdiaopt == 2 ) THEN ! test against maximum mass diameter xvbarmax = xvmx(il) /((3. + alpha(mgs,il))**3/((3. + alpha(mgs,il))*(2. + alpha(mgs,il))*(1. + alpha(mgs,il)))) @@ -8393,7 +8514,16 @@ subroutine ziegfall1d(nx,ny,nz,nor,norz,na,dtp,jgs,ixcol, & IF ( tmp < cx(mgs,il) ) THEN ! breakup g1 = 36.*(6.0 + alpha(mgs,il))*(5.0 + alpha(mgs,il))*(4.0 + alpha(mgs,il))/ & & ((3.0 + alpha(mgs,il))*(2.0 + alpha(mgs,il))*(1.0 + alpha(mgs,il))*pi**2) - zx(mgs,il) = zx(mgs,il) + g1*(rho0(mgs)/xdn(mgs,il))**2*( (qx(mgs,il)/tmp)**2 * (tmp-cx(mgs,il)) ) + ! check if incoming zx is consistent + ! Z from incoming cx, qx, and alpha + tmpz = g1/(pi/6.*xdn(mgs,il))**2 * ((rho0(mgs)*qx(mgs,il))**2)/tmp + IF ( tmpz > zx(mgs,il) ) THEN + ! find cx that gives zx + tmpc = g1/(pi/6.*xdn(mgs,il))**2 * ((rho0(mgs)*qx(mgs,il))**2)/zx(mgs,il) + cx(mgs,il) = Max(cx(mgs,il), tmpc) + ENDIF + zx(mgs,il) = g1/(pi/6.*xdn(mgs,il))**2 * ((rho0(mgs)*qx(mgs,il))**2)/cx(mgs,il) +! zx(mgs,il) = zx(mgs,il) + g1*(rho0(mgs)/xdn(mgs,il))**2*( (qx(mgs,il)/tmp)**2 * (tmp-cx(mgs,il)) ) an(igs(mgs),jgs,kgs(mgs),lz(il)) = zx(mgs,il) chw = cx(mgs,il) @@ -8467,9 +8597,9 @@ subroutine ziegfall1d(nx,ny,nz,nor,norz,na,dtp,jgs,ixcol, & call setvtz(ngscnt,qx,qxmin,qxw,cx,rho0,rhovt,xdia,cno,cnostmp, & & xmas,vtxbar,xdn,xvmn,xvmx,xv,cdx,cdxgs, & - & ipconc,ndebugzf,ngs,nz,kgs,fadvisc, & + & ipconc,ndebugzf,ngs,nz,igs,kgs,fadvisc, & & cwmasn,cwmasx,cwradn,cnina,cimn,cimx, & - & itype1,itype2,temcg,infdo,alpha,ildo,axx,bxx) + & itype1,itype2,temcg,infdo,alpha,axx,bxx,ildo) ! & itype1,itype2,temcg,infdo,alpha,ildo,axh,bxh,axhl,bxhl) @@ -8673,12 +8803,12 @@ subroutine radardd02(nx,ny,nz,nor,na,an,temk, & real dtmp (nx,nz) real tmp - real*8 dtmps, dtmpr, dtmph, dtmphl, g1, zx, ze, x + double precision :: dtmps, dtmpr, dtmph, dtmphl, g1, zx, ze, x integer i,j,k,ix,jy,kz,ihcnt - real*8 xcnoh, xcnos, dadh, dads, zhdryc, zsdryc, zhwetc,zswetc - real*8 dadr + double precision :: xcnoh, xcnos, dadh, dads, zhdryc, zsdryc, zhwetc,zswetc + double precision :: dadr real dbzmax,dbzmin parameter ( dbzmin = 0 ) @@ -9700,7 +9830,7 @@ SUBROUTINE NUCOND & ! - real ccnc(ngs), ccna(ngs), cnuc(ngs), cwnccn(ngs) + real ccnc(ngs), ccna(ngs), cnuc(ngs), cwnccn(ngs), ccnaco(ngs), ccnanu(ngs) real :: ccnc_nu(ngs), ccnc_ac(ngs), ccnc_co(ngs) real ccncuf(ngs) real sscb ! 'cloud base' SS threshold @@ -9724,7 +9854,7 @@ SUBROUTINE NUCOND & real volb, t2s real, parameter :: aa1 = 9.44e15, aa2 = 5.78e3 ! a1 in Ziegler - real ec0, ex1, ft, rhoinv(ngs) + real rhoinv(ngs) real chw, g1, rd1 @@ -9750,6 +9880,7 @@ SUBROUTINE NUCOND & real dcrit real cn(ngs), cnuf(ngs) real :: ccwmax + integer ltemq @@ -9828,16 +9959,35 @@ SUBROUTINE NUCOND & integer, parameter :: iunit = 0 - real :: frac, hwdn, tmpg + real :: frac, hwdn, tmpg, xdia1, xdia3, cwch,xvol real :: cvm,cpm,rmm real, parameter :: cpv = 1885.0 ! specific heat of water vapor at constant pressure - + real, parameter :: Mair = 0.0284 ! MOLECULAR WEIGHT OF 'AIR' (KG/MOL) + + integer :: kstag integer :: count +! Addtion T.Iguchi Y2021 Update + real, parameter :: mwwater = 0.01801528 ! Molecular weight of water (kg/mol) + real, parameter :: rhowater = 997.0 ! Density of liquid water (kg/m3) + real, parameter :: gasconst = 8.3144598 ! Gas constant (m2 kg s-2 K-1 mol-1) + real :: sswater ! unit change supersaturation from percentage to n/a + real :: sigvl, aact + + real :: alpha_ar, gamma_ar, G_ar, evs, zeta, smax + real :: f_ac, g_ac, eta_ac + real :: f_nu, g_nu, eta_nu + real :: f_co, g_co, eta_co + + real :: sm_nu, sm_ac, sm_co, ss_ac, ss_nu, ss_co + real :: uu_nu, uu_ac, uu_co + + real :: cn_ac, cn_co, cn_nu + ! ------------------------------------------------------------------------------- itile = nxi jtile = ny @@ -9877,7 +10027,12 @@ SUBROUTINE NUCOND & ltemq = Int( (temp1-163.15)/fqsat+1.5 ) ltemq = Min( nqsat, Max(1,ltemq) ) - c1 = t00(ix,jy,kz)*tabqvs(ltemq) +! c1 = t00(ix,jy,kz)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + c1 = t00(ix,jy,kz)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + c1 = rdorv*esbolton*tabqvs(ltemq)/(pn(ix,jy,kz) + pb(kz) - esbolton*tabqvs(ltemq)) + ENDIF IF ( c1 > 0. ) THEN ssfilt(ix,jy,kz) = 100.*(an(ix,jy,kz,lv)/c1 - 1.0) ! from "new" values @@ -9919,6 +10074,7 @@ SUBROUTINE NUCOND & do kz = kzb,kze do ix = nxmpb,nxi + pres(1) = pn(ix,jy,kz) + pb(kz) pqs(1) = 380.0/(pn(ix,jy,kz) + pb(kz)) theta(1) = an(ix,jy,kz,lt) temg(1) = t0(ix,jy,kz) @@ -9926,7 +10082,12 @@ SUBROUTINE NUCOND & temcg(1) = temg(1) - tfr ltemq = (temg(1)-163.15)/fqsat+1.5 ltemq = Min( nqsat, Max(1,ltemq) ) - qvs(1) = pqs(1)*tabqvs(ltemq) + ! qvs(1) = pqs(1)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + qvs(1) = pqs(1)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + qvs(1) = rdorv*esbolton*tabqvs(ltemq)/(pres(1) - esbolton*tabqvs(ltemq)) + ENDIF qis(1) = pqs(1)*tabqis(ltemq) qss(1) = qvs(1) @@ -10005,11 +10166,21 @@ SUBROUTINE NUCOND & pqs(mgs) = (380.0)/(pres(mgs)) ltemq = (temg(mgs)-163.15)/fqsat+1.5 ltemq = Min( nqsat, Max(1,ltemq) ) - qvs(mgs) = pqs(mgs)*tabqvs(ltemq) +! qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + qvs(mgs) = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF qis(mgs) = pqs(mgs)*tabqis(ltemq) ! qvap(mgs) = max( (qwvp(mgs) + qv0(mgs)), 0.0 ) - es(mgs) = 6.1078e2*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + es(mgs) = 6.1078e2*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + es(mgs) = esbolton*tabqvs(ltemq) + ENDIF +! es(mgs) = 6.1078e2*tabqvs(ltemq) qss(mgs) = qvs(mgs) @@ -10082,7 +10253,37 @@ SUBROUTINE NUCOND & ccnc(mgs) = an(igs(mgs),jy,kgs(mgs),lccn) + an(igs(mgs),jy,kgs(mgs),lccnuf) ELSE ccnc(mgs) = an(igs(mgs),jy,kgs(mgs),lccn) + IF ( lccna > 1 ) THEN + cnuc(mgs) = ccnc(mgs) + ENDIF + ENDIF + IF ( lcn_nu > 1 ) THEN + ccnc_nu(mgs) = an(igs(mgs),jy,kgs(mgs),lcn_nu) + ENDIF + IF ( lcn_co > 1 ) THEN + ccnc_co(mgs) = an(igs(mgs),jy,kgs(mgs),lcn_co) + ENDIF + IF ( lccnaco > 1 ) THEN + ccnaco(mgs) = an(igs(mgs),jy,kgs(mgs),lccnaco) + ELSE + ccnaco(mgs) = 0.0 + ENDIF + IF ( lccnanu > 1 ) THEN + ccnanu(mgs) = an(igs(mgs),jy,kgs(mgs),lccnanu) + ELSE + ccnanu(mgs) = 0.0 ENDIF + ELSEIF ( lccn > 1 .and. ( ac_opt == 1 .or. ac_opt == 11 ) ) THEN + ccnc(mgs) = an(igs(mgs),jy,kgs(mgs),lccn) + ! ccnc(mgs) = ccnc_ac(mgs) + cnuc(mgs) = ccnc(mgs) + cwnccn(mgs) = cnuc(mgs) + ! write(0,*) 'ccnc_ac,mgs = ', ccnc_ac(mgs),mgs,igs(mgs),jy,kgs(mgs) + ELSEIF ( lccn > 1 .and. ( ac_opt == 2 .or. ac_opt == 22 ) ) THEN + ccnc_nu(mgs) = an(igs(mgs),jy,kgs(mgs),lcn_nu) + ccnc(mgs) = an(igs(mgs),jy,kgs(mgs),lccn) + ! ccnc(mgs) = ccnc_ac(mgs) + ccnc_co(mgs) = an(igs(mgs),jy,kgs(mgs),lcn_co) ELSE ccnc(mgs) = cwnccn(mgs) ENDIF @@ -10094,9 +10295,21 @@ SUBROUTINE NUCOND & cnuf(mgs) = 0.0 IF ( lccna > 1 ) THEN ccna(mgs) = an(igs(mgs),jy,kgs(mgs),lccna) ! predicted count of activated ccn + IF ( ac_opt == 22 ) THEN + IF ( lccnaco > 1 ) THEN + ccnaco(mgs) = an(igs(mgs),jy,kgs(mgs),lccnaco) + ELSE + ccnaco(mgs) = 0.0 + ENDIF + IF ( lccnanu > 1 ) THEN + ccnanu(mgs) = an(igs(mgs),jy,kgs(mgs),lccnanu) + ELSE + ccnanu(mgs) = 0.0 + ENDIF + ENDIF ELSE IF ( lccn > 1 ) THEN - ccna(mgs) = cwnccn(mgs) - ccnc(mgs) ! diagnose activated ccn as background value - remaining unactivated ccn + ccna(mgs) = 0.0 ! WRF driver interface already has ccw subtracted from ccnc ELSE ccna(mgs) = cx(mgs,lc) ! approximation of number of activated ccn ENDIF @@ -10113,9 +10326,9 @@ SUBROUTINE NUCOND & DO mgs = 1,ngscnt ! default value of renucfrac is 0.0 IF ( irenuc /= 6 ) THEN - cnuc(mgs) = Max(ccnc(mgs),cwnccn(mgs))*(1. - renucfrac) + ccnc(mgs)*renucfrac + cnuc(mgs) = ccnc(mgs)*(1. - renucfrac) + ccnc(mgs)*renucfrac ELSE - cnuc(mgs) = Max(ccnc(mgs),cwnccn(mgs))*(1. - renucfrac) + Max(0.0,ccnc(mgs) - ccna(mgs))*renucfrac + cnuc(mgs) = ccnc(mgs)*(1. - renucfrac) + Max(0.0,ccnc(mgs) - ccna(mgs))*renucfrac ENDIF IF ( renucfrac >= 0.999 ) THEN IF ( temg(mgs) < 265. ) THEN @@ -10526,16 +10739,27 @@ SUBROUTINE NUCOND & QEVAP= Min( qx(mgs,lc), R1*(qss(mgs)-qvap(mgs)) ) - IF ( qx(mgs,lc) <= QEVAP ) THEN ! GO TO 63 + IF ( qx(mgs,lc) <= QEVAP ) THEN !{ GO TO 63 qwvp(mgs) = qwvp(mgs) + qx(mgs,lc) thetap(mgs) = thetap(mgs) - felvcp(mgs)*qx(mgs,lc)/(pi0(mgs)) IF ( io_flag .and. nxtra > 1 ) THEN axtra(igs(mgs),jy,kgs(mgs),1) = -qx(mgs,lc)/dtp ENDIF qx(mgs,lc) = 0. - IF ( restoreccn ) THEN + IF ( restoreccn ) THEN !{ IF ( lccna > 1 ) THEN - ccna(mgs) = ccna(mgs) - restoreccnfrac*cx(mgs,lc) + tmp = restoreccnfrac*cx(mgs,lc) + IF ( lccnaco > 1 .and. lccnanu > 1 ) THEN + ! restore CCN proportionally to each type, although coarse are presumably already lost to rain + tmp2 = ccna(mgs) + ccnaco(mgs) + ccnanu(mgs) + IF ( tmp2 > 0.0 ) THEN + ccna(mgs) = ccna(mgs) - tmp*ccna(mgs)/tmp2 + ccnaco(mgs) = ccnaco(mgs) - tmp*ccnaco(mgs)/tmp2 + ccnanu(mgs) = ccnanu(mgs) - tmp*ccnanu(mgs)/tmp2 + ENDIF + ELSE + ccna(mgs) = ccna(mgs) - tmp + ENDIF ELSEIF ( irenuc <= 2 ) THEN IF ( .not. invertccn ) THEN ccnc(mgs) = Max( ccnc(mgs), Min( qccn*rho0(mgs), ccnc(mgs) + restoreccnfrac*cx(mgs,lc) ) ) @@ -10543,16 +10767,27 @@ SUBROUTINE NUCOND & ccnc(mgs) = ccnc(mgs) + restoreccnfrac*cx(mgs,lc) ENDIF ENDIF - ENDIF + ENDIF !} cx(mgs,lc) = 0. - ELSE + ELSE !} { qctmp = qx(mgs,lc) qwvp(mgs) = qwvp(mgs) + QEVAP qx(mgs,lc) = qx(mgs,lc) - QEVAP IF ( qx(mgs,lc) .le. 0. ) THEN IF ( restoreccn ) THEN IF ( lccna > 1 ) THEN - ccna(mgs) = ccna(mgs) - restoreccnfrac*cx(mgs,lc) + tmp = restoreccnfrac*cx(mgs,lc) + IF ( lccnaco > 1 .and. lccnanu > 1 ) THEN + ! restore CCN proportionally to each type, although coarse are presumably already lost to rain + tmp2 = ccna(mgs) + ccnaco(mgs) + ccnanu(mgs) + IF ( tmp2 > 0.0 ) THEN + ccna(mgs) = ccna(mgs) - tmp*ccna(mgs)/tmp2 + ccnaco(mgs) = ccnaco(mgs) - tmp*ccnaco(mgs)/tmp2 + ccnanu(mgs) = ccnanu(mgs) - tmp*ccnanu(mgs)/tmp2 + ENDIF + ELSE + ccna(mgs) = ccna(mgs) - tmp + ENDIF ELSEIF ( irenuc <= 2 ) THEN ! ccnc(mgs) = Max( ccnc(mgs), Min( qccn*rho0(mgs), ccnc(mgs) + cx(mgs,lc) ) ) ! ccnc(mgs) = ccnc(mgs) + cx(mgs,lc) @@ -10568,7 +10803,19 @@ SUBROUTINE NUCOND & tmp = 0.9*QEVAP*cx(mgs,lc)/qctmp ! let droplets get smaller but also remove some. A factor of 1.0 would maintain same size IF ( restoreccn ) THEN IF ( lccna > 1 ) THEN - ccna(mgs) = ccna(mgs) - restoreccnfrac*tmp + tmp = restoreccnfrac*tmp + IF ( lccnaco > 1 .and. lccnanu > 1 ) THEN + ! restore CCN proportionally to each type, although coarse are presumably already lost to rain + tmp2 = ccna(mgs) + ccnaco(mgs) + ccnanu(mgs) + IF ( tmp2 > 0.0 ) THEN + ccna(mgs) = ccna(mgs) - tmp*ccna(mgs)/tmp2 + ccnaco(mgs) = ccnaco(mgs) - tmp*ccnaco(mgs)/tmp2 + ccnanu(mgs) = ccnanu(mgs) - tmp*ccnanu(mgs)/tmp2 + ENDIF + ELSE + ccna(mgs) = ccna(mgs) - tmp + ENDIF + ! ccna(mgs) = ccna(mgs) - restoreccnfrac*tmp ELSEIF ( irenuc <= 2 ) THEN ! ccnc(mgs) = Max( ccnc(mgs), Min( qccn*rho0(mgs), ccnc(mgs) + tmp ) ) ! ccnc(mgs) = ccnc(mgs) + tmp @@ -10586,7 +10833,7 @@ SUBROUTINE NUCOND & axtra(igs(mgs),jy,kgs(mgs),1) = -QEVAP/dtp ENDIF - ENDIF + ENDIF !} GO TO 631 @@ -10710,7 +10957,11 @@ SUBROUTINE NUCOND & ltemq = Min( nqsat, Max(1,ltemq) ) ltemq1 = ltemq temp1 = temg(mgs) - p380 = 380.0/pres(mgs) + IF ( iqvsopt == 0 ) THEN + p380 = 380.0/pres(mgs) + ELSE + p380 = esbolton*rdorv/(pres(mgs) - es(mgs)) + ENDIF ! taus = Max( 0.05*dtp, Min(taus, 0.25*dtp ) ) ! nc = NInt(dtp/Min(1.0,0.5*taus)) @@ -10880,7 +11131,11 @@ SUBROUTINE NUCOND & temg(mgs) = theta(mgs)*f1 ltemq = (temg(mgs)-163.15)/fqsat+1.5 ltemq = Min( nqsat, Max(1,ltemq) ) - qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + qvs(mgs) = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF ! es(mgs) = 6.1078e2*tabqvs(ltemq) ! @@ -10918,7 +11173,7 @@ SUBROUTINE NUCOND & ! IF ( ssf(mgs) > ssmx .and. ssf(mgs) < 20.0 ) THEN ! test -- fails ! IF ( ssf(mgs) > ssmx .and. ssf(mgs) < 20.0 .and. ccnc(mgs) > 0.1*cwnccn(mgs)) THEN ! test -- is OK IF ( ssf(mgs) > ssmx .and. ssf(mgs) < 20.0 .and. & - ( ccnc(mgs) > 0.05*cwnccn(mgs) .or. ( ac_opt > 0 .and. ccnc_ac(mgs) - cx(mgs,lc) > 0.0 ) ) ) THEN ! test + ( ccnc(mgs) > 0.05*cwnccn(mgs) .or. ( ac_opt > 0 .and. ccnc(mgs) - cx(mgs,lc) > 0.0 ) ) ) THEN ! test ! IF ( ssf(mgs) > ssmx ) THEN ! original condition CALL QVEXCESS(ngs,mgs,qwvp,qv0,qx(1,lc),pres,thetap,theta0,dcloud, & & pi0,tabqvs,nqsat,fqsat,cbw,fcqv1,felvcp,ssmx,pk,ngscnt) @@ -10941,7 +11196,12 @@ SUBROUTINE NUCOND & ! temg(mgs) = theta2temp( theta(mgs), pres(mgs) ) ltemq = (temg(mgs)-163.15)/fqsat+1.5 ltemq = Min( nqsat, Max(1,ltemq) ) - qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + ! qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + qvs(mgs) = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF ! es(mgs) = 6.1078e2*tabqvs(ltemq) !.... S. TWOMEY (1959) @@ -10954,11 +11214,11 @@ SUBROUTINE NUCOND & IF ( .not. flag_qndrop ) THEN ! { do not calculate number of droplets if using wrf-chem - IF ( ac_opt == 0 ) THEN + ! IF ( ac_opt == 0 ) THEN cnuctmp = cnuc(mgs) - ELSE - cnuctmp = ccnc_ac(mgs) - ENDIF + ! ELSE + ! cnuctmp = ccnc(mgs) + ! ENDIF ! IF ( ssmax(mgs) .lt. sscb .and. qx(mgs,lc) .gt. qxmin(lc)) THEN IF ( dcloud .gt. qxmin(lc) .and. wvel(mgs) > 0.0) THEN @@ -10991,11 +11251,11 @@ SUBROUTINE NUCOND & ! ccnc(mgs) = 0.0 ENDIF ELSE - cn(mgs) = Min( cn(mgs), ccnc_ac(mgs) ) + cn(mgs) = Min( cn(mgs), ccnc(mgs) ) ENDIF ! cx(mgs,lc) = cx(mgs,lc) + cn(mgs) - IF ( irenuc <= 2 .and. lccna < 1 ) ccnc(mgs) = Max(0.0, ccnc(mgs) - cn(mgs)) - ccna(mgs) = ccna(mgs) + cn(mgs) + IF ( irenuc <= 2 .and. lccna < 1 ) ccnc(mgs) = Max(0.0, ccnc(mgs) - cn(mgs)) + ccna(mgs) = ccna(mgs) + cn(mgs) ENDIF ! write(91,*) 'nuc1: cn, ix, kz = ',cn(mgs),igs(mgs),kgs(mgs),wvel(mgs),cnexp,ccnc(mgs) @@ -11289,7 +11549,11 @@ SUBROUTINE NUCOND & ltemq = Int( (temp1-163.15)/fqsat+1.5 ) ltemq = Min( nqsat, Max(1,ltemq) ) - c1= pqs(mgs)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + c1 = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + c1 = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF IF ( c1 > 0. ) THEN ssf(mgs) = Max(0.0, 100.*((qv0(mgs) + qwvp(mgs))/c1 - 1.0) ) ! from "new" values ELSE @@ -11367,7 +11631,11 @@ SUBROUTINE NUCOND & ltemq = Min( nqsat, Max(1,ltemq) ) ! c1 = t00(igs(mgs),jy,kgs(mgs))*tabqvs(ltemq) - c1= pqs(mgs)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + c1 = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + c1 = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF ssf(mgs) = 0.0 IF ( c1 > 0. ) THEN @@ -11468,7 +11736,11 @@ SUBROUTINE NUCOND & ltemq = Min( nqsat, Max(1,ltemq) ) ! c1 = t00(igs(mgs),jy,kgs(mgs))*tabqvs(ltemq) - c1= pqs(mgs)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + c1 = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + c1 = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF ssf(mgs) = 0.0 IF ( c1 > 0. ) THEN @@ -11502,14 +11774,24 @@ SUBROUTINE NUCOND & ENDIF + ELSEIF ( irenuc == 9 .or. irenuc == 10 ) THEN ! } { + + write(0,*) 'irenuc=9 requires nuwrfmods=1' + + ELSEIF ( irenuc == 11 ) THEN ! } { + + write(0,*) 'irenuc=11 requires nuwrfmods=1' ENDIF ! } + ccna(mgs) = ccna(mgs) + cn(mgs) + + ENDIF ! irenuc >= 0 .and. .not. flag_qndrop - IF( cx(mgs,lc) .GT. 0. .AND. qx(mgs,lc) .LE. qxmin(lc)) cx(mgs,lc)=0. + ! IF( cx(mgs,lc) .GT. 0. .AND. qx(mgs,lc) .LE. qxmin(lc)) cx(mgs,lc)=0. GO TO 631 !.... NUCLEATION ON CLOUD INFLOW BOUNDARY POINT @@ -11542,39 +11824,39 @@ SUBROUTINE NUCOND & IF ( qvex .gt. 0.0 ) THEN - thetap(mgs) = thetap(mgs) + felvcp(mgs)*qvex/(pi0(mgs)) - IF ( io_flag .and. nxtra > 1 ) THEN - axtra(igs(mgs),jy,kgs(mgs),1) = axtra(igs(mgs),jy,kgs(mgs),1) + qvex/dtp - ENDIF - qwvp(mgs) = qwvp(mgs) - qvex - qx(mgs,lc) = qx(mgs,lc) + qvex - IF ( .not. flag_qndrop) THEN - IF ( imaxsupopt == 1 ) THEN - cn(mgs) = Min( Max(ccnc(mgs),cwnccn(mgs)), rho0(mgs)*qvex/Max( cwmasn5, xmas(mgs,lc) ) ) - ELSEIF ( imaxsupopt == 2 ) THEN - cn(mgs) = Min( Max(ccnc(mgs),cwnccn(mgs)), rho0(mgs)*qvex/Max( cwmasn5, Max(cwmas30,xmas(mgs,lc)) ) ) - ELSEIF ( imaxsupopt == 3 ) THEN - cn(mgs) = Min( Max(ccnc(mgs),cwnccn(mgs)), rho0(mgs)*qvex/Max( cwmasn5, Max(cwmasx,xmas(mgs,lc)) ) ) -! cn(mgs) = 1.5*cxmin - ELSEIF ( imaxsupopt == 4 ) THEN - cn(mgs) = Min( Max(ccnc(mgs),cwnccn(mgs)), rho0(mgs)*qvex/Max( cwmasn5, Max(cwmas20,xmas(mgs,lc)) ) ) + thetap(mgs) = thetap(mgs) + felvcp(mgs)*qvex/(pi0(mgs)) + IF ( io_flag .and. nxtra > 1 ) THEN + axtra(igs(mgs),jy,kgs(mgs),1) = axtra(igs(mgs),jy,kgs(mgs),1) + qvex/dtp ENDIF - IF ( lccna > 1 ) THEN - ccna(mgs) = ccna(mgs) + cn(mgs) - ELSE - ccnc(mgs) = Max( 0.0, ccnc(mgs) - cn(mgs) ) - ENDIF - cx(mgs,lc) = cx(mgs,lc) + cn(mgs) - ENDIF - -! write(iunit,*) 'theta = ',theta0(mgs) + thetap(mgs) + qwvp(mgs) = qwvp(mgs) - qvex + qx(mgs,lc) = qx(mgs,lc) + qvex + IF ( .not. flag_qndrop) THEN + IF ( imaxsupopt == 1 ) THEN + cn(mgs) = Min( Max(ccnc(mgs),cwnccn(mgs)), rho0(mgs)*qvex/Max( cwmasn5, xmas(mgs,lc) ) ) + ELSEIF ( imaxsupopt == 2 ) THEN + cn(mgs) = Min( Max(ccnc(mgs),cwnccn(mgs)), rho0(mgs)*qvex/Max( cwmasn5, Max(cwmas30,xmas(mgs,lc)) ) ) + ELSEIF ( imaxsupopt == 3 ) THEN + cn(mgs) = Min( Max(ccnc(mgs),cwnccn(mgs)), rho0(mgs)*qvex/Max( cwmasn5, Max(cwmasx,xmas(mgs,lc)) ) ) +! cn(mgs) = 1.5*cxmin + ELSEIF ( imaxsupopt == 4 ) THEN + cn(mgs) = Min( Max(ccnc(mgs),cwnccn(mgs)), rho0(mgs)*qvex/Max( cwmasn5, Max(cwmas20,xmas(mgs,lc)) ) ) + ENDIF + + IF ( lccna > 1 ) THEN + !IF ( ac_opt == 0 ) THEN + ccna(mgs) = ccna(mgs) + cn(mgs) + !ENDIF + ELSE + ccnc(mgs) = Max( 0.0, ccnc(mgs) - cn(mgs) ) + ENDIF -! temg(mgs) = theta(mgs)*( pres(mgs) / poo ) ** cap + cx(mgs,lc) = cx(mgs,lc) + cn(mgs) - ENDIF + ENDIF ! flag_qndrop - - ENDIF + ENDIF ! ( qvex .gt. 0.0 ) + + ENDIF ! ( qv1 .gt. (ssmx*qvs1) ) ! ! Calculate droplet volume and check if it is within bounds. @@ -11594,7 +11876,6 @@ SUBROUTINE NUCOND & xmas(mgs,lc) = Max( xmas(mgs,lc), cwmasn ) cx(mgs,lc) = rho0(mgs)*qx(mgs,lc)/xmas(mgs,lc) ! IF ( cx(mgs,lc) > tmp*1.1 ) THEN -! write(0,*) 'nucond: kgs, ccw1,2 = ',kgs(mgs),tmp,cx(mgs,lc) ! ENDIF ENDIF ENDIF @@ -11676,11 +11957,27 @@ SUBROUTINE NUCOND & IF ( ipconc .ge. 2 ) THEN an(igs(mgs),jy,kgs(mgs),lnc) = Max(cx(mgs,lc) , 0.0) + ! IF ( ac_opt > 10 .and. (cx(mgs,lc) > 0. .or. ccna(mgs) > 0. ) ) THEN + ! write(0,*) 'i,k final cx/cna = ',igs(mgs),kgs(mgs),cx(mgs,lc),ccna(mgs) + ! ENDIF + IF ( lss > 1 ) an(igs(mgs),jy,kgs(mgs),lss) = Max( 0.0, ssmax(mgs) ) IF ( ac_opt == 0 ) THEN IF ( lccn .gt. 1 .and. lccna .lt. 1 ) THEN an(igs(mgs),jy,kgs(mgs),lccn) = Max(0.0, ccnc(mgs) ) ENDIF + ELSEIF ( ac_opt == 1 .and. lccn > 1) THEN + an(igs(mgs),jy,kgs(mgs),lccn) = Max( 0.0, ccnc(mgs) ) ! cn are depleted for ac_opt=1 or 2 + ELSEIF ( ac_opt == 11 .and. lccna > 1) THEN + ! an(igs(mgs),jy,kgs(mgs),lccna) = Max( 0.0, ccna(mgs) ) ! done below + ELSEIF ( ac_opt == 2 .and. lccn > 1) THEN + an(igs(mgs),jy,kgs(mgs),lccn) = Max( 0.0, ccnc(mgs) ) + an(igs(mgs),jy,kgs(mgs),lcn_nu) = Max( 0.0, ccnc_nu(mgs) ) + an(igs(mgs),jy,kgs(mgs),lcn_co) = Max( 0.0, ccnc_co(mgs) ) + ELSEIF ( ac_opt == 22 .and. lccna > 1) THEN + ! an(igs(mgs),jy,kgs(mgs),lccna) = Max( 0.0, ccna(mgs) ) ! done below + an(igs(mgs),jy,kgs(mgs),lccnanu) = Max( 0.0, ccnanu(mgs) ) + an(igs(mgs),jy,kgs(mgs),lccnaco) = Max( 0.0, ccnaco(mgs) ) ENDIF IF ( lccnuf .gt. 1 .and. .not. ( lccna .gt. 1 .and. i_uf_or_ccn > 0 ) ) THEN an(igs(mgs),jy,kgs(mgs),lccnuf) = Max(0.0, ccncuf(mgs) ) @@ -11688,7 +11985,7 @@ SUBROUTINE NUCOND & IF ( lccna .gt. 1 ) THEN an(igs(mgs),jy,kgs(mgs),lccna) = Max(0.0, ccna(mgs) ) ENDIF - ENDIF + ENDIF ! ipconc >= 2 IF ( ipconc .ge. 3 .and. rcond == 2 ) THEN an(igs(mgs),jy,kgs(mgs),lnr) = Max(cx(mgs,lr) , 0.0) ENDIF @@ -11721,49 +12018,108 @@ SUBROUTINE NUCOND & ! end of gather scatter (for this jy slice) -!#ifdef COMMAS -! GOTO 9999 -!#endif - ! Redistribute inappreciable cloud particles and charge ! ! Redistribution everywhere in the domain... ! - IF ( .true. ) THEN - - frac = 1.0 ! 0.25 ! 1.0 ! 0.2 +! moved to separate subroutine (below) ! -! alternate test version for ipconc .ge. 3 -! just vaporize stuff to prevent noise in the number concentrations + + + 9999 RETURN + + END SUBROUTINE NUCOND - do kz = 1,nz -! do jy = 1,1 - do ix = 1,nxi - - t0(ix,jy,kz) = an(ix,jy,kz,lt)*t77(ix,jy,kz) - - zerocx(:) = .false. - DO il = lc,lhab - IF ( iresetmoments == 1 .or. iresetmoments == il ) THEN - IF ( ln(il) > 1 ) zerocx(il) = ( an(ix,jy,kz,ln(il)) < cxmin ) - IF ( lz(il) > 1 ) zerocx(il) = ( zerocx(il) .or. an(ix,jy,kz,lz(il)) < zxmin ) - ELSE - IF ( il == lc ) THEN - IF ( ln(il) > 1 ) zerocx(il) = ( an(ix,jy,kz,ln(il)) <= 0 ) .and. .not. flag_qndrop ! do not reset if progn=1 (WRF-CHEM) - ELSE - IF ( ln(il) > 1 ) zerocx(il) = ( an(ix,jy,kz,ln(il)) <= 0 ) - ENDIF - ENDIF - ENDDO - - IF ( lhl .gt. 1 ) THEN - - IF ( lzhl .gt. 1 ) THEN - an(ix,jy,kz,lzhl) = Max(0.0, an(ix,jy,kz,lzhl) ) - - IF ( an(ix,jy,kz,lhl) .ge. frac*qxmin(lhl) .and. rescale_low_alpha ) THEN ! check 6th moment +! ##################################################################### +! ##################################################################### +! Clean up tiny values of mixing ratio +! Redistribute inappreciable cloud particles and charge +! +! Redistribution everywhere in the domain... +! + subroutine smallvalues & + & (nx,ny,nz,na,jyslab & + & ,nor,norz,dtp,nxi & + & ,t0 & + & ,an,dn, w & + & ,t77,flag_qndrop & + & ) + + + implicit none + + integer :: nx,ny,nz,na,nxi + integer :: nor,norz, jyslab ! ,nht,ngt,igsr + real :: dtp ! time step + logical,intent(in) :: flag_qndrop + +! +! external temporary arrays +! + real t77(-nor+1:nx+nor,-nor+1:ny+nor,-norz+1:nz+norz) + real t0(-nor+1:nx+nor,-nor+1:ny+nor,-norz+1:nz+norz) + real an(-nor+1:nx+nor,-nor+1:ny+nor,-norz+1:nz+norz,na) + real dn(-nor+1:nx+nor,-nor+1:ny+nor,-norz+1:nz+norz) + real w(-nor+1:nx+nor,-nor+1:ny+nor,-norz+1:nz+norz) + + ! local + + + logical zerocx(lc:lqmx) + + real :: frac, hwdn, tmpg, xdia1, xdia3, cwch,xvol + + integer ix,kz,i,n, km1 + integer :: il + integer :: jy, jgs + real :: chw, g1, z1, tmp, tmp2, fw, tmpmx, qr + + +! Redistribute inappreciable cloud particles and charge +! +! Redistribution everywhere in the domain... +! + jy = 1 + + frac = 1.0 ! 0.25 ! 1.0 ! 0.2 + + cwch = ((3. + alphah)*(2. + alphah)*(1.0 + alphah))**(-1./3.) +! +! alternate test version for ipconc .ge. 3 +! just vaporize stuff to prevent noise in the number concentrations + + + do kz = 1,nz +! do jy = 1,1 + do ix = 1,nxi + + t0(ix,jy,kz) = an(ix,jy,kz,lt)*t77(ix,jy,kz) + + zerocx(:) = .false. + DO il = lc,lhab + IF ( iresetmoments == 1 .or. iresetmoments == il ) THEN + IF ( ln(il) > 1 ) zerocx(il) = ( an(ix,jy,kz,ln(il)) < cxmin ) + IF ( lz(il) > 1 ) zerocx(il) = ( zerocx(il) .or. (an(ix,jy,kz,lz(il)) < zxmin) ) + ELSE + IF ( il == lc ) THEN + IF ( ln(il) > 1 ) THEN + zerocx(il) = ( an(ix,jy,kz,ln(il)) <= 0.0 ) .and. .not. flag_qndrop ! do not reset if progn=1 (WRF-CHEM) + ENDIF + ELSE + IF ( ln(il) > 1 ) zerocx(il) = ( an(ix,jy,kz,ln(il)) <= 0.0 ) + ENDIF + ENDIF + ENDDO + + IF ( lhl .gt. 1 ) THEN + + IF ( lzhl .gt. 1 ) THEN + + an(ix,jy,kz,lzhl) = Max(0.0, an(ix,jy,kz,lzhl) ) + + IF ( an(ix,jy,kz,lhl) .ge. frac*qxmin(lhl) .and. rescale_low_alpha ) THEN ! check 6th moment IF ( an(ix,jy,kz,lnhl) .gt. 0.0 ) THEN @@ -11796,7 +12152,7 @@ SUBROUTINE NUCOND & ENDIF !lzhl - if ( an(ix,jy,kz,lhl) .lt. frac*qxmin(lhl) .or. zerocx(lhl) ) then + if ( (an(ix,jy,kz,lhl) .lt. frac*qxmin(lhl)) .or. zerocx(lhl) ) then ! IF ( an(ix,jy,kz,lhl) .gt. 0 ) THEN an(ix,jy,kz,lv) = an(ix,jy,kz,lv) + an(ix,jy,kz,lhl) @@ -11870,6 +12226,26 @@ SUBROUTINE NUCOND & ENDIF + IF ( lvhl .gt. 1 ) THEN + IF ( an(ix,jy,kz,lvhl) .gt. 0.0 ) THEN + hwdn = dn(ix,jy,kz)*an(ix,jy,kz,lhl)/an(ix,jy,kz,lvhl) + ELSE + hwdn = xdn0(lhl) + ENDIF + hwdn = Max( xdnmn(lhl), hwdn ) + ELSE + hwdn = xdn0(lhl) + ENDIF + + qr = an(ix,jy,kz,lhl) + xvol = dn(ix,jy,kz)*an(ix,jy,kz,lhl)/(hwdn*an(ix,jy,kz,lnhl)) + chw = an(ix,jy,kz,lnhl) + + IF ( xvol .lt. xvmn(lhl) .or. xvol .gt. xvmx(lhl) ) THEN + xvol = Min( xvmx(lhl), Max( xvmn(lhl),xvol ) ) + chw = dn(ix,jy,kz)*an(ix,jy,kz,lhl)/(xvol*hwdn) + an(ix,jy,kz,lnhl) = chw + ENDIF ! CHECK INTERCEPT IF ( ipconc == 5 .and. an(ix,jy,kz,lhl) .gt. qxmin(lhl) .and. alphahl .le. 0.1 .and. lnhl .gt. 1 .and. lzhl == 0 ) THEN @@ -11880,9 +12256,9 @@ SUBROUTINE NUCOND & hwdn = xdn0(lhl) ENDIF tmp = (hwdn*an(ix,jy,kz,lnhl))/(dn(ix,jy,kz)*an(ix,jy,kz,lhl)) - tmpg = an(ix,jy,kz,lnhl)*(tmp*(3.14159))**(1./3.) + tmpg = an(ix,jy,kz,lnhl)*(tmp*pi)**(1./3.) IF ( tmpg .lt. cnohlmn ) THEN - tmp = ( (hwdn)/(dn(ix,jy,kz)*an(ix,jy,kz,lhl))*(3.14159))**(1./3.) + tmp = ( (hwdn)/(dn(ix,jy,kz)*an(ix,jy,kz,lhl))*pi)**(1./3.) an(ix,jy,kz,lnhl) = (cnohlmn/tmp)**(3./4.) ENDIF @@ -11932,7 +12308,7 @@ SUBROUTINE NUCOND & ENDIF - if ( an(ix,jy,kz,lh) .lt. frac*qxmin(lh) .or. zerocx(lh) ) then + if ( (an(ix,jy,kz,lh) .lt. frac*qxmin(lh)) .or. zerocx(lh) ) then ! IF ( an(ix,jy,kz,lh) .gt. 0 ) THEN an(ix,jy,kz,lv) = an(ix,jy,kz,lv) + an(ix,jy,kz,lh) @@ -12006,9 +12382,6 @@ SUBROUTINE NUCOND & ENDIF -! CHECK INTERCEPT - IF ( ipconc == 5 .and. an(ix,jy,kz,lh) .gt. qxmin(lh) .and. alphah .le. 0.1 .and. lnh .gt. 1 .and. lzh == 0 ) THEN - IF ( lvh .gt. 1 ) THEN IF ( an(ix,jy,kz,lvh) .gt. 0.0 ) THEN hwdn = dn(ix,jy,kz)*an(ix,jy,kz,lh)/an(ix,jy,kz,lvh) @@ -12019,22 +12392,49 @@ SUBROUTINE NUCOND & ELSE hwdn = xdn0(lh) ENDIF + + qr = an(ix,jy,kz,lh) + xvol = dn(ix,jy,kz)*an(ix,jy,kz,lh)/(hwdn*an(ix,jy,kz,lnh)) + chw = an(ix,jy,kz,lnh) + + IF ( xvol .lt. xvmn(lh) .or. xvol .gt. xvmx(lh) ) THEN + xvol = Min( xvmx(lh), Max( xvmn(lh),xvol ) ) + chw = dn(ix,jy,kz)*an(ix,jy,kz,lh)/(xvol*hwdn) + an(ix,jy,kz,lnh) = chw + ENDIF + +! CHECK INTERCEPT + IF ( ipconc == 5 .and. an(ix,jy,kz,lh) .gt. qxmin(lh) .and. alphah .le. 0.1 .and. lnh .gt. 1 .and. lzh == 0 ) THEN + tmp = (hwdn*an(ix,jy,kz,lnh))/(dn(ix,jy,kz)*an(ix,jy,kz,lh)) - tmpg = an(ix,jy,kz,lnh)*(tmp*(3.14159))**(1./3.) + tmpg = an(ix,jy,kz,lnh)*(tmp*pi)**(1./3.) IF ( tmpg .lt. cnohmn ) THEN ! tmpg = an(ix,jy,kz,lnh)*( (hwdn*an(ix,jy,kz,lnh))/(dn(ix,jy,kz)*an(ix,jy,kz,lh))*(3.14159))**(1./3.) ! tmpg = an(ix,jy,kz,lnh)**(4./3.)*( (hwdn)/(dn(ix,jy,kz)*an(ix,jy,kz,lh))*(3.14159))**(1./3.) - tmp = ( (hwdn)/(dn(ix,jy,kz)*an(ix,jy,kz,lh))*(3.14159))**(1./3.) + tmp = ( (hwdn)/(dn(ix,jy,kz)*an(ix,jy,kz,lh))*pi)**(1./3.) an(ix,jy,kz,lnh) = (cnohmn/tmp)**(3./4.) ENDIF ENDIF + + IF ( ipconc == 5 .and. imorrgdnglimit == 1 ) THEN + ! limit on characteristic diameter (i.e., 1/slope) + xdia3 = (xvol*6.*piinv)**(1./3.) + xdia1 = cwch*xdia3 + IF ( xdia1 > morrdnglimit ) THEN + xdia1 = morrdnglimit + xvol = pi/6.0*(xdia1/cwch)**3 + chw = dn(ix,jy,kz)*qr/(xvol*hwdn) + an(ix,jy,kz,lnh) = chw + xdia3 = (xvol*6.*piinv)**(1./3.) + ENDIF + + ENDIF end if - if ( an(ix,jy,kz,ls) .lt. frac*qxmin(ls) .or. zerocx(ls) & ! .or. an(ix,jy,kz,lns) .lt. 0.1 ! .and. - & ) then + if ( (an(ix,jy,kz,ls) .lt. frac*qxmin(ls)) .or. zerocx(ls) ) then IF ( t0(ix,jy,kz) .lt. 273.15 ) THEN ! IF ( an(ix,jy,kz,ls) .gt. 0 ) THEN an(ix,jy,kz,lv) = an(ix,jy,kz,lv) + an(ix,jy,kz,ls) @@ -12095,8 +12495,7 @@ SUBROUTINE NUCOND & an(ix,jy,kz,lzr) = Max(0.0, an(ix,jy,kz,lzr) ) ENDIF - if ( an(ix,jy,kz,lr) .lt. frac*qxmin(lr) .or. zerocx(lr) & - & ) then + if ( (an(ix,jy,kz,lr) .lt. frac*qxmin(lr)) .or. zerocx(lr) ) then an(ix,jy,kz,lv) = an(ix,jy,kz,lv) + an(ix,jy,kz,lr) an(ix,jy,kz,lr) = 0.0 IF ( ipconc .ge. 3 ) THEN @@ -12113,8 +12512,7 @@ SUBROUTINE NUCOND & ! ! for qci ! - IF ( an(ix,jy,kz,li) .le. frac*qxmin(li) .or. zerocx(li) & ! .or. an(ix,jy,kz,lni) .lt. 0.1 - & ) THEN + IF ( (an(ix,jy,kz,li) .le. frac*qxmin(li)) .or. zerocx(li) ) THEN an(ix,jy,kz,lv) = an(ix,jy,kz,lv) + an(ix,jy,kz,li) an(ix,jy,kz,li)= 0.0 IF ( ipconc .ge. 1 ) THEN @@ -12122,41 +12520,11 @@ SUBROUTINE NUCOND & ENDIF ENDIF -! -! for qis -! - IF ( lis > 1 ) THEN ! { - IF ( an(ix,jy,kz,lis) .le. frac*qxmin(lis) .or. zerocx(lis) & ! .or. an(ix,jy,kz,lni) .lt. 0.1 - & ) THEN ! { { - an(ix,jy,kz,lv) = an(ix,jy,kz,lv) + an(ix,jy,kz,lis) - an(ix,jy,kz,lis)= 0.0 - IF ( ipconc .ge. 1 ) THEN - an(ix,jy,kz,lnis) = 0.0 - ENDIF - - ELSEIF ( icespheres >= 2 ) THEN ! } { - km1 = Max(1, kz-1) - IF ( 0.5*( w(ix,jy,kz) + w(ix,jy,kz+1)) < -1.0 .or. & - & (icespheres == 3 .and. ( t0(ix,jy,kz) < 232.15 .or. an(ix,jy,kz,lc) < qxmin(lc) ) ) .or. & - & (icespheres == 5 .and. ( t0(ix,jy,kz) < 232.15 .or. & - & ( an(ix,jy,kz,lc) < qxmin(lc) .and. an(ix,jy,km1,lc) < qxmin(lc) )) ) .or. & - & (icespheres == 4 .and. ( t0(ix,jy,kz) < 235.15 )) ) THEN ! transfer to regular ice crystals in downdraft or at low temp - an(ix,jy,kz,li) = an(ix,jy,kz,li) + an(ix,jy,kz,lis) - an(ix,jy,kz,lni) = an(ix,jy,kz,lni) + an(ix,jy,kz,lnis) - an(ix,jy,kz,lis)= 0.0 - an(ix,jy,kz,lnis)= 0.0 - - ENDIF - - ENDIF ! } } - ENDIF ! } - ! ! for qcw ! - IF ( an(ix,jy,kz,lc) .le. frac*qxmin(lc) .or. zerocx(lc) & - & ) THEN + IF ( (an(ix,jy,kz,lc) .le. frac*qxmin(lc)) .or. zerocx(lc) ) THEN an(ix,jy,kz,lv) = an(ix,jy,kz,lv) + an(ix,jy,kz,lc) an(ix,jy,kz,lc)= 0.0 IF ( ipconc .ge. 2 ) THEN @@ -12164,19 +12532,36 @@ SUBROUTINE NUCOND & IF ( irenuc < 5 .and. lccna <= 1 ) THEN IF ( ac_opt == 0 ) THEN an(ix,jy,kz,lccn) = an(ix,jy,kz,lccn) + Max(0.0,an(ix,jy,kz,lnc)) + ELSEIF ( lccn > 1 ) THEN + an(ix,jy,kz,lccn) = an(ix,jy,kz,lccn) + Max(0.0,an(ix,jy,kz,lnc)) ENDIF ELSEIF ( lccna > 1 ) THEN - an(ix,jy,kz,lccna) = Max( 0.0, an(ix,jy,kz,lccna) - Max(0.0,an(ix,jy,kz,lnc)) ) + tmp = Max(0.0,an(ix,jy,kz,lnc)) + IF ( lccnaco > 1 .and. lccnanu > 1 ) THEN + ! restore CCN proportionally to each type, although coarse are presumably already lost to rain + tmp2 = an(ix,jy,kz,lccna) + an(ix,jy,kz,lccnaco) + an(ix,jy,kz,lccnanu) + IF ( tmp2 > 0.0 .and. tmp > 0.0 ) THEN + an(ix,jy,kz,lccna) = Max( 0.0, an(ix,jy,kz,lccna) - tmp*an(ix,jy,kz,lccna)/tmp2 ) + an(ix,jy,kz,lccnaco) = Max( 0.0, an(ix,jy,kz,lccnaco) - tmp*an(ix,jy,kz,lccnaco)/tmp2 ) + an(ix,jy,kz,lccnanu) = Max( 0.0, an(ix,jy,kz,lccnanu) - tmp*an(ix,jy,kz,lccnanu)/tmp2 ) + ENDIF + ELSE + an(ix,jy,kz,lccna) = Max( 0.0, an(ix,jy,kz,lccna) - tmp ) + ENDIF ENDIF ENDIF an(ix,jy,kz,lnc) = 0.0 IF ( lccn > 1 ) an(ix,jy,kz,lccn) = Max( 0.0, an(ix,jy,kz,lccn) ) - IF ( lccna > 0 .and. ac_opt == 0 ) THEN ! apply exponential decay to activated CCN to restore to environmental value +! IF ( lccna > 0 .and. ac_opt == 0 ) THEN ! apply exponential decay to activated CCN to restore to environmental value + IF ( lccna > 0 ) THEN ! apply exponential decay to activated CCN to restore to environmental value IF ( restoreccn ) THEN tmp = an(ix,jy,kz,li) + an(ix,jy,kz,ls) - - IF ( an(ix,jy,kz,lccna) > 1. .and. tmp < qxmin(li) ) an(ix,jy,kz,lccna) = an(ix,jy,kz,lccna)*Exp(-dtp/ccntimeconst) + IF ( tmp < qxmin(li) ) THEN + IF ( an(ix,jy,kz,lccna) > 1. .and. tmp < qxmin(li) ) an(ix,jy,kz,lccna) = an(ix,jy,kz,lccna)*Exp(-dtp/ccntimeconst) + IF ( lccnaco > 1 ) an(ix,jy,kz,lccnaco) = an(ix,jy,kz,lccnaco)*Exp(-dtp/ccntimeconst) + IF ( lccnanu > 1 ) an(ix,jy,kz,lccnanu) = an(ix,jy,kz,lccnanu)*Exp(-dtp/ccntimeconst) + ENDIF ENDIF ELSEIF ( lccn > 1 .and. restoreccn .and. ac_opt == 0 ) THEN ! in this case, we are treating the ccn field as ccna @@ -12203,20 +12588,8 @@ SUBROUTINE NUCOND & ! end do end do - ENDIF ! true/false - IF ( ndebug .ge. 1 ) write(6,*) 'END OF ICEZVD_DR' -! -! - - - 9999 RETURN - - END SUBROUTINE NUCOND - - -! ##################################################################### -! ##################################################################### + end subroutine smallvalues @@ -12329,7 +12702,8 @@ subroutine nssl_2mom_gs & real rainprod2d(-nor+1:nx+nor,-norz+ng1:nz+norz) real evapprod2d(-nor+1:nx+nor,-norz+ng1:nz+norz) - real alpha2d(-nor+1:nx+nor,-norz+ng1:nz+norz,3) + + real :: alpha2d(-nor+1:nx+nor,-norz+ng1:nz+norz,3) real, parameter :: tfrdry = 243.15 @@ -12401,7 +12775,7 @@ subroutine nssl_2mom_gs & integer i,j,k,i1 integer kzb,kze real slope1, slope2 - real x1, x2, x3 + real x1, x2, x3, y1 real eps,eps2 parameter (eps=1.e-20,eps2=1.e-5) ! @@ -12536,12 +12910,12 @@ subroutine nssl_2mom_gs & real, parameter :: mfrag = 1.0e-10 ! assumed ice fragment mass for collisional splintering (Schuur & Rutledge 00b) double precision cautn(ngs), rh(ngs), nh(ngs) real ex1, ft, rhoinv(ngs) - double precision ec0(ngs) + real :: ec0(ngs) - real ac1,bc, taus, c1,d1,e1,f1,p380,tmp,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,temp3 ! , sstdy, super - real :: flim + real ac1,bc, c1,d1,e1,f1,p380,tmp,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,temp3 ! , sstdy, super + real :: flim, xmass real dw,dwr - double precision :: tmpz, tmpzmlt + double precision :: tmpc, tmpz, tmpzmlt real ratio, delx, dely real dbigg,volt real chgtmp,fac,mixedphasefac @@ -12712,7 +13086,7 @@ subroutine nssl_2mom_gs & real :: zxmxd(ngs,lr:lhab) real :: g1x(ngs,lr:lhab) - + real :: g1xmax,g1xmin real :: qsimxdep(ngs) ! max sublimation of qi+qs+qis real :: qsimxsub(ngs) ! max depositionof qi+qs+qis logical,parameter :: DoSublimationFix = .true. @@ -12735,9 +13109,10 @@ subroutine nssl_2mom_gs & real :: qhgt10mm ! mass greater than 10mm real :: qhgt20mm ! mass greater than 20mm real :: fwmhtmp - real, parameter :: fwmhtmptem = -15. ! temperature at which fwmhtmp fully switches to liquid water only being on large particles +! real, parameter :: fwmhtmptem = -15. ! temperature at which fwmhtmp fully switches to liquid water only being on large particles real, parameter :: d1t = (6.0 * 0.268e-3/(917.* pi))**(1./3.) ! d1t is the diameter of the ice sphere with the mass (0.268e-3 kg) of an 8mm spherical drop real, parameter :: srasheym = 0.1389 ! slope fraction from Rasmussen and Heymsfield + real :: dtmp ! real swvent(ngs),hwvent(ngs),rwvent(ngs),hlvent(ngs),hwventy(ngs),hlventy(ngs),rwventz(ngs) real hxventtmp @@ -12757,6 +13132,7 @@ subroutine nssl_2mom_gs & real qhlmlr0, qhlmlr05, qhlmlr1, qhlmlr2,qhlmlr3, qhlmlr12, qhlmlr23 real qxd1, cxd1, zxd1 ! mass and number up to mltdiam1 real qxd05, cxd05 ! mass and number up to mltdiam1/2 + real :: qrbreak, crbreaksmall, crbreak, zrbreak, breakbin real :: qxd(ndiam+4), cxd(ndiam+4), qhml(ndiam+4), qhml0(ndiam+4) real :: dqxd(ndiam+4), dcxd(ndiam+4), dqhml(ndiam+4) @@ -12887,6 +13263,7 @@ subroutine nssl_2mom_gs & ! real qsaci(ngs) real qsacis(ngs) + real csacis(ngs) real qhaci(ngs) real qhacs(ngs) @@ -12895,6 +13272,7 @@ subroutine nssl_2mom_gs & real :: chacis0(ngs) real :: csaci0(ngs) ! collision rate only + real :: csacis0(ngs) ! collision rate only real :: chaci0(ngs) ! collision rate only real :: chacs0(ngs) ! collision rate only real :: chlaci0(ngs) @@ -13045,7 +13423,7 @@ subroutine nssl_2mom_gs & real ehxw(ngs),ehlw(ngs),egmw(ngs),ehw(ngs) real err(ngs),esr(ngs),eglr(ngs),eghr(ngs),efr(ngs) real ehxr(ngs),ehlr(ngs),egmr(ngs) - real eri(ngs),esi(ngs),egli(ngs),eghi(ngs),efi(ngs),efis(ngs) + real eri(ngs),esi(ngs),esis(ngs),egli(ngs),eghi(ngs),efi(ngs),efis(ngs) real ehxi(ngs),ehli(ngs),egmi(ngs),ehi(ngs),ehis(ngs),ehlis(ngs) real ers(ngs),ess(ngs),egls(ngs),eghs(ngs),efs(ngs),ehs(ngs),ehsfac(ngs) real ehscnv(ngs) @@ -13054,7 +13432,7 @@ subroutine nssl_2mom_gs & real ehsclsn(ngs),ehiclsn(ngs),ehisclsn(ngs) real efsclsn(ngs),eficlsn(ngs),efisclsn(ngs) real ehlsclsn(ngs),ehliclsn(ngs),ehlisclsn(ngs) - real esiclsn(ngs) + real esiclsn(ngs),esisclsn(ngs) real :: ehs_collsn = 0.5, ehi_collsn = 1.0 real :: efs_collsn = 0.5, efi_collsn = 1.0 @@ -13634,6 +14012,7 @@ subroutine nssl_2mom_gs & do ix = nxmpb,itile pqs(1) = t00(ix,jy,kz) + pres(1) = pn(ix,jy,kz) + pb(kz) theta(1) = an(ix,jy,kz,lt) temg(1) = t0(ix,jy,kz) @@ -13641,7 +14020,12 @@ subroutine nssl_2mom_gs & tqvcon = temg(1)-cbw ltemq = (temg(1)-163.15)/fqsat + 1.5 ltemq = Min( nqsat, Max(1,ltemq) ) - qvs(1) = pqs(1)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + qvs(1) = pqs(1)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + qvs(1) = rdorv*esbolton*tabqvs(ltemq)/(pres(1) - esbolton*tabqvs(ltemq)) + ENDIF + IF ( iqis0 == 1 .or. temg(1) <= tfr+0.5 ) THEN qis(1) = pqs(1)*tabqis(ltemq) ELSE @@ -13727,7 +14111,13 @@ subroutine nssl_2mom_gs & pqs(mgs) = (380.0)/(pres(mgs)) ltemq = (temg(mgs)-163.15)/fqsat+1.5 ltemq = Min( nqsat, Max(1,ltemq) ) - qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + + IF ( iqvsopt == 0 ) THEN + qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + qvs(mgs) = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF + IF ( iqis0 == 1 .or. temg(mgs) <= tfr+0.5 ) THEN qis(mgs) = pqs(mgs)*tabqis(ltemq) ELSE @@ -13960,6 +14350,12 @@ subroutine nssl_2mom_gs & IF ( ipconc .ge. 6 ) THEN + tmp = alphamax - 1.0 + g1xmax = (6.0 + tmp)*(5.0 + tmp)*(4.0 + tmp)/ & + & ((3.0 + tmp)*(2.0 + tmp)*(1.0 + tmp)) + g1xmin = (6.0 + alphamin)*(5.0 + alphamin)*(4.0 + alphamin)/ & + & ((3.0 + alphamin)*(2.0 + alphamin)*(1.0 + alphamin)) + IF ( lz(lr) .lt. 1 ) THEN g1x(:,lr) = (6.0 + alphar)*(5.0 + alphar)*(4.0 + alphar)/ & & ((3.0 + alphar)*(2.0 + alphar)*(1.0 + alphar)) @@ -13984,6 +14380,16 @@ subroutine nssl_2mom_gs & ENDIF + IF ( ipconc == 5 ) THEN + ! set up factors for ihlcnh=3 conversion + g1x(:,lr) = (6.0 + alphar)*(5.0 + alphar)*(4.0 + alphar)/ & + & ((3.0 + alphar)*(2.0 + alphar)*(1.0 + alphar)) + g1x(:,lh) = (6.0 + alphah)*(5.0 + alphah)*(4.0 + alphah)/ & + & ((3.0 + alphah)*(2.0 + alphah)*(1.0 + alphah)) + g1x(:,lhl) = (6.0 + alphahl)*(5.0 + alphahl)*(4.0 + alphahl)/ & + & ((3.0 + alphahl)*(2.0 + alphahl)*(1.0 + alphahl)) + ENDIF + scx(:,:) = 0.0 ! ! set shape parameters @@ -14098,7 +14504,6 @@ subroutine nssl_2mom_gs & tmp = qx(mgs,li)+qx(mgs,ls)+qx(mgs,lh) IF ( lhl > 1 ) tmp = tmp + qx(mgs,lhl) - IF ( lf > 1 ) tmp = tmp + qx(mgs,lf) cvm = cv+cvv*qx(mgs,lv)+cpl*(qx(mgs,lc)+qx(mgs,lr)) & +cpigb*(tmp) @@ -14259,7 +14664,7 @@ subroutine nssl_2mom_gs & ! alpha(mgs,lr) = Min(alphamax, c1r*tanh(c2r*(xdia(mgs,lr,3)*1000. - c3r)) + c4r) ! IF ( igs(mgs) == 19 ) write(0,*) 'imy: i,k,alpr,xdia = ',igs(mgs),kgs(mgs),alpha(mgs,lr),xdia(mgs,lr,3)*1000. - ! M&M-C 2010: + ! Milbrandt & M-C 2010: tmp = 4. + alphar i = Int(dgami*(tmp)) del = tmp - dgam*i @@ -14280,7 +14685,7 @@ subroutine nssl_2mom_gs & xdia(mgs,lh,3) = (xv(mgs,lh)*6.*piinv)**(1./3.) ! mwfac*xdia(mgs,lh,1) ! (xv(mgs,lh)*cwc0*6.0)**(1./3.) ! alpha(mgs,lh) = Min(alphamax, c1h*tanh(c2h*(xdia(mgs,lh,3)*1000. - c3h)) + c4h) - ! M&M-C 2010: + ! Milbrandt & M-C 2010: tmp = 4. + dnu(lh) i = Int(dgami*(tmp)) del = tmp - dgam*i @@ -14802,7 +15207,7 @@ subroutine nssl_2mom_gs & ! check for artificial breakup (graupel/hail larger than allowed max size) - IF ( imaxdiaopt == 1 ) THEN + IF ( imaxdiaopt == 1 .or. il /= lr ) THEN xvbarmax = xvmx(il) ELSEIF ( imaxdiaopt == 2 ) THEN ! test against maximum mass diameter xvbarmax = xvmx(il) /((3. + alpha(mgs,il))**3/((3. + alpha(mgs,il))*(2. + alpha(mgs,il))*(1. + alpha(mgs,il)))) @@ -14829,12 +15234,21 @@ subroutine nssl_2mom_gs & IF ( tmp < cx(mgs,il) ) THEN ! artificial breakup has happened, so need to adjust reflectivity and find new shape parameter g1 = 36.*(6.0 + alpha(mgs,il))*(5.0 + alpha(mgs,il))*(4.0 + alpha(mgs,il))/ & & ((3.0 + alpha(mgs,il))*(2.0 + alpha(mgs,il))*(1.0 + alpha(mgs,il))*pi**2) - zx(mgs,il) = zx(mgs,il) + g1*(rho0(mgs)/xdn(mgs,il))**2*( (qx(mgs,il)/tmp)**2 * (tmp-cx(mgs,il)) ) +! zx(mgs,il) = zx(mgs,il) + g1*(rho0(mgs)/xdn(mgs,il))**2*( (qx(mgs,il)/tmp)**2 * (tmp-cx(mgs,il)) ) + ! check if incoming zx is consistent + ! Z from incoming cx, qx, and alpha + tmpz = g1/(pi/6.*xdn(mgs,il))**2 * ((rho0(mgs)*qx(mgs,il))**2)/tmp + IF ( tmpz > zx(mgs,il) ) THEN + tmpc = g1/(pi/6.*xdn(mgs,il))**2 * ((rho0(mgs)*qx(mgs,il))**2)/zx(mgs,il) + cx(mgs,il) = Max(cx(mgs,il), tmpc) + ! find cx that gives zx + ENDIF + zx(mgs,il) = g1/(pi/6.*xdn(mgs,il))**2 * ((rho0(mgs)*qx(mgs,il))**2)/cx(mgs,il) an(igs(mgs),jgs,kgs(mgs),lz(il)) = zx(mgs,il) - chw = cx(mgs,il) - qr = qx(mgs,il) - z = zx(mgs,il) + qr = qx(mgs,il) + chw = cx(mgs,il) + z = zx(mgs,il) rdi = z*(pi/6.*xdn(mgs,il))**2*chw/((rho0(mgs)*qr)**2) alp = (6.0+alpha(mgs,il))*(5.0+alpha(mgs,il))*(4.0+alpha(mgs,il))/ & @@ -14948,12 +15362,12 @@ subroutine nssl_2mom_gs & gf1palp(mgs) = y IF ( iferwisventr == 2 ) THEN +! ventrn = Gamma(alphar + 2.5 + br/2.)/Gamma(alphar + 1.) ! adapted from Wisner et al. 1972 tmp = alpha(mgs,lr) + 2.5 + br/2. i = Int(dgami*(tmp)) del = tmp - dgam*i x = gmoi(i) + (gmoi(i+1) - gmoi(i))*del*dgami -! ventrx(mgs) = Gamma_sp(alpha(mgs,lr) + 1.5 + br/6.)/Gamma_sp(alpha(mgs,lr) + 1.) ventrxn(mgs) = x/y @@ -15128,9 +15542,9 @@ subroutine nssl_2mom_gs & call setvtz(ngscnt,qx,qxmin,qxw,cx,rho0,rhovt,xdia,cno,cnostmp, & & xmas,vtxbar,xdn,xvmn,xvmx,xv,cdx,cdxgs, & - & ipconc,ndebug,ngs,nz,kgs,fadvisc, & + & ipconc,ndebug,ngs,nz,igs,kgs,fadvisc, & & cwmasn,cwmasx,cwradn,cnina,cimn,cimx, & - & itype1,itype2,temcg,infdo,alpha,0,axx,bxx) ! ,cdh,cdhl) + & itype1,itype2,temcg,infdo,alpha,axx,bxx,0) ! ,cdh,cdhl) ! & itype1,itype2,temcg,infdo,alpha,0,axh,bxh,axhl,bxhl) ! ,cdh,cdhl) @@ -15336,8 +15750,9 @@ subroutine nssl_2mom_gs & IF ( qx(mgs,il) > qxmin(il) .and. ivshdgs > 0 ) THEN ! tmpdiam is weighted diameter of d^(shedalp-1), so for shedalp=3, this is the area-weighted diameter or maximum mass diameter. - tmpdiam = (shedalp+alpha(mgs,il))*xdia(mgs,il,1) ! *( xdn(mgs,il)/917. )**(1./3.) ! erm added density factor for equiv. solid ice sphere 10.12.2015 - + ! tmpdiam = (shedalp+alpha(mgs,il))*xdia(mgs,il,1) ! *( xdn(mgs,il)/917. )**(1./3.) ! erm added density factor for equiv. solid ice sphere 10.12.2015 + tmpdiam = (shedalp+0.0)*xdia(mgs,il,1) ! *( xdn(mgs,il)/917. )**(1./3.) ! erm added density factor for equiv. solid ice sphere 10.12.2015 + ! imltshddmr IF ( tmpdiam > sheddiam0 ) THEN vshdgs(mgs,il) = 0.523599*(1.5e-3)**3/massfacshr ! 1.5mm drops from very large ice ELSEIF ( tmpdiam > sheddiam ) THEN ! intermediate size @@ -15497,7 +15912,7 @@ subroutine nssl_2mom_gs & if (xdia(mgs,lc,1).gt.ewi_dcmin .and. xdia(mgs,li,1).gt.ewi_dimin) then ! erm 5/10/2007 test following change: ! if (xdia(mgs,lc,1).gt.12.0e-06 .and. xdia(mgs,li,1).gt.50.0e-06) then - eiw(mgs) = 0.5 + eiw(mgs) = eiw0 end if if ( temg(mgs) .ge. 273.15 ) eiw(mgs) = 0.0 end if @@ -15672,6 +16087,7 @@ subroutine nssl_2mom_gs & ! ENDIF if ( temg(mgs) .gt. 273.15 ) esi(mgs) = 0.0 end if + ! ! ! @@ -16370,24 +16786,6 @@ subroutine nssl_2mom_gs & end do - IF ( lis > 1 .and. ipconc >= 5 ) THEN - do mgs = 1,ngscnt - qhacis(mgs) = 0.0 - qhacis0(mgs) = 0.0 - IF ( ehis(mgs) .gt. 0.0 ) THEN - - vt = Sqrt((vtxbar(mgs,lh,1)-vtxbar(mgs,lis,1))**2 + & - & 0.04*vtxbar(mgs,lh,1)*vtxbar(mgs,lis,1) ) - - qhacis0(mgs) = 0.25*pi*ehisclsn(mgs)*cx(mgs,lh)*qx(mgs,lis)*vt* & - & ( da0lh(mgs)*xdia(mgs,lh,3)**2 + & - & dab1lh(mgs,lh,lis)*xdia(mgs,lh,3)*xdia(mgs,lis,3) + & - & da1(li)*xdia(mgs,lis,3)**2 ) - qhacis(mgs) = Min( ehis(mgs)*qhacis0(mgs), qxmxd(mgs,lis) ) - ENDIF - end do - ENDIF - ! ! do mgs = 1,ngscnt @@ -16638,6 +17036,9 @@ subroutine nssl_2mom_gs & end do ENDIF ! + qhlacis(:) = 0.0 + qhlacis0(:) = 0.0 + qhlacs(:) = 0.0 qhlacs0(:) = 0.0 IF ( lhl .gt. 1 ) THEN @@ -16759,7 +17160,7 @@ subroutine nssl_2mom_gs & ni = ni + cx(mgs,li)*Exp(- (40.e-6/xdia(mgs,li,1))**3 ) ENDIF IF ( imurain == 1 ) THEN ! gamma of diameter - IF ( iacrsize /= 4 ) THEN + IF ( iacrsize /= 4 ) THEN IF ( iacrsize .eq. 1 ) THEN ratio = 500.e-6/xdia(mgs,lr,1) ELSEIF ( iacrsize .eq. 2 ) THEN @@ -16793,10 +17194,10 @@ subroutine nssl_2mom_gs & qr = (tmp1 + dely*dqiacralphainv*(tmp2 - tmp1))*qx(mgs,lr) - ELSE ! iacrsize == 4 : use all - nr = cx(mgs,lr) - qr = qx(mgs,lr) - ENDIF + ELSE ! iacrsize == 4 : use all + nr = cx(mgs,lr) + qr = qx(mgs,lr) + ENDIF vt = Sqrt((vtxbar(mgs,lr,1)-vtxbar(mgs,li,1))**2 + & & 0.04*vtxbar(mgs,lr,1)*vtxbar(mgs,li,1) ) @@ -16805,9 +17206,9 @@ subroutine nssl_2mom_gs & & ( da0(li)*xdia(mgs,li,3)**2 + & & dab1lh(mgs,li,lr)*xdia(mgs,lh,3)*xdia(mgs,li,3) + & & da1(lr)*xdia(mgs,lr,3)**2 ) - + qiacr(mgs) = Min( qrmxd(mgs), qiacr(mgs) ) - + ciacr(mgs) = 0.25*pi*eri(mgs)*ni*nr*vt* & & ( da0(li)*xdia(mgs,li,3)**2 + & @@ -16895,7 +17296,13 @@ subroutine nssl_2mom_gs & ! ave. diam of freezing drops in microns IF ( qiacr(mgs)*dtp > qxmin(lh) .and. ciacr(mgs) > 1.e-3 ) THEN tmpdiam = 1.e6*( 6.*qiacr(mgs)/(1000.*pi*ciacr(mgs) ) )**(1./3.) ! avg. diameter of newly frozen drops in microns - csplinter(mgs) = lawson_splinter_fac*tmpdiam**4*ciacr(mgs) + fac = 1.0 + IF ( nsplinter .eq. 1001 ) THEN + ! fac = 0.2/sqrt(2.0*pi*10.**2)*Exp(-0.5*((258.-temg(mgs))/10.)**2 ) ! temperature dependence from Sullivan et al. 2018 ACP + ! ELSE + fac = 0.2*Exp(-0.5*((258.-temg(mgs))/10.)**2 ) ! temperature dependence from Sullivan et al. 2018 ACP + ENDIF + csplinter(mgs) = fac*lawson_splinter_fac*tmpdiam**4*ciacr(mgs) ENDIF ELSEIF ( nsplinter .ge. 0 ) THEN csplinter(mgs) = nsplinter*ciacr(mgs) @@ -16961,7 +17368,7 @@ subroutine nssl_2mom_gs & if ( ipconc .ge. 2 .or. ipelec .ge. 9 ) then do mgs = 1,ngscnt ciacw(mgs) = 0.0 - IF ( eiw(mgs) .gt. 0.0 ) THEN + IF ( eiw(mgs) .gt. 0.0 .and. xmas(mgs,lc) > 0.0 ) THEN ciacw(mgs) = qiacw(mgs)*rho0(mgs)/xmas(mgs,lc) ciacw(mgs) = min(ciacw(mgs),ccmxd(mgs)) ENDIF @@ -16972,6 +17379,7 @@ subroutine nssl_2mom_gs & if (ndebug .gt. 0 ) write(0,*) 'ICEZVD_GS: conc 18' if ( ipconc .ge. 2 .or. ipelec .ge. 1 ) then do mgs = 1,ngscnt + tmp1 = 0.0 cracw(mgs) = 0.0 cracr(mgs) = 0.0 ec0(mgs) = 1.e9 @@ -16987,7 +17395,7 @@ subroutine nssl_2mom_gs & & + 2.0*gf2*xdia(mgs,lc,1)*xdia(mgs,lr,1) & & + gf3*xdia(mgs,lr,2) ) ENDIF - ELSE ! IF ( ipconc .ge. 3 .and. + ELSE ! IF ( ipconc .ge. 3 .and. ) IF ( dmrauto <= 0 .or. rho0(mgs)*qx(mgs,lr) > 1.2*xl2p(mgs) ) THEN !{ IF ( 0.5*xdia(mgs,lr,3) .gt. rh(mgs) ) THEN ! { .or. cx(mgs,lr) .gt. nh(mgs) ! IF ( qx(mgs,lc) .gt. qxmin(lc) .and. qx(mgs,lr) .gt. qxmin(lr) ) THEN @@ -17016,7 +17424,7 @@ subroutine nssl_2mom_gs & ! Rain self collection (cracr) and break-up (factor of ec0) ! ! - ec0(mgs) = 2.e9 + ec0(mgs) = 1.0 ! 2.e9 IF ( qx(mgs,lr) .gt. qxmin(lr) ) THEN rwrad = 0.5*xdia(mgs,lr,3) @@ -17032,38 +17440,120 @@ subroutine nssl_2mom_gs & tmp = xdia(mgs,lr,3) - 0.1e-3 ENDIF -! IF ( xdia(mgs,lr,3) .gt. 2.0e-3 .or. icracr <= 0 ) THEN - IF ( tmp .gt. 1.9e-3 .or. icracr <= 0 ) THEN +! Using collection efficiency factor ec0 to simulate break-up that off-sets self-collection (Zieger 1985; Cohard & Pinty 2000) +! ec0 is 1 for rain diameter < 600 microns and then drop off toward zero until diameter of 2mm to represent passive breakup +! ec0 does not go negative here (i.e., does not follow other versions that create extra breakup at large rain diameter) + IF ( ( tmp .gt. 1.9e-3 .and. irainbreak /= 10 .and. irainbreak /= 20 ) .or. icracr <= 0 ) THEN ec0(mgs) = 0.0 cracr(mgs) = 0.0 + IF ( ibincracr == 3 ) THEN + tmp1 = aa1*(cx(mgs,lr)*xv(mgs,lr))**2* & + & (alpha(mgs,lr) + 6.)*(alpha(mgs,lr) + 5.)*(alpha(mgs,lr) + 4.)/ & + & ((alpha(mgs,lr) + 3.)*(alpha(mgs,lr) + 2.)*(alpha(mgs,lr) + 1.)) + ENDIF ELSE IF ( dmrauto <= 0 .or. rho0(mgs)*qx(mgs,lr) > 1.2*xl2p(mgs) ) THEN - IF ( xdia(mgs,lr,3) .lt. 6.1e-4 ) THEN + + IF ( xdia(mgs,lr,3) .lt. 6.1e-4 .or. irainbreak == 10 ) THEN ec0(mgs) = 1.0 ELSE - ec0(mgs) = Exp(-50.0*(50.0*(xdia(mgs,lr,3) - 6.0e-4))) + ec0(mgs) = Exp( -2500.0*(xdia(mgs,lr,3) - 6.0e-4) ) ENDIF + IF ( rwrad .ge. 50.e-6 ) THEN - cracr(mgs) = ec0(mgs)*aa2*cx(mgs,lr)**2*xv(mgs,lr) + tmp1 = aa2*cx(mgs,lr)**2*xv(mgs,lr) + cracr(mgs) = ec0(mgs)*tmp1 + IF ( irainbreak == 20 ) THEN + cracr(mgs) = tmp1 + ENDIF ELSE IF ( imurain == 3 ) THEN cracr(mgs) = ec0(mgs)*aa1*(cx(mgs,lr)*xv(mgs,lr))**2* & & (alpha(mgs,lr) + 2.)/(alpha(mgs,lr) + 1.) ELSE ! imurain == 1 - cracr(mgs) = ec0(mgs)*aa1*(cx(mgs,lr)*xv(mgs,lr))**2* & + tmp1 = aa1*(cx(mgs,lr)*xv(mgs,lr))**2* & & (alpha(mgs,lr) + 6.)*(alpha(mgs,lr) + 5.)*(alpha(mgs,lr) + 4.)/ & & ((alpha(mgs,lr) + 3.)*(alpha(mgs,lr) + 2.)*(alpha(mgs,lr) + 1.)) - + cracr(mgs) = ec0(mgs)*tmp1 + IF ( irainbreak == 20 ) THEN + cracr(mgs) = tmp1 + ENDIF ENDIF - ENDIF + ENDIF ! rwrad > 50 ! cracr(mgs) = Min(cracr(mgs),crmxd(mgs)) - ENDIF - ENDIF - ENDIF + ENDIF ! dmrauto <= 0 + ENDIF ! tmp > 1.9e-3 + + IF ( irainbreak == 100 ) THEN ! Morrison breakup + ec0(mgs) = 1.0 + IF ( xdia(mgs,lr,1) > 300.e-6 ) THEN + ec0(mgs) = 2. - Exp(2300.*(xdia(mgs,lr,1)-300.e-6)) + ENDIF + cracr(mgs) = 5.78*ec0(mgs)*cx(mgs,lr)*qx(mgs,lr) + ENDIF + + ENDIF ! ( qx(mgs,lr) .gt. qxmin(lr) ) + + ! active breakup option + crbreak = 0.0 + IF ( irainbreak == 1 .or. irainbreak == 10 ) THEN + crbreak = Max( 0.0, rainbreakfac* (rho0(mgs)*qx(mgs,lr))**2 ) ! hand fit to lower range of wkqss output + cracr(mgs) = cracr(mgs) - crbreak ! cracr is subtracted, so negative value for breakup + ELSEIF ( irainbreak == 2 .or. irainbreak == 20 ) THEN + ! irainbreak == 20 does not work as intended + crbreak = Max( 0.0, rainbreakfac*(1. - ec0(mgs))*(rho0(mgs)*qx(mgs,lr))**2 ) ! hand fit to lower range of wkqss output +! crbreak = Max(0.0, -0.18 + 1.139e6 * (rho0(mgs)*qx(mgs,lr) + 0.00038106)**2) + cracr(mgs) = cracr(mgs) - crbreak ! cracr is subtracted, so negative value for breakup + ELSEIF ( irainbreak == 11 .and. rho0(mgs)*qx(mgs,lr) > qrbrthresh1 .and. ipconc >= 5 ) THEN + + ! Ad hoc method to break up drops in the DSD tail (D > draintail) + + ratio = Min( maxratiolu, draintail/xdia(mgs,lr,1) ) + ! mass + tmp2 = gaminterp(ratio,alpha(mgs,lr),4,1) + qxd1 = qx(mgs,lr)*(tmp2) + qrbreak = dtpinv*qxd1 + + crbreaksmall = rho0(mgs)*qrbreak/(xdn(mgs,lr)*pi/6.*drsmall**3) + IF ( ( qxd1 > qxmin(lr)) ) THEN + + ! number + tmp = gaminterp(ratio,alpha(mgs,lr),1,1) + IF ( ipconc == 5 ) THEN + ! tmp = Min( 0.2, tmp ) + ENDIF + cxd1 = cx(mgs,lr)*( tmp) + IF ( rho0(mgs)*qx(mgs,lr) > qrbrthresh2 ) THEN + flim = 1.0 + ELSE + flim = (rho0(mgs)*qx(mgs,lr) - qrbrthresh1)/(qrbrthresh2 - qrbrthresh1) + ENDIF + crbreak = flim*(crbreaksmall - dtpinv*cxd1) + +! IF ( kgs(mgs) == 1 .and. qx(mgs,lr) > 0.1e-3 ) THEN +! write(0,*) 'crbreak: ',crbreak,crbreaksmall,dtpinv*cxd1,cx(mgs,lr),cracr(mgs) - crbreak +! ENDIF + cracr(mgs) = cracr(mgs) - crbreak ! cracr is subtracted, so negative value for breakup + + ! reflectivity -- not used yet: goes into zracr +! IF ( ipconc >= 6 .and. lzr > 1 ) THEN +! tmp3 = gaminterp(ratio,alpha(mgs,lr),11,1) +! zxd1 = zx(mgs,lr)*(tmp3) +! zrbreak = dtpinv*zxd1 +! ELSE +! zxd1 = 0 +! ENDIF +! zrbreak = Max(0.0, zrbreak - crbreaksmall*drsmall**6) + ELSEIF ( irainbreak == 12 ) THEN + crbreak = Max( 0.0, 3.8098 * (rho0(mgs)*qx(mgs,lr))**1.9416 ) ! best fit to lower range of wkqss (collision only) output + cracr(mgs) = cracr(mgs) - crbreak ! cracr is subtracted, so negative value for breakup + ENDIF + ENDIF ! cracw(mgs) = min(cracw(mgs),cxmxd(mgs,lc)) + end do end if @@ -17144,24 +17634,6 @@ subroutine nssl_2mom_gs & end if - chacis(:) = 0.0 - if ( lis > 1 .and. ipconc .ge. 5 .or. ipelec .ge. 1 ) then - do mgs = 1,ngscnt - IF ( ehis(mgs) .gt. 0.0 .or. ( ehisclsn(mgs) > 0.0 .and. ipelec > 0 )) THEN - - vt = Sqrt((vtxbar(mgs,lh,1)-vtxbar(mgs,lis,1))**2 + & - & 0.04*vtxbar(mgs,lh,1)*vtxbar(mgs,lis,1) ) - - chacis0(mgs) = 0.25*pi*ehisclsn(mgs)*cx(mgs,lh)*cx(mgs,lis)*vt* & - & ( da0lh(mgs)*xdia(mgs,lh,3)**2 + & - & dab0lh(mgs,lh,lis)*xdia(mgs,lh,3)*xdia(mgs,lis,3) + & - & da0(lis)*xdia(mgs,lis,3)**2 ) - - - chacis(mgs) = min(ehis(mgs)*chacis0(mgs),cxmxd(mgs,lis)) - ENDIF - end do - end if ! ! if (ndebug .gt. 0 ) write(0,*) 'ICEZVD_GS: conc 22nn' @@ -17269,28 +17741,6 @@ subroutine nssl_2mom_gs & end if - IF ( lis > 1 .and. ipconc .ge. 5) THEN - - if (ndebug .gt. 0 ) write(0,*) 'ICEZVD_GS: conc 22kk' - chlacis(:) = 0.0 - chlacis0(:) = 0.0 - do mgs = 1,ngscnt - IF ( lhl .gt. 1 .and. ( ehlis(mgs) .gt. 0.0 .or. (ipelec > 0 .and. ehlisclsn(mgs) > 0.0) ) ) THEN - - vt = Sqrt((vtxbar(mgs,lhl,1)-vtxbar(mgs,lis,1))**2 + & - & 0.04*vtxbar(mgs,lhl,1)*vtxbar(mgs,lis,1) ) - - chlacis0(mgs) = 0.25*pi*ehlisclsn(mgs)*cx(mgs,lhl)*cx(mgs,lis)*vt* & - & ( da0lhl(mgs)*xdia(mgs,lhl,3)**2 + & - & dab0(lhl,lis)*xdia(mgs,lhl,3)*xdia(mgs,lis,3) + & - & da0(lis)*xdia(mgs,lis,3)**2 ) - - - chlacis(mgs) = min(ehlis(mgs)*chlacis0(mgs),cxmxd(mgs,lis)) - ENDIF - end do - ENDIF - ! ! if (ndebug .gt. 0 ) write(0,*) 'ICEZVD_GS: conc 22jj' @@ -17659,16 +18109,45 @@ subroutine nssl_2mom_gs & ELSE !{ - - IF ( ipconc >= 6 .and. lzr > 1 ) THEN + + IF ( ipconc >= 5 .or. lzr > 1 ) THEN + + cxd1 = crfrz(mgs)*dtp + qxd1 = qrfrz(mgs)*dtp + ! interpolate along x, i.e., ratio; tmp1 = ziacrratio(i,j) + delx*dqiacrratioinv*(ziacrratio(ip1,j) - ziacrratio(i,j)) tmp2 = ziacrratio(i,jp1) + delx*dqiacrratioinv*(ziacrratio(ip1,jp1) - ziacrratio(i,jp1)) ! interpolate along alpha; - zrfrz(mgs) = (tmp1 + dely*dqiacralphainv*(tmp2 - tmp1))*zx(mgs,lr)*dtpinv + IF ( ipconc >= 6 .and. lzr > 1 ) THEN + zxd1 = (tmp1 + dely*dqiacralphainv*(tmp2 - tmp1))*zx(mgs,lr) + ! Do the correction for alphamax + zrfrz(mgs) = zxd1*dtpinv + ! tmp4 is the Z from the converted particles assuming shape of alphamax + tmp3 = g1xmax*(rho0(mgs)*qxd1)**2/((pi*xdn(mgs,lh)/6.0)**2) + tmp4 = tmp3/cxd1 + IF ( tmp4 > zxd1 ) THEN ! calculate new graupel/fd number to match zxd1 + ! increase cxd1 to make z,q,c rates consistent + ! cxd1 = g1xmax*(rho0(mgs)*qxd1)**2/(zxd1*(pi*xdn(mgs,lh)/6.0)**2) + cxd1 = tmp3/zxd1 + crfrzf(mgs) = dtpinv*cxd1 + ENDIF + ELSE + ! tmp5 is rain reflectivity moment + tmp5 = g1x(mgs,lr)*(rho0(mgs)*qx(mgs,lr))**2/((pi*xdn(mgs,lr)/6.)**2*cx(mgs,lr)) + zxd1 = (tmp1 + dely*dqiacralphainv*(tmp2 - tmp1))*tmp5 + ! tmp4 is the reflectivity of the newly-converted graupel particles (use g1x(lh) for loss term) + ! which we want to match zxd1 to prevent spurious increase in total reflectivity + tmp3 = g1x(mgs,lr)*(rho0(mgs)*qxd1)**2/((pi*xdn(mgs,lr)/6.0)**2) + tmp4 = tmp3/cxd1 + IF ( tmp4 > zxd1 ) THEN ! calculate new FD number to match zxd1 + crfrzf(mgs) = tmp3/zxd1*dtpinv + ENDIF + ENDIF ENDIF + IF ( ibiggsmallrain > 0 .and. xv(mgs,lr) < 2.*xvmn(lr) .and. ( ibiggsnow == 1 .or. ibiggsnow == 3 ) ) THEN ! IF ( ibiggsmallrain > 0 .and. xv(mgs,lr) < xvbiggsnow .and. ( ibiggsnow == 1 .or. ibiggsnow == 3 ) ) THEN @@ -17706,7 +18185,7 @@ subroutine nssl_2mom_gs & ! j = Int(Max(0.0,Min(15.,alpha(mgs,lr)))*dqiacralphainv) ! j = Int(Max(alphamin,Min(alphamax,alpha(mgs,lr)))*dqiacralphainv) IF ( alp0flag ) THEN - j = Int(Max(0.0,Min(15.,alpha(mgs,lr)))*dqiacralphainv) + j = Int(Max(0.0,Min(alphamax,alpha(mgs,lr)))*dqiacralphainv) ELSE j = Int(Max(minalphalu,Min(maxalphalu,alpha(mgs,lr)))*dqiacralphainv) ENDIF @@ -17892,7 +18371,13 @@ subroutine nssl_2mom_gs & tmp = 0 IF ( qrfrz(mgs)*dtp > qxmin(lh) .and. crfrz(mgs) > 1.e-3 ) THEN tmpdiam = 1.e6*( 6.*qrfrz(mgs)/(1000.*pi*crfrz(mgs) ))**(1./3.) ! avg. diameter of newly frozen drops in microns - tmp = lawson_splinter_fac*tmpdiam**4*crfrz(mgs) + fac = 1.0 + IF ( nsplinter .eq. 1001 ) THEN + ! fac = 0.2/sqrt(2.0*pi*10.**2)*Exp(-0.5*((258.-temg(mgs))/10.)**2 ) ! temperature dependence from Sullivan et al. 2018 ACP + ! ELSE + fac = 0.2*Exp(-0.5*((258.-temg(mgs))/10.)**2 ) ! temperature dependence from Sullivan et al. 2018 ACP + ENDIF + tmp = fac*lawson_splinter_fac*tmpdiam**4*crfrz(mgs) ENDIF ELSEIF ( nsplinter .gt. 0 ) THEN tmp = nsplinter*crfrz(mgs) @@ -18358,6 +18843,7 @@ subroutine nssl_2mom_gs & IF ( ipconc >= 7 ) THEN + ! vent coeff. for reflectivity rate from evaporation alpr = Min(alpharmax,alpha(mgs,lr) ) tmp = alpr + 5.5 + br/2. @@ -19119,11 +19605,11 @@ subroutine nssl_2mom_gs & felscptmp = (fels(mgs)-rw*temg(mgs))/cvm ENDIF - IF ( eqtset > 2 ) THEN - pipert(mgs) = pipert(mgs) + (0 & - & +felspi(mgs)*dqci(mgs) & - & +felvpi(mgs)*dqcw(mgs))*dtp - ENDIF +! IF ( eqtset > 2 ) THEN +! pipert(mgs) = pipert(mgs) + (0 & +! & +felspi(mgs)*dqci(mgs) & +! & +felvpi(mgs)*dqcw(mgs)) ! *dtp +! ENDIF ! ! @@ -19141,7 +19627,13 @@ subroutine nssl_2mom_gs & tqvcon = temgtmp-cbw ltemq = (temgtmp-163.15)/fqsat+1.5 ltemq = Min( nqsat, Max(1,ltemq) ) - qvstmp = pqs(mgs)*tabqvs(ltemq) + + IF ( iqvsopt == 0 ) THEN + qvstmp = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + qvstmp = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF + qisstmp = pqs(mgs)*tabqis(ltemq) qctmp(mgs) = max( 0.0, qctmp(mgs) ) qitmp(mgs) = max( 0.0, qitmp(mgs) ) @@ -19365,8 +19857,8 @@ subroutine nssl_2mom_gs & ! ELSE ! cscnis(mgs) = 0.0 ! ENDIF + ! write(91,*) 'qi,qscni = ',igs(mgs),kgs(mgs),qx(mgs,li),qscni(mgs),cscnis(mgs),qidpv(mgs) ENDIF - IF ( iscni .ne. 4 ) THEN ! crystal aggregation to become snow ! erm 9/5/08 commented second line and added xv to 1st line (zrnic et al 1993) @@ -19407,7 +19899,133 @@ subroutine nssl_2mom_gs & end if end do + IF ( incwet < 1 ) THEN + dhwet(:) = d1t + dhlwet(:) = d1t + dfwet(:) = d1t + ENDIF + IF ( incwet >= 1 ) THEN + ! 'incwet' = incomplete gamma for wet growth + ! Find diameter where wet growth starts, then compute dry and wet growth + ! over [dwet,infinity]. Subtract dry growth from qxacw etc. to get total + ! dry growth part + dhwet(:) = dg0thresh + 0.0001 + dhlwet(:) = dg0thresh + 0.0001 + dfwet(:) = dg0thresh + 0.0001 + + DO mgs = 1,ngscnt + + sqrtrhovt = Sqrt( rhovt(mgs) ) + fventh = sqrtrhovt*(fpndl(mgs)**(1./3.)) * (fakvisc(mgs))**(-0.5) + fventm = sqrtrhovt*(fschm(mgs)**(1./3.)) * (fakvisc(mgs))**(-0.5) + ltemq = (tfr-163.15)/fqsat+1.5 + qvs0 = pqs(mgs)*tabqvs(ltemq) + denomdp = felf(mgs) + fcw(mgs)*temcg(mgs) + denominvdp = 1.d0/(felf(mgs) + fcw(mgs)*temcg(mgs)) + + IF (((qhacw(mgs) + qhacr(mgs))*dtp > qxmin(lh) .and. qx(mgs,lh) > hlcnhqmin .and. & + temg(mgs) .le. tfr + wetgrthtoffset .and. temg(mgs) .ge. 243.15 ) ) THEN +! dw = 0.01*( Exp( -temcg(mgs)/( 1.1e4 * rho0(mgs)*ehw(mgs)*qx(mgs,lc) - 1.3e3*rho0(mgs)*qx(mgs,li) + 1.0 ) ) - 1.0 ) +! dwr = 0.01*( Exp( -temcg(mgs)/( 1.1e4 * rho0(mgs)*(ehw(mgs)*qx(mgs,lc)+ehr(mgs)*qx(mgs,lr)) - & +! 1.3e3*rho0(mgs)*qx(mgs,li) + 1.0 ) ) - 1.0 ) + x = 1.1e4 * rho0(mgs)*(ehw(mgs)*qx(mgs,lc)+ehr(mgs)*qx(mgs,lr)) - & + 1.3e3*rho0(mgs)*qx(mgs,li) + 1.0 + IF ( x > 1.e-20 ) THEN + arg = Min(70.0, (-temcg(mgs)/x )) ! prevent overflow of the exp function in 32 bit + dwr = 0.01*(exp(arg) - 1.0) + ELSE + dwr = 1.e30 + ENDIF + d = dwr + + IF ( dwr < 0.2 .and. dwr > 0.0 .and. rho0(mgs)*(qx(mgs,lc)+qx(mgs,lr)) > 1.e-4 ) THEN + + h1 = ( -ftka(mgs)*temcg(mgs) - felv(mgs)*fwvdf(mgs)*rho0(mgs)*(qx(mgs,lv) - qvs0) ) + h2 = ehi(mgs)*qx(mgs,li)*rho0(mgs)*fci(mgs)*temcg(mgs) + h3 = Max(dwehwmin, ehw(mgs))*qx(mgs,lc) + h4 = ehr(mgs)* qx(mgs,lr) + ! iterate to find minimum diameter for wet growth. Start with value of dwr + DO n = 1,10 + d = Max(d, 1.e-4) + dold = d + vth = axx(mgs,lh)*d**bxx(mgs,lh) + x2 = fventh*sqrtrhovt*Sqrt(d*vth) + IF ( x2 > 1.4 ) THEN + ah = 0.78 + 0.308*x2 ! heat ventillation + ELSE + ah = 1.0 + 0.108*x2**2 ! mass ventillation (Beard and Pruppacher 1971, eq. 9) + ENDIF + + + d = 8.*ah*h1/ & + ( ( Max(0.001,vth - vtxbar(mgs,lc,1))*h3 + & + Max(0.001,vth - vtxbar(mgs,lr,1))*h4) *rho0(mgs)*denomdp + & + Max(0.001,vth - vtxbar(mgs,li,1))*h2) + + IF ( Abs(dold - d)/dold < 0.05 .or. ( n > 3 .and. d > dg0thresh ) ) EXIT + + ENDDO + ENDIF + + dhwet(mgs) = Min(dg0thresh + 0.0001, Max( d, dwetmin )) + ELSE + dhwet(mgs) = dg0thresh + 0.0001 + ENDIF + + IF (((qhlacw(mgs) + qhlacr(mgs))*dtp > qxmin(lhl) .and. qx(mgs,lhl) > 0.01e-3 & + .and. temg(mgs) .le. tfr + wetgrthtoffset .and. temg(mgs) .ge. 243.15 ) ) THEN +! dw = 0.01*( Exp( -temcg(mgs)/( 1.1e4 * rho0(mgs)*ehlw(mgs)*qx(mgs,lc) - 1.3e3*rho0(mgs)*qx(mgs,li) + 1.0 ) ) - 1.0 ) +! dwr = 0.01*( Exp( -temcg(mgs)/( 1.1e4 * rho0(mgs)*(ehlw(mgs)*qx(mgs,lc)+ehlr(mgs)*qx(mgs,lr)) - & +! 1.3e3*rho0(mgs)*qx(mgs,li) + 1.0 ) ) - 1.0 ) + x = 1.1e4 * rho0(mgs)*(ehlw(mgs)*qx(mgs,lc)+ehlr(mgs)*qx(mgs,lr)) - & + 1.3e3*rho0(mgs)*qx(mgs,li) + 1.0 + IF ( x > 1.e-20 ) THEN + arg = Min(70.0, (-temcg(mgs)/x )) ! prevent overflow of the exp function in 32 bit + dwr = 0.01*(exp(arg) - 1.0) + ELSE + dwr = 1.e30 + ENDIF + d = dwr + IF ( dwr < 0.2 .and. dwr > 0.0 .and. rho0(mgs)*(qx(mgs,lc)+qx(mgs,lr)) > 1.e-4 ) THEN + +! write(91,*) 'dw,dwr,temcg = ',100.*dw,100.*dwr,temcg(mgs) + h1 = ( -ftka(mgs)*temcg(mgs) - felv(mgs)*fwvdf(mgs)*rho0(mgs)*(qx(mgs,lv) - qvs0) ) + h2 = ehi(mgs)*qx(mgs,li)*rho0(mgs)*fci(mgs)*temcg(mgs) + h3 = Max(dwehwmin, ehlw(mgs))*qx(mgs,lc) + h4 = ehlr(mgs)* qx(mgs,lr) + ! iterate to find minimum diameter for wet growth. Start with value of dwr + DO n = 1,10 + d = Max(d, 1.e-4) + dold = d + vth = axx(mgs,lhl)*d**bxx(mgs,lhl) + x2 = fventh*sqrtrhovt*Sqrt(d*vth) + IF ( x2 > 1.4 ) THEN + ah = 0.78 + 0.308*x2 ! heat ventillation + ELSE + ah = 1.0 + 0.108*x2**2 ! mass ventillation (Beard and Pruppacher 1971, eq. 9) + ENDIF + + + d = 8.*ah*h1/ & + ( ( Max(0.001,vth - vtxbar(mgs,lc,1))*h3 + & + Max(0.001,vth - vtxbar(mgs,lr,1))*h4) *rho0(mgs)*denomdp + & + Max(0.001,vth - vtxbar(mgs,li,1))*h2) + + IF ( Abs(dold - d)/dold < 0.05 .or. ( n > 3 .and. d > dg0thresh ) ) EXIT + + ENDDO + ENDIF + + dhlwet(mgs) = Min(dg0thresh + 0.0001, Max( d, dwetmin ) ) + ELSE + dhlwet(mgs) = dg0thresh + 0.0001 + ENDIF + + + ENDDO + + ENDIF ! incwet @@ -19448,12 +20066,88 @@ subroutine nssl_2mom_gs & ! IF ( dnu(lh) .ne. 0. ) THEN ! qhwet(mgs) = qhdry(mgs) ! ELSE - IF ( incwet == 0 ) THEN + ! IF ( incwet == 0 ) THEN qhwet(mgs) = & & ( xdia(mgs,lh,1)*hwvent(mgs)*cx(mgs,lh)*fwet1(mgs) & & + fwet2(mgs)*(qhaci(mgs) + qhacs(mgs)) ) - qhwet(mgs) = max( 0.0, qhwet(mgs)) - ELSE + qhwet(mgs) = max( 0.0, qhwet(mgs)) + + IF ( incwet == 1 .and. qhwet(mgs) < qhdry(mgs) .and. dhwet(mgs) < dg0thresh ) THEN + ! ELSE + ! IF ( dhwet(mgs) < dg0thresh ) THEN + ! find portion of qc and qr collection that are dry/wet growth for d > dwet + + ratio = Min( maxratiolu, dhwet(mgs)/xdia(mgs,lh,1) ) + + tmp1 = gaminterp(ratio,alpha(mgs,lh),13,1) ! alpha + 3 + tmp2 = gaminterp(ratio,alpha(mgs,lh),12,1) ! alpha + 2 + tmp3 = gaminterp(ratio,alpha(mgs,lh), 9,1) ! alpha + 1 + + IF ( qhacw(mgs)*dtp > qxmin(lh) ) THEN + vt = abs(vtxbar(mgs,lh,1)-vtxbar(mgs,lc,1)) + + qxacwtmp = 0.25*pi*ehw(mgs)*cx(mgs,lh)*(qx(mgs,lc)-qcwresv(mgs))*vt* & + & ( tmp1*da0lh(mgs)*xdia(mgs,lh,3)**2 + & + & tmp2*dab1lh(mgs,lh,lc)*xdia(mgs,lh,3)*xdia(mgs,lc,3) + & + & tmp3*da1lc(mgs)*xdia(mgs,lc,3)**2 ) + ELSE + qxacwtmp = 0.0 + ENDIF + + IF ( qhacr(mgs)*dtp > qxmin(lh) ) THEN + + vt = Sqrt((vtxbar(mgs,lh,1)-vtxbar(mgs,lr,1))**2 + & + & 0.04*vtxbar(mgs,lh,1)*vtxbar(mgs,lr,1) ) + + qxacrtmp = 0.25*pi*ehr(mgs)*cx(mgs,lh)*qx(mgs,lr)*vt* & + & ( tmp1*da0lh(mgs)*xdia(mgs,lh,3)**2 + & + & tmp2*dab1lh(mgs,lh,lr)*xdia(mgs,lh,3)*xdia(mgs,lr,3) + & + & tmp3*da1lr(mgs)*xdia(mgs,lr,3)**2 ) + ELSE + qxacrtmp = 0.0 + ENDIF + + ! hwvent is where the size dependency is, so hxventtmp gives the portion for d > dwet + x = gaminterp(ratio,alpha(mgs,lh),9,1) ! alpha + 1 + y = gaminterp(ratio,alpha(mgs,lh),3,1) ! alpha + b/2 + 5/2 + + hxventtmp = 0.78*x + y*hwventy(mgs) ! & + + ! find the ice and snow collection for d > dwet + qxacitmp = 0.0 + IF ( qhaci(mgs)*dtp > qxmin(lh) ) THEN + vt = abs(vtxbar(mgs,lh,1)-vtxbar(mgs,li,1)) + + qxacitmp = 0.25*pi*ehiclsn(mgs)*cx(mgs,lh)*qx(mgs,li)*vt* & + & ( tmp1*da0lh(mgs)*xdia(mgs,lh,3)**2 + & + & tmp2*dab1lh(mgs,lh,li)*xdia(mgs,lh,3)*xdia(mgs,li,3) + & + & tmp3*da1(li)*xdia(mgs,li,3)**2 ) + ENDIF + + qxacstmp = 0.0 + IF ( qhacs(mgs)*dtp > qxmin(lh) ) THEN + vt = abs(vtxbar(mgs,lh,1)-vtxbar(mgs,ls,1)) + + qxacstmp = 0.25*pi*ehsclsn(mgs)*cx(mgs,lh)*qx(mgs,ls)*vt* & + & ( da0lh(mgs)*xdia(mgs,lh,3)**2 + & + & dab1lh(mgs,lh,ls)*xdia(mgs,lh,3)*xdia(mgs,ls,3) + & + & da1(ls)*xdia(mgs,ls,3)**2 ) + ENDIF + + qxwettmp = & + & xdia(mgs,lh,1)*hxventtmp*cx(mgs,lh)*fwet1(mgs) & + & + fwet2(mgs)*(qxacitmp + qxacstmp) + + ! as dry growth but subtract part for D > Dw and add wet growth for D > Dw + qhwet(mgs) = qhacw(mgs) + qhacr(mgs) + qhaci(mgs) + qhacs(mgs) & + - ehi(mgs)*qxacitmp - ehs(mgs)*qxacstmp & + - qxacwtmp - qxacrtmp + qxwettmp + + ! qhacw(mgs) = Min( qhacw(mgs), 0.5*qx(mgs,lc)*dtpinv ) + + ! ELSE ! for dwet > 15cm, just assume dry growth + ! qhwet(mgs) = qhdry(mgs) + ! ENDIF ENDIF ! ENDIF @@ -19461,13 +20155,88 @@ subroutine nssl_2mom_gs & qhlwet(mgs) = 0.0 IF ( lhl .gt. 1 ) THEN - IF ( incwet == 0 ) THEN + !IF ( incwet == 0 ) THEN qhlwet(mgs) = & & ( xdia(mgs,lhl,1)*hlvent(mgs)*cx(mgs,lhl)*fwet1(mgs) & & + fwet2(mgs)*(qhlaci(mgs) + qhlacs(mgs)) ) qhlwet(mgs) = max( 0.0, qhlwet(mgs)) - ELSE + IF ( incwet == 1 .and. qhlwet(mgs) < qhldry(mgs) .and. dhlwet(mgs) < dg0thresh ) THEN + !ELSE +!! || defined (WRFEXTRAS) + ! IF ( dhlwet(mgs) < dg0thresh ) THEN + ! find portion of qc and qr collection that are dry/wet growth for d > dwet + + ratio = Min( maxratiolu, dhlwet(mgs)/xdia(mgs,lhl,1) ) + + tmp1 = gaminterp(ratio,alpha(mgs,lhl),13,2) ! alpha + 3 + tmp2 = gaminterp(ratio,alpha(mgs,lhl),12,2) ! alpha + 2 + tmp3 = gaminterp(ratio,alpha(mgs,lhl), 9,2) ! alpha + 1 + + IF ( qhlacw(mgs)*dtp > qxmin(lhl) ) THEN + vt = abs(vtxbar(mgs,lhl,1)-vtxbar(mgs,lc,1)) + + qxacwtmp = 0.25*pi*ehlw(mgs)*cx(mgs,lhl)*(qx(mgs,lc)-qcwresv(mgs))*vt* & + & ( tmp1*da0lhl(mgs)*xdia(mgs,lhl,3)**2 + & + & tmp2*dab1lh(mgs,lhl,lc)*xdia(mgs,lhl,3)*xdia(mgs,lc,3) + & + & tmp3*da1lc(mgs)*xdia(mgs,lc,3)**2 ) + ELSE + qxacwtmp = 0.0 + ENDIF + + IF ( qhlacr(mgs)*dtp > qxmin(lhl) ) THEN + + vt = Sqrt((vtxbar(mgs,lhl,1)-vtxbar(mgs,lr,1))**2 + & + & 0.04*vtxbar(mgs,lhl,1)*vtxbar(mgs,lr,1) ) + + qxacrtmp = 0.25*pi*ehlr(mgs)*cx(mgs,lhl)*qx(mgs,lr)*vt* & + & ( tmp1*da0lhl(mgs)*xdia(mgs,lhl,3)**2 + & + & tmp2*dab1lh(mgs,lhl,lr)*xdia(mgs,lhl,3)*xdia(mgs,lr,3) + & + & tmp3*da1lr(mgs)*xdia(mgs,lr,3)**2 ) + ELSE + qxacrtmp = 0.0 + ENDIF + + x = gaminterp(ratio,alpha(mgs,lhl),9,2) ! alpha + 1 + y = gaminterp(ratio,alpha(mgs,lhl),3,2) ! alpha + b/2 + 5/2 + + hxventtmp = 0.78*x + y*hlventy(mgs) ! & + + qxacitmp = 0.0 + IF ( qhlaci(mgs)*dtp > qxmin(lhl) ) THEN + vt = abs(vtxbar(mgs,lhl,1)-vtxbar(mgs,li,1)) + + qxacitmp = 0.25*pi*ehliclsn(mgs)*cx(mgs,lhl)*qx(mgs,li)*vt* & + & ( tmp1*da0lhl(mgs)*xdia(mgs,lhl,3)**2 + & + & tmp2*dab1lh(mgs,lhl,li)*xdia(mgs,lhl,3)*xdia(mgs,li,3) + & + & tmp3*da1(li)*xdia(mgs,li,3)**2 ) + ENDIF + + qxacstmp = 0.0 + IF ( qhlacs(mgs)*dtp > qxmin(lhl) ) THEN + vt = abs(vtxbar(mgs,lhl,1)-vtxbar(mgs,ls,1)) + + qxacstmp = 0.25*pi*ehlsclsn(mgs)*cx(mgs,lhl)*qx(mgs,ls)*vt* & + & ( da0lhl(mgs)*xdia(mgs,lhl,3)**2 + & + & dab1lh(mgs,lhl,ls)*xdia(mgs,lhl,3)*xdia(mgs,ls,3) + & + & da1(ls)*xdia(mgs,ls,3)**2 ) + ENDIF + + qxwettmp = & + & xdia(mgs,lhl,1)*hxventtmp*cx(mgs,lhl)*fwet1(mgs) & + & + fwet2(mgs)*(qxacitmp + qxacstmp) + + ! qhlacw(mgs) + qhlacr(mgs) - qxacwtmp - qxacrtmp is the 'dry' growth + ! at smaller diameters +! qhlwet(mgs) = qhlacw(mgs) + qhlacr(mgs) - qxacwtmp - qxacrtmp + qxwettmp + ! as dry growth but subtract part for D > Dw and add wet growth for D > Dw + qhlwet(mgs) = qhlacw(mgs) + qhlacr(mgs) + qhlaci(mgs) + qhlacs(mgs) & + - ehli(mgs)*qxacitmp - ehls(mgs)*qxacstmp & + - qxacwtmp - qxacrtmp + qxwettmp + + ! ELSE + ! qhlwet(mgs) = qhldry(mgs) + ! ENDIF ENDIF ! incwet ENDIF @@ -19876,14 +20645,20 @@ subroutine nssl_2mom_gs & ltest = xdia(mgs,lh,1)*(4. + alpha(mgs,lh)) > Abs( hlcnhdia ) ! test on mass-weighted diameter ENDIF + ! if incwet > 0, then should use dhwet here to avoid calculating again IF ( iusedw == 0 .and. ihlcnh == 1 ) THEN dg0(mgs) = -1. ELSE - IF (((qhacw(mgs) + qhacr(mgs))*dtp > qxmin(lh) .and. qx(mgs,lh) > hlcnhqmin .and. temg(mgs) .le. tfr-2.0 & - .and. temg(mgs) .gt. dwtempmin ) .or. ( wetgrowth(mgs) .and. qx(mgs,lh) > hlcnhqmin ) ) THEN + IF ( temg(mgs) .le. tfr+hailcnvtoffset .and. & + (( (qhacw(mgs) + qhacr(mgs))*dtp > qxmin(lh) .and. qx(mgs,lh) > hlcnhqmin & + .and. temg(mgs) .gt. dwtempmin ) .or. ( wetgrowth(mgs) .and. qx(mgs,lh) > hlcnhqmin )) ) THEN ! dw = 0.01*( Exp( -temcg(mgs)/( 1.1e4 * rho0(mgs)*ehw(mgs)*qx(mgs,lc) - 1.3e3*rho0(mgs)*qx(mgs,li) + 1.0 ) ) - 1.0 ) ! dwr = 0.01*( Exp( -temcg(mgs)/( 1.1e4 * rho0(mgs)*(ehw(mgs)*qx(mgs,lc)+ehr(mgs)*qx(mgs,lr)) - & ! 1.3e3*rho0(mgs)*qx(mgs,li) + 1.0 ) ) - 1.0 ) + IF ( incwet > 0 ) THEN + d = dhwet(mgs) + ELSE + ! First guess for dwet (not that good, but it is something) x = 1.1e4 * rho0(mgs)*(ehw(mgs)*qx(mgs,lc)+ehr(mgs)*qx(mgs,lr)) - & 1.3e3*rho0(mgs)*qx(mgs,li) + 1.0 IF ( x > 1.e-20 ) THEN @@ -19892,7 +20667,7 @@ subroutine nssl_2mom_gs & ELSE dwr = 1.e30 ENDIF - d = dwr + d = Min(dwr, dg0thresh + 0.0001) IF ( dwr < 0.2 .and. dwr > 0.0 .and. rho0(mgs)*(qx(mgs,lc)+qx(mgs,lr)) > 1.e-4 ) THEN sqrtrhovt = Sqrt( rhovt(mgs) ) fventh = sqrtrhovt*(fpndl(mgs)**(1./3.)) * (fakvisc(mgs))**(-0.5) @@ -19945,21 +20720,26 @@ subroutine nssl_2mom_gs & IF ( Abs(dold - d)/dold < 0.05 .or. ( n > 3 .and. d > dg0thresh ) ) EXIT ENDDO - ENDIF + + d = Min( d, dg0thresh + 0.0001 ) + ENDIF ! dwr < 0.2 .and. dwr > 0.0 + ENDIF ! incwet - dg0(mgs) = Min( dwmax, Max( d, dwmin ) ) + ! dg0(mgs) = Min( dwmax, Max( d, dwmin ) ) + dg0(mgs) = Max( d, dwmin ) ELSE - IF ( qx(mgs,lh) > qxmin(lh) .and. qx(mgs,lh) > hlcnhqmin .and. temg(mgs) .le. tfr-2.0 ) THEN - dg0(mgs) = dwmax - ELSE + ! IF ( qx(mgs,lh) > qxmin(lh) .and. qx(mgs,lh) > hlcnhqmin .and. temg(mgs) .le. tfr+hailcnvtoffset ) THEN + ! dg0(mgs) = dwmax + ! ELSE dg0(mgs) = dg0thresh + 0.0001 - ENDIF + ! ENDIF ENDIF IF ( ihlcnh == 3 .and. (qhacw(mgs) + qhacr(mgs))*dtp > qxmin(lh) .and. qx(mgs,lh) > hlcnhqmin & - .and. temg(mgs) .le. tfr-2.0 ) THEN + .and. temg(mgs) .le. tfr+hailcnvtoffset .and. temg(mgs) > 238.0 ) THEN ! set a secondary condition on to capture large graupel that is riming but not in wet growth - dg0(mgs) = Min( dg0(mgs), dg0thresh - 0.0001 ) +! dg0(mgs) = Min( dg0(mgs), dg0thresh - 0.0001 ) + dg0(mgs) = Min( dg0(mgs), dwmax ) ENDIF ENDIF @@ -19973,7 +20753,7 @@ subroutine nssl_2mom_gs & & ltest .and. qx(mgs,lh) .gt. hlcnhqmin ) .or. wtest ) THEN ! { ! : xdia(mgs,lh,3) .gt. 2.e-3 .and. qx(mgs,lh) .gt. 1.0e-3 THEN ! 0823.2008 erm test ! IF ( xdia(mgs,lh,3) .gt. 1.e-3 ) THEN - IF ( qhacw(mgs) .gt. 0.0 .and. qhacw(mgs) .gt. qhaci(mgs) .and. temg(mgs) .le. tfr-2.0 ) THEN ! { + IF ( qhacw(mgs) .gt. 0.0 .and. qhacw(mgs) .gt. qhaci(mgs) .and. temg(mgs) .le. tfr+hailcnvtoffset ) THEN ! { ! dh0 is the diameter dividing wet growth from dry growth (Ziegler 1985), modified by MY05 ! dh0 = 0.01*(exp(temcg(mgs)/(1.1e4*(qx(mgs,lc)+qx(mgs,lr)) - ! : 1.3e3*qx(mgs,li) + 1.0e-3 ) ) - 1.0) @@ -19987,6 +20767,7 @@ subroutine nssl_2mom_gs & ELSE dh0 = 1.e30 ENDIF + dg0(mgs) = Min(dh0, dg0thresh + 0.0001) ENDIF ! wtest ! dh0 = Max( dh0, 5.e-3 ) @@ -20021,7 +20802,7 @@ subroutine nssl_2mom_gs & IF ( wtest .and. & - ( qhacw(mgs)*dtp > qxmin(lh) .and. temg(mgs) .lt. tfr-2. .and. qx(mgs,lh) > hlcnhqmin ) ) THEN + ( qhacw(mgs)*dtp > qxmin(lh) .and. temg(mgs) .lt. tfr+hailcnvtoffset .and. qx(mgs,lh) > hlcnhqmin ) ) THEN ! convert number, mass, and reflectivity for d > dw IF ( ipconc == 5 ) THEN ! dg0(mgs) = Min( dg0(mgs), hldia1 ) @@ -20075,10 +20856,42 @@ subroutine nssl_2mom_gs & tmp3 = gaminterp(ratio,alpha(mgs,lh),11,1) zxd1 = flim*zx(mgs,lh)*(tmp3) zhlcnh(mgs) = dtpinv*zxd1 + + ! tmp4 is the Z from the converted particles assuming shape of alphamax + tmp3 = g1xmax*(rho0(mgs)*qxd1)**2/((pi*xdn(mgs,lh)/6.0)**2) + tmp4 = tmp3/cxd1 + IF ( tmp4 > zxd1 ) THEN ! calculate new hail number to match zxd1 + ! increase cxd1 to make z,q,c rates consistent + ! cxd1 = g1xmax*(rho0(mgs)*qxd1)**2/(zxd1*(pi*xdn(mgs,lh)/6.0)**2) + cxd1 = tmp3/zxd1 + chlcnhhl(mgs) = dtpinv*cxd1 + ENDIF ELSE zxd1 = 0 ENDIF + IF ( ipconc == 5 ) THEN ! Adjust cxd1 by reflectivity removed from graupel + tmp3 = gaminterp(ratio,alpha(mgs,lh),11,1) + ! tmp5 is graupel reflectivity moment + tmp5 = g1x(mgs,lh)*(rho0(mgs)*qx(mgs,lh))**2/((pi*xdn(mgs,lh)/6.)**2*cx(mgs,lh)) + zxd1 = flim*(tmp3)*tmp5 + ! tmp4 is the reflectivity of the newly-converted graupel particles (use g1x(lh) for loss term) + ! which we want to match zxd1 to prevent spurious increase in total reflectivity + tmp3 = g1x(mgs,lh)*(rho0(mgs)*qxd1)**2/((pi*xdn(mgs,lh)/6.0)**2) + tmp4 = tmp3/cxd1 + IF ( tmp4 > zxd1 ) THEN ! calculate new hail number to match zxd1 + ! cxd1 = g1x(mgs,lhl)*(rho0(mgs)*qxd1)**2/(zxd1*pi*xdn(mgs,lh)/6.0) ! trial form results in tiny hail + ! want the adjust size of the new hail so that Z is conserved, so increase number of + ! particles to make qxd1,zxd1, and C consistent. + ! want zxd1 = g1x(mgs,lh)*(rho0(mgs)*qxd1)**2/(c*(pi*xdn(mgs,lh)/6.0)**2) + ! Use g1x(mgs,lh) here instead of g1x(mgs,lhl) because rzxhlh will then multiply + ! by g1x(mgs,lhl)/g1x(mgs,lh) + ! cxd1 = g1x(mgs,lh)*(rho0(mgs)*qxd1)**2/(zxd1*(pi*xdn(mgs,lh)/6.0)**2) + cxd1 = tmp3/zxd1 + chlcnhhl(mgs) = dtpinv*cxd1 ! multiplied later by rzxhlh(mgs) + ENDIF + ENDIF + ELSE qhlcnh(mgs) = 0.0 ENDIF @@ -20104,7 +20917,7 @@ subroutine nssl_2mom_gs & ! convert number, mass, and reflectivity for d > hldia1, ! regardless of wet growth status, but as long as riming > 0 DO mgs = 1,ngscnt - IF ( qhacw(mgs)*dtp > qxmin(lh) .and. temg(mgs) .lt. tfr-2. .and. qx(mgs,lh) > qxmin(lh) ) THEN + IF ( qhacw(mgs)*dtp > qxmin(lh) .and. temg(mgs) .lt. tfr+hailcnvtoffset .and. qx(mgs,lh) > qxmin(lh) ) THEN ratio = Min( maxratiolu, hldia1/xdia(mgs,lh,1) ) ! number @@ -20872,6 +21685,7 @@ subroutine nssl_2mom_gs & ! ! ! + if (ndebug .gt. 0 ) write(0,*) 'dbg = 8a' ! ! @@ -20907,6 +21721,7 @@ subroutine nssl_2mom_gs & ! Cloud ice ! ! IF ( ipconc .ge. 1 ) THEN + if (ndebug .gt. 0 ) write(0,*) 'cloud ice sum' IF ( warmonly < 0.5 ) THEN IF ( ffrzs < 1.0 ) THEN @@ -20952,7 +21767,7 @@ subroutine nssl_2mom_gs & & +chlmul1(mgs) & & + csplinter(mgs) + csplinter2(mgs) & & +csmul(mgs) - + pccii(mgs) = pccii(mgs)*(1. - ffrzs) pccid(mgs) = & ! & il5(mgs)*(-cscni(mgs) - cscnvi(mgs) & ! - cwaci(mgs) & @@ -21073,6 +21888,7 @@ subroutine nssl_2mom_gs & pccii(mgs) = pccii(mgs) - (1.-frac)*il5(mgs)*(cwfrzc(mgs)+cwctfzc(mgs))*(1. - ffrzs) IF ( lhl .gt. 1 ) chlacw(mgs) = frac*chlacw(mgs) + ! STOP ENDIF @@ -21097,16 +21913,16 @@ subroutine nssl_2mom_gs & ! & -csmlr(mgs)/rzxs(mgs) & & -csmlrr(mgs) & & - cimlr(mgs) ) & + & - Min(0.0,cracr(mgs)) & ! cracr is negative if there is enough breakup & -crshr(mgs) !null at this point when wet snow/graupel included pcrwd(mgs) = & & il5(mgs)*(-ciacr(mgs) - crfrz(mgs) ) & ! - cipacr(mgs)) ! > -csacr(mgs) & & - chacr(mgs) - chlacr(mgs) & & +crcev(mgs) & - & - cracr(mgs) + & - Max(0.0,cracr(mgs)) ! > -il5(mgs)*ciracr(mgs) - ELSEIF ( warmonly < 0.8 ) THEN pcrwi(mgs) = & & crcnw(mgs) & @@ -21982,14 +22798,6 @@ subroutine nssl_2mom_gs & zhacr(mgs) = g1*(6.*rho0(mgs)/(pi*xdn(mgs,lh)))**2*( 2.*( qx(mgs,lh)/cx(mgs,lh)) * qhacr(mgs) ) ! zhacrf(mgs) = g1*zhacr - -! z = g1*(6.*rho0(mgs)/(pi*1000.))**2*( (qx(mgs,lh)+dtp*qhacr(mgs))**2)/(cx(mgs,lh)) - - IF ( z > zx(mgs,lh) ) THEN -! zhacr(mgs) = (z - zx(mgs,lh))*dtpinv - ELSE -! zhacr(mgs) = 0.0 - ENDIF ENDIF ! zhacr(mgs) = g1*(6.*rho0(mgs)/(pi*1000.))**2*( 2.*( tmp ) * qhacr(mgs) ) @@ -22000,12 +22808,7 @@ subroutine nssl_2mom_gs & ! : ((3.0 + alp)*(2.0 + alp)*(1.0 + alp)) IF ( qhacw(mgs) .gt. 0.0 ) THEN ! zhacw(mgs) = g1*(6.*rho0(mgs)/(pi*1000.))**2*( 2.*( qx(mgs,lh)/cx(mgs,lh)) * qhacw(mgs) ) - zhacw(mgs) = g1*(6.*rho0(mgs)/(pi*xdn(mgs,lh)))**2*( 2.*( qx(mgs,lh)/cx(mgs,lh)) * qhacw(mgs) ) - -! z = g1*(6.*rho0(mgs)/(pi*1000.))**2*( (qx(mgs,lh)+dtp*(qhacw(mgs)-qhmul1(mgs)))**2)/(cx(mgs,lh)) - IF ( z > zx(mgs,lh) ) THEN -! zhacw(mgs) = (z - zx(mgs,lh))*dtpinv - ENDIF + zhacw(mgs) = g1*(6.*rho0(mgs)/(pi*xdn(mgs,lh)))**2*( 2.*( qx(mgs,lh)/cx(mgs,lh)) * qhacw(mgs) ) ENDIF ELSE ! } { ! this is not used because of the 'true' above @@ -22081,18 +22884,18 @@ subroutine nssl_2mom_gs & r = rho0(mgs)*qhcns(mgs)/vhcns(mgs) ! density of new graupel particles IF ( imusnow == 3 ) THEN zhcns(mgs) = 3.6476*rho0(mgs)**2*(alpha(mgs,ls)+2.)/(r**2*(alpha(mgs,ls)+1.)) * & - & ( 2.*tmp * qhcns(mgs) - tmp**2 * chcns(mgs) ) + & ( 2.*tmp * qhcns(mgs) - tmp**2 * chcnsh(mgs) ) ELSE write(0,*) 'Value of imusnow not valid. Must be 3 (fix me for =1). imusnow = ',imusnow STOP ENDIF ENDIF - IF ( qhcni(mgs) > 0.0 .and. chcni(mgs) > 0.0 .and. cx(mgs,li) > cxmin .and. vhcni(mgs) > 0 ) THEN + IF ( qhcni(mgs) > 0.0 .and. chcnih(mgs) > 0.0 .and. cx(mgs,li) > cxmin .and. vhcni(mgs) > 0 ) THEN tmp = qx(mgs,li)/cx(mgs,li) r = rho0(mgs)*qhcni(mgs)/vhcni(mgs) ! density of new graupel particles zhcni(mgs) = 3.6476*rho0(mgs)**2*(alpha(mgs,li)+2.)/(r**2*(alpha(mgs,li)+1.)) * & - & ( 2.*tmp * qhcni(mgs) - tmp**2 * chcni(mgs) ) + & ( 2.*tmp * qhcni(mgs) - tmp**2 * chcnih(mgs) ) ENDIF @@ -22115,14 +22918,6 @@ subroutine nssl_2mom_gs & & - il5(mgs)*zhlcnh(mgs) - IF ( igs(mgs) == 44 .and. kgs(mgs) == 23 .or. dtp*( pqhwi(mgs) + pqhwd(mgs) ) > qxmin(lh) ) THEN -! write(0,*) 'i,k,time = ',igs(mgs),kgs(mgs),time_real -! write(0,*) 'pzhwi,d = ',pzhwi(mgs),pzhwd(mgs),dtp*( pzhwi(mgs) + pzhwd(mgs) ),zx(mgs,lh) -! write(0,*) 'pqhwi,d = ',pqhwi(mgs),pqhwd(mgs),dtp*( pqhwi(mgs) + pqhwd(mgs) ),qx(mgs,lh) -! write(0,*) 'pchwi,d = ',pchwi(mgs),pchwd(mgs),dtp*( pchwi(mgs) + pchwd(mgs) ),cx(mgs,lh) - ENDIF - - ! IF ( zhcnhl(mgs) < 0.0 ) THEN ! write(0,*) 'Problem with zhcnhl! zhcnhl,qhcnhl,chcnhl = ',zhcnhl(mgs),qhcnhl(mgs),chcnhl(mgs) ! write(0,*) 'g1,tmp = ',g1x(mgs,lhl),tmp @@ -22384,7 +23179,7 @@ subroutine nssl_2mom_gs & zracw(mgs) = g1x(mgs,lr)*(6.*rho0(mgs)/(pi*1000.))**2*( 2.*tmp * qracw(mgs) ) ENDIF - IF ( cracr(mgs) > 0.0 .and. cx(mgs,lr) > 0.0 ) THEN + IF ( cracr(mgs) /= 0.0 .and. cx(mgs,lr) > 0.0 ) THEN zracr(mgs) = g1x(mgs,lr)*(6.*rho0(mgs)/(pi*1000.))**2*( tmp**2 * cracr(mgs) ) ENDIF @@ -22399,7 +23194,7 @@ subroutine nssl_2mom_gs & IF ( iferwisventr == 2 ) THEN vent1 = Min(0.0, (12./(pii*xdn(mgs,lr)))*xdia(mgs,lr,1)**3*fvce(mgs)*rwcap(mgs)*rwventz(mgs)) - zrcev(mgs) = Max( zrcev(mgs), vent1 ) + zrcev(mgs) = Max( dble(zrcev(mgs)), vent1 ) ENDIF ! IF ( ny == 2 .and. igs(mgs) == 20 ) THEN ! write(0,*) 'k,zrcevold,new,maxdep : ',kgs(mgs),zrcev(mgs),vent1,-zxmxd(mgs,lr),alpha(mgs,lr),cx(mgs,lr) @@ -22528,7 +23323,8 @@ subroutine nssl_2mom_gs & ! > + rho0(mgs)*qhshr(mgs)/xdn(mgs,lh) !xdn(mgs,lr) ! ENDIF - IF ( lzh > 1 .and. qx(mgs,lh) > qxmin(lh) ) THEN + IF ( lzh > 1 .and. qx(mgs,lh) > qxmin(lh) .and. & + vx(mgs,lh) + dtp*(pvhwi(mgs) + pvhwd(mgs)) > rho0(mgs)*qxmin(lh)/900. ) THEN ! Calculate change in reflectivity due to density changes xdn_new = rho0(mgs)*(qx(mgs,lh) + dtp*(pqhwi(mgs) + pqhwd(mgs) ))/ & @@ -22637,7 +23433,8 @@ subroutine nssl_2mom_gs & & + rho0(mgs)*(1-il5(mgs))*vhlmlr(mgs)/xdn(mgs,lhl) & & + vhlshdr(mgs) - vhlsoak(mgs) - IF ( lzhl > 1 .and. qx(mgs,lhl) > qxmin(lhl) ) THEN + IF ( lzhl > 1 .and. qx(mgs,lhl) > qxmin(lhl) .and. & + vx(mgs,lhl) + dtp*(pvhli(mgs) + pvhld(mgs)) > rho0(mgs)*qxmin(lhl)/900. ) THEN ! Calculate change in reflectivity due to density changes xdn_new = rho0(mgs)*(qx(mgs,lhl) + dtp*(pqhli(mgs) + pqhld(mgs) ))/ & @@ -22874,12 +23671,17 @@ subroutine nssl_2mom_gs & ! write(iunit,*) il5(mgs)*qscni(mgs), qscnvi(mgs) write(iunit,*) il5(mgs)*qsaci(mgs) - write(iunit,*) il5(mgs)*qrfrzs(mgs) + write(iunit,*) il5(mgs)*qrfrzs(mgs), qiacrs(mgs) write(iunit,*) il5(mgs)*qiacrs(mgs),il3(mgs)*(qiacrf(mgs)+qracif(mgs)),il3(mgs),qiacrf(mgs),qracif(mgs) write(iunit,*) il5(mgs)*qsdpv(mgs), qscev(mgs) - write(iunit,*) qsacw(mgs) + write(iunit,*) qsacw(mgs),qwfrzc(mgs), qwctfzc(mgs), qicichr(mgs) write(iunit,*) qsacr(mgs), qscnh(mgs) - write(iunit,*) 'pqswi = ',pqswi(mgs) + write(iunit,*) il2(mgs)*qsacr(mgs) + write(iunit,*) il5(mgs)*qicicnt(mgs)*ffrzs + write(iunit,*) il3(mgs)*(qiacrf(mgs)+qracif(mgs)) ! only applies for ipconc <= 3 + write(iunit,*) Max(0.0, qscev(mgs)) + write(iunit,*) qsacw(mgs) + qscnh(mgs) + write(iunit,*) 'pqswi = ',pqswi(mgs) write(iunit,*) -qhcns(mgs) write(iunit,*) -qracs(mgs) write(iunit,*) -qhacs(mgs) @@ -23054,6 +23856,7 @@ subroutine nssl_2mom_gs & qwvp(mgs) = qwvp(mgs) + & & dtp*(pqwvi(mgs)+pqwvd(mgs)) + ! qcwresv(mgs) = qx(mgs,lc) ! temporary save of old qc value qx(mgs,lc) = qx(mgs,lc) + & & dtp*(pqcwi(mgs)+pqcwd(mgs)) qx(mgs,lr) = qx(mgs,lr) + & @@ -23358,7 +24161,11 @@ subroutine nssl_2mom_gs & ltemq = (temg(mgs)-163.15)/fqsat+1.5 ltemq = Min( nqsat, Max(1,ltemq) ) - qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + qvs(mgs) = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF IF ( ( qvap(mgs) > qvs(mgs) .or. qx(mgs,lc) > qxmin(lc) ) .and. temg(mgs) > tfrh ) THEN tmp = (qvap(mgs) - qvs(mgs))/(1. + qvs(mgs)*felv(mgs)**2/(cp*rw*temg(mgs)**2) ) @@ -23402,7 +24209,11 @@ subroutine nssl_2mom_gs & ltemq = (temg(mgs)-163.15)/fqsat+1.5 ltemq = Min( nqsat, Max(1,ltemq) ) - qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + IF ( iqvsopt == 0 ) THEN + qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + qvs(mgs) = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF qis(mgs) = pqs(mgs)*tabqis(ltemq) qss(mgs) = qvs(mgs) if ( temg(mgs) .lt. tfr ) then @@ -23470,7 +24281,7 @@ subroutine nssl_2mom_gs & IF ( eqtset > 2 ) THEN pipert(mgs) = pipert(mgs) & & +(felspi(mgs)*dqci(mgs) & - & +felvpi(mgs)*dqcw(mgs))*dtp + & +felvpi(mgs)*dqcw(mgs)) ! *dtp (remove dtp since dqxx are not rates) ENDIF end if ! dqwv(mgs) .lt. 0. (end of evap/sublim) @@ -23540,7 +24351,7 @@ subroutine nssl_2mom_gs & IF ( eqtset > 2 ) THEN pipert(mgs) = pipert(mgs) + (0 & & +felspi(mgs)*dqci(mgs) & - & +felvpi(mgs)*dqcw(mgs))*dtp + & +felvpi(mgs)*dqcw(mgs)) ! *dtp (remove dtp since dqxx are not rates) ENDIF qwvp(mgs) = qwvp(mgs) - ( dqvcnd(mgs) ) @@ -23564,7 +24375,12 @@ subroutine nssl_2mom_gs & tqvcon = temg(mgs)-cbw ltemq = (temg(mgs)-163.15)/fqsat+1.5 ltemq = Min( nqsat, Max(1,ltemq) ) - qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + + IF ( iqvsopt == 0 ) THEN + qvs(mgs) = pqs(mgs)*tabqvs(ltemq) + ELSEIF ( iqvsopt == 1 ) THEN + qvs(mgs) = rdorv*esbolton*tabqvs(ltemq)/(pres(mgs) - esbolton*tabqvs(ltemq)) + ENDIF qis(mgs) = pqs(mgs)*tabqis(ltemq) qx(mgs,lc) = max( 0.0, qx(mgs,lc) ) qitmp(mgs) = max( 0.0, qitmp(mgs) ) @@ -23718,7 +24534,7 @@ subroutine nssl_2mom_gs & IF ( qx(mgs,il) .le. 0.0 ) THEN cx(mgs,il) = 0.0 ELSE !{ - IF ( cx(mgs,il) .gt. cxmin ) THEN !{ + IF ( cx(mgs,il) .gt. cxmin .and. qx(mgs,il) > qxmin(il) ) THEN !{ only do this if mass is sufficient ! xv(mgs,il) = rho0(mgs)*qx(mgs,il)/(xdn(mgs,il)*Max(1.0e-9,cx(mgs,il))) ! xv(mgs,il) = rho0(mgs)*qx(mgs,il)/(xdn(mgs,il)*Max(cxmin,cx(mgs,il))) xv(mgs,il) = rho0(mgs)*qx(mgs,il)/(xdn(mgs,il)*cx(mgs,il)) @@ -23729,7 +24545,8 @@ subroutine nssl_2mom_gs & ! 8/26/2015 erm: apply imaxdiaopt for 2-moment also IF ( imaxdiaopt == 1 .or. il == lc .or. il == li .or. (il == lr .and. imurain == 3) .or. & - & (il == ls .and. imusnow == 3 ) ) THEN + & (il == ls .and. imusnow == 3 ) .or. ( il >= lh .and. lh > 0 ) ) THEN +! IF ( imaxdiaopt == 1 .or. (il == lr .and. imurain == 3) .or. .not. (il == lr .and. imurain == 1) ) THEN xvbarmax = xvmx(il) ELSEIF ( imaxdiaopt == 2 ) THEN ! test against maximum mass diameter xvbarmax = xvmx(il) /((3. + alpha(mgs,il))**3/((3. + alpha(mgs,il))*(2. + alpha(mgs,il))*(1. + alpha(mgs,il)))) @@ -24053,8 +24870,8 @@ subroutine nssl_2mom_gs & qr = qx(mgs,il) z = zx(mgs,il) - IF ( zx(mgs,il) .gt. 0. ) THEN !{ - + IF ( zx(mgs,il) .gt. zxmin .and. qr > qxmin(il) .and. chw > cxmin ) THEN !{ + ! rdi = z*(pi/6.*1000.)**2*chw/((rho0(mgs)*qr)**2) rdi = z*(pi/6.*xdn(mgs,il))**2*chw/((rho0(mgs)*qr)**2) diff --git a/phys/module_physics_init.F b/phys/module_physics_init.F index 6189996ea7..47d058184b 100644 --- a/phys/module_physics_init.F +++ b/phys/module_physics_init.F @@ -4669,11 +4669,12 @@ SUBROUTINE mp_init(RAINNC,SNOWNC,GRAUPELNC,config_flags,restart,warm_rain, ccn_conc = config_flags%nssl_cccn/1.225 ! set this to have correct boundary conditions ENDIF CALL nssl_2mom_init(nssl_params=nssl_params,ipctmp=nssl_ipconc,mixphase=0, & - nssl_density_on=(config_flags%nssl_density_on > 0), & - nssl_hail_on=config_flags%nssl_hail_on > 0, & - nssl_ccn_on=(config_flags%nssl_ccn_on > 0), & - nssl_icdx=config_flags%nssl_icdx, & - nssl_icdxhl=config_flags%nssl_icdxhl,ccn_is_ccna=config_flags%nssl_ccn_is_ccna) + nssl_density_on=(config_flags%nssl_density_on > 0), & + nssl_hail_on=config_flags%nssl_hail_on > 0, & + nssl_ccn_on=( config_flags%nssl_ccn_on > 0 ), & + nssl_icdx=config_flags%nssl_icdx, & + nssl_icdxhl=config_flags%nssl_icdxhl, & + ccn_is_ccna=config_flags%nssl_ccn_is_ccna) #endif #if (EM_CORE==1) CASE (CAMMGMPSCHEME) ! CAM5's microphysics From 93c73337aa5a2ad77edb411d8a24f66aebe7ab2b Mon Sep 17 00:00:00 2001 From: Ted Mansell Date: Sun, 23 Feb 2025 16:29:22 -0600 Subject: [PATCH 2/8] module_mp_nssl_2mom.F : Added missing final fixes --- phys/module_mp_nssl_2mom.F | 44 ++++++++++++++++++++++---------------- 1 file changed, 25 insertions(+), 19 deletions(-) diff --git a/phys/module_mp_nssl_2mom.F b/phys/module_mp_nssl_2mom.F index 980106d236..6a58e2507c 100644 --- a/phys/module_mp_nssl_2mom.F +++ b/phys/module_mp_nssl_2mom.F @@ -1,6 +1,6 @@ !WRF:MODEL_LAYER:PHYSICS -! prepocessed on "Feb 19 2025" at "16:07:57" +! prepocessed on "Feb 23 2025" at "16:22:22" @@ -1262,6 +1262,7 @@ SUBROUTINE nssl_2mom_init( & & nssl_alphahl, & & nssl_alphar, & & nssl_density_on, nssl_hail_on, nssl_ccn_on, nssl_icecrystals_on, ccn_is_ccna, & + & nssl_ccn_opt, & & infileunit, & & myrank, mpiroot & ) @@ -1283,8 +1284,8 @@ SUBROUTINE nssl_2mom_init( & & nssl_icdx, & & nssl_icdxhl, myrank, mpiroot, & & nssl_ufccn, & - & nssl_ccn_on - logical, intent(in), optional :: nssl_density_on, nssl_hail_on, nssl_icecrystals_on + & nssl_ccn_opt + logical, intent(in), optional :: nssl_density_on, nssl_ccn_on, nssl_hail_on, nssl_icecrystals_on integer, intent(inout), optional :: ccn_is_ccna integer, intent(in),optional :: infileunit @@ -1482,8 +1483,11 @@ SUBROUTINE nssl_2mom_init( & ENDIF IF ( present( nssl_ccn_on ) ) THEN - IF ( nssl_ccn_on > 0 ) THEN + IF ( nssl_ccn_on ) THEN ccn_on = 1 + IF ( present( nssl_ccn_opt ) ) THEN + IF ( nssl_ccn_opt > 10 ) ac_opt = 22 + ENDIF ELSE ccn_on = 0 irenuc = 2 @@ -1493,8 +1497,8 @@ SUBROUTINE nssl_2mom_init( & IF ( irenuc >= 5 ) THEN turn_on_ccna = .true. IF ( present( nssl_ccn_on ) ) THEN - IF ( nssl_ccn_on <= 0 ) THEN - write(0,*) 'NSSL_MP Error: Must have nssl_ccn_on=1 for irenuc >= 5!' + IF ( .not. nssl_ccn_on ) THEN + write(0,*) 'NSSL_MP Error: Must have nssl_ccn_on=1/true for irenuc >= 5!' STOP ENDIF ENDIF @@ -3999,6 +4003,19 @@ SUBROUTINE GAMMADP(X,GA) integer :: k,m1,m double precision :: G(26) + + DATA G/1.0D0,0.5772156649015329D0, & + & -0.6558780715202538D0, -0.420026350340952D-1, & + & 0.1665386113822915D0,-.421977345555443D-1, & + & -.96219715278770D-2, .72189432466630D-2, & + & -.11651675918591D-2, -.2152416741149D-3, & + & .1280502823882D-3, -.201348547807D-4, & + & -.12504934821D-5, .11330272320D-5, & + & -.2056338417D-6, .61160950D-8, & + & .50020075D-8, -.11812746D-8, & + & .1043427D-9, .77823D-11, & + & -.36968D-11, .51D-12, & + & -.206D-13, -.54D-14, .14D-14, .1D-15/ IF (X.EQ.INT(X)) THEN IF (X.GT.0.0D0) THEN @@ -4022,18 +4039,6 @@ SUBROUTINE GAMMADP(X,GA) ELSE Z=X ENDIF - DATA G/1.0D0,0.5772156649015329D0, & - & -0.6558780715202538D0, -0.420026350340952D-1, & - & 0.1665386113822915D0,-.421977345555443D-1, & - & -.96219715278770D-2, .72189432466630D-2, & - & -.11651675918591D-2, -.2152416741149D-3, & - & .1280502823882D-3, -.201348547807D-4, & - & -.12504934821D-5, .11330272320D-5, & - & -.2056338417D-6, .61160950D-8, & - & .50020075D-8, -.11812746D-8, & - & .1043427D-9, .77823D-11, & - & -.36968D-11, .51D-12, & - & -.206D-13, -.54D-14, .14D-14, .1D-15/ GR=G(26) DO K=25,1,-1 GR=GR*Z+G(K) @@ -12915,7 +12920,8 @@ subroutine nssl_2mom_gs & real ac1,bc, c1,d1,e1,f1,p380,tmp,tmp1,tmp2,tmp3,tmp4,tmp5,tmp6,temp3 ! , sstdy, super real :: flim, xmass real dw,dwr - double precision :: tmpc, tmpz, tmpzmlt + double precision :: tmpz, tmpzmlt + real :: tmpc real ratio, delx, dely real dbigg,volt real chgtmp,fac,mixedphasefac From 695ff560d697cc1c5c5c0638d7987028b6a63dbc Mon Sep 17 00:00:00 2001 From: Ted Mansell Date: Mon, 24 Feb 2025 15:53:36 -0600 Subject: [PATCH 3/8] Changed nssl_hail_on from size of max_dom to 1 because it has to be the same on all domains. Also deprecate nssl_ccn_on=0 (i.e., strongly recommend nssl_ccn_on=1) --- Registry/Registry.EM_COMMON | 2 +- share/module_check_a_mundo.F | 32 +++++++++++++++++--------------- 2 files changed, 18 insertions(+), 16 deletions(-) diff --git a/Registry/Registry.EM_COMMON b/Registry/Registry.EM_COMMON index aacb32b5a6..7ba37e4f58 100644 --- a/Registry/Registry.EM_COMMON +++ b/Registry/Registry.EM_COMMON @@ -2419,7 +2419,7 @@ rconfig real nssl_rho_qhl namelist,physics 1 900 rconfig real nssl_rho_qs namelist,physics 1 100. rh "Snow particle density" "" "" rconfig integer nssl_icdx namelist,physics 1 6 rh "NSSL Graupel fall speed option" "" "" rconfig integer nssl_icdxhl namelist,physics 1 6 rh "NSSL Hail fall speed option" "" "" -rconfig integer nssl_hail_on namelist,physics max_domains -1 rh "NSSL Hail flag" "" "" +rconfig integer nssl_hail_on namelist,physics 1 -1 rh "NSSL Hail flag" "" "" rconfig integer nssl_ccn_on namelist,physics 1 -1 rh "NSSL CCN flag" "" "" rconfig integer nssl_ccn_is_ccna namelist,physics 1 0 rh "NSSL flag that CCN is CCNA" "" "" rconfig integer nssl_2moment_on namelist,physics 1 -1 rh "NSSL 2-moment flag" "" "" diff --git a/share/module_check_a_mundo.F b/share/module_check_a_mundo.F index c655474244..5436cde591 100644 --- a/share/module_check_a_mundo.F +++ b/share/module_check_a_mundo.F @@ -773,7 +773,7 @@ END FUNCTION bep_bem_ngr_u !----------------------------------------------------------------------- ! Check that all values of mp_physics are the same for all domains !----------------------------------------------------------------------- - +! NOTE: This is redundant because frame/module_configure.F has already set all values to the value of mp_physics(max_dom) DO i = 2, model_config_rec % max_dom IF ( .NOT. model_config_rec % grid_allowed(i) ) CYCLE IF ( model_config_rec % mp_physics(i) .NE. & @@ -3382,27 +3382,27 @@ SUBROUTINE set_physics_rconfigs IF ( model_config_rec % mp_physics(i) .EQ. 22 ) THEN model_config_rec % mp_physics(i) = NSSL_2MOM model_config_rec % nssl_2moment_on = 1 - model_config_rec % nssl_hail_on(i) = 0 - model_config_rec % nssl_ccn_on = 0 + model_config_rec % nssl_hail_on = 0 + model_config_rec % nssl_ccn_on = 1 model_config_rec % nssl_density_on = 1 ! set graupel density WRITE (wrf_err_message, FMT='(A)') ' **CAUTION** mp_physics = 22 has been deprecated. '// & - 'Instead you can use mp_physics=18, nssl_hail_on=0, nssl_ccn_on=0' + 'Instead you can use mp_physics=18, nssl_hail_on=0, nssl_ccn_on=1' CALL wrf_debug ( 0, wrf_err_message ) ELSEIF ( model_config_rec % mp_physics(i) .EQ. 17 ) THEN model_config_rec % mp_physics(i) = NSSL_2MOM model_config_rec % nssl_2moment_on = 1 - model_config_rec % nssl_hail_on(i) = 1 - model_config_rec % nssl_ccn_on = 0 + model_config_rec % nssl_hail_on = 1 + model_config_rec % nssl_ccn_on = 1 model_config_rec % nssl_density_on = 2 ! set graupel+hail density WRITE (wrf_err_message, FMT='(A)') ' **CAUTION** mp_physics = 17 has been deprecated. '// & - 'Instead you can use mp_physics=18, nssl_ccn_on=0' + 'Please use mp_physics=18, Note that nssl_ccn_on=0 is not recommended and default set to 1' ! print statement for deprecated option CALL wrf_debug ( 0, wrf_err_message ) ELSEIF ( model_config_rec % mp_physics(i) .EQ. 19 ) THEN ! single-moment with hail + graupel density model_config_rec % mp_physics(i) = NSSL_2MOM model_config_rec % nssl_2moment_on = 0 - model_config_rec % nssl_hail_on(i) = 2 + model_config_rec % nssl_hail_on = 2 model_config_rec % nssl_density_on = 1 ! set graupel density ! print statement for deprecated option WRITE (wrf_err_message, FMT='(A)') ' **CAUTION** mp_physics = 19 has been deprecated. '// & @@ -3412,7 +3412,7 @@ SUBROUTINE set_physics_rconfigs ! single-moment without model_config_rec % mp_physics(i) = NSSL_2MOM model_config_rec % nssl_2moment_on = 0 - model_config_rec % nssl_hail_on(i) = 0 + model_config_rec % nssl_hail_on = 0 model_config_rec % nssl_density_on = 0 ! set graupel density ! print statement for deprecated option WRITE (wrf_err_message, FMT='(A)') ' **CAUTION** mp_physics = 21 has been deprecated. '// & @@ -3423,7 +3423,7 @@ SUBROUTINE set_physics_rconfigs IF ( model_config_rec % mp_physics(i) /= NSSL_2MOM ) THEN ! If not NSSL-MP, make sure extra fields are turned off (in case of stray namelist settings) model_config_rec % nssl_2moment_on = 0 - model_config_rec % nssl_hail_on(i) = 0 + model_config_rec % nssl_hail_on = 0 model_config_rec % nssl_density_on = 0 ! set graupel density model_config_rec % nssl_3moment = 0 model_config_rec % nssl_ccn_on = 0 @@ -3431,23 +3431,25 @@ SUBROUTINE set_physics_rconfigs ELSE ! make sure settings are consistent IF ( model_config_rec % nssl_ccn_on < 0 ) THEN + ! nssl_ccn_on = -1, so set default model_config_rec % nssl_ccn_on = 1 ENDIF IF ( model_config_rec % nssl_2moment_on < 0 ) THEN ! turn on number concentrations + ! nssl_2moment_on = -1 so set default model_config_rec % nssl_2moment_on = 1 ENDIF - IF ( model_config_rec % nssl_hail_on(i) < 0 ) THEN + IF ( model_config_rec % nssl_hail_on < 0 ) THEN IF ( model_config_rec % nssl_2moment_on == 0 ) THEN - model_config_rec % nssl_hail_on(i) = 2 + model_config_rec % nssl_hail_on = 2 ELSE - model_config_rec % nssl_hail_on(i) = 1 + model_config_rec % nssl_hail_on = 1 ENDIF ENDIF IF ( model_config_rec % nssl_density_on < 0 ) THEN - IF ( model_config_rec % nssl_hail_on(i) == 1 ) THEN + IF ( model_config_rec % nssl_hail_on == 1 ) THEN model_config_rec % nssl_density_on = 2 ! set default of graupel+hail density ELSE model_config_rec % nssl_density_on = 1 ! set graupel density (hail off) @@ -3456,7 +3458,7 @@ SUBROUTINE set_physics_rconfigs IF ( model_config_rec % nssl_3moment == 1 ) THEN model_config_rec % nssl_2moment_on = 1 - IF ( model_config_rec % nssl_hail_on(i) == 1 ) THEN + IF ( model_config_rec % nssl_hail_on == 1 ) THEN model_config_rec % nssl_3moment = 2 ! 3mom rain, graupel and hail ELSE model_config_rec % nssl_3moment = 1 ! 3mom rain and graupel (no hail) From 0a72336086eb89aa4634177d35b337d84fe6cdae Mon Sep 17 00:00:00 2001 From: Ted Mansell Date: Thu, 27 Feb 2025 19:01:06 -0600 Subject: [PATCH 4/8] module_mp_nssl_2mom.F : Fixes for single moment and no-hail options. --- doc/README.NSSLmp | 12 ++++++------ phys/module_mp_nssl_2mom.F | 14 ++++++++++---- share/module_check_a_mundo.F | 1 + 3 files changed, 17 insertions(+), 10 deletions(-) diff --git a/doc/README.NSSLmp b/doc/README.NSSLmp index 31ef9409aa..915b63cea8 100644 --- a/doc/README.NSSLmp +++ b/doc/README.NSSLmp @@ -22,16 +22,16 @@ Basic options in physics namelist: CCN concentration + options The legacy options (17,19,21,22) still behave as before (for now), but going - forward one should use mp_physics=18 with modifier flags: + forward one should use mp_physics=18 with modifier flags. 2025 Update, however, sets nssl_ccn_on=1 by default (keeps supersaturation much more reasonable; except for single moment). mp_physics = 22 ! NSSL scheme (2-moment) without hail - Equivalent: mp=18, nssl_hail_on=0, nssl_ccn_on=0 - = 17 ! NSSL scheme (2-moment) with hail with constant background CCN - concentration - Equivalent: mp=18, nssl_ccn_on=0 + Equivalent: mp=18, nssl_hail_on=0, nssl_ccn_on=1 + = 17 ! NSSL scheme (2-moment) with hail is now the same as mp=18 + Equivalent: mp=18, nssl_ccn_on=1 <- must explicitly set nssl_ccn_on=0 to + get old behavior = 19, NSSL 1-moment (7 class: qv,qc,qr,qi,qs,qg,qh; predicts graupel density) - Equivalent: mp=18, nssl_2moment_on=0, nssl_ccn_on=0 (do no set nssl_hail_on) + Equivalent: mp=18, nssl_2moment_on=0, nssl_ccn_on=1 (do no set nssl_hail_on) = 21, NSSL 1-moment, (6-class), very similar to Gilmore et al. 2004 Equivalent: mp=18, nssl_2moment_on=0, nssl_hail_on=0, nssl_ccn_on=0, nssl_density_on=0 diff --git a/phys/module_mp_nssl_2mom.F b/phys/module_mp_nssl_2mom.F index 6a58e2507c..5b35709c66 100644 --- a/phys/module_mp_nssl_2mom.F +++ b/phys/module_mp_nssl_2mom.F @@ -1,6 +1,6 @@ !WRF:MODEL_LAYER:PHYSICS -! prepocessed on "Feb 23 2025" at "16:22:22" +! prepocessed on "Feb 27 2025" at "18:52:58" @@ -1498,7 +1498,7 @@ SUBROUTINE nssl_2mom_init( & turn_on_ccna = .true. IF ( present( nssl_ccn_on ) ) THEN IF ( .not. nssl_ccn_on ) THEN - write(0,*) 'NSSL_MP Error: Must have nssl_ccn_on=1/true for irenuc >= 5!' + write(0,*) 'NSSL_MP Error: Must have nssl_ccn_on=1/true for irenuc >= 5!' STOP ENDIF ENDIF @@ -8063,7 +8063,7 @@ subroutine ziegfall1d(nx,ny,nz,nor,norz,na,dtp,jgs,ixcol, & ! Vaporize tiny values DO il = l1,l2 - IF ( lz(il) < 1 ) THEN + IF ( lz(il) < 1 .and. ln(il) > 1 ) THEN do mgs = 1,ngscnt IF ( cx(mgs,il) <= cxmin .or. qx(mgs,il) < qxmin(il) ) THEN cx(mgs,il) = 0.0 @@ -12242,6 +12242,7 @@ subroutine smallvalues & hwdn = xdn0(lhl) ENDIF + IF ( ipconc >= 5 .and. an(ix,jy,kz,lhl) .gt. qxmin(lhl) ) THEN qr = an(ix,jy,kz,lhl) xvol = dn(ix,jy,kz)*an(ix,jy,kz,lhl)/(hwdn*an(ix,jy,kz,lnhl)) chw = an(ix,jy,kz,lnhl) @@ -12251,6 +12252,7 @@ subroutine smallvalues & chw = dn(ix,jy,kz)*an(ix,jy,kz,lhl)/(xvol*hwdn) an(ix,jy,kz,lnhl) = chw ENDIF + ENDIF ! CHECK INTERCEPT IF ( ipconc == 5 .and. an(ix,jy,kz,lhl) .gt. qxmin(lhl) .and. alphahl .le. 0.1 .and. lnhl .gt. 1 .and. lzhl == 0 ) THEN @@ -12398,6 +12400,7 @@ subroutine smallvalues & hwdn = xdn0(lh) ENDIF + IF ( ipconc >= 5 .and. an(ix,jy,kz,lh) .gt. qxmin(lh) ) THEN qr = an(ix,jy,kz,lh) xvol = dn(ix,jy,kz)*an(ix,jy,kz,lh)/(hwdn*an(ix,jy,kz,lnh)) chw = an(ix,jy,kz,lnh) @@ -12407,6 +12410,7 @@ subroutine smallvalues & chw = dn(ix,jy,kz)*an(ix,jy,kz,lh)/(xvol*hwdn) an(ix,jy,kz,lnh) = chw ENDIF + ENDIF ! CHECK INTERCEPT IF ( ipconc == 5 .and. an(ix,jy,kz,lh) .gt. qxmin(lh) .and. alphah .le. 0.1 .and. lnh .gt. 1 .and. lzh == 0 ) THEN @@ -14392,8 +14396,10 @@ subroutine nssl_2mom_gs & & ((3.0 + alphar)*(2.0 + alphar)*(1.0 + alphar)) g1x(:,lh) = (6.0 + alphah)*(5.0 + alphah)*(4.0 + alphah)/ & & ((3.0 + alphah)*(2.0 + alphah)*(1.0 + alphah)) - g1x(:,lhl) = (6.0 + alphahl)*(5.0 + alphahl)*(4.0 + alphahl)/ & + IF ( lhl > 0 ) THEN + g1x(:,lhl) = (6.0 + alphahl)*(5.0 + alphahl)*(4.0 + alphahl)/ & & ((3.0 + alphahl)*(2.0 + alphahl)*(1.0 + alphahl)) + ENDIF ENDIF scx(:,:) = 0.0 diff --git a/share/module_check_a_mundo.F b/share/module_check_a_mundo.F index 5436cde591..3e8f1d12b5 100644 --- a/share/module_check_a_mundo.F +++ b/share/module_check_a_mundo.F @@ -3414,6 +3414,7 @@ SUBROUTINE set_physics_rconfigs model_config_rec % nssl_2moment_on = 0 model_config_rec % nssl_hail_on = 0 model_config_rec % nssl_density_on = 0 ! set graupel density + model_config_rec % nssl_ccn_on = 0 ! print statement for deprecated option WRITE (wrf_err_message, FMT='(A)') ' **CAUTION** mp_physics = 21 has been deprecated. '// & 'Instead you can use mp_physics=18, nssl_2moment_on=0, nssl_ccn_on=0, nssl_hail_on=0' From d32870f1510184936ac5fced6b1b8dcf32f0c8b5 Mon Sep 17 00:00:00 2001 From: Ted Mansell Date: Tue, 11 Mar 2025 10:07:02 -0500 Subject: [PATCH 5/8] Splits a couple long lines that gfortran might complain about and fixes hail production for the unlikely case of a user choosing single moment with hail and without graupel density and graupel density set less than 700. --- phys/module_mp_nssl_2mom.F | 45 +++++++++++++++++++++----------------- 1 file changed, 25 insertions(+), 20 deletions(-) diff --git a/phys/module_mp_nssl_2mom.F b/phys/module_mp_nssl_2mom.F index 5b35709c66..fea2d5db13 100644 --- a/phys/module_mp_nssl_2mom.F +++ b/phys/module_mp_nssl_2mom.F @@ -1,6 +1,6 @@ !WRF:MODEL_LAYER:PHYSICS -! prepocessed on "Feb 27 2025" at "18:52:58" +! prepocessed on "Mar 11 2025" at "09:41:58" @@ -1419,12 +1419,11 @@ SUBROUTINE nssl_2mom_init( & ipconc = ipctmp - IF ( ipconc < 5 ) THEN - ihlcnh = 0 - ENDIF IF ( ihlcnh <= 0 ) THEN - IF ( ipconc == 5 ) THEN + IF ( ipconc < 5 ) THEN + ihlcnh = 0 + ELSEIF ( ipconc == 5 ) THEN ihlcnh = 3 ELSEIF ( ipconc >= 6 ) THEN ihlcnh = 3 @@ -1747,9 +1746,15 @@ SUBROUTINE nssl_2mom_init( & IF ( ipconc == 0 ) THEN IF ( hail_on == 1 ) THEN ! turn on graupel density for 1-moment scheme - lvh = 9 - ltmp = 9 - denscale(lvh) = 1 + IF ( density_on >= 1 ) THEN ! turn on graupel density for 1-moment scheme + lvh = 9 + ltmp = 9 + denscale(lvh) = 1 + ELSE + ltmp = lhab + lvh = 0 + lvhl = 0 + ENDIF ELSE ! no hail, 'LFO' scheme ltmp = lhab lhl = 0 @@ -2924,7 +2929,7 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw an(ix,1,kz,lnh) = chw(ix,kz,jy) IF ( lhl > 1 ) an(ix,1,kz,lnhl) = chl(ix,kz,jy) ENDIF - IF ( lvh > 0 ) an(ix,1,kz,lvh) = vhw(ix,kz,jy) + IF ( lvh > 0 .and. present( vhw ) ) an(ix,1,kz,lvh) = vhw(ix,kz,jy) IF ( lvhl > 0 .and. present( vhl ) ) an(ix,1,kz,lvhl) = vhl(ix,kz,jy) IF ( ipconc >= 6 ) THEN @@ -3518,7 +3523,7 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw - IF ( lvh > 0 ) vhw(ix,kz,jy) = an(ix,1,kz,lvh) + IF ( lvh > 0 .and. present( vhw ) ) vhw(ix,kz,jy) = an(ix,1,kz,lvh) IF ( lvhl > 0 .and. present( vhl ) ) vhl(ix,kz,jy) = an(ix,1,kz,lvhl) #if ( WRF_CHEM == 1 ) @@ -12579,7 +12584,8 @@ subroutine smallvalues & ! write(0,*) 'restore: k, qccn,exp = ',kz,qccn,dn(ix,jy,kz)*qccn,Exp(-dtp/ccntimeconst) ! write(0,*) 'ccn1,ccn2 = ',an(ix,jy,kz,lccn),dn(ix,jy,kz)*qccn - Max(0.0 , dn(ix,jy,kz)*qccn - an(ix,jy,kz,lccn))*Exp(-dtp/ccntimeconst) ! ENDIF - IF ( an(ix,jy,kz,lccn) > 1. .and. tmp < qxmin(li) .and. ( an(ix,jy,kz,lccn) < dn(ix,jy,kz)*qccn .or. .not. invertccn ) ) THEN + IF ( an(ix,jy,kz,lccn) > 1. .and. tmp < qxmin(li) .and. & + ( an(ix,jy,kz,lccn) < dn(ix,jy,kz)*qccn .or. .not. invertccn ) ) THEN ! an(ix,jy,kz,lccn) = & ! an(ix,jy,kz,lccn) + Max(0.0 , dn(ix,jy,kz)*qccn - an(ix,jy,kz,lccn))*(1.0 - Exp(-dtp/ccntimeconst)) ! Equivalent form after expanding last term: @@ -13092,10 +13098,9 @@ subroutine nssl_2mom_gs & real :: alpha(ngs,lc:lhab) real :: dab0lh(ngs,lc:lhab,lc:lhab) real :: dab1lh(ngs,lc:lhab,lc:lhab) - real :: zx(ngs,lr:lhab) - real :: zxmxd(ngs,lr:lhab) - real :: g1x(ngs,lr:lhab) - + real :: zx(ngs,lr:lhab) + real :: zxmxd(ngs,lr:lhab) + real :: g1x(ngs,lr:lhab) real :: g1xmax,g1xmin real :: qsimxdep(ngs) ! max sublimation of qi+qs+qis real :: qsimxsub(ngs) ! max depositionof qi+qs+qis @@ -17842,7 +17847,8 @@ subroutine nssl_2mom_gs & tmp = crcnw(mgs) tmp2 = qrcnw(mgs)*cx(mgs,lr)/qx(mgs,lr) ! try mass*diameter-weighted average of old and new Dmr (using full qc mass) - crcnw(mgs) = (tmp*xdia(mgs,lc,3)*qx(mgs,lc)+tmp2*xdia(mgs,lr,3)*qx(mgs,lr))/(xdia(mgs,lc,3)*qx(mgs,lc)+xdia(mgs,lr,3)*qx(mgs,lr)) + crcnw(mgs) = (tmp*xdia(mgs,lc,3)*qx(mgs,lc)+tmp2*xdia(mgs,lr,3)*qx(mgs,lr))/ & + (xdia(mgs,lc,3)*qx(mgs,lc)+xdia(mgs,lr,3)*qx(mgs,lr)) ELSEIF ( ( dmropt == 7 ) .and. qx(mgs,lr) > qxmin(lr) ) THEN tmp = crcnw(mgs) tmp2 = qrcnw(mgs)*cx(mgs,lr)/qx(mgs,lr) @@ -17852,7 +17858,8 @@ subroutine nssl_2mom_gs & tmp = crcnw(mgs) tmp2 = qrcnw(mgs)*cx(mgs,lr)/qx(mgs,lr) ! try sqrt(diameter)-weighted average of old and new Dmr - crcnw(mgs) = (tmp*sqrt(xdia(mgs,lc,3))+tmp2*sqrt(xdia(mgs,lr,3)))/(sqrt(xdia(mgs,lc,3))+sqrt(xdia(mgs,lr,3))) + crcnw(mgs) = (tmp*sqrt(xdia(mgs,lc,3))+tmp2*sqrt(xdia(mgs,lr,3)))/ & + (sqrt(xdia(mgs,lc,3))+sqrt(xdia(mgs,lr,3))) ENDIF ELSEIF ( dmrauto == 1 .and. cx(mgs,lr) > cxmin) THEN IF ( qx(mgs,lr) > qxmin(lr) ) THEN @@ -19165,8 +19172,6 @@ subroutine nssl_2mom_gs & ELSEIF ( ibinhlmlr == 1 ) THEN ! use incomplete gamma functions to approximate the bin results -! #ifdef 1 -! #if (defined 1) && defined( COMMAS ) || defined( COMMASTMP ) ELSEIF ( ibinhlmlr == -1 ) THEN ! OLD VERSION use incomplete gamma functions to approximate the bin results @@ -20964,7 +20969,7 @@ subroutine nssl_2mom_gs & ! qhlcnh(mgs) = 0.0 ! chlcnh(mgs) = 0.0 if ( wetgrowth(mgs) .and. temg(mgs) .lt. tfr-5. .and. qx(mgs,lh) > qxmin(lh) ) then - if ( qhacw(mgs).gt.1.e-6 .and. xdn(mgs,lh) > 700. ) then + if ( qhacw(mgs).gt.1.e-6 .and. ( xdn(mgs,lh) > 700. .or. lvh == 0 ) ) then qhlcnh(mgs) = & ((pi*xdn(mgs,lh)*cx(mgs,lh)) / (6.0*rho0(mgs)*dtp)) & *exp(-hldia1/xdia(mgs,lh,1)) & From d2426166a7f37eb5c1a5e27c9bc19876c4bee1b2 Mon Sep 17 00:00:00 2001 From: Ted Mansell Date: Thu, 13 Mar 2025 12:15:37 -0500 Subject: [PATCH 6/8] solve_em.F and module_microphysics_driver.F : Removed new flags for ssat/ssati and use config_flags instead module_check_a_mundo.F : filled out a comment --- dyn_em/solve_em.F | 2 -- phys/module_microphysics_driver.F | 8 +++----- share/module_check_a_mundo.F | 2 +- 3 files changed, 4 insertions(+), 8 deletions(-) diff --git a/dyn_em/solve_em.F b/dyn_em/solve_em.F index cafc08b37a..e29c86a358 100644 --- a/dyn_em/solve_em.F +++ b/dyn_em/solve_em.F @@ -3730,8 +3730,6 @@ END SUBROUTINE CMAQ_DRIVER & ,RHO=grid%rho ,SPEC_ZONE=grid%spec_zone & & ,SR=grid%sr ,TH=th_phy & & ,ssat=grid%ssat, ssati=grid%ssati & - & ,f_ssat=( SIZE(grid%ssat) .GT. 1 ) & - & ,f_ssati=( SIZE(grid%ssati) .GT. 1 ) & & ,refl_10cm=grid%refl_10cm & ! hm, 9/22/09 for refl & ,vmi3d=grid%vmi3d & ! for P3 & ,di3d=grid%di3d & ! for P3 diff --git a/phys/module_microphysics_driver.F b/phys/module_microphysics_driver.F index ab0345b8c4..56fc82fd63 100644 --- a/phys/module_microphysics_driver.F +++ b/phys/module_microphysics_driver.F @@ -112,7 +112,6 @@ SUBROUTINE microphysics_driver( & ,qnwfa2d, qnifa2d, qnbca2d & ! for water/ice-friendly/black carbon aerosols ,qnocbb2d, qnbcbb2d & ! for biomass burning aerosols ,ssat,ssati & - ,f_ssat,f_ssati & ,refl_10cm & ! HM, 9/22/09, add for refl ,vmi3d & ! for P3 ,di3d & ! for P3 @@ -730,8 +729,7 @@ SUBROUTINE microphysics_driver( & ,f_qvoli,f_qaoli & ! for Jensen ISHMAEL ,f_qvoli2,f_qaoli2 & ! for Jensen ISHMAEL ,f_qi3,f_qni3,f_qvoli3,f_qaoli3 & ! for Jensen ISHMAEL - ,f_qnwfa, f_qnifa, f_qnbca & ! Added by G. Thompson - ,f_ssat,f_ssati + ,f_qnwfa, f_qnifa, f_qnbca ! Added by G. Thompson LOGICAL, OPTIONAL, INTENT(IN) :: diagflag @@ -2079,8 +2077,8 @@ SUBROUTINE microphysics_driver( & GRPLNCV = GRAUPELNCV, & SR=SR, & dbz = refl_10cm, & - ssat3d = ssat, f_ssat=f_ssat, & - ssati = ssati, f_ssati=f_ssati, & + ssat3d = ssat, f_ssat =(config_flags%nssl_ssat_output >= 1), & + ssati = ssati, f_ssati=(config_flags%nssl_ssat_output >= 2), & #if ( WRF_CHEM == 1 ) WETSCAV_ON = config_flags%wetscav_onoff == 1, & EVAPPROD=evapprod,RAINPROD=rainprod, & diff --git a/share/module_check_a_mundo.F b/share/module_check_a_mundo.F index 3e8f1d12b5..ebf57e437c 100644 --- a/share/module_check_a_mundo.F +++ b/share/module_check_a_mundo.F @@ -3409,7 +3409,7 @@ SUBROUTINE set_physics_rconfigs 'Instead you can use mp_physics=18, nssl_2moment_on=0, nssl_ccn_on=0' CALL wrf_debug ( 0, wrf_err_message ) ELSEIF ( model_config_rec % mp_physics(i) .EQ. 21 ) THEN - ! single-moment without + ! single-moment without graupel volume, no hail model_config_rec % mp_physics(i) = NSSL_2MOM model_config_rec % nssl_2moment_on = 0 model_config_rec % nssl_hail_on = 0 From 89cbe02e64e681e0582d3aa49fe76f69afe0b097 Mon Sep 17 00:00:00 2001 From: Ted Mansell Date: Tue, 1 Apr 2025 16:47:20 -0500 Subject: [PATCH 7/8] module_microphysics_driver.F : Pass namelist flag to NSSLmp instead of creating flags module_mp_nssl_2mom.F : Use ssat_output input parameter instead of logicals; also suggested changes from G. Firl --- phys/module_microphysics_driver.F | 5 ++-- phys/module_mp_nssl_2mom.F | 42 +++++++++++++++---------------- 2 files changed, 23 insertions(+), 24 deletions(-) diff --git a/phys/module_microphysics_driver.F b/phys/module_microphysics_driver.F index 56fc82fd63..dcc44c116a 100644 --- a/phys/module_microphysics_driver.F +++ b/phys/module_microphysics_driver.F @@ -2077,8 +2077,9 @@ SUBROUTINE microphysics_driver( & GRPLNCV = GRAUPELNCV, & SR=SR, & dbz = refl_10cm, & - ssat3d = ssat, f_ssat =(config_flags%nssl_ssat_output >= 1), & - ssati = ssati, f_ssati=(config_flags%nssl_ssat_output >= 2), & + ssat3d = ssat, & + ssati = ssati, & + nssl_ssat_output = config_flags%nssl_ssat_output, & #if ( WRF_CHEM == 1 ) WETSCAV_ON = config_flags%wetscav_onoff == 1, & EVAPPROD=evapprod,RAINPROD=rainprod, & diff --git a/phys/module_mp_nssl_2mom.F b/phys/module_mp_nssl_2mom.F index fea2d5db13..82b597fb9c 100644 --- a/phys/module_mp_nssl_2mom.F +++ b/phys/module_mp_nssl_2mom.F @@ -2387,7 +2387,7 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw elec_physics, & induc,elecz,scion,sciona,f_scion,f_sciona, & noninduc,noninducp,noninducn, & - ssat3d,ssati,f_ssat,f_ssati, & + ssat3d,ssati,nssl_ssat_output, & pcc2, pre2, depsubr, & mnucf2, melr2, ctr2, & rim1_2, rim2_2,rim3_2, & @@ -2435,7 +2435,7 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw qi,qhl,ccw,crw,cci,csw,chw,chl,vhw,vhl integer, optional, intent(in) :: is_theta_or_temp logical, optional, intent(in) :: f_zrw, f_zhw, f_zhl, f_vhw, f_vhl ! not used yet - logical, optional, intent(in) :: f_ssat,f_ssati + integer, optional, intent(in) :: nssl_ssat_output real, dimension(ims:ime, kms:kme, jms:jme), optional, intent(inout):: dbz, vzf, cn, cna, cni, cnuf real, dimension(ims:ime, kms:kme, jms:jme), optional, intent(inout):: cn_nu, cn_ac, cn_co, cinp, cna_co, cna_nu logical, optional, intent(in) :: f_cnnu, f_cnac, f_cnco, f_cinp, f_cnaco, f_cnanu @@ -3280,8 +3280,7 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw ENDDO ENDIF - IF ( ( ( present( ssat3d ) .and. present( f_ssat ) ) .or. & - ( present( ssati ) .and. present( f_ssati ) ) ) .and. makediag ) THEN + IF ( ( present( ssat3d ) .and. present( nssl_ssat_output ) ) .and. makediag ) THEN DO kz = kts,kte DO ix = its,ite @@ -3290,8 +3289,7 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw ltemq = Int( (temp1-163.15)/fqsat+1.5 ) ltemq = Min( nqsat, Max(1,ltemq) ) - IF ( present( ssat3d ) .and. present( f_ssat ) ) THEN - IF ( f_ssat ) THEN + IF ( present( ssat3d ) .and. nssl_ssat_output >= 1 ) THEN ! c1 = t00(ix,1,kz)*tabqvs(ltemq) IF ( iqvsopt == 0 ) THEN @@ -3304,15 +3302,12 @@ SUBROUTINE nssl_2mom_driver(qv, qc, qr, qi, qs, qh, qhl, ccw, crw, cci, csw, chw ssat3d(ix,kz,jy) = 100.*(an(ix,1,kz,lv)/c1 - 1.0) ! from "new" values ENDIF - ENDIF ENDIF - IF ( present( ssati ) .and. present( f_ssati ) ) THEN - IF ( f_ssati ) THEN + IF ( present( ssati ) .and. nssl_ssat_output >= 2 ) THEN t9s = (380.0/pn(ix,1,kz))*tabqis(ltemq) !saturation mixing ratio wrt ice ssati(ix,kz,jy) = 100.*(an(ix,1,kz,lv)/t9s - 1.0) ! Min(t8s,max(an(ix,1,kz,lv),0.0))/t9s ! qv/qvi ENDIF - ENDIF ENDDO ENDDO @@ -5941,6 +5936,9 @@ SUBROUTINE calc_eff_radius & double precision :: numh, numhl,denomh,denomhl logical :: flag_t4, flag_t5, flag_t6 + + real, parameter :: qmin = 1.e-8 + real, parameter :: volmin = 1.e-30 ! ------------------------------------------------------------------------------- @@ -6055,7 +6053,7 @@ SUBROUTINE calc_eff_radius & ENDIF IF ( present( t4 ) .and.( ( present(qrw) .and. present(crw) ) .or. flag_t4 ) ) THEN - IF ( qx(mgs,lr) > Max(1.e-8,qxmin(lr)) .and. cx(mgs,lr) > cxmin ) THEN + IF ( qx(mgs,lr) > Max(qmin,qxmin(lr)) .and. cx(mgs,lr) > cxmin ) THEN IF ( imurain == 1 ) THEN ! gamma-diameter ! Lambda for rain lam_r = factor_r *((xdn0(lr)*cx(mgs,lr))/(qx(mgs,lr)*rho0(mgs)))**(1./3.) @@ -6074,11 +6072,11 @@ SUBROUTINE calc_eff_radius & IF ( lhl < 1 .or. flag_t6 ) THEN ! graupel only - IF ( qx(mgs,lh) > Max(1.e-8,qxmin(lh)) ) THEN + IF ( qx(mgs,lh) > Max(qmin,qxmin(lh)) ) THEN ! Lambda for graupel hwdn = xdn0(lh) IF ( lvh > 1 ) THEN ! variable density - IF ( an(ix,jy,kz,lvh) > 1.e-30 ) THEN + IF ( an(ix,jy,kz,lvh) > volmin ) THEN hwdn = rho0(mgs)*qx(mgs,lh)/an(ix,jy,kz,lvh) ENDIF ENDIF @@ -6089,11 +6087,11 @@ SUBROUTINE calc_eff_radius & ELSE ! have hail, too, but do not have t6 array - IF ( qx(mgs,lh) > Max(1.e-8,qxmin(lh)) .and. qx(mgs,lhl) < Max(1.e-8,qxmin(lhl)) ) THEN + IF ( qx(mgs,lh) > Max(qmin,qxmin(lh)) .and. qx(mgs,lhl) < Max(qmin,qxmin(lhl)) ) THEN ! Lambda for graupel hwdn = xdn0(lh) IF ( lvh > 1 ) THEN ! variable density - IF ( an(ix,jy,kz,lvh) > 1.e-30 ) THEN + IF ( an(ix,jy,kz,lvh) > volmin ) THEN hwdn = rho0(mgs)*qx(mgs,lh)/an(ix,jy,kz,lvh) ENDIF ENDIF @@ -6101,11 +6099,11 @@ SUBROUTINE calc_eff_radius & lam_h = factor_h *((hwdn*cx(mgs,lh))/(qx(mgs,lh)*rho0(mgs)))**(1./3.) t5(ix,jy,kz) = 0.5*(alphah+3.)/lam_h - ELSEIF ( qx(mgs,lh) < Max(1.e-8,qxmin(lh)) .and. qx(mgs,lhl) > Max(1.e-8,qxmin(lhl)) ) THEN + ELSEIF ( qx(mgs,lh) < Max(qmin,qxmin(lh)) .and. qx(mgs,lhl) > Max(qmin,qxmin(lhl)) ) THEN ! Lambda for hail hldn = xdn0(lhl) IF ( lvhl > 1 ) THEN ! variable density - IF ( an(ix,jy,kz,lvhl) > 1.e-30 ) THEN + IF ( an(ix,jy,kz,lvhl) > volmin ) THEN hldn = rho0(mgs)*qx(mgs,lhl)/an(ix,jy,kz,lvhl) ENDIF ENDIF @@ -6113,19 +6111,19 @@ SUBROUTINE calc_eff_radius & lam_hl = factor_hl *((hldn*cx(mgs,lhl))/(qx(mgs,lhl)*rho0(mgs)))**(1./3.) t5(ix,jy,kz) = 0.5*(alphahl+3.)/lam_hl - ELSEIF ( qx(mgs,lh) > Max(1.e-8,qxmin(lh)) .and. qx(mgs,lhl) > Max(1.e-8,qxmin(lhl)) ) THEN + ELSEIF ( qx(mgs,lh) > Max(qmin,qxmin(lh)) .and. qx(mgs,lhl) > Max(qmin,qxmin(lhl)) ) THEN ! r_eff graupel and hail combined hldn = xdn0(lhl) IF ( lvhl > 1 ) THEN ! variable density - IF ( an(ix,jy,kz,lvhl) > 1.e-30 ) THEN + IF ( an(ix,jy,kz,lvhl) > volmin ) THEN hldn = rho0(mgs)*qx(mgs,lhl)/an(ix,jy,kz,lvhl) ENDIF ENDIF hwdn = xdn0(lh) IF ( lvh > 1 ) THEN ! variable density - IF ( an(ix,jy,kz,lvh) > 1.e-30 ) THEN + IF ( an(ix,jy,kz,lvh) > volmin ) THEN hwdn = rho0(mgs)*qx(mgs,lh)/an(ix,jy,kz,lvh) ENDIF ENDIF @@ -6150,11 +6148,11 @@ SUBROUTINE calc_eff_radius & IF ( present(t6) .and. flag_t6 .and. lhl > 1 ) THEN - IF ( qx(mgs,lhl) > Max(1.e-8,qxmin(lhl)) ) THEN + IF ( qx(mgs,lhl) > Max(qmin,qxmin(lhl)) ) THEN ! Lambda for hail hldn = xdn0(lhl) IF ( lvhl > 1 ) THEN ! variable density - IF ( an(ix,jy,kz,lvhl) > 1.e-30 ) THEN + IF ( an(ix,jy,kz,lvhl) > volmin ) THEN hldn = rho0(mgs)*qx(mgs,lhl)/an(ix,jy,kz,lvhl) ENDIF ENDIF From 84075e92bde7a14f505926662310996dc4a6475e Mon Sep 17 00:00:00 2001 From: Ted Mansell Date: Wed, 9 Apr 2025 12:05:49 -0500 Subject: [PATCH 8/8] module_mp_nssl_2mom.F : Bug fix for old droplet nucleation (irenuc=2) --- phys/module_mp_nssl_2mom.F | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/phys/module_mp_nssl_2mom.F b/phys/module_mp_nssl_2mom.F index 82b597fb9c..562564f528 100644 --- a/phys/module_mp_nssl_2mom.F +++ b/phys/module_mp_nssl_2mom.F @@ -10334,7 +10334,11 @@ SUBROUTINE NUCOND & DO mgs = 1,ngscnt ! default value of renucfrac is 0.0 IF ( irenuc /= 6 ) THEN - cnuc(mgs) = ccnc(mgs)*(1. - renucfrac) + ccnc(mgs)*renucfrac + IF ( irenuc == 2 ) THEN + cnuc(mgs) = Max(ccnc(mgs),cwnccn(mgs))*(1. - renucfrac) + ccnc(mgs)*renucfrac + ELSE + cnuc(mgs) = ccnc(mgs)*(1. - renucfrac) + ccnc(mgs)*renucfrac + ENDIF ELSE cnuc(mgs) = ccnc(mgs)*(1. - renucfrac) + Max(0.0,ccnc(mgs) - ccna(mgs))*renucfrac ENDIF