Skip to content

Thompson MP: new one-moment Graupel y-intercept #965

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

Merged
davegill merged 3 commits into from
Jan 28, 2020
Merged

Thompson MP: new one-moment Graupel y-intercept #965

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
3eb2f2c
large change to graupel Y-intercept parameter (diagnosed), attempting…
gthompsnWRF Jul 22, 2019
0284082
tweak the intercept diagnosis to new final settings.
gthompsnWRF Jul 29, 2019
1a189bb
increment to version2 lookup tables because of expanded number of ele…
gthompsnWRF Aug 2, 2019
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
10 changes: 1 addition & 9 deletions phys/module_diag_misc.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -806,13 +806,6 @@ SUBROUTINE diagnostic_output_calc( &
CASE (THOMPSON, THOMPSONAERO)

scheme_has_graupel = .true.
xmu_g = 1.
cge(1) = xbm_g + 1.
cge(2) = xmu_g + 1.
cge(3) = xbm_g + xmu_g + 1.
do n = 1, 3
cgg(n) = WGAMMA(cge(n))
enddo

! !$OMP PARALLEL DO &
! !$OMP PRIVATE ( ij )
Expand All @@ -821,7 +814,7 @@ SUBROUTINE diagnostic_output_calc( &
DO i=i_start(ij),i_end(ij)
DO k=kme-1, kms, -1
if (temp_qg(i,k,j) .LT. 1.E-6) CYCLE
zans1 = (2.5 + 2./7. * (ALOG10(temp_qg(i,k,j))+7.))
zans1 = (2.5 + 2.5/7. * (ALOG10(temp_qg(i,k,j))+7.))
zans1 = MAX(2., MIN(zans1, 7.))
N0exp = 10.**zans1
lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1))
Expand All @@ -834,7 +827,6 @@ SUBROUTINE diagnostic_output_calc( &
ENDDO
! !$OMP END PARALLEL DO


! CASE (MORR_MILB_P3)
! scheme_has_graupel = .true.
! either Hugh or Jason need to implement.
Expand Down
111 changes: 40 additions & 71 deletions phys/module_mp_thompson.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -107,8 +107,8 @@ MODULE module_mp_thompson
!.. mixing ratio. Also, when mu_g is non-zero, these become equiv
!.. y-intercept for an exponential distrib and proper values are
!.. computed based on same mixing ratio and total number concentration.
REAL, PARAMETER, PRIVATE:: gonv_min = 1.E4
REAL, PARAMETER, PRIVATE:: gonv_max = 3.E6
REAL, PARAMETER, PRIVATE:: gonv_min = 1.E3
REAL, PARAMETER, PRIVATE:: gonv_max = 9.E6

!..Mass power law relations: mass = am*D**bm
!.. Snow from Field et al. (2005), others assume spherical form.
Expand Down Expand Up @@ -211,9 +211,9 @@ MODULE module_mp_thompson
INTEGER, PARAMETER, PRIVATE:: ntb_c = 37
INTEGER, PARAMETER, PRIVATE:: ntb_i = 64
INTEGER, PARAMETER, PRIVATE:: ntb_r = 37
INTEGER, PARAMETER, PRIVATE:: ntb_s = 28
INTEGER, PARAMETER, PRIVATE:: ntb_g = 28
INTEGER, PARAMETER, PRIVATE:: ntb_g1 = 28
INTEGER, PARAMETER, PRIVATE:: ntb_s = 37
INTEGER, PARAMETER, PRIVATE:: ntb_g = 37
INTEGER, PARAMETER, PRIVATE:: ntb_g1 = 37
INTEGER, PARAMETER, PRIVATE:: ntb_r1 = 37
INTEGER, PARAMETER, PRIVATE:: ntb_i1 = 55
INTEGER, PARAMETER, PRIVATE:: ntb_t = 9
Expand Down Expand Up @@ -264,14 +264,16 @@ MODULE module_mp_thompson

!..Lookup tables for graupel content (kg/m**3).
REAL, DIMENSION(ntb_g), PARAMETER, PRIVATE:: &
r_g = (/1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
r_g = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, &
1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, &
1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, &
1.e-2/)

!..Lookup tables for snow content (kg/m**3).
REAL, DIMENSION(ntb_s), PARAMETER, PRIVATE:: &
r_s = (/1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
r_s = (/1.e-6,2.e-6,3.e-6,4.e-6,5.e-6,6.e-6,7.e-6,8.e-6,9.e-6, &
1.e-5,2.e-5,3.e-5,4.e-5,5.e-5,6.e-5,7.e-5,8.e-5,9.e-5, &
1.e-4,2.e-4,3.e-4,4.e-4,5.e-4,6.e-4,7.e-4,8.e-4,9.e-4, &
1.e-3,2.e-3,3.e-3,4.e-3,5.e-3,6.e-3,7.e-3,8.e-3,9.e-3, &
1.e-2/)
Expand All @@ -286,7 +288,8 @@ MODULE module_mp_thompson

!..Lookup tables for graupel y-intercept parameter (/m**4).
REAL, DIMENSION(ntb_g1), PARAMETER, PRIVATE:: &
N0g_exp = (/1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
N0g_exp = (/1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, &
1.e4,2.e4,3.e4,4.e4,5.e4,6.e4,7.e4,8.e4,9.e4, &
1.e5,2.e5,3.e5,4.e5,5.e5,6.e5,7.e5,8.e5,9.e5, &
1.e6,2.e6,3.e6,4.e6,5.e6,6.e6,7.e6,8.e6,9.e6, &
1.e7/)
Expand Down Expand Up @@ -1833,23 +1836,12 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
!+---+-----------------------------------------------------------------+
!..Calculate y-intercept, slope values for graupel.
!+---+-----------------------------------------------------------------+
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
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))
ygra1 = alog10(max(1.E-9, rg(k)))
zans1 = (2.5 + 2.5/7. * (ygra1+7.))
zans1 = MAX(2., MIN(zans1, 7.))
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 @@ -2903,23 +2895,12 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
!+---+-----------------------------------------------------------------+
!..Calculate y-intercept, slope values for graupel.
!+---+-----------------------------------------------------------------+
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
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))
ygra1 = alog10(max(1.E-9, rg(k)))
zans1 = (2.5 + 2.5/7. * (ygra1+7.))
zans1 = MAX(2., MIN(zans1, 7.))
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 @@ -3586,10 +3567,10 @@ subroutine qr_acr_qg

