Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

After testing with UFS/GFS/FV-3, some tuning knob changes to Thompson-MP and icloud3 (cloud fraction) scheme #1626

Merged
merged 3 commits into from
Jan 20, 2022
Merged

After testing with UFS/GFS/FV-3, some tuning knob changes to Thompson-MP and icloud3 (cloud fraction) scheme #1626

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
8e1b772
a few minor tuning knob changes to Thompson-MP plus tuning of icloud=…
gthompsnWRF Jan 8, 2022
12b6c04
Add RRTMG SW and LW mods to SWF and LWF
davegill Jan 20, 2022
a535d21
52 18:32
davegill Jan 20, 2022
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 8 additions & 8 deletions phys/module_mp_thompson.F
View file
Open in desktop
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
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12432,9 +12432,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
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -16139,9 +16139,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
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -10849,9 +10849,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)))
davegill marked this conversation as resolved.
Show resolved Hide resolved
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
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -12114,9 +12114,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
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -3843,13 +3843,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 @@ -3858,10 +3861,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 @@ -3870,31 +3873,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 @@ -4081,9 +4077,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 @@ -4136,13 +4132,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 @@ -4152,7 +4149,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 @@ -4177,13 +4174,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 @@ -4193,8 +4191,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