Add assimilation of GLM flash-extent density (FED) observations to En…

…KF (NOAA-EMC#654)

**Description**

Fixes NOAA-EMC#653

The proposed code changes
(1) add a new "fed" observation type to the EnKF
(2) add localization parameters, with namelist control, for FED
observations
(3) output basic statistics for FED observations.

In the RRFS, the FED observations will be assimilated together with
radar-reflectivity observations. The localization parameters for the
reflectivity observations in RRFS are corrlength=18 and lnsigcutoff=0.5.
However, these tight localization distances don't work well for the
sparse FED observations. Therefore, localization parameters for FED
observations, with namelist control, were added to allow the FED
observations to influence the model state over longer distances. The
default localization parameters for FED observations (corrlength=30 and
lnsigcutoff=2.0) were determined through experimentation with WRF and
FV3 convection-allowing (3-km horizontal grid spacing) ensembles.

**Type of change**

- [X] New feature (non-breaking change which adds functionality)

**How Has This Been Tested?**

Hourly cycling with simultaneous EnKF assimilation of FED and
radar-reflectivity observations has been tested for a CONUS version of
the prototype RRFSv1 for two short (1-2 days) retrospective periods, one
in July 2022 and the other in August 2023. The impacts of the FED
observations on the analyses are greatest over the oceans far from land,
where the radar network does not provide observations.

src/enkf/enkf.f90

@@ -51,7 +51,7 @@ module enkf
 !  NH, tropics and SH, and in the horizontal, vertical and time dimensions,
 !  using the namelist parameters  corrlengthnh, corrlengthtr, corrlengthsh,
 !  lnsigcutoffnh, lnsigcutofftr, lnsigcutoffsh (lnsigcutoffsatnh,
-!  lnsigcutoffsattr, lnsigcutoffsatsh for satellite obs, similar for ps obs)
+!  lnsigcutoffsattr, lnsigcutoffsatsh for satellite obs, similar for ps and fed obs)
 !  obtimelnh, obtimeltr, obtimelsh. The length scales should be given in km for the
 !  horizontal, hours for time, and 'scale heights' (units of -log(p/pref)) in the
 !  vertical. The function used for localization (function taper)

src/enkf/enkf_obs_sensitivity.f90

@@ -36,6 +36,7 @@ module enkf_obs_sensitivity
 use params, only: efsoi_flag,latbound,nlevs,nanals,datestring, &
                   lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh, &
                   lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh, &
+                  lnsigcutofffednh,lnsigcutofffedtr,lnsigcutofffedsh, &
                   lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh, &
                   corrlengthnh,corrlengthtr,corrlengthsh, &
                   obtimelnh,obtimeltr,obtimelsh,letkf_flag, &
@@ -292,6 +293,8 @@ subroutine read_ob_sens
         lnsigl(nob) = latval(deglat,lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh)
      else if (obtype(nob)(1:3) == ' ps') then
         lnsigl(nob) = latval(deglat,lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh)
+     else if (obtype(nob)(1:3) == 'fed') then
+        lnsigl(nob) = latval(deglat,lnsigcutofffednh,lnsigcutofffedtr,lnsigcutofffedsh)
      else
         lnsigl(nob)=latval(deglat,lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh)
      end if

src/enkf/enkf_obsmod.f90

@@ -109,6 +109,8 @@ module enkf_obsmod
       lnsigcutoffnh, lnsigcutoffsh, lnsigcutofftr, corrlengthnh,&
       corrlengthtr, corrlengthsh, obtimelnh, obtimeltr, obtimelsh,&
       lnsigcutoffsatnh, lnsigcutoffsatsh, lnsigcutoffsattr,&
+      lnsigcutofffednh, lnsigcutofffedsh, lnsigcutofffedtr,&
+      corrlengthfednh, corrlengthfedtr, corrlengthfedsh,   &
       varqc, huber, zhuberleft, zhuberright, modelspace_vloc, &
       lnsigcutoffpsnh, lnsigcutoffpssh, lnsigcutoffpstr, neigv, &
       lnsigcutoffrdrnh, lnsigcutoffrdrsh, lnsigcutoffrdrtr,&
@@ -276,6 +278,8 @@ subroutine readobs()
       lnsigl(nob) = latval(deglat,lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh)
    else if (obtype(nob)(1:3) == ' ps') then
       lnsigl(nob) = latval(deglat,lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh)
