Lightning threat indexes (#628)

* add IFI to post

* pass w to model & fixes to ccpp/physics changes

* send lightning threat indexes to post

* upp: remove 444 ltng and add 704

* ccpp/physics: fill ltg*_max with 0 when model is hydrostatic

* atmos_cubed_sphere: only send w to model in non-hydrostatic mode, and only if model wants it

* disable lightning_threat when hydrostatic, and initialize wgrs

* turn off lightning_threat by default since it allocates wgrs and passes w

* missing `active = (lightning_threat_indices_enabled)` for wgrd

* UPP submodule: correct ifi load on acorn

* remove unneeded initialization of ltg*_max in post_fv3

* remove unneeded hydrostatic check from maximum_hourly_diagnostics

* fix out-of-bounds access in upp calslr_roebbr

* fix another out-of-bounds access and move a message to stdout

.gitignore

@@ -3,7 +3,7 @@
 *.pyc
 
 *.sw[a-p]
-~
+*~
 
 build/
 install/

ccpp/data/GFS_typedefs.F90

@@ -162,6 +162,7 @@ module GFS_typedefs
     real (kind=kind_phys), pointer :: pgr  (:)     => null()  !< surface pressure (Pa) real
     real (kind=kind_phys), pointer :: ugrs (:,:)   => null()  !< u component of layer wind
     real (kind=kind_phys), pointer :: vgrs (:,:)   => null()  !< v component of layer wind
+    real (kind=kind_phys), pointer :: wgrs (:,:)   => null()  !< w component of layer wind
     real (kind=kind_phys), pointer :: vvl  (:,:)   => null()  !< layer mean vertical velocity in pa/sec
     real (kind=kind_phys), pointer :: tgrs (:,:)   => null()  !< model layer mean temperature in k
     real (kind=kind_phys), pointer :: qgrs (:,:,:) => null()  !< layer mean tracer concentration
@@ -1459,6 +1460,9 @@ module GFS_typedefs
     real(kind=kind_phys) :: rhcmax          ! maximum critical relative humidity, replaces rhc_max in physcons.F90
     real(kind=kind_phys) :: huge            !< huge fill value
 
+!--- lightning threat and diagsnostics
+    logical              :: lightning_threat !< report lightning threat indices
+
     contains
       procedure :: init            => control_initialize
       procedure :: init_chemistry  => control_chemistry_initialize
@@ -1979,6 +1983,11 @@ module GFS_typedefs
     real (kind=kind_phys), pointer :: aux2d(:,:)  => null()    !< auxiliary 2d arrays in output (for debugging)
     real (kind=kind_phys), pointer :: aux3d(:,:,:)=> null()    !< auxiliary 2d arrays in output (for debugging)
 
+    !--- Lightning threat indices
+    real (kind=kind_phys), pointer :: ltg1_max(:)        => null()  !
+    real (kind=kind_phys), pointer :: ltg2_max(:)        => null()  !
+    real (kind=kind_phys), pointer :: ltg3_max(:)        => null()  !
+
     contains
       procedure :: create    => diag_create
       procedure :: rad_zero  => diag_rad_zero
@@ -2057,13 +2066,20 @@ subroutine statein_create (Statein, IM, Model)
     allocate (Statein%pgr    (IM))
     allocate (Statein%ugrs   (IM,Model%levs))
     allocate (Statein%vgrs   (IM,Model%levs))
+    if(Model%lightning_threat) then
+      allocate (Statein%wgrs   (IM,Model%levs))
+    endif
     allocate (Statein%qgrs   (IM,Model%levs,Model%ntrac))
 
     Statein%qgrs   = clear_val
     Statein%pgr    = clear_val
     Statein%ugrs   = clear_val
     Statein%vgrs   = clear_val
 
+    if(Model%lightning_threat) then
+      Statein%wgrs = clear_val
+    endif
+
     !--- soil state variables - for soil SPPT - sfc-perts, mgehne
     allocate (Statein%smc  (IM,Model%lsoil))
     allocate (Statein%stc  (IM,Model%lsoil))
@@ -3539,6 +3555,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     logical :: mix_chem = .false.         ! tracer mixing option by MYNN PBL
     logical :: fire_turb = .false.        ! enh vertmix option by MYNN PBL
 
+!-- Lightning threat index
+    logical :: lightning_threat = .false.
+
 !--- aerosol scavenging factors
     integer, parameter :: max_scav_factors = 183
     character(len=40)  :: fscav_aero(max_scav_factors)
@@ -3683,7 +3702,9 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
                                addsmoke_flag, fire_turb, mix_chem,                          &
                           !--- (DFI) time ranges with radar-prescribed microphysics tendencies
                           !          and (maybe) convection suppression
-                               fh_dfi_radar, radar_tten_limits, do_cap_suppress
+                               fh_dfi_radar, radar_tten_limits, do_cap_suppress,            &
+                          !--- GSL lightning threat indices
+                               lightning_threat
 
 !--- other parameters
     integer :: nctp    =  0                !< number of cloud types in CS scheme
@@ -3757,6 +3778,8 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     Model%flag_for_scnv_generic_tend = .true.
     Model%flag_for_dcnv_generic_tend = .true.
 
+    Model%lightning_threat = lightning_threat
+
     Model%fh_dfi_radar     = fh_dfi_radar
     Model%num_dfi_radar    = 0
     Model%dfi_radar_max_intervals = dfi_radar_max_intervals ! module-level parameter, top of file
@@ -5089,6 +5112,13 @@ subroutine control_initialize (Model, nlunit, fn_nml, me, master, &
     Model%restart          = restart
     Model%lsm_cold_start   = .not. restart
     Model%hydrostatic      = hydrostatic
+
+    if(Model%hydrostatic .and. Model%lightning_threat) then
+      write(0,*) 'Turning off lightning threat index for hydrostatic run.'
+      Model%lightning_threat = .false.
+      lightning_threat = .false.
+    endif
+
     Model%jdat(1:8)        = jdat(1:8)
     allocate(Model%si(Model%levs+1))
     !--- Define sigma level for radiation initialization
@@ -5862,6 +5892,7 @@ subroutine control_print(Model)
       print *, ' thermodyn_id      : ', Model%thermodyn_id
       print *, ' sfcpress_id       : ', Model%sfcpress_id
       print *, ' gen_coord_hybrid  : ', Model%gen_coord_hybrid
+      print *, ' hydrostatic       : ', Model%hydrostatic
       print *, ' '
       print *, 'grid extent parameters'
       print *, ' isc               : ', Model%isc
@@ -6367,7 +6398,9 @@ subroutine control_print(Model)
       print *, ' first_time_step   : ', Model%first_time_step
       print *, ' restart           : ', Model%restart
       print *, ' lsm_cold_start    : ', Model%lsm_cold_start
-      print *, ' hydrostatic       : ', Model%hydrostatic
+      print *, ' '
+      print *, 'lightning threat indexes'
+      print *, ' lightning_threat  : ', Model%lightning_threat
     endif
 
   end subroutine control_print
@@ -6961,6 +6994,15 @@ subroutine diag_create (Diag, IM, Model)
       Diag%old_pgr = clear_val
     endif
 
+    if(Model%lightning_threat) then
+       allocate (Diag%ltg1_max(IM))
+       allocate (Diag%ltg2_max(IM))
+       allocate (Diag%ltg3_max(IM))
+       Diag%ltg1_max = zero
+       Diag%ltg2_max = zero
+       Diag%ltg3_max = zero
+    endif
+
     !--- Radiation
     allocate (Diag%fluxr   (IM,Model%nfxr))
     allocate (Diag%topfsw  (IM))
@@ -7525,6 +7567,13 @@ subroutine diag_phys_zero (Diag, Model, linit, iauwindow_center)
       Diag%tsnowp  = zero
     endif
 
+! GSL lightning threat indexes
+    if(Model%lightning_threat) then
+       Diag%ltg1_max = zero
+       Diag%ltg2_max = zero
+       Diag%ltg3_max = zero
+    endif
+
   end subroutine diag_phys_zero
 
 end module GFS_typedefs

ccpp/data/GFS_typedefs.meta

@@ -111,6 +111,14 @@
   dimensions = (horizontal_loop_extent)
   type = real
   kind = kind_phys
+[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
+  active = (do_lightning_threat_index_calculations)
 [vvl]
   standard_name = lagrangian_tendency_of_air_pressure
   long_name = layer mean vertical velocity
@@ -6083,6 +6091,12 @@
   units = flag
   dimensions = ()
   type = logical
+[lightning_threat]
+  standard_name = do_lightning_threat_index_calculations
+  long_name = enables the lightning threat index calculations
+  units = flag
+  dimensions = ()
+  type = logical
 [ipt]
   standard_name = index_of_horizontal_gridpoint_for_debug_output
   long_name = horizontal index for point used for diagnostic printout
@@ -8877,6 +8891,33 @@
   dimensions = (horizontal_loop_extent)
   type = real
   kind = kind_phys
+[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
+  active = (do_lightning_threat_index_calculations)
+[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
+  active = (do_lightning_threat_index_calculations)
+[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
+  active = (do_lightning_threat_index_calculations)
 
 ########################################################################
 [ccpp-table-properties]

ccpp/driver/GFS_diagnostics.F90

@@ -4247,6 +4247,41 @@ subroutine GFS_externaldiag_populate (ExtDiag, Model, Statein, Stateout, Sfcprop
        enddo
     enddo
 
+    if(Model%lightning_threat) then
+      idx = idx + 1
+      ExtDiag(idx)%axes = 2
+      ExtDiag(idx)%name = 'ltg1_max'
+      ExtDiag(idx)%desc = 'Max Lightning Threat 1'
+      ExtDiag(idx)%unit = 'flashes/(5 min)'
+      ExtDiag(idx)%mod_name = 'gfs_sfc'
+      allocate (ExtDiag(idx)%data(nblks))
+      do nb = 1,nblks
+        ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%ltg1_max
+      enddo
+
+      idx = idx + 1
+      ExtDiag(idx)%axes = 2
+      ExtDiag(idx)%name = 'ltg2_max'
+      ExtDiag(idx)%desc = 'Max Lightning Threat 2'
+      ExtDiag(idx)%unit = 'flashes/(5 min)'
+      ExtDiag(idx)%mod_name = 'gfs_sfc'
+      allocate (ExtDiag(idx)%data(nblks))
+      do nb = 1,nblks
+        ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%ltg2_max
+      enddo
+
+      idx = idx + 1
+      ExtDiag(idx)%axes = 2
+      ExtDiag(idx)%name = 'ltg3_max'
+      ExtDiag(idx)%desc = 'Max Lightning Threat 3'
+      ExtDiag(idx)%unit = 'flashes/(5 min)'
+      ExtDiag(idx)%mod_name = 'gfs_sfc'
+      allocate (ExtDiag(idx)%data(nblks))
+      do nb = 1,nblks
+        ExtDiag(idx)%data(nb)%var2 => IntDiag(nb)%ltg3_max
+      enddo
+    endif
+
     ! Cloud effective radii from Microphysics
     if (Model%imp_physics == Model%imp_physics_thompson .or. Model%imp_physics == Model%imp_physics_fer_hires .or. &
         Model%imp_physics == Model%imp_physics_nssl  ) then

io/post_fv3.F90

@@ -496,7 +496,8 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
       use vrbls3d,     only: t, q, uh, vh, wh, alpint, dpres, zint, zmid, o3,  &
                              qqr, qqs, cwm, qqi, qqw, qqg, omga, cfr, pmid,    &
                              q2, rlwtt, rswtt, tcucn, tcucns, train, el_pbl,   &
-                             pint, exch_h, ref_10cm, qqni, qqnr, qqnwfa, qqnifa
+                             pint, exch_h, ref_10cm, qqni, qqnr, qqnwfa,       &
+                             qqnifa
       use vrbls2d,     only: f, pd, sigt4, fis, pblh, ustar, z0, ths, qs, twbs,&
                              qwbs, avgcprate, cprate, avgprec, prec, lspa, sno,&
                              cldefi, th10, q10, tshltr, pshltr, albase,        &
@@ -528,7 +529,8 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
                              refdm10c_max, u10max, v10max, wspd10max, sfcuxi,  &
                              sfcvxi, t10m, t10avg, psfcavg, akhsavg, akmsavg,  &
                              albedo, tg, prate_max, pwat, snow_acm, snow_bkt,  &
-                             acgraup, graup_bucket, acfrain, frzrn_bucket 
+                             acgraup, graup_bucket, acfrain, frzrn_bucket,     &
+                             ltg1_max, ltg2_max, ltg3_max
       use soil,        only: sldpth, sh2o, smc, stc
       use masks,       only: lmv, lmh, htm, vtm, gdlat, gdlon, dx, dy, hbm2, sm, sice
       use ctlblk_mod,  only: im, jm, lm, lp1, jsta, jend, jsta_2l, jend_2u, jsta_m,jend_m, &
@@ -732,7 +734,8 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
 !                     snow phase change heat flux, snopcx
 ! GFS does not use total momentum flux,sfcuvx
 !$omp parallel do default(none),private(i,j),shared(jsta,jend,im,spval,ista,iend), &
-!$omp& shared(acfrcv,ncfrcv,acfrst,ncfrst,bgroff,rswin,rswinc,rswout,snopcx,sfcuvx)
+!$omp& shared(acfrcv,ncfrcv,acfrst,ncfrst,bgroff,rswin,rswinc,rswout,snopcx,sfcuvx,&
+!$omp& ltg1_max,ltg2_max,ltg3_max)
       do j=jsta,jend
         do i=ista,iend
           acfrcv(i,j) = spval
@@ -1038,6 +1041,39 @@ subroutine set_postvars_fv3(wrt_int_state,grid_id,mype,mpicomp)
               enddo
             endif
 
+            ! Lightning threat index 1
+            if(trim(fieldname)=='ltg1_max') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,ltg1_max,arrayr42d,fillValue,spval)
+              do j=jsta,jend
+                do i=ista, iend
+                  ltg1_max(i,j)=arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillValue) < small) ltg1_max(i,j)=spval
+                enddo
+              enddo
+            endif
+
+            ! Lightning threat index 2
+            if(trim(fieldname)=='ltg2_max') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,ltg2_max,arrayr42d,fillValue,spval)
+              do j=jsta,jend
+                do i=ista, iend
+                  ltg2_max(i,j)=arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillValue) < small) ltg2_max(i,j)=spval
+                enddo
+              enddo
+            endif
+
+            ! Lightning threat index 3
+            if(trim(fieldname)=='ltg3_max') then
+              !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,ltg3_max,arrayr42d,fillValue,spval)
+              do j=jsta,jend
+                do i=ista, iend
+                  ltg3_max(i,j)=arrayr42d(i,j)
+                  if(abs(arrayr42d(i,j)-fillValue) < small) ltg3_max(i,j)=spval
+                enddo
+              enddo
+            endif
+
             ! frictional velocity
             if(trim(fieldname)=='fricv') then
               !$omp parallel do default(none) private(i,j) shared(jsta,jend,ista,iend,ustar,arrayr42d,fillValue,spval)

io/post_nems_routines.F90

@@ -226,6 +226,7 @@ subroutine read_postnmlt(kpo,kth,kpv,po,th,pv,post_namelist)
                              lsmdef,ALSL,me,d3d_on,gocart_on,hyb_sigp,&
                              pthresh,novegtype,ivegsrc,icu_physics,   &
                              isf_surface_physics,modelname,submodelname
+      use upp_ifi_mod, only: write_ifi_debug_files
 !
 !    revision history:
 !    Jul 2019 Jun Wang: read post namelist
@@ -241,7 +242,7 @@ subroutine read_postnmlt(kpo,kth,kpv,po,th,pv,post_namelist)
       integer l,k
 
       namelist/nampgb/kpo,po,kth,th,kpv,pv,popascal,d3d_on,gocart_on,  &
-                      hyb_sigp
+                      hyb_sigp,write_ifi_debug_files
       namelist/model_inputs/modelname,submodelname
 !---------------------------------------------------------------------
 !