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

RRFS_dev: fix problems in reflectivity diagnosis for Thompson scheme. #91

Merged
merged 1 commit into from
May 21, 2021
Merged

RRFS_dev: fix problems in reflectivity diagnosis for Thompson scheme. #91

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
85 changes: 73 additions & 12 deletions physics/module_mp_thompson.F90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2048,11 +2048,27 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
!+---+-----------------------------------------------------------------+
!> - Calculate y-intercept, slope values for graupel.
!+---+-----------------------------------------------------------------+
! Ming Hu: go back to old version for Spring experiment 2021
N0_min = gonv_max
k_0 = kts
do k = kte, kts, -1
if (temp(k).ge.270.65) k_0 = MAX(k_0, k)
enddo
do k = kte, kts, -1
ygra1 = alog10(max(1.E-9, rg(k)))
zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1
if (k.gt.k_0 .and. L_qr(k) .and. mvd_r(k).gt.100.E-6) then
xslw1 = 4.01 + alog10(mvd_r(k))
else
xslw1 = 0.01
endif
ygra1 = 4.31 + alog10(max(5.E-5, rg(k)))
zans1 = (3.1 +(100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1))) + rand1
if (rand1 .ne. 0.0) then
zans1 = MAX(2., MIN(zans1, 7.))
endif
N0_exp = 10.**(zans1)
N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max)))
N0_min = MIN(N0_exp, N0_min)
N0_exp = N0_min
lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
ilamg(k) = 1./lamg
Expand Down Expand Up @@ -3124,11 +3140,28 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
!+---+-----------------------------------------------------------------+
!> - Calculate y-intercept, slope values for graupel.
!+---+-----------------------------------------------------------------+
! Ming Hu: go back to old version for Spring experiment 2021

N0_min = gonv_max
k_0 = kts
do k = kte, kts, -1
ygra1 = alog10(max(1.E-9, rg(k)))
zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1
if (temp(k).ge.270.65) k_0 = MAX(k_0, k)
enddo
do k = kte, kts, -1
if (k.gt.k_0 .and. L_qr(k) .and. mvd_r(k).gt.100.E-6) then
xslw1 = 4.01 + alog10(mvd_r(k))
else
xslw1 = 0.01
endif
ygra1 = 4.31 + alog10(max(5.E-5, rg(k)))
zans1 = (3.1 +(100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1))) + rand1
if (rand1 .ne. 0.0) then
zans1 = MAX(2., MIN(zans1, 7.))
endif
N0_exp = 10.**(zans1)
N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max)))
N0_min = MIN(N0_exp, N0_min)
N0_exp = N0_min
lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
ilamg(k) = 1./lamg
Expand Down Expand Up @@ -3483,10 +3516,14 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
t4_vts = Kap1*Mrat**mu_s*csg(7)*ils2**cse(7)
vts = rhof(k)*av_s * (t1_vts+t2_vts)/(t3_vts+t4_vts)
if (temp(k).gt. (T_0+0.1)) then
! vtsk(k) = MAX(vts*vts_boost(k), &
! & vts*((vtrk(k)-vts*vts_boost(k))/(temp(k)-T_0)))
SR = rs(k)/(rs(k)+rr(k))
vtsk(k) = vts*SR + (1.-SR)*vtrk(k)
!Ming Hu: go back to old version for Spring experiment 2021
vtsk(k) = MAX(vts*vts_boost(k), &
& vts*((vtrk(k)-vts*vts_boost(k))/(temp(k)-T_0))) !
! DH* The version below is supposed to be a better formulation,
! but gave worse results in RAPv5/HRRRv4 than the line above.
! this formulation for RAPv5/HRRRv4, reverted 20 Feb 2020
! SR = rs(k)/(rs(k)+rr(k))
! vtsk(k) = vts*SR + (1.-SR)*vtrk(k)
else
vtsk(k) = vts*vts_boost(k)
endif
Expand Down Expand Up @@ -5295,7 +5332,7 @@ end subroutine calc_effectRad
!! 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, &
t1d, p1d, dBZ, rand1, kts, kte, ii, jj, melti_org, &
vt_dBZ, first_time_step)

IMPLICIT NONE
Expand Down Expand Up @@ -5330,7 +5367,8 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
DOUBLE PRECISION:: fmelt_s, fmelt_g

INTEGER:: i, k, k_0, kbot, n
LOGICAL, INTENT(IN):: melti
LOGICAL, INTENT(IN):: melti_org
LOGICAL :: melti
LOGICAL, DIMENSION(kts:kte):: L_qr, L_qs, L_qg

DOUBLE PRECISION:: cback, x, eta, f_d
Expand All @@ -5352,6 +5390,11 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
allow_wet_graupel = .true.
endif

!Ming Hu hardwired for Spring Experiment testing
allow_wet_snow = .true.
allow_wet_graupel = .false.
melti=.true.

do k = kts, kte
dBZ(k) = -35.0
enddo
Expand Down Expand Up @@ -5463,16 +5506,34 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
!+---+-----------------------------------------------------------------+

if (ANY(L_qg .eqv. .true.)) then
! Ming Hu: go back to old version for Spring experiment 2021

N0_min = gonv_max
k_0 = kts
do k = kte, kts, -1
ygra1 = alog10(max(1.E-9, rg(k)))
zans1 = 3.0 + 2./7.*(ygra1+8.) + rand1
if (temp(k).ge.270.65) k_0 = MAX(k_0, k)
enddo
do k = kte, kts, -1
if (k.gt.k_0 .and. L_qr(k) .and. mvd_r(k).gt.100.E-6) then
xslw1 = 4.01 + alog10(mvd_r(k))
else
xslw1 = 0.01
endif
ygra1 = 4.31 + alog10(max(5.E-5, rg(k)))
zans1 = (3.1 +(100./(300.*xslw1*ygra1/(10./xslw1+1.+0.25*ygra1)+30.+10.*ygra1))) + rand1
if (rand1 .ne. 0.0) then
zans1 = MAX(2., MIN(zans1, 7.))
endif
N0_exp = 10.**(zans1)
N0_exp = MAX(DBLE(gonv_min), MIN(N0_exp, DBLE(gonv_max)))
N0_min = MIN(N0_exp, N0_min)
N0_exp = N0_min
lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
ilamg(k) = 1./lamg
N0_g(k) = N0_exp/(cgg(2)*lam_exp) * lamg**cge(2)
enddo

endif

!+---+-----------------------------------------------------------------+
Expand Down