diff --git a/physics/GFS_MP_generic_pre.F90 b/physics/GFS_MP_generic_pre.F90 index 0910f9cd2..32b097e1d 100644 --- a/physics/GFS_MP_generic_pre.F90 +++ b/physics/GFS_MP_generic_pre.F90 @@ -9,14 +9,14 @@ module GFS_MP_generic_pre !> \section arg_table_GFS_MP_generic_pre_run Argument Table !! \htmlinclude GFS_MP_generic_pre_run.html !! - subroutine GFS_MP_generic_pre_run(im, levs, ldiag3d, qdiag3d, do_aw, ntcw, nncl, & + subroutine GFS_MP_generic_pre_run(im, levs, ldiag3d, qdiag3d, do_aw, progsigma, ntcw, nncl, & ntrac, gt0, gq0, save_t, save_q, num_dfi_radar, errmsg, errflg) ! use machine, only: kind_phys implicit none integer, intent(in) :: im, levs, ntcw, nncl, ntrac, num_dfi_radar - logical, intent(in) :: ldiag3d, qdiag3d, do_aw + logical, intent(in) :: ldiag3d, qdiag3d, do_aw, progsigma real(kind=kind_phys), dimension(:,:), intent(in) :: gt0 real(kind=kind_phys), dimension(:,:,:), intent(in) :: gq0 @@ -39,7 +39,7 @@ subroutine GFS_MP_generic_pre_run(im, levs, ldiag3d, qdiag3d, do_aw, ntcw, nncl, enddo enddo endif - if (ldiag3d .or. do_aw) then + if (ldiag3d .or. do_aw .or. progsigma) then if(qdiag3d) then do n=1,ntrac do k=1,levs @@ -48,7 +48,7 @@ subroutine GFS_MP_generic_pre_run(im, levs, ldiag3d, qdiag3d, do_aw, ntcw, nncl, enddo enddo enddo - else if(do_aw) then + else if(do_aw .or. progsigma) then ! if qdiag3d, all q are saved already save_q(1:im,:,1) = gq0(1:im,:,1) do n=ntcw,ntcw+nncl-1 @@ -59,4 +59,4 @@ subroutine GFS_MP_generic_pre_run(im, levs, ldiag3d, qdiag3d, do_aw, ntcw, nncl, end subroutine GFS_MP_generic_pre_run - end module GFS_MP_generic_pre \ No newline at end of file + end module GFS_MP_generic_pre diff --git a/physics/GFS_MP_generic_pre.meta b/physics/GFS_MP_generic_pre.meta index ac0393917..a2a4947ef 100644 --- a/physics/GFS_MP_generic_pre.meta +++ b/physics/GFS_MP_generic_pre.meta @@ -42,6 +42,13 @@ dimensions = () type = logical intent = in +[progsigma] + standard_name = do_prognostic_updraft_area_fraction + long_name = flag for prognostic area fraction in cumulus convection + units = flag + dimensions = () + type = logical + intent = in [ntcw] standard_name = index_of_cloud_liquid_water_mixing_ratio_in_tracer_concentration_array long_name = tracer index for cloud condensate (or liquid water) diff --git a/physics/lsm_ruc.F90 b/physics/lsm_ruc.F90 index 4db2c9894..99b6c2b41 100644 --- a/physics/lsm_ruc.F90 +++ b/physics/lsm_ruc.F90 @@ -1043,6 +1043,20 @@ subroutine lsm_ruc_run & ! inputs z0_lnd(i,j) = z0rl_lnd(i)/100. znt_lnd(i,j) = z0rl_lnd(i)/100. + ! Workaround needed for subnormal numbers. This should be + ! done after all other sanity checks, in case a sanity check + ! results in subnormal numbers. + ! + ! This bug was caught by the UFS gfortran debug-mode + ! regression tests, and the fix is necessary to pass those + ! tests. + if(abs(snowh_lnd(i,j))<1e-20) then + snowh_lnd(i,j)=0 + endif + if(abs(sneqv_lnd(i,j))<1e-20) then + sneqv_lnd(i,j)=0 + endif + if(debug_print) then if(me==0 ) then write (0,*)'before LSMRUC for land' @@ -1360,6 +1374,19 @@ subroutine lsm_ruc_run & ! inputs z0_ice(i,j) = z0rl_ice(i)/100. znt_ice(i,j) = z0rl_ice(i)/100. + ! Workaround needed for subnormal numbers. This should be + ! done after all other sanity checks, in case a sanity check + ! results in subnormal numbers. + ! + ! Although this bug has not been triggered yet, it is expected + ! to be, like the _lnd variants many lines up from here. + if(abs(snowh_ice(i,j))<1e-20) then + snowh_ice(i,j)=0 + endif + if(abs(sneqv_ice(i,j))<1e-20) then + sneqv_ice(i,j)=0 + endif + !> - Call RUC LSM lsmruc() for ice. call lsmruc( & & delt, flag_init, lsm_cold_start, kdt, iter, nsoil, & diff --git a/physics/module_mp_thompson.F90 b/physics/module_mp_thompson.F90 index b71c8f1f7..fc5a0a2fb 100644 --- a/physics/module_mp_thompson.F90 +++ b/physics/module_mp_thompson.F90 @@ -980,7 +980,7 @@ END SUBROUTINE thompson_init !> @{ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & nwfa, nifa, nwfa2d, nifa2d, & - aero_ind_fdb, tt, th, pii, & + tt, th, pii, & p, w, dz, dt_in, dt_inner, & sedi_semi, decfl, & RAINNC, RAINNCV, & @@ -994,7 +994,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & vt_dbz_wt, first_time_step, & re_cloud, re_ice, re_snow, & has_reqc, has_reqi, has_reqs, & - rand_perturb_on, & + aero_ind_fdb, rand_perturb_on, & kme_stoch, & rand_pert, spp_prt_list, spp_var_list, & spp_stddev_cutoff, n_var_spp, & @@ -1037,7 +1037,6 @@ 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 - LOGICAL, OPTIONAL, INTENT(IN):: aero_ind_fdb REAL, DIMENSION(ims:ime, kms:kme, jms:jme), OPTIONAL, INTENT(INOUT):: & re_cloud, re_ice, re_snow REAL, DIMENSION(ims:ime, kms:kme, jms:jme), INTENT(INOUT):: pfils, pflls @@ -1071,6 +1070,7 @@ SUBROUTINE mp_gt_driver(qv, qc, qr, qi, qs, qg, ni, nr, nc, & LOGICAL, INTENT (IN) :: reset_dBZ ! Extended diagnostics, array pointers only associated if ext_diag flag is .true. LOGICAL, INTENT (IN) :: ext_diag + LOGICAL, OPTIONAL, INTENT(IN):: aero_ind_fdb REAL, DIMENSION(:,:,:), INTENT(INOUT):: & !vts1, txri, txrc, & prw_vcdc, & @@ -1483,10 +1483,15 @@ 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 - if ( .not. aero_ind_fdb) then - nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt - nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt - endif + if ( PRESENT (aero_ind_fdb) ) then + if ( .not. aero_ind_fdb) then + nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt + nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt + endif + else + nwfa1d(kts) = nwfa1d(kts) + nwfa2d(i,j)*dt + nifa1d(kts) = nifa1d(kts) + nifa2d(i,j)*dt + end if do k = kts, kte nc(i,k,j) = nc1d(k) @@ -2220,7 +2225,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & ni(k) = MAX(R2, ni1d(k)*rho(k)) if (ni(k).le. R2) 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(4999.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 @@ -2228,7 +2233,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(4999.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 @@ -2434,7 +2439,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ do k = kte, kts, -1 ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1 + zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 @@ -2462,12 +2467,9 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & do k = kts, kte !> - Rain self-collection follows Seifert, 1994 and drop break-up -!! follows Verlinde and Cotton, 1993. RAIN2M +!! follows Verlinde and Cotton, 1993. Updated after Saleeby et al 2022. RAIN2M if (L_qr(k) .and. mvd_r(k).gt. D0r) then -!-GT Ef_rr = 1.0 -!-GT if (mvd_r(k) .gt. 1500.0E-6) then - Ef_rr = 1.0 - EXP(2300.0*(mvd_r(k)-1950.0E-6)) -!-GT endif + Ef_rr = MAX(-0.1, 1.0 - EXP(2300.0*(mvd_r(k)-1950.0E-6))) pnr_rcr(k) = Ef_rr * 2.0*nr(k)*rr(k) endif @@ -2933,7 +2935,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 .or. merra2_aerosol_aware) .AND. homogIce .AND. (xni.le.499.E3) & + if ((is_aerosol_aware .or. merra2_aerosol_aware) .AND. homogIce .AND. (xni.le.4999.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 @@ -3269,7 +3271,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(4999.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 @@ -3280,8 +3282,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.4999.E3) & + niten(k) = (4999.E3-ni1d(k)*rho(k))*odts*orho !> - Rain tendency qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) & @@ -3510,7 +3512,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & !+---+-----------------------------------------------------------------+ do k = kte, kts, -1 ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1 + zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 @@ -3959,7 +3961,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pfll1(k) = pfll1(k) + sed_r(k)*DT*onstep(1) enddo - if (rr(kts).gt.R1*10.) & + if (rr(kts).gt.R1*1000.) & pptrain = pptrain + sed_r(kts)*DT*onstep(1) enddo else !if(.not. sedi_semi) @@ -4053,7 +4055,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pfil1(k) = pfil1(k) + sed_i(k)*DT*onstep(2) enddo - if (ri(kts).gt.R1*10.) & + if (ri(kts).gt.R1*1000.) & pptice = pptice + sed_i(kts)*DT*onstep(2) enddo endif @@ -4082,7 +4084,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pfil1(k) = pfil1(k) + sed_s(k)*DT*onstep(3) enddo - if (rs(kts).gt.R1*10.) & + if (rs(kts).gt.R1*1000.) & pptsnow = pptsnow + sed_s(kts)*DT*onstep(3) enddo endif @@ -4112,7 +4114,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & pfil1(k) = pfil1(k) + sed_g(k)*DT*onstep(4) enddo - if (rg(kts).gt.R1*10.) & + if (rg(kts).gt.R1*1000.) & pptgraul = pptgraul + sed_g(kts)*DT*onstep(4) enddo else ! if(.not. sedi_semi) then @@ -4137,7 +4139,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, & vtg = 0. if (rg(k).gt. R1) then ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1 + zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 @@ -4237,7 +4239,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)) + 4999.D3/rho(k)) endif qr1d(k) = qr1d(k) + qrten(k)*DT nr1d(k) = MAX(R2/rho(k), nr1d(k) + nrten(k)*DT) @@ -5702,7 +5704,7 @@ end FUNCTION iceKoop !+---+-----------------------------------------------------------------+ !>\ingroup aathompson -!! Helper routine for Phillips et al (2008) ice nucleation. Trude +!! Helper routine for Phillips et al (2008) ice nucleation. REAL FUNCTION delta_p (yy, y1, y2, aa, bb) IMPLICIT NONE @@ -5745,6 +5747,7 @@ END FUNCTION delta_p !! schemes. Since only the smallest snowflakes should impact !! radiation, compute from first portion of complicated Field number !! distribution, not the second part, which is the larger sizes. + subroutine calc_effectRad (t1d, p1d, qv1d, qc1d, nc1d, qi1d, ni1d, qs1d, & & re_qc1d, re_qi1d, re_qs1d, kts, kte) @@ -5860,6 +5863,7 @@ end subroutine calc_effectRad !! library of routines. The meltwater fraction is simply the amount !! of frozen species remaining from what initially existed at the !! melting level interface. + subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & t1d, p1d, dBZ, rand1, kts, kte, ii, jj, melti, & vt_dBZ, first_time_step) @@ -5892,7 +5896,7 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & REAL, DIMENSION(kts:kte):: ze_rain, ze_snow, ze_graupel DOUBLE PRECISION:: N0_exp, N0_min, lam_exp, lamr, lamg - REAL:: a_, b_, loga_, tc0 + REAL:: a_, b_, loga_, tc0, SR DOUBLE PRECISION:: fmelt_s, fmelt_g INTEGER:: i, k, k_0, kbot, n @@ -5915,7 +5919,7 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & else do_vt_dBZ = .false. allow_wet_snow = .true. - allow_wet_graupel = .true. + allow_wet_graupel = .false. endif do k = kts, kte @@ -6031,7 +6035,7 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & if (ANY(L_qg .eqv. .true.)) then do k = kte, kts, -1 ygra1 = alog10(max(1.E-9, rg(k))) - zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1 + zans1 = 3.4 + 2./7.*(ygra1+8.) + rand1 N0_exp = 10.**(zans1) N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max))) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 @@ -6085,7 +6089,8 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & !..Reflectivity contributed by melting snow if (allow_wet_snow .and. L_qs(k) .and. L_qs(k_0) ) then - fmelt_s = MAX(0.05d0, MIN(1.0d0-rs(k)/rs(k_0), 0.99d0)) + SR = MAX(0.01, MIN(1.0 - rs(k)/(rs(k) + rr(k)), 0.99)) + fmelt_s = DBLE(SR*SR) eta = 0.d0 oM3 = 1./smoc(k) M0 = (smob(k)*oM3) @@ -6108,7 +6113,8 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, & !..Reflectivity contributed by melting graupel if (allow_wet_graupel .and. L_qg(k) .and. L_qg(k_0) ) then - fmelt_g = MAX(0.05d0, MIN(1.0d0-rg(k)/rg(k_0), 0.99d0)) + SR = MAX(0.01, MIN(1.0 - rg(k)/(rg(k) + rr(k)), 0.99)) + fmelt_g = DBLE(SR*SR) eta = 0.d0 lamg = 1./ilamg(k) do n = 1, nrbins @@ -6214,7 +6220,7 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) real zsum,qsum,dim,dip,con1,fa1,fa2 real allold, decfl real dz(km), ww(km), qq(km) - real wi(km+1), zi(km+1), za(km+2) !hmhj + real wi(km+1), zi(km+1), za(km+2) real qn(km) real dza(km+1), qa(km+1), qmi(km+1), qpi(km+1) real net_flx(km) @@ -6265,12 +6271,12 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) enddo wi(km) = 0.5*(ww(km)+ww(km-1)) wi(km+1) = ww(km) -! + ! terminate of top of raingroup do k=2,km if( ww(k).eq.0.0 ) wi(k)=ww(k-1) enddo -! + ! diffusivity of wi con1 = 0.05 do k=km,1,-1 @@ -6283,18 +6289,18 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) do k=1,km+1 za(k) = zi(k) - wi(k)*dt enddo - za(km+2) = zi(km+1) !hmhj -! - do k=1,km+1 !hmhj + za(km+2) = zi(km+1) + + do k=1,km+1 dza(k) = za(k+1)-za(k) enddo -! + ! computer deformation at arrival point do k=1,km qa(k) = qq(k)*dz(k)/dza(k) enddo qa(km+1) = 0.0 -! + ! estimate values at arrival cell interface with monotone do k=2,km dip=(qa(k+1)-qa(k))/(dza(k+1)+dza(k)) @@ -6315,7 +6321,7 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) qmi(1)=qa(1) qmi(km+1)=qa(km+1) qpi(km+1)=qa(km+1) -! + ! interpolation to regular point qn = 0.0 kb=1 @@ -6335,7 +6341,7 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) cycle find_kb endif enddo find_kb - find_kt : do kk=kt,km+2 !hmhj + find_kt : do kk=kt,km+2 if( zi(k+1).le.za(kk) ) then kt = kk exit find_kt @@ -6378,24 +6384,21 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) endif cycle intp endif -! + enddo intp -! + ! rain out sum_precip: do k=1,km if( za(k).lt.0.0 .and. za(k+1).le.0.0 ) then -!hmhj precip = precip + qa(k)*dza(k) net_flx(k) = qa(k)*dza(k) cycle sum_precip else if ( za(k).lt.0.0 .and. za(k+1).gt.0.0 ) then -!hmhj -!hmhj precip(i) = precip(i) + qa(k)*(0.0-za(k)) - th = (0.0-za(k))/dza(k) !hmhj - th2 = th*th !hmhj - qqd = 0.5*(qpi(k)-qmi(k)) !hmhj - qqh = qqd*th2+qmi(k)*th !hmhj - precip = precip + qqh*dza(k) !hmhj + th = (0.0-za(k))/dza(k) + th2 = th*th + qqd = 0.5*(qpi(k)-qmi(k)) + qqh = qqd*th2+qmi(k)*th + precip = precip + qqh*dza(k) net_flx(k) = qqh*dza(k) exit sum_precip endif @@ -6413,11 +6416,9 @@ SUBROUTINE semi_lagrange_sedim(km,dzl,wwl,rql,precip,pfsan,dt,R1) ! ! replace the new values rql(:) = max(qn(:),R1) -! -! ---------------------------------- -! + END SUBROUTINE semi_lagrange_sedim -!+---+-----------------------------------------------------------------+ + !+---+-----------------------------------------------------------------+ !+---+-----------------------------------------------------------------+ END MODULE module_mp_thompson diff --git a/physics/mp_thompson.F90 b/physics/mp_thompson.F90 index 80762a5f1..ffe1a03d6 100644 --- a/physics/mp_thompson.F90 +++ b/physics/mp_thompson.F90 @@ -686,7 +686,6 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & if (is_aerosol_aware .or. merra2_aerosol_aware) then 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, & - aero_ind_fdb=aero_ind_fdb, & tt=tgrs, p=prsl, w=w, dz=dz, dt_in=dtstep, dt_inner=dt_inner, & sedi_semi=sedi_semi, decfl=decfl, & rainnc=rain_mp, rainncv=delta_rain_mp, & @@ -696,7 +695,8 @@ subroutine mp_thompson_run(ncol, nlev, con_g, con_rd, & refl_10cm=refl_10cm, & diagflag=diagflag, do_radar_ref=do_radar_ref_mp, & has_reqc=has_reqc, has_reqi=has_reqi, has_reqs=has_reqs, & - rand_perturb_on=spp_mp_opt, kme_stoch=kme_stoch, & + aero_ind_fdb=aero_ind_fdb, rand_perturb_on=spp_mp_opt, & + kme_stoch=kme_stoch, & rand_pert=spp_wts_mp, spp_var_list=spp_var_list, & spp_prt_list=spp_prt_list, n_var_spp=n_var_spp, & spp_stddev_cutoff=spp_stddev_cutoff, & @@ -887,6 +887,29 @@ subroutine get_niwfa(aerfld, nifa, nwfa, ncol, nlev) ! mass. NIFA is mainly summarized over five dust bins and NWFA over the ! other 10 bins. The parameters besides each bins are carefully tuned ! for a good performance of the scheme. + ! + ! The fields for the last index of the aerfld array + ! are specified as below. + ! 1: dust bin 1, 0.1 to 1.0 micrometers + ! 2: dust bin 2, 1.0 to 1.8 micrometers + ! 3: dust bin 3, 1.8 to 3.0 micrometers + ! 4: dust bin 4, 3.0 to 6.0 micrometers + ! 5: dust bin 5, 6.0 to 10.0 micrometers + ! 6: sea salt bin 1, 0.03 to 0.1 micrometers + ! 7: sea salt bin 2, 0.1 to 0.5 micrometers + ! 8: sea salt bin 3, 0.5 to 1.5 micrometers + ! 9: sea salt bin 4, 1.5 to 5.0 micrometers + ! 10: sea salt bin 5, 5.0 to 10.0 micrometers + ! 11: Sulfate, 0.35 (mean) micrometers + ! 15: water-friendly organic carbon, 0.35 (mean) micrometers + ! + ! Bin densities are as follows: + ! 1: dust bin 1: 2500 kg/m2 + ! 2-5: dust bin 2-5: 2650 kg/m2 + ! 6-10: sea salt bins 6-10: 2200 kg/m2 + ! 11: sulfate: 1700 kg/m2 + ! 15: organic carbon: 1800 kg/m2 + implicit none integer, intent(in)::ncol, nlev real (kind=kind_phys), dimension(:,:,:), intent(in) :: aerfld diff --git a/physics/progsigma_calc.f90 b/physics/progsigma_calc.f90 index 49a5e2a4f..f83514f27 100644 --- a/physics/progsigma_calc.f90 +++ b/physics/progsigma_calc.f90 @@ -14,8 +14,8 @@ !> @{ subroutine progsigma_calc (im,km,flag_init,flag_restart, & - del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & - delt,prevsq,q,kbcon1,ktcon,cnvflg,sigmain,sigmaout, & + flag_shallow,del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap, & + delt,prevsq,q,kbcon1,ktcon,cnvflg,sigmain,sigmaout, & sigmab,errmsg,errflg) ! ! @@ -30,7 +30,7 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & real(kind=kind_phys), intent(in) :: prevsq(im,km), q(im,km),del(im,km), & qmicro(im,km),tmf(im,km),dbyo1(im,km),zdqca(im,km), & omega_u(im,km),zeta(im,km) - logical, intent(in) :: flag_init,flag_restart,cnvflg(im) + logical, intent(in) :: flag_init,flag_restart,cnvflg(im),flag_shallow real(kind=kind_phys), intent(in) :: sigmain(im,km) ! intent out @@ -48,14 +48,15 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & real(kind=kind_phys) :: gcvalmx,epsilon,ZZ,cvg,mcon,buy2, & fdqb,dtdyn,dxlim,rmulacvg,tem, & - DEN,betascu,dp1,invdelt + DEN,betascu,betadcu,dp1,invdelt !Parameters gcvalmx = 0.1 rmulacvg=10. epsilon=1.E-11 km1=km-1 - betascu = 3.0 + betadcu = 2.0 + betascu = 3.6 invdelt = 1./delt !Initialization 2D @@ -212,12 +213,19 @@ subroutine progsigma_calc (im,km,flag_init,flag_restart, & !Reduce area fraction before coupling back to mass-flux computation. !This tuning could be addressed in updraft velocity equation instead. - do i= 1, im - if(cnvflg(i)) then - sigmab(i)=sigmab(i)/betascu - endif - enddo - + if(flag_shallow)then + do i= 1, im + if(cnvflg(i)) then + sigmab(i)=sigmab(i)/betascu + endif + enddo + else + do i= 1, im + if(cnvflg(i)) then + sigmab(i)=sigmab(i)/betadcu + endif + enddo + endif end subroutine progsigma_calc diff --git a/physics/radiation_clouds.f b/physics/radiation_clouds.f index edfa94439..f7e017849 100644 --- a/physics/radiation_clouds.f +++ b/physics/radiation_clouds.f @@ -2487,7 +2487,7 @@ subroutine progcld_thompson & do i = 1, IX cwp(i,k) = max(0.0, clw(i,k,ntcw) * dz(i,k)*1.E6) crp(i,k) = 0.0 - snow_mass_factor = 0.99 + snow_mass_factor = 0.90 cip(i,k) = max(0.0, (clw(i,k,ntiw) & & + (1.0-snow_mass_factor)*clw(i,k,ntsw))*dz(i,k)*1.E6) if (re_snow(i,k) .gt. snow_max_radius)then @@ -3378,7 +3378,7 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & REAL:: RH_00L, RH_00O, RH_00 REAL:: entrmnt=0.5 INTEGER:: k - REAL:: TC, qvsi, qvsw, RHUM, delz + REAL:: TC, qvsi, qvsw, RHUM, delz, var_temp REAL, DIMENSION(kts:kte):: qvs, rh, rhoa integer:: ndebug = 0 @@ -3438,7 +3438,8 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & CLDFRA(K) = 1.0 elseif (((qc(k)+qi(k)).gt.1.E-10) .and. & & ((qc(k)+qi(k)).lt.1.E-6)) then - CLDFRA(K) = MIN(0.99, 0.1*(11.0 + log10(qc(k)+qi(k)))) + var_temp = MIN(0.99, 0.1*(11.0 + log10(qc(k)+qi(k)))) + CLDFRA(K) = var_temp*var_temp else IF ((XLAND-1.5).GT.0.) THEN !--- Ocean @@ -3448,24 +3449,27 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & ENDIF tc = MAX(-80.0, t(k) - 273.15) - if (tc .lt. -12.0) RH_00 = RH_00L + if (tc .lt. -24.0) RH_00 = RH_00L if (tc .gt. 20.0) then CLDFRA(K) = 0.0 - elseif (tc .ge. -12.0) then + elseif (tc .ge. -24.0) then RHUM = MIN(rh(k), 1.0) - CLDFRA(K) = MAX(0., 1.0-SQRT((1.001-RHUM)/(1.001-RH_00))) + var_temp = SQRT(SQRT((1.001-RHUM)/(1.001-RH_00))) + CLDFRA(K) = MAX(0., 1.0-var_temp) else if (max_relh.gt.1.12 .or. (.NOT.(modify_qvapor)) ) then !..For HRRR model, the following look OK. RHUM = MIN(rh(k), 1.45) - RH_00 = RH_00 + (1.45-RH_00)*(-12.0-tc)/(-12.0+85.) - CLDFRA(K) = MAX(0.,1.0-SQRT((1.46-RHUM)/(1.46-RH_00))) + RH_00 = RH_00 + (1.45-RH_00)*(-24.0-tc)/(-24.0+85.) + var_temp = SQRT(SQRT((1.46-RHUM)/(1.46-RH_00))) + CLDFRA(K) = MAX(0., 1.0-var_temp) else !..but for the GFS model, RH is way lower. RHUM = MIN(rh(k), 1.05) - RH_00 = RH_00 + (1.05-RH_00)*(-12.0-tc)/(-12.0+85.) - CLDFRA(K) = MAX(0.,1.0-SQRT((1.06-RHUM)/(1.06-RH_00))) + RH_00 = RH_00 + (1.05-RH_00)*(-24.0-tc)/(-24.0+85.) + var_temp = SQRT(SQRT((1.06-RHUM)/(1.06-RH_00))) + CLDFRA(K) = MAX(0., 1.0-var_temp) endif endif if (CLDFRA(K).gt.0.) CLDFRA(K)=MAX(0.01,MIN(CLDFRA(K),0.99)) @@ -3784,6 +3788,8 @@ SUBROUTINE adjust_cloudFinal(cfr, qc, qi, Rho,dz, kts,kte) END SUBROUTINE adjust_cloudFinal +!+---+-----------------------------------------------------------------+ + subroutine cloud_fraction_XuRandall & & ( IX, NLAY, plyr, clwf, rhly, qstl, & ! --- inputs & cldtot ) & ! --- outputs @@ -3808,7 +3814,6 @@ subroutine cloud_fraction_XuRandall & do k = 1, NLAY do i = 1, IX clwt = 1.0e-6 * (plyr(i,k)*0.001) -! clwt = 2.0e-6 * (plyr(i,k)*0.001) if (clwf(i,k) > clwt) then @@ -3818,8 +3823,6 @@ subroutine cloud_fraction_XuRandall & tem1 = min(max(sqrt(sqrt(onemrh*qstl(i,k))),0.0001),1.0) tem1 = 2000.0 / tem1 -! tem1 = 1000.0 / tem1 - value = max( min( tem1*(clwf(i,k)-clwm), 50.0 ), 0.0 ) tem2 = sqrt( sqrt(rhly(i,k)) ) @@ -3830,6 +3833,8 @@ subroutine cloud_fraction_XuRandall & end subroutine cloud_fraction_XuRandall +!+---+-----------------------------------------------------------------+ + subroutine cloud_fraction_mass_flx_1 & & ( IX, NLAY, lmfdeep2, xrc3, plyr, clwf, rhly, qstl, & ! --- inputs & cldtot ) & ! --- outputs @@ -3879,6 +3884,8 @@ subroutine cloud_fraction_mass_flx_1 & end subroutine cloud_fraction_mass_flx_1 +!+---+-----------------------------------------------------------------+ + subroutine cloud_fraction_mass_flx_2 & & ( IX, NLAY, lmfdeep2, xrc3, plyr, clwf, rhly, qstl, & ! --- inputs & cldtot ) & ! --- outputs diff --git a/physics/samfdeepcnv.f b/physics/samfdeepcnv.f index 2552ac622..2bef2d729 100644 --- a/physics/samfdeepcnv.f +++ b/physics/samfdeepcnv.f @@ -212,6 +212,7 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & real(kind=kind_phys) omega_u(im,km),zdqca(im,km),qlks(im,km), & omegac(im),zeta(im,km),dbyo1(im,km),sigmab(im) real(kind=kind_phys) gravinv + logical flag_shallow c physical parameters ! parameter(grav=grav,asolfac=0.958) ! parameter(elocp=hvap/cp,el2orc=hvap*hvap/(rv*cp)) @@ -2882,7 +2883,8 @@ subroutine samfdeepcnv_run (im,km,first_time_step,restart, & !> - From Bengtsson et al. (2022) Prognostic closure scheme, equation 8, compute updraft area fraction based on a moisture budget if(progsigma)then - call progsigma_calc(im,km,first_time_step,restart, + flag_shallow = .false. + call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & prevsq,q,kbcon1,ktcon,cnvflg, & sigmain,sigmaout,sigmab,errmsg,errflg) diff --git a/physics/samfshalcnv.f b/physics/samfshalcnv.f index c4a6b8685..2031c0ae9 100644 --- a/physics/samfshalcnv.f +++ b/physics/samfshalcnv.f @@ -160,7 +160,7 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & & omegac(im),zeta(im,km),dbyo1(im,km), & sigmab(im) real(kind=kind_phys) gravinv,dxcrtas - + logical flag_shallow c physical parameters ! parameter(g=grav,asolfac=0.89) ! parameter(g=grav) @@ -1929,7 +1929,8 @@ subroutine samfshalcnv_run(im,km,itc,ntc,cliq,cp,cvap, & c c Prognostic closure if(progsigma)then - call progsigma_calc(im,km,first_time_step,restart, + flag_shallow = .true. + call progsigma_calc(im,km,first_time_step,restart,flag_shallow, & del,tmf,qmicro,dbyo1,zdqca,omega_u,zeta,hvap,delt, & prevsq,q,kbcon1,ktcon,cnvflg, & sigmain,sigmaout,sigmab,errmsg,errflg)