+   else if (obtype(nob)(1:3) == 'fed') then
+      lnsigl(nob) = latval(deglat,lnsigcutofffednh,lnsigcutofffedtr,lnsigcutofffedsh)
    else if ( (obtype(nob)(1:3) == 'dbz' .or. obtype(nob)(1:3) == ' rw') .and. l_use_enkf_directZDA ) then
       lnsigl(nob) = latval(deglat,lnsigcutoffrdrnh,lnsigcutoffrdrtr,lnsigcutoffrdrsh)
    else
@@ -293,6 +297,9 @@ subroutine readobs()
    if ( (obtype(nob)(1:3) == 'dbz' .or. obtype(nob)(1:3) == ' rw') .and. l_use_enkf_directZDA ) then
        corrlengthsq(nob)=latval(deglat,corrlengthrdrnh,corrlengthrdrtr,corrlengthrdrsh)**2
    end if
+   if (obtype(nob)(1:3) == 'fed') then
+       corrlengthsq(nob)=latval(deglat,corrlengthfednh,corrlengthfedtr,corrlengthfedsh)**2
+   end if
    obtimel(nob)=latval(deglat,obtimelnh,obtimeltr,obtimelsh)
 end do
 

src/enkf/innovstats.f90

@@ -45,6 +45,7 @@ subroutine print_innovstats(obfit,obsprd)
  nobsspd_nh,nobsspd_sh,nobsspd_tr,&
  nobsgps_nh,nobsgps_sh,nobsgps_tr,&
  nobsdbz_nh,nobsdbz_sh,nobsdbz_tr,&
+ nobsfed_nh,nobsfed_sh,nobsfed_tr,&
  nobsrw_nh,nobsrw_sh,nobsrw_tr,&
  nobsq_nh,nobsq_sh,nobsq_tr,nobswnd_nh,nobswnd_sh,nobswnd_tr,&
  nobsoz_nh,nobsoz_sh,nobsoz_tr,nobsps_sh,nobsps_nh,nobsps_tr,nob
@@ -67,6 +68,9 @@ subroutine print_innovstats(obfit,obsprd)
  sumdbz_nh,biasdbz_nh,sumdbz_spread_nh,sumdbz_oberr_nh,&
  sumdbz_sh,biasdbz_sh,sumdbz_spread_sh,sumdbz_oberr_sh,&
  sumdbz_tr,biasdbz_tr,sumdbz_spread_tr,sumdbz_oberr_tr,&
+ sumfed_nh,biasfed_nh,sumfed_spread_nh,sumfed_oberr_nh,&
+ sumfed_sh,biasfed_sh,sumfed_spread_sh,sumfed_oberr_sh,&
+ sumfed_tr,biasfed_tr,sumfed_spread_tr,sumfed_oberr_tr,&
  sumrw_nh,biasrw_nh,sumrw_spread_nh,sumrw_oberr_nh,&
  sumrw_sh,biasrw_sh,sumrw_spread_sh,sumrw_oberr_sh,&
  sumrw_tr,biasrw_tr,sumrw_spread_tr,sumrw_oberr_tr,&
@@ -112,6 +116,9 @@ subroutine print_innovstats(obfit,obsprd)
   nobsdbz_nh = 0
   nobsdbz_sh = 0
   nobsdbz_tr = 0
+  nobsfed_nh = 0
+  nobsfed_sh = 0
+  nobsfed_tr = 0
   nobsrw_nh = 0
   nobsrw_sh = 0
   nobsrw_tr = 0
@@ -168,6 +175,12 @@ subroutine print_innovstats(obfit,obsprd)
                  sumdbz_nh,biasdbz_nh,sumdbz_spread_nh,sumdbz_oberr_nh,nobsdbz_nh,&
                  sumdbz_sh,biasdbz_sh,sumdbz_spread_sh,sumdbz_oberr_sh,nobsdbz_sh,&
                  sumdbz_tr,biasdbz_tr,sumdbz_spread_tr,sumdbz_oberr_tr,nobsdbz_tr)
