Add FED EnVar DA Capability

src/gsi/gsi_fedOper.F90

@@ -9,6 +9,7 @@ module gsi_fedOper
 !   2023-07-10  D. Dowell - created new module for FED (flash extent
 !                           density); gsi_dbzOper.F90 code used as a
 !                           starting point for developing this new module
+!   2023-08-24  H. Wang   - Turned on intfed and stpfed
 !
 !   input argument list: see Fortran 90 style document below
 !
@@ -128,6 +129,7 @@ subroutine setup_(self, lunin, mype, is, nobs, init_pass, last_pass)
   end subroutine setup_
 
   subroutine intjo1_(self, ibin, rval,sval, qpred,sbias)
+    use intfedmod, only: intjo => intfed
     use gsi_bundlemod  , only: gsi_bundle
     use bias_predictors, only: predictors
     use m_obsNode , only: obsNode
@@ -145,9 +147,14 @@ subroutine intjo1_(self, ibin, rval,sval, qpred,sbias)
     character(len=*),parameter:: myname_=myname//"::intjo1_"
     class(obsNode),pointer:: headNode
 
+    headNode => obsLList_headNode(self%obsLL(ibin))
+    call intjo(headNode, rval,sval)
+    headNode => null()
+
   end subroutine intjo1_
 
   subroutine stpjo1_(self, ibin, dval,xval,pbcjo,sges,nstep,dbias,xbias)
+    use stpfedmod, only: stpjo => stpfed
     use gsi_bundlemod, only: gsi_bundle
     use bias_predictors, only: predictors
     use m_obsNode , only: obsNode
@@ -169,6 +176,9 @@ subroutine stpjo1_(self, ibin, dval,xval,pbcjo,sges,nstep,dbias,xbias)
     character(len=*),parameter:: myname_=myname//"::stpjo1_"
     class(obsNode),pointer:: headNode
 
+    headNode => obsLList_headNode(self%obsLL(ibin))
+    call stpjo(headNode,dval,xval,pbcjo(:),sges,nstep)
+    headNode => null()
   end subroutine stpjo1_
 
 end module gsi_fedOper

src/gsi/gsi_files.cmake

@@ -274,6 +274,7 @@ intaod.f90
 intcldch.f90
 intco.f90
 intdbz.f90
+intfed.f90
 intdw.f90
 intgps.f90
 intgust.f90
@@ -594,6 +595,7 @@ stpcalc.f90
 stpcldch.f90
 stpco.f90
 stpdbz.f90
+stpfed.f90
 stpdw.f90
 stpgps.f90
 stpgust.f90

src/gsi/gsimod.F90

@@ -23,12 +23,13 @@ module gsimod
   use gsi_dbzOper, only: diag_radardbz
   use gsi_fedOper, only: diag_fed
 
-  use obsmod, only: doradaroneob,oneoblat,oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,&
+  use obsmod, only: doradaroneob,dofedoneob,oneoblat,oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,&
      radar_no_thinning,ens_hx_dbz_cut,static_gsi_nopcp_dbz,rmesh_dbz,&
-     rmesh_vr,zmesh_dbz,zmesh_vr,if_vterminal, if_model_dbz,if_vrobs_raw,if_use_w_vr,&
+     rmesh_vr,zmesh_dbz,zmesh_vr,if_vterminal, if_model_dbz,if_model_fed,innov_use_model_fed,if_vrobs_raw,if_use_w_vr,&
      minobrangedbz,maxobrangedbz,maxobrangevr,maxtiltvr,missing_to_nopcp,&
      ntilt_radarfiles,whichradar,&
-     minobrangevr,maxtiltdbz,mintiltvr,mintiltdbz,l2rwthin,hurricane_radar 
+     minobrangevr,maxtiltdbz,mintiltvr,mintiltdbz,l2rwthin,hurricane_radar,&
+     r_hgt_fed
 
   use obsmod, only: lwrite_predterms, &
      lwrite_peakwt,use_limit,lrun_subdirs,l_foreaft_thin,lobsdiag_forenkf,&
@@ -202,7 +203,7 @@ module gsimod
   use gsi_nstcouplermod, only: gsi_nstcoupler_init_nml
   use gsi_nstcouplermod, only: nst_gsi,nstinfo,zsea1,zsea2,fac_dtl,fac_tsl
   use ncepnems_io, only: init_nems,imp_physics,lupp
