Add FED EnVar DA Capability

src/gsi/control2state.f90

@@ -282,6 +282,7 @@ subroutine control2state(xhat,sval,bval)
          endif
       enddo
    end if
+
    call gsi_bundlegetpointer (sval(jj),'ps'  ,sv_ps,  istatus)
    call gsi_bundlegetvar ( wbundle, 'ps' , sv_ps,  istatus )
 
@@ -730,6 +731,7 @@ subroutine control2state_ad(rval,bval,grad)
          endif
       enddo
    end if
+
 !  Calculate sensible temperature
    if(do_tv_to_tsen) call tv_to_tsen_ad(cv_t,rv_q,rv_tsen)
 

src/gsi/ensctl2model_ad.f90

@@ -190,7 +190,6 @@ subroutine ensctl2model_ad(eval,mval,grad)
           call gsi_bundleputvar     (wbundle_c, clouds(ic),rv_rank3,istatus)
       endif
    enddo
-
 !  Convert RHS calculations for u,v to st/vp
    if (do_getuv) then
       if(uv_hyb_ens) then

src/gsi/ensctl2state.f90

@@ -236,7 +236,6 @@ subroutine ensctl2state(xhat,mval,eval)
       enddo
    endif
 
-
 !$omp section
 
 !  Get pointers to required state variables

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,9 +23,9 @@ 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 
@@ -770,12 +770,12 @@ 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,&

src/gsi/intfed.f90

@@ -0,0 +1,189 @@
+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 gridmod, only: wrf_mass_regional, fv3_regional
+  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,valqr,valqs,valqg,valfed,valqnr
+  real(r_kind) cg_fed,p0,grad,wnotgross,wgross,pg_fed
+  real(r_kind) qrtl,qstl, qgtl, qnrtl
+  real(r_kind),pointer,dimension(:) :: sqr,sqs,sqg,sfed,sqnr
+  real(r_kind),pointer,dimension(:) :: rqr,rqs,rqg,rfed,rqnr
+  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

src/gsi/obsmod.F90

@@ -473,7 +473,7 @@ module obsmod
   ! ==== DBZ DA ===
   public :: ntilt_radarfiles
   public :: whichradar
-  public :: vr_dealisingopt, if_vterminal, if_model_dbz, inflate_obserr, if_vrobs_raw, if_use_w_vr, l2rwthin
+  public :: vr_dealisingopt, if_vterminal, if_model_dbz,if_model_fed, innov_use_model_fed, inflate_obserr, if_vrobs_raw, if_use_w_vr, l2rwthin
 
   public :: doradaroneob,oneoblat,oneoblon
   public :: oneobddiff,oneobvalue,oneobheight,oneobradid
@@ -488,6 +488,7 @@ module obsmod
   ! --- DBZ DA ---
 
   public :: iout_fed, mype_fed  
+  public :: dofedoneob
 
   public :: obsmod_init_instr_table
   public :: obsmod_final_instr_table
@@ -621,8 +622,8 @@ module obsmod
   integer(i_kind) ntilt_radarfiles,tcp_posmatch,tcp_box
 
   logical ::  ta2tb
-  logical ::  doradaroneob
-  logical :: vr_dealisingopt, if_vterminal, if_model_dbz, inflate_obserr, if_vrobs_raw, if_use_w_vr, l2rwthin
+  logical ::  doradaroneob,dofedoneob
+  logical :: vr_dealisingopt, if_vterminal, if_model_dbz,if_model_fed, innov_use_model_fed,inflate_obserr, if_vrobs_raw, if_use_w_vr, l2rwthin
   character(4) :: whichradar,oneobradid
   real(r_kind) :: oneoblat,oneoblon,oneobddiff,oneobvalue,oneobheight
   logical :: radar_no_thinning
@@ -755,11 +756,14 @@ subroutine init_obsmod_dflts
     if_vrobs_raw=.false.
     if_use_w_vr=.true.
     if_model_dbz=.false.
+    if_model_fed=.false.
+    innov_use_model_fed=.false.
     inflate_obserr=.false.
     whichradar="KKKK"
 
     oneobradid="KKKK"
     doradaroneob=.false.
+    dofedoneob=.false.
     oneoblat=-999_r_kind
     oneoblon=-999_r_kind
     oneobddiff=-999_r_kind

src/gsi/read_dbz_netcdf.f90

@@ -526,7 +526,8 @@ subroutine read_dbz_mrms_netcdf(nread,ndata,nodata,infile,obstype,lunout,sis,nob
 
 		    if(thislon>=r360) thislon=thislon-r360
                     if(thislon<zero ) thislon=thislon+r360
-
+                    if(thislon>=r360 .or. thislat >90.0_r_kind) cycle
+
 		   !-Convert back to radians                 
 
 		    thislat = thislat*deg2rad