From 9dc68cadcc5516dae5d81e86e321fde250a6be11 Mon Sep 17 00:00:00 2001 From: gthompsnWRF <35609171+gthompsnWRF@users.noreply.github.com> Date: Wed, 19 Jan 2022 19:16:19 -0700 Subject: [PATCH] After testing with UFS/GFS/FV-3, some tuning knob changes to Thompson-MP and icloud3 (cloud fraction) scheme (#1626) TYPE: enhancement KEYWORDS: microphysics, clouds, cloud fraction, Thompson microphysics SOURCE: Greg Thompson (UCAR/Joint Center for Satellite Data Assimilation) DESCRIPTION OF CHANGES: The purpose of these changes is to produce better cloud coverage and radiation amounts per tuning done with GFS global model and comparisons to observations. Problem: The Thompson-MP code (from ccpp-physics) has been used within UFS/GFS/FV-3 global model for multiple day simulations. Anning Cheng and Ruiyu Sun (NCEP/EMC) have performed most of the simulations and analysis and iterated with me to make adjustments to the Thompson-MP. At first, EMC found that high clouds were not present enough, particularly in the tropics. Once switching to a larger upper limit of cloud ice before it becomes snow, the cloud amounts improved but some clouds and radiation balances were still not matching results from FV3's GFDL microphysics (which had already been tuned to observations). Therefore an effort to account for subgrid-scale clouds using the cloud fraction scheme (icloud=3) was adopted to GFS and tested with various parameter settings. Solution: The biggest change is in how cloud ice converts to snow at a threshold size of 300 microns (up from 200 microns) and for rimed snow to convert to graupel (changed from 250 to 350 microns). A change to allowed max size of ice means that the fall velocity constant for ice was changed to keep it aligned with snow at the same cut-over size (of 300 microns). In addition, the cloud fraction scheme (icloud=3) can further improve the clouds and radiation together with Thompson-MP due to the overall under-prediction of clouds. The icloud3 option was also making far too many clouds when compared to observations so its tuning knobs were adjusted until attaining the following overall improvements compared to cloud and radiation global climatologies that EMC uses: cloud amounts of low, middle, high, and total cloud coverage, longwave radiation outgoing at top-of-atmosphere, and shortwave radiation reaching the ground. ISSUE: There are corresponding issues or pull requests in the ccpp-physics repo, [778](https://github.com/NCAR/ccpp-physics/pull/778) [781](https://github.com/NCAR/ccpp-physics/pull/781) [809](https://github.com/NCAR/ccpp-physics/pull/809) LIST OF MODIFIED FILES: M module_mp_thompson.F M module_radiation_driver.F M module_ra_rrtmg_lw.F M module_ra_rrtmg_lwf.F M module_ra_rrtmg_sw.F M module_ra_rrtmg_swf.F TESTS CONDUCTED: 1. A series of tests in GFS including 5, 7, 16, and 30-day long simulations compared to known cloud and radiation climatologies. The plots attached here in this comment were created by Anning Cheng. The numbers at the top of each panel represent global 3-day (days 3, 4, 5) average. The first plot is the high/mid/low/total cloud amount. The tunings in this PR reduced the high cloud amount from over 55% down to the low 40s, which matches observations of global cloud coverage pretty well. And, another comparison is without the cloud fraction scheme, the mid-level clouds (and low/high) are simply less than observations, which is something well published in IPCC reports of global model clouds, especially the Southern Ocean and east sides of ocean basins (low stratus). ![gfnl_cf_20190611](https://user-images.githubusercontent.com/35609171/149175845-585914ce-84d4-4d72-a7e2-e96a5cb9d228.png) The next plot's panel (a) shows the outgoing longwave radiation. The target value global average is about 240 W/m2. Without the cloud fraction scheme and changes to ```D0s``` and ```D0g```, the result would be too much outgoing longwave since it will not contain enough clouds at lower temperatures. Also, the older version of the cloud fraction scheme with its excessive amount of high clouds would make the outgoing radiation closer to 225 W/m2. So the tuning of the scheme brings the results closer to observations. ![gfnl_prcp_20190611](https://user-images.githubusercontent.com/35609171/149176400-41f8b655-3d58-4cc4-adb0-6f6037beed52.png) Lastly, the final plot shows the downward longwave reaching the ground (panel c) and downward shortwave reaching the ground (panel d). ![gfnl_sfc_SWLW_20190611](https://user-images.githubusercontent.com/35609171/149177060-964ce8bd-02dc-4359-ad5d-37502b2d0f70.png) 2. Jenkins tests are passing. RELEASE NOTE: Update of the Thompson microphysics scheme and cloud fraction scheme (icloud=3) to match the observations better. The modifications include updates to RRTMG LW and SW, and RRTMG fast LW and SW. --- phys/module_mp_thompson.F | 16 +++++----- phys/module_ra_rrtmg_lw.F | 6 ++-- phys/module_ra_rrtmg_lwf.F | 6 ++-- phys/module_ra_rrtmg_sw.F | 6 ++-- phys/module_ra_rrtmg_swf.F | 6 ++-- phys/module_radiation_driver.F | 54 ++++++++++++++++------------------ 6 files changed, 50 insertions(+), 44 deletions(-) diff --git a/phys/module_mp_thompson.F b/phys/module_mp_thompson.F index fd56262c0b..ecc4055a1d 100644 --- a/phys/module_mp_thompson.F +++ b/phys/module_mp_thompson.F @@ -198,8 +198,8 @@ MODULE module_mp_thompson REAL, PARAMETER, PRIVATE:: xm0i = 1.E-12 REAL, PARAMETER, PRIVATE:: D0c = 1.E-6 REAL, PARAMETER, PRIVATE:: D0r = 50.E-6 - REAL, PARAMETER, PRIVATE:: D0s = 200.E-6 - REAL, PARAMETER, PRIVATE:: D0g = 250.E-6 + REAL, PARAMETER, PRIVATE:: D0s = 300.E-6 + REAL, PARAMETER, PRIVATE:: D0g = 350.E-6 REAL, PRIVATE:: D0i, xm0s, xm0g !..Lookup table dimensions @@ -1641,7 +1641,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(9999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + ni(k) = MIN(999.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 @@ -1649,7 +1649,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(9999.D3, cig(1)*oig2*ri(k)/am_i*lami**bm_i) + ni(k) = MIN(999.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 @@ -2678,7 +2678,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(9999.D3, cig(1)*oig2*xri/am_i*lami**bm_i) + xni = MIN(999.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 @@ -2689,8 +2689,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.9999.E3) & - niten(k) = (9999.E3-ni1d(k)*rho(k))*odts*orho + if (xni.gt.999.E3) & + niten(k) = (999.E3-ni1d(k)*rho(k))*odts*orho !..Rain tendency qrten(k) = qrten(k) + (prr_wau(k) + prr_rcw(k) & @@ -3535,7 +3535,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, & - 9999.D3/rho(k)) + 999.D3/rho(k)) endif qr1d(k) = qr1d(k) + qrten(k)*DT nr1d(k) = MAX(R2/rho(k), nr1d(k) + nrten(k)*DT) diff --git a/phys/module_ra_rrtmg_lw.F b/phys/module_ra_rrtmg_lw.F index 718e673f9d..b75f36d3af 100644 --- a/phys/module_ra_rrtmg_lw.F +++ b/phys/module_ra_rrtmg_lw.F @@ -12512,9 +12512,11 @@ SUBROUTINE RRTMG_LWRAD( & if(iceflglw.eq.5)then do k = kts, kte - snow_mass_factor = 1.0 + snow_mass_factor = 0.99 ! Assume 1% of snow overlaps the cloud ice category + gicewp = gicewp + (qs1d(k)*(1.0-snow_mass_factor) * pdel(ncol,k)*100.0 / gravmks * 1000.0) if (resnow1d(ncol,k) .gt. 130.)then - snow_mass_factor = (130.0/resnow1d(ncol,k))*(130.0/resnow1d(ncol,k)) + snow_mass_factor = MIN(snow_mass_factor, & + & (130.0/resnow1d(ncol,k))*(130.0/resnow1d(ncol,k))) resnow1d(ncol,k) = 130.0 IF ( wrf_dm_on_monitor() ) THEN WRITE(message,*)'RRTMG: reducing snow mass (cloud path) to ', & diff --git a/phys/module_ra_rrtmg_lwf.F b/phys/module_ra_rrtmg_lwf.F index 68bcb14de6..ae4b1dea7c 100644 --- a/phys/module_ra_rrtmg_lwf.F +++ b/phys/module_ra_rrtmg_lwf.F @@ -16126,9 +16126,11 @@ SUBROUTINE RRTMG_LWRAD_FAST( & if(iceflglw.eq.5)then do k = kts, kte - snow_mass_factor = 1.0 + snow_mass_factor = 0.99 ! Assume 1% of snow overlaps the cloud ice category + gicewp = gicewp + (qs1d(k)*(1.0-snow_mass_factor) * pdel(ncol,k)*100.0 / gravmks * 1000.0) if (resnow1d(icol,k) .gt. 130.)then - snow_mass_factor = (130.0/resnow1d(icol,k))*(130.0/resnow1d(icol,k)) + snow_mass_factor = MIN(snow_mass_factor, & + & (130.0/resnow1d(ncol,k))*(130.0/resnow1d(ncol,k))) resnow1d(icol,k) = 130.0 IF ( wrf_dm_on_monitor() ) THEN WRITE(message,*)'RRTMG: reducing snow mass (cloud path) to ', & diff --git a/phys/module_ra_rrtmg_sw.F b/phys/module_ra_rrtmg_sw.F index de37e20f71..f8244075c9 100644 --- a/phys/module_ra_rrtmg_sw.F +++ b/phys/module_ra_rrtmg_sw.F @@ -10943,9 +10943,11 @@ SUBROUTINE RRTMG_SWRAD( & if(iceflgsw.eq.5)then do k = kts, kte - snow_mass_factor = 1.0 + snow_mass_factor = 0.99 ! Assume 1% of snow overlaps the cloud ice category + gicewp = gicewp + (qs1d(k)*(1.0-snow_mass_factor) * pdel(ncol,k)*100.0 / gravmks * 1000.0) if (resnow1d(ncol,k) .gt. 130.)then - snow_mass_factor = (130.0/resnow1d(ncol,k))*(130.0/resnow1d(ncol,k)) + snow_mass_factor = MIN(snow_mass_factor, & + & (130.0/resnow1d(ncol,k))*(130.0/resnow1d(ncol,k))) resnow1d(ncol,k) = 130.0 endif gsnowp = qs1d(k) * snow_mass_factor * pdel(ncol,k)*100.0 / gravmks * 1000.0 ! Grid box snow water path. diff --git a/phys/module_ra_rrtmg_swf.F b/phys/module_ra_rrtmg_swf.F index 6440d6389f..ad4f454a7c 100644 --- a/phys/module_ra_rrtmg_swf.F +++ b/phys/module_ra_rrtmg_swf.F @@ -12117,9 +12117,11 @@ SUBROUTINE RRTMG_SWRAD_FAST( & if(iceflgsw.eq.5)then do k = kts, kte - snow_mass_factor = 1.0 + snow_mass_factor = 0.99 ! Assume 1% of snow overlaps the cloud ice category + gicewp = gicewp + (qs1d(k)*(1.0-snow_mass_factor) * pdel(ncol,k)*100.0 / gravmks * 1000.0) if (resnow1d(icol,k) .gt. 130.)then - snow_mass_factor = (130.0/resnow1d(icol,k))*(130.0/resnow1d(icol,k)) + snow_mass_factor = MIN(snow_mass_factor, & + & (130.0/resnow1d(ncol,k))*(130.0/resnow1d(ncol,k))) resnow1d(icol,k) = 130.0 endif gsnowp = qs1d(k) * snow_mass_factor * pdel(icol,k)*100.0 / gravmks * 1000.0 ! Grid box snow water path. diff --git a/phys/module_radiation_driver.F b/phys/module_radiation_driver.F index a80b39b2fa..c63dbee151 100644 --- a/phys/module_radiation_driver.F +++ b/phys/module_radiation_driver.F @@ -3752,13 +3752,16 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & DO k = kts,kte delz = MAX(100., dz(k)) - RH_00L = 0.53 + MIN(0.46,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) - RH_00O = 0.86 + MIN(0.13,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) + RH_00L = 0.77 + MIN(0.22,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) + RH_00O = 0.85 + MIN(0.14,SQRT(1./(50.0+gridkm*gridkm*delz*0.01))) RHUM = rh(k) - if (qc(k).gt.1.E-6 .or. qi(k).ge.1.E-6 & - & .or. (qs(k).gt.1.E-5 .and. t(k).lt.273.)) then + if (qc(k).gt.1.E-6 .or. qi(k).ge.1.E-7 & + & .or. (qs(k).gt.1.E-6 .and. t(k).lt.273.)) then 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)))) else IF ((XLAND-1.5).GT.0.) THEN !--- Ocean @@ -3767,10 +3770,10 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & RH_00 = RH_00L ENDIF - tc = t(k) - 273.15 + tc = MAX(-80.0, t(k) - 273.15) if (tc .lt. -12.0) RH_00 = RH_00L - if (tc .ge. 29.0) then + if (tc .ge. 25.0) then CLDFRA(K) = 0.0 elseif (tc .ge. -12.0) then RHUM = MIN(rh(k), 1.0) @@ -3779,24 +3782,16 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & 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+112.) - if (RH_00 .ge. 1.5) then - WRITE (dbg_msg,*) ' FATAL: RH_00 too large (1.5): ', RH_00, RH_00L, tc - CALL wrf_error_fatal (dbg_msg) - endif + 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))) 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+112.) - if (RH_00 .ge. 1.05) then - WRITE (dbg_msg,*) ' FATAL: RH_00 too large (1.05): ', RH_00, RH_00L, tc - CALL wrf_error_fatal (dbg_msg) - endif + 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))) endif endif - if (CLDFRA(K).gt.0.) CLDFRA(K) = MAX(0.01, MIN(CLDFRA(K),0.95)) + if (CLDFRA(K).gt.0.) CLDFRA(K) = MAX(0.01, MIN(CLDFRA(K),0.99)) if (debug_flag) then WRITE (dbg_msg,*) 'DEBUG-GT: cloud fraction (k,RH_00, RHUM, CF): ',k,RH_00,RHUM,CLDFRA(K) @@ -3804,6 +3799,7 @@ SUBROUTINE cal_cldfra3(CLDFRA, qv, qc, qi, qs, dz, & endif endif + if (cldfra(k).gt.0.0 .and. p(k).lt.7000.0) CLDFRA(K) = 0.0 ENDDO @@ -3990,9 +3986,9 @@ SUBROUTINE find_cloudLayers(qvs1d, cfr1d, T1d, P1d, Dz1d, entrmnt,& ENDDO - k_cldb = k_m12C + 5 + k_cldb = k_m12C + 3 in_cloud = .false. - k = k_m12C + 4 + k = k_m12C + 2 DO WHILE (.not. in_cloud .AND. k.gt.kbot) k_cldt = 0 if (cfr1d(k).ge.0.01) then @@ -4045,13 +4041,14 @@ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte) do k = k1, k2 tdz = tdz + dz(k) enddo - max_iwc = ABS(qvs(k2)-qvs(k1)) +! max_iwc = ABS(qvs(k2)-qvs(k1)) + max_iwc = MAX(0.0, qvs(k1)-qvs(k2)) ! print*, ' max_iwc = ', max_iwc, ' over DZ=',tdz do k = k1, k2 - max_iwc = MAX(1.E-5, max_iwc - (qi(k)+qs(k))) + max_iwc = MAX(1.E-6, max_iwc - (qi(k)+qs(k))) enddo - max_iwc = MIN(2.E-3, max_iwc) + max_iwc = MIN(1.E-4, max_iwc) this_dz = 0.0 do k = k1, k2 @@ -4061,7 +4058,7 @@ SUBROUTINE adjust_cloudIce(cfr,qi,qs,qvs,T,dz,entr, k1,k2,kts,kte) this_dz = this_dz + dz(k) endif this_iwc = max_iwc*this_dz/tdz - iwc = MAX(5.E-6, this_iwc*(1.-entr)) + iwc = MAX(1.E-6, this_iwc*(1.-entr)) if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.203.16) then qi(k) = qi(k) + cfr(k)*cfr(k)*iwc endif @@ -4086,13 +4083,14 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte) do k = k1, k2 tdz = tdz + dz(k) enddo - max_lwc = ABS(qvs(k2)-qvs(k1)) +! max_lwc = ABS(qvs(k2)-qvs(k1)) + max_lwc = MAX(0.0, qvs(k1)-qvs(k2)) ! print*, ' max_lwc = ', max_lwc, ' over DZ=',tdz do k = k1, k2 - max_lwc = MAX(1.E-5, max_lwc - qc(k)) + max_lwc = MAX(1.E-6, max_lwc - qc(k)) enddo - max_lwc = MIN(2.E-3, max_lwc) + max_lwc = MIN(1.E-4, max_lwc) this_dz = 0.0 do k = k1, k2 @@ -4102,8 +4100,8 @@ SUBROUTINE adjust_cloudH2O(cfr, qc, qvs,T,dz,entr, k1,k2,kts,kte) this_dz = this_dz + dz(k) endif this_lwc = max_lwc*this_dz/tdz - lwc = MAX(5.E-6, this_lwc*(1.-entr)) - if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.253.16) then + lwc = MAX(1.E-6, this_lwc*(1.-entr)) + if (cfr(k).gt.0.0.and.cfr(k).lt.1.0.and.T(k).ge.258.16) then qc(k) = qc(k) + cfr(k)*cfr(k)*lwc endif enddo