+         else if (obtype(nob)(1:3) == 'fed') then
+            call obstats(obfit(nob),oberrvar_orig(nob),&
+                 obsprd(nob),obloclat(nob),&
+                 sumfed_nh,biasfed_nh,sumfed_spread_nh,sumfed_oberr_nh,nobsfed_nh,&
+                 sumfed_sh,biasfed_sh,sumfed_spread_sh,sumfed_oberr_sh,nobsfed_sh,&
+                 sumfed_tr,biasfed_tr,sumfed_spread_tr,sumfed_oberr_tr,nobsfed_tr)
          else if (obtype(nob)(1:3) == ' rw') then
             call obstats(obfit(nob),oberrvar_orig(nob),&
                  obsprd(nob),obloclat(nob),&
@@ -216,6 +229,9 @@ subroutine print_innovstats(obfit,obsprd)
    call printstats('  all dbz',sumdbz_nh,biasdbz_nh,sumdbz_spread_nh,sumdbz_oberr_nh,nobsdbz_nh,&
         sumdbz_sh,biasdbz_sh,sumdbz_spread_sh,sumdbz_oberr_sh,nobsdbz_sh,&
         sumdbz_tr,biasdbz_tr,sumdbz_spread_tr,sumdbz_oberr_tr,nobsdbz_tr)
+   call printstats('  all fed',sumfed_nh,biasfed_nh,sumfed_spread_nh,sumfed_oberr_nh,nobsfed_nh,&
+        sumfed_sh,biasfed_sh,sumfed_spread_sh,sumfed_oberr_sh,nobsfed_sh,&
+        sumfed_tr,biasfed_tr,sumfed_spread_tr,sumfed_oberr_tr,nobsfed_tr)
    call printstats('  all  rw',sumrw_nh,biasq_nh,sumrw_spread_nh,sumrw_oberr_nh,nobsrw_nh,&
         sumrw_sh,biasrw_sh,sumrw_spread_sh,sumrw_oberr_sh,nobsrw_sh,&
         sumrw_tr,biasrw_tr,sumrw_spread_tr,sumrw_oberr_tr,nobsrw_tr)

src/enkf/params.f90

@@ -124,6 +124,8 @@ module params
 real(r_single),public ::  lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh,&
                lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh,&
                lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh
+real(r_single),public ::  corrlengthfednh,corrlengthfedtr,corrlengthfedsh, &
+               lnsigcutofffednh,lnsigcutofffedtr,lnsigcutofffedsh
 real(r_single),public ::  corrlengthrdrnh,corrlengthrdrtr,corrlengthrdrsh, &
                lnsigcutoffrdrnh,lnsigcutoffrdrtr,lnsigcutoffrdrsh
 real(r_single),public :: analpertwtnh,analpertwtsh,analpertwttr,sprd_tol,saterrfact
@@ -261,6 +263,8 @@ module params
                    lnsigcutoffnh,lnsigcutofftr,lnsigcutoffsh,&
                    lnsigcutoffsatnh,lnsigcutoffsattr,lnsigcutoffsatsh,&
                    lnsigcutoffpsnh,lnsigcutoffpstr,lnsigcutoffpssh,&
+                   corrlengthfednh,corrlengthfedsh,corrlengthfedtr,&
+                   lnsigcutofffednh,lnsigcutofffedsh,lnsigcutofffedtr,&
                    fgfileprefixes,fgsfcfileprefixes,anlfileprefixes, &
                    anlsfcfileprefixes,incfileprefixes,incsfcfileprefixes,&
                    statefileprefixes,statesfcfileprefixes, &
@@ -317,6 +321,10 @@ subroutine read_namelist()
 corrlengthrdrnh = 10
 corrlengthrdrtr = 10
 corrlengthrdrsh = 10
+! corrlength (km) for GLM flash extent density
+corrlengthfednh = 30_r_single
+corrlengthfedtr = 30_r_single
+corrlengthfedsh = 30_r_single
 ! read in localization length scales from an external file.
 readin_localization = .false.
 ! min and max inflation.
@@ -341,6 +349,9 @@ subroutine read_namelist()
 lnsigcutoffrdrnh = 0.2_r_single  ! value for radar
 lnsigcutoffrdrtr = 0.2_r_single  ! value for radar
 lnsigcutoffrdrsh = 0.2_r_single  ! value for radar
+lnsigcutofffednh = 2._r_single    ! value for GLM flash extent density
+lnsigcutofffedtr = 2._r_single    ! value for GLM flash extent density 
+lnsigcutofffedsh = 2._r_single    ! value for GLM flash extent density
 ! ob time localization
 obtimelnh = 1.e10_r_single
 obtimeltr = 1.e10_r_single
@@ -813,6 +824,10 @@ subroutine read_namelist()
 corrlengthrdrnh = corrlengthrdrnh * 1.e3_r_single/rearth
 corrlengthrdrtr = corrlengthrdrtr * 1.e3_r_single/rearth
 corrlengthrdrsh = corrlengthrdrsh * 1.e3_r_single/rearth
+! rescale covariance localization length for GLM FED
+corrlengthfednh = corrlengthfednh * 1.e3_r_single/rearth
+corrlengthfedtr = corrlengthfedtr * 1.e3_r_single/rearth
+corrlengthfedsh = corrlengthfedsh * 1.e3_r_single/rearth
 
 ! convert targe area boundary into radians
 tar_minlat = tar_minlat * deg2rad

src/enkf/readconvobs.f90

@@ -42,9 +42,9 @@ module readconvobs
 
 
 !> observation types to read from netcdf files
-integer(i_kind), parameter :: nobtype = 11
+integer(i_kind), parameter :: nobtype = 12
 character(len=3), dimension(nobtype), parameter :: obtypes = (/'  t', '  q', ' ps', ' uv', 'tcp', &
-                                                               'gps', 'spd', ' pw', ' dw', ' rw', 'dbz' /)
+                                                               'gps', 'spd', ' pw', ' dw', ' rw', 'dbz', 'fed' /)
 
 contains
 
