From 2353eaa7ed753db84c8e3fbaad6927cd54c2f893 Mon Sep 17 00:00:00 2001 From: daviddowellNOAA <72174157+daviddowellNOAA@users.noreply.github.com> Date: Thu, 7 Dec 2023 10:02:04 -0600 Subject: [PATCH] Add assimilation of GLM flash-extent density (FED) observations to EnKF (#654) **Description** Fixes https://github.com/NOAA-EMC/GSI/issues/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 | 2 +- src/enkf/enkf_obs_sensitivity.f90 | 3 +++ src/enkf/enkf_obsmod.f90 | 7 +++++++ src/enkf/innovstats.f90 | 16 ++++++++++++++++ src/enkf/params.f90 | 15 +++++++++++++++ src/enkf/readconvobs.f90 | 12 +++++++++--- 6 files changed, 51 insertions(+), 4 deletions(-) diff --git a/src/enkf/enkf.f90 b/src/enkf/enkf.f90 index d35613b585..479f60c019 100644 --- a/src/enkf/enkf.f90 +++ b/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) diff --git a/src/enkf/enkf_obs_sensitivity.f90 b/src/enkf/enkf_obs_sensitivity.f90 index 6c37936f31..72296d5934 100644 --- a/src/enkf/enkf_obs_sensitivity.f90 +++ b/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 diff --git a/src/enkf/enkf_obsmod.f90 b/src/enkf/enkf_obsmod.f90 index ea8f6446fb..eb4f9c8e58 100644 --- a/src/enkf/enkf_obsmod.f90 +++ b/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 diff --git a/src/enkf/innovstats.f90 b/src/enkf/innovstats.f90 index 68668218fc..853532c9b9 100644 --- a/src/enkf/innovstats.f90 +++ b/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) diff --git a/src/enkf/params.f90 b/src/enkf/params.f90 index 62701c24a7..36b0c9c207 100644 --- a/src/enkf/params.f90 +++ b/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 diff --git a/src/enkf/readconvobs.f90 b/src/enkf/readconvobs.f90 index d1f4ec3ff8..a5383069a1 100644 --- a/src/enkf/readconvobs.f90 +++ b/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.