Skip to content

Alter graupel intercept #1324

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

Closed
gthompsnWRF wants to merge 9 commits into from
Closed

Alter graupel intercept #1324

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
9 commits
Select commit Hold shift + click to select a range
c38fcce
numerous changes to icloud3 option; improvements after testing with M…
gthompsnWRF Nov 7, 2020
d589ee6
add new namelist option to insert initial cloud fraction and cloud co…
gthompsnWRF Nov 7, 2020
a23deb4
change cldfra to i0rh
gthompsnWRF Nov 7, 2020
2d7f50d
utilize the icloud=3 option as a way to insert initial cloud fraction…
gthompsnWRF Nov 7, 2020
dc593dc
relatively minor tunings of graupel Y-intercept parameter
gthompsnWRF Nov 8, 2020
b1a43a2
another adjustment of constants for Y-intercept graupel, and update t…
gthompsnWRF Nov 9, 2020
ac7f056
another adjustment to intercept parameter
gthompsnWRF Nov 12, 2020
5d686c4
drop in a bug fix for 2 potentially uninitialized vars that hiccup us…
gthompsnWRF Nov 30, 2020
26a57ba
more cases tested, adjusted one value throughout
gthompsnWRF Dec 10, 2020
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
3 changes: 2 additions & 1 deletion Registry/Registry.EM_COMMON
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1551,7 +1551,7 @@ state integer STEPCU - misc 1 - r "S
state real RTHRATEN ikj misc 1 - rd "RTHRATEN" "THETA TENDENCY DUE TO RADIATION" "K s-1"
state real RTHRATENLW ikj misc 1 - r "RTHRATLW" "UNCOUPLED THETA TENDENCY DUE TO LONG WAVE RADIATION" "K s-1"
state real RTHRATENSW ikj misc 1 - r "RTHRATSW" "UNCOUPLED THETA TENDENCY DUE TO SHORT WAVE RADIATION" "K s-1"
state real CLDFRA ikj misc 1 - rh "CLDFRA" "CLOUD FRACTION" ""
state real CLDFRA ikj misc 1 - i0rh "CLDFRA" "CLOUD FRACTION" ""
state real CONVCLD ij misc 1 - r "CONVCLD" "BMJ CONVECTIVE CLOUD" "kg m-2"
state real CCLDFRA ikj misc 1 - r "CCLDFRA" "CONVECTIVE CLOUD FRACTION" ""
state real CLDFRA_OLD ikj misc 1 - r "CLDFRA_OLD" "previous time level cldfra" ""
Expand Down Expand Up @@ -2458,6 +2458,7 @@ rconfig integer tracer_pblmix namelist,physics max_domains 1
rconfig logical use_aero_icbc namelist,physics 1 .false. rh "use_aero_icbc" "Use GOCART climo 3D aerosols IC/BC data in Thompson-MP-Aero" "logical flag"
rconfig logical use_rap_aero_icbc namelist,physics 1 .false. r "use_rap_aero_icbc" "Use GOCART climo 3D aerosols IC/BC data in Thompson-MP-Aero from RAP" "logical flag"
rconfig integer use_mp_re namelist,physics 1 1 h "use_mp_re" "use effective radii computed in some mp schemes in RRTMG" "flag"
rconfig logical insert_init_cloud namelist,physics 1 .false. irh "insert_init_cloud" "Insert initial cloud using icloud=3 method (Thompson)" "logical flag"

# The following two options are hooked into various microphysics schemes to allow for ensemble perturbations of CCN and hail/graupel PSDs - GAC (AFWA)
rconfig real ccn_conc namelist,physics 1 1.0E8 h "ccn_conc" "CCN concentration" "# m-3"
Expand Down
128 changes: 128 additions & 0 deletions dyn_em/module_initialize_real.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -69,6 +69,8 @@ SUBROUTINE init_domain_rk ( grid &
)

USE module_optional_input
USE module_radiation_driver
USE module_dm, ONLY : wrf_dm_max_real
IMPLICIT NONE