good = 0
IF ( wrf_dm_on_monitor() ) THEN
INQUIRE(FILE="qr_acr_qg.dat",EXIST=lexist)
INQUIRE(FILE="qr_acr_qgV2.dat",EXIST=lexist)
IF ( lexist ) THEN
CALL wrf_message("ThompMP: read qr_acr_qg.dat instead of computing")
OPEN(63,file="qr_acr_qg.dat",form="unformatted",err=1234)
CALL wrf_message("ThompMP: read qr_acr_qgV2.dat instead of computing")
OPEN(63,file="qr_acr_qgV2.dat",form="unformatted",err=1234)
READ(63,err=1234) tcg_racg
READ(63,err=1234) tmr_racg
READ(63,err=1234) tcr_gacr
Expand All @@ -3602,18 +3583,18 @@ subroutine qr_acr_qg
INQUIRE(63,opened=lopen)
IF (lopen) THEN
IF( force_read_thompson ) THEN
CALL wrf_error_fatal("Error reading qr_acr_qg.dat. Aborting because force_read_thompson is .true.")
CALL wrf_error_fatal("Error reading qr_acr_qgV2.dat. Aborting because force_read_thompson is .true.")
ENDIF
CLOSE(63)
ELSE
IF( force_read_thompson ) THEN
CALL wrf_error_fatal("Error opening qr_acr_qg.dat. Aborting because force_read_thompson is .true.")
CALL wrf_error_fatal("Error opening qr_acr_qgV2.dat. Aborting because force_read_thompson is .true.")
ENDIF
ENDIF
ENDIF
ELSE
IF( force_read_thompson ) THEN
CALL wrf_error_fatal("Non-existent qr_acr_qg.dat. Aborting because force_read_thompson is .true.")
CALL wrf_error_fatal("Non-existent qr_acr_qgV2.dat. Aborting because force_read_thompson is .true.")
ENDIF
ENDIF
ENDIF
Expand Down Expand Up @@ -3724,8 +3705,8 @@ subroutine qr_acr_qg
#endif

