From a7f9269c42f34f48372ed21bd0bcf55e75f9d7a0 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Wed, 13 Jun 2018 07:17:37 -0600 Subject: [PATCH] Initial version of aerosol-aware Thompson MP scheme with modifications for faster initialization (physics/module_mp_thompson_hrrr.F90), dependencies (physics/module_mp_thompson_hrrr_radar.F90) and CCPP entry points (physics/mp_thompson_hrrr.F90) --- physics/module_mp_thompson_hrrr.F90 | 5998 +++++++++++++++++++++ physics/module_mp_thompson_hrrr_radar.F90 | 612 +++ physics/mp_thompson_hrrr.F90 | 557 ++ 3 files changed, 7167 insertions(+) create mode 100644 physics/module_mp_thompson_hrrr.F90 create mode 100644 physics/module_mp_thompson_hrrr_radar.F90 create mode 100644 physics/mp_thompson_hrrr.F90 diff --git a/physics/module_mp_thompson_hrrr.F90 b/physics/module_mp_thompson_hrrr.F90 new file mode 100644 index 000000000..0fb3b214d --- /dev/null +++ b/physics/module_mp_thompson_hrrr.F90 @@ -0,0 +1,5998 @@ +#define DEBUG_AEROSOLS + +!+---+-----------------------------------------------------------------+ +!.. This subroutine computes the moisture tendencies of water vapor, +!.. cloud droplets, rain, cloud ice (pristine), snow, and graupel. +!.. Prior to WRFv2.2 this code was based on Reisner et al (1998), but +!.. few of those pieces remain. A complete description is now found in +!.. Thompson, G., P. R. Field, R. M. Rasmussen, and W. D. Hall, 2008: +!.. Explicit Forecasts of winter precipitation using an improved bulk +!.. microphysics scheme. Part II: Implementation of a new snow +!.. parameterization. Mon. Wea. Rev., 136, 5095-5115. +!.. Prior to WRFv3.1, this code was single-moment rain prediction as +!.. described in the reference above, but in v3.1 and higher, the +!.. scheme is two-moment rain (predicted rain number concentration). +!.. +!.. Beginning with WRFv3.6, this is also the "aerosol-aware" scheme as +!.. described in Thompson, G. and T. Eidhammer, 2014: A study of +!.. aerosol impacts on clouds and precipitation development in a large +!.. winter cyclone. J. Atmos. Sci., 71, 3636-3658. Setting WRF +!.. namelist option mp_physics=8 utilizes the older one-moment cloud +!.. water with constant droplet concentration set as Nt_c (found below) +!.. while mp_physics=28 uses double-moment cloud droplet number +!.. concentration, which is not permitted to exceed Nt_c_max below. +!.. +!.. Most importantly, users may wish to modify the prescribed number of +!.. cloud droplets (Nt_c; see guidelines mentioned below). Otherwise, +!.. users may alter the rain and graupel size distribution parameters +!.. to use exponential (Marshal-Palmer) or generalized gamma shape. +!.. The snow field assumes a combination of two gamma functions (from +!.. Field et al. 2005) and would require significant modifications +!.. throughout the entire code to alter its shape as well as accretion +!.. rates. Users may also alter the constants used for density of rain, +!.. graupel, ice, and snow, but the latter is not constant when using +!.. Paul Field's snow distribution and moments methods. Other values +!.. users can modify include the constants for mass and/or velocity +!.. power law relations and assumed capacitances used in deposition/ +!.. sublimation/evaporation/melting. +!.. Remaining values should probably be left alone. +!.. +!..Author: Greg Thompson, NCAR-RAL, gthompsn@ucar.edu, 303-497-2805 +!..Last modified: 01 Aug 2016 Aerosol additions to v3.5.1 code 9/2013 +!.. Cloud fraction additions 11/2014 part of pre-v3.7 +!..Imported in CCPP by: Dom Heinzeller, NOAA/ESRL/GS, dom.heinzeller@noaa.gov +!..Last modified: 22 May 2018 Initial import of HRRR operational version +!+---+-----------------------------------------------------------------+ +!wrft:model_layer:physics +!+---+-----------------------------------------------------------------+ +! +MODULE module_mp_thompson_hrrr + + USE machine, ONLY : kind_phys + + USE module_mp_thompson_hrrr_radar + + IMPLICIT NONE + + LOGICAL, PARAMETER, PRIVATE:: iiwarm = .false. + ! CCPP + !LOGICAL, PRIVATE:: is_aerosol_aware = .false. + LOGICAL, PARAMETER, PRIVATE:: dustyIce = .true. + LOGICAL, PARAMETER, PRIVATE:: homogIce = .true. +!tgs - test +! LOGICAL, PARAMETER, PRIVATE:: dustyIce = .false. +! LOGICAL, PARAMETER, PRIVATE:: homogIce = .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, PARAMETER, PRIVATE:: rho_s = 100.0 + REAL, PARAMETER, PRIVATE:: rho_g = 500.0 + REAL, PARAMETER, PRIVATE:: rho_i = 890.0 + +!..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 + +!..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.E4 + REAL, PARAMETER, PRIVATE:: gonv_max = 3.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, PARAMETER, PRIVATE:: am_g = PI*rho_g/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 = 442.0 + REAL, PARAMETER, PRIVATE:: bv_g = 0.89 + REAL, PARAMETER, PRIVATE:: av_i = 1847.5 + 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 = 200.E-6 + REAL, PARAMETER, PRIVATE:: D0g = 250.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 = 28 + INTEGER, PARAMETER, PRIVATE:: ntb_g = 28 + INTEGER, PARAMETER, PRIVATE:: ntb_g1 = 28 + 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-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-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.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,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, & + 1.e7/) + +!..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(18), PRIVATE:: cse, csg + REAL, PRIVATE:: oams, obms, ocms + REAL, DIMENSION(12), PRIVATE:: cge, cgg + REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, oamg, obmg, 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 thompson_init(hgt, nwfa2d, nwfa, nifa, dx, dy, & + is_start, & + ids, ide, jds, jde, kds, kde, & + ims, ime, jms, jme, kms, kme, & + its, ite, jts, jte, kts, kte, & + mpicomm, mpirank, mpiroot, threads) + + 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 aerosol-aware scheme + + REAL, DIMENSION(ims:ime,kms:kme,jms:jme), OPTIONAL, INTENT(INOUT) :: nwfa, nifa + REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(INOUT) :: nwfa2d + REAL, DIMENSION(ims:ime,jms:jme), OPTIONAL, INTENT(IN) :: DX, DY + LOGICAL, OPTIONAL, INTENT(IN) :: is_start + INTEGER, INTENT(IN) :: mpicomm, mpirank, mpiroot + INTEGER, INTENT(IN) :: threads + + + INTEGER:: i, j, k, l, m, n + REAL:: h_01, niIN3, niCCN3, max_test + LOGICAL:: is_aerosol_aware, micro_init, has_CCN, has_IN + real :: stime, etime +#ifdef SION + INTEGER :: ierr + LOGICAL :: precomputed_tables +#else + LOGICAL, PARAMETER :: precomputed_tables = .FALSE. +#endif + + is_aerosol_aware = .FALSE. + micro_init = .FALSE. + has_CCN = .FALSE. + has_IN = .FALSE. + + if (PRESENT(nwfa2d) .AND. PRESENT(nwfa) .AND. PRESENT(nifa)) is_aerosol_aware = .TRUE. + + if_aerosol_aware: if (is_aerosol_aware) then + + if (.NOT.PRESENT(is_start)) then + write(*,*) "thompson_init: optional argument is_start required if is_aerosol_aware is true" + return + end if + +!..Check for existing aerosol data, both CCN and IN aerosols. If missing +!.. fill in just a basic vertical profile, somewhat boundary-layer following. + + max_test = MAXVAL ( nwfa(its:ite,:,jts:jte) ) + + if (max_test .lt. eps) then + write(*,*) ' Apparently there are no initial CCN aerosols.' + write(*,*) ' checked column at point (i,j) = ', its,jts +! do j = jts, min(jde-1,jte) +! do i = its, min(ide-1,ite) +! tgs for FIM + do j = jts, jte + do i = its, 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) + 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(*,*) ' Apparently initial CCN aerosols are present.' + write(*,*) ' column sum at point (i,j) = ', its,jts, SUM(nwfa(its,:,jts)) + endif + + +!tgs for FIM + max_test = MAXVAL ( nifa(its:ite,:,jts:jte) ) + + if (max_test .lt. eps) then + write(*,*) ' Apparently there are no initial IN aerosols.' + write(*,*) ' checked column at point (i,j) = ', its,jts +!tgs for FIM +! do j = jts, min(jde-1,jte) +! do i = its, min(ide-1,ite) + do j = jts, jte + do i = its, 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(*,*) ' Apparently initial IN aerosols are present.' + write(*,*) ' column sum at point (i,j) = ', its,jts, SUM(nifa(its,:,jts)) + endif + +!..Capture initial state lowest level CCN aerosol data in 2D array. + +! do j = jts, min(jde-1,jte) +! do i = its, min(ide-1,ite) +! nwfa2d(i,j) = nwfa(i,kts,j) +! enddo +! enddo + +! DH* +! Allow for external aerosol surface emission rates, fallback to original +! approach by Greg if none are provided. +! *DH +!..Scale the lowest level aerosol data into an emissions rate. This is +!.. very far from ideal, but need higher emissions where larger amount +!.. of existing and lesser emissions where not already lots of aerosols +!.. for first-order simplistic approach. Later, proper connection to +!.. emission inventory would be better, but, for now, scale like this: +!.. where: Nwfa=50 per cc, emit 0.875E4 aerosols per kg per second +!.. Nwfa=500 per cc, emit 0.875E5 aerosols per kg per second +!.. Nwfa=5000 per cc, emit 0.875E6 aerosols per kg per second +!.. for a grid with 20km spacing and scale accordingly for other spacings. + + max_test = MAXVAL ( nwfa2d(its:ite,jts:jte) ) + + if (is_start .and. max_test.lt.eps) then + write(*,*) ' Apparently CCN emission rates are not specified.' + write(*,*) ' checked column at point (i,j) = ', its,jts + write(*,*) ' calculate emission rates from initial CCN distribution' + !if (SQRT(DX*DY)/20000.0 .ge. 1.0) then + ! h_01 = 0.875 + !else + ! h_01 = (0.875 + 0.125*((20000.-SQRT(DX*DY))/16000.)) * SQRT(DX*DY)/20000. + !endif + !write(*,*) ' aerosol surface flux emission scale factor is: ', h_01 +!tgs - for FIM +! do j = jts, min(jde-1,jte) +! do i = its, min(ide-1,ite) + do j = jts, jte + do i = its, ite + if (SQRT(DX(i,j)*DY(i,j))/20000.0 .ge. 1.0) then + h_01 = 0.875 + else + h_01 = (0.875 + 0.125*((20000.-SQRT(DX(i,j)*DY(i,j)))/16000.)) * SQRT(DX(i,j)*DY(i,j))/20000. + endif + nwfa2d(i,j) = 10.0**(LOG10(nwfa(i,kts,j)*1.E-6)-3.69897) + nwfa2d(i,j) = nwfa2d(i,j)*h_01 * 1.E6 + enddo + enddo + endif + + endif if_aerosol_aware + + +!..Allocate space for lookup tables (J. Michalakes 2009Jun08). + + if (.NOT. ALLOCATED(tcg_racg) ) then + ALLOCATE(tcg_racg(ntb_g1,ntb_g,ntb_r1,ntb_r)) + micro_init = .TRUE. + endif + + if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tcr_gacr)) ALLOCATE(tcr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tmg_gacr)) ALLOCATE(tmg_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tnr_racg)) ALLOCATE(tnr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r)) + if (.NOT. ALLOCATED(tnr_gacr)) ALLOCATE(tnr_gacr(ntb_g1,ntb_g,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: 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 * 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) = cse(16) + mu_s + 1. + cse(18) = bv_s + mu_s + 3. + do n = 1, 18 + 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(5) = bm_g*2. + mu_g + bv_g + 1. + cge(6) = bm_g + mu_g + bv_g + 1. + cge(7) = bm_g + mu_g + bv_g + 2. + cge(8) = bm_g + mu_g + bv_g + 3. + cge(9) = mu_g + bv_g + 3. + cge(10) = mu_g + 2. + cge(11) = 0.5*(bv_g + 5. + 2.*mu_g) + cge(12) = 0.5*(bv_g + 5.) + mu_g + do n = 1, 12 + cgg(n) = WGAMMA(cge(n)) + enddo + oamg = 1./am_g + obmg = 1./bm_g + ocmg = oamg**obmg + oge1 = 1./cge(1) + ogg1 = 1./cgg(1) + ogg2 = 1./cgg(2) + ogg3 = 1./cgg(3) + +!+---+-----------------------------------------------------------------+ +!..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) + 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) + 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) = 5.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. +!+---+-----------------------------------------------------------------+ + +#ifdef SION + call cpu_time(stime) + call readwrite_tables("read", mpicomm, mpirank, mpiroot, ierr) + call cpu_time(etime) + if (mpirank==mpiroot) print '("Reading and broadcasting precomputed Thompson tables took ",f10.3," seconds.")', etime-stime + if (ierr==0) then + precomputed_tables = .true. + if (mpirank==mpiroot) write(0,*) "Read Thompson tables from disk, skip calculation" + else + precomputed_tables = .false. + if (mpirank==mpiroot) write(0,*) "An error occurred reading Thompson tables from disk, recalculate" + end if +#endif + + precomputed_tables_1: if (.not.precomputed_tables) then + + call cpu_time(stime) + + do m = 1, ntb_r + do k = 1, ntb_r1 + do j = 1, ntb_g + do i = 1, ntb_g1 + tcg_racg(i,j,k,m) = 0.0d0 + tmr_racg(i,j,k,m) = 0.0d0 + tcr_gacr(i,j,k,m) = 0.0d0 + tmg_gacr(i,j,k,m) = 0.0d0 + tnr_racg(i,j,k,m) = 0.0d0 + tnr_gacr(i,j,k,m) = 0.0d0 + 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 + + if (mpirank==mpiroot) WRITE (*,*)'CREATING MICROPHYSICS LOOKUP TABLES ... ' + if (mpirank==mpiroot) WRITE (*,*) ' using: mu_c=',mu_c,' mu_i=',mu_i, & + ' mu_r=',mu_r,' mu_g=',mu_g + +!..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. + ! This computation is cheap compared to the others below, and + ! doing it always ensures that the correct data is in the SIONlib + ! file containing the precomputed tables *DH + !if (is_aerosol_aware) then + if (mpirank==mpiroot) WRITE (*,*) ' calling table_ccnAct routine' + call table_ccnAct + !endif + +!..Collision efficiency between rain/snow and cloud water. + if (mpirank==mpiroot) WRITE (*,*)' creating qc collision eff tables' + call table_Efrw + call table_Efsw + +!..Drop evaporation. + WRITE(*,*) ' creating rain evap table' + call table_dropEvap + + call cpu_time(etime) + if (mpirank==mpiroot) print '("Calculating Thompson tables part 1 took ",f10.3," seconds.")', etime-stime + + end if precomputed_tables_1 + +!..Initialize various constants for computing radar reflectivity. + call cpu_time(stime) + 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 + xbm_g = bm_g + xmu_g = mu_g + call radar_init + call cpu_time(etime) + if (mpirank==mpiroot) print '("Calling radar_init took ",f10.3," seconds.")', etime-stime + + + if_not_iiwarm: if (.not. iiwarm) then + + precomputed_tables_2: if (.not.precomputed_tables) then + + call cpu_time(stime) + +!$OMP parallel num_threads(threads) + +!$OMP sections + +!$OMP section +!..Rain collecting graupel & graupel collecting rain. + WRITE (*,*) ' creating rain collecting graupel table' + call qr_acr_qg + +!$OMP section +!..Rain collecting snow & snow collecting rain. + WRITE (*,*) ' creating rain collecting snow table' + call qr_acr_qs + +!$OMP section +!..Cloud water and rain freezing (Bigg, 1953). + WRITE (*,*) ' creating freezing of water drops table' + call freezeH2O + +!$OMP section +!..Conversion of some ice mass into snow category. + WRITE (*,*) ' creating ice converting to snow table' + call qi_aut_qs + +!$OMP end sections + +!$OMP end parallel + + call cpu_time(etime) + if (mpirank==mpiroot) print '("Calculating Thompson tables part 2 took ",f10.3," seconds.")', etime-stime + +#ifdef SION + call cpu_time(stime) + call readwrite_tables("write", mpicomm, mpirank, mpiroot, ierr) + if (ierr/=0) then + write(0,*) "An error occurred writing Thompson tables to disk" + stop 1 + end if + call cpu_time(etime) + if (mpirank==mpiroot) print '("Writing Thompson tables took ",f10.3," seconds.")', etime-stime +#endif + + end if precomputed_tables_2 + + endif if_not_iiwarm + + WRITE (*,*) ' ... DONE microphysical lookup tables' + + endif if_micro_init + + END SUBROUTINE thompson_init +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!..This is a wrapper routine designed to transfer values from 3D to 1D. +!+---+-----------------------------------------------------------------+ + SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & + nwfa, nifa, nwfa2d, & + tt, th, pii, & + p, w, dz, dt_in, itimestep, & + RAINNC, RAINNCV, & + SNOWNC, SNOWNCV, & + ICENC, ICENCV, & + GRAUPELNC, GRAUPELNCV, SR, & +#if ( WRF_CHEM == 1 ) + rainprod, evapprod, & +#endif + refl_10cm, vt_dbz_wt, & ! dcd + diagflag, 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, & + errmsg, errflg) ! tile dims +#ifdef DEBUG_AEROSOLS + use mpi +#endif + + 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 + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & + tt, th + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(IN):: & + pii + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & + nc, nwfa, nifa + REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(OUT):: & + 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 + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: & + 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, & + ICENC, ICENCV, & + GRAUPELNC, GRAUPELNCV +! dcd added vt_dbz_wt + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & + refl_10cm, vt_dbz_wt + REAL, INTENT(IN):: dt_in + INTEGER, INTENT(IN):: itimestep + +!..Local variables + REAL, DIMENSION(kts:kte):: & + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + nr1d, nc1d, nwfa1d, nifa1d, & + t1d, p1d, w1d, dz1d, rho, dBZ, vt_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 + REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max + REAL:: nwfa1 + INTEGER:: i, j, k + 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 + LOGICAL :: first_time_step ! dcd + ! CCPP + ! CCPP error handling + character(len=*), optional, intent( out) :: errmsg + integer, optional, intent( out) :: errflg + LOGICAL :: is_aerosol_aware +#ifdef DEBUG_AEROSOLS + integer :: mpirank, mpisize, impi, ierr + real :: nc1d_in_debug, nwfa1d_in_debug, nifa1d_in_debug, nwfa1_in_debug + real :: nc1d_out_debug, nwfa1d_out_debug, nifa1d_out_debug, nwfa1_out_debug + + call MPI_COMM_RANK(MPI_COMM_WORLD,mpirank,ierr) + call MPI_COMM_SIZE(MPI_COMM_WORLD,mpisize,ierr) + + nc1d_in_debug = 0.0 + nwfa1d_in_debug = 0.0 + nifa1d_in_debug = 0.0 + nwfa1_in_debug = 0.0 + nc1d_out_debug = 0.0 + nwfa1d_out_debug = 0.0 + nifa1d_out_debug = 0.0 + nwfa1_out_debug = 0.0 +#endif + + ! CCPP + if (present(errmsg)) errmsg = '' + if (present(errflg)) errflg = 0 + + if ( (present(tt) .and. (present(th) .or. present(pii))) .or. & + (.not.present(tt) .and. .not.(present(th) .and. present(pii))) ) then + if (present(errmsg)) then + write(errmsg, '(a)') 'Logic error in mp_gt_driver: provide either tt or th+pii' + else + write(*,'(a)') 'Logic error in mp_gt_driver: provide either tt or th+pii' + end if + if (present(errflg)) then + errflg = 1 + return + else + stop + end if + end if + + if (present(nc) .and. present(nwfa) .and. present(nifa) .and. present(nwfa2d)) then + is_aerosol_aware = .true. + else + is_aerosol_aware = .false. + end if + +!+---+ + first_time_step = (itimestep .eq. 1) ! dcd + i_start = its + j_start = jts +!tgs - for FIM + i_end = ite + j_end = jte +! 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 + + if (.NOT. is_aerosol_aware .AND. PRESENT(nc) .AND. PRESENT(nwfa) & + .AND. PRESENT(nifa) .AND. PRESENT(nwfa2d)) then + write(*,*)'WARNING, nc-nwfa-nifa-nwfa2d present but is_aerosol_aware is FALSE' + endif + + j_loop: do j = j_start, j_end + i_loop: do i = i_start, i_end + + pptrain = 0. + pptsnow = 0. + pptgraul = 0. + pptice = 0. + RAINNCV(i,j) = 0. + IF ( PRESENT (snowncv) ) THEN + SNOWNCV(i,j) = 0. + ENDIF + IF ( PRESENT (icencv) ) THEN + ICENCV(i,j) = 0. + ENDIF + IF ( PRESENT (graupelncv) ) THEN + GRAUPELNCV(i,j) = 0. + ENDIF + SR(i,j) = 0. + + do k = kts, kte + if (present(tt)) then + t1d(k) = tt(i,k,j) + else + t1d(k) = th(i,k,j)*pii(i,k,j) + end if + 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) + 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) + enddo + nwfa1 = nwfa2d(i,j) + else + do k = kts, kte + rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) + nc1d(k) = Nt_c/rho(k) + nwfa1d(k) = 11.1E6/rho(k) + nifa1d(k) = naIN1*0.01/rho(k) + enddo + nwfa1 = 11.1E6 + endif +#ifdef DEBUG_AEROSOLS + nc1d_in_debug = nc1d_in_debug + sum(nc1d(:)) + nwfa1d_in_debug = nwfa1d_in_debug + sum(nwfa1d(:)) + nifa1d_in_debug = nifa1d_in_debug + sum(nifa1d(:)) + nwfa1_in_debug = nwfa1_in_debug + nwfa1 +#endif + call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dz1d, & + pptrain, pptsnow, pptgraul, pptice, & +#if ( WRF_CHEM == 1 ) + rainprod1d, evapprod1d, & +#endif + kts, kte, dt, i, j, is_aerosol_aware) + + 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 + RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice + IF ( PRESENT(snowncv) .AND. PRESENT(snownc) ) THEN + ! Add ice to snow if separate ice not present (as in WRF) + IF ( .NOT.PRESENT(icencv) .OR. .NOT.PRESENT(icenc) ) THEN + SNOWNCV(i,j) = pptsnow + pptice + SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + pptice + ELSE + SNOWNCV(i,j) = pptsnow + SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + ENDIF + ENDIF + ! Use separate ice if present (as in FV3) + IF ( PRESENT(icencv) .AND. PRESENT(icenc) ) THEN + ICENCV(i,j) = pptice + ICENC(i,j) = ICENC(i,j) + 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) + + + +!..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 +!-GT nwfa1d(kts) = nwfa1 + nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt_in + + do k = kts, kte + nc(i,k,j) = nc1d(k) + nwfa(i,k,j) = nwfa1d(k) + nifa(i,k,j) = nifa1d(k) + enddo + endif +#ifdef DEBUG_AEROSOLS + nc1d_out_debug = nc1d_out_debug + sum(nc1d(:)) + nwfa1d_out_debug = nwfa1d_out_debug + sum(nwfa1d(:)) + nifa1d_out_debug = nifa1d_out_debug + sum(nifa1d(:)) +#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) + if (present(tt)) then + tt(i,k,j) = t1d(k) + else + th(i,k,j) = t1d(k)/pii(i,k,j) + end if +#if ( WRF_CHEM == 1 ) + rainprod(i,k,j) = rainprod1d(k) + evapprod(i,k,j) = evapprod1d(k) +#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(*,*) 'WARNING, negative qc ', qc1d(k), & + ' at i,j,k=', i,j,k + 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(*,*) 'WARNING, negative qr ', qr1d(k), & + ' at i,j,k=', i,j,k + 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(*,*) 'WARNING, negative nr ', nr1d(k), & + ' at i,j,k=', i,j,k + 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(*,*) 'WARNING, negative qs ', qs1d(k), & + ' at i,j,k=', i,j,k + 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(*,*) 'WARNING, negative qi ', qi1d(k), & + ' at i,j,k=', i,j,k + 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(*,*) 'WARNING, negative qg ', qg1d(k), & + ' at i,j,k=', i,j,k + 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(*,*) 'WARNING, negative ni ', ni1d(k), & + ' at i,j,k=', i,j,k + endif + if (qv1d(k) .lt. 0.0) then + write(*,*) 'WARNING, negative qv ', qv1d(k), & + ' at i,j,k=', i,j,k + if (k.lt.kte-2 .and. k.gt.kts+1) then + write(*,*) ' below and above are: ', qv(i,k-1,j), qv(i,k+1,j) + 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 + ! dcd added vt_dBZ, first_time_step + call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & + t1d, p1d, dBZ, vt_dBZ, kts, kte, i, j, first_time_step) + do k = kts, kte + refl_10cm(i,k,j) = MAX(-35., dBZ(k)) + vt_dbz_wt(i,k,j) = vt_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) = 2.49E-6 + re_qi1d(k) = 4.99E-6 + re_qs1d(k) = 9.99E-6 + enddo + call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & + re_qc1d, re_qi1d, re_qs1d, kts, kte, is_aerosol_aware) + do k = kts, kte + re_cloud(i,k,j) = MAX(2.49E-6, MIN(re_qc1d(k), 50.E-6)) + re_ice(i,k,j) = MAX(4.99E-6, MIN(re_qi1d(k), 125.E-6)) + re_snow(i,k,j) = MAX(9.99E-6, MIN(re_qs1d(k), 999.E-6)) + enddo + ENDIF + + enddo i_loop + enddo j_loop +!tgs +! print *,'Thomp MP min/max TH(:,:,:)',minval(th(1,:,:)),maxval(th(1,:,:)) +! print *,'Thomp MP TH(1,:,9960)',th(1,:,9960) + +! DEBUG - GT +! write(*,'(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, ')' +! END DEBUG - GT +#ifdef DEBUG_AEROSOLS + nc1d_in_debug = nc1d_in_debug / real((j_end-j_start+1)*(i_end-i_start+1)*(kte-kts+1)) + nwfa1d_in_debug = nwfa1d_in_debug / real((j_end-i_start+1)*(i_end-i_start+1)*(kte-kts+1)) + nifa1d_in_debug = nifa1d_in_debug / real((j_end-j_start+1)*(i_end-i_start+1)*(kte-kts+1)) + nwfa1_in_debug = nwfa1_in_debug / real((j_end-j_start+1)*(i_end-i_start+1)) + nc1d_out_debug = nc1d_out_debug / real((j_end-j_start+1)*(i_end-i_start+1)*(kte-kts+1)) + nwfa1d_out_debug = nwfa1d_out_debug / real((j_end-i_start+1)*(i_end-i_start+1)*(kte-kts+1)) + nifa1d_out_debug = nifa1d_out_debug / real((j_end-j_start+1)*(i_end-i_start+1)*(kte-kts+1)) + call MPI_BARRIER(MPI_COMM_WORLD,ierr) + do impi=0,mpisize-1 + if (impi==mpirank) then + write(0,'(a,i3,a,l,e16.7)') "MPI rank ", mpirank, & + & ": is_aero, mean(nwfa1) IN:", is_aerosol_aware, nwfa1_in_debug + write(0,'(a,i3,a,3e16.7)') "MPI rank ", mpirank, & + & ": mean(nc1d_in), mean(nwfa1d_in), mean(nifa1d_in):", & + & nc1d_in_debug, nwfa1d_in_debug, nifa1d_in_debug + write(0,'(a,i3,a,3e16.7)') "MPI rank ", mpirank, & + & ": mean(nc1d_out), mean(nwfa1d_out), mean(nifa1d_out):", & + & nc1d_out_debug, nwfa1d_out_debug, nifa1d_out_debug + end if + call MPI_BARRIER(MPI_COMM_WORLD,ierr) + end do +#endif + + END SUBROUTINE mp_gt_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_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + nr1d, nc1d, nwfa1d, nifa1d, t1d, p1d, w1d, dzq, & + pptrain, pptsnow, pptgraul, pptice, & +#if ( WRF_CHEM == 1 ) + rainprod, evapprod, & +#endif + kts, kte, dt, ii, jj, is_aerosol_aware) + + implicit none + +!..Sub arguments + INTEGER, INTENT(IN):: kts, kte, ii, jj + REAL, DIMENSION(kts:kte), INTENT(INOUT):: & + qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + nr1d, nc1d, nwfa1d, nifa1d, t1d + REAL, DIMENSION(kts:kte), INTENT(IN):: p1d, w1d, dzq + REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice + REAL, INTENT(IN):: dt +#if ( WRF_CHEM == 1 ) + REAL, DIMENSION(kts:kte), INTENT(INOUT):: & + rainprod, evapprod +#endif + LOGICAL, INTENT(IN) :: is_aerosol_aware + +!..Local variables + REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, & + qrten, qsten, qgten, niten, nrten, ncten, nwfaten, nifaten + + DOUBLE PRECISION, DIMENSION(kts:kte):: prw_vcd + + DOUBLE PRECISION, DIMENSION(kts:kte):: pnc_wcd, pnc_wau, pnc_rcw, & + pnc_scw, pnc_gcw + + DOUBLE PRECISION, DIMENSION(kts:kte):: pna_rca, pna_sca, pna_gca, & + pnd_rcd, pnd_scd, pnd_gcd + + DOUBLE PRECISION, DIMENSION(kts:kte):: prr_wau, prr_rcw, prr_rcs, & + prr_rcg, prr_sml, prr_gml, & + prr_rci, prv_rev, & + pnr_wau, 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 + + DOUBLE PRECISION, PARAMETER:: zeroD0 = 0.0d0 + + REAL, DIMENSION(kts:kte):: temp, pres, qv + REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni, nr, nc, nwfa, nifa + 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):: ilamr, ilamg, N0_r, N0_g + REAL, DIMENSION(kts:kte):: mvd_r, mvd_c + REAL, DIMENSION(kts:kte):: smob, smo2, smo1, smo0, & + smoc, smod, smoe, smof + + REAL, DIMENSION(kts:kte):: sed_r, sed_s, sed_g, sed_i, sed_n,sed_c + + 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, ilamc + 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, DIMENSION(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk, & + vtck, vtnck + REAL, DIMENSION(kts:kte):: vts_boost + 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:: xnc, xri, xni, xmi, oxmi, xrc, xrr, xnr + 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 + 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 + INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq + 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 + + LOGICAL:: melti, no_micro + LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg + LOGICAL:: debug_flag + ! DH* INTEGER:: idx_lo ! land and ocean + INTEGER:: nu_c + +!+---+ + + debug_flag = .false. +! if (ii.eq.901 .and. jj.eq.379) debug_flag = .true. + if(debug_flag) then + write(*, *) 'DEBUG INFO, mp_thompson at (i,j) ', ii, ', ', jj + endif + + no_micro = .true. + dtsave = dt + odt = 1./dt + odts = 1./dtsave + iexfrq = 1 + + !DH* used in mp_thompson_gfs implementation in FV3 + !if(islmski == 1) then + ! idx_lo = 1 + !else + ! idx_lo = 2 + !DH* endif + +!+---+-----------------------------------------------------------------+ +!.. 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. + niten(k) = 0. + nrten(k) = 0. + ncten(k) = 0. + nwfaten(k) = 0. + nifaten(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. + + pna_rca(k) = 0. + pna_sca(k) = 0. + pna_gca(k) = 0. + + pnd_rcd(k) = 0. + pnd_scd(k) = 0. + pnd_gcd(k) = 0. + enddo +#if ( WRF_CHEM == 1 ) + do k = kts, kte + rainprod(k) = 0. + evapprod(k) = 0. + enddo +#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. + 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(rho(k) < 0.) then +! print *,'rho < 0. at k=',k,rho(k),temp(k),qv(k) +! stop +! endif + 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 (qc1d(k) .gt. R1) then + no_micro = .false. + rc(k) = qc1d(k)*rho(k) + nc(k) = MAX(2., nc1d(k)*rho(k)) + L_qc(k) = .true. + nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) + lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr + xDc = (bm_r + nu_c + 1.) / lamc + if (xDc.lt. D0c) then + lamc = cce(2,nu_c)/D0c + elseif (xDc.gt. D0r*2.) then + lamc = cce(2,nu_c)/(D0r*2.) + endif + nc(k) = MIN( DBLE(Nt_c_max), ccg(1,nu_c)*ocg2(nu_c)*rc(k) & + / am_r*lamc**bm_r) + 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 + + 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)/25.E-6 + ni(k) = MIN(499.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(499.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. + rg(k) = qg1d(k)*rho(k) + L_qg(k) = .true. + else + qg1d(k) = 0.0 + rg(k) = R1 + L_qg(k) = .false. + endif + enddo + +!+---+-----------------------------------------------------------------+ +! if (debug_flag) then +! do k = kts, kte +! write(*, '(a,i3,f8.2,1x,f7.2,1x, 11(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), nwfa(k), nifa(k) +! 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. +!+---+-----------------------------------------------------------------+ + 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) + 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 + enddo + +!+---+-----------------------------------------------------------------+ +!..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 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_ + + enddo + +!+---+-----------------------------------------------------------------+ +!..Calculate y-intercept, slope values for graupel. +!+---+-----------------------------------------------------------------+ + N0_min = gonv_max + k_0 = kts + do k = kte, kts, -1 + if (temp(k).ge.270.65) k_0 = MAX(k_0, k) + enddo + do k = kte, kts, -1 + if (k.gt.k_0 .and. L_qr(k) .and. mvd_r(k).gt.100.E-6) then + xslw1 = 4.01 + alog10(mvd_r(k)) + else + xslw1 = 0.01 + endif + ygra1 = 4.31 + alog10(max(5.E-5, rg(k))) + zans1 = 3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1)) + N0_exp = 10.**(zans1) + N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) + N0_min = MIN(N0_exp, N0_min) + N0_exp = N0_min + lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 + lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg + ilamg(k) = 1./lamg + N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) + 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. RAIN2M + if (L_qr(k) .and. mvd_r(k).gt. D0r) then +!-GT Ef_rr = 1.0 +!-GT if (mvd_r(k) .gt. 1500.0E-6) then + Ef_rr = 1.0 - EXP(2300.0*(mvd_r(k)-1950.0E-6)) +!-GT endif + pnr_rcr(k) = Ef_rr * 2.0*nr(k)*rr(k) + endif + + mvd_c(k) = D0c + if (L_qc(k)) then + nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) + xDc = MAX(D0c*1.E6, ((rc(k)/(am_r*nc(k)))**obmr) * 1.E6) +! print *,'nc(k),am_r,nu_c,rc(k),obmr,k ',nc(k),am_r,nu_c,rc(k),obmr,k + lamc = (nc(k)*am_r* ccg(2,nu_c) * ocg1(nu_c) / rc(k))**obmr + mvd_c(k) = (3.0+nu_c+0.672) / lamc +!DH* used in mp_thompson_gfs implementation in FV3 +! if(islmski == 1) then +! mvd_c(k) = (3.0+mu_cl+0.672) / lamc +! else +! mvd_c(k) = (3.0+mu_co+0.672) / lamc +! endif + endif + +!..Autoconversion follows Berry & Reinhardt (1974) with characteristic +!.. diameters correctly computed from gamma distrib of cloud droplets. + if (rc(k).gt. 0.01e-3) then + Dc_g = ((ccg(3,nu_c)*ocg2(nu_c))**obmr / lamc) * 1.E6 +!DH* used in mp_thompson_gfs implementation in FV3 +! if(islmski == 1) then +! Dc_g = ((ccg(3,1)*ocg2(1))**obmr / lamc) * 1.E6 +! else +! Dc_g = ((ccg(3,2)*ocg2(2))**obmr / lamc) * 1.E6 +! endif + Dc_b = (xDc*xDc*xDc*Dc_g*Dc_g*Dc_g - xDc*xDc*xDc*xDc*xDc*xDc) & + **(1./6.) + zeta1 = 0.5*((6.25E-6*xDc*Dc_b*Dc_b*Dc_b - 0.4) & + + abs(6.25E-6*xDc*Dc_b*Dc_b*Dc_b - 0.4)) + zeta = 0.027*rc(k)*zeta1 + taud = 0.5*((0.5*Dc_b - 7.5) + abs(0.5*Dc_b - 7.5)) + R1 + tau = 3.72/(rc(k)*taud) + prr_wau(k) = zeta/tau + prr_wau(k) = MIN(DBLE(rc(k)*odts), prr_wau(k)) + pnr_wau(k) = prr_wau(k) / (am_r*nu_c*D0r*D0r*D0r) ! RAIN2M +!DH* used in mp_thompson_gfs implementation in FV3 +! if(islmski == 1) then +! pnr_wau(k) = prr_wau(k) / (am_r*mu_cl*D0r*D0r*D0r) ! RAIN2M +! else +! pnr_wau(k) = prr_wau(k) / (am_r*mu_co*D0r*D0r*D0r) ! RAIN2M +! endif + pnc_wau(k) = MIN(DBLE(nc(k)*odts), prr_wau(k) & + / (am_r*mvd_c(k)*mvd_c(k)*mvd_c(k))) ! Qc2M + endif + +!..Rain collecting cloud water. In CE, assume Dc<kte-3.and.jj==31) then +! print *,'nu_c=',ii,jj,k,nu_c +! endif + lamc = (am_r*ccg(2,nu_c)*ocg1(nu_c)*nc1d(k)/qc1d(k))**obmr + xDc = (bm_r + nu_c + 1.) / lamc + if (xDc.lt. D0c) then + lamc = cce(2,nu_c)/D0c + elseif (xDc.gt. D0r*2.) then + lamc = cce(2,nu_c)/(D0r*2.) + endif + nc1d(k) = MIN(ccg(1,nu_c)*ocg2(nu_c)*qc1d(k)/am_r*lamc**bm_r,& + DBLE(Nt_c_max)/rho(k)) + 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, & + 499.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 + if (qg1d(k) .le. R1) qg1d(k) = 0.0 + enddo + + end subroutine mp_thompson +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Creation of the lookup tables and support functions found below here. +!+---+-----------------------------------------------------------------+ +!..Rain collecting graupel (and inverse). Explicit CE integration. +!+---+-----------------------------------------------------------------+ + + subroutine qr_acr_qg + + implicit none + +!..Local variables + INTEGER:: i, j, k, m, n, n2 + INTEGER:: km, km_s, km_e + DOUBLE PRECISION, DIMENSION(nbg):: vg, N_g + DOUBLE PRECISION, DIMENSION(nbr):: vr, N_r + DOUBLE PRECISION:: N0_r, N0_g, lam_exp, lamg, lamr + DOUBLE PRECISION:: massg, massr, dvg, dvr, t1, t2, z1, z2, y1, y2 + LOGICAL force_read_thompson, write_thompson_tables + LOGICAL lexist,lopen + INTEGER good + + force_read_thompson = .false. + write_thompson_tables = .false. +!+---+ + +! CALL nl_get_force_read_thompson(1,force_read_thompson) +! CALL nl_get_write_thompson_tables(1,write_thompson_tables) + + good = 0 +! IF ( wrf_dm_on_monitor() ) THEN + INQUIRE(FILE="qr_acr_qg.dat",EXIST=lexist) + IF ( lexist ) THEN + write(0,*) "ThompMP: read qr_acr_qg.dat stead of computing" + OPEN(63,file="qr_acr_qg.dat",form="unformatted",err=1234) +!sms$serial begin + READ(63,err=1234) tcg_racg + READ(63,err=1234) tmr_racg + READ(63,err=1234) tcr_gacr + READ(63,err=1234) tmg_gacr + READ(63,err=1234) tnr_racg + READ(63,err=1234) tnr_gacr +!sms$serial end + + good = 1 + 1234 CONTINUE + IF ( good .NE. 1 ) THEN + INQUIRE(63,opened=lopen) + IF (lopen) THEN + IF( force_read_thompson ) THEN + write(0,*) "Error reading qr_acr_qg.dat. Aborting because force_read_thompson is .true." ! DH* errmsg + return + ENDIF + CLOSE(63) + ELSE + IF( force_read_thompson ) THEN + write(0,*) "Error opening qr_acr_qg.dat. Aborting because force_read_thompson is .true." ! DH* errmsg + return + ENDIF + ENDIF + ELSE + INQUIRE(63,opened=lopen) + IF (lopen) THEN + CLOSE(63) + ENDIF + ENDIF + ELSE + IF( force_read_thompson ) THEN + write(0,*) "Non-existent qr_acr_qg.dat. Aborting because force_read_thompson is .true." ! DH* errmsg + return + ENDIF + ENDIF +! ENDIF +!#if defined(DM_PARALLEL) && !defined(STUBMPI) +! CALL wrf_dm_bcast_integer(good,1) +!#endif + + IF ( good .EQ. 1 ) THEN +!#if defined(DM_PARALLEL) && !defined(STUBMPI) +! CALL wrf_dm_bcast_double(tcg_racg,SIZE(tcg_racg)) +! CALL wrf_dm_bcast_double(tmr_racg,SIZE(tmr_racg)) +! CALL wrf_dm_bcast_double(tcr_gacr,SIZE(tcr_gacr)) +! CALL wrf_dm_bcast_double(tmg_gacr,SIZE(tmg_gacr)) +! CALL wrf_dm_bcast_double(tnr_racg,SIZE(tnr_racg)) +! CALL wrf_dm_bcast_double(tnr_gacr,SIZE(tnr_gacr)) +!#endif + ELSE + write(0,*) "ThompMP: computing qr_acr_qg" + do n2 = 1, nbr +! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2)) + vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) & + + 0.07934E9*Dr(n2)*Dr(n2)*Dr(n2) & + - 0.002362E12*Dr(n2)*Dr(n2)*Dr(n2)*Dr(n2) + enddo + do n = 1, nbg + vg(n) = av_g*Dg(n)**bv_g + enddo + +!..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for +!.. fortran indices. J. Michalakes, 2009Oct30. + +#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e ) +#else + km_s = 0 + km_e = ntb_r*ntb_r1 - 1 +#endif + + do km = km_s, km_e + m = km / ntb_r1 + 1 + k = mod( km , ntb_r1 ) + 1 + + lam_exp = (N0r_exp(k)*am_r*crg(1)/r_r(m))**ore1 + lamr = lam_exp * (crg(3)*org2*org1)**obmr + N0_r = N0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2) + do n2 = 1, nbr + N_r(n2) = N0_r*Dr(n2)**mu_r *DEXP(-lamr*Dr(n2))*dtr(n2) + enddo + + do j = 1, ntb_g + do i = 1, ntb_g1 + lam_exp = (N0g_exp(i)*am_g*cgg(1)/r_g(j))**oge1 + lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg + N0_g = N0g_exp(i)/(cgg(2)*lam_exp) * lamg**cge(2) + do n = 1, nbg + N_g(n) = N0_g*Dg(n)**mu_g * DEXP(-lamg*Dg(n))*dtg(n) + enddo + + t1 = 0.0d0 + t2 = 0.0d0 + z1 = 0.0d0 + z2 = 0.0d0 + y1 = 0.0d0 + y2 = 0.0d0 + do n2 = 1, nbr + massr = am_r * Dr(n2)**bm_r + do n = 1, nbg + massg = am_g * Dg(n)**bm_g + + dvg = 0.5d0*((vr(n2) - vg(n)) + DABS(vr(n2)-vg(n))) + dvr = 0.5d0*((vg(n) - vr(n2)) + DABS(vg(n)-vr(n2))) + + t1 = t1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & + *dvg*massg * N_g(n)* N_r(n2) + z1 = z1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & + *dvg*massr * N_g(n)* N_r(n2) + y1 = y1+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & + *dvg * N_g(n)* N_r(n2) + + t2 = t2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & + *dvr*massr * N_g(n)* N_r(n2) + y2 = y2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & + *dvr * N_g(n)* N_r(n2) + z2 = z2+ PI*.25*Ef_rg*(Dg(n)+Dr(n2))*(Dg(n)+Dr(n2)) & + *dvr*massg * N_g(n)* N_r(n2) + enddo + 97 continue + enddo + tcg_racg(i,j,k,m) = t1 + tmr_racg(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0) + tcr_gacr(i,j,k,m) = t2 + tmg_gacr(i,j,k,m) = z2 + tnr_racg(i,j,k,m) = y1 + tnr_gacr(i,j,k,m) = y2 + enddo + enddo + enddo + +!..Note wrf_dm_gatherv expects zero-based km_s, km_e (J. Michalakes, 2009Oct30). + +#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + CALL wrf_dm_gatherv(tcg_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tmr_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tcr_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tmg_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tnr_racg, ntb_g*ntb_g1, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tnr_gacr, ntb_g*ntb_g1, km_s, km_e, R8SIZE) +#endif + + IF ( write_thompson_tables ) THEN + write(0,*) "Writing qr_acr_qg.dat in Thompson MP init" + OPEN(63,file="qr_acr_qg.dat",form="unformatted",err=9234) + WRITE(63,err=9234) tcg_racg + WRITE(63,err=9234) tmr_racg + WRITE(63,err=9234) tcr_gacr + WRITE(63,err=9234) tmg_gacr + WRITE(63,err=9234) tnr_racg + WRITE(63,err=9234) tnr_gacr + CLOSE(63) + RETURN ! ----- RETURN + 9234 CONTINUE + write(0,*) "Error writing qr_acr_qg.dat" ! DH* errmsg + return + ENDIF + ENDIF + + end subroutine qr_acr_qg +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Rain collecting snow (and inverse). Explicit CE integration. +!+---+-----------------------------------------------------------------+ + + subroutine qr_acr_qs + + implicit none + +!..Local variables + INTEGER:: i, j, k, m, n, n2 + INTEGER:: km, km_s, km_e + DOUBLE PRECISION, DIMENSION(nbr):: vr, D1, N_r + DOUBLE PRECISION, DIMENSION(nbs):: vs, N_s + DOUBLE PRECISION:: loga_, a_, b_, second, M0, M2, M3, Mrat, oM3 + DOUBLE PRECISION:: N0_r, lam_exp, lamr, slam1, slam2 + DOUBLE PRECISION:: dvs, dvr, masss, massr + DOUBLE PRECISION:: t1, t2, t3, t4, z1, z2, z3, z4 + DOUBLE PRECISION:: y1, y2, y3, y4 + LOGICAL force_read_thompson, write_thompson_tables + LOGICAL lexist,lopen + INTEGER good +! LOGICAL, EXTERNAL :: wrf_dm_on_monitor + +!+---+ + + force_read_thompson = .false. + write_thompson_tables = .false. + +! CALL nl_get_force_read_thompson(1,force_read_thompson) +! CALL nl_get_write_thompson_tables(1,write_thompson_tables) + + good = 0 +! IF ( wrf_dm_on_monitor() ) THEN + INQUIRE(FILE="qr_acr_qs.dat",EXIST=lexist) + IF ( lexist ) THEN + write(0,*) "ThompMP: read qr_acr_qs.dat instead of computing" + OPEN(63,file="qr_acr_qs.dat",form="unformatted",err=1234) +!sms$serial begin + READ(63,err=1234)tcs_racs1 + READ(63,err=1234)tmr_racs1 + READ(63,err=1234)tcs_racs2 + READ(63,err=1234)tmr_racs2 + READ(63,err=1234)tcr_sacr1 + READ(63,err=1234)tms_sacr1 + READ(63,err=1234)tcr_sacr2 + READ(63,err=1234)tms_sacr2 + READ(63,err=1234)tnr_racs1 + READ(63,err=1234)tnr_racs2 + READ(63,err=1234)tnr_sacr1 + READ(63,err=1234)tnr_sacr2 +!sms$serial end + + good = 1 + 1234 CONTINUE + IF ( good .NE. 1 ) THEN + INQUIRE(63,opened=lopen) + IF (lopen) THEN + IF( force_read_thompson ) THEN + write(0,*) "Error reading qr_acr_qs.dat. Aborting because force_read_thompson is .true." ! DH* errmsg + return + ENDIF + CLOSE(63) + ELSE + IF( force_read_thompson ) THEN + write(0,*) "Error opening qr_acr_qs.dat. Aborting because force_read_thompson is .true." ! DH* errmsg + return + ENDIF + ENDIF + ENDIF + ELSE + IF( force_read_thompson ) THEN + write(0,*) "Non-existent qr_acr_qs.dat. Aborting because force_read_thompson is .true." ! DH* errmsg + return + ENDIF + ENDIF +! ENDIF +!#if defined(DM_PARALLEL) && !defined(STUBMPI) +! CALL wrf_dm_bcast_integer(good,1) +!#endif + + IF ( good .EQ. 1 ) THEN +!#if defined(DM_PARALLEL) && !defined(STUBMPI) +! CALL wrf_dm_bcast_double(tcs_racs1,SIZE(tcs_racs1)) +! CALL wrf_dm_bcast_double(tmr_racs1,SIZE(tmr_racs1)) +! CALL wrf_dm_bcast_double(tcs_racs2,SIZE(tcs_racs2)) +! CALL wrf_dm_bcast_double(tmr_racs2,SIZE(tmr_racs2)) +! CALL wrf_dm_bcast_double(tcr_sacr1,SIZE(tcr_sacr1)) +! CALL wrf_dm_bcast_double(tms_sacr1,SIZE(tms_sacr1)) +! CALL wrf_dm_bcast_double(tcr_sacr2,SIZE(tcr_sacr2)) +! CALL wrf_dm_bcast_double(tms_sacr2,SIZE(tms_sacr2)) +! CALL wrf_dm_bcast_double(tnr_racs1,SIZE(tnr_racs1)) +! CALL wrf_dm_bcast_double(tnr_racs2,SIZE(tnr_racs2)) +! CALL wrf_dm_bcast_double(tnr_sacr1,SIZE(tnr_sacr1)) +! CALL wrf_dm_bcast_double(tnr_sacr2,SIZE(tnr_sacr2)) +!#endif + ELSE + write(0,*) "ThompMP: computing qr_acr_qs" + do n2 = 1, nbr +! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2)) + vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) & + + 0.07934E9*Dr(n2)*Dr(n2)*Dr(n2) & + - 0.002362E12*Dr(n2)*Dr(n2)*Dr(n2)*Dr(n2) + D1(n2) = (vr(n2)/av_s)**(1./bv_s) + enddo + do n = 1, nbs + vs(n) = 1.5*av_s*Ds(n)**bv_s * DEXP(-fv_s*Ds(n)) + enddo + +!..Note values returned from wrf_dm_decomp1d are zero-based, add 1 for +!.. fortran indices. J. Michalakes, 2009Oct30. + +#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + CALL wrf_dm_decomp1d ( ntb_r*ntb_r1, km_s, km_e ) +#else + km_s = 0 + km_e = ntb_r*ntb_r1 - 1 +#endif + + do km = km_s, km_e + m = km / ntb_r1 + 1 + k = mod( km , ntb_r1 ) + 1 + + lam_exp = (N0r_exp(k)*am_r*crg(1)/r_r(m))**ore1 + lamr = lam_exp * (crg(3)*org2*org1)**obmr + N0_r = N0r_exp(k)/(crg(2)*lam_exp) * lamr**cre(2) + do n2 = 1, nbr + N_r(n2) = N0_r*Dr(n2)**mu_r * DEXP(-lamr*Dr(n2))*dtr(n2) + enddo + + do j = 1, ntb_t + do i = 1, ntb_s + +!..From the bm_s moment, compute plus one moment. If we are not +!.. using bm_s=2, then we must transform to the pure 2nd moment +!.. (variable called "second") and then to the bm_s+1 moment. + + M2 = r_s(i)*oams *1.0d0 + if (bm_s.gt.2.0-1.E-3 .and. bm_s.lt.2.0+1.E-3) then + loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*bm_s & + + sa(4)*Tc(j)*bm_s + sa(5)*Tc(j)*Tc(j) & + + sa(6)*bm_s*bm_s + sa(7)*Tc(j)*Tc(j)*bm_s & + + sa(8)*Tc(j)*bm_s*bm_s + sa(9)*Tc(j)*Tc(j)*Tc(j) & + + sa(10)*bm_s*bm_s*bm_s + a_ = 10.0**loga_ + b_ = sb(1) + sb(2)*Tc(j) + sb(3)*bm_s & + + sb(4)*Tc(j)*bm_s + sb(5)*Tc(j)*Tc(j) & + + sb(6)*bm_s*bm_s + sb(7)*Tc(j)*Tc(j)*bm_s & + + sb(8)*Tc(j)*bm_s*bm_s + sb(9)*Tc(j)*Tc(j)*Tc(j) & + + sb(10)*bm_s*bm_s*bm_s + second = (M2/a_)**(1./b_) + else + second = M2 + endif + + loga_ = sa(1) + sa(2)*Tc(j) + sa(3)*cse(1) & + + sa(4)*Tc(j)*cse(1) + sa(5)*Tc(j)*Tc(j) & + + sa(6)*cse(1)*cse(1) + sa(7)*Tc(j)*Tc(j)*cse(1) & + + sa(8)*Tc(j)*cse(1)*cse(1) + sa(9)*Tc(j)*Tc(j)*Tc(j) & + + sa(10)*cse(1)*cse(1)*cse(1) + a_ = 10.0**loga_ + b_ = sb(1)+sb(2)*Tc(j)+sb(3)*cse(1) + sb(4)*Tc(j)*cse(1) & + + sb(5)*Tc(j)*Tc(j) + sb(6)*cse(1)*cse(1) & + + sb(7)*Tc(j)*Tc(j)*cse(1) + sb(8)*Tc(j)*cse(1)*cse(1) & + + sb(9)*Tc(j)*Tc(j)*Tc(j)+sb(10)*cse(1)*cse(1)*cse(1) + M3 = a_ * second**b_ + + oM3 = 1./M3 + Mrat = M2*(M2*oM3)*(M2*oM3)*(M2*oM3) + M0 = (M2*oM3)**mu_s + slam1 = M2 * oM3 * Lam0 + slam2 = M2 * oM3 * Lam1 + + do n = 1, nbs + N_s(n) = Mrat*(Kap0*DEXP(-slam1*Ds(n)) & + + Kap1*M0*Ds(n)**mu_s * DEXP(-slam2*Ds(n)))*dts(n) + enddo + + t1 = 0.0d0 + t2 = 0.0d0 + t3 = 0.0d0 + t4 = 0.0d0 + z1 = 0.0d0 + z2 = 0.0d0 + z3 = 0.0d0 + z4 = 0.0d0 + y1 = 0.0d0 + y2 = 0.0d0 + y3 = 0.0d0 + y4 = 0.0d0 + do n2 = 1, nbr + massr = am_r * Dr(n2)**bm_r + do n = 1, nbs + masss = am_s * Ds(n)**bm_s + + dvs = 0.5d0*((vr(n2) - vs(n)) + DABS(vr(n2)-vs(n))) + dvr = 0.5d0*((vs(n) - vr(n2)) + DABS(vs(n)-vr(n2))) + + if (massr .gt. 1.5*masss) then + t1 = t1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvs*masss * N_s(n)* N_r(n2) + z1 = z1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvs*massr * N_s(n)* N_r(n2) + y1 = y1+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvs * N_s(n)* N_r(n2) + else + t3 = t3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvs*masss * N_s(n)* N_r(n2) + z3 = z3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvs*massr * N_s(n)* N_r(n2) + y3 = y3+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvs * N_s(n)* N_r(n2) + endif + + if (massr .gt. 1.5*masss) then + t2 = t2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvr*massr * N_s(n)* N_r(n2) + y2 = y2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvr * N_s(n)* N_r(n2) + z2 = z2+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvr*masss * N_s(n)* N_r(n2) + else + t4 = t4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvr*massr * N_s(n)* N_r(n2) + y4 = y4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvr * N_s(n)* N_r(n2) + z4 = z4+ PI*.25*Ef_rs*(Ds(n)+Dr(n2))*(Ds(n)+Dr(n2)) & + *dvr*masss * N_s(n)* N_r(n2) + endif + + enddo + enddo + tcs_racs1(i,j,k,m) = t1 + tmr_racs1(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0) + tcs_racs2(i,j,k,m) = t3 + tmr_racs2(i,j,k,m) = z3 + tcr_sacr1(i,j,k,m) = t2 + tms_sacr1(i,j,k,m) = z2 + tcr_sacr2(i,j,k,m) = t4 + tms_sacr2(i,j,k,m) = z4 + tnr_racs1(i,j,k,m) = y1 + tnr_racs2(i,j,k,m) = y3 + tnr_sacr1(i,j,k,m) = y2 + tnr_sacr2(i,j,k,m) = y4 + enddo + enddo + enddo + +!..Note wrf_dm_gatherv expects zero-based km_s, km_e (J. Michalakes, 2009Oct30). + +#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) ) + CALL wrf_dm_gatherv(tcs_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tmr_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tcs_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tmr_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tcr_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tms_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tcr_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tms_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tnr_racs1, ntb_s*ntb_t, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tnr_racs2, ntb_s*ntb_t, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tnr_sacr1, ntb_s*ntb_t, km_s, km_e, R8SIZE) + CALL wrf_dm_gatherv(tnr_sacr2, ntb_s*ntb_t, km_s, km_e, R8SIZE) +#endif + + IF ( write_thompson_tables ) THEN + write(0,*) "Writing qr_acr_qs.dat in Thompson MP init" + OPEN(63,file="qr_acr_qs.dat",form="unformatted",err=9234) + WRITE(63,err=9234)tcs_racs1 + WRITE(63,err=9234)tmr_racs1 + WRITE(63,err=9234)tcs_racs2 + WRITE(63,err=9234)tmr_racs2 + WRITE(63,err=9234)tcr_sacr1 + WRITE(63,err=9234)tms_sacr1 + WRITE(63,err=9234)tcr_sacr2 + WRITE(63,err=9234)tms_sacr2 + WRITE(63,err=9234)tnr_racs1 + WRITE(63,err=9234)tnr_racs2 + WRITE(63,err=9234)tnr_sacr1 + WRITE(63,err=9234)tnr_sacr2 + CLOSE(63) + RETURN ! ----- RETURN + 9234 CONTINUE + write(0,*) "Error writing qr_acr_qs.dat" ! DH* errmsg + ENDIF + ENDIF + + end subroutine qr_acr_qs +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!..This is a literal adaptation of Bigg (1954) probability of drops of +!..a particular volume freezing. Given this probability, simply freeze +!..the proportion of drops summing their masses. +!+---+-----------------------------------------------------------------+ + + subroutine freezeH2O + + implicit none + +!..Local variables + INTEGER:: i, j, k, m, n, n2 + DOUBLE PRECISION, DIMENSION(nbr):: N_r, massr + DOUBLE PRECISION, DIMENSION(nbc):: N_c, massc + DOUBLE PRECISION:: sum1, sum2, sumn1, sumn2, & + prob, vol, Texp, orho_w, & + lam_exp, lamr, N0_r, lamc, N0_c, y + INTEGER:: nu_c + REAL:: T_adjust + LOGICAL force_read_thompson, write_thompson_tables + LOGICAL lexist,lopen + INTEGER good +! LOGICAL, EXTERNAL :: wrf_dm_on_monitor + +!+---+ +! CALL nl_get_force_read_thompson(1,force_read_thompson) +! CALL nl_get_write_thompson_tables(1,write_thompson_tables) + force_read_thompson = .false. + write_thompson_tables = .false. + + + good = 0 +! IF ( wrf_dm_on_monitor() ) THEN + INQUIRE(FILE="freezeH2O.dat",EXIST=lexist) + IF ( lexist ) THEN + write(0,*) "ThompMP: read freezeH2O.dat stead of computing" + OPEN(63,file="freezeH2O.dat",form="unformatted",err=1234) +!sms$serial begin + READ(63,err=1234)tpi_qrfz + READ(63,err=1234)tni_qrfz + READ(63,err=1234)tpg_qrfz + READ(63,err=1234)tnr_qrfz + READ(63,err=1234)tpi_qcfz + READ(63,err=1234)tni_qcfz +!sms$serial end + + good = 1 + 1234 CONTINUE + IF ( good .NE. 1 ) THEN + INQUIRE(63,opened=lopen) + IF (lopen) THEN + IF( force_read_thompson ) THEN + write(0,*) "Error reading freezeH2O.dat. Aborting because force_read_thompson is .true." ! DH* errmsg + return + ENDIF + CLOSE(63) + ELSE + IF( force_read_thompson ) THEN + write(0,*) "Error opening freezeH2O.dat. Aborting because force_read_thompson is .true." ! DH* errmsg + return + ENDIF + ENDIF + ENDIF + ELSE + IF( force_read_thompson ) THEN + write(0,*) "Non-existent freezeH2O.dat. Aborting because force_read_thompson is .true." ! DH* errmsg + return + ENDIF + ENDIF +! ENDIF + +!#if defined(DM_PARALLEL) && !defined(STUBMPI) +! CALL wrf_dm_bcast_integer(good,1) +!#endif + + IF ( good .EQ. 1 ) THEN +!#if defined(DM_PARALLEL) && !defined(STUBMPI) +! CALL wrf_dm_bcast_double(tpi_qrfz,SIZE(tpi_qrfz)) +! CALL wrf_dm_bcast_double(tni_qrfz,SIZE(tni_qrfz)) +! CALL wrf_dm_bcast_double(tpg_qrfz,SIZE(tpg_qrfz)) +! CALL wrf_dm_bcast_double(tnr_qrfz,SIZE(tnr_qrfz)) +! CALL wrf_dm_bcast_double(tpi_qcfz,SIZE(tpi_qcfz)) +! CALL wrf_dm_bcast_double(tni_qcfz,SIZE(tni_qcfz)) +!#endif + ELSE + write(0,*) "ThompMP: computing freezeH2O" + + orho_w = 1./rho_w + + do n2 = 1, nbr + massr(n2) = am_r*Dr(n2)**bm_r + enddo + do n = 1, nbc + massc(n) = am_r*Dc(n)**bm_r + enddo + +!..Freeze water (smallest drops become cloud ice, otherwise graupel). + do m = 1, ntb_IN + T_adjust = MAX(-3.0, MIN(3.0 - ALOG10(Nt_IN(m)), 3.0)) + do k = 1, 45 +! print*, ' Freezing water for temp = ', -k + Texp = DEXP( DFLOAT(k) - T_adjust*1.0D0 ) - 1.0D0 + do j = 1, ntb_r1 + do i = 1, ntb_r + lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(i))**ore1 + lamr = lam_exp * (crg(3)*org2*org1)**obmr + N0_r = N0r_exp(j)/(crg(2)*lam_exp) * lamr**cre(2) + sum1 = 0.0d0 + sum2 = 0.0d0 + sumn1 = 0.0d0 + sumn2 = 0.0d0 + do n2 = nbr, 1, -1 + N_r(n2) = N0_r*Dr(n2)**mu_r*DEXP(-lamr*Dr(n2))*dtr(n2) + vol = massr(n2)*orho_w + prob = 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp) + if (massr(n2) .lt. xm0g) then + sumn1 = sumn1 + prob*N_r(n2) + sum1 = sum1 + prob*N_r(n2)*massr(n2) + else + sumn2 = sumn2 + prob*N_r(n2) + sum2 = sum2 + prob*N_r(n2)*massr(n2) + endif + if ((sum1+sum2).ge.r_r(i)) EXIT + enddo + tpi_qrfz(i,j,k,m) = sum1 + tni_qrfz(i,j,k,m) = sumn1 + tpg_qrfz(i,j,k,m) = sum2 + tnr_qrfz(i,j,k,m) = sumn2 + enddo + enddo + + do j = 1, nbc + nu_c = MIN(15, NINT(1000.E6/t_Nc(j)) + 2) + do i = 1, ntb_c + lamc = (t_Nc(j)*am_r* ccg(2,nu_c) * ocg1(nu_c) / r_c(i))**obmr + N0_c = t_Nc(j)*ocg1(nu_c) * lamc**cce(1,nu_c) + sum1 = 0.0d0 + sumn2 = 0.0d0 + do n = nbc, 1, -1 + vol = massc(n)*orho_w + prob = 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp) + N_c(n) = N0_c*Dc(n)**nu_c*EXP(-lamc*Dc(n))*dtc(n) + sumn2 = MIN(t_Nc(j), sumn2 + prob*N_c(n)) + sum1 = sum1 + prob*N_c(n)*massc(n) + if (sum1 .ge. r_c(i)) EXIT + enddo + tpi_qcfz(i,j,k,m) = sum1 + tni_qcfz(i,j,k,m) = sumn2 + enddo + enddo + enddo + enddo + + IF ( write_thompson_tables ) THEN + write(0,*) "Writing freezeH2O.dat in Thompson MP init" + OPEN(63,file="freezeH2O.dat",form="unformatted",err=9234) + WRITE(63,err=9234)tpi_qrfz + WRITE(63,err=9234)tni_qrfz + WRITE(63,err=9234)tpg_qrfz + WRITE(63,err=9234)tnr_qrfz + WRITE(63,err=9234)tpi_qcfz + WRITE(63,err=9234)tni_qcfz + CLOSE(63) + RETURN ! ----- RETURN + 9234 CONTINUE + write(0,*) "Error writing freezeH2O.dat" ! DH* errmsg + return + ENDIF + ENDIF + + end subroutine freezeH2O +!+---+-----------------------------------------------------------------+ +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Cloud ice converting to snow since portion greater than min snow +!.. size. Given cloud ice content (kg/m**3), number concentration +!.. (#/m**3) and gamma shape parameter, mu_i, break the distrib into +!.. bins and figure out the mass/number of ice with sizes larger than +!.. D0s. Also, compute incomplete gamma function for the integration +!.. of ice depositional growth from diameter=0 to D0s. Amount of +!.. ice depositional growth is this portion of distrib while larger +!.. diameters contribute to snow growth (as in Harrington et al. 1995). +!+---+-----------------------------------------------------------------+ + + subroutine qi_aut_qs + + implicit none + +!..Local variables + INTEGER:: i, j, n2 + DOUBLE PRECISION, DIMENSION(nbi):: N_i + DOUBLE PRECISION:: N0_i, lami, Di_mean, t1, t2 + REAL:: xlimit_intg + +!+---+ + + do j = 1, ntb_i1 + do i = 1, ntb_i + lami = (am_i*cig(2)*oig1*Nt_i(j)/r_i(i))**obmi + Di_mean = (bm_i + mu_i + 1.) / lami + N0_i = Nt_i(j)*oig1 * lami**cie(1) + t1 = 0.0d0 + t2 = 0.0d0 + if (SNGL(Di_mean) .gt. 5.*D0s) then + t1 = r_i(i) + t2 = Nt_i(j) + tpi_ide(i,j) = 0.0D0 + elseif (SNGL(Di_mean) .lt. D0i) then + t1 = 0.0D0 + t2 = 0.0D0 + tpi_ide(i,j) = 1.0D0 + else + xlimit_intg = lami*D0s + tpi_ide(i,j) = GAMMP(mu_i+2.0, xlimit_intg) * 1.0D0 + do n2 = 1, nbi + N_i(n2) = N0_i*Di(n2)**mu_i * DEXP(-lami*Di(n2))*dti(n2) + if (Di(n2).ge.D0s) then + t1 = t1 + N_i(n2) * am_i*Di(n2)**bm_i + t2 = t2 + N_i(n2) + endif + enddo + endif + tps_iaus(i,j) = t1 + tni_iaus(i,j) = t2 + enddo + enddo + + end subroutine qi_aut_qs +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Variable collision efficiency for rain collecting cloud water using +!.. method of Beard and Grover, 1974 if a/A less than 0.25; otherwise +!.. uses polynomials to get close match of Pruppacher & Klett Fig 14-9. +!+---+-----------------------------------------------------------------+ + + subroutine table_Efrw + + implicit none + +!..Local variables + DOUBLE PRECISION:: vtr, stokes, reynolds, Ef_rw + DOUBLE PRECISION:: p, yc0, F, G, H, z, K0, X + INTEGER:: i, j + + do j = 1, nbc + do i = 1, nbr + Ef_rw = 0.0 + p = Dc(j)/Dr(i) + if (Dr(i).lt.50.E-6 .or. Dc(j).lt.3.E-6) then + t_Efrw(i,j) = 0.0 + elseif (p.gt.0.25) then + X = Dc(j)*1.D6 + if (Dr(i) .lt. 75.e-6) then + Ef_rw = 0.026794*X - 0.20604 + elseif (Dr(i) .lt. 125.e-6) then + Ef_rw = -0.00066842*X*X + 0.061542*X - 0.37089 + elseif (Dr(i) .lt. 175.e-6) then + Ef_rw = 4.091e-06*X*X*X*X - 0.00030908*X*X*X & + + 0.0066237*X*X - 0.0013687*X - 0.073022 + elseif (Dr(i) .lt. 250.e-6) then + Ef_rw = 9.6719e-5*X*X*X - 0.0068901*X*X + 0.17305*X & + - 0.65988 + elseif (Dr(i) .lt. 350.e-6) then + Ef_rw = 9.0488e-5*X*X*X - 0.006585*X*X + 0.16606*X & + - 0.56125 + else + Ef_rw = 0.00010721*X*X*X - 0.0072962*X*X + 0.1704*X & + - 0.46929 + endif + else + vtr = -0.1021 + 4.932E3*Dr(i) - 0.9551E6*Dr(i)*Dr(i) & + + 0.07934E9*Dr(i)*Dr(i)*Dr(i) & + - 0.002362E12*Dr(i)*Dr(i)*Dr(i)*Dr(i) + stokes = Dc(j)*Dc(j)*vtr*rho_w/(9.*1.718E-5*Dr(i)) + reynolds = 9.*stokes/(p*p*rho_w) + + F = DLOG(reynolds) + G = -0.1007D0 - 0.358D0*F + 0.0261D0*F*F + K0 = DEXP(G) + z = DLOG(stokes/(K0+1.D-15)) + H = 0.1465D0 + 1.302D0*z - 0.607D0*z*z + 0.293D0*z*z*z + yc0 = 2.0D0/PI * ATAN(H) + Ef_rw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p)) + + endif + + t_Efrw(i,j) = MAX(0.0, MIN(SNGL(Ef_rw), 0.95)) + + enddo + enddo + + end subroutine table_Efrw +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Variable collision efficiency for snow collecting cloud water using +!.. method of Wang and Ji, 2000 except equate melted snow diameter to +!.. their "effective collision cross-section." +!+---+-----------------------------------------------------------------+ + + subroutine table_Efsw + + implicit none + +!..Local variables + DOUBLE PRECISION:: Ds_m, vts, vtc, stokes, reynolds, Ef_sw + DOUBLE PRECISION:: p, yc0, F, G, H, z, K0 + INTEGER:: i, j + + do j = 1, nbc + vtc = 1.19D4 * (1.0D4*Dc(j)*Dc(j)*0.25D0) + do i = 1, nbs + vts = av_s*Ds(i)**bv_s * DEXP(-fv_s*Ds(i)) - vtc + Ds_m = (am_s*Ds(i)**bm_s / am_r)**obmr + p = Dc(j)/Ds_m + if (p.gt.0.25 .or. Ds(i).lt.D0s .or. Dc(j).lt.6.E-6 & + .or. vts.lt.1.E-3) then + t_Efsw(i,j) = 0.0 + else + stokes = Dc(j)*Dc(j)*vts*rho_w/(9.*1.718E-5*Ds_m) + reynolds = 9.*stokes/(p*p*rho_w) + + F = DLOG(reynolds) + G = -0.1007D0 - 0.358D0*F + 0.0261D0*F*F + K0 = DEXP(G) + z = DLOG(stokes/(K0+1.D-15)) + H = 0.1465D0 + 1.302D0*z - 0.607D0*z*z + 0.293D0*z*z*z + yc0 = 2.0D0/PI * ATAN(H) + Ef_sw = (yc0+p)*(yc0+p) / ((1.+p)*(1.+p)) + + t_Efsw(i,j) = MAX(0.0, MIN(SNGL(Ef_sw), 0.95)) + endif + + enddo + enddo + + end subroutine table_Efsw +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Function to compute collision efficiency of collector species (rain, +!.. snow, graupel) of aerosols. Follows Wang et al, 2010, ACP, which +!.. follows Slinn (1983). +!+---+-----------------------------------------------------------------+ + + real function Eff_aero(D, Da, visc,rhoa,Temp,species) + + implicit none + real:: D, Da, visc, rhoa, Temp + character(LEN=1):: species + real:: aval, Cc, diff, Re, Sc, St, St2, vt, Eff + real, parameter:: boltzman = 1.3806503E-23 + real, parameter:: meanPath = 0.0256E-6 + + vt = 1. + if (species .eq. 'r') then + vt = -0.1021 + 4.932E3*D - 0.9551E6*D*D & + + 0.07934E9*D*D*D - 0.002362E12*D*D*D*D + elseif (species .eq. 's') then + vt = av_s*D**bv_s + elseif (species .eq. 'g') then + vt = av_g*D**bv_g + endif + + Cc = 1. + 2.*meanPath/Da *(1.257+0.4*exp(-0.55*Da/meanPath)) + diff = boltzman*Temp*Cc/(3.*PI*visc*Da) + + Re = 0.5*rhoa*D*vt/visc + Sc = visc/(rhoa*diff) + + St = Da*Da*vt*1000./(9.*visc*D) + aval = 1.+LOG(1.+Re) + St2 = (1.2 + 1./12.*aval)/(1.+aval) + + Eff = 4./(Re*Sc) * (1. + 0.4*SQRT(Re)*Sc**0.3333 & + + 0.16*SQRT(Re)*SQRT(Sc)) & + + 4.*Da/D * (0.02 + Da/D*(1.+2.*SQRT(Re))) + + if (St.gt.St2) Eff = Eff + ( (St-St2)/(St-St2+0.666667))**1.5 + Eff_aero = MAX(1.E-5, MIN(Eff, 1.0)) + + end function Eff_aero + +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Integrate rain size distribution from zero to D-star to compute the +!.. number of drops smaller than D-star that evaporate in a single +!.. timestep. Drops larger than D-star dont evaporate entirely so do +!.. not affect number concentration. +!+---+-----------------------------------------------------------------+ + + subroutine table_dropEvap + + implicit none + +!..Local variables + INTEGER:: i, j, k, n + DOUBLE PRECISION, DIMENSION(nbc):: N_c, massc + DOUBLE PRECISION:: summ, summ2, lamc, N0_c + INTEGER:: nu_c +! DOUBLE PRECISION:: Nt_r, N0, lam_exp, lam +! REAL:: xlimit_intg + + do n = 1, nbc + massc(n) = am_r*Dc(n)**bm_r + enddo + + do k = 1, nbc + nu_c = MIN(15, NINT(1000.E6/t_Nc(k)) + 2) + do j = 1, ntb_c + lamc = (t_Nc(k)*am_r* ccg(2,nu_c)*ocg1(nu_c) / r_c(j))**obmr + N0_c = t_Nc(k)*ocg1(nu_c) * lamc**cce(1,nu_c) + do i = 1, nbc +!-GT tnc_wev(i,j,k) = GAMMP(nu_c+1., SNGL(Dc(i)*lamc))*t_Nc(k) + N_c(i) = N0_c* Dc(i)**nu_c*EXP(-lamc*Dc(i))*dtc(i) +! if(j.eq.18 .and. k.eq.50) print*, ' N_c = ', N_c(i) + summ = 0. + summ2 = 0. + do n = 1, i + summ = summ + massc(n)*N_c(n) + summ2 = summ2 + N_c(n) + enddo +! if(j.eq.18 .and. k.eq.50) print*, ' DEBUG-TABLE: ', r_c(j), t_Nc(k), summ2, summ + tpc_wev(i,j,k) = summ + tnc_wev(i,j,k) = summ2 + enddo + enddo + enddo + +! +!..To do the same thing for rain. +! +! do k = 1, ntb_r +! do j = 1, ntb_r1 +! lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(k))**ore1 +! lam = lam_exp * (crg(3)*org2*org1)**obmr +! N0 = N0r_exp(j)/(crg(2)*lam_exp) * lam**cre(2) +! Nt_r = N0 * crg(2) / lam**cre(2) +! do i = 1, nbr +! xlimit_intg = lam*Dr(i) +! tnr_rev(i,j,k) = GAMMP(mu_r+1.0, xlimit_intg) * Nt_r +! enddo +! enddo +! enddo + +! TO APPLY TABLE ABOVE +!..Rain lookup table indexes. +! Dr_star = DSQRT(-2.D0*DT * t1_evap/(2.*PI) & +! * 0.78*4.*diffu(k)*xsat*rvs/rho_w) +! idx_d = NINT(1.0 + FLOAT(nbr) * DLOG(Dr_star/D0r) & +! / DLOG(Dr(nbr)/D0r)) +! idx_d = MAX(1, MIN(idx_d, nbr)) +! +! nir = NINT(ALOG10(rr(k))) +! do nn = nir-1, nir+1 +! n = nn +! if ( (rr(k)/10.**nn).ge.1.0 .and. & +! (rr(k)/10.**nn).lt.10.0) goto 154 +! enddo +!154 continue +! idx_r = INT(rr(k)/10.**n) + 10*(n-nir2) - (n-nir2) +! idx_r = MAX(1, MIN(idx_r, ntb_r)) +! +! lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr +! lam_exp = lamr * (crg(3)*org2*org1)**bm_r +! N0_exp = org1*rr(k)/am_r * lam_exp**cre(1) +! nir = NINT(DLOG10(N0_exp)) +! do nn = nir-1, nir+1 +! n = nn +! if ( (N0_exp/10.**nn).ge.1.0 .and. & +! (N0_exp/10.**nn).lt.10.0) goto 155 +! enddo +!155 continue +! idx_r1 = INT(N0_exp/10.**n) + 10*(n-nir3) - (n-nir3) +! idx_r1 = MAX(1, MIN(idx_r1, ntb_r1)) +! +! pnr_rev(k) = MIN(nr(k)*odts, SNGL(tnr_rev(idx_d,idx_r1,idx_r) & ! RAIN2M +! * odts)) + + end subroutine table_dropEvap +! +!ctrlL +!+---+-----------------------------------------------------------------+ +!..Fill the table of CCN activation data created from parcel model run +!.. by Trude Eidhammer with inputs of aerosol number concentration, +!.. vertical velocity, temperature, lognormal mean aerosol radius, and +!.. hygroscopicity, kappa. The data are read from external file and +!.. contain activated fraction of CCN for given conditions. +!+---+-----------------------------------------------------------------+ + + subroutine table_ccnAct + +! jwb USE module_domain +! jwb USE module_dm + implicit none + +! LOGICAL, EXTERNAL:: wrf_dm_on_monitor + +!..Local variables + INTEGER:: iunit_mp_th1, i + LOGICAL:: opened + CHARACTER*64 errmess + + iunit_mp_th1 = -1 +! IF ( wrf_dm_on_monitor() ) THEN + DO i = 20,99 + INQUIRE ( i , OPENED = opened ) + IF ( .NOT. opened ) THEN + iunit_mp_th1 = i + GOTO 2010 + ENDIF + ENDDO + 2010 CONTINUE +! ENDIF +!#if defined(DM_PARALLEL) && !defined(STUBMPI) +! CALL wrf_dm_bcast_bytes ( iunit_mp_th1 , IWORDSIZE ) +!#endif + IF ( iunit_mp_th1 < 0 ) THEN + write(0,*) 'module_mp_thompson: table_ccnAct: '// & + 'Can not find unused fortran unit to read in lookup table.' ! DH* errmsg + return + ENDIF + +! IF ( wrf_dm_on_monitor() ) THEN + WRITE(*, '(A,I2)') 'module_mp_thompson: opening CCN_ACTIVATE.BIN on unit ',iunit_mp_th1 + OPEN(iunit_mp_th1,FILE='CCN_ACTIVATE.BIN', & + FORM='UNFORMATTED',STATUS='OLD',ERR=9009) +! ENDIF + +!!#define DM_BCAST_MACRO(A) CALL wrf_dm_bcast_bytes(A, size(A)*R4SIZE) + +!sms$serial begin + READ(iunit_mp_th1,ERR=9010) tnccn_act +!sms$serial end + +!#if defined(DM_PARALLEL) && !defined(STUBMPI) +! DM_BCAST_MACRO(tnccn_act) +!#endif + + + RETURN + 9009 CONTINUE + WRITE( errmess , '(A,I2)' ) 'module_mp_thompson: error opening CCN_ACTIVATE.BIN on unit ',iunit_mp_th1 + write(0,*) errmess ! DH* errmsg + RETURN + 9010 CONTINUE + WRITE( errmess , '(A,I2)' ) 'module_mp_thompson: error reading CCN_ACTIVATE.BIN on unit ',iunit_mp_th1 + write(0,*) errmess ! DH* errmsg + RETURN + + end subroutine table_ccnAct +!^L +!+---+-----------------------------------------------------------------+ +!..Retrieve fraction of CCN that gets activated given the model temp, +!.. vertical velocity, and available CCN concentration. The lookup +!.. table (read from external file) has CCN concentration varying the +!.. quickest, then updraft, then temperature, then mean aerosol radius, +!.. and finally hygroscopicity, kappa. +!.. TO_DO ITEM: For radiation cooling producing fog, in which case the +!.. updraft velocity could easily be negative, we could use the temp +!.. and its tendency to diagnose a pretend postive updraft velocity. +!+---+-----------------------------------------------------------------+ + real function activ_ncloud(Tt, Ww, NCCN) + + implicit none + REAL, INTENT(IN):: Tt, Ww, NCCN + REAL:: n_local, w_local + INTEGER:: i, j, k, l, m, n + REAL:: A, B, C, D, t, u, x1, x2, y1, y2, nx, wy, fraction + + +! ta_Na = (/10.0, 31.6, 100.0, 316.0, 1000.0, 3160.0, 10000.0/) ntb_arc +! ta_Ww = (/0.01, 0.0316, 0.1, 0.316, 1.0, 3.16, 10.0, 31.6, 100.0/) ntb_arw +! ta_Tk = (/243.15, 253.15, 263.15, 273.15, 283.15, 293.15, 303.15/) ntb_art +! ta_Ra = (/0.01, 0.02, 0.04, 0.08, 0.16/) ntb_arr +! ta_Ka = (/0.2, 0.4, 0.6, 0.8/) ntb_ark + + n_local = NCCN * 1.E-6 + w_local = Ww + + if (n_local .ge. ta_Na(ntb_arc)) then + n_local = ta_Na(ntb_arc) - 1.0 + elseif (n_local .le. ta_Na(1)) then + n_local = ta_Na(1) + 1.0 + endif + do n = 2, ntb_arc + if (n_local.ge.ta_Na(n-1) .and. n_local.lt.ta_Na(n)) goto 8003 + enddo + 8003 continue + i = n + x1 = LOG(ta_Na(i-1)) + x2 = LOG(ta_Na(i)) + + if (w_local .ge. ta_Ww(ntb_arw)) then + w_local = ta_Ww(ntb_arw) - 1.0 + elseif (w_local .le. ta_Ww(1)) then + w_local = ta_Ww(1) + 0.001 + endif + do n = 2, ntb_arw + if (w_local.ge.ta_Ww(n-1) .and. w_local.lt.ta_Ww(n)) goto 8005 + enddo + 8005 continue + j = n + y1 = LOG(ta_Ww(j-1)) + y2 = LOG(ta_Ww(j)) + + k = MAX(1, MIN( NINT( (Tt - ta_Tk(1))*0.1) + 1, ntb_art)) + +!..The next two values are indexes of mean aerosol radius and +!.. hygroscopicity. Currently these are constant but a future version +!.. should implement other variables to allow more freedom such as +!.. at least simple separation of tiny size sulfates from larger +!.. sea salts. + l = 3 + m = 2 + + A = tnccn_act(i-1,j-1,k,l,m) + B = tnccn_act(i,j-1,k,l,m) + C = tnccn_act(i,j,k,l,m) + D = tnccn_act(i-1,j,k,l,m) + nx = LOG(n_local) + wy = LOG(w_local) + + t = (nx-x1)/(x2-x1) + u = (wy-y1)/(y2-y1) + +! t = (n_local-ta(Na(i-1))/(ta_Na(i)-ta_Na(i-1)) +! u = (w_local-ta_Ww(j-1))/(ta_Ww(j)-ta_Ww(j-1)) + + fraction = (1.0-t)*(1.0-u)*A + t*(1.0-u)*B + t*u*C + (1.0-t)*u*D + +! if (NCCN*fraction .gt. 0.75*Nt_c_max) then +! write(*,*) ' DEBUG-GT ', n_local, w_local, Tt, i, j, k +! endif + + activ_ncloud = NCCN*fraction + + end function activ_ncloud + +!+---+-----------------------------------------------------------------+ +!+---+-----------------------------------------------------------------+ + SUBROUTINE GCF(GAMMCF,A,X,GLN) +! --- RETURNS THE INCOMPLETE GAMMA FUNCTION Q(A,X) EVALUATED BY ITS +! --- CONTINUED FRACTION REPRESENTATION AS GAMMCF. ALSO RETURNS +! --- LN(GAMMA(A)) AS GLN. THE CONTINUED FRACTION IS EVALUATED BY +! --- A MODIFIED LENTZ METHOD. +! --- USES GAMMLN + IMPLICIT NONE + INTEGER, PARAMETER:: ITMAX=100 + REAL, PARAMETER:: gEPS=3.E-7 + REAL, PARAMETER:: FPMIN=1.E-30 + REAL, INTENT(IN):: A, X + REAL:: GAMMCF,GLN + INTEGER:: I + REAL:: AN,B,C,D,DEL,H + GLN=GAMMLN(A) + B=X+1.-A + C=1./FPMIN + D=1./B + H=D + DO 11 I=1,ITMAX + AN=-I*(I-A) + B=B+2. + D=AN*D+B + IF(ABS(D).LT.FPMIN)D=FPMIN + C=B+AN/C + IF(ABS(C).LT.FPMIN)C=FPMIN + D=1./D + DEL=D*C + H=H*DEL + IF(ABS(DEL-1.).LT.gEPS)GOTO 1 + 11 CONTINUE + PRINT *, 'A TOO LARGE, ITMAX TOO SMALL IN GCF' + 1 GAMMCF=EXP(-X+A*LOG(X)-GLN)*H + END SUBROUTINE GCF +! (C) Copr. 1986-92 Numerical Recipes Software 2.02 +!+---+-----------------------------------------------------------------+ + SUBROUTINE GSER(GAMSER,A,X,GLN) +! --- RETURNS THE INCOMPLETE GAMMA FUNCTION P(A,X) EVALUATED BY ITS +! --- ITS SERIES REPRESENTATION AS GAMSER. ALSO RETURNS LN(GAMMA(A)) +! --- AS GLN. +! --- USES GAMMLN + IMPLICIT NONE + INTEGER, PARAMETER:: ITMAX=100 + REAL, PARAMETER:: gEPS=3.E-7 + REAL, INTENT(IN):: A, X + REAL:: GAMSER,GLN + INTEGER:: N + REAL:: AP,DEL,SUM + GLN=GAMMLN(A) + IF(X.LE.0.)THEN + IF(X.LT.0.) PRINT *, 'X < 0 IN GSER' + GAMSER=0. + RETURN + ENDIF + AP=A + SUM=1./A + DEL=SUM + DO 11 N=1,ITMAX + AP=AP+1. + DEL=DEL*X/AP + SUM=SUM+DEL + IF(ABS(DEL).LT.ABS(SUM)*gEPS)GOTO 1 + 11 CONTINUE + PRINT *,'A TOO LARGE, ITMAX TOO SMALL IN GSER' + 1 GAMSER=SUM*EXP(-X+A*LOG(X)-GLN) + END SUBROUTINE GSER +! (C) Copr. 1986-92 Numerical Recipes Software 2.02 +!+---+-----------------------------------------------------------------+ + REAL FUNCTION GAMMLN(XX) +! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 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) + +! print *,'rslfmp',p,t,x +! 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/max(1.e-4,(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 = 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, is_aerosol_aware) + + 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 + LOGICAL, INTENT(IN) :: is_aerosol_aware +!..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/) + + has_qc = .false. + has_qi = .false. + has_qs = .false. + + 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(R2, nc1d(k)*rho(k)) + 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 + + if (has_qc) then + do k = kts, kte + if (rc(k).le.R1 .or. nc(k).le.R2) CYCLE + if (nc(k).lt.100) then + inu_c = 15 + elseif (nc(k).gt.1.E10) then + inu_c = 2 + else + inu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) + endif + lamc = (nc(k)*am_r*g_ratio(inu_c)/rc(k))**obmr + re_qc1d(k) = MAX(2.51E-6, MIN(SNGL(0.5D0 * DBLE(3.+inu_c)/lamc), 50.E-6)) + enddo + endif + + 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(5.01E-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(10.E-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. +!+---+-----------------------------------------------------------------+ + +! dcd added vt_dBZ, first_time_step + subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & + t1d, p1d, dBZ, vt_dBZ, kts, kte, ii, jj, first_time_step) + + IMPLICIT NONE + +!..Sub arguments + INTEGER, INTENT(IN):: kts, kte, ii, jj + REAL, DIMENSION(kts:kte), INTENT(IN):: & + qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, t1d, p1d + REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ + REAL, DIMENSION(kts:kte), INTENT(INOUT):: vt_dBZ + LOGICAL, INTENT(IN) :: first_time_step ! dcd + +!..Local variables + LOGICAL :: allow_wet_graupel ! dcd + LOGICAL :: allow_wet_snow ! dcd + REAL, DIMENSION(kts:kte):: temp, pres, qv, rho, rhof + REAL, DIMENSION(kts:kte):: rc, rr, nr, rs, rg + + 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 + + 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 + LOGICAL:: melti + LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg + + DOUBLE PRECISION:: cback, x, eta, f_d + REAL:: xslw1, ygra1, zans1 + +!+---+ + +! dcd + if (first_time_step) then +! no bright banding, to be consistent with hydrometeor retrieval in GSI + allow_wet_snow = .false. + else + allow_wet_snow = .true. + endif + allow_wet_graupel = .false. +! print*, 'allow_wet_snow = ', allow_wet_snow +! print*, 'allow_wet_graupel = ', allow_wet_graupel +! dcd + + 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) + L_qg(k) = .true. + else + rg(k) = R1 + L_qg(k) = .false. + endif + enddo + +!+---+-----------------------------------------------------------------+ +!..Calculate y-intercept, slope, and useful moments for snow. +!+---+-----------------------------------------------------------------+ + do k = kts, kte + 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 + +!+---+-----------------------------------------------------------------+ +!..Calculate y-intercept, slope values for graupel. +!+---+-----------------------------------------------------------------+ + + N0_min = gonv_max + k_0 = kts + do k = kte, kts, -1 + if (temp(k).ge.270.65) k_0 = MAX(k_0, k) + enddo + do k = kte, kts, -1 + if (k.gt.k_0 .and. L_qr(k) .and. mvd_r(k).gt.100.E-6) then + xslw1 = 4.01 + alog10(mvd_r(k)) + else + xslw1 = 0.01 + endif + ygra1 = 4.31 + alog10(max(5.E-5, rg(k))) + zans1 = 3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1)) + N0_exp = 10.**(zans1) + N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) + N0_min = MIN(N0_exp, N0_min) + N0_exp = N0_min + lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 + lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg + ilamg(k) = 1./lamg + N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2) + enddo + +!+---+-----------------------------------------------------------------+ +!..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 + +!+---+-----------------------------------------------------------------+ +!..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, 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/900.0)*(am_g/900.0) & + & * N0_g(k)*cgg(4)*ilamg(k)**cge(4) + 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_0-1, kts, -1 + +!..Reflectivity contributed by melting snow +! dcd added allow_wet_snow + if (allow_wet_snow .and. 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 +! dcd added allow_wet_graupel + if (allow_wet_graupel .and. 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 * xxDg(n)**bm_g + call rayleigh_soak_wetgraupel (x, DBLE(ocmg), 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 = kte, 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*cgg(5)*lamg**(-cge(5)) & +! & / (cgg(4)*lamg**(-cge(4))) +! 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 +! + +#ifdef SION + subroutine readwrite_tables(mode, mpicomm, mpirank, mpiroot, ierr) + +#ifdef MPI + use mpi +#endif + use sion_f90 + + implicit none + + ! Interface variables + character(len=*), intent(in) :: mode + integer, intent(in) :: mpicomm + integer, intent(in) :: mpirank + integer, intent(in) :: mpiroot + integer, intent(out) :: ierr + +#ifdef MPI + ! MPI variables + integer :: mpierr +#endif + + ! SIONlib variables + integer :: SIONLIB_fsblksize + integer :: SIONLIB_numfiles + character*2 :: SIONLIB_filemode + ! + integer :: nprocs + integer, dimension(:), allocatable :: procs + integer*8, dimension(:), allocatable :: chunksizes + ! + integer*8 :: brw + integer :: sid + integer :: f_endian, s_endian + logical :: exists + integer*8 :: tables_size + real*8 :: checksum + character(len=*), parameter :: filename = 'thompson_tables_precomp.sl' + + integer :: i + + continue + + ierr = 0 + + ! Test if SIONlib file containing pre-computed tables exists + inquire(file=trim(filename), exist=exists) + if (trim(mode)=="read") then + SIONLIB_filemode = "rb" + if (.not.exists) then + if (mpirank==mpiroot) write(0,*) "SIONlib file " // trim(filename) // & + " with precomputed Thompson MP tables not found" + ierr = 1 + return + end if + else if (trim(mode)=="write") then + SIONLIB_filemode = "wb" + SIONLIB_numfiles = 1 + if (exists) then + if (mpirank==mpiroot) write(0,*) "SIONlib file " // trim(filename) // & + " with precomputed Thompson MP tables already exists" + ierr = 1 + return + end if + end if + +#ifdef MPI + ! To avoid that MPI master task creates the file before + ! other tasks pass the inquire test above + call MPI_BARRIER(mpicomm, mpierr) +#endif + + mpi_master_io_only: if (mpirank==mpiroot) then + tables_size = sizeof(tcg_racg) + tables_size = tables_size + sizeof(tmr_racg) + tables_size = tables_size + sizeof(tcr_gacr) + tables_size = tables_size + sizeof(tmg_gacr) + tables_size = tables_size + sizeof(tnr_racg) + tables_size = tables_size + sizeof(tnr_gacr) + tables_size = tables_size + sizeof(tcs_racs1) + tables_size = tables_size + sizeof(tmr_racs1) + tables_size = tables_size + sizeof(tcs_racs2) + tables_size = tables_size + sizeof(tmr_racs2) + tables_size = tables_size + sizeof(tcr_sacr1) + tables_size = tables_size + sizeof(tms_sacr1) + tables_size = tables_size + sizeof(tcr_sacr2) + tables_size = tables_size + sizeof(tms_sacr2) + tables_size = tables_size + sizeof(tnr_racs1) + tables_size = tables_size + sizeof(tnr_racs2) + tables_size = tables_size + sizeof(tnr_sacr1) + tables_size = tables_size + sizeof(tnr_sacr2) + tables_size = tables_size + sizeof(tpi_qcfz) + tables_size = tables_size + sizeof(tni_qcfz) + tables_size = tables_size + sizeof(tpi_qrfz) + tables_size = tables_size + sizeof(tpg_qrfz) + tables_size = tables_size + sizeof(tni_qrfz) + tables_size = tables_size + sizeof(tnr_qrfz) + tables_size = tables_size + sizeof(tps_iaus) + tables_size = tables_size + sizeof(tni_iaus) + tables_size = tables_size + sizeof(tpi_ide) + tables_size = tables_size + sizeof(t_Efrw) + tables_size = tables_size + sizeof(t_Efsw) + tables_size = tables_size + sizeof(tnr_rev) + tables_size = tables_size + sizeof(tpc_wev) + tables_size = tables_size + sizeof(tnc_wev) + tables_size = tables_size + sizeof(tnccn_act) + + ! Autodetect SIONlib filesystem block size + SIONLIB_fsblksize = -1 + + nprocs = 1 + allocate (procs(1:nprocs)) + allocate (chunksizes(1:nprocs)) + do i=1,nprocs + procs(i) = i + chunksizes(i) = sizeof(checksum) + tables_size + end do + + write(0,'(a)') "Opening file " // trim(filename) + call fsion_open(trim(filename), SIONLIB_filemode, nprocs, SIONLIB_numfiles, chunksizes(1), SIONLIB_fsblksize, procs(1), sid) + if (sid<0) write(0,'(a)') "Error opening " // trim(filename) // " in " // trim(mode) // " mode" + + call fsion_seek(sid, mpirank, SION_CURRENT_BLK, SION_CURRENT_POS, ierr) + ! fsion_seek returns ierr=1 if cursor could be positioned as requested and 0 otherwise + if (ierr==1) ierr=0 + + if (trim(mode)=="read") then + ! Check that file endianness is identical to system endianness + call fsion_get_file_endianness(sid, f_endian) + call fsion_get_endianess(s_endian) + if (f_endian .ne. s_endian) then + write(0,'(a)') "Error, endianness of SIONlib file " // trim(filename) // " differs " // & + "from filesystem endianness; please delete file and recalculate tables!" + ierr = 1 + end if + if (ierr==0) then + ! Read checksum + call fsion_read(checksum, int(kind(checksum),8), int(1,8), sid, brw) + ! Read arrays tcg_racg through tnccn_act + call fsion_read(tcg_racg(1,1,1,1), int(kind(tcg_racg(1,1,1,1)),8), int(size(tcg_racg),8), sid, brw) + call fsion_read(tmr_racg(1,1,1,1), int(kind(tmr_racg(1,1,1,1)),8), int(size(tmr_racg),8), sid, brw) + call fsion_read(tcr_gacr(1,1,1,1), int(kind(tcr_gacr(1,1,1,1)),8), int(size(tcr_gacr),8), sid, brw) + call fsion_read(tmg_gacr(1,1,1,1), int(kind(tmg_gacr(1,1,1,1)),8), int(size(tmg_gacr),8), sid, brw) + call fsion_read(tnr_racg(1,1,1,1), int(kind(tnr_racg(1,1,1,1)),8), int(size(tnr_racg),8), sid, brw) + call fsion_read(tnr_gacr(1,1,1,1), int(kind(tnr_gacr(1,1,1,1)),8), int(size(tnr_gacr),8), sid, brw) + call fsion_read(tcs_racs1(1,1,1,1), int(kind(tcs_racs1(1,1,1,1)),8), int(size(tcs_racs1),8), sid, brw) + call fsion_read(tmr_racs1(1,1,1,1), int(kind(tmr_racs1(1,1,1,1)),8), int(size(tmr_racs1),8), sid, brw) + call fsion_read(tcs_racs2(1,1,1,1), int(kind(tcs_racs2(1,1,1,1)),8), int(size(tcs_racs2),8), sid, brw) + call fsion_read(tmr_racs2(1,1,1,1), int(kind(tmr_racs2(1,1,1,1)),8), int(size(tmr_racs2),8), sid, brw) + call fsion_read(tcr_sacr1(1,1,1,1), int(kind(tcr_sacr1(1,1,1,1)),8), int(size(tcr_sacr1),8), sid, brw) + call fsion_read(tms_sacr1(1,1,1,1), int(kind(tms_sacr1(1,1,1,1)),8), int(size(tms_sacr1),8), sid, brw) + call fsion_read(tcr_sacr2(1,1,1,1), int(kind(tcr_sacr2(1,1,1,1)),8), int(size(tcr_sacr2),8), sid, brw) + call fsion_read(tms_sacr2(1,1,1,1), int(kind(tms_sacr2(1,1,1,1)),8), int(size(tms_sacr2),8), sid, brw) + call fsion_read(tnr_racs1(1,1,1,1), int(kind(tnr_racs1(1,1,1,1)),8), int(size(tnr_racs1),8), sid, brw) + call fsion_read(tnr_racs2(1,1,1,1), int(kind(tnr_racs2(1,1,1,1)),8), int(size(tnr_racs2),8), sid, brw) + call fsion_read(tnr_sacr1(1,1,1,1), int(kind(tnr_sacr1(1,1,1,1)),8), int(size(tnr_sacr1),8), sid, brw) + call fsion_read(tnr_sacr2(1,1,1,1), int(kind(tnr_sacr2(1,1,1,1)),8), int(size(tnr_sacr2),8), sid, brw) + call fsion_read(tpi_qcfz(1,1,1,1), int(kind(tpi_qcfz(1,1,1,1)),8), int(size(tpi_qcfz),8), sid, brw) + call fsion_read(tni_qcfz(1,1,1,1), int(kind(tni_qcfz(1,1,1,1)),8), int(size(tni_qcfz),8), sid, brw) + call fsion_read(tpi_qrfz(1,1,1,1), int(kind(tpi_qrfz(1,1,1,1)),8), int(size(tpi_qrfz),8), sid, brw) + call fsion_read(tpg_qrfz(1,1,1,1), int(kind(tpg_qrfz(1,1,1,1)),8), int(size(tpg_qrfz),8), sid, brw) + call fsion_read(tni_qrfz(1,1,1,1), int(kind(tni_qrfz(1,1,1,1)),8), int(size(tni_qrfz),8), sid, brw) + call fsion_read(tnr_qrfz(1,1,1,1), int(kind(tnr_qrfz(1,1,1,1)),8), int(size(tnr_qrfz),8), sid, brw) + call fsion_read(tps_iaus(1,1), int(kind(tps_iaus(1,1)),8), int(size(tps_iaus),8), sid, brw) + call fsion_read(tni_iaus(1,1), int(kind(tni_iaus(1,1)),8), int(size(tni_iaus),8), sid, brw) + call fsion_read(tpi_ide(1,1), int(kind(tpi_ide(1,1)),8), int(size(tpi_ide),8), sid, brw) + call fsion_read(t_Efrw(1,1), int(kind(t_Efrw(1,1)),8), int(size(t_Efrw),8), sid, brw) + call fsion_read(t_Efsw(1,1), int(kind(t_Efsw(1,1)),8), int(size(t_Efsw),8), sid, brw) + call fsion_read(tnr_rev(1,1,1), int(kind(tnr_rev(1,1,1)),8), int(size(tnr_rev),8), sid, brw) + call fsion_read(tpc_wev(1,1,1), int(kind(tpc_wev(1,1,1)),8), int(size(tpc_wev),8), sid, brw) + call fsion_read(tnc_wev(1,1,1), int(kind(tnc_wev (1,1,1)),8), int(size(tnc_wev),8), sid, brw) + call fsion_read(tnccn_act(1,1,1,1,1), int(kind(tnccn_act(1,1,1,1,1)),8), int(size(tnccn_act),8), sid, brw) + else + ! Wrong endianness (ierr/=0) will force checksum match to fail + checksum = -1 + end if + else if (trim(mode)=="write") then + ! Calculate and write checksum + checksum = calculate_checksum() + call fsion_write(checksum, int(kind(checksum),8), int(1,8), sid, brw) + ! Write arrays tcg_racg through tnccn_act + call fsion_write(tcg_racg(1,1,1,1), int(kind(tcg_racg(1,1,1,1)),8), int(size(tcg_racg),8), sid, brw) + call fsion_write(tmr_racg(1,1,1,1), int(kind(tmr_racg(1,1,1,1)),8), int(size(tmr_racg),8), sid, brw) + call fsion_write(tcr_gacr(1,1,1,1), int(kind(tcr_gacr(1,1,1,1)),8), int(size(tcr_gacr),8), sid, brw) + call fsion_write(tmg_gacr(1,1,1,1), int(kind(tmg_gacr(1,1,1,1)),8), int(size(tmg_gacr),8), sid, brw) + call fsion_write(tnr_racg(1,1,1,1), int(kind(tnr_racg(1,1,1,1)),8), int(size(tnr_racg),8), sid, brw) + call fsion_write(tnr_gacr(1,1,1,1), int(kind(tnr_gacr(1,1,1,1)),8), int(size(tnr_gacr),8), sid, brw) + call fsion_write(tcs_racs1(1,1,1,1), int(kind(tcs_racs1(1,1,1,1)),8), int(size(tcs_racs1),8), sid, brw) + call fsion_write(tmr_racs1(1,1,1,1), int(kind(tmr_racs1(1,1,1,1)),8), int(size(tmr_racs1),8), sid, brw) + call fsion_write(tcs_racs2(1,1,1,1), int(kind(tcs_racs2(1,1,1,1)),8), int(size(tcs_racs2),8), sid, brw) + call fsion_write(tmr_racs2(1,1,1,1), int(kind(tmr_racs2(1,1,1,1)),8), int(size(tmr_racs2),8), sid, brw) + call fsion_write(tcr_sacr1(1,1,1,1), int(kind(tcr_sacr1(1,1,1,1)),8), int(size(tcr_sacr1),8), sid, brw) + call fsion_write(tms_sacr1(1,1,1,1), int(kind(tms_sacr1(1,1,1,1)),8), int(size(tms_sacr1),8), sid, brw) + call fsion_write(tcr_sacr2(1,1,1,1), int(kind(tcr_sacr2(1,1,1,1)),8), int(size(tcr_sacr2),8), sid, brw) + call fsion_write(tms_sacr2(1,1,1,1), int(kind(tms_sacr2(1,1,1,1)),8), int(size(tms_sacr2),8), sid, brw) + call fsion_write(tnr_racs1(1,1,1,1), int(kind(tnr_racs1(1,1,1,1)),8), int(size(tnr_racs1),8), sid, brw) + call fsion_write(tnr_racs2(1,1,1,1), int(kind(tnr_racs2(1,1,1,1)),8), int(size(tnr_racs2),8), sid, brw) + call fsion_write(tnr_sacr1(1,1,1,1), int(kind(tnr_sacr1(1,1,1,1)),8), int(size(tnr_sacr1),8), sid, brw) + call fsion_write(tnr_sacr2(1,1,1,1), int(kind(tnr_sacr2(1,1,1,1)),8), int(size(tnr_sacr2),8), sid, brw) + call fsion_write(tpi_qcfz(1,1,1,1), int(kind(tpi_qcfz(1,1,1,1)),8), int(size(tpi_qcfz),8), sid, brw) + call fsion_write(tni_qcfz(1,1,1,1), int(kind(tni_qcfz(1,1,1,1)),8), int(size(tni_qcfz),8), sid, brw) + call fsion_write(tpi_qrfz(1,1,1,1), int(kind(tpi_qrfz(1,1,1,1)),8), int(size(tpi_qrfz),8), sid, brw) + call fsion_write(tpg_qrfz(1,1,1,1), int(kind(tpg_qrfz(1,1,1,1)),8), int(size(tpg_qrfz),8), sid, brw) + call fsion_write(tni_qrfz(1,1,1,1), int(kind(tni_qrfz(1,1,1,1)),8), int(size(tni_qrfz),8), sid, brw) + call fsion_write(tnr_qrfz(1,1,1,1), int(kind(tnr_qrfz(1,1,1,1)),8), int(size(tnr_qrfz),8), sid, brw) + call fsion_write(tps_iaus(1,1), int(kind(tps_iaus(1,1)),8), int(size(tps_iaus),8), sid, brw) + call fsion_write(tni_iaus(1,1), int(kind(tni_iaus(1,1)),8), int(size(tni_iaus),8), sid, brw) + call fsion_write(tpi_ide(1,1), int(kind(tpi_ide(1,1)),8), int(size(tpi_ide),8), sid, brw) + call fsion_write(t_Efrw(1,1), int(kind(t_Efrw(1,1)),8), int(size(t_Efrw),8), sid, brw) + call fsion_write(t_Efsw(1,1), int(kind(t_Efsw(1,1)),8), int(size(t_Efsw),8), sid, brw) + call fsion_write(tnr_rev(1,1,1), int(kind(tnr_rev(1,1,1)),8), int(size(tnr_rev),8), sid, brw) + call fsion_write(tpc_wev(1,1,1), int(kind(tpc_wev(1,1,1)),8), int(size(tpc_wev),8), sid, brw) + call fsion_write(tnc_wev(1,1,1), int(kind(tnc_wev (1,1,1)),8), int(size(tnc_wev),8), sid, brw) + call fsion_write(tnccn_act(1,1,1,1,1), int(kind(tnccn_act(1,1,1,1,1)),8), int(size(tnccn_act),8), sid, brw) + end if + + write(0,'(a)') "Closing file " // trim(filename) + call fsion_close(sid, ierr) + + ierr = 0 + ! Test if checksum matches, this fails if wrong endianness (checksum=-1, see above) + if (trim(mode)=="read" .and. checksum/=calculate_checksum()) then + write(0,'(2(a,e20.9))') "Checksum mismatch, expected", calculate_checksum(), "but got", checksum + ierr = 1 + end if + + deallocate (procs) + deallocate (chunksizes) + + else + + ierr = 0 + + end if mpi_master_io_only + +#ifdef MPI + if (trim(mode)=="read") then + ! After reading the tables, broadcast the information to all MPI tasks. + ! First, broadcast the current error code from MPI master (0 = success) + call MPI_BCAST(ierr, 1, MPI_INTEGER, mpiroot, mpicomm, mpierr) + if (ierr/=0) return + call MPI_BCAST(tcg_racg, size(tcg_racg), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tmr_racg, size(tmr_racg), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tcr_gacr, size(tcr_gacr), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tmg_gacr, size(tmg_gacr), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tnr_racg, size(tnr_racg), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tnr_gacr, size(tnr_gacr), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tcs_racs1, size(tcs_racs1), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tmr_racs1, size(tmr_racs1), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tcs_racs2, size(tcs_racs2), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tmr_racs2, size(tmr_racs2), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tcr_sacr1, size(tcr_sacr1), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tms_sacr1, size(tms_sacr1), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tcr_sacr2, size(tcr_sacr2), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tms_sacr2, size(tms_sacr2), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tnr_racs1, size(tnr_racs1), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tnr_racs2, size(tnr_racs2), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tnr_sacr1, size(tnr_sacr1), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tnr_sacr2, size(tnr_sacr2), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tpi_qcfz, size(tpi_qcfz), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tni_qcfz, size(tni_qcfz), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tpi_qrfz, size(tpi_qrfz), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tpg_qrfz, size(tpg_qrfz), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tni_qrfz, size(tni_qrfz), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tnr_qrfz, size(tnr_qrfz), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tps_iaus, size(tps_iaus), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tni_iaus, size(tni_iaus), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tpi_ide, size(tpi_ide), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(t_Efrw, size(t_Efrw), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(t_Efsw, size(t_Efsw), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tnr_rev, size(tnr_rev), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tpc_wev, size(tpc_wev), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tnc_wev, size(tnc_wev), MPI_DOUBLE_PRECISION, mpiroot, mpicomm, mpierr) + call MPI_BCAST(tnccn_act, size(tnccn_act), MPI_REAL, mpiroot, mpicomm, mpierr) + else if (trim(mode)=="write") then + call MPI_BARRIER(mpicomm, mpierr) + end if +#endif + + return + + contains + + function calculate_checksum() result(checksum) + real*8 :: checksum + checksum = real(tables_size,8)*sum(tcg_racg) + end function calculate_checksum + + end subroutine readwrite_tables +#endif + +!+---+-----------------------------------------------------------------+ +!+---+-----------------------------------------------------------------+ +END MODULE module_mp_thompson_hrrr +!+---+-----------------------------------------------------------------+ diff --git a/physics/module_mp_thompson_hrrr_radar.F90 b/physics/module_mp_thompson_hrrr_radar.F90 new file mode 100644 index 000000000..98f5f3e78 --- /dev/null +++ b/physics/module_mp_thompson_hrrr_radar.F90 @@ -0,0 +1,612 @@ +!+---+-----------------------------------------------------------------+ +!..This set of routines facilitates computing radar reflectivity. +!.. This module is more library code whereas the individual microphysics +!.. schemes contains specific details needed for the final computation, +!.. so refer to location within each schemes calling the routine named +!.. rayleigh_soak_wetgraupel. +!.. The bulk of this code originated from Ulrich Blahak (Germany) and +!.. was adapted to WRF by G. Thompson. This version of code is only +!.. intended for use when Rayleigh scattering principles dominate and +!.. is not intended for wavelengths in which Mie scattering is a +!.. significant portion. Therefore, it is well-suited to use with +!.. 5 or 10 cm wavelength like USA NEXRAD radars. +!.. This code makes some rather simple assumptions about water +!.. coating on outside of frozen species (snow/graupel). Fraction of +!.. meltwater is simply the ratio of mixing ratio below melting level +!.. divided by mixing ratio at level just above highest T>0C. Also, +!.. immediately 90% of the melted water exists on the ice's surface +!.. and 10% is embedded within ice. No water is "shed" at all in these +!.. assumptions. The code is quite slow because it does the reflectivity +!.. calculations based on 50 individual size bins of the distributions. +!+---+-----------------------------------------------------------------+ + + MODULE module_mp_thompson_hrrr_radar + + PUBLIC :: rayleigh_soak_wetgraupel + PUBLIC :: radar_init + PRIVATE :: m_complex_water_ray + PRIVATE :: m_complex_ice_maetzler + PRIVATE :: m_complex_maxwellgarnett + PRIVATE :: get_m_mix_nested + PRIVATE :: get_m_mix + PRIVATE :: WGAMMA + PRIVATE :: GAMMLN + + + INTEGER, PARAMETER, PUBLIC:: nrbins = 50 + DOUBLE PRECISION, DIMENSION(nrbins+1), PUBLIC:: xxDx + DOUBLE PRECISION, DIMENSION(nrbins), PUBLIC:: xxDs,xdts,xxDg,xdtg + DOUBLE PRECISION, PARAMETER, PUBLIC:: lamda_radar = 0.10 ! in meters + DOUBLE PRECISION, PUBLIC:: K_w, PI5, lamda4 + COMPLEX*16, PUBLIC:: m_w_0, m_i_0 + DOUBLE PRECISION, DIMENSION(nrbins+1), PUBLIC:: simpson + DOUBLE PRECISION, DIMENSION(3), PARAMETER, PUBLIC:: basis = & + (/1.d0/3.d0, 4.d0/3.d0, 1.d0/3.d0/) + REAL, DIMENSION(4), PUBLIC:: xcre, xcse, xcge, xcrg, xcsg, xcgg + REAL, PUBLIC:: xam_r, xbm_r, xmu_r, xobmr + REAL, PUBLIC:: xam_s, xbm_s, xmu_s, xoams, xobms, xocms + REAL, PUBLIC:: xam_g, xbm_g, xmu_g, xoamg, xobmg, xocmg + REAL, PUBLIC:: xorg2, xosg2, xogg2 + + INTEGER, PARAMETER, PUBLIC:: slen = 20 + CHARACTER(len=slen), PUBLIC:: & + mixingrulestring_s, matrixstring_s, inclusionstring_s, & + hoststring_s, hostmatrixstring_s, hostinclusionstring_s, & + mixingrulestring_g, matrixstring_g, inclusionstring_g, & + hoststring_g, hostmatrixstring_g, hostinclusionstring_g + +!..Single melting snow/graupel particle 90% meltwater on external sfc + DOUBLE PRECISION, PARAMETER:: melt_outside_s = 0.9d0 + DOUBLE PRECISION, PARAMETER:: melt_outside_g = 0.9d0 + + CONTAINS + +!+---+-----------------------------------------------------------------+ +!+---+-----------------------------------------------------------------+ +!+---+-----------------------------------------------------------------+ + + subroutine radar_init + + IMPLICIT NONE + INTEGER:: n + PI5 = 3.14159*3.14159*3.14159*3.14159*3.14159 + lamda4 = lamda_radar*lamda_radar*lamda_radar*lamda_radar + m_w_0 = m_complex_water_ray (lamda_radar, 0.0d0) + m_i_0 = m_complex_ice_maetzler (lamda_radar, 0.0d0) + K_w = (ABS( (m_w_0*m_w_0 - 1.0) /(m_w_0*m_w_0 + 2.0) ))**2 + + do n = 1, nrbins+1 + simpson(n) = 0.0d0 + enddo + do n = 1, nrbins-1, 2 + simpson(n) = simpson(n) + basis(1) + simpson(n+1) = simpson(n+1) + basis(2) + simpson(n+2) = simpson(n+2) + basis(3) + enddo + + do n = 1, slen + mixingrulestring_s(n:n) = char(0) + matrixstring_s(n:n) = char(0) + inclusionstring_s(n:n) = char(0) + hoststring_s(n:n) = char(0) + hostmatrixstring_s(n:n) = char(0) + hostinclusionstring_s(n:n) = char(0) + mixingrulestring_g(n:n) = char(0) + matrixstring_g(n:n) = char(0) + inclusionstring_g(n:n) = char(0) + hoststring_g(n:n) = char(0) + hostmatrixstring_g(n:n) = char(0) + hostinclusionstring_g(n:n) = char(0) + enddo + + mixingrulestring_s = 'maxwellgarnett' + hoststring_s = 'air' + matrixstring_s = 'water' + inclusionstring_s = 'spheroidal' + hostmatrixstring_s = 'icewater' + hostinclusionstring_s = 'spheroidal' + + mixingrulestring_g = 'maxwellgarnett' + hoststring_g = 'air' + matrixstring_g = 'water' + inclusionstring_g = 'spheroidal' + hostmatrixstring_g = 'icewater' + hostinclusionstring_g = 'spheroidal' + +!..Create bins of snow (from 100 microns up to 2 cm). + xxDx(1) = 100.D-6 + xxDx(nrbins+1) = 0.02d0 + do n = 2, nrbins + xxDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nrbins) & + *DLOG(xxDx(nrbins+1)/xxDx(1)) +DLOG(xxDx(1))) + enddo + do n = 1, nrbins + xxDs(n) = DSQRT(xxDx(n)*xxDx(n+1)) + xdts(n) = xxDx(n+1) - xxDx(n) + enddo + +!..Create bins of graupel (from 100 microns up to 5 cm). + xxDx(1) = 100.D-6 + xxDx(nrbins+1) = 0.05d0 + do n = 2, nrbins + xxDx(n) = DEXP(DFLOAT(n-1)/DFLOAT(nrbins) & + *DLOG(xxDx(nrbins+1)/xxDx(1)) +DLOG(xxDx(1))) + enddo + do n = 1, nrbins + xxDg(n) = DSQRT(xxDx(n)*xxDx(n+1)) + xdtg(n) = xxDx(n+1) - xxDx(n) + enddo + + +!..The calling program must set the m(D) relations and gamma shape +!.. parameter mu for rain, snow, and graupel. Easily add other types +!.. based on the template here. For majority of schemes with simpler +!.. exponential number distribution, mu=0. + + xcre(1) = 1. + xbm_r + xcre(2) = 1. + xmu_r + xcre(3) = 1. + xbm_r + xmu_r + xcre(4) = 1. + 2.*xbm_r + xmu_r + do n = 1, 4 + xcrg(n) = WGAMMA(xcre(n)) + enddo + xorg2 = 1./xcrg(2) + + xcse(1) = 1. + xbm_s + xcse(2) = 1. + xmu_s + xcse(3) = 1. + xbm_s + xmu_s + xcse(4) = 1. + 2.*xbm_s + xmu_s + do n = 1, 4 + xcsg(n) = WGAMMA(xcse(n)) + enddo + xosg2 = 1./xcsg(2) + + xcge(1) = 1. + xbm_g + xcge(2) = 1. + xmu_g + xcge(3) = 1. + xbm_g + xmu_g + xcge(4) = 1. + 2.*xbm_g + xmu_g + do n = 1, 4 + xcgg(n) = WGAMMA(xcge(n)) + enddo + xogg2 = 1./xcgg(2) + + xobmr = 1./xbm_r + xoams = 1./xam_s + xobms = 1./xbm_s + xocms = xoams**xobms + xoamg = 1./xam_g + xobmg = 1./xbm_g + xocmg = xoamg**xobmg + + + end subroutine radar_init + +!+---+-----------------------------------------------------------------+ +!+---+-----------------------------------------------------------------+ + + COMPLEX*16 FUNCTION m_complex_water_ray(lambda,T) + +! Complex refractive Index of Water as function of Temperature T +! [deg C] and radar wavelength lambda [m]; valid for +! lambda in [0.001,1.0] m; T in [-10.0,30.0] deg C +! after Ray (1972) + + IMPLICIT NONE + DOUBLE PRECISION, INTENT(IN):: T,lambda + DOUBLE PRECISION:: epsinf,epss,epsr,epsi + DOUBLE PRECISION:: alpha,lambdas,sigma,nenner + COMPLEX*16, PARAMETER:: i = (0d0,1d0) + DOUBLE PRECISION, PARAMETER:: PIx=3.1415926535897932384626434d0 + + epsinf = 5.27137d0 + 0.02164740d0 * T - 0.00131198d0 * T*T + epss = 78.54d+0 * (1.0 - 4.579d-3 * (T - 25.0) & + + 1.190d-5 * (T - 25.0)*(T - 25.0) & + - 2.800d-8 * (T - 25.0)*(T - 25.0)*(T - 25.0)) + alpha = -16.8129d0/(T+273.16) + 0.0609265d0 + lambdas = 0.00033836d0 * exp(2513.98d0/(T+273.16)) * 1e-2 + + nenner = 1.d0+2.d0*(lambdas/lambda)**(1d0-alpha)*sin(alpha*PIx*0.5) & + + (lambdas/lambda)**(2d0-2d0*alpha) + epsr = epsinf + ((epss-epsinf) * ((lambdas/lambda)**(1d0-alpha) & + * sin(alpha*PIx*0.5)+1d0)) / nenner + epsi = ((epss-epsinf) * ((lambdas/lambda)**(1d0-alpha) & + * cos(alpha*PIx*0.5)+0d0)) / nenner & + + lambda*1.25664/1.88496 + + m_complex_water_ray = SQRT(CMPLX(epsr,-epsi)) + + END FUNCTION m_complex_water_ray + +!+---+-----------------------------------------------------------------+ + + COMPLEX*16 FUNCTION m_complex_ice_maetzler(lambda,T) + +! complex refractive index of ice as function of Temperature T +! [deg C] and radar wavelength lambda [m]; valid for +! lambda in [0.0001,30] m; T in [-250.0,0.0] C +! Original comment from the Matlab-routine of Prof. Maetzler: +! Function for calculating the relative permittivity of pure ice in +! the microwave region, according to C. Maetzler, "Microwave +! properties of ice and snow", in B. Schmitt et al. (eds.) Solar +! System Ices, Astrophys. and Space Sci. Library, Vol. 227, Kluwer +! Academic Publishers, Dordrecht, pp. 241-257 (1998). Input: +! TK = temperature (K), range 20 to 273.15 +! f = frequency in GHz, range 0.01 to 3000 + + IMPLICIT NONE + DOUBLE PRECISION, INTENT(IN):: T,lambda + DOUBLE PRECISION:: f,c,TK,B1,B2,b,deltabeta,betam,beta,theta,alfa + + c = 2.99d8 + TK = T + 273.16 + f = c / lambda * 1d-9 + + B1 = 0.0207 + B2 = 1.16d-11 + b = 335.0d0 + deltabeta = EXP(-10.02 + 0.0364*(TK-273.16)) + betam = (B1/TK) * ( EXP(b/TK) / ((EXP(b/TK)-1)**2) ) + B2*f*f + beta = betam + deltabeta + theta = 300. / TK - 1. + alfa = (0.00504d0 + 0.0062d0*theta) * EXP(-22.1d0*theta) + m_complex_ice_maetzler = 3.1884 + 9.1e-4*(TK-273.16) + m_complex_ice_maetzler = m_complex_ice_maetzler & + + CMPLX(0.0d0, (alfa/f + beta*f)) + m_complex_ice_maetzler = SQRT(CONJG(m_complex_ice_maetzler)) + + END FUNCTION m_complex_ice_maetzler + +!+---+-----------------------------------------------------------------+ + + subroutine rayleigh_soak_wetgraupel (x_g, a_geo, b_geo, fmelt, & + meltratio_outside, m_w, m_i, lambda, C_back, & + mixingrule,matrix,inclusion, & + host,hostmatrix,hostinclusion) + + IMPLICIT NONE + + DOUBLE PRECISION, INTENT(in):: x_g, a_geo, b_geo, fmelt, lambda, & + meltratio_outside + DOUBLE PRECISION, INTENT(out):: C_back + COMPLEX*16, INTENT(in):: m_w, m_i + CHARACTER(len=*), INTENT(in):: mixingrule, matrix, inclusion, & + host, hostmatrix, hostinclusion + + COMPLEX*16:: m_core, m_air + DOUBLE PRECISION:: D_large, D_g, rhog, x_w, xw_a, fm, fmgrenz, & + volg, vg, volair, volice, volwater, & + meltratio_outside_grenz, mra + INTEGER:: error + DOUBLE PRECISION, PARAMETER:: PIx=3.1415926535897932384626434d0 + +! refractive index of air: + m_air = (1.0d0,0.0d0) + +! Limiting the degree of melting --- for safety: + fm = DMAX1(DMIN1(fmelt, 1.0d0), 0.0d0) +! Limiting the ratio of (melting on outside)/(melting on inside): + mra = DMAX1(DMIN1(meltratio_outside, 1.0d0), 0.0d0) + +! ! The relative portion of meltwater melting at outside should increase +! ! from the given input value (between 0 and 1) +! ! to 1 as the degree of melting approaches 1, +! ! so that the melting particle "converges" to a water drop. +! ! Simplest assumption is linear: + mra = mra + (1.0d0-mra)*fm + + x_w = x_g * fm + + D_g = a_geo * x_g**b_geo + + if (D_g .ge. 1d-12) then + + vg = PIx/6. * D_g**3 + rhog = DMAX1(DMIN1(x_g / vg, 900.0d0), 10.0d0) + vg = x_g / rhog + + meltratio_outside_grenz = 1.0d0 - rhog / 1000. + + if (mra .le. meltratio_outside_grenz) then + !..In this case, it cannot happen that, during melting, all the + !.. air inclusions within the ice particle get filled with + !.. meltwater. This only happens at the end of all melting. + volg = vg * (1.0d0 - mra * fm) + + else + !..In this case, at some melting degree fm, all the air + !.. inclusions get filled with meltwater. + fmgrenz=(900.0-rhog)/(mra*900.0-rhog+900.0*rhog/1000.) + + if (fm .le. fmgrenz) then + !.. not all air pockets are filled: + volg = (1.0 - mra * fm) * vg + else + !..all air pockets are filled with meltwater, now the + !.. entire ice sceleton melts homogeneously: + volg = (x_g - x_w) / 900.0 + x_w / 1000. + endif + + endif + + D_large = (6.0 / PIx * volg) ** (1./3.) + volice = (x_g - x_w) / (volg * 900.0) + volwater = x_w / (1000. * volg) + volair = 1.0 - volice - volwater + + !..complex index of refraction for the ice-air-water mixture + !.. of the particle: + m_core = get_m_mix_nested (m_air, m_i, m_w, volair, volice, & + volwater, mixingrule, host, matrix, inclusion, & + hostmatrix, hostinclusion, error) + if (error .ne. 0) then + C_back = 0.0d0 + return + endif + + !..Rayleigh-backscattering coefficient of melting particle: + C_back = (ABS((m_core**2-1.0d0)/(m_core**2+2.0d0)))**2 & + * PI5 * D_large**6 / lamda4 + + else + C_back = 0.0d0 + endif + + end subroutine rayleigh_soak_wetgraupel + +!+---+-----------------------------------------------------------------+ + + complex*16 function get_m_mix_nested (m_a, m_i, m_w, volair, & + volice, volwater, mixingrule, host, matrix, & + inclusion, hostmatrix, hostinclusion, cumulerror) + + IMPLICIT NONE + + DOUBLE PRECISION, INTENT(in):: volice, volair, volwater + COMPLEX*16, INTENT(in):: m_a, m_i, m_w + CHARACTER(len=*), INTENT(in):: mixingrule, host, matrix, & + inclusion, hostmatrix, hostinclusion + INTEGER, INTENT(out):: cumulerror + + DOUBLE PRECISION:: vol1, vol2 + COMPLEX*16:: mtmp + INTEGER:: error + + !..Folded: ( (m1 + m2) + m3), where m1,m2,m3 could each be + !.. air, ice, or water + + cumulerror = 0 + get_m_mix_nested = CMPLX(1.0d0,0.0d0) + + if (host .eq. 'air') then + + if (matrix .eq. 'air') then + write(*,*) 'GET_M_MIX_NESTED: bad matrix: ', matrix + cumulerror = cumulerror + 1 + else + vol1 = volice / MAX(volice+volwater,1d-10) + vol2 = 1.0d0 - vol1 + mtmp = get_m_mix (m_a, m_i, m_w, 0.0d0, vol1, vol2, & + mixingrule, matrix, inclusion, error) + cumulerror = cumulerror + error + + if (hostmatrix .eq. 'air') then + get_m_mix_nested = get_m_mix (m_a, mtmp, 2.0*m_a, & + volair, (1.0d0-volair), 0.0d0, mixingrule, & + hostmatrix, hostinclusion, error) + cumulerror = cumulerror + error + elseif (hostmatrix .eq. 'icewater') then + get_m_mix_nested = get_m_mix (m_a, mtmp, 2.0*m_a, & + volair, (1.0d0-volair), 0.0d0, mixingrule, & + 'ice', hostinclusion, error) + cumulerror = cumulerror + error + else + write(*,*) 'GET_M_MIX_NESTED: bad hostmatrix: ', & + hostmatrix + cumulerror = cumulerror + 1 + endif + endif + + elseif (host .eq. 'ice') then + + if (matrix .eq. 'ice') then + write(*,*) 'GET_M_MIX_NESTED: bad matrix: ', matrix + cumulerror = cumulerror + 1 + else + vol1 = volair / MAX(volair+volwater,1d-10) + vol2 = 1.0d0 - vol1 + mtmp = get_m_mix (m_a, m_i, m_w, vol1, 0.0d0, vol2, & + mixingrule, matrix, inclusion, error) + cumulerror = cumulerror + error + + if (hostmatrix .eq. 'ice') then + get_m_mix_nested = get_m_mix (mtmp, m_i, 2.0*m_a, & + (1.0d0-volice), volice, 0.0d0, mixingrule, & + hostmatrix, hostinclusion, error) + cumulerror = cumulerror + error + elseif (hostmatrix .eq. 'airwater') then + get_m_mix_nested = get_m_mix (mtmp, m_i, 2.0*m_a, & + (1.0d0-volice), volice, 0.0d0, mixingrule, & + 'air', hostinclusion, error) + cumulerror = cumulerror + error + else + write(*,*) 'GET_M_MIX_NESTED: bad hostmatrix: ', & + hostmatrix + cumulerror = cumulerror + 1 + endif + endif + + elseif (host .eq. 'water') then + + if (matrix .eq. 'water') then + write(*,*) 'GET_M_MIX_NESTED: bad matrix: ', matrix + cumulerror = cumulerror + 1 + else + vol1 = volair / MAX(volice+volair,1d-10) + vol2 = 1.0d0 - vol1 + mtmp = get_m_mix (m_a, m_i, m_w, vol1, vol2, 0.0d0, & + mixingrule, matrix, inclusion, error) + cumulerror = cumulerror + error + + if (hostmatrix .eq. 'water') then + get_m_mix_nested = get_m_mix (2*m_a, mtmp, m_w, & + 0.0d0, (1.0d0-volwater), volwater, mixingrule, & + hostmatrix, hostinclusion, error) + cumulerror = cumulerror + error + elseif (hostmatrix .eq. 'airice') then + get_m_mix_nested = get_m_mix (2*m_a, mtmp, m_w, & + 0.0d0, (1.0d0-volwater), volwater, mixingrule, & + 'ice', hostinclusion, error) + cumulerror = cumulerror + error + else + write(*,*) 'GET_M_MIX_NESTED: bad hostmatrix: ', & + hostmatrix + cumulerror = cumulerror + 1 + endif + endif + + elseif (host .eq. 'none') then + + get_m_mix_nested = get_m_mix (m_a, m_i, m_w, & + volair, volice, volwater, mixingrule, & + matrix, inclusion, error) + cumulerror = cumulerror + error + + else + write(*,*) 'GET_M_MIX_NESTED: unknown matrix: ', host + cumulerror = cumulerror + 1 + endif + + IF (cumulerror .ne. 0) THEN + write(*,*) 'GET_M_MIX_NESTED: error encountered' + get_m_mix_nested = CMPLX(1.0d0,0.0d0) + endif + + end function get_m_mix_nested + +!+---+-----------------------------------------------------------------+ + + COMPLEX*16 FUNCTION get_m_mix (m_a, m_i, m_w, volair, volice, & + volwater, mixingrule, matrix, inclusion, error) + + IMPLICIT NONE + + DOUBLE PRECISION, INTENT(in):: volice, volair, volwater + COMPLEX*16, INTENT(in):: m_a, m_i, m_w + CHARACTER(len=*), INTENT(in):: mixingrule, matrix, inclusion + INTEGER, INTENT(out):: error + + error = 0 + get_m_mix = CMPLX(1.0d0,0.0d0) + + if (mixingrule .eq. 'maxwellgarnett') then + if (matrix .eq. 'ice') then + get_m_mix = m_complex_maxwellgarnett(volice, volair, volwater, & + m_i, m_a, m_w, inclusion, error) + elseif (matrix .eq. 'water') then + get_m_mix = m_complex_maxwellgarnett(volwater, volair, volice, & + m_w, m_a, m_i, inclusion, error) + elseif (matrix .eq. 'air') then + get_m_mix = m_complex_maxwellgarnett(volair, volwater, volice, & + m_a, m_w, m_i, inclusion, error) + else + write(*,*) 'GET_M_MIX: unknown matrix: ', matrix + error = 1 + endif + + else + write(*,*) 'GET_M_MIX: unknown mixingrule: ', mixingrule + error = 2 + endif + + if (error .ne. 0) then + write(*,*) 'GET_M_MIX: error encountered' + endif + + END FUNCTION get_m_mix + +!+---+-----------------------------------------------------------------+ + + COMPLEX*16 FUNCTION m_complex_maxwellgarnett(vol1, vol2, vol3, & + m1, m2, m3, inclusion, error) + + IMPLICIT NONE + + COMPLEX*16 :: m1, m2, m3 + DOUBLE PRECISION :: vol1, vol2, vol3 + CHARACTER(len=*) :: inclusion + + COMPLEX*16 :: beta2, beta3, m1t, m2t, m3t + INTEGER, INTENT(out) :: error + + error = 0 + + if (DABS(vol1+vol2+vol3-1.0d0) .gt. 1d-6) then + write(*,*) 'M_COMPLEX_MAXWELLGARNETT: sum of the ', & + 'partial volume fractions is not 1...ERROR' + m_complex_maxwellgarnett=CMPLX(-999.99d0,-999.99d0) + error = 1 + return + endif + + m1t = m1**2 + m2t = m2**2 + m3t = m3**2 + + if (inclusion .eq. 'spherical') then + beta2 = 3.0d0*m1t/(m2t+2.0d0*m1t) + beta3 = 3.0d0*m1t/(m3t+2.0d0*m1t) + elseif (inclusion .eq. 'spheroidal') then + beta2 = 2.0d0*m1t/(m2t-m1t) * (m2t/(m2t-m1t)*LOG(m2t/m1t)-1.0d0) + beta3 = 2.0d0*m1t/(m3t-m1t) * (m3t/(m3t-m1t)*LOG(m3t/m1t)-1.0d0) + else + write(*,*) 'M_COMPLEX_MAXWELLGARNETT: ', & + 'unknown inclusion: ', inclusion + m_complex_maxwellgarnett=DCMPLX(-999.99d0,-999.99d0) + error = 1 + return + endif + + m_complex_maxwellgarnett = & + SQRT(((1.0d0-vol2-vol3)*m1t + vol2*beta2*m2t + vol3*beta3*m3t) / & + (1.0d0-vol2-vol3+vol2*beta2+vol3*beta3)) + + END FUNCTION m_complex_maxwellgarnett + +!+---+-----------------------------------------------------------------+ + REAL FUNCTION GAMMLN(XX) +! --- RETURNS THE VALUE LN(GAMMA(XX)) FOR XX > 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 WGAMMA(y) + + IMPLICIT NONE + REAL, INTENT(IN):: y + + WGAMMA = EXP(GAMMLN(y)) + + END FUNCTION WGAMMA + +!+---+-----------------------------------------------------------------+ + END MODULE module_mp_thompson_hrrr_radar +!+---+-----------------------------------------------------------------+ diff --git a/physics/mp_thompson_hrrr.F90 b/physics/mp_thompson_hrrr.F90 new file mode 100644 index 000000000..ddbf478e4 --- /dev/null +++ b/physics/mp_thompson_hrrr.F90 @@ -0,0 +1,557 @@ +! CCPP license goes here, as well as further documentation + +#define DEBUG_AEROSOLS + +module mp_thompson_hrrr + + use machine, only : kind_phys + + use module_mp_thompson_hrrr, only : thompson_init, mp_gt_driver + + implicit none + + contains + +#if 0 +!! \section arg_table_mp_thompson_hrrr_init Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |-----------------|---------------------------------------------------------------|--------------------------------------------------------|------------|------|-----------|-----------|--------|----------| +!! | ncol | horizontal_loop_extent | horizontal loop extent | count | 0 | integer | | in | F | +!! | nlev | vertical_dimension | number of vertical levels | count | 0 | integer | | in | F | +!! | con_g | gravitational_acceleration | gravitational acceleration | m s-2 | 0 | real | kind_phys | in | F | +!! | con_rd | gas_constant_dry_air | ideal gas constant for dry air | J kg-1 K-1 | 0 | real | kind_phys | in | F | +!! | phil | geopotential | geopotential at model layer centers | m2 s-2 | 2 | real | kind_phys | in | F | +!! | prsl | air_pressure | mean layer pressure | Pa | 2 | real | kind_phys | in | F | +!! | tgrs | air_temperature | model layer mean temperature | K | 2 | real | kind_phys | in | F | +!! | is_aerosol_aware| flag_for_aerosol_physics | flag for aerosol-aware physics | flag | 0 | logical | | in | F | +!! | nwfa2d | tendency_of_water_friendly_surface_aerosols_at_surface | instantaneous fake surface aerosol source | kg-1 s-1 | 1 | real | kind_phys | inout | T | +!! | nwfa | water_friendly_aerosol_number_concentration | number concentration of water-friendly aerosols | kg-1 | 2 | real | kind_phys | inout | T | +!! | nifa | ice_friendly_aerosol_number_concentration | number concentration of ice-friendly aerosols | kg-1 | 2 | real | kind_phys | inout | T | +!! | area | cell_area | area of the grid cell | m2 | 1 | real | kind_phys | in | F | +!! | mpicomm | mpi_comm | MPI communicator | index | 0 | integer | | in | F | +!! | mpirank | mpi_rank | current MPI-rank | index | 0 | integer | | in | F | +!! | mpiroot | mpi_root | master MPI-rank | index | 0 | integer | | in | F | +!! | threads | omp_threads | number of OpenMP threads available to scheme | count | 0 | integer | | in | F | +!! | errmsg | error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! +#endif + subroutine mp_thompson_hrrr_init(ncol, nlev, con_g, con_rd, & + phil, prsl, tgrs, & + is_aerosol_aware, & + nwfa2d, nwfa, nifa, & + area, & + mpicomm, mpirank, mpiroot, & + threads, errmsg, errflg) + + implicit none + + ! Interface variables + integer, intent(in) :: ncol + integer, intent(in) :: nlev + + real(kind_phys), intent(in) :: con_g + real(kind_phys), intent(in) :: con_rd + real(kind_phys), intent(in) :: phil(1:ncol,1:nlev) + real(kind_phys), intent(in) :: prsl(1:ncol,1:nlev) + real(kind_phys), intent(in) :: tgrs(1:ncol,1:nlev) + logical, intent(in) :: is_aerosol_aware + real(kind_phys), optional, intent(inout) :: nwfa2d(1:ncol) + real(kind_phys), optional, intent(inout) :: nwfa(1:ncol,1:nlev) + real(kind_phys), optional, intent(inout) :: nifa(1:ncol,1:nlev) + real(kind_phys), intent(in) :: area(1:ncol) + integer, intent(in) :: mpicomm + integer, intent(in) :: mpirank + integer, intent(in) :: mpiroot + integer, intent(in) :: threads + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg + + ! Local variables + real (kind=kind_phys) :: hgt(1:ncol,1:nlev) + real (kind=kind_phys) :: rho(1:ncol,1:nlev) + real(kind_phys) :: dx(1:ncol) + logical, parameter :: is_start = .true. + ! Dimensions used in thompson_init + integer :: ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + ! Geopotential height in m2 s-2 to height in m + hgt = phil/con_g + + ! Set internal dimensions + ids = 1 + ims = 1 + its = 1 + ide = ncol + ime = ncol + ite = ncol + jds = 1 + jms = 1 + jts = 1 + jde = 1 + jme = 1 + jte = 1 + kds = 1 + kms = 1 + kts = 1 + kde = nlev + kme = nlev + kte = nlev + + if (is_aerosol_aware .and. present(nwfa2d) .and. present(nwfa) .and. present(nifa)) then +#ifdef DEBUG_AEROSOLS + if (mpirank==mpiroot) then + write(0,'(a,3e16.7)') "DH DEBUG mp thompson init before: nwfa2d min/mean/max =", & + & minval(nwfa2d), sum(nwfa2d)/real(size(nwfa2d)), maxval(nwfa2d) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson init before: nwfa min/mean/max =", & + & minval(nwfa), sum(nwfa)/real(size(nwfa)), maxval(nwfa) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson init before: nifa min/mean/max =", & + & minval(nifa), sum(nifa)/real(size(nifa)), maxval(nifa) + end if +#endif + ! Grid cell size + dx = sqrt(area) + ! Call init + call thompson_init(hgt=hgt, & + nwfa2d=nwfa2d, nwfa=nwfa, nifa=nifa, & + dx=dx, dy=dx, is_start=is_start, & + ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & + ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & + its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & + mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & + threads=threads) + if (errflg /= 0) return +#ifdef DEBUG_AEROSOLS + if (mpirank==mpiroot) then + write(0,'(a,3e16.7)') "DH DEBUG mp thompson init after: nwfa2d min/mean/max =", & + & minval(nwfa2d), sum(nwfa2d)/real(size(nwfa2d)), maxval(nwfa2d) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson init after: nwfa min/mean/max =", & + & minval(nwfa), sum(nwfa)/real(size(nwfa)), maxval(nwfa) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson init after: nifa min/mean/max =", & + & minval(nifa), sum(nifa)/real(size(nifa)), maxval(nifa) + end if +#endif + else if (is_aerosol_aware) then + write(errmsg,fmt='(*(a))') 'Logic error in mp_thompson_hrrr_init:', & + ' aerosol-aware microphysics require all of the', & + ' following optional arguments: nwfa2d, nwfa, nifa' + errflg = 1 + return + else + call thompson_init(hgt=hgt, & + ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & + ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & + its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & + mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & + threads=threads) + if (errflg /= 0) return + end if + + end subroutine mp_thompson_hrrr_init + +#if 0 +!! \section arg_table_mp_thompson_hrrr_run Argument Table +!! | local_name | standard_name | long_name | units | rank | type | kind | intent | optional | +!! |-----------------|---------------------------------------------------------------|--------------------------------------------------------|------------|------|-----------|-----------|--------|----------| +!! | ncol | horizontal_loop_extent | horizontal loop extent | count | 0 | integer | | in | F | +!! | nlev | vertical_dimension | number of vertical levels | count | 0 | integer | | in | F | +!! | con_g | gravitational_acceleration | gravitational acceleration | m s-2 | 0 | real | kind_phys | in | F | +!! | con_rd | gas_constant_dry_air | ideal gas constant for dry air | J kg-1 K-1 | 0 | real | kind_phys | in | F | +!! | spechum | water_vapor_specific_humidity_updated_by_physics | water vapor specific humidity | kg kg-1 | 2 | real | kind_phys | inout | F | +!! | qc | cloud_condensed_water_mixing_ratio_updated_by_physics | cloud water mixing ratio wrt dry+vapor (no condensates)| kg kg-1 | 2 | real | kind_phys | inout | F | +!! | qr | rain_water_mixing_ratio_updated_by_physics | rain water mixing ratio wrt dry+vapor (no condensates) | kg kg-1 | 2 | real | kind_phys | inout | F | +!! | qi | ice_water_mixing_ratio_updated_by_physics | ice water mixing ratio wrt dry+vapor (no condensates) | kg kg-1 | 2 | real | kind_phys | inout | F | +!! | qs | snow_water_mixing_ratio_updated_by_physics | snow water mixing ratio wrt dry+vapor (no condensates) | kg kg-1 | 2 | real | kind_phys | inout | F | +!! | qg | graupel_mixing_ratio_updated_by_physics | graupel mixing ratio wrt dry+vapor (no condensates) | kg kg-1 | 2 | real | kind_phys | inout | F | +!! | ni | ice_number_concentration_updated_by_physics | ice number concentration | kg-1 | 2 | real | kind_phys | inout | F | +!! | nr | rain_number_concentration_updated_by_physics | rain number concentration | kg-1 | 2 | real | kind_phys | inout | F | +!! | is_aerosol_aware| flag_for_aerosol_physics | flag for aerosol-aware physics | flag | 0 | logical | | in | F | +!! | nc | cloud_droplet_number_concentration_updated_by_physics | cloud droplet number concentration | kg-1 | 2 | real | kind_phys | inout | T | +!! | nwfa | water_friendly_aerosol_number_concentration_updated_by_physics| number concentration of water-friendly aerosols | kg-1 | 2 | real | kind_phys | inout | T | +!! | nifa | ice_friendly_aerosol_number_concentration_updated_by_physics | number concentration of ice-friendly aerosols | kg-1 | 2 | real | kind_phys | inout | T | +!! | nwfa2d | tendency_of_water_friendly_surface_aerosols_at_surface | instantaneous fake surface aerosol source | kg-1 s-1 | 1 | real | kind_phys | in | T | +!! | tgrs | air_temperature_updated_by_physics | model layer mean temperature | K | 2 | real | kind_phys | inout | F | +!! | prsl | air_pressure | mean layer pressure | Pa | 2 | real | kind_phys | in | F | +!! | phii | geopotential_at_interface | geopotential at model layer interfaces | m2 s-2 | 2 | real | kind_phys | in | F | +!! | omega | omega | layer mean vertical velocity | Pa s-1 | 2 | real | kind_phys | in | F | +!! | dtp | time_step_for_physics | physics timestep | s | 0 | real | kind_phys | in | F | +!! | kdt | index_of_time_step | current forecast iteration | index | 0 | integer | | in | F | +!! | rain | lwe_thickness_of_precipitation_amount_on_dynamics_timestep | total rain at this time step | m | 1 | real | kind_phys | inout | F | +!! | rainst | lwe_thickness_of_stratiform_precipitation_amount | stratiform rainfall amount on physics timestep | m | 1 | real | kind_phys | inout | F | +!! | snow | lwe_thickness_of_snow_amount_on_dynamics_timestep | snow fall at this time step | m | 1 | real | kind_phys | inout | F | +!! | ice | lwe_thickness_of_ice_amount_on_dynamics_timestep | ice fall at this time step | m | 1 | real | kind_phys | inout | F | +!! | graupel | lwe_thickness_of_graupel_amount_on_dynamics_timestep | graupel fall at this time step | m | 1 | real | kind_phys | inout | F | +!! | sr | ratio_of_snowfall_to_rainfall | ratio of snowfall to large-scale rainfall | frac | 1 | real | kind_phys | out | F | +!! | islmsk | sea_land_ice_mask | sea/land/ice mask (=0/1/2) | flag | 1 | integer | | in | F | +!! | refl_10cm | radar_reflectivity_10cm | instantaneous refl_10cm | dBZ | 2 | real | kind_phys | out | F | +!! | do_radar_ref | flag_for_radar_reflectivity | flag for radar reflectivity | flag | 0 | logical | | in | F | +!! | re_cloud | mean_effective_radius_for_liquid_cloud | mean effective radius for liquid cloud | micron | 2 | real | kind_phys | out | T | +!! | re_ice | mean_effective_radius_for_ice_cloud | mean effective radius for ice cloud | micron | 2 | real | kind_phys | out | T | +!! | re_snow | mean_effective_radius_for_snow_flake | mean effective radius for snow flake | micron | 2 | real | kind_phys | out | T | +!! | mpicomm | mpi_comm | MPI communicator | index | 0 | integer | | in | F | +!! | mpirank | mpi_rank | current MPI-rank | index | 0 | integer | | in | F | +!! | mpiroot | mpi_root | master MPI-rank | index | 0 | integer | | in | F | +!! | errmsg | error_message | error message for error handling in CCPP | none | 0 | character | len=* | out | F | +!! | errflg | error_flag | error flag for error handling in CCPP | flag | 0 | integer | | out | F | +!! +#endif + subroutine mp_thompson_hrrr_run(ncol, nlev, con_g, con_rd, & + spechum, qc, qr, qi, qs, qg, ni, nr, & + is_aerosol_aware, nc, nwfa, nifa, nwfa2d, & + tgrs, prsl, phii, omega, dtp, kdt, & + rain, rainst, snow, ice, graupel, sr, & + islmsk, & + refl_10cm, do_radar_ref, & + re_cloud, re_ice, re_snow, & + mpicomm, mpirank, mpiroot, & + errmsg, errflg) + + implicit none + + ! Interface variables + + ! Dimensions and constants + integer, intent(in ) :: ncol + integer, intent(in ) :: nlev + real(kind_phys), intent(in ) :: con_g + real(kind_phys), intent(in ) :: con_rd + ! Hydrometeors + real(kind_phys), intent(inout) :: spechum(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: qc(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: qr(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: qi(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: qs(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: qg(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: ni(1:ncol,1:nlev) + real(kind_phys), intent(inout) :: nr(1:ncol,1:nlev) + ! Aerosols + logical, intent(in) :: is_aerosol_aware + real(kind_phys), optional, intent(inout) :: nc(1:ncol,1:nlev) + real(kind_phys), optional, intent(inout) :: nwfa(1:ncol,1:nlev) + real(kind_phys), optional, intent(inout) :: nifa(1:ncol,1:nlev) + real(kind_phys), optional, intent(in ) :: nwfa2d(1:ncol) + ! State variables and timestep information + real(kind_phys), intent(inout) :: tgrs(1:ncol,1:nlev) + real(kind_phys), intent(in ) :: prsl(1:ncol,1:nlev) + real(kind_phys), intent(in ) :: phii(1:ncol,1:nlev+1) + real(kind_phys), intent(in ) :: omega(1:ncol,1:nlev) + real(kind_phys), intent(in ) :: dtp + integer, intent(in ) :: kdt + ! Rain/snow/graupel fall amounts and fraction of frozen precip + real(kind_phys), intent(inout) :: rain(1:ncol) + real(kind_phys), intent(inout) :: rainst(1:ncol) + real(kind_phys), intent(inout) :: snow(1:ncol) + real(kind_phys), intent(inout) :: ice(1:ncol) + real(kind_phys), intent(inout) :: graupel(1:ncol) + real(kind_phys), intent( out) :: sr(1:ncol) + ! Sea/land/ice mask (currently not used) + integer, intent(in ) :: islmsk(1:ncol) + ! Radar reflectivity + real(kind_phys), intent( out) :: refl_10cm(1:ncol,1:nlev) + logical, optional, intent(in ) :: do_radar_ref + ! Cloud effective radii + real(kind_phys), optional, intent( out) :: re_cloud(1:ncol,1:nlev) + real(kind_phys), optional, intent( out) :: re_ice(1:ncol,1:nlev) + real(kind_phys), optional, intent( out) :: re_snow(1:ncol,1:nlev) + ! MPI information + integer, intent(in) :: mpicomm + integer, intent(in) :: mpirank + integer, intent(in) :: mpiroot + ! CCPP error handling + character(len=*), intent( out) :: errmsg + integer, intent( out) :: errflg + + ! Local variables + + ! Air density + real(kind_phys) :: rho(1:ncol,1:nlev) ! kg m-3 + ! Hydrometeors + real(kind_phys) :: qv_mp(1:ncol,1:nlev) ! kg kg-1 (dry mixing ratio) + real(kind_phys) :: qc_mp(1:ncol,1:nlev) ! kg kg-1 (dry mixing ratio) + real(kind_phys) :: qr_mp(1:ncol,1:nlev) ! kg kg-1 (dry mixing ratio) + real(kind_phys) :: qi_mp(1:ncol,1:nlev) ! kg kg-1 (dry mixing ratio) + real(kind_phys) :: qs_mp(1:ncol,1:nlev) ! kg kg-1 (dry mixing ratio) + real(kind_phys) :: qg_mp(1:ncol,1:nlev) ! kg kg-1 (dry mixing ratio) + ! Vertical velocity and level width + real(kind_phys) :: w(1:ncol,1:nlev) ! m s-1 + real(kind_phys) :: dz(1:ncol,1:nlev) ! m + ! Rain/snow/graupel fall amounts + real(kind_phys) :: rain_mp(1:ncol) ! mm, dummy, not used + real(kind_phys) :: snow_mp(1:ncol) ! mm, dummy, not used + real(kind_phys) :: ice_mp(1:ncol) ! mm, dummy, not used + real(kind_phys) :: graupel_mp(1:ncol) ! mm, dummy, not used + real(kind_phys) :: delta_rain_mp(1:ncol) ! mm + real(kind_phys) :: delta_snow_mp(1:ncol) ! mm + real(kind_phys) :: delta_ice_mp(1:ncol) ! mm + real(kind_phys) :: delta_graupel_mp(1:ncol) ! mm + ! Radar reflectivity + real(kind_phys) :: vt_dbz_wt(1:ncol,1:nlev) ! dummy for reflectivity-weighted terminal velocity, which is + ! in argument list of mp_gt_driver but is not calculated + logical :: diagflag ! must be true if do_radar_ref is true, not used otherwise + integer :: do_radar_ref_mp ! integer instead of logical do_radar_ref + ! Effective cloud radii + logical :: do_effective_radii + real(kind_phys) :: re_cloud_mp(1:ncol,1:nlev) ! m + real(kind_phys) :: re_ice_mp(1:ncol,1:nlev) ! m + real(kind_phys) :: re_snow_mp(1:ncol,1:nlev) ! m + integer :: has_reqc + integer :: has_reqi + integer :: has_reqs + ! Dimensions used in mp_gt_driver + integer :: ids,ide, jds,jde, kds,kde, & + ims,ime, jms,jme, kms,kme, & + its,ite, jts,jte, kts,kte + + ! Initialize the CCPP error handling variables + errmsg = '' + errflg = 0 + + ! Convert specific humidity/moist mixing ratios to dry mixing ratios + qv_mp = spechum/(1.0_kind_phys-spechum) + qc_mp = qc/(1.0_kind_phys-spechum) + qr_mp = qr/(1.0_kind_phys-spechum) + qi_mp = qi/(1.0_kind_phys-spechum) + qs_mp = qs/(1.0_kind_phys-spechum) + qg_mp = qg/(1.0_kind_phys-spechum) + + ! Density of air in kg m-3 + rho = prsl/(con_rd*tgrs) + + if (is_aerosol_aware .and. .not. (present(nc) .and. present(nwfa) .and. present(nifa) .and. present(nwfa2d))) then + write(errmsg,fmt='(*(a))') 'Logic error in mp_thompson_hrrr_run:', & + ' aerosol-aware microphysics require all of the', & + ' following optional arguments: nc, nwfa, nifa, nwfa2d' + errflg = 1 + return + end if + + ! Convert omega in Pa s-1 to vertical velocity w in m s-1 + w = -omega/(rho*con_g) + + ! Layer width in m from geopotential in m2 s-2 + dz = (phii(:,2:nlev+1) - phii(:,1:nlev)) / con_g + + ! Accumulated values inside Thompson scheme, not used; + ! only use delta and add to inout variables (different units) + rain_mp = 0 + snow_mp = 0 + ice_mp = 0 + graupel_mp = 0 + delta_rain_mp = 0 + delta_snow_mp = 0 + delta_ice_mp = 0 + delta_graupel_mp = 0 + + ! Flags for calculating radar reflectivity; diagflag is redundant + if (do_radar_ref) then + diagflag = .true. + do_radar_ref_mp = 1 + else + diagflag = .false. + do_radar_ref_mp = 0 + end if + + if (present(re_cloud) .and. present(re_ice) .and. present(re_snow)) then + do_effective_radii = .true. + has_reqc = 1 + has_reqi = 1 + has_reqs = 1 + else if (.not.present(re_cloud) .and. .not.present(re_ice) .and. .not.present(re_snow)) then + do_effective_radii = .false. + has_reqc = 0 + has_reqi = 0 + has_reqs = 0 + else + write(errmsg,fmt='(*(a))') 'Logic error in mp_thompson_hrrr_run:', & + ' all or none of the following optional', & + ' arguments are required: re_cloud, re_ice, re_snow' + errflg = 1 + return + end if + ! Initialize to zero, intent(out) variables + re_cloud_mp = 0 + re_ice_mp = 0 + re_snow_mp = 0 + + ! Set internal dimensions + ids = 1 + ims = 1 + its = 1 + ide = ncol + ime = ncol + ite = ncol + jds = 1 + jms = 1 + jts = 1 + jde = 1 + jme = 1 + jte = 1 + kds = 1 + kms = 1 + kts = 1 + kde = nlev + kme = nlev + kte = nlev + +#ifdef DEBUG_AEROSOLS + if (mpirank==mpiroot) then + write(0,*) "DH DEBUG: called mp_thompson_hrrr_run, is_aerosol_aware=", is_aerosol_aware, & + & ", do_effective_radii=", do_effective_radii, ", do_radar_ref=", do_radar_ref, & + & ", diagflag=", diagflag, ", do_radar_ref_mp=", do_radar_ref_mp + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: prsl min/mean/max =", & + & minval(prsl), sum(prsl)/real(size(prsl)), maxval(prsl) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: tgrs min/mean/max =", & + & minval(tgrs), sum(tgrs)/real(size(tgrs)), maxval(tgrs) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: rho min/mean/max =", & + & minval(rho), sum(rho)/real(size(rho)), maxval(rho) + if (is_aerosol_aware) then + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: nwfa2d min/mean/max =", & + & minval(nwfa2d), sum(nwfa2d)/real(size(nwfa2d)), maxval(nwfa2d) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: nwfa min/mean/max =", & + & minval(nwfa), sum(nwfa)/real(size(nwfa)), maxval(nwfa) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: nifa min/mean/max =", & + & minval(nifa), sum(nifa)/real(size(nifa)), maxval(nifa) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: nc min/mean/max =", & + & minval(nc), sum(nc)/real(size(nc)), maxval(nc) + end if + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: ni min/mean/max =", & + & minval(ni), sum(ni)/real(size(ni)), maxval(ni) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: nr min/mean/max =", & + & minval(nr), sum(nr)/real(size(nr)), maxval(nr) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: omega min/mean/max =", & + & minval(omega), sum(omega)/real(size(omega)), maxval(omega) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: w min/mean/max =", & + & minval(w), sum(w)/real(size(w)), maxval(w) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: re_cloud_mp min/mean/max =", & + & minval(re_cloud_mp), sum(re_cloud_mp)/real(size(re_cloud_mp)), maxval(re_cloud_mp) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: re_ice_mp min/mean/max =", & + & minval(re_ice_mp), sum(re_ice_mp)/real(size(re_ice_mp)), maxval(re_ice_mp) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: re_snow_mp min/mean/max =", & + & minval(re_snow_mp), sum(re_snow_mp)/real(size(re_snow_mp)), maxval(re_snow_mp) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: delta_rain_mp min/mean/max =", & + & minval(delta_rain_mp), sum(delta_rain_mp)/real(size(delta_rain_mp)), maxval(delta_rain_mp) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: delta_snow_mp min/mean/max =", & + & minval(delta_snow_mp), sum(delta_snow_mp)/real(size(delta_snow_mp)), maxval(delta_snow_mp) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: delta_ice_mp min/mean/max =", & + & minval(delta_ice_mp), sum(delta_ice_mp)/real(size(delta_ice_mp)), maxval(delta_ice_mp) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run before: delta_graupel_mp min/mean/max =", & + & minval(delta_graupel_mp), sum(delta_graupel_mp)/real(size(delta_graupel_mp)), maxval(delta_graupel_mp) + end if +#endif + + ! Call Thompson MP with or without aerosols + if (is_aerosol_aware) then + call mp_gt_driver(qv=qv_mp, qc=qc_mp, qr=qr_mp, qi=qi_mp, qs=qs_mp, qg=qg_mp, & + ni=ni, nr=nr, nc=nc, & + nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, itimestep=kdt, & + rainnc=rain_mp, rainncv=delta_rain_mp, & + snownc=snow_mp, snowncv=delta_snow_mp, & + icenc=ice_mp, icencv=delta_ice_mp, & + graupelnc=graupel_mp, graupelncv=delta_graupel_mp, sr=sr, & + refl_10cm=refl_10cm, vt_dbz_wt=vt_dbz_wt, & + diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & + re_cloud=re_cloud_mp, re_ice=re_ice_mp, re_snow=re_snow_mp, & + has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & + ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & + ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & + its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & + errmsg=errmsg, errflg=errflg) + + else + call mp_gt_driver(qv=qv_mp, qc=qc_mp, qr=qr_mp, qi=qi_mp, qs=qs_mp, qg=qg_mp, & + ni=ni, nr=nr, nc=nc, & + tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, itimestep=kdt, & + rainnc=rain_mp, rainncv=delta_rain_mp, & + snownc=snow_mp, snowncv=delta_snow_mp, & + icenc=ice_mp, icencv=delta_ice_mp, & + graupelnc=graupel_mp, graupelncv=delta_graupel_mp, sr=sr, & + refl_10cm=refl_10cm, vt_dbz_wt=vt_dbz_wt, & + diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & + re_cloud=re_cloud_mp, re_ice=re_ice_mp, re_snow=re_snow_mp, & + has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & + ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & + ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & + its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & + errmsg=errmsg, errflg=errflg) + end if + if (errflg/=0) return + + ! convert dry mixing ratios to specific humidity/moist mixing ratios + spechum = qv_mp/(1.0_kind_phys+qv_mp) + qc = qc_mp/(1.0_kind_phys+qv_mp) + qr = qr_mp/(1.0_kind_phys+qv_mp) + qi = qi_mp/(1.0_kind_phys+qv_mp) + qs = qs_mp/(1.0_kind_phys+qv_mp) + qg = qg_mp/(1.0_kind_phys+qv_mp) + +#ifdef DEBUG_AEROSOLS + if (mpirank==mpiroot) then + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: prsl min/mean/max =", & + & minval(prsl), sum(prsl)/real(size(prsl)), maxval(prsl) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: tgrs min/mean/max =", & + & minval(tgrs), sum(tgrs)/real(size(tgrs)), maxval(tgrs) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: rho min/mean/max =", & + & minval(rho), sum(rho)/real(size(rho)), maxval(rho) + if (is_aerosol_aware) then + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: nwfa2d min/mean/max =", & + & minval(nwfa2d), sum(nwfa2d)/real(size(nwfa2d)), maxval(nwfa2d) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: nwfa min/mean/max =", & + & minval(nwfa), sum(nwfa)/real(size(nwfa)), maxval(nwfa) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: nifa min/mean/max =", & + & minval(nifa), sum(nifa)/real(size(nifa)), maxval(nifa) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: nc min/mean/max =", & + & minval(nc), sum(nc)/real(size(nc)), maxval(nc) + end if + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: ni min/mean/max =", & + & minval(ni), sum(ni)/real(size(ni)), maxval(ni) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: nr min/mean/max =", & + & minval(nr), sum(nr)/real(size(nr)), maxval(nr) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: omega min/mean/max =", & + & minval(omega), sum(omega)/real(size(omega)), maxval(omega) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: w min/mean/max =", & + & minval(w), sum(w)/real(size(w)), maxval(w) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: re_cloud_mp min/mean/max =", & + & minval(re_cloud_mp), sum(re_cloud_mp)/real(size(re_cloud_mp)), maxval(re_cloud_mp) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: re_ice_mp min/mean/max =", & + & minval(re_ice_mp), sum(re_ice_mp)/real(size(re_ice_mp)), maxval(re_ice_mp) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: re_snow_mp min/mean/max =", & + & minval(re_snow_mp), sum(re_snow_mp)/real(size(re_snow_mp)), maxval(re_snow_mp) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: delta_rain_mp min/mean/max =", & + & minval(delta_rain_mp), sum(delta_rain_mp)/real(size(delta_rain_mp)), maxval(delta_rain_mp) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: delta_snow_mp min/mean/max =", & + & minval(delta_snow_mp), sum(delta_snow_mp)/real(size(delta_snow_mp)), maxval(delta_snow_mp) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: delta_ice_mp min/mean/max =", & + & minval(delta_ice_mp), sum(delta_ice_mp)/real(size(delta_ice_mp)), maxval(delta_ice_mp) + write(0,'(a,3e16.7)') "DH DEBUG mp thompson run after: delta_graupel_mp min/mean/max =", & + & minval(delta_graupel_mp), sum(delta_graupel_mp)/real(size(delta_graupel_mp)), maxval(delta_graupel_mp) + end if +#endif + + ! Convert rainfall from mm to m and add deltas to inout variables + rain = rain + delta_rain_mp/1000.0_kind_phys + rainst = rainst + delta_rain_mp/1000.0_kind_phys + snow = snow + delta_snow_mp/1000.0_kind_phys + ice = ice + delta_ice_mp/1000.0_kind_phys + graupel = graupel + delta_graupel_mp/1000.0_kind_phys + + if (do_effective_radii) then + ! Convert m to micron + re_cloud = re_cloud_mp*1.0E6_kind_phys + re_ice = re_ice_mp*1.0E6_kind_phys + re_snow = re_snow_mp*1.0E6_kind_phys + end if + + end subroutine mp_thompson_hrrr_run + + ! DH* do we need to deallocate stuff? which function to call in module_mp_thompson_hrrr.F90? + subroutine mp_thompson_hrrr_finalize() + end subroutine mp_thompson_hrrr_finalize + +end module mp_thompson_hrrr