@@ -79,7 +79,7 @@ subroutine get_num_convobs_bin(obspath,datestring,num_obs_tot,num_obs_totdiag,id
     integer(i_kind)  :: iunit, nchar, nreal, ii, mype, ios, idate, i, ipe, ioff0
     integer(i_kind),dimension(2) :: nn,nobst, nobsps, nobsq, nobsuv, nobsgps, &
          nobstcp,nobstcx,nobstcy,nobstcz,nobssst, nobsspd, nobsdw, nobsrw, nobspw, &
-         nobsdbz
+         nobsdbz, nobsfed
     character(8),allocatable,dimension(:):: cdiagbuf
     real(r_single),allocatable,dimension(:,:)::rdiagbuf
     real(r_kind) :: errorlimit,errorlimit2,error,pres,obmax
@@ -104,6 +104,7 @@ subroutine get_num_convobs_bin(obspath,datestring,num_obs_tot,num_obs_totdiag,id
     nobspw = 0
     nobsgps = 0
     nobsdbz = 0
+    nobsfed = 0
     nobstcp = 0; nobstcx = 0; nobstcy = 0; nobstcz = 0
     init_pass = .true.
     peloop: do ipe=0,npefiles
@@ -187,6 +188,9 @@ subroutine get_num_convobs_bin(obspath,datestring,num_obs_tot,num_obs_totdiag,id
        else if (obtype == 'dbz') then
           nobsdbz = nobsdbz + nn
           num_obs_tot = num_obs_tot + nn(2)
+       else if (obtype == 'fed') then
+          nobsfed = nobsfed + nn
+          num_obs_tot = num_obs_tot + nn(2)
        else if (obtype == 'gps') then
           nobsgps = nobsgps + nn
           num_obs_tot = num_obs_tot + nn(2)
@@ -231,6 +235,7 @@ subroutine get_num_convobs_bin(obspath,datestring,num_obs_tot,num_obs_totdiag,id
           write(6,100) 'dw',nobsdw(1),nobsdw(2)
           write(6,100) 'rw',nobsrw(1),nobsrw(2)
           write(6,100) 'dbz',nobsdbz(1),nobsdbz(2)
+          write(6,100) 'fed',nobsfed(1),nobsfed(2)
           write(6,100) 'tcp',nobstcp(1),nobstcp(2)
           if (nobstcx(2) .gt. 0) then
              write(6,100) 'tcx',nobstcx(1),nobstcx(2)
@@ -1075,6 +1080,7 @@ subroutine get_convobs_data_bin(obspath, datestring, nobs_max, nobs_maxdiag,   &
     if (obtype == '  t' .or. obtype == ' uv' .or. obtype == ' ps' .or. &
         obtype == 'tcp' .or. obtype == '  q' .or. obtype == 'spd' .or. &
         obtype == 'sst' .or. obtype == ' rw' .or. obtype == 'dbz' .or. &
+        obtype == 'fed' .or.                                           &
         obtype == 'gps' .or. obtype == ' dw' .or. obtype == ' pw')  then
 
 !   direct reflectivitiy DA has a different routine for dbz obs.