From 5a34c967a5e28c8bd0ab8de02f65427023d38f7b Mon Sep 17 00:00:00 2001 From: Samuel Trahan Date: Tue, 31 Jan 2023 23:03:04 +0000 Subject: [PATCH 01/11] GSL lightning threat index --- physics/maximum_hourly_diagnostics.F90 | 82 +++++++++++++++++++++++-- physics/maximum_hourly_diagnostics.meta | 71 +++++++++++++++++++++ 2 files changed, 149 insertions(+), 4 deletions(-) diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index fb3a400e6..373c09467 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -29,11 +29,13 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, gt0, refl_10cm, refdmax, refdmax263k, u10m, v10m, & u10max, v10max, spd10max, pgr, t2m, q2m, t02max, & t02min, rh02max, rh02min, dtp, rain, pratemax, & + lightning_threat, ltg1_max,ltg2_max,ltg3_max, & + vvl, phii, qgraupel, qsnowwat, qicewat, & errmsg, errflg) ! Interface variables integer, intent(in) :: im, levs - logical, intent(in) :: reset, lradar + logical, intent(in) :: reset, lradar, lightning_threat integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_fer_hires, & imp_physics_nssl real(kind_phys), intent(in ) :: con_g @@ -55,20 +57,28 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, real(kind_phys), intent(inout) :: rh02max(:) real(kind_phys), intent(inout) :: rh02min(:) real(kind_phys), intent(in ) :: dtp - real(kind_phys), intent(in ) :: rain(im) - real(kind_phys), intent(inout) :: pratemax(im) + real(kind_phys), intent(in ) :: rain(:) + real(kind_phys), intent(inout) :: pratemax(:) + + real(kind_phys), intent(in), dimension(:,:) :: phii, qgraupel, qsnowwat, qicewat, vvl + real(kind_phys), intent(inout), dimension(:) :: ltg1_max, ltg2_max, ltg3_max character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg ! Local variables real(kind_phys), dimension(:), allocatable :: refd, refd263k - real(kind_phys) :: tem, pshltr, QCQ, rh02 + real(kind_phys) :: tem, pshltr, QCQ, rh02, dP, Q integer :: i ! Initialize CCPP error handling variables errmsg = '' errflg = 0 +!Lightning threat indices + if (lightning_threat) then + call lightning_threat_indices + endif + !Calculate hourly max 1-km agl and -10C reflectivity if (lradar .and. (imp_physics == imp_physics_gfdl .or. & imp_physics == imp_physics_thompson .or. & @@ -134,6 +144,70 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, pratemax(i) = max(pratemax(i),(3.6E6/dtp)*rain(i)) enddo + contains + + subroutine lightning_threat_indices + implicit none + REAL(kind_phys), PARAMETER :: clim1=1.50 + REAL(kind_phys), PARAMETER :: clim2=0.40*1.22 + REAL(kind_phys), PARAMETER :: clim3=0.02*1.22 + ! coef1 and coef2 are modified from the values given + ! in McCaul et al. + ! coef1 is x 1000 x 1.22 + ! coef2 is x 1.22 + ! are these tuning factors, scale factors?? + ! McCaul et al. used a 2-km WRF simulation + REAL(kind_phys), PARAMETER :: coef1=0.042*1000.*1.22 + REAL(kind_phys), PARAMETER :: coef2=0.20*1.22 + + REAL(kind_phys) :: totice_colint(im), msft(im), ltg1, ltg2 + LOGICAL :: ltg1_calc(im) + integer :: k, i + + totice_colint = 0 + ltg1_calc = .false. + msft = 1. + ! get area (m^2) in units of km^2 + ! msft = 1.E-6*area + do k=levs,2,-1 + do i=1,im + dP = phii(i,k) - phii(i,k+1) + Q = qgraupel(i,k) + qsnowwat(i,k) + qicewat(i,k) + totice_colint(i) = totice_colint(i) + Q * dP / con_g + + IF ( .not.ltg1_calc(i) ) THEN + IF ( 0.5*(phii(i,k-1) - phii(i,k+1)) < 258.15 ) THEN + ltg1_calc(i) = .true. + + ltg1 = coef1*vvl(i,k)* & + (( qgraupel(i,k-1) + qgraupel(i,k) )*0.5 )/msft(i) + + IF ( ltg1 .LT. clim1 ) ltg1 = 0. + + IF ( ltg1 .GT. ltg1_max(i) ) THEN + ltg1_max(i) = ltg1 + ENDIF + ENDIF + ENDIF + enddo + enddo + + do i=1,im + ltg2 = coef2 * totice_colint(i) / msft(i) + + IF ( ltg2 .LT. clim2 ) ltg2 = 0. + + IF ( ltg2 .GT. ltg2_max(i) ) THEN + ltg2_max(i) = ltg2 + ENDIF + + ltg3_max(i) = 0.95 * ltg1_max(i) + 0.05 * ltg2_max(i) + + IF ( ltg3_max(i) .LT. clim3 ) ltg3_max(i) = 0. + enddo + + end subroutine lightning_threat_indices + end subroutine maximum_hourly_diagnostics_run subroutine max_fields(phil,ref3D,grav,im,levs,refd,tk,refd263k) diff --git a/physics/maximum_hourly_diagnostics.meta b/physics/maximum_hourly_diagnostics.meta index 391dbde52..33f7eb8f0 100644 --- a/physics/maximum_hourly_diagnostics.meta +++ b/physics/maximum_hourly_diagnostics.meta @@ -238,6 +238,77 @@ type = real kind = kind_phys intent = inout +[vvl] + standard_name = lagrangian_tendency_of_air_pressure + long_name = layer mean vertical velocity + units = Pa s-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[qgraupel] + standard_name = graupel_mixing_ratio + long_name = ratio of mass of graupel to mass of dry air plus vapor (without condensates) + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[qsnowwat] + standard_name = snow_mixing_ratio + long_name = ratio of mass of snow water to mass of dry air plus vapor (without condensates) + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[qicewat] + standard_name = cloud_ice_mixing_ratio + long_name = ratio of mass of ice water to mass of dry air plus vapor (without condensates) + units = kg kg-1 + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in +[lightning_threat] + standard_name = lightning_threat_indices_enabled + long_name = lightning threat indices enabled + units = flag + dimensions = () + type = logical + intent = in +[ltg1_max] + standard_name = gsl_lightning_threat_index_1 + long_name = GSL lightning threat index 1 + units = flashes 5 min-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[ltg2_max] + standard_name = gsl_lightning_threat_index_2 + long_name = GSL lightning threat index 2 + units = flashes 5 min-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[ltg3_max] + standard_name = gsl_lightning_threat_index_3 + long_name = GSL lightning threat index 3 + units = flashes 5 min-1 + dimensions = (horizontal_loop_extent) + type = real + kind = kind_phys + intent = inout +[phii] + standard_name = geopotential_at_interface + long_name = geopotential at model layer interfaces + units = m2 s-2 + dimensions = (horizontal_loop_extent,vertical_interface_dimension) + type = real + kind = kind_phys + intent = in [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 07159470f9254c2dd2a33c46011a5e40aa3a15e9 Mon Sep 17 00:00:00 2001 From: Samuel Trahan Date: Wed, 1 Feb 2023 20:40:32 +0000 Subject: [PATCH 02/11] many changes; code almost works --- physics/maximum_hourly_diagnostics.F90 | 39 +++++++++++++++++++------ physics/maximum_hourly_diagnostics.meta | 15 +++++++--- 2 files changed, 41 insertions(+), 13 deletions(-) diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index 373c09467..f03ad330f 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -30,11 +30,11 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, u10max, v10max, spd10max, pgr, t2m, q2m, t02max, & t02min, rh02max, rh02min, dtp, rain, pratemax, & lightning_threat, ltg1_max,ltg2_max,ltg3_max, & - vvl, phii, qgraupel, qsnowwat, qicewat, & - errmsg, errflg) + wgrs, phii, qgraupel, qsnowwat, qicewat, & + kdt, errmsg, errflg) ! Interface variables - integer, intent(in) :: im, levs + integer, intent(in) :: im, levs, kdt logical, intent(in) :: reset, lradar, lightning_threat integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_fer_hires, & imp_physics_nssl @@ -60,7 +60,7 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, real(kind_phys), intent(in ) :: rain(:) real(kind_phys), intent(inout) :: pratemax(:) - real(kind_phys), intent(in), dimension(:,:) :: phii, qgraupel, qsnowwat, qicewat, vvl + real(kind_phys), intent(in), dimension(:,:) :: phii, qgraupel, qsnowwat, qicewat, wgrs real(kind_phys), intent(inout), dimension(:) :: ltg1_max, ltg2_max, ltg3_max character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -76,6 +76,7 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, !Lightning threat indices if (lightning_threat) then + print *,'call lightning_threat_indices' call lightning_threat_indices endif @@ -160,27 +161,38 @@ subroutine lightning_threat_indices REAL(kind_phys), PARAMETER :: coef1=0.042*1000.*1.22 REAL(kind_phys), PARAMETER :: coef2=0.20*1.22 - REAL(kind_phys) :: totice_colint(im), msft(im), ltg1, ltg2 + REAL(kind_phys) :: totice_colint(im), msft(im), ltg1, ltg2, high_ltg1, high_wgrs, high_graupel LOGICAL :: ltg1_calc(im) - integer :: k, i + integer :: k, i, count + + count = 0 + high_ltg1 = 0 + high_wgrs = 0 + high_graupel = 0 totice_colint = 0 ltg1_calc = .false. msft = 1. ! get area (m^2) in units of km^2 ! msft = 1.E-6*area - do k=levs,2,-1 + do k=2,levs do i=1,im - dP = phii(i,k) - phii(i,k+1) + dP = phii(i,k+1) - phii(i,k) Q = qgraupel(i,k) + qsnowwat(i,k) + qicewat(i,k) totice_colint(i) = totice_colint(i) + Q * dP / con_g IF ( .not.ltg1_calc(i) ) THEN IF ( 0.5*(phii(i,k-1) - phii(i,k+1)) < 258.15 ) THEN + count = count + 1 ltg1_calc(i) = .true. - ltg1 = coef1*vvl(i,k)* & + ltg1 = coef1*wgrs(i,k)* & (( qgraupel(i,k-1) + qgraupel(i,k) )*0.5 )/msft(i) + high_ltg1 = max(high_ltg1, ltg1) + high_graupel = max(high_graupel, qgraupel(i,k)) + if(abs(wgrs(i,k)) > high_wgrs) then + high_wgrs = wgrs(i,k) + endif IF ( ltg1 .LT. clim1 ) ltg1 = 0. @@ -192,6 +204,15 @@ subroutine lightning_threat_indices enddo enddo + if(count > 0) then + if(high_ltg1 == 0 .and. high_wgrs == 0 .and. high_graupel == 0) then + print *, 'high ltg1, wgrs, and graupel are all 0' + else +183 format('high_ltg1 = ',F30.23,' high_wgrs = ',F30.23,' high_graupel = ',F30.23) + print 183, high_ltg1, high_wgrs, high_graupel + endif + endif + do i=1,im ltg2 = coef2 * totice_colint(i) / msft(i) diff --git a/physics/maximum_hourly_diagnostics.meta b/physics/maximum_hourly_diagnostics.meta index 33f7eb8f0..107281a48 100644 --- a/physics/maximum_hourly_diagnostics.meta +++ b/physics/maximum_hourly_diagnostics.meta @@ -238,10 +238,10 @@ type = real kind = kind_phys intent = inout -[vvl] - standard_name = lagrangian_tendency_of_air_pressure - long_name = layer mean vertical velocity - units = Pa s-1 +[wgrs] + standard_name = z_wind + long_name = vertical wind + units = m s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real kind = kind_phys @@ -309,6 +309,13 @@ type = real kind = kind_phys intent = in +[kdt] + standard_name = index_of_timestep + long_name = current forecast iteration + units = index + dimensions = () + type = integer + intent = in [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 4e66f9f23ccc2948ff7eab983a0671dba1e7f332 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Wed, 1 Feb 2023 22:55:33 +0000 Subject: [PATCH 03/11] revise a debug print --- physics/maximum_hourly_diagnostics.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index f03ad330f..08f793aad 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -205,8 +205,8 @@ subroutine lightning_threat_indices enddo if(count > 0) then - if(high_ltg1 == 0 .and. high_wgrs == 0 .and. high_graupel == 0) then - print *, 'high ltg1, wgrs, and graupel are all 0' + if(abs(high_wgrs) < 0.1 .or. high_graupel < 1e-4) then + !print *, 'low wgrs or graupel' else 183 format('high_ltg1 = ',F30.23,' high_wgrs = ',F30.23,' high_graupel = ',F30.23) print 183, high_ltg1, high_wgrs, high_graupel From 59cc88e9177e872639acc0397941f32c334faa73 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Fri, 10 Feb 2023 00:59:10 +0000 Subject: [PATCH 04/11] wrong vars & units --- physics/maximum_hourly_diagnostics.F90 | 26 +++++++++++--------- physics/maximum_hourly_diagnostics.meta | 32 +++++++++++++++++++++---- 2 files changed, 43 insertions(+), 15 deletions(-) diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index 08f793aad..19f110767 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -30,8 +30,8 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, u10max, v10max, spd10max, pgr, t2m, q2m, t02max, & t02min, rh02max, rh02min, dtp, rain, pratemax, & lightning_threat, ltg1_max,ltg2_max,ltg3_max, & - wgrs, phii, qgraupel, qsnowwat, qicewat, & - kdt, errmsg, errflg) + wgrs, prsi, qgraupel, qsnowwat, qicewat, tgrs, con_rd,& + prsl, kdt, errmsg, errflg) ! Interface variables integer, intent(in) :: im, levs, kdt @@ -39,6 +39,7 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_fer_hires, & imp_physics_nssl real(kind_phys), intent(in ) :: con_g + real(kind_phys), intent(in ) :: con_rd real(kind_phys), intent(in ) :: phil(:,:) real(kind_phys), intent(in ) :: gt0(:,:) real(kind_phys), intent(in ) :: refl_10cm(:,:) @@ -58,9 +59,11 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, real(kind_phys), intent(inout) :: rh02min(:) real(kind_phys), intent(in ) :: dtp real(kind_phys), intent(in ) :: rain(:) + real(kind_phys), intent(in ) :: tgrs(:,:) + real(kind_phys), intent(in ) :: prsl(:,:) real(kind_phys), intent(inout) :: pratemax(:) - real(kind_phys), intent(in), dimension(:,:) :: phii, qgraupel, qsnowwat, qicewat, wgrs + real(kind_phys), intent(in), dimension(:,:) :: prsi, qgraupel, qsnowwat, qicewat, wgrs real(kind_phys), intent(inout), dimension(:) :: ltg1_max, ltg2_max, ltg3_max character(len=*), intent(out) :: errmsg integer, intent(out) :: errflg @@ -161,7 +164,7 @@ subroutine lightning_threat_indices REAL(kind_phys), PARAMETER :: coef1=0.042*1000.*1.22 REAL(kind_phys), PARAMETER :: coef2=0.20*1.22 - REAL(kind_phys) :: totice_colint(im), msft(im), ltg1, ltg2, high_ltg1, high_wgrs, high_graupel + REAL(kind_phys) :: totice_colint(im), msft(im), ltg1, ltg2, high_ltg1, high_wgrs, high_graupel, rho LOGICAL :: ltg1_calc(im) integer :: k, i, count @@ -175,19 +178,20 @@ subroutine lightning_threat_indices msft = 1. ! get area (m^2) in units of km^2 ! msft = 1.E-6*area - do k=2,levs + do k=1,levs-1 do i=1,im - dP = phii(i,k+1) - phii(i,k) + dP = prsi(i,k) - prsi(i,k+1) Q = qgraupel(i,k) + qsnowwat(i,k) + qicewat(i,k) - totice_colint(i) = totice_colint(i) + Q * dP / con_g + rho = prsl(i,k) / (con_rd * tgrs(i,k)) + totice_colint(i) = totice_colint(i) + Q * rho * dP / con_g IF ( .not.ltg1_calc(i) ) THEN - IF ( 0.5*(phii(i,k-1) - phii(i,k+1)) < 258.15 ) THEN + IF ( 0.5*(tgrs(i,k+1) + tgrs(i,k)) < 258.15 ) THEN count = count + 1 ltg1_calc(i) = .true. ltg1 = coef1*wgrs(i,k)* & - (( qgraupel(i,k-1) + qgraupel(i,k) )*0.5 )/msft(i) + (( qgraupel(i,k+1) + qgraupel(i,k) )*0.5 )/msft(i) high_ltg1 = max(high_ltg1, ltg1) high_graupel = max(high_graupel, qgraupel(i,k)) if(abs(wgrs(i,k)) > high_wgrs) then @@ -205,8 +209,8 @@ subroutine lightning_threat_indices enddo if(count > 0) then - if(abs(high_wgrs) < 0.1 .or. high_graupel < 1e-4) then - !print *, 'low wgrs or graupel' + if(high_ltg1 < .01 .and. (abs(high_wgrs) < 0.1 .or. high_graupel < 1e-4)) then + ! Nothing to look at else 183 format('high_ltg1 = ',F30.23,' high_wgrs = ',F30.23,' high_graupel = ',F30.23) print 183, high_ltg1, high_wgrs, high_graupel diff --git a/physics/maximum_hourly_diagnostics.meta b/physics/maximum_hourly_diagnostics.meta index 107281a48..98d30dc19 100644 --- a/physics/maximum_hourly_diagnostics.meta +++ b/physics/maximum_hourly_diagnostics.meta @@ -270,6 +270,22 @@ type = real kind = kind_phys intent = in +[con_rd] + standard_name = gas_constant_of_dry_air + long_name = ideal gas constant for dry air + units = J kg-1 K-1 + dimensions = () + type = real + kind = kind_phys + intent = in +[tgrs] + standard_name = air_temperature + long_name = model layer mean temperature + units = K + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in [lightning_threat] standard_name = lightning_threat_indices_enabled long_name = lightning threat indices enabled @@ -301,14 +317,22 @@ type = real kind = kind_phys intent = inout -[phii] - standard_name = geopotential_at_interface - long_name = geopotential at model layer interfaces - units = m2 s-2 +[prsi] + standard_name = air_pressure_at_interface + long_name = air pressure at model layer interfaces + units = Pa dimensions = (horizontal_loop_extent,vertical_interface_dimension) type = real kind = kind_phys intent = in +[prsl] + standard_name = air_pressure + long_name = mean layer pressure + units = Pa + dimensions = (horizontal_loop_extent,vertical_layer_dimension) + type = real + kind = kind_phys + intent = in [kdt] standard_name = index_of_timestep long_name = current forecast iteration From 1a9d3d2e2522ef6d6b3697be6d1aeee3521a2cd7 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Fri, 10 Feb 2023 18:43:06 +0000 Subject: [PATCH 05/11] remove msft and tweak print statements --- physics/maximum_hourly_diagnostics.F90 | 23 ++++++++++++----------- 1 file changed, 12 insertions(+), 11 deletions(-) diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index 19f110767..969af8dcd 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -164,7 +164,7 @@ subroutine lightning_threat_indices REAL(kind_phys), PARAMETER :: coef1=0.042*1000.*1.22 REAL(kind_phys), PARAMETER :: coef2=0.20*1.22 - REAL(kind_phys) :: totice_colint(im), msft(im), ltg1, ltg2, high_ltg1, high_wgrs, high_graupel, rho + REAL(kind_phys) :: totice_colint(im), ltg1, ltg2, high_ltg1, high_wgrs, high_graupel, rho LOGICAL :: ltg1_calc(im) integer :: k, i, count @@ -175,9 +175,6 @@ subroutine lightning_threat_indices totice_colint = 0 ltg1_calc = .false. - msft = 1. - ! get area (m^2) in units of km^2 - ! msft = 1.E-6*area do k=1,levs-1 do i=1,im dP = prsi(i,k) - prsi(i,k+1) @@ -191,11 +188,15 @@ subroutine lightning_threat_indices ltg1_calc(i) = .true. ltg1 = coef1*wgrs(i,k)* & - (( qgraupel(i,k+1) + qgraupel(i,k) )*0.5 )/msft(i) - high_ltg1 = max(high_ltg1, ltg1) - high_graupel = max(high_graupel, qgraupel(i,k)) - if(abs(wgrs(i,k)) > high_wgrs) then - high_wgrs = wgrs(i,k) + (( qgraupel(i,k+1) + qgraupel(i,k) )*0.5 ) + if(ltg1 > 0.01) then +184 format('Found ltg1=',F20.13,' with w=',F20.13,' Qg=',F20.13) + print 184, ltg1, wgrs(i,k), ( qgraupel(i,k+1) + qgraupel(i,k) )*0.5 + endif + if(ltg1 > high_ltg1) then + high_ltg1 = ltg1 + high_graupel = qgraupel(i,k) + high_wgrs = wgrs(i,k) endif IF ( ltg1 .LT. clim1 ) ltg1 = 0. @@ -212,13 +213,13 @@ subroutine lightning_threat_indices if(high_ltg1 < .01 .and. (abs(high_wgrs) < 0.1 .or. high_graupel < 1e-4)) then ! Nothing to look at else -183 format('high_ltg1 = ',F30.23,' high_wgrs = ',F30.23,' high_graupel = ',F30.23) +183 format('Max ltg1=',F20.13,' has w=',F20.13,' Qg=',F20.13) print 183, high_ltg1, high_wgrs, high_graupel endif endif do i=1,im - ltg2 = coef2 * totice_colint(i) / msft(i) + ltg2 = coef2 * totice_colint(i) IF ( ltg2 .LT. clim2 ) ltg2 = 0. From 5d6055ae89e4dbe6f17ed096005e5daf8eafa322 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Thu, 16 Feb 2023 23:13:41 +0000 Subject: [PATCH 06/11] remove prints and climate limits --- physics/maximum_hourly_diagnostics.F90 | 11 +++++------ 1 file changed, 5 insertions(+), 6 deletions(-) diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index 969af8dcd..267dd0d94 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -79,7 +79,6 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, !Lightning threat indices if (lightning_threat) then - print *,'call lightning_threat_indices' call lightning_threat_indices endif @@ -189,7 +188,7 @@ subroutine lightning_threat_indices ltg1 = coef1*wgrs(i,k)* & (( qgraupel(i,k+1) + qgraupel(i,k) )*0.5 ) - if(ltg1 > 0.01) then + if(.false.) then ! ltg1 > 0.01) then 184 format('Found ltg1=',F20.13,' with w=',F20.13,' Qg=',F20.13) print 184, ltg1, wgrs(i,k), ( qgraupel(i,k+1) + qgraupel(i,k) )*0.5 endif @@ -199,7 +198,7 @@ subroutine lightning_threat_indices high_wgrs = wgrs(i,k) endif - IF ( ltg1 .LT. clim1 ) ltg1 = 0. + !IF ( ltg1 .LT. clim1 ) ltg1 = 0. IF ( ltg1 .GT. ltg1_max(i) ) THEN ltg1_max(i) = ltg1 @@ -209,7 +208,7 @@ subroutine lightning_threat_indices enddo enddo - if(count > 0) then + if(.false.) then ! count > 0) then if(high_ltg1 < .01 .and. (abs(high_wgrs) < 0.1 .or. high_graupel < 1e-4)) then ! Nothing to look at else @@ -221,7 +220,7 @@ subroutine lightning_threat_indices do i=1,im ltg2 = coef2 * totice_colint(i) - IF ( ltg2 .LT. clim2 ) ltg2 = 0. + !IF ( ltg2 .LT. clim2 ) ltg2 = 0. IF ( ltg2 .GT. ltg2_max(i) ) THEN ltg2_max(i) = ltg2 @@ -229,7 +228,7 @@ subroutine lightning_threat_indices ltg3_max(i) = 0.95 * ltg1_max(i) + 0.05 * ltg2_max(i) - IF ( ltg3_max(i) .LT. clim3 ) ltg3_max(i) = 0. + !IF ( ltg3_max(i) .LT. clim3 ) ltg3_max(i) = 0. enddo end subroutine lightning_threat_indices From 88a7ebb5f1ab56166212ce10ecf0912092329f3b Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Mon, 27 Feb 2023 22:30:09 +0000 Subject: [PATCH 07/11] put clim limits back in --- physics/maximum_hourly_diagnostics.F90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index 267dd0d94..17679d3c6 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -220,7 +220,7 @@ subroutine lightning_threat_indices do i=1,im ltg2 = coef2 * totice_colint(i) - !IF ( ltg2 .LT. clim2 ) ltg2 = 0. + IF ( ltg2 .LT. clim2 ) ltg2 = 0. IF ( ltg2 .GT. ltg2_max(i) ) THEN ltg2_max(i) = ltg2 @@ -228,7 +228,7 @@ subroutine lightning_threat_indices ltg3_max(i) = 0.95 * ltg1_max(i) + 0.05 * ltg2_max(i) - !IF ( ltg3_max(i) .LT. clim3 ) ltg3_max(i) = 0. + IF ( ltg3_max(i) .LT. clim3 ) ltg3_max(i) = 0. enddo end subroutine lightning_threat_indices From 2942b9f6eec9b96bd504022bb4fce1d83045061d Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Mon, 27 Feb 2023 22:40:54 +0000 Subject: [PATCH 08/11] remove unintended changes --- physics/maximum_hourly_diagnostics.F90 | 15 +-------------- 1 file changed, 1 insertion(+), 14 deletions(-) diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index 17679d3c6..fbbcc86b3 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -188,17 +188,13 @@ subroutine lightning_threat_indices ltg1 = coef1*wgrs(i,k)* & (( qgraupel(i,k+1) + qgraupel(i,k) )*0.5 ) - if(.false.) then ! ltg1 > 0.01) then -184 format('Found ltg1=',F20.13,' with w=',F20.13,' Qg=',F20.13) - print 184, ltg1, wgrs(i,k), ( qgraupel(i,k+1) + qgraupel(i,k) )*0.5 - endif if(ltg1 > high_ltg1) then high_ltg1 = ltg1 high_graupel = qgraupel(i,k) high_wgrs = wgrs(i,k) endif - !IF ( ltg1 .LT. clim1 ) ltg1 = 0. + IF ( ltg1 .LT. clim1 ) ltg1 = 0. IF ( ltg1 .GT. ltg1_max(i) ) THEN ltg1_max(i) = ltg1 @@ -208,15 +204,6 @@ subroutine lightning_threat_indices enddo enddo - if(.false.) then ! count > 0) then - if(high_ltg1 < .01 .and. (abs(high_wgrs) < 0.1 .or. high_graupel < 1e-4)) then - ! Nothing to look at - else -183 format('Max ltg1=',F20.13,' has w=',F20.13,' Qg=',F20.13) - print 183, high_ltg1, high_wgrs, high_graupel - endif - endif - do i=1,im ltg2 = coef2 * totice_colint(i) From 1e5bfd90d1ebf323bd65e80516ef8e41c092fac9 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Tue, 28 Feb 2023 19:11:05 +0000 Subject: [PATCH 09/11] fill ltg*_max with 0 when model is hydrostatic --- physics/maximum_hourly_diagnostics.F90 | 12 +++++++++--- physics/maximum_hourly_diagnostics.meta | 7 +++++++ 2 files changed, 16 insertions(+), 3 deletions(-) diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index fbbcc86b3..df8d9202f 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -31,11 +31,11 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, t02min, rh02max, rh02min, dtp, rain, pratemax, & lightning_threat, ltg1_max,ltg2_max,ltg3_max, & wgrs, prsi, qgraupel, qsnowwat, qicewat, tgrs, con_rd,& - prsl, kdt, errmsg, errflg) + prsl, kdt, hydrostatic, errmsg, errflg) ! Interface variables integer, intent(in) :: im, levs, kdt - logical, intent(in) :: reset, lradar, lightning_threat + logical, intent(in) :: reset, lradar, lightning_threat, hydrostatic integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_fer_hires, & imp_physics_nssl real(kind_phys), intent(in ) :: con_g @@ -79,7 +79,13 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, !Lightning threat indices if (lightning_threat) then - call lightning_threat_indices + if(hydrostatic) then + ltg1_max = 0 + ltg2_max = 0 + ltg3_max = 0 + else + call lightning_threat_indices + endif endif !Calculate hourly max 1-km agl and -10C reflectivity diff --git a/physics/maximum_hourly_diagnostics.meta b/physics/maximum_hourly_diagnostics.meta index 98d30dc19..afe533375 100644 --- a/physics/maximum_hourly_diagnostics.meta +++ b/physics/maximum_hourly_diagnostics.meta @@ -340,6 +340,13 @@ dimensions = () type = integer intent = in +[hydrostatic] + standard_name = flag_for_hydrostatic_solver + long_name = flag for hydrostatic solver from dynamics + units = flag + dimensions = () + type = logical + intent = in [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP From 78568907b5237bb7057a20923fea6506de7c1016 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Mon, 6 Mar 2023 22:36:10 +0000 Subject: [PATCH 10/11] resolve reviewer comments in lightning code --- physics/maximum_hourly_diagnostics.meta | 20 ++++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/physics/maximum_hourly_diagnostics.meta b/physics/maximum_hourly_diagnostics.meta index afe533375..9fa33a667 100644 --- a/physics/maximum_hourly_diagnostics.meta +++ b/physics/maximum_hourly_diagnostics.meta @@ -239,8 +239,8 @@ kind = kind_phys intent = inout [wgrs] - standard_name = z_wind - long_name = vertical wind + standard_name = unsmoothed_nonhydrostatic_upward_air_velocity + long_name = unsmoothed non-hydrostatic upward air velocity units = m s-1 dimensions = (horizontal_loop_extent,vertical_layer_dimension) type = real @@ -287,31 +287,31 @@ kind = kind_phys intent = in [lightning_threat] - standard_name = lightning_threat_indices_enabled - long_name = lightning threat indices enabled + standard_name = do_lightning_threat_index_calculations + long_name = enables the lightning threat index calculations units = flag dimensions = () type = logical intent = in [ltg1_max] - standard_name = gsl_lightning_threat_index_1 - long_name = GSL lightning threat index 1 + standard_name = lightning_threat_index_1 + long_name = lightning threat index 1 units = flashes 5 min-1 dimensions = (horizontal_loop_extent) type = real kind = kind_phys intent = inout [ltg2_max] - standard_name = gsl_lightning_threat_index_2 - long_name = GSL lightning threat index 2 + standard_name = lightning_threat_index_2 + long_name = lightning threat index 2 units = flashes 5 min-1 dimensions = (horizontal_loop_extent) type = real kind = kind_phys intent = inout [ltg3_max] - standard_name = gsl_lightning_threat_index_3 - long_name = GSL lightning threat index 3 + standard_name = lightning_threat_index_3 + long_name = lightning threat index 3 units = flashes 5 min-1 dimensions = (horizontal_loop_extent) type = real From c886b46390752237671c540b7b19d58e03c08534 Mon Sep 17 00:00:00 2001 From: "samuel.trahan" Date: Tue, 14 Mar 2023 20:53:35 +0000 Subject: [PATCH 11/11] remove unneeded hydrostatic check from maximum_hourly_diagnostics --- physics/maximum_hourly_diagnostics.F90 | 12 +++--------- physics/maximum_hourly_diagnostics.meta | 7 ------- 2 files changed, 3 insertions(+), 16 deletions(-) diff --git a/physics/maximum_hourly_diagnostics.F90 b/physics/maximum_hourly_diagnostics.F90 index df8d9202f..cd1016053 100644 --- a/physics/maximum_hourly_diagnostics.F90 +++ b/physics/maximum_hourly_diagnostics.F90 @@ -31,11 +31,11 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, t02min, rh02max, rh02min, dtp, rain, pratemax, & lightning_threat, ltg1_max,ltg2_max,ltg3_max, & wgrs, prsi, qgraupel, qsnowwat, qicewat, tgrs, con_rd,& - prsl, kdt, hydrostatic, errmsg, errflg) + prsl, kdt, errmsg, errflg) ! Interface variables integer, intent(in) :: im, levs, kdt - logical, intent(in) :: reset, lradar, lightning_threat, hydrostatic + logical, intent(in) :: reset, lradar, lightning_threat integer, intent(in) :: imp_physics, imp_physics_gfdl, imp_physics_thompson, imp_physics_fer_hires, & imp_physics_nssl real(kind_phys), intent(in ) :: con_g @@ -79,13 +79,7 @@ subroutine maximum_hourly_diagnostics_run(im, levs, reset, lradar, imp_physics, !Lightning threat indices if (lightning_threat) then - if(hydrostatic) then - ltg1_max = 0 - ltg2_max = 0 - ltg3_max = 0 - else - call lightning_threat_indices - endif + call lightning_threat_indices endif !Calculate hourly max 1-km agl and -10C reflectivity diff --git a/physics/maximum_hourly_diagnostics.meta b/physics/maximum_hourly_diagnostics.meta index 9fa33a667..e9d0876d2 100644 --- a/physics/maximum_hourly_diagnostics.meta +++ b/physics/maximum_hourly_diagnostics.meta @@ -340,13 +340,6 @@ dimensions = () type = integer intent = in -[hydrostatic] - standard_name = flag_for_hydrostatic_solver - long_name = flag for hydrostatic solver from dynamics - units = flag - dimensions = () - type = logical - intent = in [errmsg] standard_name = ccpp_error_message long_name = error message for error handling in CCPP