Lightning threat indexes

physics/maximum_hourly_diagnostics.F90

@@ -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, hydrostatic, 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, 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
+       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(:,:)
@@ -55,20 +58,36 @@ 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
+         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
        if (lradar .and. (imp_physics == imp_physics_gfdl .or. &
           imp_physics == imp_physics_thompson .or.    &
@@ -134,6 +153,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)

physics/maximum_hourly_diagnostics.meta

@@ -238,6 +238,115 @@
   type = real
   kind = kind_phys
   intent = inout
+[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
+  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 = 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
+[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
+[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