Add FED EnVar DA Capability (NOAA-EMC#632)

- This PR supports RRFS_B GSI FED assimilation. - This PR adds a new GSI EnVar FED assimilation capability. The summary of the changes: - Read FED background and ensemble from restart phy files - Add new control/state variable of fed ( in anavinfo, section: metguess, state and control variable) - Create intfed.f90 and sfpfed.f90 for minimization. - Other related codes. For example, update hydrometers when either dbz or fed is assimilated, or both are assimilated. Previously the update of hydrometers is done only when dbz is assimilated. Please see Fixes NOAA-EMC#622 This PR was tested with: - One FED obs DA test - Real FED DA with pseudo ensemble for code development and debug - Real FED DA with real ensemble --------- Co-authored-by: David Dowell <David.Dowell@noaa.gov>
shoyokota · Oct 31, 2023 · acfe56d · acfe56d
1 parent 2cb0f5b
commit acfe56d
Show file tree

Hide file tree

Showing 13 changed files with 997 additions and 487 deletions.
diff --git a/src/gsi/cplr_get_fv3_regional_ensperts.f90 b/src/gsi/cplr_get_fv3_regional_ensperts.f90
diff --git a/src/gsi/gsi_fedOper.F90 b/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
diff --git a/src/gsi/gsi_files.cmake b/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

diff --git a/src/gsi/gsi_rfv3io_mod.f90 b/src/gsi/gsi_rfv3io_mod.f90
diff --git a/src/gsi/gsimod.F90 b/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,dbz_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 and metguess vectors. 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,20 @@ 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',330)
+  end if
+
+  if (.not. (miter == 0 .or. lobserver) .and. if_model_fed .and. .not. fed_exist) then
+     if(mype==0) write(6,*)' GSIMOD: .not. (miter == 0 .or. lobserver) and if_model_fed=.true. but fed is not in anavinfo file'
+     call die(myname_,'Please check namelist parameters and/or add fed in anavinfo (contro/state_vector and met_guess) when miter > 0 and if_model_fed=.true.',332)
+  end if
+
+  if (.not. (miter == 0 .or. lobserver) .and. if_model_dbz .and. .not. dbz_exist) then
+     if(mype==0) write(6,*)' GSIMOD: .not. (miter == 0 .or. lobserver) and if_model_dbz=.true. but dbz is not in anavinfo file'
+     call die(myname_,'Please check namelist parameters and/or add dbz 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

diff --git a/src/gsi/intfed.f90 b/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
diff --git a/src/gsi/m_berror_stats_reg.f90 b/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