!+---+-----------------------------------------------------------------+ !.. 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). !.. !.. 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: 30 Jan 2009 (WRF version) ! 01 Oct 2009 (CM1 version) ! 12 January 2011: fixed bug with mass- and energy-conserving ! equations (neweqts=1,2) in cm1r14 !+---+-----------------------------------------------------------------+ !wrft:model_layer:physics !+---+-----------------------------------------------------------------+ ! MODULE module_mp_thompson ! USE module_wrf_error ! USE module_utility, ONLY: WRFU_Clock, WRFU_Alarm ! USE module_domain, ONLY : HISTORY_ALARM, Is_alarm_tstep IMPLICIT NONE LOGICAL, PARAMETER, PRIVATE:: iiwarm = .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 = 400.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. ! setting to Continental by default: GHB, 091001 REAL, PARAMETER, PRIVATE:: Nt_c = 300.E6 !..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 = 5.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 = 125.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 !..Capacitance of sphere and plates/aggregates: D**3, D**2 REAL, PARAMETER, PRIVATE:: C_cube = 0.5 REAL, PARAMETER, PRIVATE:: C_sqrd = 0.3 !..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.95 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-8 REAL, PARAMETER, PRIVATE:: eps = 1.E-29 !..Constants in Cooper curve relation for cloud ice number. REAL, PARAMETER, PRIVATE:: TNO = 5.0 * 0.2 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:: cv = Cp - R REAL, PARAMETER, PRIVATE:: cvv = 1408.5 REAL, PARAMETER, PRIVATE:: cpl = 4190.0 REAL, PARAMETER, PRIVATE:: cpi = 2106.0 !..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_r1 = 37 INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55 INTEGER, PARAMETER, PRIVATE:: ntb_t = 9 INTEGER, PRIVATE:: nic2, nii2, nii3, nir2, nir3, nis2, nig2 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 !..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 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/) !..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. ! These are now allocatable, 20090612, JM DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:):: & tcg_racg, tmr_racg, tcr_gacr, tmg_gacr, & tnr_racg, tnr_gacr DOUBLE PRECISION, 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 DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:):: & tpi_qcfz, tni_qcfz DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:,:):: & tpi_qrfz, tpg_qrfz, tni_qrfz, tnr_qrfz DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:,:):: & tps_iaus, tni_iaus, tpi_ide REAL, ALLOCATABLE, DIMENSION(:,:):: t_Efrw REAL, ALLOCATABLE, DIMENSION(:,:):: t_Efsw DOUBLE PRECISION, ALLOCATABLE, DIMENSION(:, :, :):: tnr_rev !..Variables holding a bunch of exponents and gamma values (cloud water, !.. cloud ice, rain, snow, then graupel). REAL, DIMENSION(3), PRIVATE:: cce, ccg REAL, PRIVATE:: ocg1, ocg2 REAL, DIMENSION(6), PRIVATE:: cie, cig REAL, PRIVATE:: oig1, oig2, obmi REAL, DIMENSION(12), 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 CHARACTER*256:: mp_debug !..Various radar related variables DOUBLE PRECISION, PARAMETER, PRIVATE:: lamda_radar = 0.10 ! in meters DOUBLE PRECISION, PRIVATE:: K_w, PI5, lamda4 COMPLEX*16, PRIVATE:: m_w_0, m_i_0 DOUBLE PRECISION, DIMENSION(nbins+1), PRIVATE:: simpson DOUBLE PRECISION, DIMENSION(3), PARAMETER, PRIVATE:: basis = & (/1.d0/3.d0, 4.d0/3.d0, 1.d0/3.d0/) INTEGER, PARAMETER, PRIVATE:: slen = 20 CHARACTER(len=slen), PRIVATE:: & mixingrulestring_s, matrixstring_s, inclusionstring_s, & hoststring_s, hostmatrixstring_s, hostinclusionstring_s, & mixingrulestring_g, matrixstring_g, inclusionstring_g, & hoststring_g, hostmatrixstring_g, hostinclusionstring_g !+---+ !+---+-----------------------------------------------------------------+ !..END DECLARATIONS !+---+-----------------------------------------------------------------+ !+---+ ! CONTAINS SUBROUTINE thompson_init IMPLICIT NONE INTEGER :: i, j, k, m, n LOGICAL :: do_init ! !jm allocate the lookup tables ! ! use whether the first lookup table has been allocated to also determine whether ! to initialize them all. do_init = .FALSE. IF ( .NOT. ALLOCATED( tcg_racg ) ) THEN ALLOCATE( tcg_racg(ntb_g,ntb_r1,ntb_r) ) do_init = .TRUE. ENDIF IF ( .NOT. ALLOCATED( tmr_racg ) ) ALLOCATE( tmr_racg(ntb_g,ntb_r1,ntb_r) ) IF ( .NOT. ALLOCATED( tcr_gacr ) ) ALLOCATE( tcr_gacr(ntb_g,ntb_r1,ntb_r) ) IF ( .NOT. ALLOCATED( tmg_gacr ) ) ALLOCATE( tmg_gacr(ntb_g,ntb_r1,ntb_r) ) IF ( .NOT. ALLOCATED( tnr_racg ) ) ALLOCATE( tnr_racg(ntb_g,ntb_r1,ntb_r) ) IF ( .NOT. ALLOCATED( tnr_gacr ) ) ALLOCATE( tnr_gacr(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,45) ) IF ( .NOT. ALLOCATED( tni_qcfz ) ) ALLOCATE( tni_qcfz(ntb_c,45) ) IF ( .NOT. ALLOCATED( tpi_qrfz) ) ALLOCATE( tpi_qrfz(ntb_r,ntb_r1,45) ) IF ( .NOT. ALLOCATED( tpg_qrfz) ) ALLOCATE( tpg_qrfz(ntb_r,ntb_r1,45) ) IF ( .NOT. ALLOCATED( tni_qrfz) ) ALLOCATE( tni_qrfz(ntb_r,ntb_r1,45) ) IF ( .NOT. ALLOCATED( tnr_qrfz) ) ALLOCATE( tnr_qrfz(ntb_r,ntb_r1,45) ) 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 ( do_init ) THEN ! !jm end mods 20090608 ! !..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. 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. cce(1) = mu_c + 1. cce(2) = bm_r + mu_c + 1. cce(3) = bm_r + mu_c + 4. ccg(1) = WGAMMA(cce(1)) ccg(2) = WGAMMA(cce(2)) ccg(3) = WGAMMA(cce(3)) ocg1 = 1./ccg(1) ocg2 = 1./ccg(2) 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 + bv_i 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)) 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. do n = 1, 12 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 + bv_s + 2. cse(6) = bm_s + bv_s + 3. 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 + mu_s + bv_s + 2. 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 + mu_g + 3. 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))) !..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 lookup tables for most costly calculations. !+---+-----------------------------------------------------------------+ DO k = 1, ntb_r DO j = 1, ntb_r1 DO i = 1, ntb_g tcg_racg(i,j,k) = 0.0d0 tmr_racg(i,j,k) = 0.0d0 tcr_gacr(i,j,k) = 0.0d0 tmg_gacr(i,j,k) = 0.0d0 tnr_racg(i,j,k) = 0.0d0 tnr_gacr(i,j,k) = 0.0d0 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 k = 1, 45 DO j = 1, ntb_r1 DO i = 1, ntb_r tpi_qrfz(i,j,k) = 0.0d0 tni_qrfz(i,j,k) = 0.0d0 tpg_qrfz(i,j,k) = 0.0d0 tnr_qrfz(i,j,k) = 0.0d0 ENDDO ENDDO DO i = 1, ntb_c tpi_qcfz(i,k) = 0.0d0 tni_qcfz(i,k) = 0.0d0 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 ! CALL wrf_debug(150, 'CREATING MICROPHYSICS LOOKUP TABLES ... ') ! WRITE (wrf_err_message, '(a, f5.2, a, f5.2, a, f5.2, a, f5.2)') & ! ' using: mu_c=',mu_c,' mu_i=',mu_i,' mu_r=',mu_r,' mu_g=',mu_g ! CALL wrf_debug(150, wrf_err_message) !..Collision efficiency between rain/snow and cloud water. ! CALL wrf_debug(200, ' creating qc collision eff tables') CALL table_Efrw CALL table_Efsw !..Drop evaporation. ! CALL wrf_debug(200, ' creating rain evap table') ! CALL table_dropEvap !..Initialize various constants for computing radar reflectivity. CALL radar_init IF (.not. iiwarm) THEN !..Rain collecting graupel & graupel collecting rain. ! CALL wrf_debug(200, ' creating rain collecting graupel table') call qr_acr_qg !..Rain collecting snow & snow collecting rain. ! CALL wrf_debug(200, ' creating rain collecting snow table') call qr_acr_qs !..Cloud water and rain freezing (Bigg, 1953). ! CALL wrf_debug(200, ' creating freezing of water drops table') call freezeH2O !..Conversion of some ice mass into snow category. ! CALL wrf_debug(200, ' creating ice converting to snow table') CALL qi_aut_qs ENDIF ! CALL wrf_debug(150, ' ... DONE microphysical lookup tables') END IF ! do_init ! jm END SUBROUTINE thompson_init !+---+-----------------------------------------------------------------+ ! !+---+-----------------------------------------------------------------+ !..This is a wrapper routine designed to transfer values from 3D to 1D. !+---+-----------------------------------------------------------------+ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, nci, ncr, & th0, t3d, pii, p, dzz, dt_in, RAIN, & tevar,train, & ruh,rvh,rmh,rr,dbz3d,getdbz) implicit none include 'input.incl' include 'timestat.incl' !..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(ib:ie, jb:je, kb:ke), INTENT(INOUT):: & qv, qc, qr, qi, qs, qg, nci, ncr, t3d, dbz3d REAL, DIMENSION(ib:ie, jb:je, kb:ke), INTENT(IN):: & th0, pii, p, dzz, rmh, rr REAL, DIMENSION(ib:ie, jb:je, nrain), INTENT(INOUT):: RAIN ! REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: & ! SNOWNC, SNOWNCV, GRAUPELNC, GRAUPELNCV ! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & ! refl_10cm REAL, INTENT(IN):: dt_in real*8, intent(inout) :: tevar,train real, intent(in), dimension(ib:ie) :: ruh real, intent(in), dimension(jb:je) :: rvh logical, INTENT(IN):: getdbz ! INTEGER, INTENT(IN):: itimestep ! TYPE (WRFU_Clock):: grid_clock ! TYPE (WRFU_Alarm), POINTER:: grid_alarms(:) !..Local variables REAL, DIMENSION(nk):: & qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & nr1d, t1d, p1d, dz1d REAL, DIMENSION(nk):: dBZ !!! REAL, DIMENSION(ni,nj,nk):: refl_10cm ! 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 INTEGER:: i, j, k, n 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:: dBZ_tstep DOUBLE PRECISION, DIMENSION(1:nk):: prv_rev real :: tem1,tem2 real*8, dimension(nj) :: bud1,bud2 !+---+ tem1 = dx*dy*dz tem2 = dx*dy dt = dt_in !$omp parallel do default(shared) & !$omp private(j) do j=1,nj bud1(j)=0.0d0 bud2(j)=0.0d0 enddo dBZ_tstep = getdbz ! if ( Is_alarm_tstep(grid_clock, grid_alarms(HISTORY_ALARM)) ) then ! dBZ_tstep = .true. ! endif ! i_start = 1 ! j_start = 1 ! i_end = ni ! j_end = nj ! if ( (ie-ib+1).gt.4 .and. (je-jb+1).lt.4) then ! i_start = ib + 2 ! i_end = ie - 1 ! j_start = jb ! j_end = je ! elseif ( (ie-ib+1).lt.4 .and. (je-jb+1).gt.4) then ! i_start = ib ! i_end = ie ! j_start = jb + 2 ! j_end = je - 1 ! endif ! qc_max = 0. ! qr_max = 0. ! qs_max = 0. ! qi_max = 0. ! qg_max = 0 ! ni_max = 0. ! nr_max = 0. ! imax_qc = 0 ! imax_qr = 0 ! imax_qi = 0 ! imax_qs = 0 ! imax_qg = 0 ! imax_ni = 0 ! imax_nr = 0 ! jmax_qc = 0 ! jmax_qr = 0 ! jmax_qi = 0 ! jmax_qs = 0 ! jmax_qg = 0 ! jmax_ni = 0 ! jmax_nr = 0 ! kmax_qc = 0 ! kmax_qr = 0 ! kmax_qi = 0 ! kmax_qs = 0 ! kmax_qg = 0 ! kmax_ni = 0 ! kmax_nr = 0 ! do i = 1, 256 ! mp_debug(i:i) = char(0) ! enddo !$omp parallel do default(shared) & !$omp private(i,j,k,n,pptrain,pptsnow,pptgraul,pptice, & !$omp t1d,p1d,dz1d,qv1d,qc1d,qi1d,qr1d,qs1d,qg1d,ni1d,nr1d,prv_rev,dbz) j_loop: do j = 1,nj i_loop: do i = 1,ni pptrain = 0. pptsnow = 0. pptgraul = 0. pptice = 0. ! RAINNCV(i,j) = 0. ! SNOWNCV(i,j) = 0. ! GRAUPELNCV(i,j) = 0. ! SR(i,j) = 0. do k = 1, nk !!! t1d(k) = (th0(i,j,k)+th(i,j,k))*pii(i,j,k) t1d(k) = t3d(i,j,k) p1d(k) = p(i,j,k) dz1d(k) = dzz(i,j,k) qv1d(k) = qv(i,j,k) qc1d(k) = qc(i,j,k) qi1d(k) = qi(i,j,k) qr1d(k) = qr(i,j,k) qs1d(k) = qs(i,j,k) qg1d(k) = qg(i,j,k) ni1d(k) = nci(i,j,k) nr1d(k) = ncr(i,j,k) enddo prv_rev = 0.0 call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & nr1d, t1d, p1d, dz1d, & pptrain, pptsnow, pptgraul, pptice, & 1, nk, dt, i, j, prv_rev, neweqts ) do n = 1,nrain RAIN(i,j,n) = RAIN(i,j,n) + 0.1*( pptrain + pptsnow + pptgraul + pptice ) bud1(j) = bud1(j) + ( pptrain + pptsnow + pptgraul + pptice )*ruh(i)*rvh(j)*tem2 enddo do k=1,nk bud2(j) = bud2(j) + rr(i,j,k)*prv_rev(k)*ruh(i)*rvh(j)*rmh(i,j,k)*tem1 enddo ! 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 ! SNOWNCV(i,j) = pptsnow + pptice ! SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + pptice ! GRAUPELNCV(i,j) = pptgraul ! GRAUPELNC(i,j) = GRAUPELNC(i,j) + pptgraul ! SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12) do k = 1, nk qv(i,j,k) = qv1d(k) qc(i,j,k) = qc1d(k) qi(i,j,k) = qi1d(k) qr(i,j,k) = qr1d(k) qs(i,j,k) = qs1d(k) qg(i,j,k) = qg1d(k) nci(i,j,k) = ni1d(k) ncr(i,j,k) = nr1d(k) !!! th(i,j,k) = th(i,j,k) + (t1d(k)-tsave(k))/pii(i,j,k) t3d(i,j,k) = t1d(k) ! if (qc1d(k) .gt. qc_max) then ! imax_qc = i ! jmax_qc = j ! kmax_qc = k ! qc_max = qc1d(k) ! elseif (qc1d(k) .lt. 0.0) then ! write(mp_debug,*) 'WARNING, negative qc ', qc1d(k), & !! ' at i,j,k=', i,j,k !! CALL wrf_debug(150, mp_debug) ! endif ! if (qr1d(k) .gt. qr_max) then ! imax_qr = i ! jmax_qr = j ! kmax_qr = k ! qr_max = qr1d(k) ! elseif (qr1d(k) .lt. 0.0) then ! write(mp_debug,*) 'WARNING, negative qr ', qr1d(k), & ! ' at i,j,k=', i,j,k !! CALL wrf_debug(150, mp_debug) ! endif ! if (nr1d(k) .gt. nr_max) then ! imax_nr = i ! jmax_nr = j ! kmax_nr = k ! nr_max = nr1d(k) ! elseif (nr1d(k) .lt. 0.0) then ! write(mp_debug,*) 'WARNING, negative nr ', nr1d(k), & ! ' at i,j,k=', i,j,k !! CALL wrf_debug(150, mp_debug) ! endif ! if (qs1d(k) .gt. qs_max) then ! imax_qs = i ! jmax_qs = j ! kmax_qs = k ! qs_max = qs1d(k) ! elseif (qs1d(k) .lt. 0.0) then ! write(mp_debug,*) 'WARNING, negative qs ', qs1d(k), & ! ' at i,j,k=', i,j,k !! CALL wrf_debug(150, mp_debug) ! endif ! if (qi1d(k) .gt. qi_max) then ! imax_qi = i ! jmax_qi = j ! kmax_qi = k ! qi_max = qi1d(k) ! elseif (qi1d(k) .lt. 0.0) then ! write(mp_debug,*) 'WARNING, negative qi ', qi1d(k), & ! ' at i,j,k=', i,j,k !! CALL wrf_debug(150, mp_debug) ! endif ! if (qg1d(k) .gt. qg_max) then ! imax_qg = i ! jmax_qg = j ! kmax_qg = k ! qg_max = qg1d(k) ! elseif (qg1d(k) .lt. 0.0) then ! write(mp_debug,*) 'WARNING, negative qg ', qg1d(k), & ! ' at i,j,k=', i,j,k !! CALL wrf_debug(150, mp_debug) ! endif ! if (ni1d(k) .gt. ni_max) then ! imax_ni = i ! jmax_ni = j ! kmax_ni = k ! ni_max = ni1d(k) ! elseif (ni1d(k) .lt. 0.0) then ! write(mp_debug,*) 'WARNING, negative ni ', ni1d(k), & ! ' at i,j,k=', i,j,k !! CALL wrf_debug(150, mp_debug) ! endif ! if (qv1d(k) .lt. 0.0) then ! if (k.lt.nk-2 .and. k.gt.1) then ! qv(i,j,k) = 0.5*(qv(i,j,k-1) + qv(i,j,k+1)) ! else ! qv(i,j,k) = 1.E-7 ! endif ! write(mp_debug,*) 'WARNING, negative qv ', qv1d(k), & ! ' at i,j,k=', i,j,k !! CALL wrf_debug(150, mp_debug) ! endif enddo if (dBZ_tstep) then call calc_refl10cm (qv1d, qr1d, nr1d, qs1d, qg1d, t1d, p1d, & dBZ, 1, nk, i, j) do k = 1, nk dbz3d(i,j,k) = MAX(-35., dBZ(k)) enddo endif enddo i_loop enddo j_loop do j=1,nj train=train+bud1(j) tevar=tevar+bud2(j) enddo time_microphy=time_microphy+mytime() ! DEBUG - GT ! write(mp_debug,'(a,7(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', & ! 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', & ! 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', & ! 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', & ! 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', & ! 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', & ! 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', & ! 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')' ! write(*,*) mp_debug ! CALL wrf_debug(150, mp_debug) ! END DEBUG - GT ! do i = 1, 256 ! mp_debug(i:i) = char(0) ! enddo END SUBROUTINE mp_gt_driver !+---+-----------------------------------------------------------------+ ! !+---+-----------------------------------------------------------------+ !+---+-----------------------------------------------------------------+ !.. 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, t1d, p1d, dzq, & pptrain, pptsnow, pptgraul, pptice, & kts, kte, dt, ii, jj, prv_rev, neweqts ) implicit none include 'timestat.incl' !..Sub arguments INTEGER, INTENT(IN):: kts, kte, ii, jj REAL, DIMENSION(kts:kte), INTENT(INOUT):: & qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & nr1d, t1d, p1d REAL, DIMENSION(kts:kte), INTENT(IN):: dzq REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice REAL, INTENT(IN):: dt integer, intent(in) :: neweqts !..Local variables REAL, DIMENSION(kts:kte):: tten, qvten, qcten, qiten, & qrten, qsten, qgten, niten, nrten DOUBLE PRECISION, DIMENSION(kts:kte):: prw_vcd 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 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 REAL, DIMENSION(kts:kte):: temp, pres, qv REAL, DIMENSION(kts:kte):: rc, ri, rr, rs, rg, ni, nr REAL, DIMENSION(kts:kte):: rho, rhof, rhof2 REAL, DIMENSION(kts:kte):: qvs, qvsi REAL, DIMENSION(kts:kte):: satw, sati, ssatw, ssati REAL, DIMENSION(kts:kte):: diffu, visco, vsc2, & tcond, lvap, ocp, lvt2 REAL, DIMENSION(kts:kte):: econd,edep,efrz real :: cvm 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 REAL:: rgvm, delta_tp, orho, lfus2 REAL, DIMENSION(4):: onstep DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamc, lamr, lamg DOUBLE PRECISION:: lami, ilami REAL:: xDc, Dc_b, Dc_g, xDi, xDr, xDs, xDg, Ds_m, Dg_m DOUBLE PRECISION:: Dr_star REAL:: zeta1, zeta, taud, tau REAL:: stoke_r, stoke_s, stoke_g, stoke_i REAL:: vti, vtr, vts, vtg REAL, DIMENSION(kts:kte+1):: vtik, vtnik, vtrk, vtnrk, vtsk, vtgk 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 REAL:: dtsave, odts, odt, odzq INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq INTEGER, DIMENSION(4):: ksed1 INTEGER:: nir, nis, nig, nii, nic INTEGER:: idx_tc,idx_t,idx_s,idx_g,idx_r1,idx_r,idx_i1,idx_i,idx_c INTEGER:: idx, idx_d LOGICAL:: melti, no_micro LOGICAL, DIMENSION(kts:kte):: L_qc, L_qi, L_qr, L_qs, L_qg LOGICAL:: debug_flag !+---+ debug_flag = .false. ! if (ii.eq.280 .and. jj.eq.1) debug_flag = .true. no_micro = .true. dtsave = dt odt = 1./dt odts = 1./dtsave iexfrq = 1 !+---+-----------------------------------------------------------------+ !.. Source/sink terms. First 2 chars: "pr" represents source/sink of !.. mass while "pn" represents source/sink of number. Next char is one !.. of "v" for water vapor, "r" for rain, "i" for cloud ice, "w" for !.. cloud water, "s" for snow, and "g" for graupel. Next chars !.. represent processes: "de" for sublimation/deposition, "ev" for !.. evaporation, "fz" for freezing, "ml" for melting, "au" for !.. autoconversion, "nu" for ice nucleation, "hm" for Hallet/Mossop !.. secondary ice production, and "c" for collection followed by the !.. character for the species being collected. ALL of these terms are !.. positive (except for deposition/sublimation terms which can switch !.. signs based on super/subsaturation) and are treated as negatives !.. where necessary in the tendency equations. !+---+-----------------------------------------------------------------+ do k = kts, kte tten(k) = 0. qvten(k) = 0. qcten(k) = 0. qiten(k) = 0. qrten(k) = 0. qsten(k) = 0. qgten(k) = 0. niten(k) = 0. nrten(k) = 0. prw_vcd(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. 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. 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 (qc1d(k) .gt. R1) then no_micro = .false. rc(k) = qc1d(k)*rho(k) L_qc(k) = .true. else qc1d(k) = 0.0 rc(k) = R1 L_qc(k) = .false. endif if (qi1d(k) .gt. R1) then no_micro = .false. ri(k) = qi1d(k)*rho(k) ni(k) = MAX(1., ni1d(k)*rho(k)) 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. 30.E-6) then lami = cie(2)/30.E-6 ni(k) = MIN(250.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) = 0.01 L_qi(k) = .false. endif if (qr1d(k) .gt. R1) then no_micro = .false. rr(k) = qr1d(k)*rho(k) nr(k) = MAX(1., nr1d(k)*rho(k)) L_qr(k) = .true. if (nr(k) .gt. 1.0) then 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 if (qr1d(k) .gt. R2) then mvd_r(k) = 2.5E-3 else mvd_r(k) = 2.5E-3 / 3.0**(ALOG10(R2)-ALOG10(qr1d(k))) endif 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) = 1.0 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 !+---+-----------------------------------------------------------------+ !..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)) 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 if(neweqts.ge.1)then cvm = cv+cvv*qv1d(k)+cpl*(qc1d(k)+qr1d(k)) & +cpi*(qi1d(k)+qs1d(k)+qg1d(k)) econd(k) = (lvap(k)-Rv*temp(k))/cvm edep(k) = (lsub -Rv*temp(k))/cvm efrz(k) = (lsub-lvap(k))/cvm else econd(k) = lvap(k)*ocp(k) edep(k) = lsub*ocp(k) efrz(k) = (lsub-lvap(k))*ocp(k) endif 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 do k = kte, kts, -1 !-GT if (.not. L_qg(k)) CYCLE !-GT N0_exp = 100.0*rho(k)/rg(k) !-GT N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) N0_exp = (gonv_max-gonv_min)*0.5D0 & * tanh((0.15E-3-rg(k))/0.15E-3) & + (gonv_max+gonv_min)*0.5D0 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 - checked against my own !.. explicit/bin scheme and appears very good). RAIN2M if (L_qr(k) .and. mvd_r(k).gt. D0r) then pnr_rcr(k) = 8.*nr(k)*rr(k) endif if (.not. L_qc(k)) CYCLE xDc = MAX(D0c*1.E6, ((rc(k)/(am_r*Nt_c))**obmr) * 1.E6) lamc = (Nt_c*am_r* ccg(2) * ocg1 / rc(k))**obmr mvd_c(k) = (3.0+mu_c+0.672) / lamc !..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)*ocg2)**obmr / lamc) * 1.E6 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*mu_c/3.*D0r*D0r*D0r/rho(k)) ! RAIN2M endif !..Rain collecting cloud water. In CE, assume Dc<1). Either way, only bother to do sedimentation below !.. 1st level that contains any sedimenting particles (k=ksed1 on down). !.. New in v3.0+ is computing separate for rain, ice, snow, and !.. graupel species thus making code faster with credit to J. Schmidt. !+---+-----------------------------------------------------------------+ nstep = 0 onstep(1) = 1.0 ksed1(1) = 0 do k = kte+1, kts, -1 vtrk(k) = 0. vtnrk(k) = 0. vtik(k) = 0. vtnik(k) = 0. vtsk(k) = 0. vtgk(k) = 0. enddo do k = kte, kts, -1 vtr = 0. rhof(k) = SQRT(RHO_NOT/rho(k)) if (rr(k).gt. R2) then lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr vtr = rhof(k)*av_r*crg(6)*org3 * lamr**cre(3) & *((lamr+fv_r)**(-cre(6))) vtrk(k) = vtr ! First below is technically correct: ! vtr = rhof(k)*av_r*crg(5)*org2 * lamr**cre(2) & ! *((lamr+fv_r)**(-cre(5))) ! Test: make number fall faster (but still slower than mass) ! Goal: less prominent size sorting vtr = rhof(k)*av_r*crg(7)/crg(12) * lamr**cre(12) & *((lamr+fv_r)**(-cre(7))) vtnrk(k) = vtr endif if (MAX(vtrk(k),vtnrk(k)) .gt. 1.E-3) then ksed1(1) = MAX(ksed1(1), k) delta_tp = dzq(k)/(MAX(vtrk(k),vtnrk(k))) nstep = MAX(nstep, INT(DT/delta_tp + 1.)) endif enddo if (ksed1(1) .eq. kte) ksed1(1) = kte-1 if (nstep .gt. 0) onstep(1) = 1./REAL(nstep) !+---+-----------------------------------------------------------------+ if (.not. iiwarm) then nstep = 0 onstep(2) = 1.0 ksed1(2) = 0 do k = kte, kts, -1 vti = 0. if (ri(k).gt. R2) then lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi ilami = 1./lami vti = rhof(k)*av_i*cig(3)*oig2 * ilami**bv_i vtik(k) = vti vti = rhof(k)*av_i*cig(4)*oig1 * ilami**bv_i vtnik(k) = vti endif if (vtik(k) .gt. 1.E-3) then ksed1(2) = MAX(ksed1(2), k) delta_tp = dzq(k)/vtik(k) nstep = MAX(nstep, INT(DT/delta_tp + 1.)) endif enddo if (ksed1(2) .eq. kte) ksed1(2) = kte-1 if (nstep .gt. 0) onstep(2) = 1./REAL(nstep) !+---+-----------------------------------------------------------------+ nstep = 0 onstep(3) = 1.0 ksed1(3) = 0 do k = kte, kts, -1 vts = 0. if (rs(k).gt. R2) then xDs = smoc(k) / smob(k) Mrat = 1./xDs ils1 = 1./(Mrat*Lam0 + fv_s) ils2 = 1./(Mrat*Lam1 + fv_s) t1_vts = Kap0*csg(4)*ils1**cse(4) t2_vts = Kap1*Mrat**mu_s*csg(10)*ils2**cse(10) ils1 = 1./(Mrat*Lam0) ils2 = 1./(Mrat*Lam1) t3_vts = Kap0*csg(1)*ils1**cse(1) t4_vts = Kap1*Mrat**mu_s*csg(7)*ils2**cse(7) vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts) if (temp(k).gt. T_0) then vtsk(k) = MAX(vts*vts_boost(k), vtrk(k)) else vtsk(k) = vts*vts_boost(k) endif endif if (vtsk(k) .gt. 1.E-3) then ksed1(3) = MAX(ksed1(3), k) delta_tp = dzq(k)/vtsk(k) nstep = MAX(nstep, INT(DT/delta_tp + 1.)) endif enddo if (ksed1(3) .eq. kte) ksed1(3) = kte-1 if (nstep .gt. 0) onstep(3) = 1./REAL(nstep) !+---+-----------------------------------------------------------------+ nstep = 0 onstep(4) = 1.0 ksed1(4) = 0 do k = kte, kts, -1 vtg = 0. if (rg(k).gt. R2) then vtg = rhof(k)*av_g*cgg(6)*ogg3 * ilamg(k)**bv_g if (temp(k).gt. T_0) then vtgk(k) = MAX(vtg, vtrk(k)) else vtgk(k) = vtg endif endif if (vtgk(k) .gt. 1.E-3) then ksed1(4) = MAX(ksed1(4), k) delta_tp = dzq(k)/vtgk(k) nstep = MAX(nstep, INT(DT/delta_tp + 1.)) endif enddo if (ksed1(4) .eq. kte) ksed1(4) = kte-1 if (nstep .gt. 0) onstep(4) = 1./REAL(nstep) endif !+---+-----------------------------------------------------------------+ !..Sedimentation of mixing ratio is the integral of v(D)*m(D)*N(D)*dD, !.. whereas neglect m(D) term for number concentration. Therefore, !.. cloud ice has proper differential sedimentation. !.. New in v3.0+ is computing separate for rain, ice, snow, and !.. graupel species thus making code faster with credit to J. Schmidt. !+---+-----------------------------------------------------------------+ nstep = NINT(1./onstep(1)) do n = 1, nstep do k = kte, kts, -1 sed_r(k) = vtrk(k)*rr(k) sed_n(k) = vtnrk(k)*nr(k) enddo k = kte odzq = 1./dzq(k) orho = 1./rho(k) qrten(k) = qrten(k) - sed_r(k)*odzq*onstep(1)*orho nrten(k) = nrten(k) - sed_n(k)*odzq*onstep(1)*orho rr(k) = MAX(R1, rr(k) - sed_r(k)*odzq*DT*onstep(1)) nr(k) = MAX(1., nr(k) - sed_n(k)*odzq*DT*onstep(1)) do k = ksed1(1), kts, -1 odzq = 1./dzq(k) orho = 1./rho(k) qrten(k) = qrten(k) + (sed_r(k+1)-sed_r(k)) & *odzq*onstep(1)*orho nrten(k) = nrten(k) + (sed_n(k+1)-sed_n(k)) & *odzq*onstep(1)*orho rr(k) = MAX(R1, rr(k) + (sed_r(k+1)-sed_r(k)) & *odzq*DT*onstep(1)) nr(k) = MAX(1., nr(k) + (sed_n(k+1)-sed_n(k)) & *odzq*DT*onstep(1)) enddo pptrain = pptrain + sed_r(kts)*DT*onstep(1) enddo !+---+-----------------------------------------------------------------+ nstep = NINT(1./onstep(2)) do n = 1, nstep do k = kte, kts, -1 sed_i(k) = vtik(k)*ri(k) sed_n(k) = vtnik(k)*ni(k) enddo k = kte odzq = 1./dzq(k) orho = 1./rho(k) qiten(k) = qiten(k) - sed_i(k)*odzq*onstep(2)*orho niten(k) = niten(k) - sed_n(k)*odzq*onstep(2)*orho ri(k) = MAX(R1, ri(k) - sed_i(k)*odzq*DT*onstep(2)) ni(k) = MAX(1., ni(k) - sed_n(k)*odzq*DT*onstep(2)) do k = ksed1(2), kts, -1 odzq = 1./dzq(k) orho = 1./rho(k) qiten(k) = qiten(k) + (sed_i(k+1)-sed_i(k)) & *odzq*onstep(2)*orho niten(k) = niten(k) + (sed_n(k+1)-sed_n(k)) & *odzq*onstep(2)*orho ri(k) = MAX(R1, ri(k) + (sed_i(k+1)-sed_i(k)) & *odzq*DT*onstep(2)) ni(k) = MAX(1., ni(k) + (sed_n(k+1)-sed_n(k)) & *odzq*DT*onstep(2)) enddo pptice = pptice + sed_i(kts)*DT*onstep(2) enddo !+---+-----------------------------------------------------------------+ nstep = NINT(1./onstep(3)) do n = 1, nstep do k = kte, kts, -1 sed_s(k) = vtsk(k)*rs(k) enddo k = kte odzq = 1./dzq(k) orho = 1./rho(k) qsten(k) = qsten(k) - sed_s(k)*odzq*onstep(3)*orho rs(k) = MAX(R1, rs(k) - sed_s(k)*odzq*DT*onstep(3)) do k = ksed1(3), kts, -1 odzq = 1./dzq(k) orho = 1./rho(k) qsten(k) = qsten(k) + (sed_s(k+1)-sed_s(k)) & *odzq*onstep(3)*orho rs(k) = MAX(R1, rs(k) + (sed_s(k+1)-sed_s(k)) & *odzq*DT*onstep(3)) enddo pptsnow = pptsnow + sed_s(kts)*DT*onstep(3) enddo !+---+-----------------------------------------------------------------+ nstep = NINT(1./onstep(4)) do n = 1, nstep do k = kte, kts, -1 sed_g(k) = vtgk(k)*rg(k) enddo k = kte odzq = 1./dzq(k) orho = 1./rho(k) qgten(k) = qgten(k) - sed_g(k)*odzq*onstep(4)*orho rg(k) = MAX(R1, rg(k) - sed_g(k)*odzq*DT*onstep(4)) do k = ksed1(4), kts, -1 odzq = 1./dzq(k) orho = 1./rho(k) qgten(k) = qgten(k) + (sed_g(k+1)-sed_g(k)) & *odzq*onstep(4)*orho rg(k) = MAX(R1, rg(k) + (sed_g(k+1)-sed_g(k)) & *odzq*DT*onstep(4)) enddo pptgraul = pptgraul + sed_g(kts)*DT*onstep(4) enddo !!! time_fall=time_fall+mytime() !+---+-----------------------------------------------------------------+ !.. Instantly melt any cloud ice into cloud water if above 0C and !.. instantly freeze any cloud water found below HGFR. !+---+-----------------------------------------------------------------+ if (.not. iiwarm) then do k = kts, kte xri = MAX(0.0, qi1d(k) + qiten(k)*DT) if ( (temp(k).gt. T_0) .and. (xri.gt. 0.0) ) then qcten(k) = qcten(k) + xri*odt qiten(k) = -qi1d(k)*odt niten(k) = -ni1d(k)*odt !!! tten(k) = tten(k) - lfus*ocp(k)*xri*odt*(1-IFDRY) tten(k) = tten(k) - efrz(k)*xri*odt*(1-IFDRY) endif xrc = MAX(0.0, qc1d(k) + qcten(k)*DT) if ( (temp(k).lt. HGFR) .and. (xrc.gt. 0.0) ) then lfus2 = lsub - lvap(k) qiten(k) = qiten(k) + xrc*odt niten(k) = niten(k) + xrc/(2.*xm0i)*odt qcten(k) = -xrc*odt !!! tten(k) = tten(k) + lfus2*ocp(k)*xrc*odt*(1-IFDRY) tten(k) = tten(k) + efrz(k)*xrc*odt*(1-IFDRY) endif enddo endif !+---+-----------------------------------------------------------------+ !.. All tendencies computed, apply and pass back final values to parent. !+---+-----------------------------------------------------------------+ do k = kts, kte t1d(k) = t1d(k) + tten(k)*DT qv1d(k) = MAX(1.E-10, qv1d(k) + qvten(k)*DT) qc1d(k) = qc1d(k) + qcten(k)*DT if (qc1d(k) .le. R1) qc1d(k) = 0.0 qi1d(k) = qi1d(k) + qiten(k)*DT ni1d(k) = ni1d(k) + niten(k)*DT if (qi1d(k) .le. R1) then qi1d(k) = 0.0 ni1d(k) = 0.0 else if (ni1d(k) .gt. 1.0) then lami = (am_i*cig(2)*oig1*ni1d(k)/qi1d(k))**obmi ilami = 1./lami xDi = (bm_i + mu_i + 1.) * ilami if (xDi.lt. 30.E-6) then lami = cie(2)/30.E-6 elseif (xDi.gt. 300.E-6) then lami = cie(2)/300.E-6 endif else lami = cie(2)/D0s endif ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, & 250.D3/rho(k)) endif qr1d(k) = qr1d(k) + qrten(k)*DT nr1d(k) = nr1d(k) + nrten(k)*DT if (qr1d(k) .le. R1) then qr1d(k) = 0.0 nr1d(k) = 0.0 else if (nr1d(k) .gt. 1.0) then 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 else if (qr1d(k) .gt. R2) then mvd_r(k) = 2.5E-3 else mvd_r(k) = 2.5E-3 / 3.0**(ALOG10(R2)-ALOG10(qr1d(k))) endif 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 !!! time_microphy=time_microphy+mytime() end subroutine mp_thompson !+---+-----------------------------------------------------------------+ ! subroutine radar_init IMPLICIT NONE INTEGER:: n PI5 = PI*PI*PI*PI*PI 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, nbins+1 simpson(n) = 0.0d0 enddo do n = 1, nbins-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' end subroutine radar_init !+---+-----------------------------------------------------------------+ !..Compute radar reflectivity assuming 10 cm wavelength radar and using !.. Rayleigh approximation. Only complication is melted snow/graupel !.. which we treat as water-coated ice spheres and use Uli Blahak's !.. library of routines. The meltwater fraction is simply the amount !.. of frozen species remaining from what initially existed at the !.. melting level interface. !+---+-----------------------------------------------------------------+ subroutine calc_refl10cm (qv1d, qr1d, nr1d, qs1d, qg1d, t1d, p1d, & dBZ, kts, kte, ii, jj) IMPLICIT NONE !..Sub arguments INTEGER, INTENT(IN):: kts, kte, ii, jj REAL, DIMENSION(kts:kte), INTENT(IN):: & qv1d, qr1d, nr1d, qs1d, qg1d, t1d, p1d REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ !..Local variables REAL, DIMENSION(kts:kte):: temp, pres, qv, rho REAL, DIMENSION(kts:kte):: 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 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 !..Single melting snow/graupel particle 70% meltwater on external sfc DOUBLE PRECISION, PARAMETER:: melt_outside_s = 0.7d0 DOUBLE PRECISION, PARAMETER:: melt_outside_g = 0.7d0 DOUBLE PRECISION:: cback, x, eta, f_d !+---+ 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)) if (qr1d(k) .gt. R2) then rr(k) = qr1d(k)*rho(k) nr(k) = nr1d(k)*rho(k) lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr ilamr(k) = 1./lamr mvd_r(k) = (3.0 + mu_r + 0.672) * ilamr(k) 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) 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) endif ilamr(k) = 1./lamr N0_r(k) = nr(k)*org2*lamr**cre(2) L_qr(k) = .true. else rr(k) = R1 nr(k) = 1.0 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 do k = kte, kts, -1 !-GT N0_exp = 200.0*rho(k)/rg(k) !-GT N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) N0_exp = (gonv_max-gonv_min)*0.5D0 & * tanh((0.15E-3-rg(k))/0.15E-3) & + (gonv_max+gonv_min)*0.5D0 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. T_0) .and. (rr(k).gt. 0.001e-3) & .and. ((rs(k+1)+rg(k+1)).gt. 0.01e-3) ) 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/rho_i)*(am_s/rho_i)*smoz(k) if (L_qg(k)) ze_graupel(k) = (0.176/0.93) * (6.0/PI)*(6.0/PI) & * (am_g/rho_i)*(am_g/rho_i) & * 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, 1, -1 !..Reflectivity contributed by melting snow fmelt_s = MAX(0.05d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0)) if (fmelt_s.gt.0.01d0 .and. fmelt_s.lt.0.999d0 .and. & L_qs(k)) then 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, nbs x = am_s * Ds(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*Ds(n)) & + Kap1*(M0*Ds(n))**mu_s * DEXP(-slam2*Ds(n))) eta = eta + f_d * CBACK * simpson(n) * dts(n) enddo ze_snow(k) = SNGL(lamda4 / (pi5 * K_w) * eta) endif !..Reflectivity contributed by melting graupel fmelt_g = MAX(0.05d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0)) if (fmelt_g.gt.0.01d0 .and. fmelt_g.lt.0.999d0 .and. & L_qg(k)) then eta = 0.d0 lamg = 1./ilamg(k) do n = 1, nbg x = am_g * Dg(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)*Dg(n)**mu_g * DEXP(-lamg*Dg(n)) eta = eta + f_d * CBACK * simpson(n) * dtg(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 end subroutine calc_refl10cm !+---+-----------------------------------------------------------------+ ! !+---+-----------------------------------------------------------------+ !..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, n, n2 DOUBLE PRECISION, DIMENSION(nbg):: vg, N_g DOUBLE PRECISION, DIMENSION(nbr):: vr, N_r DOUBLE PRECISION:: N0_exp, N0_r, N0_g, lam_exp, lamg, lamr DOUBLE PRECISION:: massg, massr, dvg, dvr, t1, t2, z1, z2, y1, y2 !+---+ do n2 = 1, nbr vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2)) enddo do n = 1, nbg vg(n) = av_g*Dg(n)**bv_g enddo do k = 1, ntb_r do j = 1, ntb_r1 lam_exp = (N0r_exp(j)*am_r*crg(1)/r_r(k))**ore1 lamr = lam_exp * (crg(3)*org2*org1)**obmr N0_r = N0r_exp(j)/(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 i = 1, ntb_g !-GT N0_exp = 100.0d0/r_g(i) !-GT N0_exp = DMAX1(gonv_min*1.d0,DMIN1(N0_exp,gonv_max*1.d0)) N0_exp = (gonv_max-gonv_min)*0.5D0 & * tanh((0.15E-3-r_g(i))/0.15E-3) & + (gonv_max+gonv_min)*0.5D0 lam_exp = (N0_exp*am_g*cgg(1)/r_g(i))**oge1 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg N0_g = N0_exp/(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) = t1 tmr_racg(i,j,k) = DMIN1(z1, r_r(k)*1.0d0) tcr_gacr(i,j,k) = t2 tmg_gacr(i,j,k) = z2 tnr_racg(i,j,k) = y1 tnr_gacr(i,j,k) = y2 enddo enddo enddo end subroutine qr_acr_qg !+---+-----------------------------------------------------------------+ ! !+---+-----------------------------------------------------------------+ !..Rain collecting snow (and inverse). Explicit CE integration. !+---+-----------------------------------------------------------------+ subroutine qr_acr_qs implicit none !..Local variables INTEGER:: i, j, k, m, n, n2 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 !+---+ do n2 = 1, nbr vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*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 do m = 1, ntb_r do k = 1, ntb_r1 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. 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. 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 enddo end subroutine qr_acr_qs !+---+-----------------------------------------------------------------+ ! !+---+-----------------------------------------------------------------+ !..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, 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 !+---+ 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 k = 1, 45 ! print*, ' Freezing water for temp = ', -k Texp = DEXP( DFLOAT(k) ) - 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 = 1, nbr 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 enddo tpi_qrfz(i,j,k) = sum1 tni_qrfz(i,j,k) = sumn1 tpg_qrfz(i,j,k) = sum2 tnr_qrfz(i,j,k) = sumn2 enddo enddo do i = 1, ntb_c lamc = 1.0D-6 * (Nt_c*am_r* ccg(2) * ocg1 / r_c(i))**obmr N0_c = 1.0D-18 * Nt_c*ocg1 * lamc**cce(1) sum1 = 0.0d0 sumn2 = 0.0d0 do n = 1, nbc y = Dc(n)*1.0D6 vol = massc(n)*orho_w prob = 1.0D0 - DEXP(-120.0D0*vol*5.2D-4 * Texp) N_c(n) = N0_c* y**mu_c * EXP(-lamc*y)*dtc(n) N_c(n) = 1.0D24 * N_c(n) sumn2 = sumn2 + prob*N_c(n) sum1 = sum1 + prob*N_c(n)*massc(n) enddo tpi_qcfz(i,k) = sum1 tni_qcfz(i,k) = sumn2 enddo enddo end subroutine freezeH2O !+---+-----------------------------------------------------------------+ ! !+---+-----------------------------------------------------------------+ !..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 !+---+ 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 tpi_ide(i,j) = GAMMP(mu_i+2.0, SNGL(lami)*D0s) * 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 ! !+---+-----------------------------------------------------------------+ !..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 #ifdef DP t_Efrw(i,j) = MAX(0.0, MIN( (Ef_rw), 0.95)) #else t_Efrw(i,j) = MAX(0.0, MIN(SNGL(Ef_rw), 0.95)) #endif enddo enddo end subroutine table_Efrw ! !+---+-----------------------------------------------------------------+ !..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)) #ifdef DP t_Efsw(i,j) = MAX(0.0, MIN( (Ef_sw), 0.95)) #else t_Efsw(i,j) = MAX(0.0, MIN(SNGL(Ef_sw), 0.95)) #endif endif enddo enddo end subroutine table_Efsw ! !+---+-----------------------------------------------------------------+ !..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 DOUBLE PRECISION:: Nt_r, N0, lam_exp, lam INTEGER:: i, j, k 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 #ifdef DP tnr_rev(i,j,k) = GAMMP(mu_r+1.0, (Dr(i)*lam)) * Nt_r #else tnr_rev(i,j,k) = GAMMP(mu_r+1.0, SNGL(Dr(i)*lam)) * Nt_r #endif enddo enddo enddo end subroutine table_dropEvap ! 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)) ! ! !+---+-----------------------------------------------------------------+ !+---+-----------------------------------------------------------------+ 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 ! include 'constants.incl' REAL, INTENT(IN):: P, T REAL:: ESL,X REAL, PARAMETER:: C0= .611583699E03 REAL, PARAMETER:: C1= .444606896E02 REAL, PARAMETER:: C2= .143177157E01 REAL, PARAMETER:: C3= .264224321E-1 REAL, PARAMETER:: C4= .299291081E-3 REAL, PARAMETER:: C5= .203154182E-5 REAL, PARAMETER:: C6= .702620698E-8 REAL, PARAMETER:: C7= .379534310E-11 REAL, PARAMETER:: C8=-.321582393E-13 ! X=MAX(-80.,T-273.16) ! ESL=612.2*EXP(17.67*X/(T-29.65)) ! ESL=C0+X*(C1+X*(C2+X*(C3+X*(C4+X*(C5+X*(C6+X*(C7+X*C8))))))) ESL=611.2 * EXP( 17.67 * ( T - 273.15 ) / ( T - 29.65 ) ) RSLF=.622*ESL/(P-ESL) END FUNCTION RSLF ! ! 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. ! Psat = 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)) ! !+---+-----------------------------------------------------------------+ ! 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=611.2 * EXP( 21.8745584 * ( T - 273.15 ) / ( T - 7.66 ) ) RSIF=.622*ESI/(P-ESI) END FUNCTION RSIF ! ! 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) ! !+---+-----------------------------------------------------------------+ 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) 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*PI*0.5) & + (lambdas/lambda)**(2d0-2d0*alpha) epsr = epsinf + ((epss-epsinf) * ((lambdas/lambda)**(1d0-alpha) & * sin(alpha*PI*0.5)+1d0)) / nenner epsi = ((epss-epsinf) * ((lambdas/lambda)**(1d0-alpha) & * cos(alpha*PI*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 ! 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 = PI/6. * D_g**3 rhog = DMAX1(DMIN1(x_g / vg, DBLE(rho_i)), 10.0d0) vg = x_g / rhog meltratio_outside_grenz = 1.0d0 - rhog / rho_w 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=(rho_i-rhog)/(mra*rho_i-rhog+rho_i*rhog/rho_w) 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) / rho_i + x_w / rho_w endif endif D_large = (6.0 / PI * volg) ** (1./3.) volice = (x_g - x_w) / (volg * rho_i) volwater = x_w / (rho_w * 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(mp_debug,*) 'GET_M_MIX_NESTED: bad matrix: ', matrix ! CALL wrf_debug(150, mp_debug) 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(mp_debug,*) 'GET_M_MIX_NESTED: bad hostmatrix: ', & hostmatrix ! CALL wrf_debug(150, mp_debug) cumulerror = cumulerror + 1 endif endif elseif (host .eq. 'ice') then if (matrix .eq. 'ice') then write(mp_debug,*) 'GET_M_MIX_NESTED: bad matrix: ', matrix ! CALL wrf_debug(150, mp_debug) 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(mp_debug,*) 'GET_M_MIX_NESTED: bad hostmatrix: ', & hostmatrix ! CALL wrf_debug(150, mp_debug) cumulerror = cumulerror + 1 endif endif elseif (host .eq. 'water') then if (matrix .eq. 'water') then write(mp_debug,*) 'GET_M_MIX_NESTED: bad matrix: ', matrix ! CALL wrf_debug(150, mp_debug) 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.0d0*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.0d0*m_a, mtmp, m_w, & 0.0d0, (1.0d0-volwater), volwater, mixingrule, & 'ice', hostinclusion, error) cumulerror = cumulerror + error else write(mp_debug,*) 'GET_M_MIX_NESTED: bad hostmatrix: ', & hostmatrix ! CALL wrf_debug(150, mp_debug) 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(mp_debug,*) 'GET_M_MIX_NESTED: unknown matrix: ', host ! CALL wrf_debug(150, mp_debug) cumulerror = cumulerror + 1 endif IF (cumulerror .ne. 0) THEN write(mp_debug,*) 'GET_M_MIX_NESTED: error encountered' ! CALL wrf_debug(150, mp_debug) 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(mp_debug,*) 'GET_M_MIX: unknown matrix: ', matrix ! CALL wrf_debug(150, mp_debug) error = 1 endif else write(mp_debug,*) 'GET_M_MIX: unknown mixingrule: ', mixingrule ! CALL wrf_debug(150, mp_debug) error = 2 endif if (error .ne. 0) then write(mp_debug,*) 'GET_M_MIX: error encountered' ! CALL wrf_debug(150, mp_debug) 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(mp_debug,*) 'M_COMPLEX_MAXWELLGARNETT: sum of the ', & 'partial volume fractions is not 1...ERROR' ! CALL wrf_debug(150, mp_debug) 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(mp_debug,*) 'M_COMPLEX_MAXWELLGARNETT: ', & 'unknown inclusion: ', inclusion ! CALL wrf_debug(150, mp_debug) 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 !+---+-----------------------------------------------------------------+ END MODULE module_mp_thompson