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

Minor cloud tunings to improve cloud fraction and radiation in UFS #1

Merged
merged 12 commits into from
Oct 18, 2022
+102 −71
Merged

Minor cloud tunings to improve cloud fraction and radiation in UFS #1

Changes from 6 commits
Commits
Show all changes
12 commits
Select commit Hold shift + click to select a range
c2ac20e
increase max ice number, adjust cloud fraction generally lower
gthompsnWRF Apr 13, 2022
bd59240
adjust/tune for slightly lower reflectivities
gthompsnWRF Apr 20, 2022
c2062bb
permit more ice nucleation, adjust constant in zans (reflectivity) calc
gthompsnWRF Apr 20, 2022
ed0a53d
fix for length of lines
gthompsnWRF Apr 26, 2022
a5eae9f
trying to finalize constant in zans calculation
gthompsnWRF May 6, 2022
14bf941
Merge branch 'small_cloud_tuning_GT' into small_cloud_tuning_GT_gjf
Sep 8, 2022
815438c
Merge branch 'main' into small_cloud_tuning_GT_gjf
Sep 9, 2022
8cfde78
move aero_ind_fdb variable in argument list for mp_gt_driver; add che…
Sep 9, 2022
e9611bd
Merge branch 'ufs/dev' into small_cloud_tuning_GT_gjf
Sep 22, 2022
04f4118
Merge branch 'ufs/dev' into small_cloud_tuning_GT_gjf
Oct 6, 2022
9009159
add additional documentation from Anning Cheng's MERRA2 implementatio…
Oct 6, 2022
353d0ff
Merge branch 'ufs/dev' into small_cloud_tuning_GT_gjf
Oct 14, 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
94 changes: 45 additions & 49 deletions physics/module_mp_thompson.F90
View file
Original file line number Diff line number Diff line change
@@ -2203,15 +2203,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(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
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(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
@@ -2417,7 +2417,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
@@ -2445,12 +2445,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

@@ -2916,7 +2913,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.499.E3) &
if (is_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
@@ -3252,7 +3249,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
@@ -3263,8 +3260,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) &
@@ -3493,7 +3490,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
@@ -3941,7 +3938,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)
@@ -4035,7 +4032,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
@@ -4064,7 +4061,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
@@ -4094,7 +4091,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
@@ -4119,7 +4116,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
@@ -4219,7 +4216,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)
@@ -5684,7 +5681,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

@@ -5727,6 +5724,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)

@@ -5842,6 +5840,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)
@@ -5874,7 +5873,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
@@ -5897,7 +5896,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
@@ -6013,7 +6012,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
@@ -6067,7 +6066,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)
@@ -6090,7 +6090,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
@@ -6196,7 +6197,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)
@@ -6247,12 +6248,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
@@ -6265,18 +6266,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))
@@ -6297,7 +6298,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
@@ -6317,7 +6318,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
@@ -6360,24 +6361,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
@@ -6395,11 +6393,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
33 changes: 20 additions & 13 deletions physics/radiation_clouds.f
View file
Original file line number Diff line number Diff line change
@@ -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