From 3eb2f2c5826814f695001b358197c715d017f07d Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Mon, 22 Jul 2019 11:22:42 -0600 Subject: [PATCH 1/3] large change to graupel Y-intercept parameter (diagnosed), attempting to follow Field et al (2018JAMC); bug fix to xmu_g=0 in diag_misc along with very minor change to diagnostic Y-intercept to align properly with mp_thompson --- phys/module_diag_misc.F | 10 +----- phys/module_mp_thompson.F | 75 ++++++++++++--------------------------- 2 files changed, 23 insertions(+), 62 deletions(-) diff --git a/phys/module_diag_misc.F b/phys/module_diag_misc.F index 1d323974ee..52cd38be6c 100644 --- a/phys/module_diag_misc.F +++ b/phys/module_diag_misc.F @@ -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 ) @@ -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 = (3.5 + 2./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)) @@ -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. diff --git a/phys/module_mp_thompson.F b/phys/module_mp_thompson.F index b228a50ceb..dadfd7804f 100644 --- a/phys/module_mp_thompson.F +++ b/phys/module_mp_thompson.F @@ -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. @@ -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 @@ -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/) @@ -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/) @@ -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 = (3.5 + 2./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 @@ -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 = (3.5 + 2./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 @@ -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 = (3.5 + 2./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 From 028408271d7f367bf2a7ffcd8cbfe536b2f035f8 Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Mon, 29 Jul 2019 10:50:13 -0600 Subject: [PATCH 2/3] tweak the intercept diagnosis to new final settings. --- phys/module_diag_misc.F | 2 +- phys/module_mp_thompson.F | 6 +++--- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/phys/module_diag_misc.F b/phys/module_diag_misc.F index 52cd38be6c..85deee14fa 100644 --- a/phys/module_diag_misc.F +++ b/phys/module_diag_misc.F @@ -814,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 = (3.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)) diff --git a/phys/module_mp_thompson.F b/phys/module_mp_thompson.F index dadfd7804f..779d342fe4 100644 --- a/phys/module_mp_thompson.F +++ b/phys/module_mp_thompson.F @@ -1839,7 +1839,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.5 + 2./7. * (ygra1+7.)) + zans1 = (2.5 + 2.5/7. * (ygra1+7.)) zans1 = MAX(2., MIN(zans1, 7.)) N0_exp = 10.**(zans1) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 @@ -2898,7 +2898,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.5 + 2./7. * (ygra1+7.)) + zans1 = (2.5 + 2.5/7. * (ygra1+7.)) zans1 = MAX(2., MIN(zans1, 7.)) N0_exp = 10.**(zans1) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 @@ -5262,7 +5262,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.5 + 2./7. * (ygra1+7.)) + zans1 = (2.5 + 2.5/7. * (ygra1+7.)) zans1 = MAX(2., MIN(zans1, 7.)) N0_exp = 10.**(zans1) lam_exp = (N0_exp*am_g*cgg(1)/rg(k))**oge1 From 1a189bb1deac154526d124160702a16da7bfdb6d Mon Sep 17 00:00:00 2001 From: Greg Thompson Date: Fri, 2 Aug 2019 13:28:35 -0600 Subject: [PATCH 3/3] increment to version2 lookup tables because of expanded number of elements of snow, graupel, and graupel Y-intercept values since the older ones are now incompatible. --- phys/module_mp_thompson.F | 36 ++++++++++++++++++------------------ 1 file changed, 18 insertions(+), 18 deletions(-) diff --git a/phys/module_mp_thompson.F b/phys/module_mp_thompson.F index 779d342fe4..744dfa936d 100644 --- a/phys/module_mp_thompson.F +++ b/phys/module_mp_thompson.F @@ -3567,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 @@ -3583,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 @@ -3705,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 @@ -3716,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 @@ -3753,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 @@ -3775,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 @@ -3791,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 @@ -3984,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 @@ -4001,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