!+---+-----------------------------------------------------------------+ ! ! RELEASED JANUARY 2025 ! ! This is the WRF module for the RCON microphysics scheme, the intent ! of which is to improve warm rain representation within the Thompson-Eidhammer scheme. ! ! RCON is based heavily on the Thompson-Eidhammer scheme with a couple significant ! changes that improve upon the code in module_mp_rcon.F to generate more realistic ! rainfall during warm rain events with additional benefits for cold rain, especially ! warm processes during cold rain events. ! ! Most deviations from the Thompson-Eidhammer scheme are marked with one of the following comments: ! !RC, !RC !\RC (for blocks) ! ! ! Among the most significant changes for rain productions are 1) the use of a wider cloud water DSD ! of lognormal shape instead of the gamma DSD used by the Thompson–Eidhammer parameterization and ! 2) enhancement of the cloud-to-rain autoconversion parameterization to accomodate the new shape. ! The changes here also allow for sedimentation of cloud water within the lowest model layer, which ! effectively creates a drizzle mode in the scheme. One new diagnostic is provided in wrfout, CLOUDNC, ! which is the accumulated cloud water at the surface, i.e. drizzle. ! ! A future version of this scheme will introduce an expanded aerosol table to increase accuracy of cloud representation. ! ! ! Further details and extensive evaluation of the rain/cloud changes is found in the following paper: ! Conrick, R., C. F. Mass, and L. McMurdie, 2023: Improving Simulations of Warm Rain in a Bulk Microphysics Scheme. ! Mon. Wea. Rev., 152, 169–185, https://doi.org/10.1175/MWR-D-23-0035.1. ! ! ! ! NOTE 1: The intent of this scheme is for the Thompson aerosol options to be used. ! Either the climotological aerosol data or the 'default' profile may be used. ! NOTE 2: The Thompson hail scheme is NOT supported at this time. ! ! ! Questions/comments/bugs? ! Contact Robert Conrick at robert.conrick@gmail.com or Clifford Mass at cmass@uw.edu ! ! ! !+---+-----------------------------------------------------------------+ !wrft:model_layer:physics !+---+-----------------------------------------------------------------+ ! MODULE module_mp_rcon USE module_wrf_error USE module_mp_radar #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) USE module_dm, ONLY : wrf_dm_max_real #endif USE module_model_constants, only : RE_QC_BG, RE_QI_BG, RE_QS_BG IMPLICIT NONE LOGICAL, PARAMETER, PRIVATE:: iiwarm = .false. LOGICAL, PRIVATE:: is_aerosol_aware = .false. LOGICAL, PARAMETER, PRIVATE:: dustyIce = .true. LOGICAL, PARAMETER, PRIVATE:: homogIce = .true. LOGICAL, PRIVATE:: is_hail_aware = .false. INTEGER, PARAMETER, PRIVATE:: IFDRY = 0 REAL, PARAMETER, PRIVATE:: T_0 = 273.15 REAL, PARAMETER, PRIVATE:: PI = 3.1415926536 !..Densities of rain, snow, graupel, and cloud ice. REAL, PARAMETER, PRIVATE:: rho_w = 1000.0 REAL, PRIVATE:: rho_s = 100.0 REAL, PARAMETER, PRIVATE:: rho_i = 890.0 INTEGER, PARAMETER, PRIVATE:: NRHG = 9 INTEGER, PARAMETER, PRIVATE:: NRHG1 = 1 INTEGER :: dimNRHG REAL, DIMENSION(NRHG), PARAMETER, PRIVATE:: & rho_g = (/50., 100., 200., 300., 400., 500., 600., 700., 800./) INTEGER, PARAMETER :: idx_bg1 = 5 ! index for rhog when mp=8 or 28 !..Prescribed number of cloud droplets. Set according to known data or !.. roughly 100 per cc (100.E6 m^-3) for Maritime cases and !.. 300 per cc (300.E6 m^-3) for Continental. Gamma shape parameter, !.. mu_c, calculated based on Nt_c is important in autoconversion !.. scheme. In 2-moment cloud water, Nt_c represents a maximum of !.. droplet concentration and nu_c is also variable depending on local !.. droplet number concentration. REAL, PARAMETER, PRIVATE:: Nt_c = 100.E6 REAL, PARAMETER, PRIVATE:: Nt_c_max = 1999.E6 !..Declaration of constants for assumed CCN/IN aerosols when none in !.. the input data. Look inside the init routine for modifications !.. due to surface land-sea points or vegetation characteristics. REAL, PARAMETER, PRIVATE:: naIN0 = 1.5E6 REAL, PARAMETER, PRIVATE:: naIN1 = 0.5E6 REAL, PARAMETER, PRIVATE:: naCCN0 = 300.0E6 REAL, PARAMETER, PRIVATE:: naCCN1 = 50.0E6 REAL, PARAMETER, PRIVATE:: naBC0 = 150.0E6 REAL, PARAMETER, PRIVATE:: naBC1 = 25.0E6 !..Generalized gamma distributions for rain, graupel and cloud ice. !.. N(D) = N_0 * D**mu * exp(-lamda*D); mu=0 is exponential. REAL, PARAMETER, PRIVATE:: mu_r = 0.0 REAL, PARAMETER, PRIVATE:: mu_g = 0.0 REAL, PARAMETER, PRIVATE:: mu_i = 0.0 REAL, PRIVATE:: mu_c !..Sum of two gamma distrib for snow (Field et al. 2005). !.. N(D) = M2**4/M3**3 * [Kap0*exp(-M2*Lam0*D/M3) !.. + Kap1*(M2/M3)**mu_s * D**mu_s * exp(-M2*Lam1*D/M3)] !.. M2 and M3 are the (bm_s)th and (bm_s+1)th moments respectively !.. calculated as function of ice water content and temperature. REAL, PARAMETER, PRIVATE:: mu_s = 0.6357 REAL, PARAMETER, PRIVATE:: Kap0 = 490.6 REAL, PARAMETER, PRIVATE:: Kap1 = 17.46 REAL, PARAMETER, PRIVATE:: Lam0 = 20.78 REAL, PARAMETER, PRIVATE:: Lam1 = 3.29 !..Y-intercept parameter for graupel is not constant and depends on !.. mixing ratio. Also, when mu_g is non-zero, these become equiv !.. y-intercept for an exponential distrib and proper values are !.. computed based on same mixing ratio and total number concentration. REAL, PARAMETER, PRIVATE:: gonv_min = 1.E2 REAL, PARAMETER, PRIVATE:: gonv_max = 1.E6 !..Mass power law relations: mass = am*D**bm !.. Snow from Field et al. (2005), others assume spherical form. REAL, PARAMETER, PRIVATE:: am_r = PI*rho_w/6.0 REAL, PARAMETER, PRIVATE:: bm_r = 3.0 REAL, PARAMETER, PRIVATE:: am_s = 0.069 REAL, PARAMETER, PRIVATE:: bm_s = 2.0 REAL, DIMENSION(NRHG), PARAMETER,PRIVATE:: am_g = (/ & & PI*rho_g(1)/6.0, PI*rho_g(2)/6.0, PI*rho_g(3)/6.0, & & PI*rho_g(4)/6.0, PI*rho_g(5)/6.0, PI*rho_g(6)/6.0, & & PI*rho_g(7)/6.0, PI*rho_g(8)/6.0, PI*rho_g(9)/6.0 /) REAL, PARAMETER, PRIVATE:: bm_g = 3.0 REAL, PARAMETER, PRIVATE:: am_i = PI*rho_i/6.0 REAL, PARAMETER, PRIVATE:: bm_i = 3.0 !..Fallspeed power laws relations: v = (av*D**bv)*exp(-fv*D) !.. Rain from Ferrier (1994), ice, snow, and graupel from !.. Thompson et al (2008). Coefficient fv is zero for graupel/ice. REAL, PARAMETER, PRIVATE:: av_r = 4854.0 REAL, PARAMETER, PRIVATE:: bv_r = 1.0 REAL, PARAMETER, PRIVATE:: fv_r = 195.0 REAL, PARAMETER, PRIVATE:: av_s = 40.0 REAL, PARAMETER, PRIVATE:: bv_s = 0.55 REAL, PARAMETER, PRIVATE:: fv_s = 100.0 REAL, PARAMETER, PRIVATE:: av_g_old = 442.0 REAL, PARAMETER, PRIVATE:: bv_g_old = 0.89 REAL, DIMENSION(NRHG), PRIVATE:: & ! Computed from A. Heymsfield: Best - Reynolds relation & av_g = (/ 45.9173813, 67.0867386, 98.0158463, 122.353378, & & 143.204224, 161.794724, 178.762115, 194.488785, & & 209.225876/) REAL, DIMENSION(NRHG), PRIVATE:: & ! Computed from A. Heymsfield: Best - Reynolds relation & bv_g = (/0.640961647, 0.640961647, 0.640961647, 0.640961647, & & 0.640961647, 0.640961647, 0.640961647, 0.640961647, & & 0.640961647/) REAL, PARAMETER, PRIVATE:: a_coeff = 0.47244157 REAL, PARAMETER, PRIVATE:: b_coeff = 0.54698726 REAL, PARAMETER, PRIVATE:: av_i = 1493.9 REAL, PARAMETER, PRIVATE:: bv_i = 1.0 REAL, PARAMETER, PRIVATE:: av_c = 0.316946E8 REAL, PARAMETER, PRIVATE:: bv_c = 2.0 !..Capacitance of sphere and plates/aggregates: D**3, D**2 REAL, PARAMETER, PRIVATE:: C_cube = 0.5 REAL, PARAMETER, PRIVATE:: C_sqrd = 0.15 !..Collection efficiencies. Rain/snow/graupel collection of cloud !.. droplets use variables (Ef_rw, Ef_sw, Ef_gw respectively) and !.. get computed elsewhere because they are dependent on stokes !.. number. REAL, PARAMETER, PRIVATE:: Ef_si = 0.05 REAL, PARAMETER, PRIVATE:: Ef_rs = 0.95 REAL, PARAMETER, PRIVATE:: Ef_rg = 0.75 REAL, PARAMETER, PRIVATE:: Ef_ri = 0.95 !..Minimum microphys values !.. R1 value, 1.E-12, cannot be set lower because of numerical !.. problems with Paul Field's moments and should not be set larger !.. because of truncation problems in snow/ice growth. REAL, PARAMETER, PRIVATE:: R1 = 1.E-12 REAL, PARAMETER, PRIVATE:: R2 = 1.E-6 REAL, PARAMETER, PRIVATE:: eps = 1.E-15 !..Constants in Cooper curve relation for cloud ice number. REAL, PARAMETER, PRIVATE:: TNO = 5.0 REAL, PARAMETER, PRIVATE:: ATO = 0.304 !..Rho_not used in fallspeed relations (rho_not/rho)**.5 adjustment. REAL, PARAMETER, PRIVATE:: rho_not = 101325.0/(287.05*298.0) !..Schmidt number REAL, PARAMETER, PRIVATE:: Sc = 0.632 REAL, PRIVATE:: Sc3 !..Homogeneous freezing temperature REAL, PARAMETER, PRIVATE:: HGFR = 235.16 !..Water vapor and air gas constants at constant pressure REAL, PARAMETER, PRIVATE:: Rv = 461.5 REAL, PARAMETER, PRIVATE:: oRv = 1./Rv REAL, PARAMETER, PRIVATE:: R = 287.04 REAL, PARAMETER, PRIVATE:: Cp = 1004.0 REAL, PARAMETER, PRIVATE:: R_uni = 8.314 ! J (mol K)-1 DOUBLE PRECISION, PARAMETER, PRIVATE:: k_b = 1.38065E-23 ! Boltzmann constant [J/K] DOUBLE PRECISION, PARAMETER, PRIVATE:: M_w = 18.01528E-3 ! molecular mass of water [kg/mol] DOUBLE PRECISION, PARAMETER, PRIVATE:: M_a = 28.96E-3 ! molecular mass of air [kg/mol] DOUBLE PRECISION, PARAMETER, PRIVATE:: N_avo = 6.022E23 ! Avogadro number [1/mol] DOUBLE PRECISION, PARAMETER, PRIVATE:: ma_w = M_w / N_avo ! mass of water molecule [kg] REAL, PARAMETER, PRIVATE:: ar_volume = 4./3.*PI*(2.5e-6)**3 ! assume radius of 0.025 micrometer, 2.5e-6 cm !..Enthalpy of sublimation, vaporization, and fusion at 0C. REAL, PARAMETER, PRIVATE:: lsub = 2.834E6 REAL, PARAMETER, PRIVATE:: lvap0 = 2.5E6 REAL, PARAMETER, PRIVATE:: lfus = lsub - lvap0 REAL, PARAMETER, PRIVATE:: olfus = 1./lfus !..Ice initiates with this mass (kg), corresponding diameter calc. !..Min diameters and mass of cloud, rain, snow, and graupel (m, kg). REAL, PARAMETER, PRIVATE:: xm0i = 1.E-12 REAL, PARAMETER, PRIVATE:: D0c = 1.E-6 REAL, PARAMETER, PRIVATE:: D0r = 50.E-6 REAL, PARAMETER, PRIVATE:: D0s = 300.E-6 REAL, PARAMETER, PRIVATE:: D0g = 350.E-6 REAL, PRIVATE:: D0i, xm0s, xm0g !..Lookup table dimensions INTEGER, PARAMETER, PRIVATE:: nbins = 100 INTEGER, PARAMETER, PRIVATE:: nbc = nbins INTEGER, PARAMETER, PRIVATE:: nbi = nbins INTEGER, PARAMETER, PRIVATE:: nbr = nbins INTEGER, PARAMETER, PRIVATE:: nbs = nbins INTEGER, PARAMETER, PRIVATE:: nbg = nbins INTEGER, PARAMETER, PRIVATE:: ntb_c = 37 INTEGER, PARAMETER, PRIVATE:: ntb_i = 64 INTEGER, PARAMETER, PRIVATE:: ntb_r = 37 INTEGER, PARAMETER, PRIVATE:: ntb_s = 37 INTEGER, PARAMETER, PRIVATE:: ntb_g = 37 INTEGER, PARAMETER, PRIVATE:: ntb_g1 = 37 INTEGER, PARAMETER, PRIVATE:: ntb_r1 = 37 INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55 INTEGER, PARAMETER, PRIVATE:: ntb_t = 9 INTEGER, PRIVATE:: nic1, nic2, nii2, nii3, nir2, nir3, nis2, nig2, nig3 INTEGER, PARAMETER, PRIVATE:: ntb_arc = 7 INTEGER, PARAMETER, PRIVATE:: ntb_arw = 9 INTEGER, PARAMETER, PRIVATE:: ntb_art = 7 INTEGER, PARAMETER, PRIVATE:: ntb_arr = 5 INTEGER, PARAMETER, PRIVATE:: ntb_ark = 4 INTEGER, PARAMETER, PRIVATE:: ntb_IN = 55 INTEGER, PRIVATE:: niIN2 DOUBLE PRECISION, DIMENSION(nbins+1):: xDx DOUBLE PRECISION, DIMENSION(nbc):: Dc, dtc DOUBLE PRECISION, DIMENSION(nbi):: Di, dti DOUBLE PRECISION, DIMENSION(nbr):: Dr, dtr DOUBLE PRECISION, DIMENSION(nbs):: Ds, dts DOUBLE PRECISION, DIMENSION(nbg):: Dg, dtg DOUBLE PRECISION, DIMENSION(nbc):: t_Nc !..Lookup tables for cloud water content (kg/m**3). REAL, DIMENSION(ntb_c), PARAMETER, PRIVATE:: & r_c = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, & 1.e-2/) !..Lookup tables for cloud ice content (kg/m**3). REAL, DIMENSION(ntb_i), PARAMETER, PRIVATE:: & r_i = (/1.e-10,2.e-10,3.e-10,4.e-10, & 5.e-10,6.e-10,7.e-10,8.e-10,9.e-10, & 1.e-9,2.e-9,3.e-9,4.e-9,5.e-9,6.e-9,7.e-9,8.e-9,9.e-9, & 1.e-8,2.e-8,3.e-8,4.e-8,5.e-8,6.e-8,7.e-8,8.e-8,9.e-8, & 1.e-7,2.e-7,3.e-7,4.e-7,5.e-7,6.e-7,7.e-7,8.e-7,9.e-7, & 1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & 1.e-3/) !..Lookup tables for rain content (kg/m**3). REAL, DIMENSION(ntb_r), PARAMETER, PRIVATE:: & r_r = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, & 1.e-2/) !..Lookup tables for graupel content (kg/m**3). REAL, DIMENSION(ntb_g), PARAMETER, PRIVATE:: & r_g = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, & 1.e-2/) !..Lookup tables for snow content (kg/m**3). REAL, DIMENSION(ntb_s), PARAMETER, PRIVATE:: & r_s = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, & 1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, & 1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, & 1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, & 1.e-2/) !..Lookup tables for rain y-intercept parameter (/m**4). REAL, DIMENSION(ntb_r1), PARAMETER, PRIVATE:: & N0r_exp = (/1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, & 1.e7,2.e7,3.e7,4.e7,5.e7,6.e7,7.e7,8.e7,9.e7, & 1.e8,2.e8,3.e8,4.e8,5.e8,6.e8,7.e8,8.e8,9.e8, & 1.e9,2.e9,3.e9,4.e9,5.e9,6.e9,7.e9,8.e9,9.e9, & 1.e10/) !..Lookup tables for graupel y-intercept parameter (/m**4). REAL, DIMENSION(ntb_g1), PARAMETER, PRIVATE:: & N0g_exp = (/1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, & 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, & 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, & 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, & 1.e6/) !..Lookup tables for ice number concentration (/m**3). REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: & Nt_i = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, & 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, & 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, & 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, & 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, & 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, & 1.e6/) !..Aerosol table parameter: Number of available aerosols, vertical !.. velocity, temperature, aerosol mean radius, and hygroscopicity. REAL, DIMENSION(ntb_arc), PARAMETER, PRIVATE:: & ta_Na = (/10.0, 31.6, 100.0, 316.0, 1000.0, 3160.0, 10000.0/) REAL, DIMENSION(ntb_arw), PARAMETER, PRIVATE:: & ta_Ww = (/0.01, 0.0316, 0.1, 0.316, 1.0, 3.16, 10.0, 31.6, 100.0/) REAL, DIMENSION(ntb_art), PARAMETER, PRIVATE:: & ta_Tk = (/243.15, 253.15, 263.15, 273.15, 283.15, 293.15, 303.15/) REAL, DIMENSION(ntb_arr), PARAMETER, PRIVATE:: & ta_Ra = (/0.01, 0.02, 0.04, 0.08, 0.16/) REAL, DIMENSION(ntb_ark), PARAMETER, PRIVATE:: & ta_Ka = (/0.2, 0.4, 0.6, 0.8/) !..Lookup tables for IN concentration (/m**3) from 0.001 to 1000/Liter. REAL, DIMENSION(ntb_IN), PARAMETER, PRIVATE:: & Nt_IN = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, & 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, & 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, & 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, & 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, & 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, & 1.e6/) !..For snow moments conversions (from Field et al. 2005) REAL, DIMENSION(10), PARAMETER, PRIVATE:: & sa = (/ 5.065339, -0.062659, -3.032362, 0.029469, -0.000285, & 0.31255, 0.000204, 0.003199, 0.0, -0.015952/) REAL, DIMENSION(10), PARAMETER, PRIVATE:: & sb = (/ 0.476221, -0.015896, 0.165977, 0.007468, -0.000141, & 0.060366, 0.000079, 0.000594, 0.0, -0.003577/) !..Temperatures (5 C interval 0 to -40) used in lookup tables. REAL, DIMENSION(ntb_t), PARAMETER, PRIVATE:: & Tc = (/-0.01, -5., -10., -15., -20., -25., -30., -35., -40./) !..Lookup tables for various accretion/collection terms. !.. ntb_x refers to the number of elements for rain, snow, graupel, !.. and temperature array indices. Variables beginning with t-p/c/m/n !.. represent lookup tables. Save compile-time memory by making !.. allocatable (2009Jun12, J. Michalakes). INTEGER, PARAMETER, PRIVATE:: R8SIZE = 8 INTEGER, PARAMETER, PRIVATE:: R4SIZE = 4 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:,:):: & tcg_racg, tmr_racg, tcr_gacr, & ! tmg_gacr tnr_racg, tnr_gacr REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: & tcs_racs1, tmr_racs1, tcs_racs2, tmr_racs2, & tcr_sacr1, tms_sacr1, tcr_sacr2, tms_sacr2, & tnr_racs1, tnr_racs2, tnr_sacr1, tnr_sacr2 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: & tpi_qcfz, tni_qcfz REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: & tpi_qrfz, tpg_qrfz, tni_qrfz, tnr_qrfz REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: & tps_iaus, tni_iaus, tpi_ide REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efrw REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efsw REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: tnr_rev REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: & tpc_wev, tnc_wev REAL (KIND=R4SIZE), ALLOCATABLE, DIMENSION(:,:,:,:,:):: tnccn_act !..Variables holding a bunch of exponents and gamma values (cloud water, !.. cloud ice, rain, snow, then graupel). REAL, DIMENSION(5,15), PRIVATE:: cce, ccg REAL, DIMENSION(15), PRIVATE:: ocg1, ocg2 REAL, DIMENSION(7), PRIVATE:: cie, cig REAL, PRIVATE:: oig1, oig2, obmi REAL, DIMENSION(13), PRIVATE:: cre, crg REAL, PRIVATE:: ore1, org1, org2, org3, obmr REAL, DIMENSION(17), PRIVATE:: cse, csg REAL, PRIVATE:: oams, obms, ocms REAL, DIMENSION(12,NRHG), PRIVATE:: cge, cgg REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, obmg REAL, DIMENSION(NRHG), PRIVATE:: oamg, ocmg !..Declaration of precomputed constants in various rate eqns. REAL:: t1_qr_qc, t1_qr_qi, t2_qr_qi, t1_qg_qc, t1_qs_qc, t1_qs_qi REAL:: t1_qr_ev, t2_qr_ev REAL:: t1_qs_sd, t2_qs_sd, t1_qg_sd, t2_qg_sd REAL:: t1_qs_me, t2_qs_me, t1_qg_me, t2_qg_me !+---+ !+---+-----------------------------------------------------------------+ !..END DECLARATIONS !+---+-----------------------------------------------------------------+ !+---+ !ctrlL CONTAINS SUBROUTINE rcon_init(hgt, orho, nwfa2d, nbca2d, ng, & nwfa, nifa, nbca, wif_input_opt, frc_urb2d, & dx, dy, is_start, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte) IMPLICIT NONE INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte REAL, DIMENSION(ims:ime,kms:kme,jms:jme), INTENT(IN):: hgt !..OPTIONAL variables that control application of hail-aware scheme REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(INOUT) :: ng !..OPTIONAL variables that control application of aerosol-aware scheme REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(INOUT) :: nwfa, nifa, nbca REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(INOUT) :: nwfa2d, nbca2d REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(IN) :: orho REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(IN) :: frc_urb2d REAL, OPTIONAL, INTENT(IN) :: DX, DY LOGICAL, OPTIONAL, INTENT(IN) :: is_start INTEGER, OPTIONAL, INTENT(IN) :: wif_input_opt CHARACTER*256:: mp_debug INTEGER:: i, j, k, l, m, n REAL:: h_01, niIN3, niCCN3, niBC3, max_test, z1 LOGICAL:: micro_init, has_CCN, has_IN !+---+ if (PRESENT(ng)) then is_hail_aware = .TRUE. dimNRHG = NRHG else av_g(idx_bg1) = av_g_old bv_g(idx_bg1) = bv_g_old dimNRHG = NRHG1 endif is_aerosol_aware = .FALSE. micro_init = .FALSE. has_CCN = .FALSE. has_IN = .FALSE. write(mp_debug,*) ' DEBUG checking column of hgt ', its+1,jts+1 CALL wrf_debug(250, mp_debug) do k = kts, kte write(mp_debug,*) ' DEBUGT k, hgt = ', k, hgt(its+1,k,jts+1) CALL wrf_debug(250, mp_debug) enddo if (PRESENT(nwfa2d) .AND. PRESENT(nwfa) .AND. PRESENT(nifa)) is_aerosol_aware = .TRUE. if (is_aerosol_aware) then !..Check for existing aerosol data, both CCN and IN aerosols. If missing !.. fill in just a basic vertical profile, somewhat boundary-layer following. #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) max_test = wrf_dm_max_real ( MAXVAL(nwfa(its:ite-1,:,jts:jte-1)) ) #else max_test = MAXVAL ( nwfa(its:ite-1,:,jts:jte-1) ) #endif if (max_test .lt. eps) then write(mp_debug,*) ' Apparently there are no initial CCN aerosols.' CALL wrf_debug(100, mp_debug) write(mp_debug,*) ' checked column at point (i,j) = ', its,jts CALL wrf_debug(100, mp_debug) do j = jts, min(jde-1,jte) do i = its, min(ide-1,ite) if (hgt(i,1,j).le.1000.0) then h_01 = 0.8 elseif (hgt(i,1,j).ge.2500.0) then h_01 = 0.01 else h_01 = 0.8*cos(hgt(i,1,j)*0.001 - 1.0) endif niCCN3 = -1.0*ALOG(naCCN1/naCCN0)/h_01 nwfa(i,1,j) = naCCN1+naCCN0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niCCN3) z1=hgt(i,2,j)-hgt(i,1,j) nwfa2d(i,j) = nwfa(i,1,j) * 0.000196 * (50./z1) do k = 2, kte nwfa(i,k,j) = naCCN1+naCCN0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niCCN3) enddo enddo enddo else has_CCN = .TRUE. write(mp_debug,*) ' Apparently initial CCN aerosols are present.' CALL wrf_debug(100, mp_debug) write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nwfa(its,:,jts)) CALL wrf_debug(100, mp_debug) endif #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) max_test = wrf_dm_max_real ( MAXVAL(nifa(its:ite-1,:,jts:jte-1)) ) #else max_test = MAXVAL ( nifa(its:ite-1,:,jts:jte-1) ) #endif if (max_test .lt. eps) then write(mp_debug,*) ' Apparently there are no initial IN aerosols.' CALL wrf_debug(100, mp_debug) write(mp_debug,*) ' checked column at point (i,j) = ', its,jts CALL wrf_debug(100, mp_debug) do j = jts, min(jde-1,jte) do i = its, min(ide-1,ite) if (hgt(i,1,j).le.1000.0) then h_01 = 0.8 elseif (hgt(i,1,j).ge.2500.0) then h_01 = 0.01 else h_01 = 0.8*cos(hgt(i,1,j)*0.001 - 1.0) endif niIN3 = -1.0*ALOG(naIN1/naIN0)/h_01 nifa(i,1,j) = naIN1+naIN0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niIN3) do k = 2, kte nifa(i,k,j) = naIN1+naIN0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niIN3) enddo enddo enddo else has_IN = .TRUE. write(mp_debug,*) ' Apparently initial IN aerosols are present.' CALL wrf_debug(100, mp_debug) write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nifa(its,:,jts)) CALL wrf_debug(100, mp_debug) endif if ( wif_input_opt .eq. 2 ) then #if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) max_test = wrf_dm_max_real ( MAXVAL(nbca(its:ite-1,:,jts:jte-1)) ) #else max_test = MAXVAL ( nbca(its:ite-1,:,jts:jte-1) ) #endif if (max_test .lt. eps) then write(mp_debug,*) ' Apparently there are no initial BC aerosols.' CALL wrf_debug(100, mp_debug) write(mp_debug,*) ' checked column at point (i,j) = ', its,jts CALL wrf_debug(100, mp_debug) do j = jts, min(jde-1,jte) do i = its, min(ide-1,ite) if (hgt(i,1,j).le.1000.0) then h_01 = 0.8 elseif (hgt(i,1,j).ge.2500.0) then h_01 = 0.01 else h_01 = 0.8*cos(hgt(i,1,j)*0.001 - 1.0) endif niBC3 = -1.0*ALOG(naBC1/naBC0)/h_01 nbca(i,1,j) = naBC1+naBC0*exp(-((hgt(i,2,j)-hgt(i,1,j))/1000.)*niBC3) z1=hgt(i,2,j)-hgt(i,1,j) nbca2d(i,j) = nbca(i,1,j) * 0.000098 * (50./z1) * (1. + frc_urb2d(i,j)) do k = 2, kte nbca(i,k,j) = naBC1+naBC0*exp(-((hgt(i,k,j)-hgt(i,1,j))/1000.)*niBC3) enddo enddo enddo else write(mp_debug,*) ' Apparently initial BC aerosols are present.' CALL wrf_debug(100, mp_debug) write(mp_debug,*) ' column sum at point (i,j) = ', its,jts, SUM(nbca(its,:,jts)) CALL wrf_debug(100, mp_debug) endif endif endif !..Allocate space for lookup tables (J. Michalakes 2009Jun08). if (.NOT. ALLOCATED(tcg_racg) ) then ALLOCATE(tcg_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) micro_init = .TRUE. endif if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tcr_gacr)) ALLOCATE(tcr_gacr(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) ! if (.NOT. ALLOCATED(tmg_gacr)) ALLOCATE(tmg_gacr(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tnr_racg)) ALLOCATE(tnr_racg(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tnr_gacr)) ALLOCATE(tnr_gacr(ntb_g1,ntb_g,dimNRHG,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tcs_racs1)) ALLOCATE(tcs_racs1(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tmr_racs1)) ALLOCATE(tmr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tcs_racs2)) ALLOCATE(tcs_racs2(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tmr_racs2)) ALLOCATE(tmr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tcr_sacr1)) ALLOCATE(tcr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tms_sacr1)) ALLOCATE(tms_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tcr_sacr2)) ALLOCATE(tcr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tms_sacr2)) ALLOCATE(tms_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tnr_racs1)) ALLOCATE(tnr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tnr_racs2)) ALLOCATE(tnr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tnr_sacr1)) ALLOCATE(tnr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tnr_sacr2)) ALLOCATE(tnr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r)) if (.NOT. ALLOCATED(tpi_qcfz)) ALLOCATE(tpi_qcfz(ntb_c,nbc,45,ntb_IN)) if (.NOT. ALLOCATED(tni_qcfz)) ALLOCATE(tni_qcfz(ntb_c,nbc,45,ntb_IN)) if (.NOT. ALLOCATED(tpi_qrfz)) ALLOCATE(tpi_qrfz(ntb_r,ntb_r1,45,ntb_IN)) if (.NOT. ALLOCATED(tpg_qrfz)) ALLOCATE(tpg_qrfz(ntb_r,ntb_r1,45,ntb_IN)) if (.NOT. ALLOCATED(tni_qrfz)) ALLOCATE(tni_qrfz(ntb_r,ntb_r1,45,ntb_IN)) if (.NOT. ALLOCATED(tnr_qrfz)) ALLOCATE(tnr_qrfz(ntb_r,ntb_r1,45,ntb_IN)) if (.NOT. ALLOCATED(tps_iaus)) ALLOCATE(tps_iaus(ntb_i,ntb_i1)) if (.NOT. ALLOCATED(tni_iaus)) ALLOCATE(tni_iaus(ntb_i,ntb_i1)) if (.NOT. ALLOCATED(tpi_ide)) ALLOCATE(tpi_ide(ntb_i,ntb_i1)) if (.NOT. ALLOCATED(t_Efrw)) ALLOCATE(t_Efrw(nbr,nbc)) if (.NOT. ALLOCATED(t_Efsw)) ALLOCATE(t_Efsw(nbs,nbc)) if (.NOT. ALLOCATED(tnr_rev)) ALLOCATE(tnr_rev(nbr, ntb_r1, ntb_r)) if (.NOT. ALLOCATED(tpc_wev)) ALLOCATE(tpc_wev(nbc,ntb_c,nbc)) if (.NOT. ALLOCATED(tnc_wev)) ALLOCATE(tnc_wev(nbc,ntb_c,nbc)) if (.NOT. ALLOCATED(tnccn_act)) & ALLOCATE(tnccn_act(ntb_arc,ntb_arw,ntb_art,ntb_arr,ntb_ark)) if (micro_init) then !..From Martin et al. (1994), assign gamma shape parameter mu for cloud !.. drops according to general dispersion characteristics (disp=~0.25 !.. for Maritime and 0.45 for Continental). !.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime !.. to 2 for really dirty air. This not used in 2-moment cloud water !.. scheme and nu_c used instead and varies from 2 to 15 (integer-only). mu_c = MIN(15., (1000.E6/Nt_c + 2.)) !..Schmidt number to one-third used numerous times. Sc3 = Sc**(1./3.) !..Compute min ice diam from mass, min snow/graupel mass from diam. D0i = (xm0i/am_i)**(1./bm_i) xm0s = am_s * D0s**bm_s xm0g = am_g(NRHG) * D0g**bm_g !..These constants various exponents and gamma() assoc with cloud, !.. rain, snow, and graupel. do n = 1, 15 cce(1,n) = n + 1. cce(2,n) = bm_r + n + 1. cce(3,n) = bm_r + n + 4. cce(4,n) = n + bv_c + 1. cce(5,n) = bm_r + n + bv_c + 1. ccg(1,n) = WGAMMA(cce(1,n)) ccg(2,n) = WGAMMA(cce(2,n)) ccg(3,n) = WGAMMA(cce(3,n)) ccg(4,n) = WGAMMA(cce(4,n)) ccg(5,n) = WGAMMA(cce(5,n)) ocg1(n) = 1./ccg(1,n) ocg2(n) = 1./ccg(2,n) enddo cie(1) = mu_i + 1. cie(2) = bm_i + mu_i + 1. cie(3) = bm_i + mu_i + bv_i + 1. cie(4) = mu_i + bv_i + 1. cie(5) = mu_i + 2. cie(6) = bm_i*0.5 + mu_i + bv_i + 1. cie(7) = bm_i*0.5 + mu_i + 1. cig(1) = WGAMMA(cie(1)) cig(2) = WGAMMA(cie(2)) cig(3) = WGAMMA(cie(3)) cig(4) = WGAMMA(cie(4)) cig(5) = WGAMMA(cie(5)) cig(6) = WGAMMA(cie(6)) cig(7) = WGAMMA(cie(7)) oig1 = 1./cig(1) oig2 = 1./cig(2) obmi = 1./bm_i cre(1) = bm_r + 1. cre(2) = mu_r + 1. cre(3) = bm_r + mu_r + 1. cre(4) = bm_r*2. + mu_r + 1. cre(5) = mu_r + bv_r + 1. cre(6) = bm_r + mu_r + bv_r + 1. cre(7) = bm_r*0.5 + mu_r + bv_r + 1. cre(8) = bm_r + mu_r + bv_r + 3. cre(9) = mu_r + bv_r + 3. cre(10) = mu_r + 2. cre(11) = 0.5*(bv_r + 5. + 2.*mu_r) cre(12) = bm_r*0.5 + mu_r + 1. cre(13) = bm_r*2. + mu_r + bv_r + 1. do n = 1, 13 crg(n) = WGAMMA(cre(n)) enddo obmr = 1./bm_r ore1 = 1./cre(1) org1 = 1./crg(1) org2 = 1./crg(2) org3 = 1./crg(3) cse(1) = bm_s + 1. cse(2) = bm_s + 2. cse(3) = bm_s*2. cse(4) = bm_s + bv_s + 1. cse(5) = bm_s*2. + bv_s + 1. cse(6) = bm_s*2. + 1. cse(7) = bm_s + mu_s + 1. cse(8) = bm_s + mu_s + 2. cse(9) = bm_s + mu_s + 3. cse(10) = bm_s + mu_s + bv_s + 1. cse(11) = bm_s*2. + mu_s + bv_s + 1. cse(12) = bm_s*2. + mu_s + 1. cse(13) = bv_s + 2. cse(14) = bm_s + bv_s cse(15) = mu_s + 1. cse(16) = 1.0 + (1.0 + bv_s)/2. cse(17) = bm_s + bv_s + 2. do n = 1, 17 csg(n) = WGAMMA(cse(n)) enddo oams = 1./am_s obms = 1./bm_s ocms = oams**obms cge(1,:) = bm_g + 1. cge(2,:) = mu_g + 1. cge(3,:) = bm_g + mu_g + 1. cge(4,:) = bm_g*2. + mu_g + 1. cge(10,:) = mu_g + 2. cge(12,:) = bm_g*0.5 + mu_g + 1. do m = 1, NRHG cge(5,m) = bm_g*2. + mu_g + bv_g(m) + 1. cge(6,m) = bm_g + mu_g + bv_g(m) + 1. cge(7,m) = bm_g*0.5 + mu_g + bv_g(m) + 1. cge(8,m) = mu_g + bv_g(m) + 1. ! not used cge(9,m) = mu_g + bv_g(m) + 3. cge(11,m) = 0.5*(bv_g(m) + 5. + 2.*mu_g) enddo do m = 1, NRHG do n = 1, 12 cgg(n,m) = WGAMMA(cge(n,m)) enddo enddo oamg = 1./am_g obmg = 1./bm_g do m = 1, NRHG oamg(m) = 1./am_g(m) ocmg(m) = oamg(m)**obmg enddo oge1 = 1./cge(1,1) ogg1 = 1./cgg(1,1) ogg2 = 1./cgg(2,1) ogg3 = 1./cgg(3,1) !+---+-----------------------------------------------------------------+ !..Simplify various rate eqns the best we can now. !+---+-----------------------------------------------------------------+ !..Rain collecting cloud water and cloud ice t1_qr_qc = PI*.25*av_r * crg(9) t1_qr_qi = PI*.25*av_r * crg(9) t2_qr_qi = PI*.25*am_r*av_r * crg(8) !..Graupel collecting cloud water ! t1_qg_qc = PI*.25*av_g * cgg(9) !..Snow collecting cloud water t1_qs_qc = PI*.25*av_s !..Snow collecting cloud ice t1_qs_qi = PI*.25*av_s !..Evaporation of rain; ignore depositional growth of rain. t1_qr_ev = 0.78 * crg(10) t2_qr_ev = 0.308*Sc3*SQRT(av_r) * crg(11) !..Sublimation/depositional growth of snow t1_qs_sd = 0.86 t2_qs_sd = 0.28*Sc3*SQRT(av_s) !..Melting of snow t1_qs_me = PI*4.*C_sqrd*olfus * 0.86 t2_qs_me = PI*4.*C_sqrd*olfus * 0.28*Sc3*SQRT(av_s) !..Sublimation/depositional growth of graupel t1_qg_sd = 0.86 * cgg(10,1) ! t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11) !..Melting of graupel t1_qg_me = PI*4.*C_cube*olfus * 0.86 * cgg(10,1) ! t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11) !..Constants for helping find lookup table indexes. nic2 = NINT(ALOG10(r_c(1))) nii2 = NINT(ALOG10(r_i(1))) nii3 = NINT(ALOG10(Nt_i(1))) nir2 = NINT(ALOG10(r_r(1))) nir3 = NINT(ALOG10(N0r_exp(1))) nis2 = NINT(ALOG10(r_s(1))) nig2 = NINT(ALOG10(r_g(1))) nig3 = NINT(ALOG10(N0g_exp(1))) niIN2 = NINT(ALOG10(Nt_IN(1))) !..Create bins of cloud water (from min diameter up to 100 microns). Dc(1) = D0c*1.0d0 dtc(1) = D0c*1.0d0 do n = 2, nbc Dc(n) = Dc(n-1) + 1.0D-6 dtc(n) = (Dc(n) - Dc(n-1)) enddo !..Create bins of cloud ice (from min diameter up to 5x min snow size). xDx(1) = D0i*1.0d0 xDx(nbi+1) = 2.0d0*D0s do n = 2, nbi xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbi) & *DLOG(xDx(nbi+1)/xDx(1)) +DLOG(xDx(1))) enddo do n = 1, nbi Di(n) = DSQRT(xDx(n)*xDx(n+1)) dti(n) = xDx(n+1) - xDx(n) enddo !..Create bins of rain (from min diameter up to 5 mm). xDx(1) = D0r*1.0d0 xDx(nbr+1) = 0.005d0 do n = 2, nbr xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbr) & *DLOG(xDx(nbr+1)/xDx(1)) +DLOG(xDx(1))) enddo do n = 1, nbr Dr(n) = DSQRT(xDx(n)*xDx(n+1)) dtr(n) = xDx(n+1) - xDx(n) enddo !..Create bins of snow (from min diameter up to 2 cm). xDx(1) = D0s*1.0d0 xDx(nbs+1) = 0.02d0 do n = 2, nbs xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbs) & *DLOG(xDx(nbs+1)/xDx(1)) +DLOG(xDx(1))) enddo do n = 1, nbs Ds(n) = DSQRT(xDx(n)*xDx(n+1)) dts(n) = xDx(n+1) - xDx(n) enddo !..Create bins of graupel (from min diameter up to 5 cm). xDx(1) = D0g*1.0d0 xDx(nbg+1) = 0.05d0 do n = 2, nbg xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbg) & *DLOG(xDx(nbg+1)/xDx(1)) +DLOG(xDx(1))) enddo do n = 1, nbg Dg(n) = DSQRT(xDx(n)*xDx(n+1)) dtg(n) = xDx(n+1) - xDx(n) enddo !..Create bins of cloud droplet number concentration (1 to 3000 per cc). xDx(1) = 1.0d0 xDx(nbc+1) = 3000.0d0 do n = 2, nbc xDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nbc) & *DLOG(xDx(nbc+1)/xDx(1)) +DLOG(xDx(1))) enddo do n = 1, nbc t_Nc(n) = DSQRT(xDx(n)*xDx(n+1)) * 1.D6 enddo nic1 = DLOG(t_Nc(nbc)/t_Nc(1)) !+---+-----------------------------------------------------------------+ !..Create lookup tables for most costly calculations. !+---+-----------------------------------------------------------------+ do m = 1, ntb_r do k = 1, ntb_r1 do n = 1, dimNRHG do j = 1, ntb_g do i = 1, ntb_g1 tcg_racg(i,j,n,k,m) = 0.0d0 tmr_racg(i,j,n,k,m) = 0.0d0 tcr_gacr(i,j,n,k,m) = 0.0d0 !tmg_gacr(i,j,n,k,m) = 0.0d0 tnr_racg(i,j,n,k,m) = 0.0d0 tnr_gacr(i,j,n,k,m) = 0.0d0 enddo enddo enddo enddo enddo do m = 1, ntb_r do k = 1, ntb_r1 do j = 1, ntb_t do i = 1, ntb_s tcs_racs1(i,j,k,m) = 0.0d0 tmr_racs1(i,j,k,m) = 0.0d0 tcs_racs2(i,j,k,m) = 0.0d0 tmr_racs2(i,j,k,m) = 0.0d0 tcr_sacr1(i,j,k,m) = 0.0d0 tms_sacr1(i,j,k,m) = 0.0d0 tcr_sacr2(i,j,k,m) = 0.0d0 tms_sacr2(i,j,k,m) = 0.0d0 tnr_racs1(i,j,k,m) = 0.0d0 tnr_racs2(i,j,k,m) = 0.0d0 tnr_sacr1(i,j,k,m) = 0.0d0 tnr_sacr2(i,j,k,m) = 0.0d0 enddo enddo enddo enddo do m = 1, ntb_IN do k = 1, 45 do j = 1, ntb_r1 do i = 1, ntb_r tpi_qrfz(i,j,k,m) = 0.0d0 tni_qrfz(i,j,k,m) = 0.0d0 tpg_qrfz(i,j,k,m) = 0.0d0 tnr_qrfz(i,j,k,m) = 0.0d0 enddo enddo do j = 1, nbc do i = 1, ntb_c tpi_qcfz(i,j,k,m) = 0.0d0 tni_qcfz(i,j,k,m) = 0.0d0 enddo enddo enddo enddo do j = 1, ntb_i1 do i = 1, ntb_i tps_iaus(i,j) = 0.0d0 tni_iaus(i,j) = 0.0d0 tpi_ide(i,j) = 0.0d0 enddo enddo do j = 1, nbc do i = 1, nbr t_Efrw(i,j) = 0.0 enddo do i = 1, nbs t_Efsw(i,j) = 0.0 enddo enddo do k = 1, ntb_r do j = 1, ntb_r1 do i = 1, nbr tnr_rev(i,j,k) = 0.0d0 enddo enddo enddo do k = 1, nbc do j = 1, ntb_c do i = 1, nbc tpc_wev(i,j,k) = 0.0d0 tnc_wev(i,j,k) = 0.0d0 enddo enddo enddo do m = 1, ntb_ark do l = 1, ntb_arr do k = 1, ntb_art do j = 1, ntb_arw do i = 1, ntb_arc tnccn_act(i,j,k,l,m) = 1.0 enddo enddo enddo enddo enddo CALL wrf_debug(150, 'CREATING MICROPHYSICS LOOKUP TABLES ... ') WRITE (wrf_err_message, '(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') & ' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g CALL wrf_debug(150, wrf_err_message) !..Read a static file containing CCN activation of aerosols. The !.. data were created from a parcel model by Feingold & Heymsfield with !.. further changes by Eidhammer and Kriedenweis. if (is_aerosol_aware) then CALL wrf_debug(200, ' calling table_ccnAct routine') call table_ccnAct endif !..Collision efficiency between rain/snow and cloud water. CALL wrf_debug(200, ' creating qc collision eff tables') call table_Efrw call table_Efsw !..Drop evaporation. CALL wrf_debug(200, ' creating rain evap table') call table_dropEvap !..Initialize various constants for computing radar reflectivity. xam_r = am_r xbm_r = bm_r xmu_r = mu_r xam_s = am_s xbm_s = bm_s xmu_s = mu_s xam_g = am_g(idx_bg1) xbm_g = bm_g xmu_g = mu_g call radar_init if (.not. iiwarm) then !..Rain collecting graupel & graupel collecting rain. CALL wrf_debug(200, ' creating rain collecting graupel table') call qr_acr_qg(dimNRHG) !..Rain collecting snow & snow collecting rain. CALL wrf_debug(200, ' creating rain collecting snow table') call qr_acr_qs !..Cloud water and rain freezing (Bigg, 1953). CALL wrf_debug(200, ' creating freezing of water drops table') call freezeH2O !..Conversion of some ice mass into snow category. CALL wrf_debug(200, ' creating ice converting to snow table') call qi_aut_qs endif CALL wrf_debug(150, ' ... DONE microphysical lookup tables') endif END SUBROUTINE rcon_init !+---+-----------------------------------------------------------------+ !ctrlL !+---+-----------------------------------------------------------------+ !..This is a wrapper routine designed to transfer values from 3D to 1D. !+---+-----------------------------------------------------------------+ SUBROUTINE mp_rcon_driver(qv, qc, qr, qi, qs, qg, qb, ni, nr, nc, ng,& nwfa, nifa, nbca, nwfa2d, nifa2d, nbca2d, & aer_init_opt, wif_input_opt, & th, pii, p, w, dz, & dt_in, itimestep, & RAINNC, RAINNCV, & CLOUDNC, & !RC SNOWNC, SNOWNCV, & GRAUPELNC, GRAUPELNCV, SR, & #if ( WRF_CHEM == 1 ) wetscav_on, rainprod, evapprod, & #endif refl_10cm, diagflag, ke_diag, do_radar_ref, & re_cloud, re_ice, re_snow, & has_reqc, has_reqi, has_reqs, & ids,ide, jds,jde, kds,kde, & ! domain dims ims,ime, jms,jme, kms,kme, & ! memory dims its,ite, jts,jte, kts,kte & ! tile dims ) !RC implicit none !..Subroutine arguments INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & qv, qc, qr, qi, qs, qg, ni, nr, th REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & nc, nwfa, nifa, nbca, qb, ng REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d, & nbca2d REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & re_cloud, re_ice, re_snow INTEGER, INTENT(IN):: has_reqc, has_reqi, has_reqs #if ( WRF_CHEM == 1 ) REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & rainprod, evapprod #endif !RC REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: CLOUDNC !\RC REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: & pii, p, w, dz REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: & RAINNC, RAINNCV, SR REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(INOUT):: & SNOWNC, SNOWNCV, GRAUPELNC, GRAUPELNCV REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & refl_10cm REAL, INTENT(IN):: dt_in INTEGER, INTENT(IN):: itimestep INTEGER, OPTIONAL, INTENT(IN):: aer_init_opt, wif_input_opt #if ( WRF_CHEM == 1 ) LOGICAL, INTENT(in) :: wetscav_on #endif !..Local variables REAL, DIMENSION(kts:kte):: & qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, & ni1d, nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d,& t1d, p1d, w1d, dz1d, rho, dBZ REAL, DIMENSION(kts:kte):: re_qc1d, re_qi1d, re_qs1d #if ( WRF_CHEM == 1 ) REAL, DIMENSION(kts:kte):: & rainprod1d, evapprod1d #endif REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic REAL:: dt, pptrain, pptsnow, pptgraul, pptice, pptcloud !RC; added pptcloud REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max REAL:: nwfa1 REAL:: ygra1, zans1 DOUBLE PRECISION:: lamg, lam_exp, lamr, N0_min, N0_exp INTEGER:: i, j, k, k_0, k_inj INTEGER:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr INTEGER:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr INTEGER:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr INTEGER:: i_start, j_start, i_end, j_end LOGICAL, OPTIONAL, INTENT(IN) :: diagflag INTEGER, OPTIONAL, INTENT(IN) :: do_radar_ref, ke_diag CHARACTER*256:: mp_debug integer :: kediagloc !+---+ i_start = its j_start = jts i_end = MIN(ite, ide-1) j_end = MIN(jte, jde-1) !..For idealized testing by developer. ! if ( (ide-ids+1).gt.4 .and. (jde-jds+1).lt.4 .and. & ! ids.eq.its.and.ide.eq.ite.and.jds.eq.jts.and.jde.eq.jte) then ! i_start = its + 2 ! i_end = ite ! j_start = jts ! j_end = jte ! endif dt = dt_in qc_max = 0. qr_max = 0. qs_max = 0. qi_max = 0. qg_max = 0 ni_max = 0. nr_max = 0. imax_qc = 0 imax_qr = 0 imax_qi = 0 imax_qs = 0 imax_qg = 0 imax_ni = 0 imax_nr = 0 jmax_qc = 0 jmax_qr = 0 jmax_qi = 0 jmax_qs = 0 jmax_qg = 0 jmax_ni = 0 jmax_nr = 0 kmax_qc = 0 kmax_qr = 0 kmax_qi = 0 kmax_qs = 0 kmax_qg = 0 kmax_ni = 0 kmax_nr = 0 do i = 1, 256 mp_debug(i:i) = char(0) enddo if (.NOT. is_aerosol_aware .AND. PRESENT(nc) .AND. PRESENT(nwfa) & .AND. PRESENT(nifa) .AND. PRESENT(nwfa2d)) then write(mp_debug,*) 'WARNING, nc-nwfa-nifa-nwfa2d present but is_aerosol_aware is FALSE' CALL wrf_debug(0, mp_debug) endif j_loop: do j = j_start, j_end i_loop: do i = i_start, i_end pptrain = 0. pptcloud = 0. !RC pptsnow = 0. pptgraul = 0. pptice = 0. RAINNCV(i,j) = 0. IF ( PRESENT (snowncv) ) THEN SNOWNCV(i,j) = 0. ENDIF IF ( PRESENT (graupelncv) ) THEN GRAUPELNCV(i,j) = 0. ENDIF SR(i,j) = 0. do k = kts, kte t1d(k) = th(i,k,j)*pii(i,k,j) p1d(k) = p(i,k,j) w1d(k) = w(i,k,j) dz1d(k) = dz(i,k,j) qv1d(k) = qv(i,k,j) qc1d(k) = qc(i,k,j) qi1d(k) = qi(i,k,j) qr1d(k) = qr(i,k,j) qs1d(k) = qs(i,k,j) qg1d(k) = qg(i,k,j) ni1d(k) = ni(i,k,j) nr1d(k) = nr(i,k,j) rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) enddo if (is_aerosol_aware) then do k = kts, kte nc1d(k) = nc(i,k,j) nwfa1d(k) = nwfa(i,k,j) nifa1d(k) = nifa(i,k,j) if (wif_input_opt .eq. 2) then nbca1d(k) = nbca(i,k,j) else nbca1d(k) = 0.0 endif enddo nwfa1 = nwfa2d(i,j) else do k = kts, kte nc1d(k) = Nt_c/rho(k) nwfa1d(k) = 11.1E6/rho(k) nifa1d(k) = naIN1*0.01/rho(k) nbca1d(k) = 5.55E6/rho(k) enddo endif !..If not the variable-density graupel-hail hybrid, then set the vol mixing !.. ratio to mass mixing ratio divided by constant density (500kg/m3) value. if (is_hail_aware) then do k = kts, kte ng1d(k) = ng(i,k,j) qb1d(k) = qb(i,k,j) enddo else do k = kte, kts, -1 if (qg1d(k).gt.R1) then ygra1 = alog10(max(1.E-9, qg1d(k)*rho(k))) zans1 = 3.0 + 2./7.*(ygra1+8.) zans1 = MAX(2., MIN(zans1, 6.)) N0_exp = 10.**(zans1) lam_exp = (N0_exp*am_g(idx_bg1)*cgg(1,1)/(rho(k)*qg1d(k)))**oge1 lamg = lam_exp * (cgg(3,1)*ogg2*ogg1)**obmg ng1d(k) = cgg(2,1)*ogg3*rho(k)*qg1d(k)*lamg**bm_g / am_g(idx_bg1) ng1d(k) = MAX(R2, ng1d(k)/rho(k)) qb1d(k) = qg1d(k)/rho_g(idx_bg1) else ng1d(k) = 0 qb1d(k) = 0 endif enddo endif call mp_rcon(aer_init_opt, wif_input_opt, & qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, ni1d, & nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d, t1d, p1d, w1d, dz1d, & pptrain, pptcloud, pptsnow, pptgraul, pptice, & #if ( WRF_CHEM == 1 ) wetscav_on, rainprod1d, evapprod1d, & #endif kts, kte, dt, i, j) !RC pcp_ra(i,j) = pptrain pcp_sn(i,j) = pptsnow pcp_gr(i,j) = pptgraul pcp_ic(i,j) = pptice RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice + pptcloud !RC RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice + pptcloud !RC CLOUDNC(i,j) = CLOUDNC(i,j) + pptcloud !RC IF ( PRESENT(snowncv) .AND. PRESENT(snownc) ) THEN SNOWNCV(i,j) = pptsnow + pptice SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + pptice ENDIF IF ( PRESENT(graupelncv) .AND. PRESENT(graupelnc) ) THEN GRAUPELNCV(i,j) = pptgraul GRAUPELNC(i,j) = GRAUPELNC(i,j) + pptgraul ENDIF SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12) !.. !..BEGIN AEROSOL EMISSIONS !.. !..Reset lowest model level to initial state aerosols (fake sfc source). !.. Changed 13 May 2013 to fake emissions in which nwfa2d is aerosol !.. number tendency (number per kg per second). if (is_aerosol_aware) then !..Add anthropogenic emissions !-GT nwfa1d(kts) = nwfa1 nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt_in nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt_in if (wif_input_opt .eq. 2) then nbca1d(kts) = nbca1d(kts) + nbca2d(i,j)*dt_in else nbca1d(kts) = 0.0 endif !.. !..END AEROSOL EMISSIONS !.. do k = kts, kte nc(i,k,j) = nc1d(k) nwfa(i,k,j) = nwfa1d(k) nifa(i,k,j) = nifa1d(k) if (wif_input_opt .eq. 2) then nbca(i,k,j) = nbca1d(k) else nbca(i,k,j) = 0.0 endif enddo endif if (is_hail_aware) then do k = kts, kte ng(i,k,j) = ng1d(k) qb(i,k,j) = qb1d(k) enddo endif do k = kts, kte qv(i,k,j) = qv1d(k) qc(i,k,j) = qc1d(k) qi(i,k,j) = qi1d(k) qr(i,k,j) = qr1d(k) qs(i,k,j) = qs1d(k) qg(i,k,j) = qg1d(k) ni(i,k,j) = ni1d(k) nr(i,k,j) = nr1d(k) th(i,k,j) = t1d(k)/pii(i,k,j) #if ( WRF_CHEM == 1 ) IF ( wetscav_on ) THEN rainprod(i,k,j) = rainprod1d(k) evapprod(i,k,j) = evapprod1d(k) ENDIF #endif if (qc1d(k) .gt. qc_max) then imax_qc = i jmax_qc = j kmax_qc = k qc_max = qc1d(k) elseif (qc1d(k) .lt. 0.0) then write(mp_debug,*) 'WARNING, negative qc ', qc1d(k), & ' at i,j,k=', i,j,k CALL wrf_debug(150, mp_debug) endif if (qr1d(k) .gt. qr_max) then imax_qr = i jmax_qr = j kmax_qr = k qr_max = qr1d(k) elseif (qr1d(k) .lt. 0.0) then write(mp_debug,*) 'WARNING, negative qr ', qr1d(k), & ' at i,j,k=', i,j,k CALL wrf_debug(150, mp_debug) endif if (nr1d(k) .gt. nr_max) then imax_nr = i jmax_nr = j kmax_nr = k nr_max = nr1d(k) elseif (nr1d(k) .lt. 0.0) then write(mp_debug,*) 'WARNING, negative nr ', nr1d(k), & ' at i,j,k=', i,j,k CALL wrf_debug(150, mp_debug) endif if (qs1d(k) .gt. qs_max) then imax_qs = i jmax_qs = j kmax_qs = k qs_max = qs1d(k) elseif (qs1d(k) .lt. 0.0) then write(mp_debug,*) 'WARNING, negative qs ', qs1d(k), & ' at i,j,k=', i,j,k CALL wrf_debug(150, mp_debug) endif if (qi1d(k) .gt. qi_max) then imax_qi = i jmax_qi = j kmax_qi = k qi_max = qi1d(k) elseif (qi1d(k) .lt. 0.0) then write(mp_debug,*) 'WARNING, negative qi ', qi1d(k), & ' at i,j,k=', i,j,k CALL wrf_debug(150, mp_debug) endif if (qg1d(k) .gt. qg_max) then imax_qg = i jmax_qg = j kmax_qg = k qg_max = qg1d(k) elseif (qg1d(k) .lt. 0.0) then write(mp_debug,*) 'WARNING, negative qg ', qg1d(k), & ' at i,j,k=', i,j,k CALL wrf_debug(150, mp_debug) endif if (ni1d(k) .gt. ni_max) then imax_ni = i jmax_ni = j kmax_ni = k ni_max = ni1d(k) elseif (ni1d(k) .lt. 0.0) then write(mp_debug,*) 'WARNING, negative ni ', ni1d(k), & ' at i,j,k=', i,j,k CALL wrf_debug(150, mp_debug) endif if (qv1d(k) .lt. 0.0) then write(mp_debug,*) 'WARNING, negative qv ', qv1d(k), & ' at i,j,k=', i,j,k CALL wrf_debug(150, mp_debug) if (k.lt.kte-2 .and. k.gt.kts+1) then write(mp_debug,*) ' below and above are: ', qv(i,k-1,j), qv(i,k+1,j) CALL wrf_debug(150, mp_debug) qv(i,k,j) = MAX(1.E-7, 0.5*(qv(i,k-1,j) + qv(i,k+1,j))) else qv(i,k,j) = 1.E-7 endif endif enddo IF ( PRESENT (diagflag) ) THEN if (diagflag .and. do_radar_ref == 1) then IF ( present(ke_diag) ) THEN kediagloc = ke_diag ELSE kediagloc = kte ENDIF call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & ng1d, qb1d, t1d, p1d, dBZ, kts, kte, i, j, kediagloc) do k = kts, kte refl_10cm(i,k,j) = MAX(-35., dBZ(k)) enddo endif ENDIF IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN do k = kts, kte re_qc1d(k) = RE_QC_BG re_qi1d(k) = RE_QI_BG re_qs1d(k) = RE_QS_BG enddo call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & re_qc1d, re_qi1d, re_qs1d, kts, kte) do k = kts, kte re_cloud(i,k,j) = MAX(RE_QC_BG, MIN(re_qc1d(k), 75.E-6)) !RC re_ice(i,k,j) = MAX(RE_QI_BG, MIN(re_qi1d(k), 125.E-6)) re_snow(i,k,j) = MAX(RE_QS_BG, MIN(re_qs1d(k), 999.E-6)) enddo ENDIF enddo i_loop enddo j_loop ! DEBUG - GT write(mp_debug,'(a,7(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', & 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', & 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', & 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', & 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', & 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', & 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', & 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')' CALL wrf_debug(150, mp_debug) ! END DEBUG - GT do i = 1, 256 mp_debug(i:i) = char(0) enddo END SUBROUTINE mp_rcon_driver !+---+-----------------------------------------------------------------+ !ctrlL !+---+-----------------------------------------------------------------+ !+---+-----------------------------------------------------------------+ !.. This subroutine computes the moisture tendencies of water vapor, !.. cloud droplets, rain, cloud ice (pristine), snow, and graupel. !.. Previously this code was based on Reisner et al (1998), but few of !.. those pieces remain. A complete description is now found in !.. Thompson et al. (2004, 2008). !+---+-----------------------------------------------------------------+ ! subroutine mp_rcon (aer_init_opt, wif_input_opt, & qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, & ni1d, nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d, & t1d, p1d, w1d, dzq, & pptrain, pptcloud, pptsnow, pptgraul, pptice, & !RC #if ( WRF_CHEM == 1 ) wetscav_on, rainprod, evapprod, & #endif kts, kte, dt, ii, jj) implicit none !..Sub arguments INTEGER, OPTIONAL, INTENT(IN):: aer_init_opt, wif_input_opt INTEGER, INTENT(IN):: kts, kte, ii, jj REAL, DIMENSION(kts:kte), INTENT(INOUT):: & qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, qb1d, & ni1d, nr1d, nc1d, ng1d, nwfa1d, nifa1d, nbca1d, t1d REAL, DIMENSION(kts:kte), INTENT(IN):: p1d, w1d, dzq REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice, pptcloud !RC REAL, INTENT(IN):: dt #if ( WRF_CHEM == 1 ) REAL, DIMENSION(kts:kte), INTENT(INOUT):: & rainprod, evapprod LOGICAL, INTENT(IN) :: wetscav_on #endif !..Local variables REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, qrten, & qsten, qgten, qbten, niten, nrten, ncten, ngten, nwfaten, nifaten, & nbcaten DOUBLE PRECISION, DIMENSION(kts:kte):: prw_vcd DOUBLE PRECISION, DIMENSION(kts:kte):: pnc_wcd, & pnc_scw, pnc_gcw DOUBLE PRECISION, DIMENSION(kts:kte):: pna_rca, pna_sca, pna_gca, & pnd_rcd, pnd_scd, pnd_gcd, pnb_rcb, pnb_scb, pnb_gcb DOUBLE PRECISION, DIMENSION(kts:kte):: prr_rcs, & prr_rcg, prr_sml, prr_gml, & prr_rci, prv_rev, & pnr_rcs, pnr_rcg, & pnr_rci, pnr_sml, pnr_gml, & pnr_rev, pnr_rcr, pnr_rfz DOUBLE PRECISION, DIMENSION(kts:kte):: pri_inu, pni_inu, pri_ihm, & pni_ihm, pri_wfz, pni_wfz, & pri_rfz, pni_rfz, pri_ide, & pni_ide, pri_rci, pni_rci, & pni_sci, pni_iau, pri_iha, pni_iha DOUBLE PRECISION, DIMENSION(kts:kte):: prs_iau, prs_sci, prs_rcs, & prs_scw, prs_sde, prs_ihm, & prs_ide DOUBLE PRECISION, DIMENSION(kts:kte):: prg_scw, prg_rfz, prg_gde, & prg_gcw, prg_rci, prg_rcs, prg_rcg, prg_ihm, & png_rcs, png_rcg, png_scw, png_gde, & pbg_scw, pbg_rfz, pbg_gcw, pbg_rci, pbg_rcs, pbg_rcg, & pbg_sml, pbg_gml DOUBLE PRECISION, PARAMETER:: zeroD0 = 0.0d0 REAL, DIMENSION(kts:kte):: temp, twet, pres, qv REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, rb REAL, DIMENSION(kts:kte):: ni, nr, nc, ns, ng, nwfa, nifa, nbca REAL, DIMENSION(kts:kte):: rho, rhof, rhof2 REAL, DIMENSION(kts:kte):: qvs, qvsi, delQvs REAL, DIMENSION(kts:kte):: satw, sati, ssatw, ssati REAL, DIMENSION(kts:kte):: diffu, visco, vsc2, & tcond, lvap, ocp, lvt2 DOUBLE PRECISION, DIMENSION(kts:kte):: ilamg, N0_r, N0_g DOUBLE PRECISION:: N0_melt REAL, DIMENSION(kts:kte):: mvd_r, mvd_c, mvd_g REAL, DIMENSION(kts:kte):: smob, smo2, smo1, smo0, & smoc, smod, smoe, smof, smog REAL, DIMENSION(kts:kte):: sed_r,sed_s,sed_g,sed_i,sed_n,sed_c,sed_b REAL:: rgvm, delta_tp, orho, lfus2 REAL, DIMENSION(5):: onstep DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg DOUBLE PRECISION:: lami, ilami REAL:: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m DOUBLE PRECISION:: Dr_star, Dc_star REAL:: zeta1, zeta, taud, tau REAL:: stoke_r, stoke_s, stoke_g, stoke_i REAL:: vti, vtr, vts, vtg, vtc REAL:: xrho_g, afall, vtg1, vtg2 REAL:: bfall = 3*b_coeff - 1 REAL, DIMENSION(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk, & vtngk, vtck, vtnck REAL, DIMENSION(kts:kte):: vts_boost REAL:: M0, slam1, slam2 REAL:: Mrat, ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts, C_snow REAL:: a_, b_, loga_, A1, A2, tf REAL:: tempc, tc0, r_mvd1, r_mvd2, xkrat REAL:: dew_t, Tlcl, The REAL:: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr, xrg, xng, xrb REAL:: xsat, rate_max, sump, ratio REAL:: clap, fcd, dfcd REAL:: otemp, rvs, rvs_p, rvs_pp, gamsc, alphsc, t1_evap, t1_subl REAL:: r_frac, g_frac, const_Ri, rime_dens REAL:: Ef_rw, Ef_sw, Ef_gw, Ef_rr REAL:: Ef_ra, Ef_sa, Ef_ga REAL:: dtsave, odts, odt, odzq, hgt_agl REAL:: xslw1, ygra1, zans1, eva_factor REAL:: SR, melt_f INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq, k_melting INTEGER, DIMENSION(5):: ksed1 INTEGER:: nir, nis, nig, nii, nic, niin INTEGER:: idx_tc, idx_t, idx_s, idx_g1, idx_g, idx_r1, idx_r, & idx_i1, idx_i, idx_c, idx, idx_d, idx_n, idx_in INTEGER, DIMENSION(kts:kte):: idx_bg LOGICAL:: melti, no_micro LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg LOGICAL:: debug_flag CHARACTER*256:: mp_debug INTEGER:: nu_c !RC; processes and microphysical quantities necessary for this scheme. !RC; included below in case someone wants to output easier. DOUBLE PRECISION, DIMENSION(kts:kte):: & prr_wau,pnr_wau,pnc_wau, & prr_rcw,pnc_rcw,ilamr,ilamc !/RC !RC; lognormal variables REAL:: ln_varx,ln_tmp1,ln_tmp2 REAL, DIMENSION(kts:kte):: ln_dn REAL, DIMENSION(kts:kte):: ln_sigma REAL:: sigma_d !/RC !+---+ debug_flag = .false. ! if (ii.eq.901 .and. jj.eq.379) debug_flag = .true. if(debug_flag) then write(mp_debug, *) 'DEBUG INFO, mp_rcon at (i,j) ', ii, ', ', jj CALL wrf_debug(550, mp_debug) endif no_micro = .true. dtsave = dt odt = 1./dt odts = 1./dtsave iexfrq = 1 !+---+-----------------------------------------------------------------+ !.. Source/sink terms. First 2 chars: "pr" represents source/sink of !.. mass while "pn" represents source/sink of number. Next char is one !.. of "v" for water vapor, "r" for rain, "i" for cloud ice, "w" for !.. cloud water, "s" for snow, and "g" for graupel. Next chars !.. represent processes: "de" for sublimation/deposition, "ev" for !.. evaporation, "fz" for freezing, "ml" for melting, "au" for !.. autoconversion, "nu" for ice nucleation, "hm" for Hallet/Mossop !.. secondary ice production, and "c" for collection followed by the !.. character for the species being collected. ALL of these terms are !.. positive (except for deposition/sublimation terms which can switch !.. signs based on super/subsaturation) and are treated as negatives !.. where necessary in the tendency equations. !+---+-----------------------------------------------------------------+ do k = kts, kte tten(k) = 0. qvten(k) = 0. qcten(k) = 0. qiten(k) = 0. qrten(k) = 0. qsten(k) = 0. qgten(k) = 0. ngten(k) = 0. qbten(k) = 0. niten(k) = 0. nrten(k) = 0. ncten(k) = 0. nwfaten(k) = 0. nifaten(k) = 0. nbcaten(k) = 0. prw_vcd(k) = 0. pnc_wcd(k) = 0. pnc_wau(k) = 0. pnc_rcw(k) = 0. pnc_scw(k) = 0. pnc_gcw(k) = 0. prv_rev(k) = 0. prr_wau(k) = 0. prr_rcw(k) = 0. prr_rcs(k) = 0. prr_rcg(k) = 0. prr_sml(k) = 0. prr_gml(k) = 0. prr_rci(k) = 0. pnr_wau(k) = 0. pnr_rcs(k) = 0. pnr_rcg(k) = 0. pnr_rci(k) = 0. pnr_sml(k) = 0. pnr_gml(k) = 0. pnr_rev(k) = 0. pnr_rcr(k) = 0. pnr_rfz(k) = 0. pri_inu(k) = 0. pni_inu(k) = 0. pri_ihm(k) = 0. pni_ihm(k) = 0. pri_wfz(k) = 0. pni_wfz(k) = 0. pri_rfz(k) = 0. pni_rfz(k) = 0. pri_ide(k) = 0. pni_ide(k) = 0. pri_rci(k) = 0. pni_rci(k) = 0. pni_sci(k) = 0. pni_iau(k) = 0. pri_iha(k) = 0. pni_iha(k) = 0. prs_iau(k) = 0. prs_sci(k) = 0. prs_rcs(k) = 0. prs_scw(k) = 0. prs_sde(k) = 0. prs_ihm(k) = 0. prs_ide(k) = 0. prg_scw(k) = 0. prg_rfz(k) = 0. prg_gde(k) = 0. prg_gcw(k) = 0. prg_rci(k) = 0. prg_rcs(k) = 0. prg_rcg(k) = 0. prg_ihm(k) = 0. ! new source/sink terms for 3-moment graupel png_scw(k) = 0. png_rcs(k) = 0. png_rcg(k) = 0. png_gde(k) = 0. pbg_scw(k) = 0. pbg_rfz(k) = 0. pbg_gcw(k) = 0. pbg_rci(k) = 0. pbg_rcs(k) = 0. pbg_rcg(k) = 0. pbg_sml(k) = 0. pbg_gml(k) = 0. pna_rca(k) = 0. pna_sca(k) = 0. pna_gca(k) = 0. pnd_rcd(k) = 0. pnd_scd(k) = 0. pnd_gcd(k) = 0. pnb_rcb(k) = 0. pnb_scb(k) = 0. pnb_gcb(k) = 0. enddo #if ( WRF_CHEM == 1 ) if ( wetscav_on ) then do k = kts, kte rainprod(k) = 0. evapprod(k) = 0. enddo endif #endif !..Bug fix (2016Jun15), prevent use of uninitialized value(s) of snow moments. do k = kts, kte smo0(k) = 0. smo1(k) = 0. smo2(k) = 0. smob(k) = 0. smoc(k) = 0. smod(k) = 0. smoe(k) = 0. smof(k) = 0. smog(k) = 0. ns(k) = 0. mvd_r(k) = 0. mvd_c(k) = 0. enddo !+---+-----------------------------------------------------------------+ !..Put column of data into local arrays. !+---+-----------------------------------------------------------------+ do k = kts, kte temp(k) = t1d(k) qv(k) = MAX(1.E-10, qv1d(k)) pres(k) = p1d(k) rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622)) if (is_aerosol_aware) then if (aer_init_opt .lt. 2) then ! Constant or climatology (e.g., GOCART) nwfa(k) = MAX(11.1E6, MIN(9999.E6, nwfa1d(k)*rho(k))) nifa(k) = MAX(naIN1*0.01, MIN(9999.E6, nifa1d(k)*rho(k))) if (wif_input_opt .eq. 2) then ! Considering BC aerosol nbca(k) = MAX(5.55E6, MIN(9999.E6, nbca1d(k)*rho(k))) else nbca(k) = 0.0 endif else ! First guess (e.g., GEOS-5) nwfa(k) = MAX(0.0, nwfa1d(k)*rho(k)) nifa(k) = MAX(0.0, nifa1d(k)*rho(k)) if (wif_input_opt .eq. 2) then ! Considering BC aerosol nbca(k) = MAX(0.0, nbca1d(k)*rho(k)) else nbca(k) = 0.0 endif endif else nwfa(k) = MAX(11.1E6, MIN(9999.E6, nwfa1d(k)*rho(k))) nifa(k) = MAX(naIN1*0.01, MIN(9999.E6, nifa1d(k)*rho(k))) nbca(k) = MAX(5.55E6, MIN(9999.E6, nbca1d(k)*rho(k))) endif ln_sigma(k) = 0.20 !RC !RC; This whole if block is written to accomodate lognormal cloud water. if (qc1d(k) .gt. R1) then no_micro = .false. rc(k) = qc1d(k)*rho(k) nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max)) L_qc(k) = .true. ! RC; set variable sigma sigma_d = (rc(k)/nc(k))**(0.3333) ln_sigma(k) = sigma_d*(-1.185E3) + 0.815 ln_sigma(k) = MAX(ln_sigma(k) , 0.2) ln_sigma(k) = MIN(ln_sigma(k) , 0.7) !/RC xDc = ( ( (1./am_r)*(rc(k)/nc(k)) )**obmr ) if (xDc.lt. D0c) then xDc = D0c elseif (xDc.gt. D0r*2.) then xDc = D0r*2. endif ln_dn(k) = xDc * EXP(-1.5 * ln_sigma(k)*ln_sigma(k)) !RC nc(k) = MAX( 2. , MIN( DBLE(Nt_c_max) , (1./am_r)*rc(k)*EXP(-4.5*ln_sigma(k)*ln_sigma(k))*(1. / ( xDc*EXP(-1.5*ln_sigma(k)*ln_sigma(k)) ))**3. ) ) !RC mvd_c(k) = MIN( D0r*2., ln_dn(k)*EXP(3.*ln_sigma(k)*ln_sigma(k)) ) !RC if (.NOT. is_aerosol_aware) nc(k) = Nt_c else qc1d(k) = 0.0 nc1d(k) = 0.0 rc(k) = R1 nc(k) = 2. L_qc(k) = .false. endif !\RC if (qi1d(k) .gt. R1) then no_micro = .false. ri(k) = qi1d(k)*rho(k) ni(k) = MAX(R2, ni1d(k)*rho(k)) if (ni(k).le. R2) then lami = cie(2)/5.E-6 ni(k) = MIN(999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) endif L_qi(k) = .true. lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi ilami = 1./lami xDi = (bm_i + mu_i + 1.) * ilami if (xDi.lt. 5.E-6) then lami = cie(2)/5.E-6 ni(k) = MIN(999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) elseif (xDi.gt. 300.E-6) then lami = cie(2)/300.E-6 ni(k) = cig(1)*oig2*ri(k)/am_i*lami**bm_i endif else qi1d(k) = 0.0 ni1d(k) = 0.0 ri(k) = R1 ni(k) = R2 L_qi(k) = .false. endif if (qr1d(k) .gt. R1) then no_micro = .false. rr(k) = qr1d(k)*rho(k) nr(k) = MAX(R2, nr1d(k)*rho(k)) if (nr(k).le. R2) then mvd_r(k) = 1.0E-3 lamr = (3.0 + mu_r + 0.672) / mvd_r(k) nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r endif L_qr(k) = .true. lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr mvd_r(k) = (3.0 + mu_r + 0.672) / lamr if (mvd_r(k) .gt. 2.5E-3) then mvd_r(k) = 2.5E-3 lamr = (3.0 + mu_r + 0.672) / mvd_r(k) nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r elseif (mvd_r(k) .lt. D0r*0.75) then mvd_r(k) = D0r*0.75 lamr = (3.0 + mu_r + 0.672) / mvd_r(k) nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r endif else qr1d(k) = 0.0 nr1d(k) = 0.0 rr(k) = R1 nr(k) = R2 L_qr(k) = .false. endif if (qs1d(k) .gt. R1) then no_micro = .false. rs(k) = qs1d(k)*rho(k) L_qs(k) = .true. else qs1d(k) = 0.0 rs(k) = R1 L_qs(k) = .false. endif if (qg1d(k) .gt. R1) then no_micro = .false. L_qg(k) = .true. rg(k) = qg1d(k)*rho(k) ng(k) = MAX(R2, ng1d(k)*rho(k)) rb(k) = MAX(qg1d(k)/rho_g(NRHG), qb1d(k)) rb(k) = MIN(qg1d(k)/rho_g(1), rb(k)) qb1d(k) = rb(k) idx_bg(k) = MAX(1,MIN(NINT(qg1d(k)/rb(k) *0.01)+1,NRHG)) if (ng(k).le. R2) then mvd_g(k) = 1.5E-3 lamg = (3.0 + mu_g + 0.672) / mvd_g(k) ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k)) endif lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg mvd_g(k) = (3.0 + mu_g + 0.672) / lamg if (mvd_g(k) .gt. 25.4E-3) then mvd_g(k) = 25.4E-3 lamg = (3.0 + mu_g + 0.672) / mvd_g(k) ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k)) elseif (mvd_g(k) .lt. D0r) then mvd_g(k) = D0r lamg = (3.0 + mu_g + 0.672) / mvd_g(k) ng(k) = cgg(2,1)*ogg3*rg(k)*lamg**bm_g / am_g(idx_bg(k)) endif else qg1d(k) = 0.0 ng1d(k) = 0.0 qb1d(k) = 0.0 idx_bg(k) = idx_bg1 rg(k) = R1 ng(k) = R2 rb(k) = R1/rho(k)/rho_g(NRHG) L_qg(k) = .false. endif if (.not. is_hail_aware) idx_bg(k) = idx_bg1 enddo !+---+-----------------------------------------------------------------+ ! if (debug_flag) then ! write(mp_debug,*) 'DEBUG-VERBOSE at (i,j) ', ii, ', ', jj ! CALL wrf_debug(550, mp_debug) ! do k = kts, kte ! write(mp_debug, '(a,i3,f8.2,1x,f7.2,1x, 13(1x,e13.6))') & ! & 'VERBOSE: ', k, pres(k)*0.01, temp(k)-273.15, qv(k), rc(k), rr(k), ri(k), rs(k), rg(k), nc(k), nr(k), ni(k), ng(k), rb(k), nwfa(k), nifa(k) ! CALL wrf_debug(550, mp_debug) ! enddo ! endif !+---+-----------------------------------------------------------------+ !+---+-----------------------------------------------------------------+ !..Derive various thermodynamic variables frequently used. !.. Saturation vapor pressure (mixing ratio) over liquid/ice comes from !.. Flatau et al. 1992; enthalpy (latent heat) of vaporization from !.. Bohren & Albrecht 1998; others from Pruppacher & Klett 1978. !+---+-----------------------------------------------------------------+ k_melting = 0 do k = kts, kte tempc = temp(k) - 273.15 rhof(k) = SQRT(RHO_NOT/rho(k)) rhof2(k) = SQRT(rhof(k)) qvs(k) = rslf(pres(k), temp(k)) delQvs(k) = MAX(0.0, rslf(pres(k), 273.15)-qv(k)) if (tempc .le. 0.0) then qvsi(k) = rsif(pres(k), temp(k)) else qvsi(k) = qvs(k) k_melting = MAX(k, k_melting) endif satw(k) = qv(k)/qvs(k) sati(k) = qv(k)/qvsi(k) ssatw(k) = satw(k) - 1. ssati(k) = sati(k) - 1. if (abs(ssatw(k)).lt. eps) ssatw(k) = 0.0 if (abs(ssati(k)).lt. eps) ssati(k) = 0.0 if (no_micro .and. ssati(k).gt. 0.0) no_micro = .false. diffu(k) = 2.11E-5*(temp(k)/273.15)**1.94 * (101325./pres(k)) if (tempc .ge. 0.0) then visco(k) = (1.718+0.0049*tempc)*1.0E-5 else visco(k) = (1.718+0.0049*tempc-1.2E-5*tempc*tempc)*1.0E-5 endif ocp(k) = 1./(Cp*(1.+0.887*qv(k))) vsc2(k) = SQRT(rho(k)/visco(k)) lvap(k) = lvap0 + (2106.0 - 4218.0)*tempc tcond(k) = (5.69 + 0.0168*tempc)*1.0E-5 * 418.936 twet(k) = temp(k) enddo if (k_melting .gt. 0) then do k = kts, k_melting if (satw(k) .lt. 0.999) then dew_t = MIN(temp(k)-0.001, t_dew(pres(k), qv(k))) Tlcl = t_lcl(temp(k), dew_t) The = theta_e(pres(k), temp(k), qv(k), Tlcl) twet(k) = MIN(temp(k), compT_fr_The(The, pres(k))) endif enddo endif !+---+-----------------------------------------------------------------+ !..If no existing hydrometeor species and no chance to initiate ice or !.. condense cloud water, just exit quickly! !+---+-----------------------------------------------------------------+ if (no_micro) return !+---+-----------------------------------------------------------------+ !..Calculate y-intercept, slope, and useful moments for snow. !+---+-----------------------------------------------------------------+ if (.not. iiwarm) then do k = kts, kte if (.not. L_qs(k)) CYCLE tc0 = MIN(-0.1, temp(k)-273.15) smob(k) = rs(k)*oams !..All other moments based on reference, 2nd moment. If bm_s.ne.2, !.. then we must compute actual 2nd moment and use as reference. if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then smo2(k) = smob(k) else loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s & + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 & + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s & + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 & + sa(10)*bm_s*bm_s*bm_s a_ = 10.0**loga_ b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s & + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 & + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s & + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 & + sb(10)*bm_s*bm_s*bm_s smo2(k) = (smob(k)/a_)**(1./b_) endif !..Calculate 0th moment. Represents snow number concentration. loga_ = sa(1) + sa(2)*tc0 + sa(5)*tc0*tc0 + sa(9)*tc0*tc0*tc0 a_ = 10.0**loga_ b_ = sb(1) + sb(2)*tc0 + sb(5)*tc0*tc0 + sb(9)*tc0*tc0*tc0 smo0(k) = a_ * smo2(k)**b_ !..Calculate 1st moment. Useful for depositional growth and melting. loga_ = sa(1) + sa(2)*tc0 + sa(3) & + sa(4)*tc0 + sa(5)*tc0*tc0 & + sa(6) + sa(7)*tc0*tc0 & + sa(8)*tc0 + sa(9)*tc0*tc0*tc0 & + sa(10) a_ = 10.0**loga_ b_ = sb(1)+ sb(2)*tc0 + sb(3) + sb(4)*tc0 & + sb(5)*tc0*tc0 + sb(6) & + sb(7)*tc0*tc0 + sb(8)*tc0 & + sb(9)*tc0*tc0*tc0 + sb(10) smo1(k) = a_ * smo2(k)**b_ !..Calculate bm_s+1 (th) moment. Useful for diameter calcs. loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & + sa(10)*cse(1)*cse(1)*cse(1) a_ = 10.0**loga_ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) smoc(k) = a_ * smo2(k)**b_ !..Calculate snow number concentration (explicit integral, not smo0) M0 = smob(k)/smoc(k) Mrat = smob(k)*M0*M0*M0 slam1 = M0 * Lam0 slam2 = M0 * Lam1 ns(k) = Mrat*Kap0/slam1 & + Mrat*Kap1*M0**mu_s*csg(15)/slam2**cse(15) !..Calculate bv_s+2 (th) moment. Useful for riming. loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(13) & + sa(4)*tc0*cse(13) + sa(5)*tc0*tc0 & + sa(6)*cse(13)*cse(13) + sa(7)*tc0*tc0*cse(13) & + sa(8)*tc0*cse(13)*cse(13) + sa(9)*tc0*tc0*tc0 & + sa(10)*cse(13)*cse(13)*cse(13) a_ = 10.0**loga_ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(13) + sb(4)*tc0*cse(13) & + sb(5)*tc0*tc0 + sb(6)*cse(13)*cse(13) & + sb(7)*tc0*tc0*cse(13) + sb(8)*tc0*cse(13)*cse(13) & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(13)*cse(13)*cse(13) smoe(k) = a_ * smo2(k)**b_ !..Calculate 1+(bv_s+1)/2 (th) moment. Useful for depositional growth. loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(16) & + sa(4)*tc0*cse(16) + sa(5)*tc0*tc0 & + sa(6)*cse(16)*cse(16) + sa(7)*tc0*tc0*cse(16) & + sa(8)*tc0*cse(16)*cse(16) + sa(9)*tc0*tc0*tc0 & + sa(10)*cse(16)*cse(16)*cse(16) a_ = 10.0**loga_ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(16) + sb(4)*tc0*cse(16) & + sb(5)*tc0*tc0 + sb(6)*cse(16)*cse(16) & + sb(7)*tc0*tc0*cse(16) + sb(8)*tc0*cse(16)*cse(16) & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(16)*cse(16)*cse(16) smof(k) = a_ * smo2(k)**b_ !..Calculate bm_s + bv_s+2 (th) moment. Useful for riming into graupel. loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(17) & + sa(4)*tc0*cse(17) + sa(5)*tc0*tc0 & + sa(6)*cse(17)*cse(17) + sa(7)*tc0*tc0*cse(17) & + sa(8)*tc0*cse(17)*cse(17) + sa(9)*tc0*tc0*tc0 & + sa(10)*cse(17)*cse(17)*cse(17) a_ = 10.0**loga_ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(17) + sb(4)*tc0*cse(17) & + sb(5)*tc0*tc0 + sb(6)*cse(17)*cse(17) & + sb(7)*tc0*tc0*cse(17) + sb(8)*tc0*cse(17)*cse(17) & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(17)*cse(17)*cse(17) smog(k) = a_ * smo2(k)**b_ enddo !+---+-----------------------------------------------------------------+ !..Calculate y-intercept, slope values for graupel. !+---+-----------------------------------------------------------------+ do k = kte, kts, -1 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg ilamg(k) = 1./lamg N0_g(k) = ng(k)*ogg2*lamg**cge(2,1) enddo endif !+---+-----------------------------------------------------------------+ !..Calculate y-intercept, slope values for rain. !+---+-----------------------------------------------------------------+ do k = kte, kts, -1 lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr ilamr(k) = 1./lamr mvd_r(k) = (3.0 + mu_r + 0.672) / lamr N0_r(k) = nr(k)*org2*lamr**cre(2) enddo !+---+-----------------------------------------------------------------+ !..Compute warm-rain process terms (except evap done later). !+---+-----------------------------------------------------------------+ do k = kts, kte !..Rain self-collection follows Seifert, 1994 and drop break-up !.. follows Verlinde and Cotton, 1993. if (L_qr(k) .and. mvd_r(k).gt. D0r) then Ef_rr = 1.0 - EXP(2300.0*(mvd_r(k)-1950.0E-6)) pnr_rcr(k) = Ef_rr * 2.0*nr(k)*rr(k) endif !RC mvd_c(k) = D0c ln_dn(k) = 0.0 if (L_qc(k)) then xDc = ( ( (1./am_r)*(rc(k)/nc(k)) )**obmr ) ! dv_bar in paper ln_dn(k) = xDc * EXP(-1.5 * ln_sigma(k)*ln_sigma(k)) mvd_c(k) = MIN( D0r*2., ln_dn(k)*EXP(3.*ln_sigma(k)*ln_sigma(k)) ) endif !/RC !RC; Autoconversion follows Nickerson et al. 1986 if (rc(k) .gt. 0.01e-3) then ln_varx = EXP(9*ln_sigma(k)*ln_sigma(k))-1. ln_tmp1 = 0.067*rc(k)*rc(k)*(1.0E16*((rc(k)/nc(k))**1.3333)*((ln_varx)**0.5) - 2.7 ) ln_tmp2 = 1.0E4*( (rc(k)*(ln_varx**0.5)/nc(k) )**obmr ) - 1.2 prr_wau(k) = MAX( 0.d0,DBLE(ln_tmp1)*DBLE(ln_tmp2)) prr_wau(k) = MIN(DBLE(rc(k)*odts), prr_wau(k)) pnr_wau(k) = prr_wau(k) / (am_r*200.*D0r*D0r*D0r) pnc_wau(k) = MIN(DBLE(nc(k)*odts), prr_wau(k) & / (am_r*mvd_c(k)*mvd_c(k)*mvd_c(k))) endif !RC; End new autoconversion !..Rain collecting cloud water. In CE, assume Dc<1). Either way, only bother to do sedimentation below !.. 1st level that contains any sedimenting particles (k=ksed1 on down). !.. New in v3.0+ is computing separate for rain, ice, snow, and !.. graupel species thus making code faster with credit to J. Schmidt. !+---+-----------------------------------------------------------------+ nstep = 0 onstep(:) = 1.0 ksed1(:) = 1 do k = kte+1, kts, -1 vtrk(k) = 0. vtnrk(k) = 0. vtik(k) = 0. vtnik(k) = 0. vtsk(k) = 0. vtgk(k) = 0. vtngk(k) = 0. vtck(k) = 0. vtnck(k) = 0. enddo if (ANY(L_qr .eqv. .true.)) then do k = kte, kts, -1 vtr = 0. rhof(k) = SQRT(RHO_NOT/rho(k)) if (rr(k).gt. R1) then lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) & *((lamr+fv_r)**(-cre(6))) vtrk(k) = vtr ! First below is technically correct: ! vtr = rhof(k)*av_r*crg(5)*org2 * lamr**cre(2) & ! *((lamr+fv_r)**(-cre(5))) ! Test: make number fall faster (but still slower than mass) ! Goal: less prominent size sorting vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) & *((lamr+fv_r)**(-cre(7))) vtnrk(k) = vtr else vtrk(k) = vtrk(k+1) vtnrk(k) = vtnrk(k+1) endif if (MAX(vtrk(k),vtnrk(k)) .gt. 1.E-3) then ksed1(1) = MAX(ksed1(1), k) delta_tp = dzq(k)/(MAX(vtrk(k),vtnrk(k))) nstep = MAX(nstep, INT(DT/delta_tp + 1.)) endif enddo if (ksed1(1) .eq. kte) ksed1(1) = kte-1 if (nstep .gt. 0) onstep(1) = 1./REAL(nstep) endif !+---+-----------------------------------------------------------------+ if (ANY(L_qc .eqv. .true.)) then hgt_agl = 0. do k = kts, kte-1 if (rc(k) .gt. R2) ksed1(5) = k hgt_agl = hgt_agl + dzq(k) if (hgt_agl .gt. 500.0) goto 151 enddo 151 continue do k = ksed1(5), kts, -1 vtc = 0. !RC; Adjustments to cloud fallspeed integrals. Allows drizzle mode. if (rc(k) .gt. R1 .and. w1d(k) .lt. 0.1E-1) then !RC; allow sedimentation only when w<0.01m/s xDc = ( ( (1./am_r)*(rc(k)/nc(k)) )**obmr ) vtc = (av_c*(xDc * EXP(-1.5 * ln_sigma(k)*ln_sigma(k)))**bv_c)*EXP(5.*ln_sigma(k)*ln_sigma(k)) vtck(k) = MIN(0.1d0,vtc) !RC; These get big, so we need to limit it to 0.1 vtc = (av_c*(xDc * EXP(-1.5 * ln_sigma(k)*ln_sigma(k)))**bv_c)*EXP(1.3*ln_sigma(k)*ln_sigma(k)) vtnck(k) = MIN(0.1d0,vtc) !\RC endif enddo endif !+---+-----------------------------------------------------------------+ if (.not. iiwarm) then if (ANY(L_qi .eqv. .true.)) then nstep = 0 do k = kte, kts, -1 vti = 0. if (ri(k).gt. R1) then lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi ilami = 1./lami vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i vtik(k) = vti ! First below is technically correct: ! vti = rhof(k)*av_i*cig(4)*oig1 * ilami**bv_i ! Goal: less prominent size sorting vti = rhof(k)*av_i*cig(6)/cig(7) * ilami**bv_i vtnik(k) = vti else vtik(k) = vtik(k+1) vtnik(k) = vtnik(k+1) endif if (vtik(k) .gt. 1.E-3) then ksed1(2) = MAX(ksed1(2), k) delta_tp = dzq(k)/vtik(k) nstep = MAX(nstep, INT(DT/delta_tp + 1.)) endif enddo if (ksed1(2) .eq. kte) ksed1(2) = kte-1 if (nstep .gt. 0) onstep(2) = 1./REAL(nstep) endif !+---+-----------------------------------------------------------------+ if (ANY(L_qs .eqv. .true.)) then nstep = 0 do k = kte, kts, -1 vts = 0. if (rs(k).gt. R1) then xDs = smoc(k) / smob(k) Mrat = 1./xDs ils1 = 1./(Mrat*Lam0 + fv_s) ils2 = 1./(Mrat*Lam1 + fv_s) t1_vts = Kap0*csg(4)*ils1**cse(4) t2_vts = Kap1*Mrat**mu_s*csg(10)*ils2**cse(10) ils1 = 1./(Mrat*Lam0) ils2 = 1./(Mrat*Lam1) t3_vts = Kap0*csg(1)*ils1**cse(1) t4_vts = Kap1*Mrat**mu_s*csg(7)*ils2**cse(7) vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts) if (prr_sml(k) .gt. 0.0) then SR = rs(k)/(rs(k)+rr(k)) vtsk(k) = vts*SR + (1.-SR)*vtrk(k) else vtsk(k) = vts*vts_boost(k) endif else vtsk(k) = vtsk(k+1) endif if (vtsk(k) .gt. 1.E-3) then ksed1(3) = MAX(ksed1(3), k) delta_tp = dzq(k)/vtsk(k) nstep = MAX(nstep, INT(DT/delta_tp + 1.)) endif enddo if (ksed1(3) .eq. kte) ksed1(3) = kte-1 if (nstep .gt. 0) onstep(3) = 1./REAL(nstep) endif !+---+-----------------------------------------------------------------+ if (ANY(L_qg .eqv. .true.)) then nstep = 0 do k = kte, kts, -1 vtg = 0. if (rg(k).gt. R1) then if (is_hail_aware) then xrho_g = MAX(rho_g(1),MIN(rg(k)/rho(k)/rb(k),rho_g(NRHG))) afall = a_coeff*((4.0*xrho_g*9.8)/(3.0*rho(k)))**b_coeff afall = afall * visco(k)**(1.0-2.0*b_coeff) else afall = av_g_old bfall = bv_g_old endif vtg = rhof(k)*afall*cgg(6,idx_bg(k))*ogg3 * ilamg(k)**bfall vtgk(k) = vtg ! Goal: less prominent size sorting ! the ELSE section below is technically (mathematically) correct: if (mu_g .eq. 0) then vtg = rhof(k)*afall*cgg(7,idx_bg(k))/cgg(12,idx_bg(k)) * ilamg(k)**bfall else vtg = rhof(k)*afall*cgg(8,idx_bg(k))*ogg2 * ilamg(k)**bfall endif vtngk(k) = vtg else vtgk(k) = vtgk(k+1) vtngk(k) = vtngk(k+1) endif if (vtgk(k) .gt. 1.E-3) then ksed1(4) = MAX(ksed1(4), k) delta_tp = dzq(k)/vtgk(k) nstep = MAX(nstep, INT(DT/delta_tp + 1.)) endif enddo if (ksed1(4) .eq. kte) ksed1(4) = kte-1 if (nstep .gt. 0) onstep(4) = 1./REAL(nstep) endif endif !+---+-----------------------------------------------------------------+ !..Sedimentation of mixing ratio is the integral of v(D)*m(D)*N(D)*dD, !.. whereas neglect m(D) term for number concentration. Therefore, !.. cloud ice has proper differential sedimentation. !+---+-----------------------------------------------------------------+ if (ANY(L_qr .eqv. .true.)) then nstep = NINT(1./onstep(1)) do n = 1, nstep do k = kte, kts, -1 sed_r(k) = vtrk(k)*rr(k) sed_n(k) = vtnrk(k)*nr(k) enddo k = kte odzq = 1./dzq(k) orho = 1./rho(k) qrten(k) = qrten(k) - sed_r(k)*odzq*onstep(1)*orho nrten(k) = nrten(k) - sed_n(k)*odzq*onstep(1)*orho rr(k) = MAX(R1, rr(k) - sed_r(k)*odzq*DT*onstep(1)) nr(k) = MAX(R2, nr(k) - sed_n(k)*odzq*DT*onstep(1)) do k = ksed1(1), kts, -1 odzq = 1./dzq(k) orho = 1./rho(k) qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) & *odzq*onstep(1)*orho nrten(k) = nrten(k) + (sed_n(k+1)-sed_n(k)) & *odzq*onstep(1)*orho rr(k) = MAX(R1, rr(k) + (sed_r(k+1)-sed_r(k)) & *odzq*DT*onstep(1)) nr(k) = MAX(R2, nr(k) + (sed_n(k+1)-sed_n(k)) & *odzq*DT*onstep(1)) enddo if (rr(kts).gt.R1*1000.) & pptrain = pptrain + sed_r(kts)*DT*onstep(1) enddo endif !+---+-----------------------------------------------------------------+ if (ANY(L_qc .eqv. .true.)) then do k = kte, kts, -1 sed_c(k) = vtck(k)*rc(k) sed_n(k) = vtnck(k)*nc(k) enddo do k = ksed1(5), kts, -1 odzq = 1./dzq(k) orho = 1./rho(k) qcten(k) = qcten(k) + (sed_c(k+1)-sed_c(k)) *odzq*orho ncten(k) = ncten(k) + (sed_n(k+1)-sed_n(k)) *odzq*orho rc(k) = MAX(R1, rc(k) + (sed_c(k+1)-sed_c(k)) *odzq*DT) nc(k) = MAX(10., nc(k) + (sed_n(k+1)-sed_n(k)) *odzq*DT) enddo !RC; next two lines we allow cloud water to contribute to surface precip. !RC; No need for the onstep() variable since cloud water falls so slowly. if (rc(kts) .gt. R1*1000.) then pptcloud = pptcloud + sed_c(kts)*DT endif endif !+---+-----------------------------------------------------------------+ if (ANY(L_qi .eqv. .true.)) then nstep = NINT(1./onstep(2)) do n = 1, nstep do k = kte, kts, -1 sed_i(k) = vtik(k)*ri(k) sed_n(k) = vtnik(k)*ni(k) enddo k = kte odzq = 1./dzq(k) orho = 1./rho(k) qiten(k) = qiten(k) - sed_i(k)*odzq*onstep(2)*orho niten(k) = niten(k) - sed_n(k)*odzq*onstep(2)*orho ri(k) = MAX(R1, ri(k) - sed_i(k)*odzq*DT*onstep(2)) ni(k) = MAX(R2, ni(k) - sed_n(k)*odzq*DT*onstep(2)) do k = ksed1(2), kts, -1 odzq = 1./dzq(k) orho = 1./rho(k) qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) & *odzq*onstep(2)*orho niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) & *odzq*onstep(2)*orho ri(k) = MAX(R1, ri(k) + (sed_i(k+1)-sed_i(k)) & *odzq*DT*onstep(2)) ni(k) = MAX(R2, ni(k) + (sed_n(k+1)-sed_n(k)) & *odzq*DT*onstep(2)) enddo if (ri(kts).gt.R1*1000.) & pptice = pptice + sed_i(kts)*DT*onstep(2) enddo endif !+---+-----------------------------------------------------------------+ if (ANY(L_qs .eqv. .true.)) then nstep = NINT(1./onstep(3)) do n = 1, nstep do k = kte, kts, -1 sed_s(k) = vtsk(k)*rs(k) enddo k = kte odzq = 1./dzq(k) orho = 1./rho(k) qsten(k) = qsten(k) - sed_s(k)*odzq*onstep(3)*orho rs(k) = MAX(R1, rs(k) - sed_s(k)*odzq*DT*onstep(3)) do k = ksed1(3), kts, -1 odzq = 1./dzq(k) orho = 1./rho(k) qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) & *odzq*onstep(3)*orho rs(k) = MAX(R1, rs(k) + (sed_s(k+1)-sed_s(k)) & *odzq*DT*onstep(3)) enddo if (rs(kts).gt.R1*1000.) & pptsnow = pptsnow + sed_s(kts)*DT*onstep(3) enddo endif !+---+-----------------------------------------------------------------+ if (ANY(L_qg .eqv. .true.)) then nstep = NINT(1./onstep(4)) do n = 1, nstep do k = kte, kts, -1 sed_g(k) = vtgk(k)*rg(k) sed_n(k) = vtngk(k)*ng(k) sed_b(k) = vtgk(k)*rb(k) enddo k = kte odzq = 1./dzq(k) orho = 1./rho(k) qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho ngten(k) = ngten(k) - sed_n(k)*odzq*onstep(4)*orho qbten(k) = qbten(k) - sed_b(k)*odzq*onstep(4) rg(k) = MAX(R1, rg(k) - sed_g(k)*odzq*DT*onstep(4)) ng(k) = MAX(R2, ng(k) - sed_n(k)*odzq*DT*onstep(4)) rb(k) = MAX(R1/rho(k)/rho_g(NRHG), rb(k) - sed_b(k)*odzq*DT*onstep(4)) do k = ksed1(4), kts, -1 odzq = 1./dzq(k) orho = 1./rho(k) qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) & *odzq*onstep(4)*orho ngten(k) = ngten(k) + (sed_n(k+1)-sed_n(k)) & *odzq*onstep(4)*orho qbten(k) = qbten(k) + (sed_b(k+1)-sed_b(k)) & *odzq*onstep(4) rg(k) = MAX(R1, rg(k) + (sed_g(k+1)-sed_g(k)) & *odzq*DT*onstep(4)) ng(k) = MAX(R2, ng(k) + (sed_n(k+1)-sed_n(k)) & *odzq*DT*onstep(4)) rb(k) = MAX(rg(k)/rho(k)/rho_g(NRHG), rb(k) + (sed_b(k+1)-sed_b(k)) & *odzq*DT*onstep(4)) enddo if (rg(kts).gt.R1*1000.) & pptgraul = pptgraul + sed_g(kts)*DT*onstep(4) enddo endif !+---+-----------------------------------------------------------------+ !.. Instantly melt any cloud ice into cloud water if above 0C and !.. instantly freeze any cloud water found below HGFR. !+---+-----------------------------------------------------------------+ if (.not. iiwarm) then do k = kts, kte xri = MAX(0.0, qi1d(k) + qiten(k)*DT) if ( (temp(k).gt. T_0) .and. (xri.gt. 0.0) ) then qcten(k) = qcten(k) + xri*odt ncten(k) = ncten(k) + ni1d(k)*odt qiten(k) = qiten(k) - xri*odt niten(k) = -ni1d(k)*odt tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-IFDRY) endif xrc = MAX(0.0, qc1d(k) + qcten(k)*DT) if ( (temp(k).lt. HGFR) .and. (xrc.gt. 0.0) ) then lfus2 = lsub - lvap(k) xnc = nc1d(k) + ncten(k)*DT qiten(k) = qiten(k) + xrc*odt niten(k) = niten(k) + xnc*odt qcten(k) = qcten(k) - xrc*odt ncten(k) = ncten(k) - xnc*odt tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-IFDRY) endif enddo endif !+---+-----------------------------------------------------------------+ !.. All tendencies computed, apply and pass back final values to parent. !+---+-----------------------------------------------------------------+ do k = kts, kte t1d(k) = t1d(k) + tten(k)*DT qv1d(k) = MAX(1.E-10, qv1d(k) + qvten(k)*DT) qc1d(k) = qc1d(k) + qcten(k)*DT nc1d(k) = MAX(2./rho(k), MIN(nc1d(k) + ncten(k)*DT, Nt_c_max)) if (is_aerosol_aware) then if (aer_init_opt .lt. 2) then nwfa1d(k) = MAX(11.1E6, MIN(9999.E6, & (nwfa1d(k)+nwfaten(k)*DT))) nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6, & (nifa1d(k)+nifaten(k)*DT))) if (wif_input_opt .eq. 2) then nbca1d(k) = MAX(5.55E6, MIN(9999.E6, & (nbca1d(k)+nbcaten(k)*DT))) else nbca1d(k) = 0.0 endif else nwfa1d(k) = MAX(0.0, (nwfa1d(k)+nwfaten(k)*DT)) nifa1d(k) = MAX(0.0, (nifa1d(k)+nifaten(k)*DT)) if (wif_input_opt .eq. 2) then nbca1d(k) = MAX(0.0, (nbca1d(k)+nbcaten(k)*DT)) else nbca1d(k) = 0.0 endif endif else nwfa1d(k) = MAX(11.1E6, MIN(9999.E6, & (nwfa1d(k)+nwfaten(k)*DT))) nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6, & (nifa1d(k)+nifaten(k)*DT))) nbca1d(k) = MAX(5.55E6, MIN(9999.E6, & (nbca1d(k)+nbcaten(k)*DT))) endif if (qc1d(k) .le. R1) then qc1d(k) = 0.0 nc1d(k) = 0.0 else !RC; Lognormal size adjustments xDc = ( ( (1./am_r)*(qc1d(k)/nc1d(k)) )**obmr ) if (xDc.lt. D0c) then xDc = D0c elseif (xDc.gt. D0r*2.) then xDc = D0r*2. endif nc1d(k) = MAX(2., MIN( DBLE(Nt_c_max/rho(k)), (1./am_r)*qc1d(k)*EXP(-4.5*ln_sigma(k)*ln_sigma(k))*(1. / ( xDc*EXP(-1.5*ln_sigma(k)*ln_sigma(k)) ))**3. )) !RC; end adjustments endif qi1d(k) = qi1d(k) + qiten(k)*DT ni1d(k) = MAX(R2/rho(k), ni1d(k) + niten(k)*DT) if (qi1d(k) .le. R1) then qi1d(k) = 0.0 ni1d(k) = 0.0 else lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi ilami = 1./lami xDi = (bm_i + mu_i + 1.) * ilami if (xDi.lt. 5.E-6) then lami = cie(2)/5.E-6 elseif (xDi.gt. 300.E-6) then lami = cie(2)/300.E-6 endif ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, & 999.D3/rho(k)) endif qr1d(k) = qr1d(k) + qrten(k)*DT nr1d(k) = MAX(R2/rho(k), nr1d(k) + nrten(k)*DT) if (qr1d(k) .le. R1) then qr1d(k) = 0.0 nr1d(k) = 0.0 else lamr = (am_r*crg(3)*org2*nr1d(k)/qr1d(k))**obmr mvd_r(k) = (3.0 + mu_r + 0.672) / lamr if (mvd_r(k) .gt. 2.5E-3) then mvd_r(k) = 2.5E-3 elseif (mvd_r(k) .lt. D0r*0.75) then mvd_r(k) = D0r*0.75 endif lamr = (3.0 + mu_r + 0.672) / mvd_r(k) nr1d(k) = crg(2)*org3*qr1d(k)*lamr**bm_r / am_r endif qs1d(k) = qs1d(k) + qsten(k)*DT if (qs1d(k) .le. R1) qs1d(k) = 0.0 qg1d(k) = qg1d(k) + qgten(k)*DT ng1d(k) = MAX(R2/rho(k), ng1d(k) + ngten(k)*DT) if (qg1d(k) .le. R1) then qg1d(k) = 0.0 ng1d(k) = 0.0 qb1d(k) = 0.0 else qb1d(k) = MAX(qg1d(k)/rho_g(NRHG), qb1d(k) + qbten(k)*DT) qb1d(k) = MIN(qg1d(k)/rho_g(1), qb1d(k)) idx_bg(k) = MAX(1,MIN(NINT(qg1d(k)/qb1d(k) *0.01)+1,NRHG)) if (.not. is_hail_aware) idx_bg(k) = idx_bg1 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng1d(k)/qg1d(k))**obmg mvd_g(k) = (3.0 + mu_g + 0.672) / lamg if (mvd_g(k) .gt. 25.4E-3) then mvd_g(k) = 25.4E-3 elseif (mvd_g(k) .lt. D0r) then mvd_g(k) = D0r endif lamg = (3.0 + mu_g + 0.672) / mvd_g(k) ng1d(k) = cgg(2,1)*ogg3*qg1d(k)*lamg**bm_g / am_g(idx_bg(k)) endif !RC; We set ilamc to 0 if qc 0. IMPLICIT NONE REAL, INTENT(IN):: XX DOUBLE PRECISION, PARAMETER:: STP = 2.5066282746310005D0 DOUBLE PRECISION, DIMENSION(6), PARAMETER:: & COF = (/76.18009172947146D0, -86.50532032941677D0, & 24.01409824083091D0, -1.231739572450155D0, & .1208650973866179D-2, -.5395239384953D-5/) DOUBLE PRECISION:: SER,TMP,X,Y INTEGER:: J X=XX Y=X TMP=X+5.5D0 TMP=(X+0.5D0)*LOG(TMP)-TMP SER=1.000000000190015D0 DO 11 J=1,6 Y=Y+1.D0 SER=SER+COF(J)/Y 11 CONTINUE GAMMLN=TMP+LOG(STP*SER/X) END FUNCTION GAMMLN ! (C) Copr. 1986-92 Numerical Recipes Software 2.02 !+---+-----------------------------------------------------------------+ REAL FUNCTION GAMMP(A,X) ! --- COMPUTES THE INCOMPLETE GAMMA FUNCTION P(A,X) ! --- SEE ABRAMOWITZ AND STEGUN 6.5.1 ! --- USES GCF,GSER IMPLICIT NONE REAL, INTENT(IN):: A,X REAL:: GAMMCF,GAMSER,GLN GAMMP = 0. IF((X.LT.0.) .OR. (A.LE.0.)) THEN PRINT *, 'BAD ARGUMENTS IN GAMMP' RETURN ELSEIF(X.LT.A+1.)THEN CALL GSER(GAMSER,A,X,GLN) GAMMP=GAMSER ELSE CALL GCF(GAMMCF,A,X,GLN) GAMMP=1.-GAMMCF ENDIF END FUNCTION GAMMP ! (C) Copr. 1986-92 Numerical Recipes Software 2.02 !+---+-----------------------------------------------------------------+ REAL FUNCTION WGAMMA(y) IMPLICIT NONE REAL, INTENT(IN):: y WGAMMA = EXP(GAMMLN(y)) END FUNCTION WGAMMA !+---+-----------------------------------------------------------------+ ! THIS FUNCTION CALCULATES THE LIQUID SATURATION VAPOR MIXING RATIO AS ! A FUNCTION OF TEMPERATURE AND PRESSURE ! REAL FUNCTION RSLF(P,T) IMPLICIT NONE REAL, INTENT(IN):: P, T REAL:: ESL,X REAL, PARAMETER:: C0= .611583699E03 REAL, PARAMETER:: C1= .444606896E02 REAL, PARAMETER:: C2= .143177157E01 REAL, PARAMETER:: C3= .264224321E-1 REAL, PARAMETER:: C4= .299291081E-3 REAL, PARAMETER:: C5= .203154182E-5 REAL, PARAMETER:: C6= .702620698E-8 REAL, PARAMETER:: C7= .379534310E-11 REAL, PARAMETER:: C8=-.321582393E-13 X=MAX(-80.,T-273.16) ! ESL=612.2*EXP(17.67*X/(T-29.65)) ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8))))))) ESL=MIN(ESL, P*0.15) ! Even with P=1050mb and T=55C, the sat. vap. pres only contributes to ~15% of total pres. RSLF=.622*ESL/(P-ESL) ! ALTERNATIVE ! ; Source: Murphy and Koop, Review of the vapour pressure of ice and ! supercooled water for atmospheric applications, Q. J. R. ! Meteorol. Soc (2005), 131, pp. 1539-1565. ! ESL = EXP(54.842763 - 6763.22 / T - 4.210 * ALOG(T) + 0.000367 * T ! + TANH(0.0415 * (T - 218.8)) * (53.878 - 1331.22 ! / T - 9.44523 * ALOG(T) + 0.014025 * T)) END FUNCTION RSLF !+---+-----------------------------------------------------------------+ ! THIS FUNCTION CALCULATES THE ICE SATURATION VAPOR MIXING RATIO AS A ! FUNCTION OF TEMPERATURE AND PRESSURE ! REAL FUNCTION RSIF(P,T) IMPLICIT NONE REAL, INTENT(IN):: P, T REAL:: ESI,X REAL, PARAMETER:: C0= .609868993E03 REAL, PARAMETER:: C1= .499320233E02 REAL, PARAMETER:: C2= .184672631E01 REAL, PARAMETER:: C3= .402737184E-1 REAL, PARAMETER:: C4= .565392987E-3 REAL, PARAMETER:: C5= .521693933E-5 REAL, PARAMETER:: C6= .307839583E-7 REAL, PARAMETER:: C7= .105785160E-9 REAL, PARAMETER:: C8= .161444444E-12 X=MAX(-80.,T-273.16) ESI=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8))))))) ESI=MIN(ESI, P*0.15) RSIF=.622*ESI/max(1.e-4,(P-ESI)) ! ALTERNATIVE ! ; Source: Murphy and Koop, Review of the vapour pressure of ice and ! supercooled water for atmospheric applications, Q. J. R. ! Meteorol. Soc (2005), 131, pp. 1539-1565. ! ESI = EXP(9.550426 - 5723.265/T + 3.53068*ALOG(T) - 0.00728332*T) END FUNCTION RSIF !+---+-----------------------------------------------------------------+ real function iceDeMott(tempc, qv, qvs, qvsi, rho, nifa) implicit none REAL, INTENT(IN):: tempc, qv, qvs, qvsi, rho, nifa !..Local vars REAL:: satw, sati, siw, p_x, si0x, dtt, dsi, dsw, dab, fc, hx REAL:: ntilde, n_in, nmax, nhat, mux, xni, nifa_cc REAL, PARAMETER:: p_c1 = 1000. REAL, PARAMETER:: p_rho_c = 0.76 REAL, PARAMETER:: p_alpha = 1.0 REAL, PARAMETER:: p_gam = 2. REAL, PARAMETER:: delT = 5. REAL, PARAMETER:: T0x = -40. REAL, PARAMETER:: Sw0x = 0.97 REAL, PARAMETER:: delSi = 0.1 REAL, PARAMETER:: hdm = 0.15 REAL, PARAMETER:: p_psi = 0.058707*p_gam/p_rho_c REAL, PARAMETER:: aap = 1. REAL, PARAMETER:: bbp = 0. REAL, PARAMETER:: y1p = -35. REAL, PARAMETER:: y2p = -25. REAL, PARAMETER:: rho_not0 = 101325./(287.05*273.15) !+---+ xni = 0.0 ! satw = qv/qvs ! sati = qv/qvsi ! siw = qvs/qvsi ! p_x = -1.0261+(3.1656e-3*tempc)+(5.3938e-4*(tempc*tempc)) & ! + (8.2584e-6*(tempc*tempc*tempc)) ! si0x = 1.+(10.**p_x) ! if (sati.ge.si0x .and. satw.lt.0.985) then ! dtt = delta_p (tempc, T0x, T0x+delT, 1., hdm) ! dsi = delta_p (sati, Si0x, Si0x+delSi, 0., 1.) ! dsw = delta_p (satw, Sw0x, 1., 0., 1.) ! fc = dtt*dsi*0.5 ! hx = min(fc+((1.-fc)*dsw), 1.) ! ntilde = p_c1*p_gam*((exp(12.96*(sati-1.1)))**0.3) / p_rho_c ! if (tempc .le. y1p) then ! n_in = ntilde ! elseif (tempc .ge. y2p) then ! n_in = p_psi*p_c1*exp(12.96*(sati-1.)-0.639) ! else ! if (tempc .le. -30.) then ! nmax = p_c1*p_gam*(exp(12.96*(siw-1.1)))**0.3/p_rho_c ! else ! nmax = p_psi*p_c1*exp(12.96*(siw-1.)-0.639) ! endif ! ntilde = MIN(ntilde, nmax) ! nhat = MIN(p_psi*p_c1*exp(12.96*(sati-1.)-0.639), nmax) ! dab = delta_p (tempc, y1p, y2p, aap, bbp) ! n_in = MIN(nhat*(ntilde/nhat)**dab, nmax) ! endif ! mux = hx*p_alpha*n_in*rho ! xni = mux*((6700.*nifa)-200.)/((6700.*5.E5)-200.) ! elseif (satw.ge.0.985 .and. tempc.gt.HGFR-273.15) then nifa_cc = MAX(0.5, nifa*RHO_NOT0*1.E-6/rho) ! xni = 3.*nifa_cc**(1.25)*exp((0.46*(-tempc))-11.6) ! [DeMott, 2015] xni = (5.94e-5*(-tempc)**3.33) & ! [DeMott, 2010] * (nifa_cc**((-0.0264*(tempc))+0.0033)) xni = xni*rho/RHO_NOT0 * 1000. ! endif iceDeMott = MAX(0., xni) end FUNCTION iceDeMott !+---+-----------------------------------------------------------------+ !..Newer research since Koop et al (2001) suggests that the freezing !.. rate should be lower than original paper, so J_rate is reduced !.. by two orders of magnitude. real function iceKoop(temp, qv, qvs, naero, dt) implicit none REAL, INTENT(IN):: temp, qv, qvs, naero, DT REAL:: mu_diff, a_w_i, delta_aw, log_J_rate, J_rate, prob_h, satw REAL:: xni xni = 0.0 satw = qv/qvs mu_diff = 210368.0 + (131.438*temp) - (3.32373E6/temp) & & - (41729.1*alog(temp)) a_w_i = exp(mu_diff/(R_uni*temp)) delta_aw = satw - a_w_i log_J_rate = -906.7 + (8502.0*delta_aw) & & - (26924.0*delta_aw*delta_aw) & & + (29180.0*delta_aw*delta_aw*delta_aw) log_J_rate = MIN(20.0, log_J_rate) J_rate = 10.**log_J_rate ! cm-3 s-1 prob_h = MIN(1.-exp(-J_rate*ar_volume*DT), 1.) if (prob_h .gt. 0.) then xni = MIN(prob_h*naero, 1000.E3) endif iceKoop = MAX(0.0, xni) end FUNCTION iceKoop !+---+-----------------------------------------------------------------+ !.. Helper routine for Phillips et al (2008) ice nucleation. Trude REAL FUNCTION delta_p (yy, y1, y2, aa, bb) IMPLICIT NONE REAL, INTENT(IN):: yy, y1, y2, aa, bb REAL:: dab, A, B, a0, a1, a2, a3 A = 6.*(aa-bb)/((y2-y1)*(y2-y1)*(y2-y1)) B = aa+(A*y1*y1*y1/6.)-(A*y1*y1*y2*0.5) a0 = B a1 = A*y1*y2 a2 = -A*(y1+y2)*0.5 a3 = A/3. if (yy.le.y1) then dab = aa else if (yy.ge.y2) then dab = bb else dab = a0+(a1*yy)+(a2*yy*yy)+(a3*yy*yy*yy) endif if (dab.lt.aa) then dab = aa endif if (dab.gt.bb) then dab = bb endif delta_p = dab END FUNCTION delta_p !+---+-----------------------------------------------------------------+ !ctrlL !+---+-----------------------------------------------------------------+ !..Compute _radiation_ effective radii of cloud water, ice, and snow. !.. These are entirely consistent with microphysics assumptions, not !.. constant or otherwise ad hoc as is internal to most radiation !.. schemes. Since only the smallest snowflakes should impact !.. radiation, compute from first portion of complicated Field number !.. distribution, not the second part, which is the larger sizes. !+---+-----------------------------------------------------------------+ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & & re_qc1d, re_qi1d, re_qs1d, kts, kte) IMPLICIT NONE !..Sub arguments INTEGER, INTENT(IN):: kts, kte REAL, DIMENSION(kts:kte), INTENT(IN):: & & t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d REAL, DIMENSION(kts:kte), INTENT(INOUT):: re_qc1d, re_qi1d, re_qs1d !..Local variables INTEGER:: k REAL, DIMENSION(kts:kte):: rho, rc, nc, ri, ni, rs REAL:: smo2, smob, smoc REAL:: tc0, loga_, a_, b_ DOUBLE PRECISION:: lamc, lami LOGICAL:: has_qc, has_qi, has_qs INTEGER:: inu_c real, dimension(15), parameter:: g_ratio = (/24,60,120,210,336, & & 504,720,990,1320,1716,2184,2730,3360,4080,4896/) DOUBLE PRECISION:: ln_sigma,ln_dn,sigma_d !RC has_qc = .false. has_qi = .false. has_qs = .false. re_qc1d(:) = RE_QC_BG re_qi1d(:) = RE_QI_BG re_qs1d(:) = RE_QS_BG do k = kts, kte rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) rc(k) = MAX(R1, qc1d(k)*rho(k)) nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max)) if (.NOT. is_aerosol_aware) nc(k) = Nt_c if (rc(k).gt.R1 .and. nc(k).gt.R2) has_qc = .true. ri(k) = MAX(R1, qi1d(k)*rho(k)) ni(k) = MAX(R2, ni1d(k)*rho(k)) if (ri(k).gt.R1 .and. ni(k).gt.R2) has_qi = .true. rs(k) = MAX(R1, qs1d(k)*rho(k)) if (rs(k).gt.R1) has_qs = .true. enddo !RC; New effective radii calculation for cloud. if (has_qc) then do k = kts, kte re_qc1d(k) = 2.49E-6 if (rc(k).le.R1 .or. nc(k).le.R2) CYCLE sigma_d = (rc(k)/nc(k))**(0.3333) ln_sigma = sigma_d*(-1.185E3) + 0.815 ln_sigma = MAX(ln_sigma , 0.2) ln_sigma = MIN(ln_sigma , 0.7) ln_dn = EXP(-3.0*0.5*ln_sigma**2)*(6.*rc(k)/(PI*nc(k)*1000.))**0.3333 re_qc1d(k) = MAX(2.51E-6, MIN(SNGL(0.5D0 * (ln_dn*EXP(2.5*ln_sigma*ln_sigma) ) ), 75.E-6)) enddo endif !/RC if (has_qi) then do k = kts, kte if (ri(k).le.R1 .or. ni(k).le.R2) CYCLE lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi re_qi1d(k) = MAX(2.51E-6, MIN(SNGL(0.5D0 * DBLE(3.+mu_i)/lami), 125.E-6)) enddo endif if (has_qs) then do k = kts, kte if (rs(k).le.R1) CYCLE tc0 = MIN(-0.1, t1d(k)-273.15) smob = rs(k)*oams !..All other moments based on reference, 2nd moment. If bm_s.ne.2, !.. then we must compute actual 2nd moment and use as reference. if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then smo2 = smob else loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s & & + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 & & + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s & & + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 & & + sa(10)*bm_s*bm_s*bm_s a_ = 10.0**loga_ b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s & & + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 & & + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s & & + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 & & + sb(10)*bm_s*bm_s*bm_s smo2 = (smob/a_)**(1./b_) endif !..Calculate bm_s+1 (th) moment. Useful for diameter calcs. loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & & + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & & + sa(10)*cse(1)*cse(1)*cse(1) a_ = 10.0**loga_ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) & & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) smoc = a_ * smo2**b_ re_qs1d(k) = MAX(5.01E-6, MIN(0.5*(smoc/smob), 999.E-6)) enddo endif end subroutine calc_effectRad !+---+-----------------------------------------------------------------+ !..Compute radar reflectivity assuming 10 cm wavelength radar and using !.. Rayleigh approximation. Only complication is melted snow/graupel !.. which we treat as water-coated ice spheres and use Uli Blahak's !.. library of routines. The meltwater fraction is simply the amount !.. of frozen species remaining from what initially existed at the !.. melting level interface. !+---+-----------------------------------------------------------------+ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & ng1d, qb1d, t1d, p1d, dBZ, kts, kte, ii, jj, ke_diag) IMPLICIT NONE !..Sub arguments INTEGER, INTENT(IN):: kts, kte, ii, jj, ke_diag REAL, DIMENSION(kts:kte), INTENT(IN):: & qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, ng1d, qb1d, t1d, p1d REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ ! REAL, DIMENSION(kts:kte), INTENT(INOUT):: vt_dBZ !..Local variables REAL, DIMENSION(kts:kte):: temp, pres, qv, rho, rhof REAL, DIMENSION(kts:kte):: rc, rr, nr, rs, rg, ng, rb DOUBLE PRECISION, DIMENSION(kts:kte):: ilamr, ilamg, N0_r, N0_g REAL, DIMENSION(kts:kte):: mvd_r REAL, DIMENSION(kts:kte):: smob, smo2, smoc, smoz REAL:: oM3, M0, Mrat, slam1, slam2, xDs REAL:: ils1, ils2, t1_vts, t2_vts, t3_vts, t4_vts REAL:: vtr_dbz_wt, vts_dbz_wt, vtg_dbz_wt REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel INTEGER, DIMENSION(kts:kte):: idx_bg DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamr, lamg REAL:: a_, b_, loga_, tc0 DOUBLE PRECISION:: fmelt_s, fmelt_g INTEGER:: i, k, k_0, kbot, n, ktop LOGICAL:: melti LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg DOUBLE PRECISION:: cback, x, eta, f_d REAL:: xslw1, ygra1, zans1 INTEGER :: k_0loop !+---+ do k = kts, kte dBZ(k) = -35.0 enddo !+---+-----------------------------------------------------------------+ !..Put column of data into local arrays. !+---+-----------------------------------------------------------------+ do k = kts, kte temp(k) = t1d(k) qv(k) = MAX(1.E-10, qv1d(k)) pres(k) = p1d(k) rho(k) = 0.622*pres(k)/(R*temp(k)*(qv(k)+0.622)) rhof(k) = SQRT(RHO_NOT/rho(k)) rc(k) = MAX(R1, qc1d(k)*rho(k)) if (qr1d(k) .gt. R1) then rr(k) = qr1d(k)*rho(k) nr(k) = MAX(R2, nr1d(k)*rho(k)) lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr ilamr(k) = 1./lamr N0_r(k) = nr(k)*org2*lamr**cre(2) mvd_r(k) = (3.0 + mu_r + 0.672) * ilamr(k) L_qr(k) = .true. else rr(k) = R1 nr(k) = R1 mvd_r(k) = 50.E-6 L_qr(k) = .false. endif if (qs1d(k) .gt. R2) then rs(k) = qs1d(k)*rho(k) L_qs(k) = .true. else rs(k) = R1 L_qs(k) = .false. endif if (qg1d(k) .gt. R2) then rg(k) = qg1d(k)*rho(k) ng(k) = MAX(R2, ng1d(k)*rho(k)) rb(k) = MAX(qg1d(k)/rho_g(NRHG), qb1d(k)) rb(k) = MIN(qg1d(k)/rho_g(1), rb(k)) idx_bg(k) = MAX(1,MIN(NINT(qg1d(k)/rb(k) *0.01)+1,NRHG)) if (.not. is_hail_aware) idx_bg(k) = idx_bg1 L_qg(k) = .true. else rg(k) = R1 ng(k) = R2 idx_bg(k) = idx_bg1 L_qg(k) = .false. endif enddo !+---+-----------------------------------------------------------------+ !..Calculate y-intercept, slope, and useful moments for snow. !+---+-----------------------------------------------------------------+ do k = kts, kte smo2(k) = 0. smob(k) = 0. smoc(k) = 0. smoz(k) = 0. enddo if ( ( ke_diag > kts .and. ANY(L_qs .eqv. .true.) ) .or. & (ke_diag == kts .and. L_qs(kts) .eqv. .true. ) ) then do k = kts, ke_diag ! kte if (.not. L_qs(k)) CYCLE tc0 = MIN(-0.1, temp(k)-273.15) smob(k) = rs(k)*oams !..All other moments based on reference, 2nd moment. If bm_s.ne.2, !.. then we must compute actual 2nd moment and use as reference. if (bm_s.gt.(2.0-1.e-3) .and. bm_s.lt.(2.0+1.e-3)) then smo2(k) = smob(k) else loga_ = sa(1) + sa(2)*tc0 + sa(3)*bm_s & & + sa(4)*tc0*bm_s + sa(5)*tc0*tc0 & & + sa(6)*bm_s*bm_s + sa(7)*tc0*tc0*bm_s & & + sa(8)*tc0*bm_s*bm_s + sa(9)*tc0*tc0*tc0 & & + sa(10)*bm_s*bm_s*bm_s a_ = 10.0**loga_ b_ = sb(1) + sb(2)*tc0 + sb(3)*bm_s & & + sb(4)*tc0*bm_s + sb(5)*tc0*tc0 & & + sb(6)*bm_s*bm_s + sb(7)*tc0*tc0*bm_s & & + sb(8)*tc0*bm_s*bm_s + sb(9)*tc0*tc0*tc0 & & + sb(10)*bm_s*bm_s*bm_s smo2(k) = (smob(k)/a_)**(1./b_) endif !..Calculate bm_s+1 (th) moment. Useful for diameter calcs. loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) & & + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 & & + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) & & + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 & & + sa(10)*cse(1)*cse(1)*cse(1) a_ = 10.0**loga_ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) & & + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) & & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) smoc(k) = a_ * smo2(k)**b_ !..Calculate bm_s*2 (th) moment. Useful for reflectivity. loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(3) & & + sa(4)*tc0*cse(3) + sa(5)*tc0*tc0 & & + sa(6)*cse(3)*cse(3) + sa(7)*tc0*tc0*cse(3) & & + sa(8)*tc0*cse(3)*cse(3) + sa(9)*tc0*tc0*tc0 & & + sa(10)*cse(3)*cse(3)*cse(3) a_ = 10.0**loga_ b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(3) + sb(4)*tc0*cse(3) & & + sb(5)*tc0*tc0 + sb(6)*cse(3)*cse(3) & & + sb(7)*tc0*tc0*cse(3) + sb(8)*tc0*cse(3)*cse(3) & & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(3)*cse(3)*cse(3) smoz(k) = a_ * smo2(k)**b_ enddo endif !+---+-----------------------------------------------------------------+ !..Calculate y-intercept, slope values for graupel. !+---+-----------------------------------------------------------------+ if (ANY(L_qg .eqv. .true.)) then do k = kte, kts, -1 lamg = (am_g(idx_bg(k))*cgg(3,1)*ogg2*ng(k)/rg(k))**obmg ilamg(k) = 1./lamg N0_g(k) = ng(k)*ogg2*lamg**cge(2,1) enddo else ilamg(:) = 0. N0_g(:) = 0. endif !+---+-----------------------------------------------------------------+ !..Locate K-level of start of melting (k_0 is level above). !+---+-----------------------------------------------------------------+ melti = .false. k_0 = kts do k = kte-1, kts, -1 if ( (temp(k).gt.273.15) .and. L_qr(k) & & .and. (L_qs(k+1).or.L_qg(k+1)) ) then k_0 = MAX(k+1, k_0) melti=.true. goto 195 endif enddo 195 continue ! Set loop limit for wet ice according to whether the full 3D field is needed or just k=1 k_0loop = Min(k_0, ke_diag+1) !+---+-----------------------------------------------------------------+ !..Do not do the so-called bright band if using the variable density !.. graupel-hail category since the density increases during melting and !.. will account for bright band behavior explicitly. !+---+-----------------------------------------------------------------+ ! if (is_hail_aware) melti = .false. !+---+-----------------------------------------------------------------+ !..Assume Rayleigh approximation at 10 cm wavelength. Rain (all temps) !.. and non-water-coated snow and graupel when below freezing are !.. simple. Integrations of m(D)*m(D)*N(D)*dD. !+---+-----------------------------------------------------------------+ do k = kts, ke_diag ! kte ze_rain(k) = 1.e-22 ze_snow(k) = 1.e-22 ze_graupel(k) = 1.e-22 if (L_qr(k)) ze_rain(k) = N0_r(k)*crg(4)*ilamr(k)**cre(4) if (L_qs(k)) ze_snow(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) & & * (am_s/900.0)*(am_s/900.0)*smoz(k) if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) & & * (am_g(idx_bg(k))/900.0)*(am_g(idx_bg(k))/900.0) & & * N0_g(k)*cgg(4,1)*ilamg(k)**cge(4,1) enddo !+---+-----------------------------------------------------------------+ !..Special case of melting ice (snow/graupel) particles. Assume the !.. ice is surrounded by the liquid water. Fraction of meltwater is !.. extremely simple based on amount found above the melting level. !.. Uses code from Uli Blahak (rayleigh_soak_wetgraupel and supporting !.. routines). !+---+-----------------------------------------------------------------+ if (.not. iiwarm .and. melti .and. k_0.ge.2) then do k = k_0loop-1, kts, -1 !..Reflectivity contributed by melting snow if (L_qs(k) .and. L_qs(k_0) ) then fmelt_s = MAX(0.05d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0)) eta = 0.d0 oM3 = 1./smoc(k) M0 = (smob(k)*oM3) Mrat = smob(k)*M0*M0*M0 slam1 = M0 * Lam0 slam2 = M0 * Lam1 do n = 1, nrbins x = am_s * xxDs(n)**bm_s call rayleigh_soak_wetgraupel (x, DBLE(ocms), DBLE(obms), & & fmelt_s, melt_outside_s, m_w_0, m_i_0, lamda_radar, & & CBACK, mixingrulestring_s, matrixstring_s, & & inclusionstring_s, hoststring_s, & & hostmatrixstring_s, hostinclusionstring_s) f_d = Mrat*(Kap0*DEXP(-slam1*xxDs(n)) & & + Kap1*(M0*xxDs(n))**mu_s * DEXP(-slam2*xxDs(n))) eta = eta + f_d * CBACK * simpson(n) * xdts(n) enddo ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta) endif !..Reflectivity contributed by melting graupel ! if (L_qg(k) .and. L_qg(k_0) ) then ! fmelt_g = MAX(0.05d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0)) ! eta = 0.d0 ! lamg = 1./ilamg(k) ! do n = 1, nrbins ! x = am_g(idx_bg(k)) * xxDg(n)**bm_g ! call rayleigh_soak_wetgraupel (x, DBLE(ocmg(idx_bg(k))), DBLE(obmg), & ! & fmelt_g, melt_outside_g, m_w_0, m_i_0, lamda_radar, & ! & CBACK, mixingrulestring_g, matrixstring_g, & ! & inclusionstring_g, hoststring_g, & ! & hostmatrixstring_g, hostinclusionstring_g) ! f_d = N0_g(k)*xxDg(n)**mu_g * DEXP(-lamg*xxDg(n)) ! eta = eta + f_d * CBACK * simpson(n) * xdtg(n) ! enddo ! ze_graupel(k) = SNGL(lamda4 / (pi5 * K_w) * eta) ! endif enddo endif do k = ke_diag, kts, -1 dBZ(k) = 10.*log10((ze_rain(k)+ze_snow(k)+ze_graupel(k))*1.d18) enddo !..Reflectivity-weighted terminal velocity (snow, rain, graupel, mix). ! do k = kte, kts, -1 ! vt_dBZ(k) = 1.E-3 ! if (rs(k).gt.R2) then ! Mrat = smob(k) / smoc(k) ! ils1 = 1./(Mrat*Lam0 + fv_s) ! ils2 = 1./(Mrat*Lam1 + fv_s) ! t1_vts = Kap0*csg(5)*ils1**cse(5) ! t2_vts = Kap1*Mrat**mu_s*csg(11)*ils2**cse(11) ! ils1 = 1./(Mrat*Lam0) ! ils2 = 1./(Mrat*Lam1) ! t3_vts = Kap0*csg(6)*ils1**cse(6) ! t4_vts = Kap1*Mrat**mu_s*csg(12)*ils2**cse(12) ! vts_dbz_wt = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts) ! if (temp(k).ge.273.15 .and. temp(k).lt.275.15) then ! vts_dbz_wt = vts_dbz_wt*1.5 ! elseif (temp(k).ge.275.15) then ! vts_dbz_wt = vts_dbz_wt*2.0 ! endif ! else ! vts_dbz_wt = 1.E-3 ! endif ! if (rr(k).gt.R1) then ! lamr = 1./ilamr(k) ! vtr_dbz_wt = rhof(k)*av_r*crg(13)*(lamr+fv_r)**(-cre(13)) & ! & / (crg(4)*lamr**(-cre(4))) ! else ! vtr_dbz_wt = 1.E-3 ! endif ! if (rg(k).gt.R2) then ! lamg = 1./ilamg(k) ! vtg_dbz_wt = rhof(k)*av_g(idx_bg(k))*cgg(5,idx_bg(k))*lamg**(-cge(5,idx_bg(k))) & ! & / (cgg(4,1)*lamg**(-cge(4,1))) ! else ! vtg_dbz_wt = 1.E-3 ! endif ! vt_dBZ(k) = (vts_dbz_wt*ze_snow(k) + vtr_dbz_wt*ze_rain(k) & ! & + vtg_dbz_wt*ze_graupel(k)) & ! & / (ze_rain(k)+ze_snow(k)+ze_graupel(k)) ! enddo end subroutine calc_refl10cm ! !+---+-----------------------------------------------------------------+ ! real function theta_e(pres_Pa,temp_K,w_non,tlcl_K) !.. !.. The following code was based on Bolton (1980) eqn #43 !.. and claims to have 0.3 K maximum error within -35 < T < 35 C !.. pres_Pa = Pressure in Pascals !.. temp_K = Temperature in Kelvin !.. w_non = mixing ratio (non-dimensional = kg/kg) !.. tlcl_K = Temperature at Lifting Condensation Level (K) !.. IMPLICIT NONE real:: pres_Pa, temp_K, w_non, tlcl_K real:: pp, tt, rr, tlc, power, xx, p1, p2 !+---+ pp = pres_Pa tt = temp_K rr = w_non + 1.e-8 tlc = tlcl_K power=(0.2854*(1.0 - (0.28*rr) )) xx = tt * (100000.0/pp)**power p1 = (3.376/tlc) - 0.00254 p2 = (rr*1000.0) * (1.0 + 0.81*rr) theta_e = xx * exp(p1*p2) return end ! !+---+-----------------------------------------------------------------+ ! real function t_lcl(temp_K,tdew_K) !.. !.. The following code was based on Bolton (1980) eqn #15 !.. and claims to have 0.1 K maximum error within -35 < T < 35 C !.. temp_K = Temperature in Kelvin !.. tdew_K = Dewpoint T at Lifting Condensation Level (K) !.. IMPLICIT NONE real:: temp_K, tdew_K real:: tt, tttd, denom !+---+ tt = temp_K tttd= tdew_K denom= ( 1.0/(tttd-56.0) ) + (log(tt/tttd)/800.) t_lcl = ( 1.0 / denom ) + 56.0 return end ! !+---+-----------------------------------------------------------------+ ! real function t_dew(pres_Pa,w_non) !.. !.. pres_Pa = Pressure in Pascals !.. w_non = mixing ratio (non-dimensional = kg/kg) !.. IMPLICIT NONE real:: pres_Pa, w_non real:: p, RR, ES, ESLN !+---+ p = pres_Pa RR=w_non+1e-8 ES=P*RR/(.622+RR) ESLN=LOG(ES) T_Dew=(35.86*ESLN-4947.2325)/(ESLN-23.6837) return end ! !+---+-----------------------------------------------------------------+ ! real function theta_wetb(thetae_K) !.. !.. Eqn below was gotten from polynomial fit to data in !.. Smithsonian Meteorological Tables showing Theta-e !.. and Theta-w !.. IMPLICIT NONE real:: thetae_K real:: x, answer !+---+ real*8 c(0:6), d(0:6) data c/-1.00922292e-10, -1.47945344e-8, -1.7303757e-6 & & ,-0.00012709, 1.15849867e-6, -3.518296861e-9 & & ,3.5741522e-12/ data d/0.00000000, -3.5223513e-10, -5.7250807e-8 & & ,-5.83975422e-6, 4.72445163e-8, -1.13402845e-10 & & ,8.729580402e-14/ x=min(475.0,thetae_K) if( x .le. 335.5 ) then answer = c(0)+x*(c(1)+x*(c(2)+x*(c(3)+x*(c(4)+x*(c(5)+ & & x*c(6) ))))) else answer = d(0)+x*(d(1)+x*(d(2)+x*(d(3)+x*(d(4)+x*(d(5)+ & & x*d(6) ))))) endif theta_wetb = answer + 273.15 return end ! !+---+-----------------------------------------------------------------+ ! real function compT_fr_The(thelcl_K,pres_Pa) !.. !.. pres_Pa = Pressure in Pascals !.. thelcl = Theta-e at LCL (units in Kelvin) !.. !.. Temperature (K) is returned given Theta-e at LCL !.. and a pressure. This describes a moist-adiabat. !.. This temperature is the parcel temp at level Pres !.. along moist adiabat described by theta-e. !.. IMPLICIT NONE real:: thelcl_K, pres_Pa real:: guess, epsilon, w1, w2, tenu, tenup, cor, thwlcl_K integer:: iter !+---+ guess= (thelcl_K - 0.5 * ( max(thelcl_K-270., 0.))**1.05) & & * (pres_Pa/100000.0)**.2 epsilon=0.01 do iter=1,100 w1 = rslf(pres_Pa,guess) w2 = rslf(pres_Pa,guess+1.) tenu = theta_e(pres_Pa,guess,w1,guess) tenup = theta_e(pres_Pa,guess+1,w2,guess+1.) cor = (thelcl_K - tenu) / (tenup - tenu) guess = guess + cor if( (cor.lt.epsilon) .and. (-cor.lt.epsilon) ) then compT_fr_The=guess return endif enddo ! print*, ' convergence not reached ' thwlcl_K=theta_wetb(thelcl_K) compT_fr_The = thwlcl_K*((pres_Pa/100000.0)**0.286) return end !+---+-----------------------------------------------------------------+ !+---+-----------------------------------------------------------------+ END MODULE module_mp_rcon !+---+-----------------------------------------------------------------+