/wrfv2_fire/phys/module_mp_thompson.F
FORTRAN Legacy | 3629 lines | 2525 code | 368 blank | 736 comment | 188 complexity | 3d748c85f72e2afd42fc1cc9f03b6927 MD5 | raw file
Possible License(s): AGPL-1.0
Large files files are truncated, but you can click here to view the full file
- !+---+-----------------------------------------------------------------+
- !.. 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: 20 Mar 2012
- !+---+-----------------------------------------------------------------+
- !wrft:model_layer:physics
- !+---+-----------------------------------------------------------------+
- !
- MODULE module_mp_thompson
- USE module_wrf_error
- ! USE module_mp_radar
- ! 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 = 500.0
- REAL, PARAMETER, PRIVATE:: rho_i = 890.0
- !..Prescribed number of cloud droplets. Set according to known data or
- !.. roughly 100 per cc (100.E6 m^-3) for Maritime cases and
- !.. 300 per cc (300.E6 m^-3) for Continental. Gamma shape parameter,
- !.. mu_c, calculated based on Nt_c is important in autoconversion
- !.. scheme.
- REAL, PARAMETER, PRIVATE:: Nt_c = 100.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 = 2.E4
- REAL, PARAMETER, PRIVATE:: gonv_max = 2.E6
- !..Mass power law relations: mass = am*D**bm
- !.. Snow from Field et al. (2005), others assume spherical form.
- REAL, PARAMETER, PRIVATE:: am_r = PI*rho_w/6.0
- REAL, PARAMETER, PRIVATE:: bm_r = 3.0
- REAL, PARAMETER, PRIVATE:: am_s = 0.069
- REAL, PARAMETER, PRIVATE:: bm_s = 2.0
- REAL, PARAMETER, PRIVATE:: am_g = PI*rho_g/6.0
- REAL, PARAMETER, PRIVATE:: bm_g = 3.0
- REAL, PARAMETER, PRIVATE:: am_i = PI*rho_i/6.0
- REAL, PARAMETER, PRIVATE:: bm_i = 3.0
- !..Fallspeed power laws relations: v = (av*D**bv)*exp(-fv*D)
- !.. Rain from Ferrier (1994), ice, snow, and graupel from
- !.. Thompson et al (2008). Coefficient fv is zero for graupel/ice.
- REAL, PARAMETER, PRIVATE:: av_r = 4854.0
- REAL, PARAMETER, PRIVATE:: bv_r = 1.0
- REAL, PARAMETER, PRIVATE:: fv_r = 195.0
- REAL, PARAMETER, PRIVATE:: av_s = 40.0
- REAL, PARAMETER, PRIVATE:: bv_s = 0.55
- REAL, PARAMETER, PRIVATE:: fv_s = 100.0
- REAL, PARAMETER, PRIVATE:: av_g = 442.0
- REAL, PARAMETER, PRIVATE:: bv_g = 0.89
- REAL, PARAMETER, PRIVATE:: av_i = 1847.5
- REAL, PARAMETER, PRIVATE:: bv_i = 1.0
- !..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.75
- REAL, PARAMETER, PRIVATE:: Ef_ri = 0.95
- !..Minimum microphys values
- !.. R1 value, 1.E-12, cannot be set lower because of numerical
- !.. problems with Paul Field's moments and should not be set larger
- !.. because of truncation problems in snow/ice growth.
- REAL, PARAMETER, PRIVATE:: R1 = 1.E-12
- REAL, PARAMETER, PRIVATE:: R2 = 1.E-6
- REAL, PARAMETER, PRIVATE:: eps = 1.E-15
- !..Constants in Cooper curve relation for cloud ice number.
- REAL, PARAMETER, PRIVATE:: TNO = 5.0
- REAL, PARAMETER, PRIVATE:: ATO = 0.304
- !..Rho_not used in fallspeed relations (rho_not/rho)**.5 adjustment.
- REAL, PARAMETER, PRIVATE:: rho_not = 101325.0/(287.05*298.0)
- !..Schmidt number
- REAL, PARAMETER, PRIVATE:: Sc = 0.632
- REAL, PRIVATE:: Sc3
- !..Homogeneous freezing temperature
- REAL, PARAMETER, PRIVATE:: HGFR = 235.16
- !..Water vapor and air gas constants at constant pressure
- REAL, PARAMETER, PRIVATE:: Rv = 461.5
- REAL, PARAMETER, PRIVATE:: oRv = 1./Rv
- REAL, PARAMETER, PRIVATE:: R = 287.04
- REAL, PARAMETER, PRIVATE:: Cp = 1004.0
- !..Enthalpy of sublimation, vaporization, and fusion at 0C.
- REAL, PARAMETER, PRIVATE:: lsub = 2.834E6
- REAL, PARAMETER, PRIVATE:: lvap0 = 2.5E6
- REAL, PARAMETER, PRIVATE:: lfus = lsub - lvap0
- REAL, PARAMETER, PRIVATE:: olfus = 1./lfus
- !..Ice initiates with this mass (kg), corresponding diameter calc.
- !..Min diameters and mass of cloud, rain, snow, and graupel (m, kg).
- REAL, PARAMETER, PRIVATE:: xm0i = 1.E-12
- REAL, PARAMETER, PRIVATE:: D0c = 1.E-6
- REAL, PARAMETER, PRIVATE:: D0r = 50.E-6
- REAL, PARAMETER, PRIVATE:: D0s = 200.E-6
- REAL, PARAMETER, PRIVATE:: D0g = 250.E-6
- REAL, PRIVATE:: D0i, xm0s, xm0g
- !..Lookup table dimensions
- INTEGER, PARAMETER, PRIVATE:: nbins = 100
- INTEGER, PARAMETER, PRIVATE:: nbc = nbins
- INTEGER, PARAMETER, PRIVATE:: nbi = nbins
- INTEGER, PARAMETER, PRIVATE:: nbr = nbins
- INTEGER, PARAMETER, PRIVATE:: nbs = nbins
- INTEGER, PARAMETER, PRIVATE:: nbg = nbins
- INTEGER, PARAMETER, PRIVATE:: ntb_c = 37
- INTEGER, PARAMETER, PRIVATE:: ntb_i = 64
- INTEGER, PARAMETER, PRIVATE:: ntb_r = 37
- INTEGER, PARAMETER, PRIVATE:: ntb_s = 28
- INTEGER, PARAMETER, PRIVATE:: ntb_g = 28
- INTEGER, PARAMETER, PRIVATE:: ntb_g1 = 28
- INTEGER, PARAMETER, PRIVATE:: ntb_r1 = 37
- INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55
- INTEGER, PARAMETER, PRIVATE:: ntb_t = 9
- INTEGER, PRIVATE:: nic2, nii2, nii3, nir2, nir3, nis2, nig2, nig3
- 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 graupel y-intercept parameter (/m**4).
- REAL, DIMENSION(ntb_g1), PARAMETER, PRIVATE:: &
- N0g_exp = (/1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
- 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
- 1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, &
- 1.e7/)
- !..Lookup tables for ice number concentration (/m**3).
- REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: &
- Nt_i = (/1.0,2.0,3.0,4.0,5.0,6.0,7.0,8.0,9.0, &
- 1.e1,2.e1,3.e1,4.e1,5.e1,6.e1,7.e1,8.e1,9.e1, &
- 1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, &
- 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, &
- 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
- 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
- 1.e6/)
- !..For snow moments conversions (from Field et al. 2005)
- REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
- sa = (/ 5.065339, -0.062659, -3.032362, 0.029469, -0.000285, &
- 0.31255, 0.000204, 0.003199, 0.0, -0.015952/)
- REAL, DIMENSION(10), PARAMETER, PRIVATE:: &
- sb = (/ 0.476221, -0.015896, 0.165977, 0.007468, -0.000141, &
- 0.060366, 0.000079, 0.000594, 0.0, -0.003577/)
- !..Temperatures (5 C interval 0 to -40) used in lookup tables.
- REAL, DIMENSION(ntb_t), PARAMETER, PRIVATE:: &
- Tc = (/-0.01, -5., -10., -15., -20., -25., -30., -35., -40./)
- !..Lookup tables for various accretion/collection terms.
- !.. ntb_x refers to the number of elements for rain, snow, graupel,
- !.. and temperature array indices. Variables beginning with t-p/c/m/n
- !.. represent lookup tables. Save compile-time memory by making
- !.. allocatable (2009Jun12, J. Michalakes).
- INTEGER, PARAMETER, PRIVATE:: R8SIZE = 8
- REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
- tcg_racg, tmr_racg, tcr_gacr, tmg_gacr, &
- tnr_racg, tnr_gacr
- REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: &
- tcs_racs1, tmr_racs1, tcs_racs2, tmr_racs2, &
- tcr_sacr1, tms_sacr1, tcr_sacr2, tms_sacr2, &
- tnr_racs1, tnr_racs2, tnr_sacr1, tnr_sacr2
- REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: &
- tpi_qcfz, tni_qcfz
- REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: &
- tpi_qrfz, tpg_qrfz, tni_qrfz, tnr_qrfz
- REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: &
- tps_iaus, tni_iaus, tpi_ide
- REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efrw
- REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:):: t_Efsw
- REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:):: tnr_rev
- !..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(7), PRIVATE:: cie, cig
- REAL, PRIVATE:: oig1, oig2, obmi
- REAL, DIMENSION(13), PRIVATE:: cre, crg
- REAL, PRIVATE:: ore1, org1, org2, org3, obmr
- REAL, DIMENSION(18), PRIVATE:: cse, csg
- REAL, PRIVATE:: oams, obms, ocms
- REAL, DIMENSION(12), PRIVATE:: cge, cgg
- REAL, PRIVATE:: oge1, ogg1, ogg2, ogg3, oamg, obmg, ocmg
- !..Declaration of precomputed constants in various rate eqns.
- REAL:: t1_qr_qc, t1_qr_qi, t2_qr_qi, t1_qg_qc, t1_qs_qc, t1_qs_qi
- REAL:: t1_qr_ev, t2_qr_ev
- REAL:: t1_qs_sd, t2_qs_sd, t1_qg_sd, t2_qg_sd
- REAL:: t1_qs_me, t2_qs_me, t1_qg_me, t2_qg_me
- CHARACTER*256:: mp_debug
- !+---+
- !+---+-----------------------------------------------------------------+
- !..END DECLARATIONS
- !+---+-----------------------------------------------------------------+
- !+---+
- !ctrlL
- CONTAINS
- SUBROUTINE thompson_init
- IMPLICIT NONE
- INTEGER:: i, j, k, m, n
- LOGICAL:: micro_init
- !..Allocate space for lookup tables (J. Michalakes 2009Jun08).
- micro_init = .FALSE.
- if (.NOT. ALLOCATED(tcg_racg) ) then
- ALLOCATE(tcg_racg(ntb_g1,ntb_g,ntb_r1,ntb_r))
- micro_init = .TRUE.
- endif
- if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tcr_gacr)) ALLOCATE(tcr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tmg_gacr)) ALLOCATE(tmg_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tnr_racg)) ALLOCATE(tnr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tnr_gacr)) ALLOCATE(tnr_gacr(ntb_g1,ntb_g,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tcs_racs1)) ALLOCATE(tcs_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tmr_racs1)) ALLOCATE(tmr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tcs_racs2)) ALLOCATE(tcs_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tmr_racs2)) ALLOCATE(tmr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tcr_sacr1)) ALLOCATE(tcr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tms_sacr1)) ALLOCATE(tms_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tcr_sacr2)) ALLOCATE(tcr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tms_sacr2)) ALLOCATE(tms_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tnr_racs1)) ALLOCATE(tnr_racs1(ntb_s,ntb_t,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tnr_racs2)) ALLOCATE(tnr_racs2(ntb_s,ntb_t,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tnr_sacr1)) ALLOCATE(tnr_sacr1(ntb_s,ntb_t,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tnr_sacr2)) ALLOCATE(tnr_sacr2(ntb_s,ntb_t,ntb_r1,ntb_r))
- if (.NOT. ALLOCATED(tpi_qcfz)) ALLOCATE(tpi_qcfz(ntb_c,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 (micro_init) then
- !..From Martin et al. (1994), assign gamma shape parameter mu for cloud
- !.. drops according to general dispersion characteristics (disp=~0.25
- !.. for Maritime and 0.45 for Continental).
- !.. disp=SQRT((mu+2)/(mu+1) - 1) so mu varies from 15 for Maritime
- !.. to 2 for really dirty air.
- 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*0.5 + mu_i + bv_i + 1.
- cie(7) = bm_i*0.5 + mu_i + 1.
- cig(1) = WGAMMA(cie(1))
- cig(2) = WGAMMA(cie(2))
- cig(3) = WGAMMA(cie(3))
- cig(4) = WGAMMA(cie(4))
- cig(5) = WGAMMA(cie(5))
- cig(6) = WGAMMA(cie(6))
- cig(7) = WGAMMA(cie(7))
- oig1 = 1./cig(1)
- oig2 = 1./cig(2)
- obmi = 1./bm_i
- cre(1) = bm_r + 1.
- cre(2) = mu_r + 1.
- cre(3) = bm_r + mu_r + 1.
- cre(4) = bm_r*2. + mu_r + 1.
- cre(5) = mu_r + bv_r + 1.
- cre(6) = bm_r + mu_r + bv_r + 1.
- cre(7) = bm_r*0.5 + mu_r + bv_r + 1.
- cre(8) = bm_r + mu_r + bv_r + 3.
- cre(9) = mu_r + bv_r + 3.
- cre(10) = mu_r + 2.
- cre(11) = 0.5*(bv_r + 5. + 2.*mu_r)
- cre(12) = bm_r*0.5 + mu_r + 1.
- cre(13) = bm_r*2. + mu_r + bv_r + 1.
- do n = 1, 13
- crg(n) = WGAMMA(cre(n))
- enddo
- obmr = 1./bm_r
- ore1 = 1./cre(1)
- org1 = 1./crg(1)
- org2 = 1./crg(2)
- org3 = 1./crg(3)
- cse(1) = bm_s + 1.
- cse(2) = bm_s + 2.
- cse(3) = bm_s*2.
- cse(4) = bm_s + bv_s + 1.
- cse(5) = bm_s*2. + bv_s + 1.
- cse(6) = bm_s*2. + 1.
- cse(7) = bm_s + mu_s + 1.
- cse(8) = bm_s + mu_s + 2.
- cse(9) = bm_s + mu_s + 3.
- cse(10) = bm_s + mu_s + bv_s + 1.
- cse(11) = bm_s*2. + mu_s + bv_s + 1.
- cse(12) = bm_s*2. + mu_s + 1.
- cse(13) = bv_s + 2.
- cse(14) = bm_s + bv_s
- cse(15) = mu_s + 1.
- cse(16) = 1.0 + (1.0 + bv_s)/2.
- cse(17) = cse(16) + mu_s + 1.
- cse(18) = bv_s + mu_s + 3.
- do n = 1, 18
- csg(n) = WGAMMA(cse(n))
- enddo
- oams = 1./am_s
- obms = 1./bm_s
- ocms = oams**obms
- cge(1) = bm_g + 1.
- cge(2) = mu_g + 1.
- cge(3) = bm_g + mu_g + 1.
- cge(4) = bm_g*2. + mu_g + 1.
- cge(5) = bm_g*2. + mu_g + bv_g + 1.
- cge(6) = bm_g + mu_g + bv_g + 1.
- cge(7) = bm_g + mu_g + bv_g + 2.
- cge(8) = bm_g + mu_g + bv_g + 3.
- cge(9) = mu_g + bv_g + 3.
- cge(10) = mu_g + 2.
- cge(11) = 0.5*(bv_g + 5. + 2.*mu_g)
- cge(12) = 0.5*(bv_g + 5.) + mu_g
- do n = 1, 12
- cgg(n) = WGAMMA(cge(n))
- enddo
- oamg = 1./am_g
- obmg = 1./bm_g
- ocmg = oamg**obmg
- oge1 = 1./cge(1)
- ogg1 = 1./cgg(1)
- ogg2 = 1./cgg(2)
- ogg3 = 1./cgg(3)
- !+---+-----------------------------------------------------------------+
- !..Simplify various rate eqns the best we can now.
- !+---+-----------------------------------------------------------------+
- !..Rain collecting cloud water and cloud ice
- t1_qr_qc = PI*.25*av_r * crg(9)
- t1_qr_qi = PI*.25*av_r * crg(9)
- t2_qr_qi = PI*.25*am_r*av_r * crg(8)
- !..Graupel collecting cloud water
- t1_qg_qc = PI*.25*av_g * cgg(9)
- !..Snow collecting cloud water
- t1_qs_qc = PI*.25*av_s
- !..Snow collecting cloud ice
- t1_qs_qi = PI*.25*av_s
- !..Evaporation of rain; ignore depositional growth of rain.
- t1_qr_ev = 0.78 * crg(10)
- t2_qr_ev = 0.308*Sc3*SQRT(av_r) * crg(11)
- !..Sublimation/depositional growth of snow
- t1_qs_sd = 0.86
- t2_qs_sd = 0.28*Sc3*SQRT(av_s)
- !..Melting of snow
- t1_qs_me = PI*4.*C_sqrd*olfus * 0.86
- t2_qs_me = PI*4.*C_sqrd*olfus * 0.28*Sc3*SQRT(av_s)
- !..Sublimation/depositional growth of graupel
- t1_qg_sd = 0.86 * cgg(10)
- t2_qg_sd = 0.28*Sc3*SQRT(av_g) * cgg(11)
- !..Melting of graupel
- t1_qg_me = PI*4.*C_cube*olfus * 0.86 * cgg(10)
- t2_qg_me = PI*4.*C_cube*olfus * 0.28*Sc3*SQRT(av_g) * cgg(11)
- !..Constants for helping find lookup table indexes.
- nic2 = NINT(ALOG10(r_c(1)))
- nii2 = NINT(ALOG10(r_i(1)))
- nii3 = NINT(ALOG10(Nt_i(1)))
- nir2 = NINT(ALOG10(r_r(1)))
- nir3 = NINT(ALOG10(N0r_exp(1)))
- nis2 = NINT(ALOG10(r_s(1)))
- nig2 = NINT(ALOG10(r_g(1)))
- nig3 = NINT(ALOG10(N0g_exp(1)))
- !..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 m = 1, ntb_r
- do k = 1, ntb_r1
- do j = 1, ntb_g
- do i = 1, ntb_g1
- tcg_racg(i,j,k,m) = 0.0d0
- tmr_racg(i,j,k,m) = 0.0d0
- tcr_gacr(i,j,k,m) = 0.0d0
- tmg_gacr(i,j,k,m) = 0.0d0
- tnr_racg(i,j,k,m) = 0.0d0
- tnr_gacr(i,j,k,m) = 0.0d0
- enddo
- enddo
- enddo
- enddo
- do m = 1, ntb_r
- do k = 1, ntb_r1
- do j = 1, ntb_t
- do i = 1, ntb_s
- tcs_racs1(i,j,k,m) = 0.0d0
- tmr_racs1(i,j,k,m) = 0.0d0
- tcs_racs2(i,j,k,m) = 0.0d0
- tmr_racs2(i,j,k,m) = 0.0d0
- tcr_sacr1(i,j,k,m) = 0.0d0
- tms_sacr1(i,j,k,m) = 0.0d0
- tcr_sacr2(i,j,k,m) = 0.0d0
- tms_sacr2(i,j,k,m) = 0.0d0
- tnr_racs1(i,j,k,m) = 0.0d0
- tnr_racs2(i,j,k,m) = 0.0d0
- tnr_sacr1(i,j,k,m) = 0.0d0
- tnr_sacr2(i,j,k,m) = 0.0d0
- enddo
- enddo
- enddo
- enddo
- do 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.
- ! xam_r = am_r
- ! xbm_r = bm_r
- ! xmu_r = mu_r
- ! xam_s = am_s
- ! xbm_s = bm_s
- ! xmu_s = mu_s
- ! xam_g = am_g
- ! xbm_g = bm_g
- ! xmu_g = mu_g
- ! call radar_init
- 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')
- endif
- END SUBROUTINE thompson_init
- !+---+-----------------------------------------------------------------+
- !ctrlL
- !+---+-----------------------------------------------------------------+
- !..This is a wrapper routine designed to transfer values from 3D to 1D.
- !+---+-----------------------------------------------------------------+
- SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, &
- th, pii, p, dz, dt_in, itimestep, &
- RAINNC, RAINNCV, &
- SNOWNC, SNOWNCV, &
- GRAUPELNC, GRAUPELNCV, SR, &
- #ifdef WRF_CHEM
- rainprod, evapprod, &
- #endif
- ! refl_10cm, grid_clock, grid_alarms, &
- ids,ide, jds,jde, kds,kde, & ! domain dims
- ims,ime, jms,jme, kms,kme, & ! memory dims
- its,ite, jts,jte, kts,kte) ! tile dims
- implicit none
- !..Subroutine arguments
- INTEGER, INTENT(IN):: ids,ide, jds,jde, kds,kde, &
- ims,ime, jms,jme, kms,kme, &
- its,ite, jts,jte, kts,kte
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
- qv, qc, qr, qi, qs, qg, ni, nr, th
- #ifdef WRF_CHEM
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
- rainprod, evapprod
- #endif
- REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(IN):: &
- pii, p, dz
- REAL, DIMENSION(ims:ime, jms:jme), INTENT(INOUT):: &
- RAINNC, RAINNCV, SR
- REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(INOUT):: &
- SNOWNC, SNOWNCV, GRAUPELNC, GRAUPELNCV
- ! REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: &
- ! refl_10cm
- REAL, INTENT(IN):: dt_in
- INTEGER, INTENT(IN):: itimestep
- ! TYPE (WRFU_Clock):: grid_clock
- ! TYPE (WRFU_Alarm), POINTER:: grid_alarms(:)
- !..Local variables
- REAL, DIMENSION(kts:kte):: &
- qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
- nr1d, t1d, p1d, dz1d, dBZ
- #ifdef WRF_CHEM
- REAL, DIMENSION(kts:kte):: &
- rainprod1d, evapprod1d
- #endif
- REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic
- REAL:: dt, pptrain, pptsnow, pptgraul, pptice
- REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max
- INTEGER:: i, j, k
- INTEGER:: imax_qc,imax_qr,imax_qi,imax_qs,imax_qg,imax_ni,imax_nr
- INTEGER:: jmax_qc,jmax_qr,jmax_qi,jmax_qs,jmax_qg,jmax_ni,jmax_nr
- INTEGER:: kmax_qc,kmax_qr,kmax_qi,kmax_qs,kmax_qg,kmax_ni,kmax_nr
- INTEGER:: i_start, j_start, i_end, j_end
- LOGICAL:: dBZ_tstep
- !+---+
- dBZ_tstep = .false.
- ! if ( Is_alarm_tstep(grid_clock, grid_alarms(HISTORY_ALARM)) ) then
- ! dBZ_tstep = .true.
- ! endif
- i_start = its
- j_start = jts
- i_end = ite
- j_end = jte
- !..For idealized testing by developer.
- ! if ( (ide-ids+1).gt.4 .and. (jde-jds+1).lt.4 .and. &
- ! ids.eq.its.and.ide.eq.ite.and.jds.eq.jts.and.jde.eq.jte) then
- ! i_start = its + 2
- ! i_end = ite
- ! j_start = jts
- ! j_end = jte
- ! endif
- dt = dt_in
-
- qc_max = 0.
- qr_max = 0.
- qs_max = 0.
- qi_max = 0.
- qg_max = 0
- ni_max = 0.
- nr_max = 0.
- imax_qc = 0
- imax_qr = 0
- imax_qi = 0
- imax_qs = 0
- imax_qg = 0
- imax_ni = 0
- imax_nr = 0
- jmax_qc = 0
- jmax_qr = 0
- jmax_qi = 0
- jmax_qs = 0
- jmax_qg = 0
- jmax_ni = 0
- jmax_nr = 0
- kmax_qc = 0
- kmax_qr = 0
- kmax_qi = 0
- kmax_qs = 0
- kmax_qg = 0
- kmax_ni = 0
- kmax_nr = 0
- do i = 1, 256
- mp_debug(i:i) = char(0)
- enddo
- j_loop: do j = j_start, j_end
- i_loop: do i = i_start, i_end
- pptrain = 0.
- pptsnow = 0.
- pptgraul = 0.
- pptice = 0.
- RAINNCV(i,j) = 0.
- IF ( PRESENT (snowncv) ) THEN
- SNOWNCV(i,j) = 0.
- ENDIF
- IF ( PRESENT (graupelncv) ) THEN
- GRAUPELNCV(i,j) = 0.
- ENDIF
- SR(i,j) = 0.
- do k = kts, kte
- t1d(k) = th(i,k,j)*pii(i,k,j)
- p1d(k) = p(i,k,j)
- dz1d(k) = dz(i,k,j)
- qv1d(k) = qv(i,k,j)
- qc1d(k) = qc(i,k,j)
- qi1d(k) = qi(i,k,j)
- qr1d(k) = qr(i,k,j)
- qs1d(k) = qs(i,k,j)
- qg1d(k) = qg(i,k,j)
- ni1d(k) = ni(i,k,j)
- nr1d(k) = nr(i,k,j)
- enddo
- call mp_thompson(qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
- nr1d, t1d, p1d, dz1d, &
- pptrain, pptsnow, pptgraul, pptice, &
- #ifdef WRF_CHEM
- rainprod1d, evapprod1d, &
- #endif
- kts, kte, dt, i, j)
- pcp_ra(i,j) = pptrain
- pcp_sn(i,j) = pptsnow
- pcp_gr(i,j) = pptgraul
- pcp_ic(i,j) = pptice
- RAINNCV(i,j) = pptrain + pptsnow + pptgraul + pptice
- RAINNC(i,j) = RAINNC(i,j) + pptrain + pptsnow + pptgraul + pptice
- IF ( PRESENT(snowncv) .AND. PRESENT(snownc) ) THEN
- SNOWNCV(i,j) = pptsnow + pptice
- SNOWNC(i,j) = SNOWNC(i,j) + pptsnow + pptice
- ENDIF
- IF ( PRESENT(graupelncv) .AND. PRESENT(graupelnc) ) THEN
- GRAUPELNCV(i,j) = pptgraul
- GRAUPELNC(i,j) = GRAUPELNC(i,j) + pptgraul
- ENDIF
- SR(i,j) = (pptsnow + pptgraul + pptice)/(RAINNCV(i,j)+1.e-12)
- do k = kts, kte
- qv(i,k,j) = qv1d(k)
- qc(i,k,j) = qc1d(k)
- qi(i,k,j) = qi1d(k)
- qr(i,k,j) = qr1d(k)
- qs(i,k,j) = qs1d(k)
- qg(i,k,j) = qg1d(k)
- ni(i,k,j) = ni1d(k)
- nr(i,k,j) = nr1d(k)
- #ifdef WRF_CHEM
- rainprod(i,k,j) = rainprod1d(k)
- evapprod(i,k,j) = evapprod1d(k)
- #endif
- th(i,k,j) = t1d(k)/pii(i,k,j)
- 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.kte-2 .and. k.gt.kts+1) then
- qv(i,k,j) = 0.5*(qv(i,k-1,j) + qv(i,k+1,j))
- else
- qv(i,k,j) = 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, qc1d, qr1d, nr1d, qs1d, qg1d, &
- ! t1d, p1d, dBZ, kts, kte, i, j)
- ! do k = kts, kte
- ! refl_10cm(i,k,j) = MAX(-35., dBZ(k))
- ! enddo
- ! endif
- enddo i_loop
- enddo j_loop
- ! DEBUG - GT
- write(mp_debug,'(a,7(a,e13.6,1x,a,i3,a,i3,a,i3,a,1x))') 'MP-GT:', &
- 'qc: ', qc_max, '(', imax_qc, ',', jmax_qc, ',', kmax_qc, ')', &
- 'qr: ', qr_max, '(', imax_qr, ',', jmax_qr, ',', kmax_qr, ')', &
- 'qi: ', qi_max, '(', imax_qi, ',', jmax_qi, ',', kmax_qi, ')', &
- 'qs: ', qs_max, '(', imax_qs, ',', jmax_qs, ',', kmax_qs, ')', &
- 'qg: ', qg_max, '(', imax_qg, ',', jmax_qg, ',', kmax_qg, ')', &
- 'ni: ', ni_max, '(', imax_ni, ',', jmax_ni, ',', kmax_ni, ')', &
- 'nr: ', nr_max, '(', imax_nr, ',', jmax_nr, ',', kmax_nr, ')'
- CALL wrf_debug(150, mp_debug)
- ! END DEBUG - GT
- do i = 1, 256
- mp_debug(i:i) = char(0)
- enddo
- END SUBROUTINE mp_gt_driver
- !+---+-----------------------------------------------------------------+
- !ctrlL
- !+---+-----------------------------------------------------------------+
- !+---+-----------------------------------------------------------------+
- !.. This subroutine computes the moisture tendencies of water vapor,
- !.. cloud droplets, rain, cloud ice (pristine), snow, and graupel.
- !.. Previously this code was based on Reisner et al (1998), but few of
- !.. those pieces remain. A complete description is now found in
- !.. Thompson et al. (2004, 2008).
- !+---+-----------------------------------------------------------------+
- !
- subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
- nr1d, t1d, p1d, dzq, &
- pptrain, pptsnow, pptgraul, pptice, &
- #ifdef WRF_CHEM
- rainprod, evapprod, &
- #endif
- kts, kte, dt, ii, jj)
- implicit none
- !..Sub arguments
- INTEGER, INTENT(IN):: kts, kte, ii, jj
- REAL, DIMENSION(kts:kte), INTENT(INOUT):: &
- qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
- nr1d, t1d, p1d
- #ifdef WRF_CHEM
- REAL, DIMENSION(kts:kte), INTENT(INOUT):: &
- rainprod, evapprod
- #endif
- REAL, DIMENSION(kts:kte), INTENT(IN):: dzq
- REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice
- REAL, INTENT(IN):: dt
- !..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
- DOUBLE PRECISION, PARAMETER:: zeroD0 = 0.0d0
- 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
- 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, Ef_rr
- REAL:: dtsave, odts, odt, odzq
- REAL:: xslw1, ygra1, zans1
- 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_g1, idx_g, idx_r1, idx_r, &
- idx_i1, idx_i, idx_c, 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.315 .and. jj.eq.2) 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
- #ifdef WRF_CHEM
- do k = kts, kte
- rainprod(k) = 0.
- evapprod(k) = 0.
- enddo
- #endif
- !+---+-----------------------------------------------------------------+
- !..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(R2, 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. 20.E-6) then
- lami = cie(2)/20.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) = R2
- L_qi(k) = .false.
- endif
- if (qr1d(k) .gt. R1) then
- no_micro = .false.
- rr(k) = qr1d(k)*rho(k)
- nr(k) = MAX(R2, nr1d(k)*rho(k))
- L_qr(k) = .true.
- lamr = (am_r*crg(3)*org2*nr(k)/rr(k))**obmr
- mvd_r(k) = (3.0 + mu_r + 0.672) / lamr
- if (mvd_r(k) .gt. 2.5E-3) then
- mvd_r(k) = 2.5E-3
- lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
- nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
- elseif (mvd_r(k) .lt. D0r*0.75) then
- mvd_r(k) = D0r*0.75
- lamr = (3.0 + mu_r + 0.672) / mvd_r(k)
- nr(k) = crg(2)*org3*rr(k)*lamr**bm_r / am_r
- endif
- else
- qr1d(k) = 0.0
- nr1d(k) = 0.0
- rr(k) = R1
- nr(k) = R2
- L_qr(k) = .false.
- endif
- if (qs1d(k) .gt. R1) then
- no_micro = .false.
- rs(k) = qs1d(k)*rho(k)
- L_qs(k) = .true.
- else
- qs1d(k) = 0.0
- rs(k) = R1
- L_qs(k) = .false.
- endif
- if (qg1d(k) .gt. R1) then
- no_micro = .false.
- rg(k) = qg1d(k)*rho(k)
- L_qg(k) = .true.
- else
- qg1d(k) = 0.0
- rg(k) = R1
- L_qg(k) = .false.
- endif
- enddo
- !+---+-----------------------------------------------------------------+
- !..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
- 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.
- !+---+-------…
Large files files are truncated, but you can click here to view the full file