From def46ffe8f4a785444e72b28e8ec49064bb5acd3 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 4 Jun 2020 15:59:08 -0600 Subject: [PATCH 1/5] Implement option to roll back Thompson MP to WRFV3.8.1 used in RAPv5/HRRRv4 --- physics/module_mp_thompson.F90 | 187 +++++++++++++++++++++++++++++---- 1 file changed, 167 insertions(+), 20 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index b3ccb7412..191070b62 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1,6 +1,12 @@ !>\file module_mp_thompson.F90 !! This file contains the entity of GSD Thompson MP scheme. +! DH* 2020-06-05 +! Use the following preprocessor directive to roll back +! to the WRFv3.8.1, used in RAPv5/HRRRv4 for more reasonable +! representation of mesoscale storms and reflectivity values +!#define WRF381 + !>\ingroup aathompson !! This module computes the moisture tendencies of water vapor, @@ -43,9 +49,16 @@ !!\author Greg Thompson, NCAR-RAL, gthompsn@ucar.edu, 303-497-2805 !! !! - Last modified: 24 Jan 2018 Aerosol additions to v3.5.1 code 9/2013 -!! Cloud fraction additions 11/2014 part of pre-v3.7 +!! Cloud fraction additions 11/2014 part of pre-v3.7 !! - Imported in CCPP by: Dom Heinzeller, NOAA/ESRL/GSD, dom.heinzeller@noaa.gov !! - Last modified: 6 Aug 2018 Update of initial import to WRFV4.0 +!! - Last modified: 13 Mar 2020 Add logic to turtn on/off the calculation +!! of melting layer in radar reflectivity routine +!! - Last modified: 2 Jun 2020 Add option to rollback to version 3.8.1 +!! used in RAPv5/HRRRv4, include stochastic physics +!! perturbations to the graupel intercept parameter, +!! the cloud water shape parameter, and the number +!! concentration of nucleated aerosols. MODULE module_mp_thompson USE machine, only : kind_phys @@ -450,6 +463,13 @@ SUBROUTINE thompson_init(nwfa2d, nifa2d, nwfa, nifa, & if (.NOT. ALLOCATED(tcg_racg) ) then ALLOCATE(tcg_racg(ntb_g1,ntb_g,ntb_r1,ntb_r)) micro_init = .TRUE. + if (mpirank==mpiroot) then +#ifdef WRF381 + write(0,*) "Using Thompson MP from WRFv3.8.1 (RAPv5/HRRRv4)" +#else + write(0,*) "Using Thompson MP from WRFv4.0+" +#endif + endif endif if (.NOT. ALLOCATED(tmr_racg)) ALLOCATE(tmr_racg(ntb_g1,ntb_g,ntb_r1,ntb_r)) @@ -961,14 +981,6 @@ SUBROUTINE thompson_init(nwfa2d, nifa2d, nwfa, nifa, & end if precomputed_tables_2 - ! DH* TEMPORARY GUARD 20181203 - if (minval(tnccn_act)==maxval(tnccn_act)) then - write(0,*) "TEMPORARY GUARD: abort model because table_ccnact seems to be faulty." - call sleep(5) - stop - end if - ! *DH - endif if_not_iiwarm if (mpirank==mpiroot) write(0,*) ' ... DONE microphysical lookup tables' @@ -997,6 +1009,9 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & vt_dbz_wt, first_time_step, & re_cloud, re_ice, re_snow, & has_reqc, has_reqi, has_reqs, & + rand_perturb_on, & + kme_stoch, & + rand_pert, & ids,ide, jds,jde, kds,kde, & ! domain dims ims,ime, jms,jme, kms,kme, & ! memory dims its,ite, jts,jte, kts,kte, & ! tile dims @@ -1019,6 +1034,10 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(OUT):: & re_cloud, re_ice, re_snow + INTEGER, INTENT(IN) :: rand_perturb_on, kme_stoch + REAL, DIMENSION(ims:ime,kms:kme_stoch,jms:jme), INTENT(IN), OPTIONAL:: & + rand_pert + INTEGER, INTENT(IN):: has_reqc, has_reqi, has_reqs #if ( WRF_CHEM == 1 ) REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: & @@ -1054,7 +1073,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & REAL:: dt, pptrain, pptsnow, pptgraul, pptice REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max REAL:: nwfa1 - INTEGER:: i, j, k + REAL:: rand1, rand2, rand3, min_rand + INTEGER:: i, j, k, m 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 @@ -1160,6 +1180,32 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & j_loop: do j = j_start, j_end i_loop: do i = i_start, i_end +!+---+-----------------------------------------------------------------+ +!..Introduce stochastic parameter perturbations by creating as many scalar rand1, rand2, ... +!.. variables as needed to perturb different pieces of microphysics. gthompsn 21Mar2018 +! Setting spp_mp to 1 gives graupel Y-intercept pertubations (2^0) +! 2 gives cloud water distribution gamma shape parameter perturbations (2^1) +! 4 gives CCN & IN activation perturbations (2^2) +! 3 gives both 1+2 +! 5 gives both 1+4 +! 6 gives both 2+4 +! 7 gives all 1+2+4 +! For now (22Mar2018), standard deviation should be only 0.25 and cut-off at 1.5 +! in order to constrain the various perturbations from being too extreme. +!+---+-----------------------------------------------------------------+ + rand1 = 0.0 + rand2 = 0.0 + rand3 = 0.0 + if (rand_perturb_on .ne. 0) then + if (MOD(rand_perturb_on,2) .ne. 0) rand1 = rand_pert(i,1,j) + m = RSHIFT(ABS(rand_perturb_on),1) + if (MOD(m,2) .ne. 0) rand2 = rand_pert(i,1,j)*2. + m = RSHIFT(ABS(rand_perturb_on),2) + if (MOD(m,2) .ne. 0) rand3 = 0.1*(rand_pert(i,1,j)+ABS(min_rand)) + m = RSHIFT(ABS(rand_perturb_on),3) + endif +!+---+-----------------------------------------------------------------+ + pptrain = 0. pptsnow = 0. pptgraul = 0. @@ -1218,6 +1264,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & #if ( WRF_CHEM == 1 ) rainprod1d, evapprod1d, & #endif + rand1, rand2, rand3, & kts, kte, dt, i, j) pcp_ra(i,j) = pptrain @@ -1485,6 +1532,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & #if ( WRF_CHEM == 1 ) rainprod, evapprod, & #endif + rand1, rand2, rand3, & kts, kte, dt, ii, jj) #ifdef MPI use mpi @@ -1499,6 +1547,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & REAL, DIMENSION(kts:kte), INTENT(IN):: p1d, w1d, dzq REAL, INTENT(INOUT):: pptrain, pptsnow, pptgraul, pptice REAL, INTENT(IN):: dt + REAL, INTENT(IN):: rand1, rand2, rand3 + #if ( WRF_CHEM == 1 ) REAL, DIMENSION(kts:kte), INTENT(INOUT):: & rainprod, evapprod @@ -1735,7 +1785,12 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & rc(k) = qc1d(k)*rho(k) nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max)) L_qc(k) = .true. - nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) + if (rand2 .eq. 0.0) then + nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) + else + nu_c = NINT(1000.E6/nc(k)) + 2 + nu_c = MAX(2, MIN(nu_c+NINT(rand2), 15)) + endif lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr xDc = (bm_r + nu_c + 1.) / lamc if (xDc.lt. D0c) then @@ -1984,7 +2039,10 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & xslw1 = 0.01 endif ygra1 = 4.31 + alog10(max(5.E-5, rg(k))) - zans1 = 3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1)) + zans1 = (3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1))) + rand1 + if (rand1 .ne. 0.0) then + zans1 = MAX(2., MIN(zans1, 7.)) + endif N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) N0_min = MIN(N0_exp, N0_min) @@ -2025,7 +2083,12 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & mvd_c(k) = D0c if (L_qc(k)) then - nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) + if (rand2 .eq. 0.0) then + nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) + else + nu_c = NINT(1000.E6/nc(k)) + 2 + nu_c = MAX(2, MIN(nu_c+NINT(rand2), 15)) + endif xDc = MAX(D0c*1.E6, ((rc(k)/(am_r*nc(k)))**obmr) * 1.E6) lamc = (nc(k)*am_r* ccg(2,nu_c) * ocg1(nu_c) / rc(k))**obmr mvd_c(k) = (3.0+nu_c+0.672) / lamc @@ -2427,6 +2490,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & .and. temp(k).lt.253.15) ) then if (dustyIce .AND. is_aerosol_aware) then xnc = iceDeMott(tempc,qv(k),qvs(k),qvsi(k),rho(k),nifa(k)) + xnc = xnc*(1.0 + 3.*rand3) else xnc = MIN(250.E3, TNO*EXP(ATO*(T_0-temp(k)))) endif @@ -2633,7 +2697,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !! supersat again. sump = pri_inu(k) + pri_ide(k) + prs_ide(k) & + prs_sde(k) + prg_gde(k) + pri_iha(k) +! DH* 2020-06-02 I believe that the WRF381 version +! is wrong, because the units do not match. +#ifdef WRF381 + rate_max = (qv(k)-qvsi(k))*odts*0.999 +#else rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999 +#endif if ( (sump.gt. eps .and. sump.gt. rate_max) .or. & (sump.lt. -eps .and. sump.lt. rate_max) ) then ratio = rate_max/sump @@ -2765,7 +2835,12 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & xrc=MAX(R1, (qc1d(k) + qcten(k)*dtsave)*rho(k)) xnc=MAX(2., (nc1d(k) + ncten(k)*dtsave)*rho(k)) if (xrc .gt. R1) then - nu_c = MIN(15, NINT(1000.E6/xnc) + 2) + if (rand2 .eq. 0.0) then + nu_c = MIN(15, NINT(1000.E6/xnc) + 2) + else + nu_c = NINT(1000.E6/xnc) + 2 + nu_c = MAX(2, MIN(nu_c+NINT(rand2), 15)) + endif lamc = (xnc*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr xDc = (bm_r + nu_c + 1.) / lamc if (xDc.lt. D0c) then @@ -3055,7 +3130,10 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & xslw1 = 0.01 endif ygra1 = 4.31 + alog10(max(5.E-5, rg(k))) - zans1 = 3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1)) + zans1 = (3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1))) + rand1 + if (rand1 .ne. 0.0) then + zans1 = MAX(2., MIN(zans1, 7.)) + endif N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) N0_min = MIN(N0_exp, N0_min) @@ -3103,7 +3181,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ ! DROPLET NUCLEATION if (clap .gt. eps) then if (is_aerosol_aware) then - xnc = MAX(2., activ_ncloud(temp(k), w1d(k), nwfa(k))) + xnc = MAX(2., activ_ncloud(temp(k), w1d(k)+rand3, nwfa(k))) else xnc = Nt_c endif @@ -3342,7 +3420,12 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & do k = ksed1(5), kts, -1 vtc = 0. if (rc(k) .gt. R1 .and. w1d(k) .lt. 1.E-1) then - nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) + if (rand2 .eq. 0.0) then + nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) + else + nu_c = NINT(1000.E6/nc(k)) + 2 + nu_c = MAX(2, MIN(nu_c+NINT(rand2), 15)) + endif lamc = (nc(k)*am_r*ccg(2,nu_c)*ocg1(nu_c)/rc(k))**obmr ilamc = 1./lamc vtc = rhof(k)*av_c*ccg(5,nu_c)*ocg2(nu_c) * ilamc**bv_c @@ -3408,7 +3491,12 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts) if (temp(k).gt. (T_0+0.1)) then vtsk(k) = MAX(vts*vts_boost(k), & - & vts*((vtrk(k)-vts*vts_boost(k))/(temp(k)-T_0))) + & vts*((vtrk(k)-vts*vts_boost(k))/(temp(k)-T_0))) ! +! DH* The version below is supposed to be a better formulation, +! but gave worse results in RAPv5/HRRRv4 than the line above. + ! this formulation for RAPv5/HRRRv4, reverted 20 Feb 2020 + ! SR = rs(k)/(rs(k)+rr(k)) ! bug fix from G. Thompson, 10 May 2019 + ! vtsk(k) = vts*SR + (1.-SR)*vtrk(k) else vtsk(k) = vts*vts_boost(k) endif @@ -3459,6 +3547,10 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !> - 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. +!.. Bug fix, 2013Nov01 to tendencies using rho(k+1) correction thanks to +!.. Eric Skyllingstad. !+---+-----------------------------------------------------------------+ if (ANY(L_qr .eqv. .true.)) then @@ -3488,7 +3580,11 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & *odzq*DT*onstep(1)) enddo +#ifdef WRF381 + if (rr(kts).gt.R1*10.) & +#else if (rr(kts).gt.R1*1000.) & +#endif pptrain = pptrain + sed_r(kts)*DT*onstep(1) enddo endif @@ -3539,7 +3635,11 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & *odzq*DT*onstep(2)) enddo +#ifdef WRF381 + if (ri(kts).gt.R1*10.) & +#else if (ri(kts).gt.R1*1000.) & +#endif pptice = pptice + sed_i(kts)*DT*onstep(2) enddo endif @@ -3566,7 +3666,11 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & *odzq*DT*onstep(3)) enddo +#ifdef WRF381 + if (rs(kts).gt.R1*10.) & +#else if (rs(kts).gt.R1*1000.) & +#endif pptsnow = pptsnow + sed_s(kts)*DT*onstep(3) enddo endif @@ -3593,7 +3697,11 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & *odzq*DT*onstep(4)) enddo +#ifdef WRF381 + if (rg(kts).gt.R1*10.) & +#else if (rg(kts).gt.R1*1000.) & +#endif pptgraul = pptgraul + sed_g(kts)*DT*onstep(4) enddo endif @@ -3634,16 +3742,31 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & qv1d(k) = MAX(1.E-10, qv1d(k) + qvten(k)*DT) qc1d(k) = qc1d(k) + qcten(k)*DT nc1d(k) = MAX(2./rho(k), MIN(nc1d(k) + ncten(k)*DT, Nt_c_max)) +! DH* 2020-06-05 I believe WRF381 is wrong in terms of units; +! dividing by rho turns number concentration per volume into +! number concentration per mass. +#ifdef WRF381 + nwfa1d(k) = MAX(11.1E6/rho(k), MIN(9999.E6/rho(k), & + (nwfa1d(k)+nwfaten(k)*DT))) + nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6/rho(k), & + (nifa1d(k)+nifaten(k)*DT))) +#else nwfa1d(k) = MAX(11.1E6, MIN(9999.E6, & (nwfa1d(k)+nwfaten(k)*DT))) nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6, & (nifa1d(k)+nifaten(k)*DT))) +#endif if (qc1d(k) .le. R1) then qc1d(k) = 0.0 nc1d(k) = 0.0 else - nu_c = MIN(15, NINT(1000.E6/(nc1d(k)*rho(k))) + 2) + if (rand2 .eq. 0.0) then + nu_c = MIN(15, NINT(1000.E6/(nc1d(k)*rho(k))) + 2) + else + nu_c = NINT(1000.E6/(nc1d(k)*rho(k))) + 2 + nu_c = MAX(2, MIN(nu_c+NINT(rand2), 15)) + endif lamc = (am_r*ccg(2,nu_c)*ocg1(nu_c)*nc1d(k)/qc1d(k))**obmr xDc = (bm_r + nu_c + 1.) / lamc if (xDc.lt. D0c) then @@ -5124,7 +5247,14 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & do k = kts, kte rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) rc(k) = MAX(R1, qc1d(k)*rho(k)) +#ifdef WRF381 + nc(k) = MAX(R2, MIN(nc1d(k)*rho(k), Nt_c_max)) +#else + ! DH* 2020-06-05 is using 2.0 instead of R2 + ! a bug in the WRFv4.0+ version of Thompson? + ! For ni(k) a few lines below, it is still R2 nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max)) +#endif if (.NOT. is_aerosol_aware) nc(k) = Nt_c if (rc(k).gt.R1 .and. nc(k).gt.R2) has_qc = .true. ri(k) = MAX(R1, qi1d(k)*rho(k)) @@ -5136,7 +5266,9 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & if (has_qc) then do k = kts, kte +#ifndef WRF381 re_qc1d(k) = 2.49E-6 +#endif if (rc(k).le.R1 .or. nc(k).le.R2) CYCLE if (nc(k).lt.100) then inu_c = 15 @@ -5152,16 +5284,24 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & if (has_qi) then do k = kts, kte +#ifndef WRF381 re_qi1d(k) = 2.49E-6 +#endif if (ri(k).le.R1 .or. ni(k).le.R2) CYCLE lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi +#ifdef WRF381 + re_qi1d(k) = MAX(5.01E-6, MIN(SNGL(0.5D0 * DBLE(3.+mu_i)/lami), 125.E-6)) +#else re_qi1d(k) = MAX(2.51E-6, MIN(SNGL(0.5D0 * DBLE(3.+mu_i)/lami), 125.E-6)) +#endif enddo endif if (has_qs) then do k = kts, kte +#ifndef WRF381 re_qs1d(k) = 4.99E-6 +#endif if (rs(k).le.R1) CYCLE tc0 = MIN(-0.1, t1d(k)-273.15) smob = rs(k)*oams @@ -5196,7 +5336,11 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) smoc = a_ * smo2**b_ +#ifdef WRF381 + re_qs1d(k) = MAX(10.E-6, MIN(0.5*(smoc/smob), 999.E-6)) +#else re_qs1d(k) = MAX(5.01E-6, MIN(0.5*(smoc/smob), 999.E-6)) +#endif enddo endif @@ -5383,7 +5527,10 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & xslw1 = 0.01 endif ygra1 = 4.31 + alog10(max(5.E-5, rg(k))) - zans1 = 3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1)) + zans1 = (3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1))) + rand1 + if (rand1 .ne. 0.0) then + zans1 = MAX(2., MIN(zans1, 7.)) + endif N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) N0_min = MIN(N0_exp, N0_min) From c41d691d400fc6f4d3132602f127127cfb3f2ccc Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 4 Jun 2020 16:09:45 -0600 Subject: [PATCH 2/5] physics/module_mp_thompson.F90: add guard to prevent running Thompson MP with the untested stochastic perturbations code --- physics/module_mp_thompson.F90 | 11 +++++++++++ 1 file changed, 11 insertions(+) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 191070b62..ce6df30e3 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1091,6 +1091,17 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & if (present(errmsg)) errmsg = '' if (present(errflg)) errflg = 0 + ! DH* 2020-06-05: The stochastic perturbations code was retrofitted + ! from a newer version of the Thompson MP scheme, but it has not been + ! tested yet. + if (rand_perturb_on .ne. 0) then + errmsg = 'Logic error in mp_gt_driver: the stochastic perturbations code ' // & + 'has not been tested yet with this version of the Thompson scheme' + errflg = 1 + return + end if + ! *DH 2020-06-05 + if ( (present(tt) .and. (present(th) .or. present(pii))) .or. & (.not.present(tt) .and. .not.(present(th) .and. present(pii))) ) then if (present(errmsg)) then From 990ffbad9219d0ef595a8d5f8f0364d88145e451 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 5 Jun 2020 14:19:47 -0600 Subject: [PATCH 3/5] Add stochastic perturbation variables to mp_thompson.F90, bugfix in module_mp_thompson.F90 --- physics/module_mp_thompson.F90 | 22 +++++++++++++++------- physics/mp_thompson.F90 | 14 ++++++++++++++ 2 files changed, 29 insertions(+), 7 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index ce6df30e3..532071a8e 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1100,6 +1100,13 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & errflg = 1 return end if + ! Activate this code when removing the guard above + !if (rand_perturb_on .ne. 0 .and. .not. present(rand_pert)) then + ! errmsg = 'Logic error in mp_gt_driver: random perturbations are on, ' // & + ! 'but optional argument rand_pert is not present' + ! errflg = 1 + ! return + !end if ! *DH 2020-06-05 if ( (present(tt) .and. (present(th) .or. present(pii))) .or. & @@ -1428,13 +1435,13 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & endif ! if (present(vt_dbz_wt) .and. present(first_time_step)) then - call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & - t1d, p1d, dBZ, kts, kte, i, j, & - melti, vt_dbz_wt(i,:,j), & + call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & + t1d, p1d, dBZ, rand1, kts, kte, i, j, & + melti, vt_dbz_wt(i,:,j), & first_time_step) else - call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & - t1d, p1d, dBZ, kts, kte, i, j, & + call calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & + t1d, p1d, dBZ, rand1, kts, kte, i, j, & melti) end if do k = kts, kte @@ -5366,13 +5373,14 @@ end subroutine calc_effectRad !! of frozen species remaining from what initially existed at the !! melting level interface. subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & - t1d, p1d, dBZ, kts, kte, ii, jj, melti, vt_dBZ, & - first_time_step) + t1d, p1d, dBZ, rand1, kts, kte, ii, jj, melti, & + vt_dBZ, first_time_step) IMPLICIT NONE !..Sub arguments INTEGER, INTENT(IN):: kts, kte, ii, jj + REAL, INTENT(IN):: rand1 REAL, DIMENSION(kts:kte), INTENT(IN):: & qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, t1d, p1d REAL, DIMENSION(kts:kte), INTENT(INOUT):: dBZ diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 3f2ee144e..1653c825d 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -472,6 +472,12 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & integer :: has_reqc integer :: has_reqi integer :: has_reqs + ! DH* 2020-06-05 hardcode these values for not using random perturbations, + ! hasn't been tested yet with this version of module_mp_thompson.F90 + integer, parameter :: rand_perturb_on = 0 + integer, parameter :: kme_stoch = 1 + !real(kind_phys) :: rand_pert(1:ncol,1:kme_stoch) + ! *DH 2020-06-05 ! Dimensions used in mp_gt_driver integer :: ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & @@ -601,6 +607,10 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & re_cloud=re_cloud, re_ice=re_ice, re_snow=re_snow, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & + rand_perturb_on=rand_perturb_on, kme_stoch=kme_stoch, & + ! DH* 2020-06-05 not passing this optional argument, see + ! comment in module_mp_thompson.F90 / mp_gt_driver + !rand_pert=rand_pert, & ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & @@ -618,6 +628,10 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & re_cloud=re_cloud, re_ice=re_ice, re_snow=re_snow, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & + rand_perturb_on=rand_perturb_on, kme_stoch=kme_stoch, & + ! DH* 2020-06-05 not passing this optional argument, see + ! comment in module_mp_thompson.F90 / mp_gt_driver + !rand_pert=rand_pert, & ids=ids, ide=ide, jds=jds, jde=jde, kds=kds, kde=kde, & ims=ims, ime=ime, jms=jms, jme=jme, kms=kms, kme=kme, & its=its, ite=ite, jts=jts, jte=jte, kts=kts, kte=kte, & From 626ec0e4c8a96b3599caaf18a28691c842bb2b9d Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 9 Jun 2020 08:15:53 -0600 Subject: [PATCH 4/5] Clean up of effective radii calculation for Thompson MP: move initialization and bounds into the calc_effectRad routine, use settings consistent with previous version of code --- physics/GFS_rrtmg_pre.F90 | 19 +++++++-------- physics/module_mp_thompson.F90 | 44 ++++++++++++++++++++++++---------- physics/mp_thompson.F90 | 15 +----------- 3 files changed, 40 insertions(+), 38 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 42411c88f..381fa159f 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -716,23 +716,20 @@ subroutine GFS_rrtmg_pre_run (Model, Grid, Sfcprop, Statein, & ! input end do ! Call Thompson's subroutine to compute effective radii do i=1,im - ! Initialize to default in units m as in module_mp_thompson.F90 - re_cloud(i,:) = 2.49E-6 - re_ice(i,:) = 4.99E-6 - re_snow(i,:) = 9.99E-6 + ! Effective radii [m] are now intent(out), bounds applied in calc_effectRad + !tgs: progclduni has different limits for ice radii (10.0-150.0) than + ! calc_effectRad (4.99-125.0 for WRFv3.8.1; 2.49-125.0 for WRFv4+) + ! it will raise the low limit from 5 to 10, but the high limit will remain 125. call calc_effectRad (tlyr(i,:), plyr(i,:), qv_mp(i,:), qc_mp(i,:), & nc_mp(i,:), qi_mp(i,:), ni_mp(i,:), qs_mp(i,:), & re_cloud(i,:), re_ice(i,:), re_snow(i,:), 1, lm ) end do - ! Scale Thompson's effective radii from meter to micron and apply bounds + ! Scale Thompson's effective radii from meter to micron do k=1,lm do i=1,im - re_cloud(i,k) = MAX(2.49, MIN(re_cloud(i,k)*1.e6, 50.)) - re_ice(i,k) = MAX(4.99, MIN(re_ice(i,k)*1.e6, 125.)) - !tgs: progclduni has different limits for ice radii: 10.0-150.0 - ! it will raise the low limit from 5 to 10, but the - ! high limit will remain 125. - re_snow(i,k) = MAX(9.99, MIN(re_snow(i,k)*1.e6, 999.)) + re_cloud(i,k) = re_cloud(i,k)*1.e6 + re_ice(i,k) = re_ice(i,k)*1.e6 + re_snow(i,k) = re_snow(i,k)*1.e6 end do end do do k=1,lm diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 532071a8e..705d245ae 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -5246,7 +5246,7 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & INTEGER, INTENT(IN):: kts, kte REAL, DIMENSION(kts:kte), INTENT(IN):: & & t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d - REAL, DIMENSION(kts:kte), INTENT(INOUT):: re_qc1d, re_qi1d, re_qs1d + REAL, DIMENSION(kts:kte), INTENT(OUT):: re_qc1d, re_qi1d, re_qs1d !..Local variables INTEGER:: k REAL, DIMENSION(kts:kte):: rho, rc, nc, ri, ni, rs @@ -5262,6 +5262,30 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & has_qi = .false. has_qs = .false. +! DH* 2020-06-08 Moved the initial values and bounds from +! the calling routines into calc_effectRad (to prevent +! multiple definitions that may be inconsistent). The +! initial values and bounds from the calling routines were +! +! re_cloud(i,k) = MAX(2.49, MIN(re_cloud(i,k)*1.e6, 50.)) +! re_ice(i,k) = MAX(4.99, MIN(re_ice(i,k)*1.e6, 125.)) +! re_snow(i,k) = MAX(9.99, MIN(re_snow(i,k)*1.e6, 999.)) +! +! independent of the version of Thompson MP. These values +! are consistent with the WRFv3.8.1 settings, but inconsistent +! with the WRFv4+ settings. In order to apply the same bounds +! as before this change, use the WRF v3.8.1 settings throughout. +#if 1 +!ifdef WRF381 + re_qc1d(:) = 2.49E-6 + re_qi1d(:) = 4.99E-6 + re_qs1d(:) = 9.99E-6 +#else + re_qc1d(:) = 2.49E-6 + re_qi1d(:) = 2.49E-6 + re_qs1d(:) = 4.99E-6 +#endif + do k = kts, kte rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) rc(k) = MAX(R1, qc1d(k)*rho(k)) @@ -5270,7 +5294,8 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & #else ! DH* 2020-06-05 is using 2.0 instead of R2 ! a bug in the WRFv4.0+ version of Thompson? - ! For ni(k) a few lines below, it is still R2 + ! For ni(k) a few lines below, it is still R2. + ! Note that R2 is defined as R2 = 1.E-6 nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max)) #endif if (.NOT. is_aerosol_aware) nc(k) = Nt_c @@ -5284,9 +5309,6 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & if (has_qc) then do k = kts, kte -#ifndef WRF381 - re_qc1d(k) = 2.49E-6 -#endif if (rc(k).le.R1 .or. nc(k).le.R2) CYCLE if (nc(k).lt.100) then inu_c = 15 @@ -5302,12 +5324,10 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & if (has_qi) then do k = kts, kte -#ifndef WRF381 - re_qi1d(k) = 2.49E-6 -#endif if (ri(k).le.R1 .or. ni(k).le.R2) CYCLE lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi -#ifdef WRF381 +#if 1 +!ifdef WRF381 re_qi1d(k) = MAX(5.01E-6, MIN(SNGL(0.5D0 * DBLE(3.+mu_i)/lami), 125.E-6)) #else re_qi1d(k) = MAX(2.51E-6, MIN(SNGL(0.5D0 * DBLE(3.+mu_i)/lami), 125.E-6)) @@ -5317,9 +5337,6 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & if (has_qs) then do k = kts, kte -#ifndef WRF381 - re_qs1d(k) = 4.99E-6 -#endif if (rs(k).le.R1) CYCLE tc0 = MIN(-0.1, t1d(k)-273.15) smob = rs(k)*oams @@ -5354,7 +5371,8 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & & + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) & & + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1) smoc = a_ * smo2**b_ -#ifdef WRF381 +#if 1 +!ifdef WRF381 re_qs1d(k) = MAX(10.E-6, MIN(0.5*(smoc/smob), 999.E-6)) #else re_qs1d(k) = MAX(5.01E-6, MIN(0.5*(smoc/smob), 999.E-6)) diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 1653c825d..ec19945b0 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -320,25 +320,12 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & ! Calculate initial cloud effective radii if requested if (present(re_cloud) .and. present(re_ice) .and. present(re_snow)) then - do i = 1, ncol - do k = 1, nlev - re_cloud(i,k) = 2.49E-6 - re_ice(i,k) = 4.99E-6 - re_snow(i,k) = 9.99E-6 - end do - end do + ! Effective radii [m] are now intent(out), bounds applied in calc_effectRad do i = 1, ncol call calc_effectRad (tgrs(i,:), prsl(i,:), qv_mp(i,:), qc_mp(i,:), & nc_mp(i,:), qi_mp(i,:), ni_mp(i,:), qs_mp(i,:), & re_cloud(i,:), re_ice(i,:), re_snow(i,:), 1, nlev) end do - do i = 1, ncol - do k = 1, nlev - re_cloud(i,k) = MAX(2.49E-6, MIN(re_cloud(i,k), 50.E-6)) - re_ice(i,k) = MAX(4.99E-6, MIN(re_ice(i,k), 125.E-6)) - re_snow(i,k) = MAX(9.99E-6, MIN(re_snow(i,k), 999.E-6)) - end do - end do !! Convert to micron: required for bit-for-bit identical restarts; !! otherwise entering mp_thompson_init and converting mu to m and !! back (without updating re_*) introduces b4b differences. From 4619424a2040d51cf63d2b92dde10c3a42cb02fe Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 9 Jun 2020 10:20:21 -0600 Subject: [PATCH 5/5] physics/module_mp_thompson.F90: update comment on possible bug in nc calculation --- physics/module_mp_thompson.F90 | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 705d245ae..304afc6d5 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -5295,7 +5295,9 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & ! DH* 2020-06-05 is using 2.0 instead of R2 ! a bug in the WRFv4.0+ version of Thompson? ! For ni(k) a few lines below, it is still R2. - ! Note that R2 is defined as R2 = 1.E-6 + ! Note that R2 is defined as R2 = 1.E-6, and is + ! used in other parts of Thompson MP for ni/nr + ! calculations (but not for nc calculations) nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max)) #endif if (.NOT. is_aerosol_aware) nc(k) = Nt_c