IF ( write_thompson_tables .AND. wrf_dm_on_monitor() ) THEN
CALL wrf_message("Writing qr_acr_qg.dat in Thompson MP init")
OPEN(63,file="qr_acr_qg.dat",form="unformatted",err=9234)
CALL wrf_message("Writing qr_acr_qgV2.dat in Thompson MP init")
OPEN(63,file="qr_acr_qgV2.dat",form="unformatted",err=9234)
WRITE(63,err=9234) tcg_racg
WRITE(63,err=9234) tmr_racg
WRITE(63,err=9234) tcr_gacr
Expand All @@ -3735,7 +3716,7 @@ subroutine qr_acr_qg
CLOSE(63)
RETURN ! ----- RETURN
9234 CONTINUE
CALL wrf_error_fatal("Error writing qr_acr_qg.dat")
CALL wrf_error_fatal("Error writing qr_acr_qgV2.dat")
ENDIF
ENDIF

Expand Down Expand Up @@ -3772,10 +3753,10 @@ subroutine qr_acr_qs

good = 0
IF ( wrf_dm_on_monitor() ) THEN
INQUIRE(FILE="qr_acr_qs.dat",EXIST=lexist)
INQUIRE(FILE="qr_acr_qsV2.dat",EXIST=lexist)
IF ( lexist ) THEN
CALL wrf_message("ThompMP: read qr_acr_qs.dat instead of computing")
OPEN(63,file="qr_acr_qs.dat",form="unformatted",err=1234)
CALL wrf_message("ThompMP: read qr_acr_qsV2.dat instead of computing")
OPEN(63,file="qr_acr_qsV2.dat",form="unformatted",err=1234)
READ(63,err=1234)tcs_racs1
READ(63,err=1234)tmr_racs1
READ(63,err=1234)tcs_racs2
Expand All @@ -3794,12 +3775,12 @@ subroutine qr_acr_qs
INQUIRE(63,opened=lopen)
IF (lopen) THEN
IF( force_read_thompson ) THEN
CALL wrf_error_fatal("Error reading qr_acr_qs.dat. Aborting because force_read_thompson is .true.")
CALL wrf_error_fatal("Error reading qr_acr_qsV2.dat. Aborting because force_read_thompson is .true.")
ENDIF
CLOSE(63)
ELSE
IF( force_read_thompson ) THEN
CALL wrf_error_fatal("Error opening qr_acr_qs.dat. Aborting because force_read_thompson is .true.")
CALL wrf_error_fatal("Error opening qr_acr_qsV2.dat. Aborting because force_read_thompson is .true.")
ENDIF
ENDIF
ELSE
Expand All @@ -3810,7 +3791,7 @@ subroutine qr_acr_qs
ENDIF
ELSE
IF( force_read_thompson ) THEN
CALL wrf_error_fatal("Non-existent qr_acr_qs.dat. Aborting because force_read_thompson is .true.")
CALL wrf_error_fatal("Non-existent qr_acr_qsV2.dat. Aborting because force_read_thompson is .true.")
ENDIF
ENDIF
ENDIF
Expand Down Expand Up @@ -4003,8 +3984,8 @@ subroutine qr_acr_qs
#endif

IF ( write_thompson_tables .AND. wrf_dm_on_monitor() ) THEN
CALL wrf_message("Writing qr_acr_qs.dat in Thompson MP init")
OPEN(63,file="qr_acr_qs.dat",form="unformatted",err=9234)
CALL wrf_message("Writing qr_acr_qsV2.dat in Thompson MP init")
OPEN(63,file="qr_acr_qsV2.dat",form="unformatted",err=9234)
WRITE(63,err=9234)tcs_racs1
WRITE(63,err=9234)tmr_racs1
WRITE(63,err=9234)tcs_racs2
Expand All @@ -4020,7 +4001,7 @@ subroutine qr_acr_qs
CLOSE(63)
RETURN ! ----- RETURN
9234 CONTINUE
CALL wrf_error_fatal("Error writing qr_acr_qs.dat")
CALL wrf_error_fatal("Error writing qr_acr_qsV2.dat")
ENDIF
ENDIF

Expand Down Expand Up @@ -5279,23 +5260,11 @@ subroutine calc_refl10cm (qv1d, qc1d, qr1d, nr1d, qs1d, qg1d, &
!+---+-----------------------------------------------------------------+

if (ANY(L_qg .eqv. .true.)) then
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
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))
ygra1 = alog10(max(1.E-9, rg(k)))
zans1 = (2.5 + 2.5/7. * (ygra1+7.))
zans1 = MAX(2., MIN(zans1, 7.))
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