-  use wrf_vars_mod, only: init_wrf_vars
+  use wrf_vars_mod, only: init_wrf_vars,fed_exist
   use gsi_rfv3io_mod,only : fv3sar_bg_opt
   use radarz_cst,            only: mphyopt, MFflg
   use radarz_iface,          only: init_mphyopt
@@ -510,6 +511,15 @@ module gsimod
 !  2023-07-30 Zhao    - added namelist options for analysis of significant wave height
 !                       (aka howv in GSI code): corp_howv, hwllp_howv
 !                       (in namelist session rapidrefresh_cldsurf)
+!  
+!  2023-09-14 H. Wang - add namelist option for FED EnVar DA. 
+!                        - if_model_fed=.true.        :  FED in background and ens. If
+!                          perform FED DA, this has to be true along with fed in
+!                          control/analysis variable list. If only run GSI observer,
+!                          it can be false.
+!                        - innov_use_model_fed=.true. :  Use FED from BG to calculate innovation.
+!                          this requires if_model_fed=.true. 
+!                          it works either an EnVar DA run or a GSI observer run.
 !
 !EOP
 !-------------------------------------------------------------------------
@@ -770,16 +780,17 @@ module gsimod
        use_sp_eqspace,lnested_loops,lsingleradob,thin4d,use_readin_anl_sfcmask,&
        luse_obsdiag,id_drifter,id_ship,verbose,print_obs_para,lsingleradar,singleradar,lnobalance, &
        missing_to_nopcp,minobrangedbz,minobrangedbz,maxobrangedbz,&
-       maxobrangevr,maxtiltvr,whichradar,doradaroneob,oneoblat,&
+       maxobrangevr,maxtiltvr,whichradar,doradaroneob,dofedoneob,oneoblat,&
        oneoblon,oneobheight,oneobvalue,oneobddiff,oneobradid,&
        rmesh_vr,zmesh_dbz,zmesh_vr, ntilt_radarfiles, whichradar,&
        radar_no_thinning,ens_hx_dbz_cut,static_gsi_nopcp_dbz,rmesh_dbz,&
        minobrangevr, maxtiltdbz, mintiltvr,mintiltdbz,if_vterminal,if_vrobs_raw,if_use_w_vr,&
-       if_model_dbz,imp_physics,lupp,netcdf_diag,binary_diag,l_wcp_cwm,aircraft_recon,diag_version,&
+       if_model_dbz,if_model_fed,innov_use_model_fed,imp_physics,lupp,netcdf_diag,binary_diag,l_wcp_cwm,aircraft_recon,diag_version,&
        write_fv3_incr,incvars_to_zero,incvars_zero_strat,incvars_efold,diag_version,&
        cao_check,lcalc_gfdl_cfrac,tau_fcst,efsoi_order,lupdqc,lqcoef,cnvw_option,l2rwthin,hurricane_radar,&
        l_reg_update_hydro_delz, l_obsprvdiag,&
-       l_use_dbz_directDA, l_use_rw_columntilt, ta2tb, optconv
+       l_use_dbz_directDA, l_use_rw_columntilt, ta2tb, optconv, &
+       r_hgt_fed
 
 ! GRIDOPTS (grid setup variables,including regional specific variables):
 !     jcap     - spectral resolution
@@ -1978,6 +1989,15 @@ subroutine gsimain_initialize
      endif
   endif
 
+  if (innov_use_model_fed .and. .not.if_model_fed) then
+     if(mype==0) write(6,*)' GSIMOD: invalid innov_use_model_fed=.true. but if_model_fed=.false.'
+     call die(myname_,'invalid innov_use_model_fed,if_model_fed, check namelist settings',332)
+  end if
+
+  if (miter > 0 .and. if_model_fed .and. .not. fed_exist) then
+     if(mype==0) write(6,*)' GSIMOD: invalid miter > 0 and if_model_fed=.true. but fed is not in anavinfo file'
+     call die(myname_,'Please add fed in anavinfo (contro/state_vector and met_guess) when miter > 0 and if_model_fed=.true.',334)
+  end if
 
 ! Ensure valid number of horizontal scales
   if (nhscrf<0 .or. nhscrf>3) then

src/gsi/intfed.f90

