Skip to content

Commit

Permalink
Merge pull request #42 from SamuelTrahanNOAA/lightning
Browse files Browse the repository at this point in the history 
Lightning threat indexes
  • Loading branch information
@grantfirl
grantfirl authored Mar 23, 2023
2 parents 630d374 + 521bc28 commit 03acf73
Show file tree
Hide file tree
Showing 2 changed files with 194 additions and 6 deletions.
98 changes: 92 additions & 6 deletions physics/maximum_hourly_diagnostics.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -29,14 +29,17 @@ 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, &
errmsg, errflg)
lightning_threat, ltg1_max,ltg2_max,ltg3_max, &
wgrs, prsi, qgraupel, qsnowwat, qicewat, tgrs, con_rd,&
prsl, kdt, errmsg, errflg)

! Interface variables
integer, intent(in) :: im, levs
logical, intent(in) :: reset, lradar
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
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(:,:)
Expand All @@ -55,20 +58,30 @@ 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(in ) :: tgrs(:,:)
real(kind_phys), intent(in ) :: prsl(:,:)
real(kind_phys), intent(inout) :: pratemax(:)

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

! 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. &
Expand Down Expand Up @@ -134,6 +147,79 @@ 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), ltg1, ltg2, high_ltg1, high_wgrs, high_graupel, rho
LOGICAL :: ltg1_calc(im)
integer :: k, i, count

count = 0
high_ltg1 = 0
high_wgrs = 0
high_graupel = 0

totice_colint = 0
ltg1_calc = .false.
do k=1,levs-1
do i=1,im
dP = prsi(i,k) - prsi(i,k+1)
Q = qgraupel(i,k) + qsnowwat(i,k) + qicewat(i,k)
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*(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 )
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 .GT. ltg1_max(i) ) THEN
ltg1_max(i) = ltg1
ENDIF
ENDIF
ENDIF
enddo
enddo

do i=1,im
ltg2 = coef2 * totice_colint(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)
Expand Down
102 changes: 102 additions & 0 deletions physics/maximum_hourly_diagnostics.meta
View file
Original file line number Diff line number Diff line change
Expand Up @@ -238,6 +238,108 @@
type = real
kind = kind_phys
intent = inout
[wgrs]
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
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
[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 = do_lightning_threat_index_calculations
long_name = enables the lightning threat index calculations
units = flag
dimensions = ()
type = logical
intent = in
[ltg1_max]
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 = 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 = lightning_threat_index_3
long_name = lightning threat index 3
units = flashes 5 min-1
dimensions = (horizontal_loop_extent)
type = real
kind = kind_phys
intent = inout
[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
units = index
dimensions = ()
type = integer
intent = in
[errmsg]
standard_name = ccpp_error_message
long_name = error message for error handling in CCPP
Expand Down

0 comments on commit 03acf73

Please sign in to comment.