1. Add oneob capablity for FED DA

2. Use 3D interplotion when FED from BG is used
3. Clean code
	modified:   src/gsi/gsimod.F90
	modified:   src/gsi/obsmod.F90
	modified:   src/gsi/read_fed.f90
	modified:   src/gsi/setupfed.f90

src/gsi/gsimod.F90

@@ -23,7 +23,7 @@ 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_model_fed,innov_use_model_fed,if_vrobs_raw,if_use_w_vr,&
      minobrangedbz,maxobrangedbz,maxobrangevr,maxtiltvr,missing_to_nopcp,&
@@ -763,7 +763,7 @@ 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,&

src/gsi/obsmod.F90

@@ -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,7 +622,7 @@ module obsmod
   integer(i_kind) ntilt_radarfiles,tcp_posmatch,tcp_box
 
   logical ::  ta2tb
-  logical ::  doradaroneob
+  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
@@ -756,12 +757,13 @@ subroutine init_obsmod_dflts
     if_use_w_vr=.true.
     if_model_dbz=.false.
     if_model_fed=.false.
-    innov_use_model_fed=.true.
+    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_fed.f90

@@ -35,7 +35,7 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs)
   use gridmod, only: tll2xy
   use mod_wrfmass_to_a, only: wrfmass_obs_to_a8
   use mpimod, only: npe
-  use obsmod, only: perturb_obs,iadatemn
+  use obsmod, only: perturb_obs,iadatemn,dofedoneob,oneoblat,oneoblon
 
   use netcdf
   implicit none
@@ -389,13 +389,21 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs)
       kint_maxloc=-1
       fed_max=-999.99
       ndata2=0 
+
+      ILOOP : &
       do i=1,numfed
         do k=1,maxlvl
           if( fed3d_column(k+2,i) >= fed_lowbnd2 .or. fed3d_column(k+2,i) == fed_lowbnd) then !Rong Kong
             dlon_earth = fed3d_column(1,i)                ! longitude (degrees) of observation
                                                        ! ilone=18    ! index of longitude (degrees)
             dlat_earth = fed3d_column(2,i)                ! latitude (degrees) of observation
                                                        ! ilate=19    ! index of latitude (degrees)
+
+           if (dofedoneob) then
+             dlat_earth=oneoblat
+             dlon_earth=oneoblon
+           endif
+
            !-Check format of longitude and correct if necessary
            if(dlon_earth>=r360) dlon_earth=dlon_earth-r360
            if(dlon_earth<zero ) dlon_earth=dlon_earth+r360
@@ -405,6 +413,12 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs)
             rlon00 = dlon_earth*deg2rad
             rlat00 = dlat_earth*deg2rad
             call tll2xy(rlon00,rlat00,dlon,dlat,outside)
+           if (dofedoneob) then
+             if (outside) then
+               write(6,*)'READ_FED: ONE OB OUTSIDE; STOP2(61) ',dlat_earth,dlon_earth
+               call stop2(61)
+             end if
+           end if
             if (outside) cycle
 
                                            !If observation is outside the domain
@@ -425,7 +439,7 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs)
             cdata_out( 4,ndata2) = hgt_fed(k)          ! obs absolute height (m) above MSL
                                                        ! ipres=4     ! index of pressure
             cdata_out( 5,ndata2) = fed3d_column(k+2,i) ! FED value
-                                                       ! idbzob=5    ! index of dbz observation
+
             cdata_out( 6,ndata2) = rstation_id         ! station id (charstring equivalent to real double)
                                                        ! id=6        ! index of station id
 
@@ -472,6 +486,7 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs)
                cdata_out(26,ndata2) = 1.0_r_kind       ! obs perturbation
                                                        ! iptrb=26    ! index of q perturbation
             end if
+
 !          print*,'cdata_out(:,ndata2)=',cdata_out(:,ndata2)
             if(fed3d_column(k+2,i) > fed_max)then
                kint_maxloc=k
@@ -480,9 +495,12 @@ subroutine read_fed(nread,ndata,nodata,infile,obstype,lunout,twind,sis,nobs)
                i_maxloc=fed3d_column(1,i)
                fed_max =fed3d_column(k+2,i)
             end if
+
+            if( dofedoneob ) exit ILOOP
+
           endif
-        enddo
-      enddo
+        enddo ! k 
+      enddo ILOOP ! i
 
 !---all looping done now print diagnostic output
   write(6,*)'READ_FED: Reached eof on FED file'