@@ -0,0 +1,187 @@
+module intfedmod
+!$$$ module documentation block
+!           .      .    .                                       .
+! module:   intfedmod    module for intfed and its tangent linear intfed_tl
+!   prgmmr:
+!
+! abstract: module for intfed and its tangent linear intfed_tl
+!
+! program history log:
+!   2023-08-24  H. Wang - add tangent linear of fed operator to directly assimilate FED 
+!
+! subroutines included:
+!   sub intfed_
+!
+! variable definitions:
+!
+! attributes:
+!   language: f90
+!   machine:
+!
+!$$$ end documentation block
+
+use m_obsNode, only: obsNode
+use m_fedNode, only: fedNode
+use m_fedNode, only: fedNode_typecast
+use m_fedNode, only: fedNode_nextcast
+use m_obsdiagNode, only: obsdiagNode_set
+implicit none
+
+PRIVATE
+PUBLIC intfed
+
+interface intfed; module procedure &
+          intfed_
+end interface
+
+contains
+
+subroutine intfed_(fedhead,rval,sval)
+!$$$  subprogram documentation block
+!                .      .    .                                       .
+! subprogram:    intfed       apply nonlin qc operator for GLM FED 
+!
+! abstract: apply observation operator for radar winds
+!             with nonlinear qc operator
+!
+! program history log:
+!   2023-08-24  H.Wang   - modified based on intdbz.f90 
+!                        - using tangent linear fed operator 
+
+!
+!   input argument list:
+!     fedhead   - obs type pointer to obs structure     
+!     sfed      - current fed solution increment
+!
+!   output argument list:
+!     rfed       - fed results from fed observation operator
+!
+! attributes:
+!   language: f90
+!   machine:  ibm RS/6000 SP
+!
+!$$$
+  use kinds, only: r_kind,i_kind
+  use constants, only: half,one,tiny_r_kind,cg_term,r3600
+  use obsmod, only: lsaveobsens,l_do_adjoint,luse_obsdiag
+  use qcmod, only: nlnqc_iter,varqc_iter
+  use jfunc, only: jiter
+  use gsi_bundlemod, only: gsi_bundle
+  use gsi_bundlemod, only: gsi_bundlegetpointer
+  use gsi_4dvar, only: ladtest_obs
+  use wrf_vars_mod, only : fed_exist
+  implicit none
+
+! Declare passed variables
+  class(obsNode), pointer, intent(in   ) :: fedhead
+  type(gsi_bundle),        intent(in   ) :: sval
+  type(gsi_bundle),        intent(inout) :: rval
+
+! Declare local variables
+  integer(i_kind) j1,j2,j3,j4,j5,j6,j7,j8,ier,istatus
+! real(r_kind) penalty
+  real(r_kind) val,w1,w2,w3,w4,w5,w6,w7,w8,valfed
+  real(r_kind) cg_fed,p0,grad,wnotgross,wgross,pg_fed
+  real(r_kind),pointer,dimension(:) :: sfed
+  real(r_kind),pointer,dimension(:) :: rfed
+  type(fedNode), pointer :: fedptr
+
+!  If no fed obs type data return
+  if(.not. associated(fedhead))return
+
+! Retrieve pointers
+! Simply return if any pointer not found
+  ier=0
+  if(fed_exist)then
+    call gsi_bundlegetpointer(sval,'fed',sfed,istatus);ier=istatus+ier
+    call gsi_bundlegetpointer(rval,'fed',rfed,istatus);ier=istatus+ier
+  else
+    return
+  end if
+
+  if(ier/=0)return
+
+
+  fedptr => fedNode_typecast(fedhead)
+  do while (associated(fedptr))
+     j1=fedptr%ij(1)
+     j2=fedptr%ij(2)
+     j3=fedptr%ij(3)
+     j4=fedptr%ij(4)
+     j5=fedptr%ij(5)
+     j6=fedptr%ij(6)
+     j7=fedptr%ij(7)
+     j8=fedptr%ij(8)
+     w1=fedptr%wij(1)
+     w2=fedptr%wij(2)
+     w3=fedptr%wij(3)
+     w4=fedptr%wij(4)
+     w5=fedptr%wij(5)
+     w6=fedptr%wij(6)
+     w7=fedptr%wij(7)
+     w8=fedptr%wij(8)
+
+
+!    Forward model
+     if( fed_exist )then
+       val = w1* sfed(j1)+w2* sfed(j2)+w3* sfed(j3)+w4* sfed(j4)+ &
+             w5* sfed(j5)+w6* sfed(j6)+w7* sfed(j7)+w8* sfed(j8)
+     end if
+
+     if(luse_obsdiag)then
+        if (lsaveobsens) then
+           grad = val*fedptr%raterr2*fedptr%err2
+           !-- fedptr%diags%obssen(jiter) = grad
+           call obsdiagNode_set(fedptr%diags,jiter=jiter,obssen=grad)
+
+        else
+           !-- if (fedptr%luse) fedptr%diags%tldepart(jiter)=val
+           if (fedptr%luse) call obsdiagNode_set(fedptr%diags,jiter=jiter,tldepart=val)
+        endif
+     endif
+
+     if (l_do_adjoint) then
+        if (.not. lsaveobsens) then
+           if( .not. ladtest_obs ) val=val-fedptr%res
+
+!          gradient of nonlinear operator
+           if (nlnqc_iter .and. fedptr%pg > tiny_r_kind .and. &
+                                fedptr%b  > tiny_r_kind) then
+              pg_fed=fedptr%pg*varqc_iter
+              cg_fed=cg_term/fedptr%b
+              wnotgross= one-pg_fed
+              wgross = pg_fed*cg_fed/wnotgross
+              p0   = wgross/(wgross+exp(-half*fedptr%err2*val**2))
+              val = val*(one-p0)
+           endif
+
+           if( ladtest_obs)  then
+              grad = val
+           else
+              grad = val*fedptr%raterr2*fedptr%err2
+           end if
+
+        endif
+
+!       Adjoint
+        if(fed_exist)then
+          valfed = grad
+          rfed(j1)=rfed(j1)+w1*valfed
+          rfed(j2)=rfed(j2)+w2*valfed
+          rfed(j3)=rfed(j3)+w3*valfed
+          rfed(j4)=rfed(j4)+w4*valfed
+          rfed(j5)=rfed(j5)+w5*valfed
+          rfed(j6)=rfed(j6)+w6*valfed
+          rfed(j7)=rfed(j7)+w7*valfed
+          rfed(j8)=rfed(j8)+w8*valfed
+        end if
+
+     endif
+
+     !fedptr => fedptr%llpoint
+     fedptr => fedNode_nextcast(fedptr)
+  end do
+  return
+end subroutine intfed_
+
+end module intfedmod

src/gsi/m_berror_stats_reg.f90

@@ -400,7 +400,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt
   integer(i_kind) :: nrf2_td2m,nrf2_mxtm,nrf2_mitm,nrf2_pmsl,nrf2_howv,nrf2_tcamt,nrf2_lcbas,nrf2_cldch
   integer(i_kind) :: nrf2_uwnd10m,nrf2_vwnd10m
   integer(i_kind) :: nrf3_sfwter,nrf3_vpwter
-  integer(i_kind) :: nrf3_dbz
+  integer(i_kind) :: nrf3_dbz,nrf3_fed
   integer(i_kind) :: nrf3_ql,nrf3_qi,nrf3_qr,nrf3_qs,nrf3_qg,nrf3_qnr,nrf3_w
   integer(i_kind) :: inerr,istat
   integer(i_kind) :: nsigstat,nlatstat,isig
@@ -624,6 +624,7 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt
   nrf3_sf =getindex(cvars3d,'sf')
   nrf3_vp =getindex(cvars3d,'vp')
   nrf3_dbz=getindex(cvars3d,'dbz')
+  nrf3_fed=getindex(cvars3d,'fed')
   nrf2_sst=getindex(cvars2d,'sst')
   nrf2_gust=getindex(cvars2d,'gust')
   nrf2_vis=getindex(cvars2d,'vis')
@@ -671,6 +672,16 @@ subroutine berror_read_wgt_reg(msig,mlat,corz,corp,hwll,hwllp,vz,rlsig,varq,qopt
     vz(:,:,nrf3_dbz)=vz(:,:,nrf3_t)
   endif
 
+  if( nrf3_fed>0 )then
+    if(.not. nrf3_t>0) then
+      write(6,*)'not as expect,stop'
+      stop
+    endif
+    corz(:,:,nrf3_fed)=10.0_r_kind
+    hwll(:,:,nrf3_fed)=hwll(:,:,nrf3_t)
+    vz(:,:,nrf3_fed)=vz(:,:,nrf3_t)
+  endif
+
   if (nrf3_oz>0) then 
      factoz = 0.0002_r_kind*r25
      corz(:,:,nrf3_oz)=factoz