From d8f4ffbfbe519ee627c6144505d029c4720bd9f3 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Sat, 13 Feb 2021 07:13:55 -0700 Subject: [PATCH 01/31] update Thompson more similar to WRF-v4.3 --- physics/module_mp_thompson.F90 | 437 +++++++++++++++------------------ 1 file changed, 194 insertions(+), 243 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index c19b594dd..fcbcb8164 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -53,6 +53,10 @@ !! perturbations to the graupel intercept parameter, !! the cloud water shape parameter, and the number !! concentration of nucleated aerosols. +!! - Last modified: 12 Feb 2021 G. Thompson updated to align more closely +!! with his WRF version, including bug fixes and designed +!! changes. + MODULE module_mp_thompson USE machine, only : kind_phys @@ -120,8 +124,8 @@ MODULE module_mp_thompson !.. mixing ratio. Also, when mu_g is non-zero, these become equiv !.. y-intercept for an exponential distrib and proper values are !.. computed based on same mixing ratio and total number concentration. - REAL, PARAMETER, PRIVATE:: gonv_min = 1.E4 - REAL, PARAMETER, PRIVATE:: gonv_max = 3.E6 + REAL, PARAMETER, PRIVATE:: gonv_min = 1.E2 + REAL, PARAMETER, PRIVATE:: gonv_max = 1.E6 !..Mass power law relations: mass = am*D**bm !.. Snow from Field et al. (2005), others assume spherical form. @@ -183,7 +187,7 @@ MODULE module_mp_thompson REAL, PRIVATE:: Sc3 !..Homogeneous freezing temperature - REAL, PARAMETER, PRIVATE:: HGFR = 235.16 !< Homogeneous freezing temperature + REAL, PARAMETER, PRIVATE:: HGFR = 235.16 !..Water vapor and air gas constants at constant pressure REAL, PARAMETER, PRIVATE:: Rv = 461.5 @@ -214,6 +218,15 @@ MODULE module_mp_thompson REAL, PARAMETER, PRIVATE:: D0g = 250.E-6 REAL, PRIVATE:: D0i, xm0s, xm0g +!..Min and max radiative effective radius of cloud water, cloud ice, and snow; +!.. performed by subroutine calc_effectRad + REAL, PARAMETER:: re_qc_min = 2.50E-6 ! 2.5 microns + REAL, PARAMETER:: re_qc_max = 50.0E-6 ! 50 microns + REAL, PARAMETER:: re_qi_min = 2.50E-6 ! 2.5 microns + REAL, PARAMETER:: re_qi_max = 125.0E-6 ! 125 microns + REAL, PARAMETER:: re_qs_min = 5.00E-6 ! 5 microns + REAL, PARAMETER:: re_qs_max = 999.0E-6 ! 999 microns (1 mm) + !..Lookup table dimensions INTEGER, PARAMETER, PRIVATE:: nbins = 100 INTEGER, PARAMETER, PRIVATE:: nbc = nbins @@ -226,7 +239,7 @@ MODULE module_mp_thompson 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_g1 = 37 INTEGER, PARAMETER, PRIVATE:: ntb_r1 = 37 INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55 INTEGER, PARAMETER, PRIVATE:: ntb_t = 9 @@ -299,10 +312,11 @@ MODULE module_mp_thompson !> 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, & + N0g_exp = (/1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, & + 1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, & + 1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, & 1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, & - 1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, & - 1.e7/) + 1.e6/) !> Lookup tables for ice number concentration (/m**3). REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: & @@ -354,6 +368,15 @@ MODULE module_mp_thompson !.. 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). + +!..To permit possible creation of new lookup tables as variables expand/change, +!.. specify a name of external file(s) including version number for pre-computed +!.. Thompson tables. + character(len=*), parameter :: thomp_table_file = 'thompson_tables_precomp_v2.sl' + character(len=*), parameter :: qr_acr_qg_file = 'qr_acr_qgV2.dat' + character(len=*), parameter :: qr_acr_qs_file = 'qr_acr_qsV2.dat' + character(len=*), parameter :: freeze_h2o_file = 'freezeH2O.dat' + INTEGER, PARAMETER, PRIVATE:: R8SIZE = 8 INTEGER, PARAMETER, PRIVATE:: R4SIZE = 4 REAL (KIND=R8SIZE), ALLOCATABLE, DIMENSION(:,:,:,:):: & @@ -748,7 +771,7 @@ SUBROUTINE thompson_init(nwfa2d, nifa2d, nwfa, nifa, & mpi_communicator = mpicomm #ifdef SION call cpu_time(stime) - call readwrite_tables("read", mpicomm, mpirank, mpiroot, ierr) + call readwrite_tables(thomp_table_file, "read", mpicomm, mpirank, mpiroot, ierr) call cpu_time(etime) if (ierr==0) then precomputed_tables = .true. @@ -897,6 +920,10 @@ SUBROUTINE thompson_init(nwfa2d, nifa2d, nwfa, nifa, & if (mpirank==mpiroot) write(0,*) ' creating rain evap table' call table_dropEvap +!> - Call qi_aut_qs() to create conversion of some ice mass into snow category + if (mpirank==mpiroot) write(0,*) ' creating ice converting to snow table' + call qi_aut_qs + call cpu_time(etime) if (mpirank==mpiroot) print '("Calculating Thompson tables part 1 took ",f10.3," seconds.")', etime-stime @@ -945,19 +972,12 @@ SUBROUTINE thompson_init(nwfa2d, nifa2d, nwfa, nifa, & call cpu_time(etime) if (mpirank==mpiroot) print '("Computing freezing of water drops table took ",f10.3," seconds.")', etime-stime -!> - Call qi_aut_qs() to create conversion of some ice mass into snow category - if (mpirank==mpiroot) write(0,*) ' creating ice converting to snow table' - call cpu_time(stime) - call qi_aut_qs - call cpu_time(etime) - if (mpirank==mpiroot) print '("Computing ice converting to snow table took ",f10.3," seconds.")', etime-stime - call cpu_time(etime) if (mpirank==mpiroot) print '("Calculating Thompson tables part 2 took ",f10.3," seconds.")', etime-stime #ifdef SION call cpu_time(stime) - call readwrite_tables("write", mpicomm, mpirank, mpiroot, ierr) + call readwrite_tables(thomp_table_file, "write", mpicomm, mpirank, mpiroot, ierr) if (ierr/=0) then write(0,*) "An error occurred writing Thompson tables to disk" stop 1 @@ -1019,7 +1039,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & nc, nwfa, nifa REAL, DIMENSION(ims:ime, jms:jme), OPTIONAL, INTENT(IN):: nwfa2d, nifa2d - REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(OUT):: & + REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & 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:: & @@ -1059,7 +1079,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & REAL, DIMENSION(its:ite, jts:jte):: pcp_ra, pcp_sn, pcp_gr, pcp_ic REAL:: dt, pptrain, pptsnow, pptgraul, pptice REAL:: qc_max, qr_max, qs_max, qi_max, qg_max, ni_max, nr_max - REAL:: nwfa1 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 @@ -1195,8 +1214,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & ! 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. +! For now (22Mar2018), standard deviation should be up to 0.75 and cut-off at 3.0 +! stddev in order to constrain the various perturbations from being too extreme. !+---+-----------------------------------------------------------------+ rand1 = 0.0 rand2 = 0.0 @@ -1206,7 +1225,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & 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)) + if (MOD(m,2) .ne. 0) rand3 = 0.25*(rand_pert(i,1,j)+ABS(min_rand)) m = RSHIFT(ABS(rand_perturb_on),3) endif !+---+-----------------------------------------------------------------+ @@ -1244,6 +1263,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & qg1d(k) = qg(i,k,j) ni1d(k) = ni(i,k,j) nr1d(k) = nr(i,k,j) + rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) enddo if (is_aerosol_aware) then do k = kts, kte @@ -1251,10 +1271,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & nwfa1d(k) = nwfa(i,k,j) nifa1d(k) = nifa(i,k,j) enddo - nwfa1 = nwfa2d(i,j) else do k = kts, kte - rho(k) = 0.622*p1d(k)/(R*t1d(k)*(qv1d(k)+0.622)) nc1d(k) = Nt_c/rho(k) nwfa1d(k) = 11.1E6/rho(k) nifa1d(k) = naIN1*0.01/rho(k) @@ -1305,7 +1323,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & !.. Changed 13 May 2013 to fake emissions in which nwfa2d is aerosol !.. number tendency (number per kg per second). if (is_aerosol_aware) then -!-GT nwfa1d(kts) = nwfa1 nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt_in nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt_in @@ -1439,17 +1456,17 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & IF (has_reqc.ne.0 .and. has_reqi.ne.0 .and. has_reqs.ne.0) THEN do k = kts, kte - re_qc1d(k) = 2.50E-6 ! 2.49E-6 - re_qi1d(k) = 5.00E-6 ! 4.99E-6 - re_qs1d(k) = 1.00E-5 ! 9.99E-6 + re_qc1d(k) = re_qc_min + re_qi1d(k) = re_qi_min + re_qs1d(k) = re_qs_min enddo !> - Call calc_effectrad() call calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & re_qc1d, re_qi1d, re_qs1d, kts, kte) do k = kts, kte - re_cloud(i,k,j) = MAX(2.50E-6, MIN(re_qc1d(k), 50.E-6)) ! MAX(2.49E-6, MIN(re_qc1d(k), 50.E-6)) - re_ice(i,k,j) = MAX(5.00E-6, MIN(re_qi1d(k), 125.E-6)) ! MAX(4.99E-6, MIN(re_qi1d(k), 125.E-6)) - re_snow(i,k,j) = MAX(1.00E-5, MIN(re_qs1d(k), 999.E-6)) ! MAX(9.99E-6, MIN(re_qs1d(k), 999.E-6)) + re_cloud(i,k,j) = MAX(re_qc_min, MIN(re_qc1d(k), re_qc_max)) + re_ice(i,k,j) = MAX(re_qi_min, MIN(re_qi1d(k), re_qi_max)) + re_snow(i,k,j) = MAX(re_qs_min, MIN(re_qs1d(k), re_qs_max)) enddo ENDIF @@ -1631,7 +1648,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & REAL:: r_frac, g_frac REAL:: Ef_rw, Ef_sw, Ef_gw, Ef_rr REAL:: Ef_ra, Ef_sa, Ef_ga - REAL:: dtsave, odts, odt, odzq, hgt_agl + REAL:: dtsave, odts, odt, odzq, hgt_agl, SR REAL:: xslw1, ygra1, zans1, eva_factor INTEGER:: i, k, k2, n, nn, nstep, k_0, kbot, IT, iexfrq INTEGER, DIMENSION(5):: ksed1 @@ -1785,14 +1802,17 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & nwfa(k) = MAX(11.1E6, MIN(9999.E6, nwfa1d(k)*rho(k))) nifa(k) = MAX(naIN1*0.01, MIN(9999.E6, nifa1d(k)*rho(k))) mvd_r(k) = D0r + mvd_c(k) = D0c if (qc1d(k) .gt. R1) then no_micro = .false. rc(k) = qc1d(k)*rho(k) nc(k) = MAX(2., MIN(nc1d(k)*rho(k), Nt_c_max)) L_qc(k) = .true. - if (rand2 .eq. 0.0) then - nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) + if (nc(k).gt.10000.E6) then + nu_c = 2 + elseif (nc(k).lt.100.) then + nu_c = 15 else nu_c = NINT(1000.E6/nc(k)) + 2 nu_c = MAX(2, MIN(nu_c+NINT(rand2), 15)) @@ -1820,8 +1840,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & ri(k) = qi1d(k)*rho(k) ni(k) = MAX(R2, ni1d(k)*rho(k)) if (ni(k).le. R2) then - lami = cie(2)/25.E-6 - ni(k) = MIN(499.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + lami = cie(2)/5.E-6 + ni(k) = MIN(9999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) endif L_qi(k) = .true. lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi @@ -1829,7 +1849,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & xDi = (bm_i + mu_i + 1.) * ilami if (xDi.lt. 5.E-6) then lami = cie(2)/5.E-6 - ni(k) = MIN(499.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + ni(k) = MIN(9999.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 @@ -2033,26 +2053,11 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ !> - Calculate y-intercept, slope values for graupel. !+---+-----------------------------------------------------------------+ - N0_min = gonv_max - k_0 = kts - do k = kte, kts, -1 - if (temp(k).ge.270.65) k_0 = MAX(k_0, k) - enddo do k = kte, kts, -1 - if (k.gt.k_0 .and. L_qr(k) .and. mvd_r(k).gt.100.E-6) then - xslw1 = 4.01 + alog10(mvd_r(k)) - else - xslw1 = 0.01 - endif - ygra1 = 4.31 + alog10(max(5.E-5, rg(k))) - zans1 = (3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1))) + rand1 - if (rand1 .ne. 0.0) then - zans1 = MAX(2., MIN(zans1, 7.)) - endif + ygra1 = alog10(max(1.E-9, rg(k))) + zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1 N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) - N0_min = MIN(N0_exp, N0_min) - N0_exp = N0_min lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg ilamg(k) = 1./lamg @@ -2087,10 +2092,11 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pnr_rcr(k) = Ef_rr * 2.0*nr(k)*rr(k) endif - mvd_c(k) = D0c if (L_qc(k)) then - if (rand2 .eq. 0.0) then - nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) + if (nc(k).gt.10000.E6) then + nu_c = 2 + elseif (nc(k).lt.100.) then + nu_c = 15 else nu_c = NINT(1000.E6/nc(k)) + 2 nu_c = MAX(2, MIN(nu_c+NINT(rand2), 15)) @@ -2098,6 +2104,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & 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 + mvd_c(k) = MAX(D0c, MIN(mvd_c(k), D0r)) endif !> - Autoconversion follows Berry & Reinhardt (1974) with characteristic @@ -2113,7 +2120,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & tau = 3.72/(rc(k)*taud) prr_wau(k) = zeta/tau prr_wau(k) = MIN(DBLE(rc(k)*odts), prr_wau(k)) - pnr_wau(k) = prr_wau(k) / (am_r*nu_c*D0r*D0r*D0r) ! RAIN2M + pnr_wau(k) = prr_wau(k) / (am_r*nu_c*200.*D0r*D0r*D0r) ! RAIN2M pnc_wau(k) = MIN(DBLE(nc(k)*odts), prr_wau(k) & / (am_r*mvd_c(k)*mvd_c(k)*mvd_c(k))) ! Qc2M endif @@ -2153,7 +2160,9 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ if (.not. iiwarm) then do k = kts, kte - vts_boost(k) = 1.5 + vts_boost(k) = 1.0 + xDs = 0.0 + if (L_qs(k)) xDs = smoc(k) / smob(k) !> - Temperature lookup table indexes. tempc = temp(k) - 273.15 @@ -2305,13 +2314,12 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !> - Snow collecting cloud water. In CE, assume Dc< - Snow and graupel collecting aerosols, wet scavenging. if (rs(k) .gt. r_s(1)) then - xDs = smoc(k) / smob(k) Ef_sa = Eff_aero(xDs,0.04E-6,visco(k),rho(k),temp(k),'s') pna_sca(k) = rhof(k)*t1_qs_qc*Ef_sa*nwfa(k)*smoe(k) pna_sca(k) = MIN(DBLE(nwfa(k)*odts), pna_sca(k)) @@ -2387,6 +2394,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + tnr_racs2(idx_s,idx_t,idx_r1,idx_r) & + tnr_sacr1(idx_s,idx_t,idx_r1,idx_r) & + tnr_sacr2(idx_s,idx_t,idx_r1,idx_r) + pnr_rcs(k) = MIN(DBLE(nr(k)*odts), pnr_rcs(k)) else prs_rcs(k) = -tcs_racs1(idx_s,idx_t,idx_r1,idx_r) & - tms_sacr1(idx_s,idx_t,idx_r1,idx_r) & @@ -2394,10 +2402,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & + tcr_sacr2(idx_s,idx_t,idx_r1,idx_r) prs_rcs(k) = MAX(DBLE(-rs(k)*odts), prs_rcs(k)) prr_rcs(k) = -prs_rcs(k) - pnr_rcs(k) = tnr_racs2(idx_s,idx_t,idx_r1,idx_r) & ! RAIN2M - + tnr_sacr2(idx_s,idx_t,idx_r1,idx_r) endif - pnr_rcs(k) = MIN(DBLE(nr(k)*odts), pnr_rcs(k)) endif !> - Rain collecting graupel. Cannot assume Wisner (1972) approximation @@ -2422,17 +2427,59 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & endif endif + if (temp(k).lt.T_0) then + rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999 + +!> - Deposition/sublimation of snow/graupel follows Srivastava & Coen (1992) + if (L_qs(k)) then + C_snow = C_sqrd + (tempc+1.5)*(C_cube-C_sqrd)/(-30.+1.5) + C_snow = MAX(C_sqrd, MIN(C_snow, C_cube)) + prs_sde(k) = C_snow*t1_subl*diffu(k)*ssati(k)*rvs & + * (t1_qs_sd*smo1(k) & + + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k)) + if (prs_sde(k).lt. 0.) then + prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k), DBLE(rate_max)) + else + prs_sde(k) = MIN(prs_sde(k), DBLE(rate_max)) + endif + endif + + if (L_qg(k) .and. ssati(k).lt. -eps) then + prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs & + * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) & + + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11)) + if (prg_gde(k).lt. 0.) then + prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k), DBLE(rate_max)) + else + prg_gde(k) = MIN(prg_gde(k), DBLE(rate_max)) + endif + endif + +!> - A portion of rimed snow converts to graupel but some remains snow. +!! Interp from 15 to 95% as riming factor increases from 5.0 to 30.0 +!! 0.028 came from (.75-.15)/(30.-5.). This remains ad-hoc and should +!! be revisited. + if (prs_scw(k).gt.5.0*prs_sde(k) .and. & + prs_sde(k).gt.eps) then + r_frac = MIN(30.0D0, prs_scw(k)/prs_sde(k)) + g_frac = MIN(0.75, 0.15 + (r_frac-5.)*.028) + vts_boost(k) = MIN(1.5, 1.1 + (r_frac-5.)*.016) + prg_scw(k) = g_frac*prs_scw(k) + prs_scw(k) = (1. - g_frac)*prs_scw(k) + endif + + endif + !+---+-----------------------------------------------------------------+ !> - Next IF block handles only those processes below 0C. !+---+-----------------------------------------------------------------+ if (temp(k).lt.T_0) then - vts_boost(k) = 1.0 rate_max = (qv(k)-qvsi(k))*rho(k)*odts*0.999 !+---+---------------- BEGIN NEW ICE NUCLEATION -----------------------+ -!> - Begin NEW ICE NUCLEATION: Freezing of supercooled water (rain or cloud) is influenced by dust +!> - Freezing of supercooled water (rain or cloud) is influenced by dust !! but still using Bigg 1953 with a temperature adjustment of a few !! degrees depending on dust concentration. A default value by way !! of idx_IN is 1.0 per Liter of air is used when dustyIce flag is @@ -2475,7 +2522,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pnr_rfz(k) = MIN(DBLE(nr(k)*odts), pnr_rfz(k)) elseif (rr(k).gt. R1 .and. temp(k).lt.HGFR) then pri_rfz(k) = rr(k)*odts - pnr_rfz(k) = nr(k)*odts ! RAIN2M pni_rfz(k) = pnr_rfz(k) endif @@ -2496,7 +2542,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) + xnc = xnc*(1.0 + 50.*rand3) else xnc = MIN(250.E3, TNO*EXP(ATO*(T_0-temp(k)))) endif @@ -2508,7 +2554,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !> - Freezing of aqueous aerosols based on Koop et al (2001, Nature) xni = smo0(k)+ni(k) + (pni_rfz(k)+pni_wfz(k)+pni_inu(k))*dtsave - if (is_aerosol_aware .AND. homogIce .AND. (xni.le.500.E3) & + if (is_aerosol_aware .AND. homogIce .AND. (xni.le.999.E3) & & .AND.(temp(k).lt.238).AND.(ssati(k).ge.0.4) ) then xnc = iceKoop(temp(k),qv(k),qvs(k),nwfa(k), dtsave) pni_iha(k) = xnc*odts @@ -2554,32 +2600,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & endif endif -!> - Deposition/sublimation of snow/graupel follows Srivastava & Coen -!! (1992). - if (L_qs(k)) then - C_snow = C_sqrd + (tempc+1.5)*(C_cube-C_sqrd)/(-30.+1.5) - C_snow = MAX(C_sqrd, MIN(C_snow, C_cube)) - prs_sde(k) = C_snow*t1_subl*diffu(k)*ssati(k)*rvs & - * (t1_qs_sd*smo1(k) & - + t2_qs_sd*rhof2(k)*vsc2(k)*smof(k)) - if (prs_sde(k).lt. 0.) then - prs_sde(k) = MAX(DBLE(-rs(k)*odts), prs_sde(k), DBLE(rate_max)) - else - prs_sde(k) = MIN(prs_sde(k), DBLE(rate_max)) - endif - endif - - if (L_qg(k) .and. ssati(k).lt. -eps) then - prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs & - * N0_g(k) * (t1_qg_sd*ilamg(k)**cge(10) & - + t2_qg_sd*vsc2(k)*rhof2(k)*ilamg(k)**cge(11)) - if (prg_gde(k).lt. 0.) then - prg_gde(k) = MAX(DBLE(-rg(k)*odts), prg_gde(k), DBLE(rate_max)) - else - prg_gde(k) = MIN(prg_gde(k), DBLE(rate_max)) - endif - endif - !> - Snow collecting cloud ice. In CE, assume Di< - A portion of rimed snow converts to graupel but some remains snow. -!! Interp from 15 to 95% as riming factor increases from 2.0 to 30.0 -!! 0.028 came from (.95-.15)/(30.-2.). This remains ad-hoc and should -!! be revisited. - if (prs_scw(k).gt.2.0*prs_sde(k) .and. & - prs_sde(k).gt.eps) then - r_frac = MIN(30.0D0, prs_scw(k)/prs_sde(k)) - g_frac = MIN(0.95, 0.15 + (r_frac-2.)*.028) - vts_boost(k) = MIN(1.5, 1.1 + (r_frac-2.)*.016) - prg_scw(k) = g_frac*prs_scw(k) - prs_scw(k) = (1. - g_frac)*prs_scw(k) - endif - else !> - Melt snow and graupel and enhance from collisions with liquid. @@ -2643,12 +2650,13 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & if (L_qs(k)) then prr_sml(k) = (tempc*tcond(k)-lvap0*diffu(k)*delQvs(k)) & * (t1_qs_me*smo1(k) + t2_qs_me*rhof2(k)*vsc2(k)*smof(k)) - prr_sml(k) = prr_sml(k) + 4218.*olfus*tempc & - * (prr_rcs(k)+prs_scw(k)) + if (prr_sml(k) .gt. 0.) then + prr_sml(k) = prr_sml(k) + 4218.*olfus*tempc & + * (prr_rcs(k)+prs_scw(k)) + endif prr_sml(k) = MIN(DBLE(rs(k)*odts), MAX(0.D0, prr_sml(k))) pnr_sml(k) = smo0(k)/rs(k)*prr_sml(k) * 10.0**(-0.25*tempc) ! RAIN2M pnr_sml(k) = MIN(DBLE(smo0(k)*odts), pnr_sml(k)) -! if (tempc.gt.3.5 .or. rs(k).lt.0.005E-3) pnr_sml(k)=0.0 if (ssati(k).lt. 0.) then prs_sde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs & @@ -2667,7 +2675,6 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & prr_gml(k) = MIN(DBLE(rg(k)*odts), MAX(0.D0, prr_gml(k))) pnr_gml(k) = N0_g(k)*cgg(2)*ilamg(k)**cge(2) / rg(k) & ! RAIN2M * prr_gml(k) * 10.0**(-0.5*tempc) -! if (tempc.gt.7.5 .or. rg(k).lt.0.005E-3) pnr_gml(k)=0.0 if (ssati(k).lt. 0.) then prg_gde(k) = C_cube*t1_subl*diffu(k)*ssati(k)*rvs & @@ -2677,7 +2684,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & endif endif -!> - This change will be required if users run adaptive time step that +!> - This change will be required if users run adaptive time step that !! results in delta-t that is generally too long to allow cloud water !! collection by snow/graupel above melting temperature. !! Credit to Bjorn-Egil Nygaard for discovering. @@ -2835,8 +2842,10 @@ 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 - if (rand2 .eq. 0.0) then - nu_c = MIN(15, NINT(1000.E6/xnc) + 2) + if (xnc.gt.10000.E6) then + nu_c = 2 + elseif (xnc.lt.100.) then + nu_c = 15 else nu_c = NINT(1000.E6/xnc) + 2 nu_c = MAX(2, MIN(nu_c+NINT(rand2), 15)) @@ -2881,7 +2890,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & xDi = (bm_i + mu_i + 1.) * ilami if (xDi.lt. 5.E-6) then lami = cie(2)/5.E-6 - xni = MIN(499.D3, cig(1)*oig2*xri/am_i*lami**bm_i) + xni = MIN(9999.D3, cig(1)*oig2*xri/am_i*lami**bm_i) niten(k) = (xni-ni1d(k)*rho(k))*odts*orho elseif (xDi.gt. 300.E-6) then lami = cie(2)/300.E-6 @@ -2905,7 +2914,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !> - Rain number tendency nrten(k) = nrten(k) + (pnr_wau(k) + pnr_sml(k) + pnr_gml(k) & - (pnr_rfz(k) + pnr_rcr(k) + pnr_rcg(k) & - + pnr_rcs(k) + pnr_rci(k)) ) & + + pnr_rcs(k) + pnr_rci(k) + pni_rfz(k)) ) & * orho !> - Rain mass/number balance; keep median volume diameter between @@ -2993,7 +3002,9 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & lvt2(k)=lvap(k)*lvap(k)*ocp(k)*oRv*otemp*otemp nwfa(k) = MAX(11.1E6, (nwfa1d(k) + nwfaten(k)*DT)*rho(k)) + enddo + do k = kts, kte if ((qc1d(k) + qcten(k)*DT) .gt. R1) then rc(k) = (qc1d(k) + qcten(k)*DT)*rho(k) nc(k) = MAX(2., MIN((nc1d(k)+ncten(k)*DT)*rho(k), Nt_c_max)) @@ -3118,26 +3129,11 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ !> - Calculate y-intercept, slope values for graupel. !+---+-----------------------------------------------------------------+ - N0_min = gonv_max - k_0 = kts - do k = kte, kts, -1 - if (temp(k).ge.270.65) k_0 = MAX(k_0, k) - enddo do k = kte, kts, -1 - if (k.gt.k_0 .and. L_qr(k) .and. mvd_r(k).gt.100.E-6) then - xslw1 = 4.01 + alog10(mvd_r(k)) - else - xslw1 = 0.01 - endif - ygra1 = 4.31 + alog10(max(5.E-5, rg(k))) - zans1 = (3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1))) + rand1 - if (rand1 .ne. 0.0) then - zans1 = MAX(2., MIN(zans1, 7.)) - endif + ygra1 = alog10(max(1.E-9, rg(k))) + zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1 N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) - N0_min = MIN(N0_exp, N0_min) - N0_exp = N0_min lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg ilamg(k) = 1./lamg @@ -3313,11 +3309,11 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !..TEST: G. Thompson 10 May 2013 !> - Reduce the rain evaporation in same places as melting graupel occurs. -!Rationale: falling and simultaneous melting graupel in subsaturated -!regions will not melt as fast because particle temperature stays -!..at 0C. Also not much shedding of the water from the graupel so -!..likely that the water-coated graupel evaporating much slower than -!..if the water was immediately shed off. +!! Rationale: falling and simultaneous melting graupel in subsaturated +!! regions will not melt as fast because particle temperature stays +!! at 0C. Also not much shedding of the water from the graupel so +!! likely that the water-coated graupel evaporating much slower than +!! if the water was immediately shed off. IF (prr_gml(k).gt.0.0) THEN eva_factor = MIN(1.0, 0.01+(0.99-0.01)*(tempc/20.0)) prv_rev(k) = prv_rev(k)*eva_factor @@ -3420,8 +3416,10 @@ 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 - if (rand2 .eq. 0.0) then - nu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) + if (nc(k).gt.10000.E6) then + nu_c = 2 + elseif (nc(k).lt.100.) then + nu_c = 15 else nu_c = NINT(1000.E6/nc(k)) + 2 nu_c = MAX(2, MIN(nu_c+NINT(rand2), 15)) @@ -3490,13 +3488,10 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & 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+0.1)) then - vtsk(k) = MAX(vts*vts_boost(k), & - & 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) +! vtsk(k) = MAX(vts*vts_boost(k), & +! & vts*((vtrk(k)-vts*vts_boost(k))/(temp(k)-T_0))) + SR = rs(k)/(rs(k)+rr(k)) + vtsk(k) = vts*SR + (1.-SR)*vtrk(k) else vtsk(k) = vts*vts_boost(k) endif @@ -3547,10 +3542,6 @@ 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 @@ -3580,11 +3571,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & *odzq*DT*onstep(1)) enddo -#if 1 if (rr(kts).gt.R1*10.) & -#else - if (rr(kts).gt.R1*1000.) & -#endif pptrain = pptrain + sed_r(kts)*DT*onstep(1) enddo endif @@ -3635,11 +3622,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & *odzq*DT*onstep(2)) enddo -#if 1 if (ri(kts).gt.R1*10.) & -#else - if (ri(kts).gt.R1*1000.) & -#endif pptice = pptice + sed_i(kts)*DT*onstep(2) enddo endif @@ -3666,11 +3649,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & *odzq*DT*onstep(3)) enddo -#if 1 if (rs(kts).gt.R1*10.) & -#else - if (rs(kts).gt.R1*1000.) & -#endif pptsnow = pptsnow + sed_s(kts)*DT*onstep(3) enddo endif @@ -3697,11 +3676,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & *odzq*DT*onstep(4)) enddo -#if 1 if (rg(kts).gt.R1*10.) & -#else - if (rg(kts).gt.R1*1000.) & -#endif pptgraul = pptgraul + sed_g(kts)*DT*onstep(4) enddo endif @@ -3742,19 +3717,21 @@ 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)) - nwfa1d(k) = MAX(11.1E6/rho(k), MIN(9999.E6/rho(k), & + nwfa1d(k) = MAX(11.1E6, MIN(9999.E6, & (nwfa1d(k)+nwfaten(k)*DT))) - nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6/rho(k), & + nifa1d(k) = MAX(naIN1*0.01, MIN(9999.E6, & (nifa1d(k)+nifaten(k)*DT))) if (qc1d(k) .le. R1) then qc1d(k) = 0.0 nc1d(k) = 0.0 else - if (rand2 .eq. 0.0) then - nu_c = MIN(15, NINT(1000.E6/(nc1d(k)*rho(k))) + 2) + if (nc1d(k)*rho(k).gt.10000.E6) then + nu_c = 2 + elseif (nc1d(k)*rho(k).lt.100.) then + nu_c = 15 else - nu_c = NINT(1000.E6/(nc1d(k)*rho(k))) + 2 - nu_c = MAX(2, MIN(nu_c+NINT(rand2), 15)) + 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 @@ -3782,7 +3759,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & lami = cie(2)/300.E-6 endif ni1d(k) = MIN(cig(1)*oig2*qi1d(k)/am_i*lami**bm_i, & - 499.D3/rho(k)) + 9999.D3/rho(k)) endif qr1d(k) = qr1d(k) + qrten(k)*DT nr1d(k) = MAX(R2/rho(k), nr1d(k) + nrten(k)*DT) @@ -3837,13 +3814,12 @@ subroutine qr_acr_qg good = 0 - INQUIRE(FILE="qr_acr_qg.dat",EXIST=lexist) + INQUIRE(FILE=qr_acr_qg_file, EXIST=lexist) #ifdef MPI call MPI_BARRIER(mpi_communicator,ierr) #endif IF ( lexist ) THEN - !write(0,*) "ThompMP: read qr_acr_qg.dat instead of computing" - OPEN(63,file="qr_acr_qg.dat",form="unformatted",err=1234) + OPEN(63,file=qr_acr_qg_file,form="unformatted",err=1234) !sms$serial begin READ(63,err=1234) tcg_racg READ(63,err=1234) tmr_racg @@ -3858,13 +3834,13 @@ subroutine qr_acr_qg INQUIRE(63,opened=lopen) IF (lopen) THEN IF( force_read_thompson ) THEN - write(0,*) "Error reading qr_acr_qg.dat. Aborting because force_read_thompson is .true." + write(0,*) "Error reading "//qr_acr_qg_file//" Aborting because force_read_thompson is .true." return ENDIF CLOSE(63) ELSE IF( force_read_thompson ) THEN - write(0,*) "Error opening qr_acr_qg.dat. Aborting because force_read_thompson is .true." + write(0,*) "Error opening "//qr_acr_qg_file//" Aborting because force_read_thompson is .true." return ENDIF ENDIF @@ -3876,7 +3852,7 @@ subroutine qr_acr_qg ENDIF ELSE IF( force_read_thompson ) THEN - write(0,*) "Non-existent qr_acr_qg.dat. Aborting because force_read_thompson is .true." + write(0,*) "Non-existent "//qr_acr_qg_file//" Aborting because force_read_thompson is .true." return ENDIF ENDIF @@ -3959,7 +3935,7 @@ subroutine qr_acr_qg tcg_racg(i,j,k,m) = t1 tmr_racg(i,j,k,m) = DMIN1(z1, r_r(m)*1.0d0) tcr_gacr(i,j,k,m) = t2 - tmg_gacr(i,j,k,m) = z2 + tmg_gacr(i,j,k,m) = DMIN1(z2, r_g(j)*1.0d0) tnr_racg(i,j,k,m) = y1 tnr_gacr(i,j,k,m) = y2 enddo @@ -3967,8 +3943,8 @@ subroutine qr_acr_qg enddo IF ( write_thompson_tables ) THEN - write(0,*) "Writing qr_acr_qg.dat in Thompson MP init" - OPEN(63,file="qr_acr_qg.dat",form="unformatted",err=9234) + write(0,*) "Writing "//qr_acr_qg_file//" in Thompson MP init" + OPEN(63,file=qr_acr_qg_file,form="unformatted",err=9234) WRITE(63,err=9234) tcg_racg WRITE(63,err=9234) tmr_racg WRITE(63,err=9234) tcr_gacr @@ -3978,7 +3954,7 @@ subroutine qr_acr_qg CLOSE(63) RETURN ! ----- RETURN 9234 CONTINUE - write(0,*) "Error writing qr_acr_qg.dat" + write(0,*) "Error writing //qr_acr_qg_file return ENDIF ENDIF @@ -4013,13 +3989,13 @@ subroutine qr_acr_qs write_thompson_tables = .false. good = 0 - INQUIRE(FILE="qr_acr_qs.dat",EXIST=lexist) + INQUIRE(FILE=qr_acr_qs_file, EXIST=lexist) #ifdef MPI call MPI_BARRIER(mpi_communicator,ierr) #endif IF ( lexist ) THEN - !write(0,*) "ThompMP: read qr_acr_qs.dat instead of computing" - OPEN(63,file="qr_acr_qs.dat",form="unformatted",err=1234) + !write(0,*) "ThompMP: read "//qr_acr_qs_file//" instead of computing" + OPEN(63,file=qr_acr_qs_file,form="unformatted",err=1234) !sms$serial begin READ(63,err=1234)tcs_racs1 READ(63,err=1234)tmr_racs1 @@ -4040,13 +4016,13 @@ subroutine qr_acr_qs INQUIRE(63,opened=lopen) IF (lopen) THEN IF( force_read_thompson ) THEN - write(0,*) "Error reading qr_acr_qs.dat. Aborting because force_read_thompson is .true." + write(0,*) "Error reading "//qr_acr_qs_file//" Aborting because force_read_thompson is .true." return ENDIF CLOSE(63) ELSE IF( force_read_thompson ) THEN - write(0,*) "Error opening qr_acr_qs.dat. Aborting because force_read_thompson is .true." + write(0,*) "Error opening "//qr_acr_qs_file//" Aborting because force_read_thompson is .true." return ENDIF ENDIF @@ -4058,7 +4034,7 @@ subroutine qr_acr_qs ENDIF ELSE IF( force_read_thompson ) THEN - write(0,*) "Non-existent qr_acr_qs.dat. Aborting because force_read_thompson is .true." + write(0,*) "Non-existent "//qr_acr_qs_file//" Aborting because force_read_thompson is .true." return ENDIF ENDIF @@ -4219,8 +4195,8 @@ subroutine qr_acr_qs enddo IF ( write_thompson_tables ) THEN - write(0,*) "Writing qr_acr_qs.dat in Thompson MP init" - OPEN(63,file="qr_acr_qs.dat",form="unformatted",err=9234) + write(0,*) "Writing "//qr_acr_qs_file//" in Thompson MP init" + OPEN(63,file=qr_acr_qs_file,form="unformatted",err=9234) WRITE(63,err=9234)tcs_racs1 WRITE(63,err=9234)tmr_racs1 WRITE(63,err=9234)tcs_racs2 @@ -4236,7 +4212,7 @@ subroutine qr_acr_qs CLOSE(63) RETURN ! ----- RETURN 9234 CONTINUE - write(0,*) "Error writing qr_acr_qs.dat" + write(0,*) "Error writing "//qr_acr_qs_file ENDIF ENDIF @@ -4274,13 +4250,13 @@ subroutine freezeH2O(threads) write_thompson_tables = .false. good = 0 - INQUIRE(FILE="freezeH2O.dat",EXIST=lexist) + INQUIRE(FILE=freeze_h2o_file",EXIST=lexist) #ifdef MPI call MPI_BARRIER(mpi_communicator,ierr) #endif IF ( lexist ) THEN - !write(0,*) "ThompMP: read freezeH2O.dat instead of computing" - OPEN(63,file="freezeH2O.dat",form="unformatted",err=1234) + !write(0,*) "ThompMP: read "//freeze_h2o_file//" instead of computing" + OPEN(63,file=freeze_h2o_file,form="unformatted",err=1234) !sms$serial begin READ(63,err=1234)tpi_qrfz READ(63,err=1234)tni_qrfz @@ -4295,13 +4271,13 @@ subroutine freezeH2O(threads) INQUIRE(63,opened=lopen) IF (lopen) THEN IF( force_read_thompson ) THEN - write(0,*) "Error reading freezeH2O.dat. Aborting because force_read_thompson is .true." + write(0,*) "Error reading "//freeze_h2o_file//" Aborting because force_read_thompson is .true." return ENDIF CLOSE(63) ELSE IF( force_read_thompson ) THEN - write(0,*) "Error opening freezeH2O.dat. Aborting because force_read_thompson is .true." + write(0,*) "Error opening "//freeze_h2o_file//" Aborting because force_read_thompson is .true." return ENDIF ENDIF @@ -4313,7 +4289,7 @@ subroutine freezeH2O(threads) ENDIF ELSE IF( force_read_thompson ) THEN - write(0,*) "Non-existent freezeH2O.dat. Aborting because force_read_thompson is .true." + write(0,*) "Non-existent "//freeze_h2o_file//" Aborting because force_read_thompson is .true." return ENDIF ENDIF @@ -4397,8 +4373,8 @@ subroutine freezeH2O(threads) enddo IF ( write_thompson_tables ) THEN - write(0,*) "Writing freezeH2O.dat in Thompson MP init" - OPEN(63,file="freezeH2O.dat",form="unformatted",err=9234) + write(0,*) "Writing "//freeze_h2o_file//" in Thompson MP init" + OPEN(63,file=freeze_h2o_file,form="unformatted",err=9234) WRITE(63,err=9234)tpi_qrfz WRITE(63,err=9234)tni_qrfz WRITE(63,err=9234)tpg_qrfz @@ -4408,7 +4384,7 @@ subroutine freezeH2O(threads) CLOSE(63) RETURN ! ----- RETURN 9234 CONTINUE - write(0,*) "Error writing freezeH2O.dat" + write(0,*) "Error writing "//freeze_h2o_file return ENDIF ENDIF @@ -5120,7 +5096,7 @@ real function iceDeMott(tempc, qv, qvs, qvsi, rho, nifa) ! mux = hx*p_alpha*n_in*rho ! xni = mux*((6700.*nifa)-200.)/((6700.*5.E5)-200.) ! elseif (satw.ge.0.985 .and. tempc.gt.HGFR-273.15) then - nifa_cc = nifa*RHO_NOT0*1.E-6/rho + nifa_cc = MAX(0.5, nifa*RHO_NOT0*1.E-6/rho) ! xni = 3.*nifa_cc**(1.25)*exp((0.46*(-tempc))-11.6) ! [DeMott, 2015] xni = (5.94e-5*(-tempc)**3.33) & ! [DeMott, 2010] * (nifa_cc**((-0.0264*(tempc))+0.0033)) @@ -5233,23 +5209,6 @@ 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. - re_qc1d(:) = 2.50E-6 ! 2.49E-6 - re_qi1d(:) = 5.00E-6 ! 4.99E-6 - re_qs1d(:) = 1.00E-5 ! 9.99E-6 - 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)) @@ -5274,7 +5233,7 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & inu_c = MIN(15, NINT(1000.E6/nc(k)) + 2) endif lamc = (nc(k)*am_r*g_ratio(inu_c)/rc(k))**obmr - re_qc1d(k) = MAX(2.51E-6, MIN(SNGL(0.5D0 * DBLE(3.+inu_c)/lamc), 50.E-6)) + re_qc1d(k) = SNGL(0.5D0 * DBLE(3.+inu_c)/lamc) enddo endif @@ -5282,7 +5241,7 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & do k = kts, kte if (ri(k).le.R1 .or. ni(k).le.R2) CYCLE lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi - re_qi1d(k) = MAX(5.01E-6, MIN(SNGL(0.5D0 * DBLE(3.+mu_i)/lami), 125.E-6)) + re_qi1d(k) = SNGL(0.5D0 * DBLE(3.+mu_i)/lami) enddo endif @@ -5322,7 +5281,7 @@ 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_ - re_qs1d(k) = MAX(1.01E-5, MIN(0.5*(smoc/smob), 999.E-6)) + re_qs1d(k) = 0.5*(smoc/smob) enddo endif @@ -5441,8 +5400,15 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & !+---+-----------------------------------------------------------------+ !..Calculate y-intercept, slope, and useful moments for snow. !+---+-----------------------------------------------------------------+ + do k = kts, kte + smo2(k) = 0. + smob(k) = 0. + smoc(k) = 0. + smoz(k) = 0. + enddo if (ANY(L_qs .eqv. .true.)) then do k = kts, kte + if (.not. L_qs(k)) CYCLE tc0 = MIN(-0.1, temp(k)-273.15) smob(k) = rs(k)*oams @@ -5498,26 +5464,11 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & !+---+-----------------------------------------------------------------+ if (ANY(L_qg .eqv. .true.)) then - N0_min = gonv_max - k_0 = kts - do k = kte, kts, -1 - if (temp(k).ge.270.65) k_0 = MAX(k_0, k) - enddo do k = kte, kts, -1 - if (k.gt.k_0 .and. L_qr(k) .and. mvd_r(k).gt.100.E-6) then - xslw1 = 4.01 + alog10(mvd_r(k)) - else - xslw1 = 0.01 - endif - ygra1 = 4.31 + alog10(max(5.E-5, rg(k))) - zans1 = (3.1 + (100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1))) + rand1 - if (rand1 .ne. 0.0) then - zans1 = MAX(2., MIN(zans1, 7.)) - endif + ygra1 = alog10(max(1.E-9, rg(k))) + zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1 N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) - N0_min = MIN(N0_exp, N0_min) - N0_exp = N0_min lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg ilamg(k) = 1./lamg @@ -5666,7 +5617,7 @@ end subroutine calc_refl10cm #ifdef SION !>\ingroup aathompson - subroutine readwrite_tables(mode, mpicomm, mpirank, mpiroot, ierr) + subroutine readwrite_tables(filename, mode, mpicomm, mpirank, mpiroot, ierr) #ifdef MPI use mpi @@ -5676,6 +5627,7 @@ subroutine readwrite_tables(mode, mpicomm, mpirank, mpiroot, ierr) implicit none ! Interface variables + character(len=*), intent(in) :: filename character(len=*), intent(in) :: mode integer, intent(in) :: mpicomm integer, intent(in) :: mpirank @@ -5702,7 +5654,6 @@ subroutine readwrite_tables(mode, mpicomm, mpirank, mpiroot, ierr) logical :: exists integer*8 :: tables_size real*8 :: checksum - character(len=*), parameter :: filename = 'thompson_tables_precomp.sl' integer :: i From 68113c9989ddaa203a7c84a1dabbc58adcc86f39 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Sat, 13 Feb 2021 11:50:04 -0700 Subject: [PATCH 02/31] simplify initialization a bit and add convert_dry_rho flag --- physics/module_mp_thompson.F90 | 17 +-- physics/mp_thompson.F90 | 245 +++++++++++---------------------- 2 files changed, 86 insertions(+), 176 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index fcbcb8164..444e867ae 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -219,7 +219,7 @@ MODULE module_mp_thompson REAL, PRIVATE:: D0i, xm0s, xm0g !..Min and max radiative effective radius of cloud water, cloud ice, and snow; -!.. performed by subroutine calc_effectRad +!.. performed by subroutine calc_effectRad. On purpose, these should stay PUBLIC. REAL, PARAMETER:: re_qc_min = 2.50E-6 ! 2.5 microns REAL, PARAMETER:: re_qc_max = 50.0E-6 ! 50 microns REAL, PARAMETER:: re_qi_min = 2.50E-6 ! 2.5 microns @@ -440,24 +440,17 @@ MODULE module_mp_thompson !! lookup tables in Thomspson scheme. !>\section gen_thompson_init thompson_init General Algorithm !> @{ - SUBROUTINE thompson_init(nwfa2d, nifa2d, nwfa, nifa, & - mpicomm, mpirank, mpiroot, & + SUBROUTINE thompson_init(mpicomm, mpirank, mpiroot, & threads, errmsg, errflg) IMPLICIT NONE -!..OPTIONAL variables that control application of aerosol-aware scheme - - REAL, DIMENSION(:,:), OPTIONAL, INTENT(IN) :: nwfa, nifa - REAL, DIMENSION(:), OPTIONAL, INTENT(IN) :: nwfa2d, nifa2d INTEGER, INTENT(IN) :: mpicomm, mpirank, mpiroot INTEGER, INTENT(IN) :: threads CHARACTER(len=*), INTENT(INOUT) :: errmsg INTEGER, INTENT(INOUT) :: errflg - INTEGER:: i, j, k, l, m, n - REAL:: h_01, airmass, niIN3, niCCN3, max_test LOGICAL:: micro_init real :: stime, etime #ifdef SION @@ -467,14 +460,8 @@ SUBROUTINE thompson_init(nwfa2d, nifa2d, nwfa, nifa, & LOGICAL, PARAMETER :: precomputed_tables = .FALSE. #endif - is_aerosol_aware = .FALSE. micro_init = .FALSE. - if (present(nwfa2d) .and. & - present(nifa2d) .and. & - present(nwfa) .and. & - present(nifa) ) is_aerosol_aware = .true. - !> - Allocate space for lookup tables (J. Michalakes 2009Jun08). if (.NOT. ALLOCATED(tcg_racg) ) then diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index ec19945b0..e22d82316 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -10,6 +10,7 @@ module mp_thompson use module_mp_thompson, only : thompson_init, mp_gt_driver, thompson_finalize, calc_effectRad use module_mp_thompson, only : naIN0, naIN1, naCCN0, naCCN1, eps, Nt_c + use module_mp_thompson, only : re_qc_min, re_qc_max, re_qi_min, re_qi_max, re_qs_min, re_qs_max use module_mp_thompson_make_number_concentrations, only: make_IceNumber, make_DropletNumber, make_RainNumber @@ -20,6 +21,7 @@ module mp_thompson private logical :: is_initialized = .False. + logical :: convert_dry_rho = .False. contains @@ -80,17 +82,8 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & character(len=*), intent( out) :: errmsg integer, intent( out) :: errflg - ! Hydrometeors - real(kind_phys) :: qv_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) - real(kind_phys) :: qc_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) - real(kind_phys) :: qr_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) - real(kind_phys) :: qi_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) - real(kind_phys) :: qs_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) - real(kind_phys) :: qg_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) - real(kind_phys) :: ni_mp(1:ncol,1:nlev) !< kg-1 - real(kind_phys) :: nr_mp(1:ncol,1:nlev) !< kg-1 - real(kind_phys) :: nc_mp(1:ncol,1:nlev) !< kg-1 ! + real(kind_phys) :: qv(1:ncol,1:nlev) ! kg kg-1 (water vapor mixing ratio) real(kind_phys) :: hgt(1:ncol,1:nlev) ! m real(kind_phys) :: rho(1:ncol,1:nlev) ! kg m-3 real(kind_phys) :: orho(1:ncol,1:nlev) ! m3 kg-1 @@ -120,16 +113,9 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & end if ! Call Thompson init - if (is_aerosol_aware) then - call thompson_init(nwfa2d=nwfa2d, nifa2d=nifa2d, nwfa=nwfa, nifa=nifa, & - mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & - threads=threads, errmsg=errmsg, errflg=errflg) - if (errflg /= 0) return - else - call thompson_init(mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & - threads=threads, errmsg=errmsg, errflg=errflg) - if (errflg /= 0) return - end if + call thompson_init(mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & + threads=threads, errmsg=errmsg, errflg=errflg) + if (errflg /= 0) return ! For restart runs, the init is done here if (restart) then @@ -137,25 +123,6 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & return end if - ! Fix initial values of hydrometeors - where(spechum<0) spechum = 0.0 - where(qc<0) qc = 0.0 - where(qr<0) qr = 0.0 - where(qi<0) qi = 0.0 - where(qs<0) qs = 0.0 - where(qg<0) qg = 0.0 - where(ni<0) ni = 0.0 - where(nr<0) nr = 0.0 - - if (is_aerosol_aware) then - ! Fix initial values of aerosols - where(nc<0) nc = 0.0 - where(nwfa<0) nwfa = 0.0 - where(nifa<0) nifa = 0.0 - where(nwfa2d<0) nwfa2d = 0.0 - where(nifa2d<0) nifa2d = 0.0 - end if - ! Geopotential height in m2 s-2 to height in m hgt = phil/con_g @@ -163,49 +130,33 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & rho = prsl/(con_rd*tgrs) orho = 1.0/rho - ! Prior to calling the functions: make_DropletNumber, make_IceNumber, make_RainNumber, - ! the incoming mixing ratios should be converted to units of mass/num per cubic meter - ! rather than per kg of air. So, to pass back to the model state variables, - ! they also need to be switched back to mass/number per kg of air, because - ! what is returned by the functions is in units of number per cubic meter. - ! They also need to be converted to dry mixing ratios. - - !> - Convert specific humidity/moist mixing ratios to dry mixing ratios - qv_mp = spechum/(1.0_kind_phys-spechum) - qc_mp = qc/(1.0_kind_phys-spechum) - qr_mp = qr/(1.0_kind_phys-spechum) - qi_mp = qi/(1.0_kind_phys-spechum) - qs_mp = qs/(1.0_kind_phys-spechum) - qg_mp = qg/(1.0_kind_phys-spechum) - - !> - Convert number concentrations from moist to dry - ni_mp = ni/(1.0_kind_phys-spechum) - nr_mp = nr/(1.0_kind_phys-spechum) - if (is_aerosol_aware) then - nc_mp = nc/(1.0_kind_phys-spechum) - end if + ! Ensure non-negative mass mixing ratios of all water variables + where(spechum<0) spechum = 1.0E-10 ! COMMENT, gthompsn, spechum should *never* be identically zero. + where(qc<0) qc = 0.0 + where(qr<0) qr = 0.0 + where(qi<0) qi = 0.0 + where(qs<0) qs = 0.0 + where(qg<0) qg = 0.0 - ! If qi is in boundary conditions but ni is not, calculate ni from qi, rho and tgrs - if (maxval(qi_mp)>0.0 .and. maxval(ni_mp)==0.0) then - ni_mp = make_IceNumber(qi_mp*rho, tgrs) * orho - end if + ! Convert specific humidity to water vapor mixing ratio + qv = spechum/(1.0_kind_phys-spechum) - ! If ni is in boundary conditions but qi is not, reset ni to zero - if (maxval(ni_mp)>0.0 .and. maxval(qi_mp)==0.0) ni_mp = 0.0 + ! Ensure we have 1st guess ice number where mass non-zero but no number. + where(qi <= 0.0) ni=0.0 + where(qi > 0 .and. ni <= 0.0) ni = make_IceNumber(qi*rho, tgrs) * orho + where(qi = 0.0 .and. ni > 0.0) ni=0.0 - ! If qr is in boundary conditions but nr is not, calculate nr from qr, rho and tgrs - if (maxval(qr_mp)>0.0 .and. maxval(nr_mp)==0.0) then - nr_mp = make_RainNumber(qr_mp*rho, tgrs) * orho - end if + ! Ensure we have 1st guess rain number where mass non-zero but no number. + where(qr <= 0.0) nr=0.0 + where(qr > 0 .and. nr <= 0.0) nr = make_RainNumber(qr*rho, tgrs) * orho + where(qr = 0.0 .and. nr > 0.0) nr=0.0 - ! If nr is in boundary conditions but qr is not, reset nr to zero - if (maxval(nr_mp)>0.0 .and. maxval(qr_mp)==0.0) nr_mp = 0.0 !..Check for existing aerosol data, both CCN and IN aerosols. If missing !.. fill in just a basic vertical profile, somewhat boundary-layer following. if (is_aerosol_aware) then - ! CCN + ! Potential cloud condensation nuclei (CCN) if (MAXVAL(nwfa) .lt. eps) then if (mpirank==mpiroot) write(*,*) ' Apparently there are no initial CCN aerosols.' do i = 1, ncol @@ -219,7 +170,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & niCCN3 = -1.0*ALOG(naCCN1/naCCN0)/h_01 nwfa(i,1) = naCCN1+naCCN0*exp(-((hgt(i,2)-hgt(i,1))/1000.)*niCCN3) airmass = 1./orho(i,1) * (hgt(i,2)-hgt(i,1))*area(i) ! kg - nwfa2d(i) = nwfa(i,1) * 0.000196 * (airmass*2.E-10) + nwfa2d(i) = nwfa(i,1) * 0.000196 * (airmass*5.E-11) do k = 2, nlev nwfa(i,k) = naCCN1+naCCN0*exp(-((hgt(i,k)-hgt(i,1))/1000.)*niCCN3) enddo @@ -227,8 +178,6 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & else if (mpirank==mpiroot) write(*,*) ' Apparently initial CCN aerosols are present.' if (MAXVAL(nwfa2d) .lt. eps) then -! Hard-coded switch between new (from WRFv4.0, top) and old (until WRFv3.9.1.1, bottom) surface emission rate calculations -#if 0 !+---+-----------------------------------------------------------------+ !..Scale the lowest level aerosol data into an emissions rate. This is !.. very far from ideal, but need higher emissions where larger amount @@ -239,41 +188,16 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & !.. that was tested as ~(20kmx20kmx50m = 2.E10 m**-3) !+---+-----------------------------------------------------------------+ if (mpirank==mpiroot) write(*,*) ' Apparently there are no initial CCN aerosol surface emission rates.' - if (mpirank==mpiroot) write(*,*) ' Use new (WRFv4+) formula to calculate CCN surface emission rates.' do i = 1, ncol airmass = 1./orho(i,1) * (hgt(i,2)-hgt(i,1))*area(i) ! kg - nwfa2d(i) = nwfa(i,1) * 0.000196 * (airmass*2.E-10) - enddo -#else - !+---+-----------------------------------------------------------------+ - !..Scale the lowest level aerosol data into an emissions rate. This is - !.. very far from ideal, but need higher emissions where larger amount - !.. of existing and lesser emissions where not already lots of aerosols - !.. for first-order simplistic approach. Later, proper connection to - !.. emission inventory would be better, but, for now, scale like this: - !.. where: Nwfa=50 per cc, emit 0.875E4 aerosols per kg per second - !.. Nwfa=500 per cc, emit 0.875E5 aerosols per kg per second - !.. Nwfa=5000 per cc, emit 0.875E6 aerosols per kg per second - !.. for a grid with 20km spacing and scale accordingly for other spacings. - !+---+-----------------------------------------------------------------+ - if (mpirank==mpiroot) write(*,*) ' Apparently there are no initial CCN aerosol surface emission rates.' - if (mpirank==mpiroot) write(*,*) ' Use old (pre WRFv4) formula to calculate CCN surface emission rates.' - do i = 1, ncol - if (SQRT(area(i))/20000.0 .ge. 1.0) then - h_01 = 0.875 - else - h_01 = (0.875 + 0.125*((20000.-SQRT(area(i)))/16000.)) * SQRT(area(i))/20000. - endif - nwfa2d(i) = 10.0**(LOG10(nwfa(i,1)*1.E-6)-3.69897) - nwfa2d(i) = nwfa2d(i)*h_01 * 1.E6 + nwfa2d(i) = nwfa(i,1) * 0.000196 * (airmass*5.E-11) enddo -#endif else if (mpirank==mpiroot) write(*,*) ' Apparently initial CCN aerosol surface emission rates are present.' endif endif - ! IN + ! Potential ice nuclei (IN) if (MAXVAL(nifa) .lt. eps) then if (mpirank==mpiroot) write(*,*) ' Apparently there are no initial IN aerosols.' do i = 1, ncol @@ -302,19 +226,20 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & endif endif - ! If qc is in boundary conditions but nc is not, calculate nc from qc, rho and nwfa - if (maxval(qc_mp)>0.0 .and. maxval(nc_mp)==0.0) then - nc_mp = make_DropletNumber(qc_mp*rho, nwfa) * orho - end if + ! Ensure we have 1st guess cloud droplet number where mass non-zero but no number. + where(qc <= 0.0) nc=0.0 + where(qc > 0 .and. nc <= 0.0) nc = make_DropletNumber(qc*rho, nwfa) * orho + where(qc = 0.0 .and. nc > 0.0) nc=0.0 - ! If nc is in boundary conditions but qc is not, reset nc to zero - if (maxval(nc_mp)>0.0 .and. maxval(qc_mp)==0.0) nc_mp = 0.0 + ! Ensure non-negative aerosol number concentrations. + where(nwfa <= 0.0) nwfa = 1.1E6 + where(nifa <= 0.0) nifa = naIN1*0.01 else ! Constant droplet concentration for single moment cloud water as in ! module_mp_thompson.F90, only needed for effective radii calculation - nc_mp = Nt_c/rho + nc = Nt_c/rho end if @@ -322,9 +247,14 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & if (present(re_cloud) .and. present(re_ice) .and. present(re_snow)) then ! 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,:), & + call calc_effectRad (tgrs(i,:), prsl(i,:), qv(i,:), qc(i,:), & + nc(i,:), qi(i,:), ni(i,:), qs(i,:), & re_cloud(i,:), re_ice(i,:), re_snow(i,:), 1, nlev) + do k = 1, nlev + re_cloud(i,k) = MAX(re_qc_min, MIN(re_cloud(i,k), re_qc_max)) + re_ice(i,k) = MAX(re_qi_min, MIN(re_ice(i,k), re_qi_max)) + re_snow(i,k) = MAX(re_qs_min, MIN(re_snow(i,k), re_qs_max)) + 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 @@ -341,13 +271,6 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & return end if - !> - Convert number concentrations from dry to moist - ni = ni_mp/(1.0_kind_phys+qv_mp) - nr = nr_mp/(1.0_kind_phys+qv_mp) - if (is_aerosol_aware) then - nc = nc_mp/(1.0_kind_phys+qv_mp) - end if - is_initialized = .true. end subroutine mp_thompson_init @@ -428,17 +351,8 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ! Air density real(kind_phys) :: rho(1:ncol,1:nlev) !< kg m-3 - ! Hydrometeors - real(kind_phys) :: qv_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) - real(kind_phys) :: qc_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) - real(kind_phys) :: qr_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) - real(kind_phys) :: qi_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) - real(kind_phys) :: qs_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) - real(kind_phys) :: qg_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) - real(kind_phys) :: ni_mp(1:ncol,1:nlev) !< kg-1 - real(kind_phys) :: nr_mp(1:ncol,1:nlev) !< kg-1 - real(kind_phys) :: nc_mp(1:ncol,1:nlev) !< kg-1 - + ! Water vapor mixing ratio (instead of specific humidity) + real(kind_phys) :: qv(1:ncol,1:nlev) !< kg kg-1 ! Vertical velocity and level width real(kind_phys) :: w(1:ncol,1:nlev) !< m s-1 real(kind_phys) :: dz(1:ncol,1:nlev) !< m @@ -494,19 +408,26 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & return end if - !> - Convert specific humidity/moist mixing ratios to dry mixing ratios - qv_mp = spechum/(1.0_kind_phys-spechum) - qc_mp = qc/(1.0_kind_phys-spechum) - qr_mp = qr/(1.0_kind_phys-spechum) - qi_mp = qi/(1.0_kind_phys-spechum) - qs_mp = qs/(1.0_kind_phys-spechum) - qg_mp = qg/(1.0_kind_phys-spechum) - - !> - Convert number concentrations from moist to dry - ni_mp = ni/(1.0_kind_phys-spechum) - nr_mp = nr/(1.0_kind_phys-spechum) - if (is_aerosol_aware) then - nc_mp = nc/(1.0_kind_phys-spechum) + !> - Convert specific humidity to water vapor mixing ratio. + !> - Also, hydrometeor variables are mass or number mixing ratio + !> - either kg of species per kg of dry air, or per kg of (dry + vapor). + + qv = spechum/(1.0_kind_phys-spechum) + + if (convert_dry_rho) then + qc = qc/(1.0_kind_phys-spechum) + qr = qr/(1.0_kind_phys-spechum) + qi = qi/(1.0_kind_phys-spechum) + qs = qs/(1.0_kind_phys-spechum) + qg = qg/(1.0_kind_phys-spechum) + + ni = ni/(1.0_kind_phys-spechum) + nr = nr/(1.0_kind_phys-spechum) + if (is_aerosol_aware) then + nc = nc/(1.0_kind_phys-spechum) + nwfa = nwfa/(1.0_kind_phys-spechum) + nifa = nifa/(1.0_kind_phys-spechum) + end if end if !> - Density of air in kg m-3 @@ -582,9 +503,8 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & !> - Call mp_gt_driver() with or without aerosols if (is_aerosol_aware) then - call mp_gt_driver(qv=qv_mp, qc=qc_mp, qr=qr_mp, qi=qi_mp, qs=qs_mp, qg=qg_mp, & - ni=ni_mp, nr=nr_mp, nc=nc_mp, & - nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & + call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & + nc=nc, nwfa=nwfa, nifa=nifa, nwfa2d=nwfa2d, nifa2d=nifa2d, & tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & @@ -604,8 +524,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & errmsg=errmsg, errflg=errflg, reset=reset) else - call mp_gt_driver(qv=qv_mp, qc=qc_mp, qr=qr_mp, qi=qi_mp, qs=qs_mp, qg=qg_mp, & - ni=ni_mp, nr=nr_mp, & + call mp_gt_driver(qv=qv, qc=qc, qr=qr, qi=qi, qs=qs, qg=qg, ni=ni, nr=nr, & tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtp, & rainnc=rain_mp, rainncv=delta_rain_mp, & snownc=snow_mp, snowncv=delta_snow_mp, & @@ -626,19 +545,23 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if if (errflg/=0) return - !> - Convert dry mixing ratios to specific humidity/moist mixing ratios - spechum = qv_mp/(1.0_kind_phys+qv_mp) - qc = qc_mp/(1.0_kind_phys+qv_mp) - qr = qr_mp/(1.0_kind_phys+qv_mp) - qi = qi_mp/(1.0_kind_phys+qv_mp) - qs = qs_mp/(1.0_kind_phys+qv_mp) - qg = qg_mp/(1.0_kind_phys+qv_mp) - - !> - Convert number concentrations from dry to moist - ni = ni_mp/(1.0_kind_phys+qv_mp) - nr = nr_mp/(1.0_kind_phys+qv_mp) - if (is_aerosol_aware) then - nc = nc_mp/(1.0_kind_phys+qv_mp) + !> - Convert water vapor mixing ratio back to specific humidity + spechum = qv/(1.0_kind_phys+qv) + + if (convert_dry_rho) then + qc = qc/(1.0_kind_phys+qv) + qr = qr/(1.0_kind_phys+qv) + qi = qi/(1.0_kind_phys+qv) + qs = qs/(1.0_kind_phys+qv) + qg = qg/(1.0_kind_phys+qv) + + ni = ni/(1.0_kind_phys+qv) + nr = nr/(1.0_kind_phys+qv) + if (is_aerosol_aware) then + nc = nc/(1.0_kind_phys+qv) + nwfa = nwfa/(1.0_kind_phys+qv) + nifa = nifa/(1.0_kind_phys+qv) + end if end if !> - Convert rainfall deltas from mm to m (on physics timestep); add to inout variables From b204497750d1fad78f9752eff8a7dabc60eb8a83 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Wed, 17 Feb 2021 10:05:14 -0700 Subject: [PATCH 03/31] fix bug to include rho(air) in nwfa --- physics/mp_thompson.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index e22d82316..ac8437262 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -228,7 +228,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & ! Ensure we have 1st guess cloud droplet number where mass non-zero but no number. where(qc <= 0.0) nc=0.0 - where(qc > 0 .and. nc <= 0.0) nc = make_DropletNumber(qc*rho, nwfa) * orho + where(qc > 0 .and. nc <= 0.0) nc = make_DropletNumber(qc*rho, nwfa*rho) * orho where(qc = 0.0 .and. nc > 0.0) nc=0.0 ! Ensure non-negative aerosol number concentrations. From be017d02cc0fc17f9c0a302693076149c9f21c8c Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 22 Feb 2021 10:28:20 -0700 Subject: [PATCH 04/31] Revert change to CODEOWNERS --- CODEOWNERS | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CODEOWNERS b/CODEOWNERS index 0d5230f89..b6c597371 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -3,7 +3,7 @@ # These owners will be the default owners for everything in the repo. #* @defunkt -* @climbfuji @llpcarson @grantfirl @JulieSchramm +* @DomHeinzeller # Order is important. The last matching pattern has the most precedence. # So if a pull request only touches javascript files, only these owners From 5f7968b2628b8ab8838d9e9d7cdb4a498f95f009 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Tue, 23 Feb 2021 01:51:19 +0000 Subject: [PATCH 05/31] 1. Fixed units inconsistency for water friendly aerosols in the computation of number concentration of liquid water. 2. Replaced dry-air density with the moist-air density for consistency with Thompson MP. --- physics/GFS_rrtmg_pre.F90 | 10 ++++++---- physics/GFS_suite_interstitial.F90 | 10 +++++----- physics/mp_thompson.F90 | 13 +++++++------ 3 files changed, 18 insertions(+), 15 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index 109df3b65..ab488c687 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -288,8 +288,6 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & plyr(i,k1) = prsl(i,k2) * 0.01 ! pa to mb (hpa) tlyr(i,k1) = tgrs(i,k2) prslk1(i,k1) = prslk(i,k2) - rho(i,k1) = plyr(i,k1)/(con_rd*tlyr(i,k1)) - orho(i,k1) = 1.0/rho(i,k1) !> - Compute relative humidity. es = min( prsl(i,k2), fpvs( tgrs(i,k2) ) ) ! fpvs and prsl in pa @@ -636,6 +634,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & do i=1,IM qvs = qgrs(i,k,ntqv) qv_mp (i,k) = qvs/(1.-qvs) + rho (i,k) = 0.622*prsl(i,k)/(con_rd*tgrs(i,k)*(qv_mp(i,k)+0.622)) + orho (i,k) = 1.0/rho(i,k) qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs) qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs) qs_mp (i,k) = tracer1(i,k,ntsw)/(1.-qvs) @@ -649,6 +649,8 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & do i=1,IM qvs = qgrs(i,k,ntqv) qv_mp (i,k) = qvs/(1.-qvs) + rho (i,k) = 0.622*prsl(i,k)/(con_rd*tgrs(i,k)*(qv_mp(i,k)+0.622)) + orho (i,k) = 1.0/rho(i,k) qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs) qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs) qs_mp (i,k) = tracer1(i,k,ntsw)/(1.-qvs) @@ -761,7 +763,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & do k=1,lm do i=1,im if (ltaerosol .and. qc_mp(i,k)>1.e-12 .and. nc_mp(i,k)<100.) then - nc_mp(i,k) = make_DropletNumber(qc_mp(i,k)*rho(i,k), nwfa(i,k)) * orho(i,k) + nc_mp(i,k) = make_DropletNumber(qc_mp(i,k)*rho(i,k), nwfa(i,k)*rho(i,k)) * orho(i,k) endif if (qi_mp(i,k)>1.e-12 .and. ni_mp(i,k)<100.) then ni_mp(i,k) = make_IceNumber(qi_mp(i,k)*rho(i,k), tlyr(i,k)) * orho(i,k) @@ -774,7 +776,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & !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,:), & + call calc_effectRad (tlyr(i,:), plyr(i,:)*100., 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 diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index 551f0e600..f9a26bb57 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -717,7 +717,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to ! local variables integer :: i,k,n,tracers - real(kind=kind_phys), dimension(im,levs) :: rho_dryair + real(kind=kind_phys), dimension(im,levs) :: rho_air real(kind=kind_phys), dimension(im,levs) :: qv_mp !< kg kg-1 (dry mixing ratio) real(kind=kind_phys), dimension(im,levs) :: qc_mp !< kg kg-1 (dry mixing ratio) real(kind=kind_phys), dimension(im,levs) :: qi_mp !< kg kg-1 (dry mixing ratio) @@ -767,16 +767,16 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to if (imp_physics == imp_physics_thompson .and. (ntlnc>0 .or. ntinc>0)) then do k=1,levs do i=1,im - !> - Density of air in kg m-3 - rho_dryair(i,k) = prsl(i,k) / (con_rd*save_tcp(i,k)) !> - Convert specific humidity to dry mixing ratio qv_mp(i,k) = spechum(i,k) / (one-spechum(i,k)) + !> - Density of air in kg m-3 + rho_air(i,k) = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+0.622)) if (ntlnc>0) then !> - Convert moist mixing ratio to dry mixing ratio qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) / (one-spechum(i,k)) !> - Convert number concentration from moist to dry nc_mp(i,k) = gq0(i,k,ntlnc) / (one-spechum(i,k)) - nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho_dryair(i,k), nwfa(i,k)) * (one/rho_dryair(i,k))) + nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho_air(i,k), nwfa(i,k)*rho_air(i,k)) * (one/rho_air(i,k))) !> - Convert number concentrations from dry to moist gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k)) endif @@ -785,7 +785,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) / (one-spechum(i,k)) !> - Convert number concentration from moist to dry ni_mp(i,k) = gq0(i,k,ntinc) / (one-spechum(i,k)) - ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho_dryair(i,k), save_tcp(i,k)) * (one/rho_dryair(i,k))) + ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho_air(i,k), save_tcp(i,k)) * (one/rho_air(i,k))) !> - Convert number concentrations from dry to moist gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k)) endif diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index ec19945b0..82ddc95be 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -159,10 +159,6 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & ! Geopotential height in m2 s-2 to height in m hgt = phil/con_g - ! Density of air in kg m-3 and inverse density of air - rho = prsl/(con_rd*tgrs) - orho = 1.0/rho - ! Prior to calling the functions: make_DropletNumber, make_IceNumber, make_RainNumber, ! the incoming mixing ratios should be converted to units of mass/num per cubic meter ! rather than per kg of air. So, to pass back to the model state variables, @@ -178,6 +174,10 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & qs_mp = qs/(1.0_kind_phys-spechum) qg_mp = qg/(1.0_kind_phys-spechum) + ! Density of air in kg m-3 and inverse density of air + rho = 0.622*prsl/(con_rd*tgrs*(qv_mp+0.622)) + orho = 1.0/rho + !> - Convert number concentrations from moist to dry ni_mp = ni/(1.0_kind_phys-spechum) nr_mp = nr/(1.0_kind_phys-spechum) @@ -304,7 +304,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & ! If qc is in boundary conditions but nc is not, calculate nc from qc, rho and nwfa if (maxval(qc_mp)>0.0 .and. maxval(nc_mp)==0.0) then - nc_mp = make_DropletNumber(qc_mp*rho, nwfa) * orho + nc_mp = make_DropletNumber(qc_mp*rho, nwfa*rho) * orho end if ! If nc is in boundary conditions but qc is not, reset nc to zero @@ -428,6 +428,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & ! Air density real(kind_phys) :: rho(1:ncol,1:nlev) !< kg m-3 + !rho = 0.622*prsl/(con_rd*tgrs*(qv_mp+0.622)) ! Hydrometeors real(kind_phys) :: qv_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) real(kind_phys) :: qc_mp(1:ncol,1:nlev) !< kg kg-1 (dry mixing ratio) @@ -510,7 +511,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if !> - Density of air in kg m-3 - rho = prsl/(con_rd*tgrs) + rho = 0.622*prsl/(con_rd*tgrs*(qv_mp+0.622)) !> - Convert omega in Pa s-1 to vertical velocity w in m s-1 w = -omega/(rho*con_g) From 9cd399d66063956ae51c52b89e0ae14972fd9040 Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Tue, 23 Feb 2021 17:40:22 +0000 Subject: [PATCH 06/31] Replace dry-air density to moist-air density in RRTMGP. Fix gthe units of water-friendly aerosols in the call to make_DropletNumber. --- physics/GFS_rrtmgp_thompsonmp_pre.F90 | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/physics/GFS_rrtmgp_thompsonmp_pre.F90 b/physics/GFS_rrtmgp_thompsonmp_pre.F90 index ea27f3d2b..de1fa3547 100644 --- a/physics/GFS_rrtmgp_thompsonmp_pre.F90 +++ b/physics/GFS_rrtmgp_thompsonmp_pre.F90 @@ -147,11 +147,11 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do ! Cloud particle sizes and number concentrations... ! First, prepare cloud mixing-ratios and number concentrations for Calc_Re - rho = p_lay(1:nCol,1:nLev)/(con_rd*t_lay(1:nCol,1:nLev)) - orho = 1./rho do iLay = 1, nLev do iCol = 1, nCol qv_mp(iCol,iLay) = q_lay(iCol,iLay)/(1.-q_lay(iCol,iLay)) + rho(iCol,iLay) = 0.622*p_lay(1:nCol,1:nLev)/(con_rd*t_lay(1:nCol,1:nLev)*(qv_mp(iCol,iLay)+0.622)) + orho(iCol,iLay) = 1./rho(iCol,iLay) qc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq) / (1.-q_lay(iCol,iLay)) qi_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice) / (1.-q_lay(iCol,iLay)) qs_mp(iCol,iLay) = tracer(iCol,iLay,i_cldsnow) / (1.-q_lay(iCol,iLay)) @@ -169,7 +169,7 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do do iLay = 1, nLev do iCol = 1, nCol if (ltaerosol .and. qc_mp(iCol,iLay) > 1.e-12 .and. nc_mp(iCol,iLay) < 100.) then - nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho(iCol,iLay), nwfa(iCol,iLay)) * orho(iCol,iLay) + nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho(iCol,iLay), nwfa(iCol,iLay)*rho(iCol,iLay)) * orho(iCol,iLay) endif if (qi_mp(iCol,iLay) > 1.e-12 .and. ni_mp(iCol,iLay) < 100.) then ni_mp(iCol,iLay) = make_IceNumber(qi_mp(iCol,iLay)*rho(iCol,iLay), t_lay(iCol,iLay)) * orho(iCol,iLay) From cae708d26b4a4fe3fc06c0191fbe72fb8af0674c Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Wed, 24 Feb 2021 11:27:23 -0700 Subject: [PATCH 07/31] fix for one more consistency check of ice number --- physics/module_mp_thompson.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 444e867ae..2e6331c33 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -2888,8 +2888,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & niten(k) = -ni1d(k)*odts endif xni=MAX(0.,(ni1d(k) + niten(k)*dtsave)*rho(k)) - if (xni.gt.499.E3) & - niten(k) = (499.E3-ni1d(k)*rho(k))*odts*orho + if (xni.gt.9999.E3) & + niten(k) = (9999.E3-ni1d(k)*rho(k))*odts*orho !> - Rain tendency qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) & From f6d5fc1352c0a93fa9876cab200e5b3af05b4ac1 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Thu, 25 Feb 2021 09:09:29 -0700 Subject: [PATCH 08/31] small compile-time bug fixes --- physics/module_mp_thompson.F90 | 7 +++---- physics/mp_thompson.F90 | 22 +++++++++++----------- 2 files changed, 14 insertions(+), 15 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 2e6331c33..ed8309789 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1264,7 +1264,6 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & nwfa1d(k) = 11.1E6/rho(k) nifa1d(k) = naIN1*0.01/rho(k) enddo - nwfa1 = 11.1E6 endif !> - Call mp_thompson() @@ -3717,7 +3716,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & elseif (nc1d(k)*rho(k).lt.100.) then nu_c = 15 else - nu_c = NINT(1000.E6/(nc1d(k)*rho(k)) + 2 + 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 @@ -3941,7 +3940,7 @@ subroutine qr_acr_qg CLOSE(63) RETURN ! ----- RETURN 9234 CONTINUE - write(0,*) "Error writing //qr_acr_qg_file + write(0,*) "Error writing "//qr_acr_qg_file return ENDIF ENDIF @@ -4237,7 +4236,7 @@ subroutine freezeH2O(threads) write_thompson_tables = .false. good = 0 - INQUIRE(FILE=freeze_h2o_file",EXIST=lexist) + INQUIRE(FILE=freeze_h2o_file,EXIST=lexist) #ifdef MPI call MPI_BARRIER(mpi_communicator,ierr) #endif diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index ac8437262..79994674f 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -142,14 +142,14 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & qv = spechum/(1.0_kind_phys-spechum) ! Ensure we have 1st guess ice number where mass non-zero but no number. - where(qi <= 0.0) ni=0.0 - where(qi > 0 .and. ni <= 0.0) ni = make_IceNumber(qi*rho, tgrs) * orho - where(qi = 0.0 .and. ni > 0.0) ni=0.0 + where(qi .LE. 0.0) ni=0.0 + where(qi .GT. 0 .and. ni .LE. 0.0) ni = make_IceNumber(qi*rho, tgrs) * orho + where(qi .EQ. 0.0 .and. ni .GT. 0.0) ni=0.0 ! Ensure we have 1st guess rain number where mass non-zero but no number. - where(qr <= 0.0) nr=0.0 - where(qr > 0 .and. nr <= 0.0) nr = make_RainNumber(qr*rho, tgrs) * orho - where(qr = 0.0 .and. nr > 0.0) nr=0.0 + where(qr .LE. 0.0) nr=0.0 + where(qr .GT. 0 .and. nr .LE. 0.0) nr = make_RainNumber(qr*rho, tgrs) * orho + where(qr .EQ. 0.0 .and. nr .GT. 0.0) nr=0.0 !..Check for existing aerosol data, both CCN and IN aerosols. If missing @@ -227,13 +227,13 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & endif ! Ensure we have 1st guess cloud droplet number where mass non-zero but no number. - where(qc <= 0.0) nc=0.0 - where(qc > 0 .and. nc <= 0.0) nc = make_DropletNumber(qc*rho, nwfa*rho) * orho - where(qc = 0.0 .and. nc > 0.0) nc=0.0 + where(qc .LE. 0.0) nc=0.0 + where(qc .GT. 0 .and. nc .LE. 0.0) nc = make_DropletNumber(qc*rho, nwfa*rho) * orho + where(qc .EQ. 0.0 .and. nc .GT. 0.0) nc=0.0 ! Ensure non-negative aerosol number concentrations. - where(nwfa <= 0.0) nwfa = 1.1E6 - where(nifa <= 0.0) nifa = naIN1*0.01 + where(nwfa .LE. 0.0) nwfa = 1.1E6 + where(nifa .LE. 0.0) nifa = naIN1*0.01 else From a2d800b1136af477f3982766e2ee8ac94a2b6fd6 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 25 Feb 2021 12:09:14 -0700 Subject: [PATCH 09/31] Update submodule pointer for physics/rte-rrtmgp --- physics/rte-rrtmgp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/rte-rrtmgp b/physics/rte-rrtmgp index 566bee9cd..33c8a984c 160000 --- a/physics/rte-rrtmgp +++ b/physics/rte-rrtmgp @@ -1 +1 @@ -Subproject commit 566bee9cd6f9977e82d75d9b4964b20b1ff6163d +Subproject commit 33c8a984c17cf41be5d4c2928242e1b4239bfc40 From 5a9d2567d2cee270cf60e9fcdf35eed64105d4ca Mon Sep 17 00:00:00 2001 From: tanyasmirnova Date: Thu, 25 Feb 2021 19:17:34 +0000 Subject: [PATCH 10/31] Update submodule pointer for physics/rte-rrtmgp --- physics/rte-rrtmgp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/rte-rrtmgp b/physics/rte-rrtmgp index 566bee9cd..33c8a984c 160000 --- a/physics/rte-rrtmgp +++ b/physics/rte-rrtmgp @@ -1 +1 @@ -Subproject commit 566bee9cd6f9977e82d75d9b4964b20b1ff6163d +Subproject commit 33c8a984c17cf41be5d4c2928242e1b4239bfc40 From 4c74498e664b912a46df829bb828893a05c12ade Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 26 Feb 2021 07:24:07 -0700 Subject: [PATCH 11/31] Bugfixes, performance improvements and formatting updates in physics/GFS_suite_interstitial.F90 and physics/GFS_rrtmgp_thompsonmp_pre.F90 --- physics/GFS_rrtmgp_thompsonmp_pre.F90 | 140 +++++++++++++------------- physics/GFS_suite_interstitial.F90 | 11 +- 2 files changed, 74 insertions(+), 77 deletions(-) diff --git a/physics/GFS_rrtmgp_thompsonmp_pre.F90 b/physics/GFS_rrtmgp_thompsonmp_pre.F90 index de1fa3547..58e1bddea 100644 --- a/physics/GFS_rrtmgp_thompsonmp_pre.F90 +++ b/physics/GFS_rrtmgp_thompsonmp_pre.F90 @@ -16,24 +16,24 @@ module GFS_rrtmgp_thompsonmp_pre make_RainNumber use rrtmgp_lw_cloud_optics, only: radliq_lwr, radliq_upr, radice_lwr, radice_upr implicit none - + ! Parameters specific to THOMPSON MP scheme. real(kind_phys), parameter :: & rerain_def = 1000.0 ! Default rain radius to 1000 microns - + public GFS_rrtmgp_thompsonmp_pre_init, GFS_rrtmgp_thompsonmp_pre_run, GFS_rrtmgp_thompsonmp_pre_finalize - -contains + +contains ! ###################################################################################### ! ###################################################################################### subroutine GFS_rrtmgp_thompsonmp_pre_init() end subroutine GFS_rrtmgp_thompsonmp_pre_init - + ! ###################################################################################### ! ###################################################################################### !! \section arg_table_GFS_rrtmgp_thompsonmp_pre_run !! \htmlinclude GFS_rrtmgp_thompsonmp_pre_run.html -!! +!! subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, doLWrad, & i_cldliq, i_cldice, i_cldrain, i_cldsnow, i_cldgrpl, i_cldtot, i_cldliq_nc, & i_cldice_nc, i_twa, effr_in, p_lev, p_lay, tv_lay, t_lay, effrin_cldliq, & @@ -42,14 +42,14 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do imfdeepcnv_gf, doGP_cldoptics_PADE, doGP_cldoptics_LUT, & cld_frac, cld_lwp, cld_reliq, cld_iwp, cld_reice, cld_swp, cld_resnow, cld_rwp, & cld_rerain, precip_frac, errmsg, errflg) - - ! Inputs + + ! Inputs integer, intent(in) :: & nCol, & ! Number of horizontal grid points nLev, & ! Number of vertical layers ncnd, & ! Number of cloud condensation types. - nTracers, & ! Number of tracers from model. - i_cldliq, & ! Index into tracer array for cloud liquid amount. + nTracers, & ! Number of tracers from model. + i_cldliq, & ! Index into tracer array for cloud liquid amount. i_cldice, & ! cloud ice amount. i_cldrain, & ! cloud rain amount. i_cldsnow, & ! cloud snow amount. @@ -61,9 +61,9 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do imfdeepcnv, & ! Choice of mass-flux deep convection scheme imfdeepcnv_gf ! Flag for Grell-Freitas deep convection scheme logical, intent(in) :: & - doSWrad, & ! Call SW radiation? - doLWrad, & ! Call LW radiation - effr_in, & ! Use cloud effective radii provided by model? + doSWrad, & ! Call SW radiation? + doLWrad, & ! Call LW radiation + effr_in, & ! Use cloud effective radii provided by model? uni_cld, & ! Use provided cloud-fraction? lmfshal, & ! Flag for mass-flux shallow convection scheme used by Xu-Randall lmfdeep2, & ! Flag for some scale-aware mass-flux convection scheme active @@ -75,7 +75,7 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do con_g, & ! Physical constant: gravitational constant con_rd ! Physical constant: gas-constant for dry air - real(kind_phys), dimension(nCol,nLev), intent(in) :: & + real(kind_phys), dimension(nCol,nLev), intent(in) :: & tv_lay, & ! Virtual temperature (K) t_lay, & ! Temperature (K) qs_lay, & ! Saturation vapor pressure (Pa) @@ -83,13 +83,13 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do relhum, & ! Relative humidity p_lay, & ! Pressure at model-layers (Pa) cld_frac_mg ! Cloud-fraction from MG scheme. WTF????? - real(kind_phys), dimension(nCol,nLev+1), intent(in) :: & + real(kind_phys), dimension(nCol,nLev+1), intent(in) :: & p_lev ! Pressure at model-level interfaces (Pa) real(kind_phys), dimension(nCol, nLev, nTracers),intent(in) :: & - tracer ! Cloud condensate amount in layer by type () - + tracer ! Cloud condensate amount in layer by type () + ! In/Outs - real(kind_phys), dimension(nCol,nLev), intent(inout) :: & + real(kind_phys), dimension(nCol,nLev), intent(inout) :: & cld_frac, & ! Total cloud fraction cld_lwp, & ! Cloud liquid water path cld_reliq, & ! Cloud liquid effective radius @@ -97,31 +97,32 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do cld_reice, & ! Cloud ice effecive radius cld_swp, & ! Cloud snow water path cld_resnow, & ! Cloud snow effective radius - cld_rwp, & ! Cloud rain water path + cld_rwp, & ! Cloud rain water path cld_rerain, & ! Cloud rain effective radius precip_frac, & ! Precipitation fraction effrin_cldliq, & ! Effective radius for liquid cloud-particles (microns) effrin_cldice, & ! Effective radius for ice cloud-particles (microns) - effrin_cldsnow ! Effective radius for snow cloud-particles (microns) - - ! Outputs + effrin_cldsnow ! Effective radius for snow cloud-particles (microns) + + ! Outputs character(len=*), intent(out) :: & errmsg ! Error message - integer, intent(out) :: & + integer, intent(out) :: & errflg ! Error flag - + ! Local variables real(kind_phys) :: alpha0, pfac, tem1, cld_mr real(kind_phys), dimension(nCol, nLev, min(4,ncnd)) :: cld_condensate integer :: iCol,iLay,l - real(kind_phys), dimension(nCol,nLev) :: deltaP, deltaZ, rho, orho, re_cloud, re_ice,& + real(kind_phys) :: rho, orho + real(kind_phys), dimension(nCol,nLev) :: deltaP, deltaZ, re_cloud, re_ice,& re_snow, qv_mp, qc_mp, qi_mp, qs_mp, nc_mp, ni_mp, nwfa logical :: top_at_1 - + ! Initialize CCPP error handling variables errmsg = '' errflg = 0 - + if (.not. (doSWrad .or. doLWrad)) return ! Cloud condensate @@ -129,29 +130,30 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do cld_condensate(1:nCol,1:nLev,2) = tracer(1:nCol,1:nLev,i_cldice) ! -ice water cld_condensate(1:nCol,1:nLev,3) = tracer(1:nCol,1:nLev,i_cldrain) ! -rain water cld_condensate(1:nCol,1:nLev,4) = tracer(1:nCol,1:nLev,i_cldsnow) + &! -snow + grapuel - tracer(1:nCol,1:nLev,i_cldgrpl) - + tracer(1:nCol,1:nLev,i_cldgrpl) + ! Cloud water path (g/m2) - deltaP = abs(p_lev(:,2:nLev+1)-p_lev(:,1:nLev))/100. + deltaP = abs(p_lev(:,2:nLev+1)-p_lev(:,1:nLev))/100. do iLay = 1, nLev do iCol = 1, nCol - ! Compute liquid/ice condensate path from mixing ratios (kg/kg)->(g/m2) + ! Compute liquid/ice condensate path from mixing ratios (kg/kg)->(g/m2) tem1 = (1.0e5/con_g) * deltaP(iCol,iLay) cld_lwp(iCol,iLay) = max(0., cld_condensate(iCol,iLay,1) * tem1) cld_iwp(iCol,iLay) = max(0., cld_condensate(iCol,iLay,2) * tem1) cld_rwp(iCol,iLay) = max(0., cld_condensate(iCol,iLay,3) * tem1) - cld_swp(iCol,iLay) = max(0., cld_condensate(iCol,iLay,4) * tem1) + cld_swp(iCol,iLay) = max(0., cld_condensate(iCol,iLay,4) * tem1) enddo - enddo - + enddo + ! Cloud particle sizes and number concentrations... - - ! First, prepare cloud mixing-ratios and number concentrations for Calc_Re + + ! Prepare cloud mixing-ratios and number concentrations for calc_effectRad, + ! and update number concentrations, consistent with sub-grid clouds do iLay = 1, nLev - do iCol = 1, nCol + do iCol = 1, nCol qv_mp(iCol,iLay) = q_lay(iCol,iLay)/(1.-q_lay(iCol,iLay)) - rho(iCol,iLay) = 0.622*p_lay(1:nCol,1:nLev)/(con_rd*t_lay(1:nCol,1:nLev)*(qv_mp(iCol,iLay)+0.622)) - orho(iCol,iLay) = 1./rho(iCol,iLay) + rho = 0.622*p_lay(iCol,iLay)/(con_rd*t_lay(iCol,iLay)*(qv_mp(iCol,iLay)+0.622)) + orho = 1./rho qc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq) / (1.-q_lay(iCol,iLay)) qi_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice) / (1.-q_lay(iCol,iLay)) qs_mp(iCol,iLay) = tracer(iCol,iLay,i_cldsnow) / (1.-q_lay(iCol,iLay)) @@ -159,24 +161,18 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do ni_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice_nc) / (1.-q_lay(iCol,iLay)) if (ltaerosol) then nwfa(iCol,iLay) = tracer(iCol,iLay,i_twa) + if (qc_mp(iCol,iLay) > 1.e-12 .and. nc_mp(iCol,iLay) < 100.) then + nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho, nwfa(iCol,iLay)*rho) * orho + endif else - nc_mp(iCol,iLay) = nt_c*orho(iCol,iLay) - endif - enddo - enddo - - ! Update number concentration, consistent with sub-grid clouds - do iLay = 1, nLev - do iCol = 1, nCol - if (ltaerosol .and. qc_mp(iCol,iLay) > 1.e-12 .and. nc_mp(iCol,iLay) < 100.) then - nc_mp(iCol,iLay) = make_DropletNumber(qc_mp(iCol,iLay)*rho(iCol,iLay), nwfa(iCol,iLay)*rho(iCol,iLay)) * orho(iCol,iLay) + nc_mp(iCol,iLay) = nt_c*orho endif if (qi_mp(iCol,iLay) > 1.e-12 .and. ni_mp(iCol,iLay) < 100.) then - ni_mp(iCol,iLay) = make_IceNumber(qi_mp(iCol,iLay)*rho(iCol,iLay), t_lay(iCol,iLay)) * orho(iCol,iLay) + ni_mp(iCol,iLay) = make_IceNumber(qi_mp(iCol,iLay)*rho, t_lay(iCol,iLay)) * orho endif enddo enddo - + ! Compute effective radii for liquid/ice/snow using subgrid scale clouds ! Call Thompson's subroutine to compute effective radii do iCol=1,nCol @@ -184,14 +180,14 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do nc_mp(iCol,:), qi_mp(iCol,:), ni_mp(iCol,:), qs_mp(iCol,:), & re_cloud(iCol,:), re_ice(iCol,:), re_snow(iCol,:), 1, nLev ) enddo - - ! Scale Thompson's effective radii from meter to micron + + ! Scale Thompson's effective radii from meter to micron effrin_cldliq(1:nCol,1:nLev) = re_cloud(1:nCol,1:nLev)*1.e6 effrin_cldice(1:nCol,1:nLev) = re_ice(1:nCol,1:nLev)*1.e6 effrin_cldsnow(1:nCol,1:nLev) = re_snow(1:nCol,1:nLev)*1.e6 - - ! Bound effective radii for RRTMGP, LUT's for cloud-optics go from - ! 2.5 - 21.5 microns for liquid clouds, + + ! Bound effective radii for RRTMGP, LUT's for cloud-optics go from + ! 2.5 - 21.5 microns for liquid clouds, ! 10 - 180 microns for ice-clouds if (doGP_cldoptics_PADE .or. doGP_cldoptics_LUT) then where(effrin_cldliq .lt. radliq_lwr) effrin_cldliq = radliq_lwr @@ -205,32 +201,32 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do cld_reice(1:nCol,1:nLev) = effrin_cldice(1:nCol,1:nLev) cld_resnow(1:nCol,1:nLev) = effrin_cldsnow(1:nCol,1:nLev) cld_rerain(1:nCol,1:nLev) = rerain_def - + ! Compute cloud-fraction. Else, use value provided - if(.not. do_mynnedmf .or. imfdeepcnv .ne. imfdeepcnv_gf ) then ! MYNN PBL or GF conv + if(.not. do_mynnedmf .or. imfdeepcnv .ne. imfdeepcnv_gf ) then ! MYNN PBL or GF conv ! Cloud-fraction if (uni_cld) then - cld_frac(1:nCol,1:nLev) = cld_frac_mg(1:nCol,1:nLev) + cld_frac(1:nCol,1:nLev) = cld_frac_mg(1:nCol,1:nLev) else if( lmfshal) alpha0 = 100. ! Default (from GATE simulations) if(.not. lmfshal) alpha0 = 2000. - ! Xu-Randall (1996) cloud-fraction + ! Xu-Randall (1996) cloud-fraction do iLay = 1, nLev do iCol = 1, nCol cld_mr = cld_condensate(iCol,iLay,1) + cld_condensate(iCol,iLay,2) + & cld_condensate(iCol,iLay,4) - cld_frac(iCol,iLay) = cld_frac_XuRandall(p_lay(iCol,iLay), & - qs_lay(iCol,iLay), relhum(iCol,iLay), cld_mr, alpha0) + cld_frac(iCol,iLay) = cld_frac_XuRandall(p_lay(iCol,iLay), & + qs_lay(iCol,iLay), relhum(iCol,iLay), cld_mr, alpha0) enddo enddo endif endif - + ! Precipitation fraction (Hack. For now use cloud-fraction) precip_frac(1:nCol,1:nLev) = cld_frac(1:nCol,1:nLev) - + end subroutine GFS_rrtmgp_thompsonmp_pre_run - + ! ###################################################################################### ! ###################################################################################### subroutine GFS_rrtmgp_thompsonmp_pre_finalize() @@ -241,11 +237,11 @@ end subroutine GFS_rrtmgp_thompsonmp_pre_finalize ! Xu-Randall(1996) A Semiempirical Cloudiness Parameterization for Use in Climate Models ! https://doi.org/10.1175/1520-0469(1996)053<3084:ASCPFU>2.0.CO;2 ! - ! cld_frac = {1-exp[-alpha*cld_mr/((1-relhum)*qs_lay)**lambda]}*relhum**P + ! cld_frac = {1-exp[-alpha*cld_mr/((1-relhum)*qs_lay)**lambda]}*relhum**P ! ! ###################################################################################### function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha) - + ! Inputs real(kind_phys), intent(in) :: & p_lay, & ! Pressure (Pa) @@ -253,7 +249,7 @@ function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha) relhum, & ! Relative humidity cld_mr, & ! Total cloud mixing ratio alpha ! Scheme parameter (default=100) - + ! Outputs real(kind_phys) :: cld_frac_XuRandall @@ -262,21 +258,21 @@ function cld_frac_XuRandall(p_lay, qs_lay, relhum, cld_mr, alpha) ! Parameters real(kind_phys) :: & - lambda = 0.50, & ! + lambda = 0.50, & ! P = 0.25 - + clwt = 1.0e-6 * (p_lay*0.001) if (cld_mr > clwt) then onemrh = max(1.e-10, 1.0 - relhum) tem1 = alpha / min(max((onemrh*qs_lay)**lambda,0.0001),1.0) tem2 = max(min(tem1*(cld_mr - clwt), 50.0 ), 0.0 ) tem3 = sqrt(sqrt(relhum)) ! This assumes "p" = 0.25. Identical, but cheaper than relhum**p - ! + ! cld_frac_XuRandall = max( tem3*(1.0-exp(-tem2)), 0.0 ) else cld_frac_XuRandall = 0.0 endif - + return end function end module GFS_rrtmgp_thompsonmp_pre diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index 464ede8cc..b51b6f9c8 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -692,7 +692,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to ! local variables integer :: i,k,n,tracers - real(kind=kind_phys), dimension(im,levs) :: rho_air + real(kind=kind_phys) :: rho, orho real(kind=kind_phys), dimension(im,levs) :: qv_mp !< kg kg-1 (dry mixing ratio) real(kind=kind_phys), dimension(im,levs) :: qc_mp !< kg kg-1 (dry mixing ratio) real(kind=kind_phys), dimension(im,levs) :: qi_mp !< kg kg-1 (dry mixing ratio) @@ -744,14 +744,15 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to do i=1,im !> - Convert specific humidity to dry mixing ratio qv_mp(i,k) = spechum(i,k) / (one-spechum(i,k)) - !> - Density of air in kg m-3 - rho_air(i,k) = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+0.622)) + !> - Density of air in kg m-3 and inverse density + rho = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+0.622)) + orho = one/rho if (ntlnc>0) then !> - Convert moist mixing ratio to dry mixing ratio qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) / (one-spechum(i,k)) !> - Convert number concentration from moist to dry nc_mp(i,k) = gq0(i,k,ntlnc) / (one-spechum(i,k)) - nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho_air(i,k), nwfa(i,k)*rho_air(i,k)) * (one/rho_air(i,k))) + nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) !> - Convert number concentrations from dry to moist gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k)) endif @@ -760,7 +761,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) / (one-spechum(i,k)) !> - Convert number concentration from moist to dry ni_mp(i,k) = gq0(i,k,ntinc) / (one-spechum(i,k)) - ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho_air(i,k), save_tcp(i,k)) * (one/rho_air(i,k))) + ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) !> - Convert number concentrations from dry to moist gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k)) endif From a13ef16702a01b971192c319afe812e46a40b3f9 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Fri, 26 Feb 2021 15:07:12 -0700 Subject: [PATCH 12/31] adopt moist air density in place of dry air density --- physics/mp_thompson.F90 | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 79994674f..84c30a18e 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -126,10 +126,6 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & ! Geopotential height in m2 s-2 to height in m hgt = phil/con_g - ! Density of air in kg m-3 and inverse density of air - rho = prsl/(con_rd*tgrs) - orho = 1.0/rho - ! Ensure non-negative mass mixing ratios of all water variables where(spechum<0) spechum = 1.0E-10 ! COMMENT, gthompsn, spechum should *never* be identically zero. where(qc<0) qc = 0.0 @@ -141,6 +137,10 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & ! Convert specific humidity to water vapor mixing ratio qv = spechum/(1.0_kind_phys-spechum) + ! Density of moist air in kg m-3 and inverse density of air + rho = 0.622*prsl/(con_rd*tgrs*(qv+0.622)) + orho = 1.0/rho + ! Ensure we have 1st guess ice number where mass non-zero but no number. where(qi .LE. 0.0) ni=0.0 where(qi .GT. 0 .and. ni .LE. 0.0) ni = make_IceNumber(qi*rho, tgrs) * orho From c6f1cab8cc644c803739047aa238c7155102c0ef Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 1 Mar 2021 07:32:32 -0700 Subject: [PATCH 13/31] Fix bug due to merge in physics/mp_thompson.F90 --- physics/mp_thompson.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 337161140..fbfbe72a9 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -431,7 +431,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if !> - Density of air in kg m-3 - rho = 0.622*prsl/(con_rd*tgrs*(qv_mp+0.622)) + rho = 0.622*prsl/(con_rd*tgrs*(qv+0.622)) !> - Convert omega in Pa s-1 to vertical velocity w in m s-1 w = -omega/(rho*con_g) From 2a2d75076ab43f60b7ea7776f5e3b0e7dc98d9dd Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Mon, 1 Mar 2021 08:40:16 -0700 Subject: [PATCH 14/31] finalze the correct min,max aerosol numbers with regard to air density --- physics/module_mp_thompson.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index ed8309789..09bb5c939 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -1261,8 +1261,8 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & else do k = kts, kte nc1d(k) = Nt_c/rho(k) - nwfa1d(k) = 11.1E6/rho(k) - nifa1d(k) = naIN1*0.01/rho(k) + nwfa1d(k) = 11.1E6 + nifa1d(k) = naIN1*0.01 enddo endif @@ -1785,8 +1785,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & 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)) - nwfa(k) = MAX(11.1E6, MIN(9999.E6, nwfa1d(k)*rho(k))) - nifa(k) = MAX(naIN1*0.01, MIN(9999.E6, nifa1d(k)*rho(k))) + nwfa(k) = MAX(11.1E6*rho(k), MIN(9999.E6*rho(k), nwfa1d(k)*rho(k))) + nifa(k) = MAX(naIN1*0.01*rho(k), MIN(9999.E6*rho(k), nifa1d(k)*rho(k))) mvd_r(k) = D0r mvd_c(k) = D0c @@ -2987,7 +2987,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & ocp(k) = 1./(Cp*(1.+0.887*qv(k))) lvt2(k)=lvap(k)*lvap(k)*ocp(k)*oRv*otemp*otemp - nwfa(k) = MAX(11.1E6, (nwfa1d(k) + nwfaten(k)*DT)*rho(k)) + nwfa(k) = MAX(11.1E6*rho(k), (nwfa1d(k) + nwfaten(k)*DT)*rho(k)) enddo do k = kts, kte From 3f2c4e5f3278e88934ebc5119c8e20707feb9e02 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Thu, 4 Mar 2021 11:56:26 -0700 Subject: [PATCH 15/31] Make convert_dry_rho an input argument, fix a bug in the non-aerosol-aware version (nc not allocated), reduce noise written to stdout --- physics/module_mp_thompson.F90 | 11 +++--- physics/mp_thompson.F90 | 62 ++++++++++++++++++++++++++++------ physics/mp_thompson.meta | 16 +++++++++ 3 files changed, 72 insertions(+), 17 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 09bb5c939..25466f412 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -423,9 +423,7 @@ MODULE module_mp_thompson !..If SIONlib isn't used, write Thompson tables with master MPI task !.. after computing them in thompson_init -#ifndef SION LOGICAL:: thompson_table_writer -#endif !+---+ !+---+-----------------------------------------------------------------+ @@ -767,7 +765,7 @@ SUBROUTINE thompson_init(mpicomm, mpirank, mpiroot, & precomputed_tables = .false. if (mpirank==mpiroot) write(0,*) "An error occurred reading Thompson tables from disk, recalculate" end if -#else +#endif ! Standard tables are only written by master MPI task; ! (physics init cannot be called by multiple threads, ! hence no need to test for a specific thread number) @@ -776,7 +774,6 @@ SUBROUTINE thompson_init(mpicomm, mpirank, mpiroot, & else thompson_table_writer = .false. end if -#endif precomputed_tables_1: if (.not.precomputed_tables) then @@ -3847,7 +3844,7 @@ subroutine qr_acr_qg #ifndef SION if (thompson_table_writer) write_thompson_tables = .true. #endif - write(0,*) "ThompMP: computing qr_acr_qg" + if (thompson_table_writer) write(0,*) "ThompMP: computing qr_acr_qg" do n2 = 1, nbr ! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2)) vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) & @@ -4029,7 +4026,7 @@ subroutine qr_acr_qs #ifndef SION if (thompson_table_writer) write_thompson_tables = .true. #endif - write(0,*) "ThompMP: computing qr_acr_qs" + if (thompson_table_writer) write(0,*) "ThompMP: computing qr_acr_qs" do n2 = 1, nbr ! vr(n2) = av_r*Dr(n2)**bv_r * DEXP(-fv_r*Dr(n2)) vr(n2) = -0.1021 + 4.932E3*Dr(n2) - 0.9551E6*Dr(n2)*Dr(n2) & @@ -4284,7 +4281,7 @@ subroutine freezeH2O(threads) #ifndef SION if (thompson_table_writer) write_thompson_tables = .true. #endif - write(0,*) "ThompMP: computing freezeH2O" + if (thompson_table_writer) write(0,*) "ThompMP: computing freezeH2O" orho_w = 1./rho_w diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index fbfbe72a9..252944b77 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -21,7 +21,6 @@ module mp_thompson private logical :: is_initialized = .False. - logical :: convert_dry_rho = .False. contains @@ -31,6 +30,7 @@ module mp_thompson !! subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & imp_physics, imp_physics_thompson, & + convert_dry_rho, & spechum, qc, qr, qi, qs, qg, ni, nr, & is_aerosol_aware, nc, nwfa2d, nifa2d, & nwfa, nifa, tgrs, prsl, phil, area, & @@ -48,6 +48,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & integer, intent(in ) :: imp_physics integer, intent(in ) :: imp_physics_thompson ! Hydrometeors + logical, intent(in ) :: convert_dry_rho real(kind_phys), intent(inout) :: spechum(:,:) real(kind_phys), intent(inout) :: qc(:,:) real(kind_phys), intent(inout) :: qr(:,:) @@ -83,10 +84,11 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & integer, intent( out) :: errflg ! - real(kind_phys) :: qv(1:ncol,1:nlev) ! kg kg-1 (water vapor mixing ratio) - real(kind_phys) :: hgt(1:ncol,1:nlev) ! m - real(kind_phys) :: rho(1:ncol,1:nlev) ! kg m-3 - real(kind_phys) :: orho(1:ncol,1:nlev) ! m3 kg-1 + real(kind_phys) :: qv(1:ncol,1:nlev) ! kg kg-1 (water vapor mixing ratio) + real(kind_phys) :: hgt(1:ncol,1:nlev) ! m + real(kind_phys) :: rho(1:ncol,1:nlev) ! kg m-3 + real(kind_phys) :: orho(1:ncol,1:nlev) ! m3 kg-1 + real(kind_phys) :: nc_local(1:ncol,1:nlev) ! needed because nc is only allocated if is_aerosol_aware is true ! real (kind=kind_phys) :: h_01, airmass, niIN3, niCCN3 integer :: i, k @@ -134,9 +136,28 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & where(qs<0) qs = 0.0 where(qg<0) qg = 0.0 - ! Convert specific humidity to water vapor mixing ratio + !> - Convert specific humidity to water vapor mixing ratio. + !> - Also, hydrometeor variables are mass or number mixing ratio + !> - either kg of species per kg of dry air, or per kg of (dry + vapor). + qv = spechum/(1.0_kind_phys-spechum) + if (convert_dry_rho) then + qc = qc/(1.0_kind_phys-spechum) + qr = qr/(1.0_kind_phys-spechum) + qi = qi/(1.0_kind_phys-spechum) + qs = qs/(1.0_kind_phys-spechum) + qg = qg/(1.0_kind_phys-spechum) + + ni = ni/(1.0_kind_phys-spechum) + nr = nr/(1.0_kind_phys-spechum) + if (is_aerosol_aware) then + nc = nc/(1.0_kind_phys-spechum) + nwfa = nwfa/(1.0_kind_phys-spechum) + nifa = nifa/(1.0_kind_phys-spechum) + end if + end if + ! Density of moist air in kg m-3 and inverse density of air rho = 0.622*prsl/(con_rd*tgrs*(qv+0.622)) orho = 1.0/rho @@ -229,17 +250,20 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & ! Ensure we have 1st guess cloud droplet number where mass non-zero but no number. where(qc .LE. 0.0) nc=0.0 where(qc .GT. 0 .and. nc .LE. 0.0) nc = make_DropletNumber(qc*rho, nwfa*rho) * orho - where(qc .EQ. 0.0 .and. nc .GT. 0.0) nc=0.0 + where(qc .EQ. 0.0 .and. nc .GT. 0.0) nc = 0.0 ! Ensure non-negative aerosol number concentrations. where(nwfa .LE. 0.0) nwfa = 1.1E6 where(nifa .LE. 0.0) nifa = naIN1*0.01 + ! Copy to local array for calculating cloud effective radii below + nc_local = nc + else ! Constant droplet concentration for single moment cloud water as in ! module_mp_thompson.F90, only needed for effective radii calculation - nc = Nt_c/rho + nc_local = Nt_c/rho end if @@ -247,8 +271,8 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & if (present(re_cloud) .and. present(re_ice) .and. present(re_snow)) then ! Effective radii [m] are now intent(out), bounds applied in calc_effectRad do i = 1, ncol - call calc_effectRad (tgrs(i,:), prsl(i,:), qv(i,:), qc(i,:), & - nc(i,:), qi(i,:), ni(i,:), qs(i,:), & + call calc_effectRad (tgrs(i,:), prsl(i,:), qv(i,:), qc(i,:), & + nc_local(i,:), qi(i,:), ni(i,:), qs(i,:), & re_cloud(i,:), re_ice(i,:), re_snow(i,:), 1, nlev) do k = 1, nlev re_cloud(i,k) = MAX(re_qc_min, MIN(re_cloud(i,k), re_qc_max)) @@ -271,6 +295,22 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & return end if + if (convert_dry_rho) then + !qc = qc/(1.0_kind_phys+qv) + !qr = qr/(1.0_kind_phys+qv) + !qi = qi/(1.0_kind_phys+qv) + !qs = qs/(1.0_kind_phys+qv) + !qg = qg/(1.0_kind_phys+qv) + + ni = ni/(1.0_kind_phys+qv) + nr = nr/(1.0_kind_phys+qv) + if (is_aerosol_aware) then + nc = nc/(1.0_kind_phys+qv) + nwfa = nwfa/(1.0_kind_phys+qv) + nifa = nifa/(1.0_kind_phys+qv) + end if + end if + is_initialized = .true. end subroutine mp_thompson_init @@ -283,6 +323,7 @@ end subroutine mp_thompson_init !>\section gen_thompson_hrrr Thompson MP General Algorithm !>@{ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & + convert_dry_rho, & spechum, qc, qr, qi, qs, qg, ni, nr, & is_aerosol_aware, nc, nwfa, nifa, & nwfa2d, nifa2d, & @@ -303,6 +344,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & real(kind_phys), intent(in ) :: con_g real(kind_phys), intent(in ) :: con_rd ! Hydrometeors + logical, intent(in ) :: convert_dry_rho real(kind_phys), intent(inout) :: spechum(1:ncol,1:nlev) real(kind_phys), intent(inout) :: qc(1:ncol,1:nlev) real(kind_phys), intent(inout) :: qr(1:ncol,1:nlev) diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index 4cfee6afc..ed54f8d02 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -65,6 +65,14 @@ type = integer intent = in optional = F +[convert_dry_rho] + standard_name = flag_for_converting_hydrometeors_from_moist_to_dry_air + long_name = flag for converting hydrometeors from moist to dry air + units = flag + dimensions = () + type = logical + intent = in + optional = F [spechum] standard_name = water_vapor_specific_humidity long_name = water vapor specific humidity @@ -341,6 +349,14 @@ kind = kind_phys intent = in optional = F +[convert_dry_rho] + standard_name = flag_for_converting_hydrometeors_from_moist_to_dry_air + long_name = flag for converting hydrometeors from moist to dry air + units = flag + dimensions = () + type = logical + intent = in + optional = F [spechum] standard_name = water_vapor_specific_humidity_updated_by_physics long_name = water vapor specific humidity From 15f1aa751a67aa2593f1700f951265f49d735947 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 5 Mar 2021 10:15:27 -0700 Subject: [PATCH 16/31] More bugfixes for Thompson MP: set is_aerosol_aware in mp_thompson_init -> init_thompson, apply bounds after calls to calc_effectRad in radiation, set intent(out) variables in calc_effectRad --- physics/GFS_rrtmg_pre.F90 | 16 ++++++++++++---- physics/GFS_rrtmgp_thompsonmp_pre.F90 | 11 +++++++++-- physics/module_mp_thompson.F90 | 20 +++++++++++++++++++- physics/mp_thompson.F90 | 9 +++++---- 4 files changed, 45 insertions(+), 11 deletions(-) diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index ab488c687..fee9b8815 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -68,10 +68,13 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & use surface_perturbation, only: cdfnor ! For Thompson MP - use module_mp_thompson, only: calc_effectRad, Nt_c - use module_mp_thompson_make_number_concentrations, only: & - make_IceNumber, & - make_DropletNumber, & + use module_mp_thompson, only: calc_effectRad, Nt_c, & + re_qc_min, re_qc_max, & + re_qi_min, re_qi_max, & + re_qs_min, re_qs_max + use module_mp_thompson_make_number_concentrations, only: & + make_IceNumber, & + make_DropletNumber, & make_RainNumber implicit none @@ -779,6 +782,11 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & call calc_effectRad (tlyr(i,:), plyr(i,:)*100., 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 ) + do k=1,lm + re_cloud(i,k) = MAX(re_qc_min, MIN(re_cloud(i,k), re_qc_max)) + re_ice(i,k) = MAX(re_qi_min, MIN(re_ice(i,k), re_qi_max)) + re_snow(i,k) = MAX(re_qs_min, MIN(re_snow(i,k), re_qs_max)) + end do end do ! Scale Thompson's effective radii from meter to micron do k=1,lm diff --git a/physics/GFS_rrtmgp_thompsonmp_pre.F90 b/physics/GFS_rrtmgp_thompsonmp_pre.F90 index 58e1bddea..08cc1b292 100644 --- a/physics/GFS_rrtmgp_thompsonmp_pre.F90 +++ b/physics/GFS_rrtmgp_thompsonmp_pre.F90 @@ -8,8 +8,10 @@ module GFS_rrtmgp_thompsonmp_pre use rrtmgp_aux, only: & check_error_msg use module_mp_thompson, only: & - calc_effectRad, & - Nt_c + calc_effectRad, Nt_c, & + re_qc_min, re_qc_max, & + re_qi_min, re_qi_max, & + re_qs_min, re_qs_max use module_mp_thompson_make_number_concentrations, only: & make_IceNumber, & make_DropletNumber, & @@ -179,6 +181,11 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do call calc_effectRad (t_lay(iCol,:), p_lay(iCol,:), qv_mp(iCol,:), qc_mp(iCol,:), & nc_mp(iCol,:), qi_mp(iCol,:), ni_mp(iCol,:), qs_mp(iCol,:), & re_cloud(iCol,:), re_ice(iCol,:), re_snow(iCol,:), 1, nLev ) + do iLay = 1, nLev + re_cloud(iCol,iLay) = MAX(re_qc_min, MIN(re_cloud(iCol,iLay), re_qc_max)) + re_ice(iCol,iLay) = MAX(re_qi_min, MIN(re_ice(iCol,iLay), re_qi_max)) + re_snow(iCol,iLay) = MAX(re_qs_min, MIN(re_snow(iCol,iLay), re_qs_max)) + end do enddo ! Scale Thompson's effective radii from meter to micron diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 25466f412..8ba269388 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -438,11 +438,13 @@ MODULE module_mp_thompson !! lookup tables in Thomspson scheme. !>\section gen_thompson_init thompson_init General Algorithm !> @{ - SUBROUTINE thompson_init(mpicomm, mpirank, mpiroot, & + SUBROUTINE thompson_init(is_aerosol_aware_in, & + mpicomm, mpirank, mpiroot, & threads, errmsg, errflg) IMPLICIT NONE + LOGICAL, INTENT(IN) :: is_aerosol_aware_in INTEGER, INTENT(IN) :: mpicomm, mpirank, mpiroot INTEGER, INTENT(IN) :: threads CHARACTER(len=*), INTENT(INOUT) :: errmsg @@ -458,6 +460,16 @@ SUBROUTINE thompson_init(mpicomm, mpirank, mpiroot, & LOGICAL, PARAMETER :: precomputed_tables = .FALSE. #endif +! Set module variable is_aerosol_aware + is_aerosol_aware = is_aerosol_aware_in + if (mpirank==mpiroot) then + if (is_aerosol_aware) then + write (0,'(a)') 'Using aerosol-aware version of Thompson microphysics' + else + write (0,'(a)') 'Using non-aerosol-aware version of Thompson microphysics' + end if + end if + micro_init = .FALSE. !> - Allocate space for lookup tables (J. Michalakes 2009Jun08). @@ -5218,6 +5230,8 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & lamc = (nc(k)*am_r*g_ratio(inu_c)/rc(k))**obmr re_qc1d(k) = SNGL(0.5D0 * DBLE(3.+inu_c)/lamc) enddo + else + re_qc1d(:) = 0.0D0 endif if (has_qi) then @@ -5226,6 +5240,8 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi re_qi1d(k) = SNGL(0.5D0 * DBLE(3.+mu_i)/lami) enddo + else + re_qi1d(:) = 0.0D0 endif if (has_qs) then @@ -5266,6 +5282,8 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & smoc = a_ * smo2**b_ re_qs1d(k) = 0.5*(smoc/smob) enddo + else + re_qs1d(:) = 0.0D0 endif end subroutine calc_effectRad diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 252944b77..5d5e631f5 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -115,8 +115,9 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & end if ! Call Thompson init - call thompson_init(mpicomm=mpicomm, mpirank=mpirank, mpiroot=mpiroot, & - threads=threads, errmsg=errmsg, errflg=errflg) + call thompson_init(is_aerosol_aware_in=is_aerosol_aware, mpicomm=mpicomm, & + mpirank=mpirank, mpiroot=mpiroot, threads=threads, & + errmsg=errmsg, errflg=errflg) if (errflg /= 0) return ! For restart runs, the init is done here @@ -276,8 +277,8 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & re_cloud(i,:), re_ice(i,:), re_snow(i,:), 1, nlev) do k = 1, nlev re_cloud(i,k) = MAX(re_qc_min, MIN(re_cloud(i,k), re_qc_max)) - re_ice(i,k) = MAX(re_qi_min, MIN(re_ice(i,k), re_qi_max)) - re_snow(i,k) = MAX(re_qs_min, MIN(re_snow(i,k), re_qs_max)) + re_ice(i,k) = MAX(re_qi_min, MIN(re_ice(i,k), re_qi_max)) + re_snow(i,k) = MAX(re_qs_min, MIN(re_snow(i,k), re_qs_max)) end do end do !! Convert to micron: required for bit-for-bit identical restarts; From 1f664cfaa1f31bc57e62d8b12af6ff7bd44f531e Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Fri, 5 Mar 2021 15:26:47 -0700 Subject: [PATCH 17/31] update GFS_phys_time_vary.scm for NoahMP init changes and UGWP v1 --- physics/GFS_phys_time_vary.scm.F90 | 481 +++++++++++++++++-- physics/GFS_phys_time_vary.scm.meta | 693 +++++++++++++++++++++++++++- 2 files changed, 1136 insertions(+), 38 deletions(-) diff --git a/physics/GFS_phys_time_vary.scm.F90 b/physics/GFS_phys_time_vary.scm.F90 index a1acc3fa0..dabd53e94 100644 --- a/physics/GFS_phys_time_vary.scm.F90 +++ b/physics/GFS_phys_time_vary.scm.F90 @@ -2,8 +2,8 @@ !! Contains code related to GFS physics suite setup (physics part of time_vary_step) !>\defgroup mod_GFS_phys_time_vary GFS Physics Time Update -!! This module contains GFS physics time vary subroutines including ozone, stratospheric water vapor, -!! aerosol, IN&CCN and surface properties updates. +!! This module contains GFS physics time vary subroutines including ozone, stratospheric water vapor, +!! aerosol, IN&CCN and surface properties updates. !> @{ module GFS_phys_time_vary @@ -23,10 +23,19 @@ module GFS_phys_time_vary use iccn_def, only : ciplin, ccnin, ci_pres use iccninterp, only : read_cidata, setindxci, ciinterpol -#if 0 + use cires_tauamf_data, only: cires_indx_ugwp, read_tau_amf, tau_amf_interp + use cires_tauamf_data, only: tau_limb, days_limb, ugwp_taulat + !--- variables needed for calculating 'sncovr' use namelist_soilveg, only: salp_data, snupx -#endif + use set_soilveg_mod, only: set_soilveg + + ! --- needed for Noah MP init + use noahmp_tables, only: laim_table,saim_table,sla_table, & + bexp_table,smcmax_table,smcwlt_table, & + dwsat_table,dksat_table,psisat_table, & + isurban_table,isbarren_table, & + isice_table,iswater_table implicit none @@ -36,9 +45,13 @@ module GFS_phys_time_vary logical :: is_initialized = .false. - real(kind=kind_phys), parameter :: con_hr = 3600.0_kind_phys - real(kind=kind_phys), parameter :: con_99 = 99.0_kind_phys - real(kind=kind_phys), parameter :: con_100 = 100.0_kind_phys + real(kind=kind_phys), parameter :: con_hr = 3600.0_kind_phys + real(kind=kind_phys), parameter :: con_99 = 99.0_kind_phys + real(kind=kind_phys), parameter :: con_100 = 100.0_kind_phys + real(kind=kind_phys), parameter :: missing_value = 9.99e20_kind_phys + real(kind=kind_phys), parameter :: drythresh = 1.e-4_kind_phys + real(kind=kind_phys), parameter :: zero = 0.0_kind_phys + real(kind=kind_phys), parameter :: one = 1.0_kind_phys contains @@ -52,7 +65,14 @@ subroutine GFS_phys_time_vary_init ( jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl, & jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm, & jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, imap, jmap, & - errmsg, errflg) + do_ugwp_v1, jindx1_tau, jindx2_tau, ddy_j1tau, ddy_j2tau, & + isot, ivegsrc, nlunit, sncovr, sncovr_ice, lsm, lsm_noahmp, lsm_ruc, min_seaice, & + fice, landfrac, vtype, weasd, lsoil, zs, dzs, lsnow_lsm_lbound, lsnow_lsm_ubound, & + tvxy, tgxy, tahxy, canicexy, canliqxy, eahxy, cmxy, chxy, fwetxy, sneqvoxy, alboldxy,& + qsnowxy, wslakexy, albdvis, albdnir, albivis, albinir, emiss, taussxy, waxy, wtxy, & + zwtxy, xlaixy, xsaixy, lfmassxy, stmassxy, rtmassxy, woodxy, stblcpxy, fastcpxy, & + smcwtdxy, deeprechxy, rechxy, snowxy, snicexy, snliqxy, tsnoxy , smoiseq, zsnsoxy, & + slc, smc, stc, tsfcl, snowd, canopy, tg3, stype, con_t0c, nthrds, errmsg, errflg) implicit none @@ -71,12 +91,85 @@ subroutine GFS_phys_time_vary_init ( integer, intent(inout) :: jindx1_ci(:), jindx2_ci(:), iindx1_ci(:), iindx2_ci(:) real(kind_phys), intent(inout) :: ddy_ci(:), ddx_ci(:) integer, intent(inout) :: imap(:), jmap(:) - + logical, intent(in) :: do_ugwp_v1 + real(kind_phys), intent(inout) :: ddy_j1tau(:), ddy_j2tau(:) + integer, intent(inout) :: jindx1_tau(:), jindx2_tau(:) + + integer, intent(in) :: isot, ivegsrc, nlunit + real(kind_phys), intent(inout) :: sncovr(:), sncovr_ice(:) + integer, intent(in) :: lsm, lsm_noahmp, lsm_ruc + real(kind_phys), intent(in) :: min_seaice, fice(:) + real(kind_phys), intent(in) :: landfrac(:), vtype(:) + real(kind_phys), intent(inout) :: weasd(:) + + ! NoahMP - only allocated when NoahMP is used + integer, intent(in) :: lsoil, lsnow_lsm_lbound, lsnow_lsm_ubound + real(kind_phys), intent(in) :: zs(:) + real(kind_phys), intent(in) :: dzs(:) + real(kind_phys), intent(inout) :: tvxy(:) + real(kind_phys), intent(inout) :: tgxy(:) + real(kind_phys), intent(inout) :: tahxy(:) + real(kind_phys), intent(inout) :: canicexy(:) + real(kind_phys), intent(inout) :: canliqxy(:) + real(kind_phys), intent(inout) :: eahxy(:) + real(kind_phys), intent(inout) :: cmxy(:) + real(kind_phys), intent(inout) :: chxy(:) + real(kind_phys), intent(inout) :: fwetxy(:) + real(kind_phys), intent(inout) :: sneqvoxy(:) + real(kind_phys), intent(inout) :: alboldxy(:) + real(kind_phys), intent(inout) :: qsnowxy(:) + real(kind_phys), intent(inout) :: wslakexy(:) + real(kind_phys), intent(inout) :: albdvis(:) + real(kind_phys), intent(inout) :: albdnir(:) + real(kind_phys), intent(inout) :: albivis(:) + real(kind_phys), intent(inout) :: albinir(:) + real(kind_phys), intent(inout) :: emiss(:) + real(kind_phys), intent(inout) :: taussxy(:) + real(kind_phys), intent(inout) :: waxy(:) + real(kind_phys), intent(inout) :: wtxy(:) + real(kind_phys), intent(inout) :: zwtxy(:) + real(kind_phys), intent(inout) :: xlaixy(:) + real(kind_phys), intent(inout) :: xsaixy(:) + real(kind_phys), intent(inout) :: lfmassxy(:) + real(kind_phys), intent(inout) :: stmassxy(:) + real(kind_phys), intent(inout) :: rtmassxy(:) + real(kind_phys), intent(inout) :: woodxy(:) + real(kind_phys), intent(inout) :: stblcpxy(:) + real(kind_phys), intent(inout) :: fastcpxy(:) + real(kind_phys), intent(inout) :: smcwtdxy(:) + real(kind_phys), intent(inout) :: deeprechxy(:) + real(kind_phys), intent(inout) :: rechxy(:) + real(kind_phys), intent(inout) :: snowxy(:) + real(kind_phys), intent(inout) :: snicexy(:,lsnow_lsm_lbound:) + real(kind_phys), intent(inout) :: snliqxy(:,lsnow_lsm_lbound:) + real(kind_phys), intent(inout) :: tsnoxy (:,lsnow_lsm_lbound:) + real(kind_phys), intent(inout) :: smoiseq(:,:) + real(kind_phys), intent(inout) :: zsnsoxy(:,lsnow_lsm_lbound:) + real(kind_phys), intent(inout) :: slc(:,:) + real(kind_phys), intent(inout) :: smc(:,:) + real(kind_phys), intent(inout) :: stc(:,:) + real(kind_phys), intent(in) :: tsfcl(:) + real(kind_phys), intent(in) :: snowd(:) + real(kind_phys), intent(in) :: canopy(:) + real(kind_phys), intent(in) :: tg3(:) + real(kind_phys), intent(in) :: stype(:) + real(kind_phys), intent(in) :: con_t0c + + integer, intent(in) :: nthrds character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! Local variables - integer :: i, j, ix + integer :: i, j, ix, vegtyp + real(kind_phys) :: rsnow + + !--- Noah MP + integer :: soiltyp, isnow, is, imn + real(kind=kind_phys) :: masslai, masssai, snd + real(kind=kind_phys) :: bexp, ddz, smcmax, smcwlt, dwsat, dksat, psisat + + real(kind=kind_phys), dimension(:), allocatable :: dzsno + real(kind=kind_phys), dimension(:), allocatable :: dzsnso ! Initialize CCPP error handling variables errmsg = '' @@ -153,6 +246,14 @@ subroutine GFS_phys_time_vary_init ( ! hardcoded in module iccn_def.F and GFS_typedefs.F90 endif +!> - Call tau_amf dats for ugwp_v1 + if (do_ugwp_v1) then + call read_tau_amf(me, master, errmsg, errflg) + endif + +!> - Initialize soil vegetation (needed for sncovr calculation further down) + call set_soilveg(me, isot, ivegsrc, nlunit) + !> - Call setindxoz() to initialize ozone data if (ntoz > 0) then call setindxoz (im, xlat_d, jindx1_o3, jindx2_o3, ddy_o3) @@ -178,6 +279,12 @@ subroutine GFS_phys_time_vary_init ( iindx1_ci, iindx2_ci, ddx_ci) endif +!> - Call cires_indx_ugwp to read monthly-mean GW-tau diagnosed from FV3GFS-runs that can resolve GWs + if (do_ugwp_v1) then + call cires_indx_ugwp (im, me, master, xlat_d, jindx1_tau, jindx2_tau, & + ddy_j1tau, ddy_j2tau) + endif + !--- initial calculation of maps local ix -> global i and j ix = 0 do j = 1,ny @@ -188,32 +295,319 @@ subroutine GFS_phys_time_vary_init ( enddo enddo -#if 0 - !Calculate sncovr if it was read in but empty (from FV3/io/FV3GFS_io.F90/sfc_prop_restart_read) - ! if (first_time_step) then - ! if (nint(Sfcprop%sncovr(1)) == -9999) then - ! !--- compute sncovr from existing variables - ! !--- code taken directly from read_fix.f - ! do ix = 1, im - ! Sfcprop%sncovr(ix) = 0.0 - ! if (Sfcprop%slmsk(ix) > 0.001) then - ! vegtyp = Sfcprop%vtype(ix) - ! if (vegtyp == 0) vegtyp = 7 - ! rsnow = 0.001*Sfcprop%weasd(ix)/snupx(vegtyp) - ! if (0.001*Sfcprop%weasd(ix) < snupx(vegtyp)) then - ! Sfcprop%sncovr(ix) = 1.0 - (exp(-salp_data*rsnow) - rsnow*exp(-salp_data)) - ! else - ! Sfcprop%sncovr(ix) = 1.0 - ! endif - ! endif - ! enddo - ! ! DH* 20201104: don't forget snocvr_ice for RUC LSM (see FV3GFS_io.F90) - ! endif - ! endif -#endif + !--- if sncovr does not exist in the restart, need to create it + if (all(sncovr < zero)) then + if (me == master ) write(0,'(a)') 'GFS_phys_time_vary_init: compute sncovr from weasd and soil vegetation parameters' + !--- compute sncovr from existing variables + !--- code taken directly from read_fix.f + sncovr(:) = zero + do ix=1,im + if (landfrac(ix) >= drythresh .or. fice(ix) >= min_seaice) then + vegtyp = vtype(ix) + if (vegtyp == 0) vegtyp = 7 + rsnow = 0.001_kind_phys*weasd(ix)/snupx(vegtyp) + if (0.001_kind_phys*weasd(ix) < snupx(vegtyp)) then + sncovr(ix) = one - (exp(-salp_data*rsnow) - rsnow*exp(-salp_data)) + else + sncovr(ix) = one + endif + endif + enddo + endif + + !--- For RUC LSM: create sncovr_ice from sncovr + if (lsm == lsm_ruc) then + if (all(sncovr_ice < zero)) then + if (me == master ) write(0,'(a)') 'GFS_phys_time_vary_init: fill sncovr_ice with sncovr for RUC LSM' + sncovr_ice(:) = sncovr(:) + endif + endif + + if (lsm == lsm_noahmp) then + if (all(tvxy < zero)) then + + allocate(dzsno (lsnow_lsm_lbound:lsnow_lsm_ubound)) + allocate(dzsnso(lsnow_lsm_lbound:lsoil) ) + dzsno(:) = missing_value + dzsnso(:) = missing_value + + tvxy(:) = missing_value + tgxy(:) = missing_value + tahxy(:) = missing_value + canicexy(:) = missing_value + canliqxy(:) = missing_value + eahxy(:) = missing_value + cmxy(:) = missing_value + chxy(:) = missing_value + fwetxy(:) = missing_value + sneqvoxy(:) = missing_value + alboldxy(:) = missing_value + qsnowxy(:) = missing_value + wslakexy(:) = missing_value + albdvis(:) = missing_value + albdnir(:) = missing_value + albivis(:) = missing_value + albinir(:) = missing_value + emiss(:) = missing_value + taussxy(:) = missing_value + waxy(:) = missing_value + wtxy(:) = missing_value + zwtxy(:) = missing_value + xlaixy(:) = missing_value + xsaixy(:) = missing_value + + lfmassxy(:) = missing_value + stmassxy(:) = missing_value + rtmassxy(:) = missing_value + woodxy(:) = missing_value + stblcpxy(:) = missing_value + fastcpxy(:) = missing_value + smcwtdxy(:) = missing_value + deeprechxy(:) = missing_value + rechxy(:) = missing_value + + snowxy (:) = missing_value + snicexy(:,:) = missing_value + snliqxy(:,:) = missing_value + tsnoxy (:,:) = missing_value + smoiseq(:,:) = missing_value + zsnsoxy(:,:) = missing_value + + do ix=1,im + if (landfrac(ix) >= drythresh) then + tvxy(ix) = tsfcl(ix) + tgxy(ix) = tsfcl(ix) + tahxy(ix) = tsfcl(ix) + + if (snowd(ix) > 0.01_kind_phys .and. tsfcl(ix) > con_t0c ) tvxy(ix) = con_t0c + if (snowd(ix) > 0.01_kind_phys .and. tsfcl(ix) > con_t0c ) tgxy(ix) = con_t0c + if (snowd(ix) > 0.01_kind_phys .and. tsfcl(ix) > con_t0c ) tahxy(ix) = con_t0c + + canicexy(ix) = 0.0_kind_phys + canliqxy(ix) = canopy(ix) + + eahxy(ix) = 2000.0_kind_phys + +! eahxy = psfc*qv/(0.622+qv); qv is mixing ratio, converted from sepcific +! humidity specific humidity /(1.0 - specific humidity) + + cmxy(ix) = zero + chxy(ix) = zero + fwetxy(ix) = zero + sneqvoxy(ix) = weasd(ix) ! mm + alboldxy(ix) = 0.65_kind_phys + qsnowxy(ix) = zero + +! if (srflag(ix) > 0.001) qsnowxy(ix) = tprcp(ix)/dtp + ! already set to 0.0 + wslakexy(ix) = zero + taussxy(ix) = zero + albdvis(ix) = 0.2_kind_phys + albdnir(ix) = 0.2_kind_phys + albivis(ix) = 0.2_kind_phys + albinir(ix) = 0.2_kind_phys + emiss(ix) = 0.95_kind_phys + + + waxy(ix) = 4900.0_kind_phys + wtxy(ix) = waxy(ix) + zwtxy(ix) = (25.0_kind_phys + 2.0_kind_phys) - waxy(ix) / 1000.0_kind_phys / 0.2_kind_phys + + vegtyp = vtype(ix) + if (vegtyp == 0) vegtyp = 7 + imn = idate(2) + + if ((vegtyp == isbarren_table) .or. (vegtyp == isice_table) .or. (vegtyp == isurban_table) .or. (vegtyp == iswater_table)) then + + xlaixy(ix) = zero + xsaixy(ix) = zero + + lfmassxy(ix) = zero + stmassxy(ix) = zero + rtmassxy(ix) = zero + + woodxy (ix) = zero + stblcpxy (ix) = zero + fastcpxy (ix) = zero + + else + + xlaixy(ix) = max(laim_table(vegtyp, imn),0.05_kind_phys) +! xsaixy(ix) = max(saim_table(vegtyp, imn),0.05) + xsaixy(ix) = max(xlaixy(ix)*0.1_kind_phys,0.05_kind_phys) + + masslai = 1000.0_kind_phys / max(sla_table(vegtyp),one) + lfmassxy(ix) = xlaixy(ix)*masslai + masssai = 1000.0_kind_phys / 3.0_kind_phys + stmassxy(ix) = xsaixy(ix)* masssai + + rtmassxy(ix) = 500.0_kind_phys + + woodxy(ix) = 500.0_kind_phys + stblcpxy(ix) = 1000.0_kind_phys + fastcpxy(ix) = 1000.0_kind_phys + + endif ! non urban ... + + if (vegtyp == isice_table) then + do is = 1,lsoil + stc(ix,is) = min(stc(ix,is),min(tg3(ix),263.15_kind_phys)) + smc(ix,is) = one + slc(ix,is) = zero + enddo + endif + + snd = snowd(ix)/1000.0_kind_phys ! go to m from snwdph + + if (weasd(ix) /= zero .and. snd == zero ) then + snd = weasd(ix)/1000.0 + endif + + if (vegtyp == 15) then ! land ice in MODIS/IGBP + if (weasd(ix) < 0.1_kind_phys) then + weasd(ix) = 0.1_kind_phys + snd = 0.01_kind_phys + endif + endif + + if (snd < 0.025_kind_phys ) then + snowxy(ix) = zero + dzsno(-2:0) = zero + elseif (snd >= 0.025_kind_phys .and. snd <= 0.05_kind_phys ) then + snowxy(ix) = -1.0_kind_phys + dzsno(0) = snd + elseif (snd > 0.05_kind_phys .and. snd <= 0.10_kind_phys ) then + snowxy(ix) = -2.0_kind_phys + dzsno(-1) = 0.5_kind_phys*snd + dzsno(0) = 0.5_kind_phys*snd + elseif (snd > 0.10_kind_phys .and. snd <= 0.25_kind_phys ) then + snowxy(ix) = -2.0_kind_phys + dzsno(-1) = 0.05_kind_phys + dzsno(0) = snd - 0.05_kind_phys + elseif (snd > 0.25_kind_phys .and. snd <= 0.45_kind_phys ) then + snowxy(ix) = -3.0_kind_phys + dzsno(-2) = 0.05_kind_phys + dzsno(-1) = 0.5_kind_phys*(snd-0.05_kind_phys) + dzsno(0) = 0.5_kind_phys*(snd-0.05_kind_phys) + elseif (snd > 0.45_kind_phys) then + snowxy(ix) = -3.0_kind_phys + dzsno(-2) = 0.05_kind_phys + dzsno(-1) = 0.20_kind_phys + dzsno(0) = snd - 0.05_kind_phys - 0.20_kind_phys + else + errmsg = 'Error in GFS_phys_time_vary.fv3.F90: Problem with the logic assigning snow layers in Noah MP initialization' + errflg = 1 + return + endif + +! Now we have the snowxy field +! snice + snliq + tsno allocation and compute them from what we have + + tsnoxy(ix,:) = zero + snicexy(ix,:) = zero + snliqxy(ix,:) = zero + zsnsoxy(ix,:) = zero + + isnow = nint(snowxy(ix))+1 ! snowxy <=0.0, dzsno >= 0.0 + + do is = isnow,0 + tsnoxy(ix,is) = tgxy(ix) + snliqxy(ix,is) = zero + snicexy(ix,is) = one * dzsno(is) * weasd(ix)/snd + enddo +! +!zsnsoxy, all negative ? +! + do is = isnow,0 + dzsnso(is) = -dzsno(is) + enddo + + do is = 1,4 + dzsnso(is) = -dzs(is) + enddo +! +! Assign to zsnsoxy +! + zsnsoxy(ix,isnow) = dzsnso(isnow) + do is = isnow+1,4 + zsnsoxy(ix,is) = zsnsoxy(ix,is-1) + dzsnso(is) + enddo +! +! smoiseq +! Init water table related quantities here +! + soiltyp = stype(ix) + if (soiltyp /= 0) then + bexp = bexp_table(soiltyp) + smcmax = smcmax_table(soiltyp) + smcwlt = smcwlt_table(soiltyp) + dwsat = dwsat_table(soiltyp) + dksat = dksat_table(soiltyp) + psisat = -psisat_table(soiltyp) + endif + + if (vegtyp == isurban_table) then + smcmax = 0.45_kind_phys + smcwlt = 0.40_kind_phys + endif + + if ((bexp > zero) .and. (smcmax > zero) .and. (-psisat > zero)) then + do is = 1, lsoil + if ( is == 1 )then + ddz = -zs(is+1) * 0.5_kind_phys + elseif ( is < lsoil ) then + ddz = ( zs(is-1) - zs(is+1) ) * 0.5_kind_phys + else + ddz = zs(is-1) - zs(is) + endif + smoiseq(ix,is) = min(max(find_eq_smc(bexp, dwsat, dksat, ddz, smcmax),1.e-4_kind_phys),smcmax*0.99_kind_phys) + enddo + else ! bexp <= 0.0 + smoiseq(ix,1:4) = smcmax + endif ! end the bexp condition + + smcwtdxy(ix) = smcmax + deeprechxy(ix) = zero + rechxy(ix) = zero + + endif + + enddo ! ix + + deallocate(dzsno) + deallocate(dzsnso) + + endif + endif !if Noah MP cold start ends is_initialized = .true. + contains + +! +! Use newton-raphson method to find eq soil moisture +! + function find_eq_smc(bexp, dwsat, dksat, ddz, smcmax) result(smc) + implicit none + real(kind=kind_phys), intent(in) :: bexp, dwsat, dksat, ddz, smcmax + real(kind=kind_phys) :: smc + real(kind=kind_phys) :: expon, aa, bb, func, dfunc, dx + integer :: iter + ! + expon = bexp + 1. + aa = dwsat / ddz + bb = dksat / smcmax ** expon + smc = 0.5 * smcmax + ! + do iter = 1,100 + func = (smc - smcmax) * aa + bb * smc ** expon + dfunc = aa + bb * expon * smc ** bexp + dx = func / dfunc + smc = smc - dx + if ( abs (dx) < 1.e-6_kind_phys) return + enddo + end function find_eq_smc + end subroutine GFS_phys_time_vary_init !! @} @@ -228,7 +622,8 @@ subroutine GFS_phys_time_vary_timestep_init ( jindx1_o3, jindx2_o3, ddy_o3, ozpl, jindx1_h, jindx2_h, ddy_h, h2opl, & jindx1_aer, jindx2_aer, ddy_aer, iindx1_aer, iindx2_aer, ddx_aer, aer_nm, & jindx1_ci, jindx2_ci, ddy_ci, iindx1_ci, iindx2_ci, ddx_ci, in_nm, ccn_nm, & - imap, jmap, prsl, seed0, rann, errmsg, errflg) + imap, jmap, prsl, seed0, rann, do_ugwp_v1, jindx1_tau, jindx2_tau, ddy_j1tau, ddy_j2tau,& + tau_amf, errmsg, errflg) implicit none @@ -252,7 +647,11 @@ subroutine GFS_phys_time_vary_timestep_init ( real(kind_phys), intent(in) :: prsl(:,:) integer, intent(in) :: seed0 real(kind_phys), intent(inout) :: rann(:,:) - ! + + logical, intent(in) :: do_ugwp_v1 + integer, intent(in) :: jindx1_tau(:), jindx2_tau(:) + real(kind_phys), intent(in) :: ddy_j1tau(:), ddy_j2tau(:) + real(kind_phys), intent(inout) :: tau_amf(:) character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -347,6 +746,13 @@ subroutine GFS_phys_time_vary_timestep_init ( levs, prsl, in_nm, ccn_nm) endif +!> - Call cires_indx_ugwp to read monthly-mean GW-tau diagnosed from FV3GFS-runs that resolve GW-activ + if (do_ugwp_v1) then + call tau_amf_interp(me, master, im, idate, fhour, & + jindx1_tau, jindx2_tau, & + ddy_j1tau, ddy_j2tau, tau_amf) + endif + ! Not needed for SCM: !> - Call gcycle() to repopulate specific time-varying surface properties for AMIP/forecast runs !if (nscyc > 0) then @@ -423,6 +829,11 @@ subroutine GFS_phys_time_vary_finalize(errmsg, errflg) if (allocated(ccnin) ) deallocate(ccnin) if (allocated(ci_pres) ) deallocate(ci_pres) + ! Deallocate UGWP-input arrays + if (allocated(ugwp_taulat)) deallocate(ugwp_taulat) + if (allocated(tau_limb )) deallocate(tau_limb) + if (allocated(days_limb )) deallocate(days_limb) + is_initialized = .false. end subroutine GFS_phys_time_vary_finalize diff --git a/physics/GFS_phys_time_vary.scm.meta b/physics/GFS_phys_time_vary.scm.meta index cf0b3afbd..1edaa32c8 100644 --- a/physics/GFS_phys_time_vary.scm.meta +++ b/physics/GFS_phys_time_vary.scm.meta @@ -1,7 +1,8 @@ [ccpp-table-properties] name = GFS_phys_time_vary type = scheme - dependencies = aerclm_def.F,aerinterp.F90,h2o_def.f,h2ointerp.f90,iccn_def.F,iccninterp.F90,machine.F,mersenne_twister.f,namelist_soilveg.f,ozinterp.f90,ozne_def.f + dependencies = aerclm_def.F,aerinterp.F90,h2o_def.f,h2ointerp.f90,iccn_def.F,iccninterp.F90,machine.F,mersenne_twister.f + dependencies = namelist_soilveg.f,set_soilveg.f,ozinterp.f90,ozne_def.f,cires_tauamf_data.F90,noahmp_tables.f90 ######################################################################## [ccpp-arg-table] @@ -307,6 +308,641 @@ type = integer intent = inout optional = F +[do_ugwp_v1] + standard_name = flag_for_ugwp_version_1 + long_name = flag to activate ver 1 CIRES UGWP + units = flag + dimensions = () + type = logical + intent = in + optional = F +[jindx1_tau] + standard_name = lower_latitude_index_of_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag_for_interpolation + long_name = index1 for weight1 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = integer + intent = inout + optional = F +[jindx2_tau] + standard_name = upper_latitude_index_of_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag_for_interpolation + long_name = index2 for weight2 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = integer + intent = inout + optional = F +[ddy_j1tau] + standard_name = latitude_interpolation_weight_complement_for_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag + long_name = interpolation weight1 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = real + intent = inout + kind = kind_phys + optional = F +[ddy_j2tau] + standard_name = latitude_interpolation_weight_for_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag + long_name = interpolation weight2 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = real + intent = inout + kind = kind_phys + optional = F +[isot] + standard_name = soil_type_dataset_choice + long_name = soil type dataset choice + units = index + dimensions = () + type = integer + intent = in + optional = F +[ivegsrc] + standard_name = vegetation_type_dataset_choice + long_name = land use dataset choice + units = index + dimensions = () + type = integer + intent = in + optional = F +[nlunit] + standard_name = iounit_namelist + long_name = fortran unit number for file opens + units = none + dimensions = () + type = integer + intent = in + optional = F +[sncovr] + standard_name = surface_snow_area_fraction_over_land + long_name = surface snow area fraction + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[sncovr_ice] + standard_name = surface_snow_area_fraction_over_ice + long_name = surface snow area fraction over ice + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[lsm] + standard_name = flag_for_land_surface_scheme + long_name = flag for land surface model + units = flag + dimensions = () + type = integer + intent = in + optional = F +[lsm_noahmp] + standard_name = flag_for_noahmp_land_surface_scheme + long_name = flag for NOAH MP land surface model + units = flag + dimensions = () + type = integer + intent = in + optional = F +[lsm_ruc] + standard_name = flag_for_ruc_land_surface_scheme + long_name = flag for RUC land surface model + units = flag + dimensions = () + type = integer + intent = in + optional = F +[min_seaice] + standard_name = sea_ice_minimum + long_name = minimum sea ice value + units = frac + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[fice] + standard_name = sea_ice_concentration + long_name = ice fraction over open water + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[landfrac] + standard_name = land_area_fraction + long_name = fraction of horizontal grid area occupied by land + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[vtype] + standard_name = vegetation_type_classification_real + long_name = vegetation type for lsm + units = index + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[weasd] + standard_name = water_equivalent_accumulated_snow_depth + long_name = water equiv of acc snow depth over land and sea ice + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[lsoil] + standard_name = soil_vertical_dimension + long_name = number of soil layers + units = count + dimensions = () + type = integer + intent = in + optional = F +[zs] + standard_name = depth_of_soil_levels_for_land_surface_model + long_name = depth of soil levels for land surface model + units = m + dimensions = (soil_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = in + optional = F +[dzs] + standard_name = thickness_of_soil_levels_for_land_surface_model + long_name = thickness of soil levels for land surface model + units = m + dimensions = (soil_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = in + optional = F +[lsnow_lsm_lbound] + standard_name = lower_bound_of_snow_vertical_dimension_for_land_surface_model + long_name = lower bound of of snow-related arrays for land surface model + units = count + dimensions = () + type = integer + intent = in + optional = F +[lsnow_lsm_ubound] + standard_name = upper_bound_of_snow_vertical_dimension_for_land_surface_model + long_name = upper bound of of snow-related arrays for land surface model + units = count + dimensions = () + type = integer + intent = in + optional = F +[tvxy] + standard_name = vegetation_temperature + long_name = vegetation temperature + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[tgxy] + standard_name = ground_temperature_for_noahmp + long_name = ground temperature for noahmp + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[tahxy] + standard_name = canopy_air_temperature + long_name = canopy air temperature + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[canicexy] + standard_name = canopy_intercepted_ice_mass + long_name = canopy intercepted ice mass + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[canliqxy] + standard_name = canopy_intercepted_liquid_water + long_name = canopy intercepted liquid water + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[eahxy] + standard_name = canopy_air_vapor_pressure + long_name = canopy air vapor pressure + units = Pa + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[cmxy] + standard_name = surface_drag_coefficient_for_momentum_for_noahmp + long_name = surface drag coefficient for momentum for noahmp + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[chxy] + standard_name = surface_drag_coefficient_for_heat_and_moisture_for_noahmp + long_name = surface exchange coeff heat & moisture for noahmp + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fwetxy] + standard_name = area_fraction_of_wet_canopy + long_name = area fraction of canopy that is wetted/snowed + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[sneqvoxy] + standard_name = snow_mass_at_previous_time_step + long_name = snow mass at previous time step + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[alboldxy] + standard_name = snow_albedo_at_previous_time_step + long_name = snow albedo at previous time step + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[qsnowxy] + standard_name = snow_precipitation_rate_at_surface + long_name = snow precipitation rate at surface + units = mm s-1 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[wslakexy] + standard_name = lake_water_storage + long_name = lake water storage + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[taussxy] + standard_name = nondimensional_snow_age + long_name = non-dimensional snow age + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[waxy] + standard_name = water_storage_in_aquifer + long_name = water storage in aquifer + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[wtxy] + standard_name = water_storage_in_aquifer_and_saturated_soil + long_name = water storage in aquifer and saturated soil + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[zwtxy] + standard_name = water_table_depth + long_name = water table depth + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[xlaixy] + standard_name = leaf_area_index + long_name = leaf area index + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[xsaixy] + standard_name = stem_area_index + long_name = stem area index + units = none + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[lfmassxy] + standard_name = leaf_mass + long_name = leaf mass + units = g m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[stmassxy] + standard_name = stem_mass + long_name = stem mass + units = g m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[rtmassxy] + standard_name = fine_root_mass + long_name = fine root mass + units = g m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[woodxy] + standard_name = wood_mass + long_name = wood mass including woody roots + units = g m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[stblcpxy] + standard_name = slow_soil_pool_mass_content_of_carbon + long_name = stable carbon in deep soil + units = g m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[fastcpxy] + standard_name = fast_soil_pool_mass_content_of_carbon + long_name = short-lived carbon in shallow soil + units = g m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[smcwtdxy] + standard_name = soil_water_content_between_soil_bottom_and_water_table + long_name = soil water content between the bottom of the soil and the water table + units = m3 m-3 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[deeprechxy] + standard_name = water_table_recharge_when_deep + long_name = recharge to or from the water table when deep + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[rechxy] + standard_name = water_table_recharge_when_shallow + long_name = recharge to or from the water table when shallow + units = m + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[albdvis] + standard_name = surface_albedo_direct_visible + long_name = direct surface albedo visible band + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[albdnir] + standard_name = surface_albedo_direct_NIR + long_name = direct surface albedo NIR band + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[albivis] + standard_name = surface_albedo_diffuse_visible + long_name = diffuse surface albedo visible band + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[albinir] + standard_name = surface_albedo_diffuse_NIR + long_name = diffuse surface albedo NIR band + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[emiss] + standard_name = surface_emissivity_lsm + long_name = surface emissivity from lsm + units = frac + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[snowxy] + standard_name = number_of_snow_layers + long_name = number of snow layers + units = count + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[snicexy] + standard_name = snow_layer_ice + long_name = snow layer ice + units = mm + dimensions = (horizontal_dimension,lower_bound_of_snow_vertical_dimension_for_land_surface_model:upper_bound_of_snow_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = inout + optional = F +[snliqxy] + standard_name = snow_layer_liquid_water + long_name = snow layer liquid water + units = mm + dimensions = (horizontal_dimension,lower_bound_of_snow_vertical_dimension_for_land_surface_model:upper_bound_of_snow_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = inout + optional = F +[tsnoxy] + standard_name = snow_temperature + long_name = snow_temperature + units = K + dimensions = (horizontal_dimension,lower_bound_of_snow_vertical_dimension_for_land_surface_model:upper_bound_of_snow_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = inout + optional = F +[smoiseq] + standard_name = equilibrium_soil_water_content + long_name = equilibrium soil water content + units = m3 m-3 + dimensions = (horizontal_dimension,soil_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = inout + optional = F +[zsnsoxy] + standard_name = layer_bottom_depth_from_snow_surface + long_name = depth from the top of the snow surface at the bottom of the layer + units = m + dimensions = (horizontal_dimension,lower_bound_of_snow_vertical_dimension_for_land_surface_model:soil_vertical_dimension_for_land_surface_model) + type = real + kind = kind_phys + intent = inout + optional = F +[slc] + standard_name = volume_fraction_of_unfrozen_soil_moisture + long_name = liquid soil moisture + units = frac + dimensions = (horizontal_dimension,soil_vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[smc] + standard_name = volume_fraction_of_soil_moisture + long_name = total soil moisture + units = frac + dimensions = (horizontal_dimension,soil_vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[stc] + standard_name = soil_temperature + long_name = soil temperature + units = K + dimensions = (horizontal_dimension,soil_vertical_dimension) + type = real + kind = kind_phys + intent = inout + optional = F +[tsfcl] + standard_name = surface_skin_temperature_over_land + long_name = surface skin temperature over land + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[snowd] + standard_name = surface_snow_thickness_water_equivalent + long_name = water equivalent snow depth + units = mm + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[canopy] + standard_name = canopy_water_amount + long_name = canopy water amount + units = kg m-2 + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[tg3] + standard_name = deep_soil_temperature + long_name = deep soil temperature + units = K + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[stype] + standard_name = soil_type_classification_real + long_name = soil type for lsm + units = index + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = in + optional = F +[con_t0c] + standard_name = temperature_at_zero_celsius + long_name = temperature at 0 degree Celsius + units = K + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F +[nthrds] + standard_name = omp_threads + long_name = number of OpenMP threads available for physics schemes + units = count + dimensions = () + type = integer + intent = in + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP @@ -442,7 +1078,7 @@ [nsswr] standard_name = number_of_timesteps_between_shortwave_radiation_calls long_name = number of timesteps between shortwave radiation calls - units = + units = dimensions = () type = integer intent = in @@ -775,6 +1411,57 @@ kind = kind_phys intent = inout optional = F +[do_ugwp_v1] + standard_name = flag_for_ugwp_version_1 + long_name = flag to activate ver 1 CIRES UGWP + units = flag + dimensions = () + type = logical + intent = in + optional = F +[jindx1_tau] + standard_name = lower_latitude_index_of_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag_for_interpolation + long_name = index1 for weight1 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = integer + intent = in + optional = F +[jindx2_tau] + standard_name = upper_latitude_index_of_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag_for_interpolation + long_name = index2 for weight2 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = integer + intent = in + optional = F +[ddy_j1tau] + standard_name = latitude_interpolation_weight_complement_for_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag + long_name = interpolation weight1 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = real + intent = in + kind = kind_phys + optional = F +[ddy_j2tau] + standard_name = latitude_interpolation_weight_for_absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag + long_name = interpolation weight2 for tau NGWs + units = none + dimensions = (horizontal_dimension) + type = real + intent = in + kind = kind_phys + optional = F +[tau_amf] + standard_name = absolute_momentum_flux_due_to_nonorographic_gravity_wave_drag + long_name = ngw_absolute_momentum_flux + units = various + dimensions = (horizontal_dimension) + type = real + kind = kind_phys + intent = inout + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP @@ -813,4 +1500,4 @@ dimensions = () type = integer intent = out - optional = F \ No newline at end of file + optional = F From 4d8f7f4117d5695c534580bdf6350c39eb95c2ad Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 5 Mar 2021 17:30:51 -0700 Subject: [PATCH 18/31] Bugfix of my own bugfix in physics/module_mp_thompson.F90: always initialize cloud effective radii --- physics/module_mp_thompson.F90 | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index 8ba269388..dfe31f375 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -5204,6 +5204,10 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & has_qi = .false. has_qs = .false. + re_qc1d(:) = 0.0D0 + re_qi1d(:) = 0.0D0 + re_qs1d(:) = 0.0D0 + 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)) @@ -5230,8 +5234,6 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & lamc = (nc(k)*am_r*g_ratio(inu_c)/rc(k))**obmr re_qc1d(k) = SNGL(0.5D0 * DBLE(3.+inu_c)/lamc) enddo - else - re_qc1d(:) = 0.0D0 endif if (has_qi) then @@ -5240,8 +5242,6 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & lami = (am_i*cig(2)*oig1*ni(k)/ri(k))**obmi re_qi1d(k) = SNGL(0.5D0 * DBLE(3.+mu_i)/lami) enddo - else - re_qi1d(:) = 0.0D0 endif if (has_qs) then @@ -5282,8 +5282,6 @@ subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & smoc = a_ * smo2**b_ re_qs1d(k) = 0.5*(smoc/smob) enddo - else - re_qs1d(:) = 0.0D0 endif end subroutine calc_effectRad From 22f45dde5467e51886c9bd622ee885e968ca17f0 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Tue, 9 Mar 2021 12:05:43 -0700 Subject: [PATCH 19/31] modify check for missing NoahMP input to initiate cold start in GFS_phys_time_vary.scm.F90 --- physics/GFS_phys_time_vary.scm.F90 | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/physics/GFS_phys_time_vary.scm.F90 b/physics/GFS_phys_time_vary.scm.F90 index dabd53e94..a54ffa7a9 100644 --- a/physics/GFS_phys_time_vary.scm.F90 +++ b/physics/GFS_phys_time_vary.scm.F90 @@ -324,8 +324,7 @@ subroutine GFS_phys_time_vary_init ( endif if (lsm == lsm_noahmp) then - if (all(tvxy < zero)) then - + if (all(tvxy <= zero)) then allocate(dzsno (lsnow_lsm_lbound:lsnow_lsm_ubound)) allocate(dzsnso(lsnow_lsm_lbound:lsoil) ) dzsno(:) = missing_value From fd8f925177daf5811e6998e1a48588611fee9191 Mon Sep 17 00:00:00 2001 From: Grant Firl Date: Tue, 9 Mar 2021 12:20:35 -0700 Subject: [PATCH 20/31] make cires_ugwpv1_module.F90/cires_ugwp1_init use the passed-in namelist --- physics/cires_ugwpv1_module.F90 | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/physics/cires_ugwpv1_module.F90 b/physics/cires_ugwpv1_module.F90 index 13b7752a5..dbbd3cd2c 100644 --- a/physics/cires_ugwpv1_module.F90 +++ b/physics/cires_ugwpv1_module.F90 @@ -177,7 +177,7 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & real(kind=kind_phys), intent (in) :: con_pi, con_rerth character(len=64), intent (in) :: fn_nml2 - character(len=64), parameter :: fn_nml='input.nml' +! character(len=64), parameter :: fn_nml='input.nml' character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -194,15 +194,15 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & ! - if (me == master) print *, trim (fn_nml), ' GW-namelist file ' - inquire (file =trim (fn_nml) , exist = exists) + if (me == master) print *, trim (fn_nml2), ' GW-namelist file ' + inquire (file =trim (fn_nml2) , exist = exists) ! if (.not. exists) then if (me == master) & - write (6, *) 'separate ugwp :: namelist file: ', trim (fn_nml), ' does not exist' + write (6, *) 'separate ugwp :: namelist file: ', trim (fn_nml2), ' does not exist' else - open (unit = nlunit, file = trim(fn_nml), action = 'read', status = 'old', iostat = ios) + open (unit = nlunit, file = trim(fn_nml2), action = 'read', status = 'old', iostat = ios) endif rewind (nlunit) read (nlunit, nml = cires_ugwp_nml) From f5733901180d8fe6919425c47b56aef5b7c44939 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 12 Mar 2021 16:38:07 -0700 Subject: [PATCH 21/31] physics/radlw_main.F90: apply local bounds for cloud effective radii instead of aborting the model run --- physics/radlw_main.F90 | 66 +++++++++++++++++++++--------------------- 1 file changed, 33 insertions(+), 33 deletions(-) diff --git a/physics/radlw_main.F90 b/physics/radlw_main.F90 index de8d9e973..7655e76d2 100644 --- a/physics/radlw_main.F90 +++ b/physics/radlw_main.F90 @@ -1250,7 +1250,7 @@ subroutine rrtmg_lw_run & endif !mz* HWRF: calculate taucmc with mcica - if (iovr == 4) then + if (iovr == 4) then call cldprmc(nlay, inflglw, iceflglw, liqflglw, & & cldfmc, ciwpmc, & & clwpmc, cswpmc, reicmc, relqmc, resnmc, & @@ -8854,25 +8854,25 @@ subroutine cldprmc(nlayers, inflag, iceflag, liqflag, cldfmc, & abscosno(ig) = 0.0_rb elseif (iceflag .eq. 0) then - if (radice .lt. 10.0_rb) stop 'ICE RADIUS TOO SMALL' - abscoice(ig) = absice0(1) + absice0(2)/radice +! if (radice .lt. 10.0_rb) stop 'ICE RADIUS TOO SMALL' + abscoice(ig) = absice0(1) + absice0(2)/max(radice,10.0_rb) abscosno(ig) = 0.0_rb elseif (iceflag .eq. 1) then - if (radice .lt. 13.0_rb .or. radice .gt. 130._rb) stop& - & 'ICE RADIUS OUT OF BOUNDS' +! if (radice .lt. 13.0_rb .or. radice .gt. 130._rb) stop& +! & 'ICE RADIUS OUT OF BOUNDS' ncbands = 5 ib = icb(ngb(ig)) - abscoice(ig) = absice1(1,ib) + absice1(2,ib)/radice + abscoice(ig) = absice1(1,ib) + absice1(2,ib)/min(max(radice,13.0_rb),130._rb) abscosno(ig) = 0.0_rb ! For iceflag=2 option, ice particle effective radius is limited to 5.0 to 131.0 microns elseif (iceflag .eq. 2) then - if (radice .lt. 5.0_rb .or. radice .gt. 131.0_rb) stop& - & 'ICE RADIUS OUT OF BOUNDS' +! if (radice .lt. 5.0_rb .or. radice .gt. 131.0_rb) stop& +! & 'ICE RADIUS OUT OF BOUNDS' ncbands = 16 - factor = (radice - 2._rb)/3._rb + factor = (min(max(radice,5.0_rb),131._rb) - 2._rb)/3._rb index = int(factor) if (index .eq. 43) index = 42 fint = factor - float(index) @@ -8885,15 +8885,15 @@ subroutine cldprmc(nlayers, inflag, iceflag, liqflag, cldfmc, & ! For iceflag=3 option, ice particle generalized effective size is limited to 5.0 to 140.0 microns elseif (iceflag .ge. 3) then - if (radice .lt. 5.0_rb .or. radice .gt. 140.0_rb) then - write(errmsg,'(a,i5,i5,f8.2,f8.2)' ) & - & 'ERROR: ICE GENERALIZED EFFECTIVE SIZE OUT OF BOUNDS' & - & ,ig, lay, ciwpmc(ig,lay), radice - errflg = 1 - return - end if +! if (radice .lt. 5.0_rb .or. radice .gt. 140.0_rb) then +! write(errmsg,'(a,i5,i5,f8.2,f8.2)' ) & +! & 'ERROR: ICE GENERALIZED EFFECTIVE SIZE OUT OF BOUNDS' & +! & ,ig, lay, ciwpmc(ig,lay), radice +! errflg = 1 +! return +! end if ncbands = 16 - factor = (radice - 2._rb)/3._rb + factor = (min(max(radice,5.0_rb),140._rb) - 2._rb)/3._rb index = int(factor) if (index .eq. 46) index = 45 fint = factor - float(index) @@ -8908,15 +8908,15 @@ subroutine cldprmc(nlayers, inflag, iceflag, liqflag, cldfmc, & !..Incorporate additional effects due to snow. if (cswpmc(ig,lay).gt.0.0_rb .and. iceflag .eq. 5) then radsno = resnmc(lay) - if (radsno .lt. 5.0_rb .or. radsno .gt. 140.0_rb) then - write(errmsg,'(a,i5,i5,f8.2,f8.2)' ) & - & 'ERROR: SNOW GENERALIZED EFFECTIVE SIZE OUT OF BOUNDS' & - & ,ig, lay, cswpmc(ig,lay), radsno - errflg = 1 - return - end if +! if (radsno .lt. 5.0_rb .or. radsno .gt. 140.0_rb) then +! write(errmsg,'(a,i5,i5,f8.2,f8.2)' ) & +! & 'ERROR: SNOW GENERALIZED EFFECTIVE SIZE OUT OF BOUNDS' & +! & ,ig, lay, cswpmc(ig,lay), radsno +! errflg = 1 +! return +! end if ncbands = 16 - factor = (radsno - 2._rb)/3._rb + factor = (min(max(radsno,5.0_rb),140.0_rb) - 2._rb)/3._rb index = int(factor) if (index .eq. 46) index = 45 fint = factor - float(index) @@ -8937,14 +8937,14 @@ subroutine cldprmc(nlayers, inflag, iceflag, liqflag, cldfmc, & elseif (liqflag .eq. 1) then radliq = relqmc(lay) - if (radliq .lt. 2.5_rb .or. radliq .gt. 60._rb) then - write(errmsg,'(a,i5,i5,f8.2,f8.2)' ) & -& 'ERROR: LIQUID EFFECTIVE SIZE OUT OF BOUNDS' & -& ,ig, lay, clwpmc(ig,lay), radliq - errflg = 1 - return - end if - index = int(radliq - 1.5_rb) +! if (radliq .lt. 2.5_rb .or. radliq .gt. 60._rb) then +! write(errmsg,'(a,i5,i5,f8.2,f8.2)' ) & +!& 'ERROR: LIQUID EFFECTIVE SIZE OUT OF BOUNDS' & +!& ,ig, lay, clwpmc(ig,lay), radliq +! errflg = 1 +! return +! end if + index = int(min(max(radliq,2.5_rb),60._rb) - 1.5_rb) if (index .eq. 0) index = 1 if (index .eq. 58) index = 57 fint = radliq - 1.5_rb - float(index) From c579871215ecfcc432ded8531f3a66defce4a9e6 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Fri, 12 Mar 2021 16:38:40 -0700 Subject: [PATCH 22/31] Bugfix for uninitialized variable in physics/module_bl_mynn.F90 --- physics/module_bl_mynn.F90 | 1 + 1 file changed, 1 insertion(+) diff --git a/physics/module_bl_mynn.F90 b/physics/module_bl_mynn.F90 index fa892eba8..d691de909 100644 --- a/physics/module_bl_mynn.F90 +++ b/physics/module_bl_mynn.F90 @@ -2947,6 +2947,7 @@ SUBROUTINE mynn_tendencies(kts,kte, & khdz(k) = rhoz(k)*dfh(k) kmdz(k) = rhoz(k)*dfm(k) ENDDO + rhoz(kte+1)=rhoz(kte) khdz(kte+1)=rhoz(kte+1)*dfh(kte) kmdz(kte+1)=rhoz(kte+1)*dfm(kte) From b23f06a4bcdafd5effc20a20333eaa125c160ace Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 15 Mar 2021 05:42:25 -0600 Subject: [PATCH 23/31] Bugfix in physics/module_sf_noahmplsm.f90 for uninitialized variable in subroutine surrad --- physics/module_sf_noahmplsm.f90 | 7 ++++++- 1 file changed, 6 insertions(+), 1 deletion(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 567f4a0cf..54bec6de5 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -2597,6 +2597,12 @@ subroutine albedo (parameters,vegtyp ,ist ,ice ,nsoil , & !in ftid(ib) = 0. ftii(ib) = 0. if (ib.eq.1) fsun = 0. + frevd = 0. + frevi = 0. + fregd = 0. + fregi = 0. + bgap = 0. + wgap = 0. end do if(cosz <= 0) goto 100 @@ -3262,7 +3268,6 @@ subroutine twostream (parameters,ib ,ic ,vegtyp ,cosz ,vai , & ! frev(ib) = freveg freg(ib) = frebar - ! flux absorbed by vegetation fab(ib) = 1. - fre(ib) - (1.-albgrd(ib))*ftd(ib) & From b67acadf85ca69e274e24097b53efd77b75add39 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 15 Mar 2021 05:44:13 -0600 Subject: [PATCH 24/31] Bugfixes in UGWP v1 for uninitialized variable azmeti(jk-1) --- physics/cires_ugwpv1_module.F90 | 25 ++++++++++++++----------- physics/cires_ugwpv1_solv2.F90 | 19 +++++++++---------- physics/ugwpv1_gsldrag.F90 | 5 +++-- 3 files changed, 26 insertions(+), 23 deletions(-) diff --git a/physics/cires_ugwpv1_module.F90 b/physics/cires_ugwpv1_module.F90 index 13b7752a5..c8108cba2 100644 --- a/physics/cires_ugwpv1_module.F90 +++ b/physics/cires_ugwpv1_module.F90 @@ -191,16 +191,19 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & integer :: k integer :: ddd_ugwp, curday_ugwp ! integer :: version - + + ! Initialize CCPP error handling variables + errmsg = '' + errflg = 0 ! if (me == master) print *, trim (fn_nml), ' GW-namelist file ' inquire (file =trim (fn_nml) , exist = exists) ! if (.not. exists) then - if (me == master) & - write (6, *) 'separate ugwp :: namelist file: ', trim (fn_nml), ' does not exist' - + write(errmsg,'(3a)') 'cires_ugwpv1_init: namelist file: ', trim (fn_nml), ' does not exist' + errflg = 1 + return else open (unit = nlunit, file = trim(fn_nml), action = 'read', status = 'old', iostat = ios) endif @@ -209,11 +212,6 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & close (nlunit) ! - ! Initialize CCPP error handling variables - errmsg = '' - errflg = 0 - - strsolver= knob_ugwp_orosolv curday_ugwp = jdat_gfs(1)*10000 + jdat_gfs(2)*100 +jdat_gfs(3) @@ -248,7 +246,7 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & ! allocate(fcor(latr), fcor2(latr) ) ! allocate( kvg(levs+1), ktg(levs+1) ) - allocate( krad(levs+1), kion(levs+1) ) + allocate( krad(levs+1), kion(levs+1) ) allocate( zkm(levs), pmb(levs) ) ! @@ -263,7 +261,12 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & ! ! find ilaunch ! - + if (knob_ugwp_palaunch gt 900.e2) then + write(errmsg,'(a,e16.7)') 'cires_ugwpv1_init: unrealistic value of knob_ugwp_palaunch', knob_ugwp_palaunch + errflg = 1 + return + endif + do k=levs, 1, -1 if (pmb(k) .gt. knob_ugwp_palaunch ) exit enddo diff --git a/physics/cires_ugwpv1_solv2.F90 b/physics/cires_ugwpv1_solv2.F90 index ee8f7bc83..afd94ff5c 100644 --- a/physics/cires_ugwpv1_solv2.F90 +++ b/physics/cires_ugwpv1_solv2.F90 @@ -267,15 +267,14 @@ subroutine cires_ugwpv1_ngw_solv2(mpi_id, master, im, levs, kdt, dtp, & !=====ksrc - aum(km2:levs) = um(jl,km2:levs) - avm(km2:levs) = vm(jl,km2:levs) - atm(km2:levs) = tm(jl,km2:levs) - aqm(km2:levs) = qm(jl,km2:levs) - azmet(km2:levs) = zmet(jl,km2:levs) - aprsi(km2:levs+1) = prsi(jl,km2:levs+1) - azmeti(km2:levs+1) = zmeti(jl,km2:levs+1) + aum(1:levs) = um(jl,1:levs) + avm(1:levs) = vm(jl,1:levs) + atm(1:levs) = tm(jl,1:levs) + aqm(1:levs) = qm(jl,1:levs) + azmet(1:levs) = zmet(jl,1:levs) + aprsi(1:levs+1) = prsi(jl,1:levs+1) + azmeti(1:levs+1) = zmeti(jl,1:levs+1) - rho_src = aprsl(ksrc)*rdi/atm(ksrc) taub_ch = max(tau_ngw(jl), tau_min) taub_src = taub_ch @@ -288,8 +287,8 @@ subroutine cires_ugwpv1_ngw_solv2(mpi_id, master, im, levs, kdt, dtp, & ! do jk = km2, levs dz_meti(jk) = azmeti(jk+1)-azmeti(jk) - dz_met(jk) = azmet(jk)-azmeti(jk-1) - enddo + dz_met(jk) = azmet(jk)-azmeti(jk-1) + enddo ! --------------------------------------------- ! interface mean flow parameters launch -> levs+1 ! --------------------------------------------- diff --git a/physics/ugwpv1_gsldrag.F90 b/physics/ugwpv1_gsldrag.F90 index 00fd42dbd..f473dc9bd 100644 --- a/physics/ugwpv1_gsldrag.F90 +++ b/physics/ugwpv1_gsldrag.F90 @@ -160,7 +160,7 @@ subroutine ugwpv1_gsldrag_init ( & if ( do_ugwp_v0_orog_only .or. do_ugwp_v0) then print *, ' ccpp do_ugwp_v0 active ', do_ugwp_v0 print *, ' ccpp do_ugwp_v1_orog_only active ', do_ugwp_v0_orog_only - write(errmsg,'(*(a))') " the CIRES CCPP-suite does not & + write(errmsg,'(*(a))') " the CIRES CCPP-suite does not & support schemes " errflg = 1 return @@ -171,7 +171,7 @@ subroutine ugwpv1_gsldrag_init ( & print *, ' do_ugwp_v1_w_gsldrag ', do_ugwp_v1_w_gsldrag print *, ' do_ugwp_v1_orog_only ', do_ugwp_v1_orog_only print *, ' do_gsl_drag_ls_bl ',do_gsl_drag_ls_bl - write(errmsg,'(*(a))') " the CIRES CCPP-suite intend to & + write(errmsg,'(*(a))') " the CIRES CCPP-suite intend to & support with but has Logic error" errflg = 1 return @@ -232,6 +232,7 @@ subroutine ugwpv1_gsldrag_init ( & call cires_ugwpv1_init (me, master, nlunit, logunit, jdat, con_pi, & con_rerth, fn_nml2, lonr, latr, levs, ak, bk, & con_p0, dtp, errmsg, errflg) + if (errflg/=0) return end if if (me == master) then From 10f532cd7feb8b8a6f9a6028d8e061b2a70b2b91 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 15 Mar 2021 15:19:12 -0600 Subject: [PATCH 25/31] Bugfix in physics/cires_ugwpv1_module.F90 --- physics/cires_ugwpv1_module.F90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/physics/cires_ugwpv1_module.F90 b/physics/cires_ugwpv1_module.F90 index c8108cba2..4746d61ff 100644 --- a/physics/cires_ugwpv1_module.F90 +++ b/physics/cires_ugwpv1_module.F90 @@ -261,7 +261,7 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & ! ! find ilaunch ! - if (knob_ugwp_palaunch gt 900.e2) then + if (knob_ugwp_palaunch > 900.e2) then write(errmsg,'(a,e16.7)') 'cires_ugwpv1_init: unrealistic value of knob_ugwp_palaunch', knob_ugwp_palaunch errflg = 1 return From 62df7ba52c8b559a7745a8e0382bf1503dd5337f Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 15 Mar 2021 15:24:48 -0600 Subject: [PATCH 26/31] Adjust two tuning parameters for RUC LSM in physics/module_sf_ruclsm.F90 and physics/module_soil_pre.F90 --- physics/module_sf_ruclsm.F90 | 3 ++- physics/module_soil_pre.F90 | 4 ++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/physics/module_sf_ruclsm.F90 b/physics/module_sf_ruclsm.F90 index a0e74ce7a..1eceaf183 100644 --- a/physics/module_sf_ruclsm.F90 +++ b/physics/module_sf_ruclsm.F90 @@ -2586,7 +2586,8 @@ SUBROUTINE SOIL (debug_print, & ! evaporation, effects sparsely vegetated areas--> cooler during the day ! fc=max(qmin,ref*0.25) ! ! For now we'll go back to ref*0.5 - fc=max(qmin,ref*0.5) +! Replace 0.5 with 0.7 2021/03/15 + fc=max(qmin,ref*0.7) fex_fc=1. if((soilmois(1)+qmin) > fc .or. (qvatm-qvg) > 0.) then soilres = 1. diff --git a/physics/module_soil_pre.F90 b/physics/module_soil_pre.F90 index 82fe23f24..8eb5a5775 100644 --- a/physics/module_soil_pre.F90 +++ b/physics/module_soil_pre.F90 @@ -42,8 +42,8 @@ SUBROUTINE init_soil_depth_3 ( zs , dzs , num_soil_levels ) IF ( num_soil_levels .EQ. 6) THEN zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 1.60 , 3.00 /) ELSEIF ( num_soil_levels .EQ. 9) THEN - !zs = (/ 0.00 , 0.01 , 0.04 , 0.10 , 0.30, 0.60, 1.00 , 1.60, 3.00 /) - zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 0.60, 1.00, 1.60 , 2.20, 3.00 /) + zs = (/ 0.00 , 0.01 , 0.04 , 0.10 , 0.30, 0.60, 1.00 , 1.60, 3.00 /) + !zs = (/ 0.00 , 0.05 , 0.20 , 0.40 , 0.60, 1.00, 1.60 , 2.20, 3.00 /) ENDIF zs2(1) = 0. From bbec192809f8e283b79eced5d032faf45060b423 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Mon, 15 Mar 2021 15:36:53 -0600 Subject: [PATCH 27/31] Optimize initializing intent(out) variables in physics/module_sf_noahmplsm.f90 --- physics/module_sf_noahmplsm.f90 | 10 ++++------ 1 file changed, 4 insertions(+), 6 deletions(-) diff --git a/physics/module_sf_noahmplsm.f90 b/physics/module_sf_noahmplsm.f90 index 54bec6de5..8041d07b4 100644 --- a/physics/module_sf_noahmplsm.f90 +++ b/physics/module_sf_noahmplsm.f90 @@ -2583,6 +2583,10 @@ subroutine albedo (parameters,vegtyp ,ist ,ice ,nsoil , & !in mpe = 1.e-06 bgap = 0. wgap = 0. + frevd = 0. + frevi = 0. + fregd = 0. + fregi = 0. ! initialize output because solar radiation only done if cosz > 0 @@ -2597,12 +2601,6 @@ subroutine albedo (parameters,vegtyp ,ist ,ice ,nsoil , & !in ftid(ib) = 0. ftii(ib) = 0. if (ib.eq.1) fsun = 0. - frevd = 0. - frevi = 0. - fregd = 0. - fregi = 0. - bgap = 0. - wgap = 0. end do if(cosz <= 0) goto 100 From e47c6acbc503e22178dcde84995dada2287d80a4 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 30 Mar 2021 08:08:21 -0600 Subject: [PATCH 28/31] Add logic to convert dry/moist mixing ratios based on flag to GFS_suite_interstitial_4_run --- physics/GFS_suite_interstitial.F90 | 85 ++++++++++++++++++----------- physics/GFS_suite_interstitial.meta | 8 +++ 2 files changed, 60 insertions(+), 33 deletions(-) diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index 4af19296a..a0361b6e3 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -155,11 +155,10 @@ end subroutine GFS_suite_interstitial_2_init subroutine GFS_suite_interstitial_2_finalize() end subroutine GFS_suite_interstitial_2_finalize -#if 0 + !> \section arg_table_GFS_suite_interstitial_2_run Argument Table !! \htmlinclude GFS_suite_interstitial_2_run.html !! -#endif subroutine GFS_suite_interstitial_2_run (im, levs, lssav, ldiag3d, lsidea, cplflx, flag_cice, shal_cnv, old_monin, mstrat, & do_shoc, frac_grid, imfshalcnv, dtf, xcosz, adjsfcdsw, adjsfcdlw, cice, pgr, ulwsfc_cice, lwhd, htrsw, htrlw, xmu, ctei_rm, & work1, work2, prsi, tgrs, prsl, qgrs_water_vapor, qgrs_cloud_water, cp, hvap, prslk, suntim, adjsfculw, adjsfculw_lnd, & @@ -475,11 +474,9 @@ end subroutine GFS_suite_interstitial_3_init subroutine GFS_suite_interstitial_3_finalize() end subroutine GFS_suite_interstitial_3_finalize -#if 0 !> \section arg_table_GFS_suite_interstitial_3_run Argument Table !! \htmlinclude GFS_suite_interstitial_3_run.html !! -#endif subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, & satmedmf, trans_trac, do_shoc, ltaerosol, ntrac, ntcw, & ntiw, ntclamt, ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, & @@ -513,7 +510,7 @@ subroutine GFS_suite_interstitial_3_run (im, levs, nn, cscnv, & real(kind=kind_phys), dimension(im, levs), intent(inout) :: rhc, save_qc ! save_qi is not allocated for Zhao-Carr MP real(kind=kind_phys), dimension(:, :), intent(inout) :: save_qi - real(kind=kind_phys), dimension(:, :), intent(inout) :: save_tcp ! ONLY ALLOCATE FOR THOMPSON! TODO + real(kind=kind_phys), dimension(:, :), intent(inout) :: save_tcp real(kind=kind_phys), dimension(im, levs, nn), intent(inout) :: clw character(len=*), intent(out) :: errmsg @@ -652,7 +649,7 @@ end subroutine GFS_suite_interstitial_4_finalize !! subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_total, ntrac, ntcw, ntiw, ntclamt, & ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & - imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, dtf, save_qc, save_qi, con_pi, & + imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, convert_dry_rho, dtf, save_qc, save_qi, con_pi, & gq0, clw, prsl, save_tcp, con_rd, nwfa, spechum, dqdti, errmsg, errflg) use machine, only: kind_phys @@ -666,7 +663,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & imp_physics_zhao_carr, imp_physics_zhao_carr_pdf - logical, intent(in) :: ltaerosol, cplchm + logical, intent(in) :: ltaerosol, cplchm, convert_dry_rho real(kind=kind_phys), intent(in) :: con_pi, dtf real(kind=kind_phys), dimension(im,levs), intent(in) :: save_qc @@ -739,33 +736,55 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to enddo if (imp_physics == imp_physics_thompson .and. (ntlnc>0 .or. ntinc>0)) then - do k=1,levs - do i=1,im - !> - Convert specific humidity to dry mixing ratio - qv_mp(i,k) = spechum(i,k) / (one-spechum(i,k)) - !> - Density of air in kg m-3 and inverse density - rho = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+0.622)) - orho = one/rho - if (ntlnc>0) then - !> - Convert moist mixing ratio to dry mixing ratio - qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) / (one-spechum(i,k)) - !> - Convert number concentration from moist to dry - nc_mp(i,k) = gq0(i,k,ntlnc) / (one-spechum(i,k)) - nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) - !> - Convert number concentrations from dry to moist - gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k)) - endif - if (ntinc>0) then - !> - Convert moist mixing ratio to dry mixing ratio - qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) / (one-spechum(i,k)) - !> - Convert number concentration from moist to dry - ni_mp(i,k) = gq0(i,k,ntinc) / (one-spechum(i,k)) - ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) - !> - Convert number concentrations from dry to moist - gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k)) - endif + if_convert_dry_rho: if (convert_dry_rho) then + do k=1,levs + do i=1,im + !> - Convert specific humidity to dry mixing ratio + qv_mp(i,k) = spechum(i,k) / (one-spechum(i,k)) + !> - Density of air in kg m-3 and inverse density + rho = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+0.622)) + orho = one/rho + if (ntlnc>0) then + !> - Convert moist mixing ratio to dry mixing ratio + qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) / (one-spechum(i,k)) + !> - Convert number concentration from moist to dry + nc_mp(i,k) = gq0(i,k,ntlnc) / (one-spechum(i,k)) + nc_mp(i,k) = max(zero, nc_mp(i,k) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) + !> - Convert number concentrations from dry to moist + gq0(i,k,ntlnc) = nc_mp(i,k) / (one+qv_mp(i,k)) + endif + if (ntinc>0) then + !> - Convert moist mixing ratio to dry mixing ratio + qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) / (one-spechum(i,k)) + !> - Convert number concentration from moist to dry + ni_mp(i,k) = gq0(i,k,ntinc) / (one-spechum(i,k)) + ni_mp(i,k) = max(zero, ni_mp(i,k) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) + !> - Convert number concentrations from dry to moist + gq0(i,k,ntinc) = ni_mp(i,k) / (one+qv_mp(i,k)) + endif + enddo enddo - enddo + else + do k=1,levs + do i=1,im + !> - Density of air in kg m-3 and inverse density + rho = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(spechum(i,k)+0.622)) + orho = one/rho + if (ntlnc>0) then + !> - Update cloud water mixing ratio + qc_mp(i,k) = (clw(i,k,2)-save_qc(i,k)) + !> - Update cloud water number concentration + gq0(i,k,ntlnc) = max(zero, gq0(i,k,ntlnc) + make_DropletNumber(qc_mp(i,k) * rho, nwfa(i,k)*rho) * orho) + endif + if (ntinc>0) then + !> - Update cloud ice mixing ratio + qi_mp(i,k) = (clw(i,k,1)-save_qi(i,k)) + !> - Update cloud ice number concentration + gq0(i,k,ntinc) = max(zero, gq0(i,k,ntinc) + make_IceNumber(qi_mp(i,k) * rho, save_tcp(i,k)) * orho) + endif + enddo + enddo + end if if_convert_dry_rho endif else diff --git a/physics/GFS_suite_interstitial.meta b/physics/GFS_suite_interstitial.meta index c09d02434..b7c1c2f67 100644 --- a/physics/GFS_suite_interstitial.meta +++ b/physics/GFS_suite_interstitial.meta @@ -1769,6 +1769,14 @@ type = integer intent = in optional = F +[convert_dry_rho] + standard_name = flag_for_converting_hydrometeors_from_moist_to_dry_air + long_name = flag for converting hydrometeors from moist to dry air + units = flag + dimensions = () + type = logical + intent = in + optional = F [dtf] standard_name = time_step_for_dynamics long_name = dynamics timestep From 92e0702a1010345e59cc5102f596af250d5910fe Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 30 Mar 2021 08:10:51 -0600 Subject: [PATCH 29/31] Revert CODEOWNERS change --- CODEOWNERS | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/CODEOWNERS b/CODEOWNERS index b6c597371..0d5230f89 100644 --- a/CODEOWNERS +++ b/CODEOWNERS @@ -3,7 +3,7 @@ # These owners will be the default owners for everything in the repo. #* @defunkt -* @DomHeinzeller +* @climbfuji @llpcarson @grantfirl @JulieSchramm # Order is important. The last matching pattern has the most precedence. # So if a pull request only touches javascript files, only these owners From 252d0d38d75330a6e267dd4cd6d08a28c3abb917 Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 30 Mar 2021 14:08:44 -0600 Subject: [PATCH 30/31] Bugfix in cires_ugwpv1_module.F90: add missing INTERNAL_FILE_NML logic --- physics/cires_ugwpv1_module.F90 | 14 ++++++++------ physics/ugwpv1_gsldrag.F90 | 4 ++-- 2 files changed, 10 insertions(+), 8 deletions(-) diff --git a/physics/cires_ugwpv1_module.F90 b/physics/cires_ugwpv1_module.F90 index c1fd53523..f5fe7f2ec 100644 --- a/physics/cires_ugwpv1_module.F90 +++ b/physics/cires_ugwpv1_module.F90 @@ -143,8 +143,8 @@ module cires_ugwpv1_module !----------------------------------------------------------------------------------- subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & - con_rerth, fn_nml2, lonr, latr, levs, ak, bk, pref, dtp, & - errmsg, errflg) + con_rerth, fn_nml2, input_nml_file, lonr, latr, levs, ak, bk, & + pref, dtp, errmsg, errflg) ! ! input_nml_file ='input.nml'=fn_nml ..... OLD_namelist and cdmvgwd(4) Corrected Bug Oct 4 ! @@ -177,7 +177,7 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & real(kind=kind_phys), intent (in) :: con_pi, con_rerth character(len=64), intent (in) :: fn_nml2 -! character(len=64), parameter :: fn_nml='input.nml' + character(len=*), intent (in) :: input_nml_file(:) character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -187,7 +187,7 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & integer :: ios logical :: exists - integer :: ncid, iernc, vid, dimid, status + integer :: ncid, iernc, vid, dimid, status integer :: k integer :: ddd_ugwp, curday_ugwp ! integer :: version @@ -196,7 +196,9 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & errmsg = '' errflg = 0 -! +#ifdef INTERNAL_FILE_NML + read (input_nml_file, nml = cires_ugwp_nml) +#else if (me == master) print *, trim (fn_nml2), ' GW-namelist file ' inquire (file =trim (fn_nml2) , exist = exists) ! @@ -210,7 +212,7 @@ subroutine cires_ugwpv1_init (me, master, nlunit, logunit, jdat_gfs, con_pi, & rewind (nlunit) read (nlunit, nml = cires_ugwp_nml) close (nlunit) -! +#endif strsolver= knob_ugwp_orosolv diff --git a/physics/ugwpv1_gsldrag.F90 b/physics/ugwpv1_gsldrag.F90 index f473dc9bd..649d994a8 100644 --- a/physics/ugwpv1_gsldrag.F90 +++ b/physics/ugwpv1_gsldrag.F90 @@ -230,8 +230,8 @@ subroutine ugwpv1_gsldrag_init ( & if ( do_ugwp_v1 ) then call cires_ugwpv1_init (me, master, nlunit, logunit, jdat, con_pi, & - con_rerth, fn_nml2, lonr, latr, levs, ak, bk, & - con_p0, dtp, errmsg, errflg) + con_rerth, fn_nml2, input_nml_file, lonr, latr, & + levs, ak, bk, con_p0, dtp, errmsg, errflg) if (errflg/=0) return end if From 5dbb259348ecd70f2e10f6aab57d42c303c11c1c Mon Sep 17 00:00:00 2001 From: Dom Heinzeller Date: Tue, 30 Mar 2021 15:39:30 -0600 Subject: [PATCH 31/31] Cleanup: use con_eps = con_rd/con_rv instead of 0.622 in many places --- physics/GFS_phys_time_vary.fv3.F90 | 3 -- physics/GFS_phys_time_vary.scm.F90 | 3 -- physics/GFS_rrtmg_pre.F90 | 10 ++--- physics/GFS_rrtmg_pre.meta | 2 +- physics/GFS_rrtmgp_thompsonmp_pre.F90 | 7 ++-- physics/GFS_rrtmgp_thompsonmp_pre.meta | 55 +++++++++++++++----------- physics/GFS_suite_interstitial.F90 | 8 ++-- physics/GFS_suite_interstitial.meta | 9 +++++ physics/gcm_shoc.F90 | 14 ++++--- physics/gcm_shoc.meta | 9 +++++ physics/gfdl_cloud_microphys.F90 | 6 +-- physics/gfdl_cloud_microphys.meta | 9 +++++ physics/gfdl_sfc_layer.F90 | 14 ++++--- physics/gfdl_sfc_layer.meta | 6 +++ physics/mp_thompson.F90 | 15 +++---- physics/mp_thompson.meta | 18 +++++++++ 16 files changed, 124 insertions(+), 64 deletions(-) diff --git a/physics/GFS_phys_time_vary.fv3.F90 b/physics/GFS_phys_time_vary.fv3.F90 index 57d253083..0cc6a66b8 100644 --- a/physics/GFS_phys_time_vary.fv3.F90 +++ b/physics/GFS_phys_time_vary.fv3.F90 @@ -446,9 +446,6 @@ subroutine GFS_phys_time_vary_init ( eahxy(ix) = 2000.0_kind_phys -! eahxy = psfc*qv/(0.622+qv); qv is mixing ratio, converted from sepcific -! humidity specific humidity /(1.0 - specific humidity) - cmxy(ix) = zero chxy(ix) = zero fwetxy(ix) = zero diff --git a/physics/GFS_phys_time_vary.scm.F90 b/physics/GFS_phys_time_vary.scm.F90 index a54ffa7a9..fb46de2bd 100644 --- a/physics/GFS_phys_time_vary.scm.F90 +++ b/physics/GFS_phys_time_vary.scm.F90 @@ -387,9 +387,6 @@ subroutine GFS_phys_time_vary_init ( eahxy(ix) = 2000.0_kind_phys -! eahxy = psfc*qv/(0.622+qv); qv is mixing ratio, converted from sepcific -! humidity specific humidity /(1.0 - specific humidity) - cmxy(ix) = zero chxy(ix) = zero fwetxy(ix) = zero diff --git a/physics/GFS_rrtmg_pre.F90 b/physics/GFS_rrtmg_pre.F90 index a622cf8f0..0c9eaf3f0 100644 --- a/physics/GFS_rrtmg_pre.F90 +++ b/physics/GFS_rrtmg_pre.F90 @@ -24,7 +24,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & imp_physics_zhao_carr_pdf, imp_physics_mg, imp_physics_wsm6, & imp_physics_fer_hires, julian, yearlen, lndp_var_list, lsswr, lslwr, & ltaerosol, lgfdlmprad, uni_cld, effr_in, do_mynnedmf, lmfshal, & - lmfdeep2, fhswr, fhlwr, solhr, sup, eps, epsm1, fvirt, & + lmfdeep2, fhswr, fhlwr, solhr, sup, con_eps, epsm1, fvirt, & rog, rocp, con_rd, xlat_d, xlat, xlon, coslat, sinlat, tsfc, slmsk, & prsi, prsl, prslk, tgrs, sfc_wts, mg_cld, effrr_in, pert_clds,sppt_wts,& sppt_amp, cnvw_in, cnvc_in, qgrs, aer_nm, dx, icloud, & !inputs from here and above @@ -103,7 +103,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & lmfshal, lmfdeep2, pert_clds real(kind=kind_phys), intent(in) :: fhswr, fhlwr, solhr, sup, julian, sppt_amp - real(kind=kind_phys), intent(in) :: eps, epsm1, fvirt, rog, rocp, con_rd + real(kind=kind_phys), intent(in) :: con_eps, epsm1, fvirt, rog, rocp, con_rd real(kind=kind_phys), dimension(:), intent(in) :: xlat_d, xlat, xlon, & coslat, sinlat, tsfc, & @@ -300,7 +300,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & !> - Compute relative humidity. es = min( prsl(i,k2), fpvs( tgrs(i,k2) ) ) ! fpvs and prsl in pa - qs = max( QMIN, eps * es / (prsl(i,k2) + epsm1*es) ) + qs = max( QMIN, con_eps * es / (prsl(i,k2) + epsm1*es) ) rhly(i,k1) = max( 0.0, min( 1.0, max(QMIN, qgrs(i,k2,ntqv))/qs ) ) qstl(i,k1) = qs enddo @@ -643,7 +643,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & do i=1,IM qvs = qgrs(i,k,ntqv) qv_mp (i,k) = qvs/(1.-qvs) - rho (i,k) = 0.622*prsl(i,k)/(con_rd*tgrs(i,k)*(qv_mp(i,k)+0.622)) + rho (i,k) = con_eps*prsl(i,k)/(con_rd*tgrs(i,k)*(qv_mp(i,k)+con_eps)) orho (i,k) = 1.0/rho(i,k) qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs) qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs) @@ -658,7 +658,7 @@ subroutine GFS_rrtmg_pre_run (im, levs, lm, lmk, lmp, n_var_lndp, & do i=1,IM qvs = qgrs(i,k,ntqv) qv_mp (i,k) = qvs/(1.-qvs) - rho (i,k) = 0.622*prsl(i,k)/(con_rd*tgrs(i,k)*(qv_mp(i,k)+0.622)) + rho (i,k) = con_eps*prsl(i,k)/(con_rd*tgrs(i,k)*(qv_mp(i,k)+con_eps)) orho (i,k) = 1.0/rho(i,k) qc_mp (i,k) = tracer1(i,k,ntcw)/(1.-qvs) qi_mp (i,k) = tracer1(i,k,ntiw)/(1.-qvs) diff --git a/physics/GFS_rrtmg_pre.meta b/physics/GFS_rrtmg_pre.meta index 14403f63d..e26cdeac1 100644 --- a/physics/GFS_rrtmg_pre.meta +++ b/physics/GFS_rrtmg_pre.meta @@ -481,7 +481,7 @@ kind = kind_phys intent = in optional = F -[eps] +[con_eps] standard_name = ratio_of_dry_air_to_water_vapor_gas_constants long_name = rd/rv units = none diff --git a/physics/GFS_rrtmgp_thompsonmp_pre.F90 b/physics/GFS_rrtmgp_thompsonmp_pre.F90 index c6661d948..b54f27d65 100644 --- a/physics/GFS_rrtmgp_thompsonmp_pre.F90 +++ b/physics/GFS_rrtmgp_thompsonmp_pre.F90 @@ -41,7 +41,7 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do i_cldliq, i_cldice, i_cldrain, i_cldsnow, i_cldgrpl, i_cldtot, i_cldliq_nc, & i_cldice_nc, i_twa, effr_in, p_lev, p_lay, tv_lay, t_lay, effrin_cldliq, & effrin_cldice, effrin_cldsnow, tracer, qs_lay, q_lay, relhum, cld_frac_mg, con_g, & - con_rd, uni_cld, lmfshal, lmfdeep2, ltaerosol, do_mynnedmf, imfdeepcnv, & + con_rd, con_eps, uni_cld, lmfshal, lmfdeep2, ltaerosol, do_mynnedmf, imfdeepcnv, & imfdeepcnv_gf, doGP_cldoptics_PADE, doGP_cldoptics_LUT, & cld_frac, cld_lwp, cld_reliq, cld_iwp, cld_reice, cld_swp, cld_resnow, cld_rwp, & cld_rerain, precip_frac, errmsg, errflg) @@ -76,7 +76,8 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do doGP_cldoptics_PADE ! (PADE approximation) real(kind_phys), intent(in) :: & con_g, & ! Physical constant: gravitational constant - con_rd ! Physical constant: gas-constant for dry air + con_rd, & ! Physical constant: gas-constant for dry air + con_eps ! Physical constant: gas constant air / gas constant H2O real(kind_phys), dimension(nCol,nLev), intent(in) :: & tv_lay, & ! Virtual temperature (K) @@ -155,7 +156,7 @@ subroutine GFS_rrtmgp_thompsonmp_pre_run(nCol, nLev, nTracers, ncnd, doSWrad, do do iLay = 1, nLev do iCol = 1, nCol qv_mp(iCol,iLay) = q_lay(iCol,iLay)/(1.-q_lay(iCol,iLay)) - rho = 0.622*p_lay(iCol,iLay)/(con_rd*t_lay(iCol,iLay)*(qv_mp(iCol,iLay)+0.622)) + rho = con_eps*p_lay(iCol,iLay)/(con_rd*t_lay(iCol,iLay)*(qv_mp(iCol,iLay)+con_eps)) orho = 1./rho qc_mp(iCol,iLay) = tracer(iCol,iLay,i_cldliq) / (1.-q_lay(iCol,iLay)) qi_mp(iCol,iLay) = tracer(iCol,iLay,i_cldice) / (1.-q_lay(iCol,iLay)) diff --git a/physics/GFS_rrtmgp_thompsonmp_pre.meta b/physics/GFS_rrtmgp_thompsonmp_pre.meta index 90ec59760..54d266b67 100644 --- a/physics/GFS_rrtmgp_thompsonmp_pre.meta +++ b/physics/GFS_rrtmgp_thompsonmp_pre.meta @@ -171,7 +171,7 @@ standard_name = effective_radius_of_stratiform_cloud_liquid_water_particle_in_um long_name = eff. radius of cloud liquid water particle in micrometer units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -180,7 +180,7 @@ standard_name = effective_radius_of_stratiform_cloud_ice_particle_in_um long_name = eff. radius of cloud ice water particle in micrometer units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -189,7 +189,7 @@ standard_name = effective_radius_of_stratiform_cloud_snow_particle_in_um long_name = effective radius of cloud snow particle in micrometers units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -247,7 +247,7 @@ standard_name = air_pressure_at_interface_for_RRTMGP_in_hPa long_name = air pressure at vertical interface for radiation calculation units = hPa - dimensions = (horizontal_dimension,vertical_dimension_plus_one) + dimensions = (horizontal_loop_extent,vertical_dimension_plus_one) type = real kind = kind_phys intent = in @@ -256,7 +256,7 @@ standard_name = air_pressure_at_layer_for_RRTMGP_in_hPa long_name = air pressure at vertical layer for radiation calculation units = hPa - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = in @@ -265,7 +265,7 @@ standard_name = virtual_temperature long_name = layer virtual temperature units = K - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = in @@ -283,7 +283,7 @@ standard_name = saturation_vapor_pressure long_name = saturation vapor pressure units = Pa - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = in @@ -292,7 +292,7 @@ standard_name = water_vapor_mixing_ratio long_name = water vaport mixing ratio units = kg/kg - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = in @@ -301,7 +301,7 @@ standard_name = relative_humidity long_name = layer relative humidity units = frac - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = in @@ -310,7 +310,7 @@ standard_name = chemical_tracers long_name = chemical tracers units = g g-1 - dimensions = (horizontal_dimension,vertical_dimension,number_of_tracers) + dimensions = (horizontal_loop_extent,vertical_dimension,number_of_tracers) type = real kind = kind_phys intent = in @@ -330,14 +330,23 @@ units = J kg-1 K-1 dimensions = () type = real - kind = kind_phys + kind = kind_phys + intent = in + optional = F +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys intent = in - optional = F + optional = F [cld_frac] standard_name = total_cloud_fraction long_name = layer total cloud fraction units = frac - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -346,7 +355,7 @@ standard_name = cloud_liquid_water_path long_name = layer cloud liquid water path units = g m-2 - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -355,7 +364,7 @@ standard_name = mean_effective_radius_for_liquid_cloud long_name = mean effective radius for liquid cloud units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -364,7 +373,7 @@ standard_name = cloud_ice_water_path long_name = layer cloud ice water path units = g m-2 - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -373,7 +382,7 @@ standard_name = mean_effective_radius_for_ice_cloud long_name = mean effective radius for ice cloud units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -382,7 +391,7 @@ standard_name = cloud_snow_water_path long_name = layer cloud snow water path units = g m-2 - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -391,7 +400,7 @@ standard_name = mean_effective_radius_for_snow_flake long_name = mean effective radius for snow cloud units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -400,7 +409,7 @@ standard_name = cloud_rain_water_path long_name = layer cloud rain water path units = g m-2 - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -409,7 +418,7 @@ standard_name = mean_effective_radius_for_rain_drop long_name = mean effective radius for rain cloud units = um - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout @@ -418,11 +427,11 @@ standard_name = precipitation_fraction_by_layer long_name = precipitation fraction in each layer units = frac - dimensions = (horizontal_dimension,vertical_dimension) + dimensions = (horizontal_loop_extent,vertical_dimension) type = real kind = kind_phys intent = inout - optional = F + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/GFS_suite_interstitial.F90 b/physics/GFS_suite_interstitial.F90 index a0361b6e3..93106d2de 100644 --- a/physics/GFS_suite_interstitial.F90 +++ b/physics/GFS_suite_interstitial.F90 @@ -650,7 +650,7 @@ end subroutine GFS_suite_interstitial_4_finalize subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_total, ntrac, ntcw, ntiw, ntclamt, & ntrw, ntsw, ntrnc, ntsnc, ntgl, ntgnc, ntlnc, ntinc, nn, imp_physics, imp_physics_gfdl, imp_physics_thompson, & imp_physics_zhao_carr, imp_physics_zhao_carr_pdf, convert_dry_rho, dtf, save_qc, save_qi, con_pi, & - gq0, clw, prsl, save_tcp, con_rd, nwfa, spechum, dqdti, errmsg, errflg) + gq0, clw, prsl, save_tcp, con_rd, con_eps, nwfa, spechum, dqdti, errmsg, errflg) use machine, only: kind_phys use module_mp_thompson_make_number_concentrations, only: make_IceNumber, make_DropletNumber @@ -673,7 +673,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to real(kind=kind_phys), dimension(im,levs,ntrac), intent(inout) :: gq0 real(kind=kind_phys), dimension(im,levs,nn), intent(inout) :: clw real(kind=kind_phys), dimension(im,levs), intent(in) :: prsl - real(kind=kind_phys), intent(in) :: con_rd + real(kind=kind_phys), intent(in) :: con_rd, con_eps real(kind=kind_phys), dimension(:,:), intent(in) :: nwfa, save_tcp real(kind=kind_phys), dimension(im,levs), intent(in) :: spechum @@ -742,7 +742,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to !> - Convert specific humidity to dry mixing ratio qv_mp(i,k) = spechum(i,k) / (one-spechum(i,k)) !> - Density of air in kg m-3 and inverse density - rho = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+0.622)) + rho = con_eps*prsl(i,k) / (con_rd*save_tcp(i,k)*(qv_mp(i,k)+con_eps)) orho = one/rho if (ntlnc>0) then !> - Convert moist mixing ratio to dry mixing ratio @@ -768,7 +768,7 @@ subroutine GFS_suite_interstitial_4_run (im, levs, ltaerosol, cplchm, tracers_to do k=1,levs do i=1,im !> - Density of air in kg m-3 and inverse density - rho = 0.622*prsl(i,k) / (con_rd*save_tcp(i,k)*(spechum(i,k)+0.622)) + rho = con_eps*prsl(i,k) / (con_rd*save_tcp(i,k)*(spechum(i,k)+con_eps)) orho = one/rho if (ntlnc>0) then !> - Update cloud water mixing ratio diff --git a/physics/GFS_suite_interstitial.meta b/physics/GFS_suite_interstitial.meta index b7c1c2f67..5b4b0dbf9 100644 --- a/physics/GFS_suite_interstitial.meta +++ b/physics/GFS_suite_interstitial.meta @@ -1858,6 +1858,15 @@ kind = kind_phys intent = in optional = F +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [nwfa] standard_name = water_friendly_aerosol_number_concentration long_name = number concentration of water-friendly aerosols diff --git a/physics/gcm_shoc.F90 b/physics/gcm_shoc.F90 index dd7791e18..af7d6db49 100644 --- a/physics/gcm_shoc.F90 +++ b/physics/gcm_shoc.F90 @@ -24,7 +24,7 @@ end subroutine shoc_finalize !! \htmlinclude shoc_run.html !! subroutine shoc_run (nx, nzm, tcr, tcrf, con_cp, con_g, con_hvap, con_hfus, con_rv, con_rd, & - con_pi, con_fvirt, dtp, prsl, delp, phii, phil, u, v, omega, rhc, & + con_pi, con_fvirt, con_eps, dtp, prsl, delp, phii, phil, u, v, omega, rhc, & supice, pcrit, cefac, cesfac, tkef1, dis_opt, hflx, evap, prnum, & gt0, gq0, ntrac, ntqv, ntcw, ntiw, ntrw, ntsw, ntgl, ntlnc, ntinc, & cld_sgs, tke, tkh, wthv_sec, errmsg, errflg) @@ -32,7 +32,8 @@ subroutine shoc_run (nx, nzm, tcr, tcrf, con_cp, con_g, con_hvap, con_hfus, con_ implicit none integer, intent(in) :: nx, nzm, ntrac, ntqv, ntcw, ntiw, ntrw, ntsw, ntgl, ntlnc, ntinc - real(kind=kind_phys), intent(in) :: tcr, tcrf, con_cp, con_g, con_hvap, con_hfus, con_rv, con_rd, con_pi, con_fvirt, & + real(kind=kind_phys), intent(in) :: tcr, tcrf, con_cp, con_g, con_hvap, con_hfus, con_rv, & + con_rd, con_pi, con_fvirt, con_eps, & dtp, supice, pcrit, cefac, cesfac, tkef1, dis_opt ! real(kind=kind_phys), intent(in), dimension(nx) :: hflx, evap @@ -118,7 +119,8 @@ subroutine shoc_run (nx, nzm, tcr, tcrf, con_cp, con_g, con_hvap, con_hfus, con_ rhc, supice, pcrit, cefac, cesfac, tkef1, dis_opt, & cld_sgs, tke, hflx, evap, prnum, tkh, wthv_sec, & ntlnc, ncpl, ncpi, & - con_cp, con_g, con_hvap, con_hfus, con_rv, con_rd, con_pi, con_fvirt) + con_cp, con_g, con_hvap, con_hfus, con_rv, con_rd, con_pi, & + con_fvirt, con_eps) if (ntiw < 0) then ! this is valid only for Zhao-Carr scheme do k=1,nzm @@ -166,13 +168,13 @@ subroutine shoc_work (ix, nx, nzm, nz, dtn, & pcrit, cefac, cesfac, tkef1, dis_opt, & cld_sgs, tke, hflx, evap, prnum, tkh, & wthv_sec, ntlnc, ncpl, ncpi, & - cp, ggr, lcond, lfus, rv, rgas, pi, epsv) + cp, ggr, lcond, lfus, rv, rgas, pi, epsv, eps) use funcphys , only : fpvsl, fpvsi, fpvs ! saturation vapor pressure for water & ice implicit none - real, intent(in) :: cp, ggr, lcond, lfus, rv, rgas, pi, epsv + real, intent(in) :: cp, ggr, lcond, lfus, rv, rgas, pi, epsv, eps integer, intent(in) :: ix ! max number of points in the physics window in the x integer, intent(in) :: nx ! Number of points in the physics window in the x @@ -219,7 +221,7 @@ subroutine shoc_work (ix, nx, nzm, nz, dtn, & real, intent(in) :: prnum (nx,nzm) ! turbulent Prandtl number real, intent(inout) :: wthv_sec (ix,nzm) ! Buoyancy flux, K*m/s - real, parameter :: zero=0.0_kp, one=1.0_kp, half=0.5_kp, two=2.0_kp, eps=0.622_kp, & + real, parameter :: zero=0.0_kp, one=1.0_kp, half=0.5_kp, two=2.0_kp, & three=3.0_kp, oneb3=one/three, twoby3=two/three, fourb3=twoby3+twoby3 real, parameter :: sqrt2 = sqrt(two), twoby15 = two / 15.0_kp, & nmin = 1.0_kp, RI_cub = 6.4e-14_kp, RL_cub = 1.0e-15_kp, & diff --git a/physics/gcm_shoc.meta b/physics/gcm_shoc.meta index 8cb03727d..047286317 100644 --- a/physics/gcm_shoc.meta +++ b/physics/gcm_shoc.meta @@ -113,6 +113,15 @@ kind = kind_phys intent = in optional = F +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [dtp] standard_name = time_step_for_physics long_name = time step for physics diff --git a/physics/gfdl_cloud_microphys.F90 b/physics/gfdl_cloud_microphys.F90 index 1ccedb956..c41323ad5 100644 --- a/physics/gfdl_cloud_microphys.F90 +++ b/physics/gfdl_cloud_microphys.F90 @@ -113,7 +113,7 @@ end subroutine gfdl_cloud_microphys_finalize !! \htmlinclude gfdl_cloud_microphys_run.html !! subroutine gfdl_cloud_microphys_run( & - levs, im, con_g, con_fvirt, con_rd, frland, garea, islmsk, & + levs, im, con_g, con_fvirt, con_rd, con_eps, frland, garea, islmsk, & gq0, gq0_ntcw, gq0_ntrw, gq0_ntiw, gq0_ntsw, gq0_ntgl, gq0_ntclamt, & gt0, gu0, gv0, vvl, prsl, phii, del, & rain0, ice0, snow0, graupel0, prcp0, sr, & @@ -134,7 +134,7 @@ subroutine gfdl_cloud_microphys_run( & ! interface variables integer, intent(in ) :: levs, im - real(kind=kind_phys), intent(in ) :: con_g, con_fvirt, con_rd + real(kind=kind_phys), intent(in ) :: con_g, con_fvirt, con_rd, con_eps real(kind=kind_phys), intent(in ), dimension(1:im) :: frland, garea integer, intent(in ), dimension(1:im) :: islmsk real(kind=kind_phys), intent(inout), dimension(1:im,1:levs) :: gq0, gq0_ntcw, gq0_ntrw, gq0_ntiw, & @@ -295,7 +295,7 @@ subroutine gfdl_cloud_microphys_run( & allocate(den(1:im,1:levs)) do k=1,levs do i=1,im - den(i,k)=0.622*prsl(i,k)/(con_rd*gt0(i,k)*(gq0(i,k)+0.622)) + den(i,k)=con_eps*prsl(i,k)/(con_rd*gt0(i,k)*(gq0(i,k)+con_eps)) enddo enddo call cloud_diagnosis (1, im, 1, levs, den(1:im,1:levs), & diff --git a/physics/gfdl_cloud_microphys.meta b/physics/gfdl_cloud_microphys.meta index 07847ed17..961a3e33f 100644 --- a/physics/gfdl_cloud_microphys.meta +++ b/physics/gfdl_cloud_microphys.meta @@ -168,6 +168,15 @@ kind = kind_phys intent = in optional = F +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [frland] standard_name = land_area_fraction_for_microphysics long_name = land area fraction used in microphysics schemes diff --git a/physics/gfdl_sfc_layer.F90 b/physics/gfdl_sfc_layer.F90 index 93e38c982..008e716e2 100644 --- a/physics/gfdl_sfc_layer.F90 +++ b/physics/gfdl_sfc_layer.F90 @@ -377,7 +377,7 @@ subroutine gfdl_sfc_layer_run (im, nsoil, km, xlat, xlon, flag_iter, lsm, & xxfh(i), ztmax(i), z0max(i), tstrc(i), & pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), & scurx(i), scury(i), pert_Cd, ens_random_seed, ens_Cdamp, upc(i), vpc(i), t1(i), q1(i), & - dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), errmsg, & + dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, & errflg) if (errflg /= 0) return @@ -526,7 +526,7 @@ subroutine gfdl_sfc_layer_run (im, nsoil, km, xlat, xlon, flag_iter, lsm, & xxfh(i), ztmax(i), z0max(i), tstrc(i), & pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), & scurx(i), scury(i), pert_Cd, ens_random_seed, ens_Cdamp, upc(i), vpc(i), t1(i), q1(i), & - dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), errmsg, & + dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, & errflg) if (errflg /= 0) return @@ -633,7 +633,7 @@ subroutine gfdl_sfc_layer_run (im, nsoil, km, xlat, xlon, flag_iter, lsm, & xxfh(i), znt_ocn(i), mznt(i), tstrc(i), & pspc(i), pkmax(i), wetc(i), slwdc(i), z1_cm(i), icoef_sf, iwavecpl, lcurr_sf, charn(i), msang(i), & scurx(i), scury(i), pert_Cd, ens_random_seed, ens_Cdamp, upc(i), vpc(i), t1(i), q1(i), & - dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), errmsg, & + dt, wind10(i), xxfh2(i), ntsflg, sfenth, tzot(i), ep2, errmsg, & errflg) if (errflg /= 0) return @@ -756,7 +756,7 @@ SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !m icoef_sf,iwavecpl,lcurr_sf,alpha,gamma,xcur,ycur, & pert_Cd, ens_random_seed, ens_Cdamp, & upc,vpc,tpc,rpc,dt,wind10,xxfh2,ntsflg,sfenth, & - tzot, errmsg, errflg) + tzot, ep2, errmsg, errflg) !------------------------------------------------------------------------ ! @@ -819,6 +819,8 @@ SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !m real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: tpc real(kind=kind_phys), intent ( in), dimension (ims :ime ) :: rpc + real(kind=kind_phys), intent ( in) :: ep2 + character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -1207,7 +1209,7 @@ SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !m if(psps1 .EQ. 0.0)then psps1 = .1 endif - rstso(i) = 0.622*estso(i)/psps1 + rstso(i) = ep2*estso(i)/psps1 vrts (i) = 1. + boycon*ecof(i)*rstso(i) enddo @@ -1735,7 +1737,7 @@ SUBROUTINE MFLUX2( fxh,fxe,fxmx,fxmy,cdm,rib,xxfh,zoc,mzoc,tstrc, & !m if(psps2 .EQ. 0.0)then psps2 = .1 endif - rstsop(i) = 0.622*estsop(i)/psps2 + rstsop(i) = ep2*estsop(i)/psps2 rdiff (i) = amin1(0.0,(rkmaxp(i) - rstsop(i))) foft(i) = tss(i) + delsrad(i)*(slwa(i) - aap(i)*tsp(i)**4 - & diff --git a/physics/gfdl_sfc_layer.meta b/physics/gfdl_sfc_layer.meta index 77024c813..a9829fec3 100644 --- a/physics/gfdl_sfc_layer.meta +++ b/physics/gfdl_sfc_layer.meta @@ -783,6 +783,8 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys + intent = inout + optional = F [qss_lnd] standard_name = surface_specific_humidity_over_land long_name = surface air saturation specific humidity over land @@ -790,6 +792,8 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys + intent = inout + optional = F [qss_ice] standard_name = surface_specific_humidity_over_ice long_name = surface air saturation specific humidity over ice @@ -797,6 +801,8 @@ dimensions = (horizontal_loop_extent) type = real kind = kind_phys + intent = inout + optional = F [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 5d5e631f5..6031af94a 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -28,9 +28,9 @@ module mp_thompson !! \section arg_table_mp_thompson_init Argument Table !! \htmlinclude mp_thompson_init.html !! - subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & - imp_physics, imp_physics_thompson, & - convert_dry_rho, & + subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, con_eps, & + restart, imp_physics, & + imp_physics_thompson, convert_dry_rho,& spechum, qc, qr, qi, qs, qg, ni, nr, & is_aerosol_aware, nc, nwfa2d, nifa2d, & nwfa, nifa, tgrs, prsl, phil, area, & @@ -43,7 +43,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & ! Interface variables integer, intent(in ) :: ncol integer, intent(in ) :: nlev - real(kind_phys), intent(in ) :: con_g, con_rd + real(kind_phys), intent(in ) :: con_g, con_rd, con_eps logical, intent(in ) :: restart integer, intent(in ) :: imp_physics integer, intent(in ) :: imp_physics_thompson @@ -160,7 +160,7 @@ subroutine mp_thompson_init(ncol, nlev, con_g, con_rd, restart, & end if ! Density of moist air in kg m-3 and inverse density of air - rho = 0.622*prsl/(con_rd*tgrs*(qv+0.622)) + rho = con_eps*prsl/(con_rd*tgrs*(qv+con_eps)) orho = 1.0/rho ! Ensure we have 1st guess ice number where mass non-zero but no number. @@ -324,7 +324,7 @@ end subroutine mp_thompson_init !>\section gen_thompson_hrrr Thompson MP General Algorithm !>@{ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & - convert_dry_rho, & + con_eps, convert_dry_rho, & spechum, qc, qr, qi, qs, qg, ni, nr, & is_aerosol_aware, nc, nwfa, nifa, & nwfa2d, nifa2d, & @@ -344,6 +344,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & integer, intent(in ) :: nlev real(kind_phys), intent(in ) :: con_g real(kind_phys), intent(in ) :: con_rd + real(kind_phys), intent(in ) :: con_eps ! Hydrometeors logical, intent(in ) :: convert_dry_rho real(kind_phys), intent(inout) :: spechum(1:ncol,1:nlev) @@ -474,7 +475,7 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & end if !> - Density of air in kg m-3 - rho = 0.622*prsl/(con_rd*tgrs*(qv+0.622)) + rho = con_eps*prsl/(con_rd*tgrs*(qv+con_eps)) !> - Convert omega in Pa s-1 to vertical velocity w in m s-1 w = -omega/(rho*con_g) diff --git a/physics/mp_thompson.meta b/physics/mp_thompson.meta index ed54f8d02..7f1e9197e 100644 --- a/physics/mp_thompson.meta +++ b/physics/mp_thompson.meta @@ -41,6 +41,15 @@ kind = kind_phys intent = in optional = F +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [restart] standard_name = flag_for_restart long_name = flag for restart (warmstart) or coldstart @@ -349,6 +358,15 @@ kind = kind_phys intent = in optional = F +[con_eps] + standard_name = ratio_of_dry_air_to_water_vapor_gas_constants + long_name = rd/rv + units = none + dimensions = () + type = real + kind = kind_phys + intent = in + optional = F [convert_dry_rho] standard_name = flag_for_converting_hydrometeors_from_moist_to_dry_air long_name = flag for converting hydrometeors from moist to dry air