src/gsi/setupfed.f90

@@ -18,9 +18,10 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa
 !                   work of Sebok and Back (2021, unpublished)
 !                 - capped maximum model FED
 !                Hongli Wang         NOAA GSL    2023-09-14
-!                 - Add option to use fed from background file to cal fed innov
-!                 - fed in BG exist (if_model_fed=.true.), and is used 
-!                    to cal innov (innov_use_model_fed=.true.)
+!                 - Add option to use fed from background file to calculate fed innov
+!                     - The bellow two namelist parameters need to be true
+!                        - if_model_fed=.true.          fed in BG exist
+!                        - innov_use_model_fed=.true.   turn on flag to use FED from BG to cal innov
 !
 !
   use mpeu_util, only: die,perr
@@ -34,7 +35,7 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa
   use obsmod, only: rmiss_single,&
        lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset                                                   
   use obsmod, only: oberror_tune
-  use obsmod, only: if_model_fed,innov_use_model_fed
+  use obsmod, only: if_model_fed,innov_use_model_fed,dofedoneob,oneobddiff,oneobvalue
   use m_obsNode, only: obsNode
   use m_fedNode, only: fedNode
   use m_fedNode, only: fedNode_appendto
@@ -235,6 +236,10 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa
      muse(i)=nint(data(iuse,i)) <= jiter
   end do
 
+  if (dofedoneob) then
+    muse=.true.
+  end if
+
   numequal=0
   numnotequal=0
 
@@ -527,9 +532,9 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa
 
      npt = 0
      FEDMdiag(i) = 0.
-     if ( if_model_fed .and. innov_use_model_fed) then
+     if (if_model_fed .and. innov_use_model_fed) then
         !use fed from background file
-        call tintrp2a11(ges_fed(:,:,1,:),FEDMdiag(i),dlat,dlon,dtime,hrdifsig,mype,nfldsig) 
+        call tintrp31(ges_fed,FEDMdiag(i),      dlat,dlon,dpres,dtime,hrdifsig,mype,nfldsig)       
      else
         call tintrp2a11(rp,FEDMdiag(i),dlat,dlon,dtime,hrdifsig,mype,nfldsig)
      end if
@@ -557,24 +562,29 @@ subroutine setupfed(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,fed_diagsa
      !--------------Calculate departure from observation----------------!
 
      ddiff = data(ifedob,i) - FEDMdiag(i)
-
-!    If requested, setup for single obs test.
-!    Note: do not use this way to run single obs test for fed in the current version. (g.zhao)
-     if (oneobtest) then
-        ddiff=maginnov
-!        if (trim(adjustl(oneob_type))=='fed') then
-!           data(ifedob,i) = maginnov
-!           ddiff = data(ifedob,i) - FEDMdiag(i)
-!        end if
+     if (dofedoneob) then
+        !use magoberr to define obs error, but oneobtest=.false.
+        if(magoberr <= zero) magoberr=1.0_r_kind 
         error=one/(magoberr)
         ratio_errors=one
-     end if
-
+        if (jiter==1) then
+          if (oneobvalue > 0_r_kind) then
+            data(ifedob,i) = oneobvalue
+            ddiff = data(ifedob,i) - FEDMdiag(i) 
+          else
+            ddiff = oneobddiff
+            data(ifedob,i) = FEDMdiag(i)+ddiff
+          endif
+          write(6,*)"FED_ONEOB: O_Val,B_Val= ",data(ifedob,i),FEDMdiag(i)
+          write(6,*)"FED_ONEOB: Innov,Error= ",ddiff,magoberr
+        else
+          ddiff = data(ifedob,i) - FEDMdiag(i)
+       end if
+     end if !oneob
 
 !    Gross error checks
      obserror = one/max(ratio_errors*error,tiny_r_kind)     
      obserrlm = max(cermin(ikx),min(cermax(ikx),obserror))
-
      residual = abs(ddiff)              != y-H(xb)      
      ratio    = residual/obserrlm       != y-H(xb)/sqrt(R)   
 
@@ -836,7 +846,6 @@ subroutine init_vars_
      !    get fed ....
          varname='fed'
          call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus)
-         print*,"FED_setupfed: ",istatus
          if (istatus==0) then
            if(allocated(ges_fed))then
             write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '