Skip to content

Commit

Permalink
After testing with UFS/GFS/FV-3, some tuning knob changes to Thompson…
Browse files Browse the repository at this point in the history
…-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](NCAR/ccpp-physics#778)
[781](NCAR/ccpp-physics#781)
[809](NCAR/ccpp-physics#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.
  • Loading branch information
gthompsnWRF authored Jan 20, 2022
1 parent 96fd889 commit 9dc68ca
Show file tree
Hide file tree
Showing 6 changed files with 50 additions and 44 deletions.
16 changes: 8 additions & 8 deletions phys/module_mp_thompson.F
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -1641,15 +1641,15 @@ 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
ilami = 1./lami
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
Expand Down Expand Up @@ -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
Expand All @@ -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) &
Expand Down Expand Up @@ -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)
Expand Down
6 changes: 4 additions & 2 deletions phys/module_ra_rrtmg_lw.F
Original file line number Diff line number Diff line change
Expand Up @@ -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 ', &
Expand Down
6 changes: 4 additions & 2 deletions phys/module_ra_rrtmg_lwf.F
Original file line number Diff line number Diff line change
Expand Up @@ -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 ', &
Expand Down
6 changes: 4 additions & 2 deletions phys/module_ra_rrtmg_sw.F
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
6 changes: 4 additions & 2 deletions phys/module_ra_rrtmg_swf.F
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down
54 changes: 26 additions & 28 deletions phys/module_radiation_driver.F
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -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)
Expand All @@ -3779,31 +3782,24 @@ 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)
CALL wrf_debug (150, dbg_msg)
endif
endif
if (cldfra(k).gt.0.0 .and. p(k).lt.7000.0) CLDFRA(K) = 0.0
ENDDO
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand All @@ -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
Expand All @@ -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
Expand Down

0 comments on commit 9dc68ca

Please sign in to comment.