! Input space and data. No gridded meteorological data has been stored, though.
Expand Down Expand Up @@ -121,6 +123,13 @@ SUBROUTINE init_domain_rk ( grid &
INTEGER :: auto_levels_opt
REAL :: max_dz, dzbot, dzstretch_s, dzstretch_u, z1, airmass

INTEGER:: i_end, j_end
REAL, ALLOCATABLE, DIMENSION(:):: temp_P, temp_T, temp_R, temp_Dz
REAL, ALLOCATABLE, DIMENSION(:):: temp_Qv, temp_Qc, temp_Qi, temp_Qs
REAL, ALLOCATABLE, DIMENSION(:):: temp_CF, temp_Nc, temp_Ni
REAL:: max_relh, temp_xland, gridkm
LOGICAL :: debug_flag = .FALSE.

REAL :: dclat

! INTEGER , PARAMETER :: nl_max = 1000
Expand Down Expand Up @@ -4370,6 +4379,125 @@ SUBROUTINE init_domain_rk ( grid &
enddo
ENDIF

!+---+-----------------------------------------------------------------+
!..We can consider that in circumstance of a 'cold start' we can make
!.. an attempt to insert some initial clouds to get a better starting
!.. radiation representation due to clouds using the icloud3 cloud fraction
!.. scheme.
!+---+-----------------------------------------------------------------+

if (config_flags%insert_init_cloud .AND. &
(P_QC .gt. 1 .AND. P_QI .gt. 1)) then

ALLOCATE(temp_P(kts:kte-1))
ALLOCATE(temp_Dz(kts:kte-1))
ALLOCATE(temp_T(kts:kte-1))
ALLOCATE(temp_R(kts:kte-1))
ALLOCATE(temp_Qv(kts:kte-1))
ALLOCATE(temp_Qc(kts:kte-1))
ALLOCATE(temp_Nc(kts:kte-1))
ALLOCATE(temp_Qi(kts:kte-1))
ALLOCATE(temp_Ni(kts:kte-1))
ALLOCATE(temp_Qs(kts:kte-1))
ALLOCATE(temp_CF(kts:kte-1))

i_end = MIN(ite,ide-1)
j_end = MIN(jte,jde-1)
max_relh = 1.5
#if ( defined( DM_PARALLEL ) && ( ! defined( STUBMPI ) ) )
max_relh = wrf_dm_max_real ( MAXVAL(grid%u_1(its:i_end,:,jts:j_end)) )
#else
max_relh = MAXVAL ( grid%u_1(its:i_end,:,jts:j_end) )
#endif
max_relh = max_relh*0.01

gridkm = SQRT(config_flags%dx*config_flags%dy)*0.001

!..As it occurs up above, temporarily utilizing the v_1 variable,
!.. to hold temperature, which it does when time_loop=0.

IF ( internal_time_loop .GT. 1 ) THEN
grid%v_1 = grid%t_2+t0
CALL theta_to_t ( grid%v_1 , grid%p_hyd , p00 , &
ids , ide , jds , jde , kds , kde , &
ims , ime , jms , jme , kms , kme , &
its , ite , jts , jte , kts , kte )
ENDIF

do j = jts, j_end
do i = its, i_end
IF ( skip_middle_points_t ( ids , ide , jds , jde , i , j , em_width , hold_ups ) ) CYCLE
debug_flag = .false.
! if (i.eq.9 .and. j.eq.9) debug_flag = .true.

temp_xland = grid%xland(i,j)
if (grid%lakemask(i,j) .eq. 1) temp_xland = 1
do k = kts, kte-1
temp_Dz(k) = (grid%ph_2(i,k+1,j)+grid%phb(i,k+1,j) - (grid%ph_2(i,k,j)+grid%phb(i,k,j)))/g
temp_P(k) = grid%p_hyd(i,k,j)
temp_T(k) = grid%v_1(i,k,j) ! Around line num 1800 v_1 used to hold temperature.
temp_R(k) = 1./grid%alt(i,k,j)
temp_Qv(k) = moist(i,k,j,P_QV)
temp_Qc(k) = MAX(0., moist(i,k,j,P_QC))
temp_Qi(k) = MAX(0., moist(i,k,j,P_QI))
if (P_QS .gt. 1) then
temp_Qs(k) = MAX(0., moist(i,k,j,P_QS))
else
temp_Qs(k) = 0.
endif
if (P_QNI .gt. 1) then
temp_Ni(k) = MAX(0., scalar(i,k,j,P_QNI))
else
temp_Ni(k) = 0.
endif
if (P_QNC .gt. 1) then
temp_Nc(k) = MAX(0., scalar(i,k,j,P_QNC))
else
temp_Nc(k) = 0.
endif
temp_CF(k) = 0.
enddo

call cal_cldfra3(temp_CF,temp_Qv,temp_Qc,temp_Qi,temp_Qs, &
& temp_Dz, temp_P, temp_T, temp_xland, gridkm, &
& config_flags%insert_init_cloud, max_relh, &
& kts, kte-1, debug_flag)

do k = kts, kte-1
grid%cldfra(i,k,j) = temp_CF(k)
enddo

if (debug_flag) then
do k = kts, kte-1
write(*,*) ' DEBUG_column: ', temp_P(k), temp_T(k), temp_Qv(k), temp_Qc(k), temp_Qi(k), moist(i,k,j,P_QC), moist(i,k,j,P_QI)
enddo
endif

if (config_flags%insert_init_cloud) then
do k = kts, kte-1
moist(i,k,j,P_QV) = MAX(temp_Qv(k), moist(i,k,j,P_QV))
moist(i,k,j,P_QC) = temp_Qc(k)
moist(i,k,j,P_QI) = temp_Qi(k)
enddo
endif

enddo
enddo

DEALLOCATE(temp_P)
DEALLOCATE(temp_Dz)
DEALLOCATE(temp_T)
DEALLOCATE(temp_R)
DEALLOCATE(temp_Qv)
DEALLOCATE(temp_Qc)
DEALLOCATE(temp_Nc)
DEALLOCATE(temp_Qi)
DEALLOCATE(temp_Ni)
DEALLOCATE(temp_Qs)
DEALLOCATE(temp_CF)

endif

!+---+-----------------------------------------------------------------+
!..Let us ensure that double-moment microphysics variables have numbers
!.. where there is mass. Currently doing this for Thompson-MP only, but
Expand Down
4 changes: 2 additions & 2 deletions phys/module_diag_nwp.F
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -606,8 +606,8 @@ SUBROUTINE diagnostic_output_nwp( &
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 = (3.5 + 2./7. * (ALOG10(temp_qg(i,k,j))+7.))
zans1 = MAX(2., MIN(zans1, 7.))
zans1 = 3.0 + 2./7.*(ALOG10(temp_qg(i,k,j))+8.)
zans1 = MAX(2., MIN(zans1, 6.))
N0exp = 10.**zans1
lam_exp = (N0exp*xam_g*cgg(1)/temp_qg(i,k,j))**(1./cge(1))
lamg = lam_exp * (cgg(3)/cgg(2)/cgg(1))**(1./xbm_g)
Expand Down
44 changes: 23 additions & 21 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.E3
REAL, PARAMETER, PRIVATE:: gonv_max = 9.E6
REAL, PARAMETER, PRIVATE:: gonv_min = 1.E2
REAL, PARAMETER, PRIVATE:: gonv_max = 1.E7

!..Mass power law relations: mass = am*D**bm
!.. Snow from Field et al. (2005), others assume spherical form.
Expand Down Expand Up @@ -288,11 +288,11 @@ MODULE module_mp_thompson

!..Lookup tables for graupel y-intercept parameter (/m**4).
REAL, DIMENSION(ntb_g1), PARAMETER, PRIVATE:: &
N0g_exp = (/1.e3,2.e3,3.e3,4.e3,5.e3,6.e3,7.e3,8.e3,9.e3, &
N0g_exp = (/1.e2,2.e2,3.e2,4.e2,5.e2,6.e2,7.e2,8.e2,9.e2, &
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/)
1.e6/)

!..Lookup tables for ice number concentration (/m**3).
REAL, DIMENSION(ntb_i1), PARAMETER, PRIVATE:: &
Expand Down Expand Up @@ -1577,7 +1577,7 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
enddo
#endif

!..Bug fix (2016Jun15), prevent use of uninitialized value(s) of snow moments.
!..Bug fix (2016Jun15 and 2020Nov30), prevent use of uninitialized value(s).
do k = kts, kte
smo0(k) = 0.
smo1(k) = 0.
Expand All @@ -1587,6 +1587,8 @@ subroutine mp_thompson (qv1d, qc1d, qi1d, qr1d, qs1d, qg1d, ni1d, &
smod(k) = 0.
smoe(k) = 0.
smof(k) = 0.
mvd_r(k) = 0.
mvd_c(k) = 0.
enddo

!+---+-----------------------------------------------------------------+
Expand Down Expand Up @@ -1848,8 +1850,8 @@ 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.5 + 2./7. * (ygra1+7.))
zans1 = MAX(2., MIN(zans1, 7.))
zans1 = 3.0 + 2./7.*(ygra1+8.)
zans1 = MAX(2., MIN(zans1, 6.))
N0_exp = 10.**(zans1)
lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
Expand Down Expand Up @@ -2909,8 +2911,8 @@ 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.5 + 2./7. * (ygra1+7.))
zans1 = MAX(2., MIN(zans1, 7.))
zans1 = 3.0 + 2./7.*(ygra1+8.)
zans1 = MAX(2., MIN(zans1, 6.))
N0_exp = 10.**(zans1)
lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
Expand Down Expand Up @@ -3578,10 +3580,10 @@ subroutine qr_acr_qg

good = 0
IF ( wrf_dm_on_monitor() ) THEN
INQUIRE(FILE="qr_acr_qgV2.dat",EXIST=lexist)
INQUIRE(FILE="qr_acr_qgV3.dat",EXIST=lexist)
IF ( lexist ) THEN
CALL wrf_message("ThompMP: read qr_acr_qgV2.dat instead of computing")
OPEN(63,file="qr_acr_qgV2.dat",form="unformatted",err=1234)
CALL wrf_message("ThompMP: read qr_acr_qgV3.dat instead of computing")
OPEN(63,file="qr_acr_qgV3.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 @@ -3594,12 +3596,12 @@ subroutine qr_acr_qg
INQUIRE(63,opened=lopen)
IF (lopen) THEN
IF( force_read_thompson ) THEN
CALL wrf_error_fatal("Error reading qr_acr_qgV2.dat. Aborting because force_read_thompson is .true.")
CALL wrf_error_fatal("Error reading qr_acr_qgV3.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_qgV2.dat. Aborting because force_read_thompson is .true.")
CALL wrf_error_fatal("Error opening qr_acr_qgV3.dat. Aborting because force_read_thompson is .true.")
ENDIF
ENDIF
ELSE
Expand All @@ -3610,7 +3612,7 @@ subroutine qr_acr_qg
ENDIF
ELSE
IF( force_read_thompson ) THEN
CALL wrf_error_fatal("Non-existent qr_acr_qgV2.dat. Aborting because force_read_thompson is .true.")
CALL wrf_error_fatal("Non-existent qr_acr_qgV3.dat. Aborting because force_read_thompson is .true.")
ENDIF
ENDIF
ENDIF
Expand Down Expand Up @@ -3722,8 +3724,8 @@ subroutine qr_acr_qg
#endif

IF ( write_thompson_tables .AND. wrf_dm_on_monitor() ) THEN
CALL wrf_message("Writing qr_acr_qgV2.dat in Thompson MP init")
OPEN(63,file="qr_acr_qgV2.dat",form="unformatted",err=9234)
CALL wrf_message("Writing qr_acr_qgV3.dat in Thompson MP init")
OPEN(63,file="qr_acr_qgV3.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 @@ -3733,7 +3735,7 @@ subroutine qr_acr_qg
CLOSE(63)
RETURN ! ----- RETURN
9234 CONTINUE
CALL wrf_error_fatal("Error writing qr_acr_qgV2.dat")
CALL wrf_error_fatal("Error writing qr_acr_qgV3.dat")
ENDIF
ENDIF

Expand Down Expand Up @@ -5294,8 +5296,8 @@ 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.5 + 2./7. * (ygra1+7.))
zans1 = MAX(2., MIN(zans1, 7.))
zans1 = 3.0 + 2./7.*(ygra1+8.)
zans1 = MAX(2., MIN(zans1, 6.))
N0_exp = 10.**(zans1)
lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1
lamg = lam_exp * (cgg(3)*ogg2*ogg1)**obmg
Expand Down
Loading