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
Merged

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

Show file tree
Hide file tree
Changes from all 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
113 changes: 57 additions & 56 deletions physics/module_mp_thompson.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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, &
Expand All @@ -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, &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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, &
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -2220,15 +2225,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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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) &
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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)

Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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
Expand All @@ -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))
Expand All @@ -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
Expand All @@ -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
Expand Down Expand Up @@ -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
Expand All @@ -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
Expand Down
Loading