diff --git a/fix b/fix
index 15ffa60307..2a86d5b1a0 160000
--- a/fix
+++ b/fix
@@ -1 +1 @@
-Subproject commit 15ffa60307bbc19746d8caeb41782de0b7604be6
+Subproject commit 2a86d5b1a09a8908618d1c8a376c68e8d9b3d02c
diff --git a/modulefiles/gsi_jet.intel.lua b/modulefiles/gsi_jet.intel.lua
index 48189ba241..b210e6efa2 100644
--- a/modulefiles/gsi_jet.intel.lua
+++ b/modulefiles/gsi_jet.intel.lua
@@ -1,7 +1,7 @@
 help([[
 ]])
 
-prepend_path("MODULEPATH", "/mnt/lfs4/HFIP/hfv3gfs/role.epic/spack-stack/spack-stack-1.6.0/envs/gsi-addon-dev-rocky8/install/modulefiles/Core")
+prepend_path("MODULEPATH", "/contrib/spack-stack/spack-stack-1.6.0/envs/gsi-addon-intel/install/modulefiles/Core")
 
 local python_ver=os.getenv("python_ver") or "3.11.6"
 local stack_intel_ver=os.getenv("stack_intel_ver") or "2021.5.0"
@@ -20,6 +20,6 @@ load(pathJoin("prod_util", prod_util_ver))
 pushenv("CFLAGS", "-axSSE4.2,AVX,CORE-AVX2")
 pushenv("FFLAGS", "-axSSE4.2,AVX,CORE-AVX2")
 
-pushenv("GSI_BINARY_SOURCE_DIR", "/mnt/lfs4/HFIP/hfv3gfs/glopara/git/fv3gfs/fix/gsi/20240208")
+pushenv("GSI_BINARY_SOURCE_DIR", "/lfs5/HFIP/hfv3gfs/glopara/FIX/fix/gsi/20240208")
 
 whatis("Description: GSI environment on Jet with Intel Compilers")
diff --git a/regression/global_enkf.sh b/regression/global_enkf.sh
index a35f8d109f..f5b4f3d1c6 100755
--- a/regression/global_enkf.sh
+++ b/regression/global_enkf.sh
@@ -76,14 +76,11 @@ export use_correlated_oberrs=${use_correlated_oberrs:-".false."}
 if [ $USE_CORRELATED_OBERRS == "YES" ]; then
     export use_correlated_oberrs=".true."
 fi
-export imp_physics=${imp_physics:-"11"}
 export lupp=${lupp:-".true."}
 export corrlength=${corrlength:-1250}
 export lnsigcutoff=${lnsigcutoff:-2.75}
 export analpertwt=${analpertwt:-0.85}
 export readin_localization_enkf=".false."
-export readin_localization_enkf=${readin_localization_enkf:-".true."}
-export reducedgrid=${reducedgrid:-".true."}
 export letkf_flag=${letkf_flag:-".true."}
 export getkf=${getkf:-".true."}
 export denkf=${denkf:-".true."}
diff --git a/regression/regression_namelists.sh b/regression/regression_namelists.sh
index a4f283f92b..a6bd78994d 100755
--- a/regression/regression_namelists.sh
+++ b/regression/regression_namelists.sh
@@ -23,12 +23,12 @@ export gsi_namelist="
    crtm_coeffs_path='./crtm_coeffs/',
    newpc4pred=.true.,adp_anglebc=.true.,angord=4,passive_bc=.true.,use_edges=.false.,
    diag_precon=.true.,step_start=1.e-3,emiss_bc=.true.,thin4d=.true.,cwoption=3,
-   verbose=.false.,imp_physics=11,lupp=.true.,
+   verbose=.false.,imp_physics=8,lupp=.true.,
    binary_diag=.false.,netcdf_diag=.true.,
    lobsdiag_forenkf=.false.,
    nhr_anal=3,6,9,nhr_obsbin=1,
    l4densvar=.true.,ens_nstarthr=3,nhr_assimilation=6,lwrite4danl=.true.,
-   optconv=0.05,cao_check=.true.,ta2tb=.false.,
+   optconv=0.05,cao_check=.true.,ta2tb=.true.,
    tzr_qc=1,sfcnst_comb=.true.,
    write_fv3_incr=.true.,incvars_to_zero= 'liq_wat_inc','icmr_inc','rwmr_inc','snmr_inc','grle_inc',
    incvars_zero_strat='sphum_inc','liq_wat_inc','icmr_inc','rwmr_inc','snmr_inc','grle_inc',incvars_efold=5, use_gfs_ncio=.true.,
@@ -962,8 +962,8 @@ export gsi_namelist="
   sprd_tol=1.e30,paoverpb_thresh=0.98,
   nlons=$LONA,nlats=$LATA,nlevs=$LEVS,nanals=$NMEM_ENKF,
   deterministic=.true.,sortinc=.true.,lupd_satbiasc=.false.,
-  reducedgrid=.true.,readin_localization=.true.,
-  use_gfs_nemsio=${use_gfs_nemsio},use_gfs_ncio=${use_gfs_ncio},imp_physics=$imp_physics,lupp=$lupp,
+  reducedgrid=.false.,readin_localization=.false.,
+  use_gfs_nemsio=${use_gfs_nemsio},use_gfs_ncio=${use_gfs_ncio},imp_physics=8,lupp=$lupp,
   univaroz=.false.,adp_anglebc=.true.,angord=4,use_edges=.false.,emiss_bc=.true.,
   letkf_flag=${letkf_flag},nobsl_max=${nobsl_max},denkf=${denkf},getkf=${getkf}.,
   nhr_anal=${IAUFHRS_ENKF},nhr_state=${IAUFHRS_ENKF},
diff --git a/regression/regression_namelists_db.sh b/regression/regression_namelists_db.sh
index b96c208070..34d80cba96 100755
--- a/regression/regression_namelists_db.sh
+++ b/regression/regression_namelists_db.sh
@@ -21,7 +21,7 @@ export gsi_namelist="
    crtm_coeffs_path='./crtm_coeffs/',
    newpc4pred=.true.,adp_anglebc=.true.,angord=4,passive_bc=.true.,use_edges=.false.,
    diag_precon=.true.,step_start=1.e-3,emiss_bc=.true.,thin4d=.true.,cwoption=3,
-   verbose=.false.,imp_physics=11,lupp=.true.,binary_diag=.false.,netcdf_diag=.true.,
+   verbose=.false.,imp_physics=8,lupp=.true.,binary_diag=.false.,netcdf_diag=.true.,
    lobsdiag_forenkf=.false.,
    nhr_anal=3,6,9,
    l4densvar=.true.,ens_nstarthr=3,nhr_obsbin=1,nhr_assimilation=6,lwrite4danl=.true., tzr_qc=1,sfcnst_comb=.true.,
diff --git a/regression/regression_param.sh b/regression/regression_param.sh
index bfc6f042fc..209762569b 100755
--- a/regression/regression_param.sh
+++ b/regression/regression_param.sh
@@ -90,8 +90,8 @@ case $regtest in
            topts[1]="0:05:00" ; popts[1]="40/3/"  ; ropts[1]="/1"
            topts[2]="0:05:00" ; popts[2]="40/5/"  ; ropts[2]="/1"
         elif [[ "$machine" = "Orion" ]]; then
-           topts[1]="0:15:00" ; popts[1]="5/4/" ; ropts[1]="/1"
-           topts[2]="0:15:00" ; popts[2]="10/4/" ; ropts[2]="/2"
+           topts[1]="0:15:00" ; popts[1]="20/6/" ; ropts[1]="/1"
+           topts[2]="0:15:00" ; popts[2]="20/12/" ; ropts[2]="/2"
         elif [[ "$machine" = "Hercules" ]]; then
            topts[1]="0:05:00" ; popts[1]="40/3/" ; ropts[1]="/1"
            topts[2]="0:05:00" ; popts[2]="40/5/" ; ropts[2]="/2"
diff --git a/src/enkf/gridio_gfs.f90 b/src/enkf/gridio_gfs.f90
index 35a0c3fbe4..6e68f1a347 100644
--- a/src/enkf/gridio_gfs.f90
+++ b/src/enkf/gridio_gfs.f90
@@ -195,9 +195,6 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, &
   sst_ind = getindex(vars2d, 'sst')
   use_full_hydro = ( ql_ind > 0 .and. qi_ind > 0  .and. &
                      qr_ind > 0 .and. qs_ind > 0 .and. qg_ind > 0 )
-  ! Currently, we do not let precipiation to affect the enkf analysis  
-  ! The following line will be removed after testing
-  use_full_hydro = .false.
 
   if (.not. isinitialized) call init_spec_vars(nlons,nlats,ntrunc,4)
 
@@ -459,9 +456,6 @@ subroutine readgriddata_pnc(vars3d,vars2d,n3d,n2d,levels,ndim,ntimes, &
   end if
 
   ! cloud derivatives
-  ! Currently, we do not let precipiation to affect the enkf analysis  
-  ! The following line will be removed after testing
-  use_full_hydro = .true.
   if (.not. use_full_hydro .and. iope==0) then
      if (ql_ind > 0 .or. qi_ind > 0) then
         do k=1,nlevs
diff --git a/src/gsi/convinfo.f90 b/src/gsi/convinfo.f90
index f6ba29c2e1..83c84d3f38 100755
--- a/src/gsi/convinfo.f90
+++ b/src/gsi/convinfo.f90
@@ -155,6 +155,8 @@ subroutine convinfo_read
 !   2009-01-22  todling - protect against non-initialized destroy call
 !   2010-05-29  todling - interface consistent w/ similar routines
 !   2014-07-10  carley  - add check to bypass blank lines in convinfo file
+!   2023-05-09  s.vetra-carvalho - extended iotype to be 8 characters long to
+!   allow for more descriptive observations names
 !
 !   input argument list:
 !     mype - mpi task id
@@ -171,7 +173,7 @@ subroutine convinfo_read
     implicit none
     
     character(len=1)cflg
-    character(len=7) iotype
+    character(len=8) iotype
     character(len=140) crecord
     integer(i_kind) lunin,i,nc,ier,istat
     integer(i_kind) nlines
@@ -188,13 +190,13 @@ subroutine convinfo_read
     nlines=0
     read1: do
        cflg=' '
-       iotype='       '
+       iotype='        '
        read(lunin,1030,iostat=istat,end=1130)cflg,iotype,crecord
-1030   format(a1,a7,2x,a140)
+1030   format(a1,a8,1x,a140)
        if (istat /= 0) exit
        nlines=nlines+1
        if(cflg == '!')cycle
-       if (cflg==' '.and.iotype=='       ') then
+       if (cflg==' '.and.iotype=='        ') then
          if(print_verbose)write(6,*) 'Encountered a blank line in convinfo file at line number: ',nlines,' skipping!'
          cycle
        end if
@@ -242,9 +244,9 @@ subroutine convinfo_read
 
     do i=1,nlines
        cflg=' '
-       iotype='       '
+       iotype='        '
        read(lunin,1030)cflg,iotype,crecord
-       if (cflg==' '.and.iotype=='       ') then
+       if (cflg==' '.and.iotype=='        ') then
          if(print_verbose)write(6,*) 'Encountered a blank line in convinfo file at line number: ',i,' skipping!'
          cycle
        end if
@@ -282,7 +284,7 @@ subroutine convinfo_read
        if(print_verbose .and. mype == 0)write(6,1031)ioctype(nc),ictype(nc),icsubtype(nc),icuse(nc),ctwind(nc),ncnumgrp(nc), &
           ncgroup(nc),ncmiter(nc),cgross(nc),cermax(nc),cermin(nc),cvar_b(nc),cvar_pg(nc), &
           ithin_conv(nc),rmesh_conv(nc),pmesh_conv(nc),idum,pmot_conv(nc),ptime_conv(nc),index_sub(nc),ibeta(nc),ikapa(nc)
-1031   format('READ_CONVINFO: ',a7,1x,i3,1x,i4,1x,i2,1x,g13.6,1x,3(I3,1x),5g13.6,i5,2g13.6,i5,2g13.6,3i5)
+1031   format('READ_CONVINFO: ',a8,1x,i3,1x,i4,1x,i2,1x,g13.6,1x,3(I3,1x),5g13.6,i5,2g13.6,i5,2g13.6,3i5)
     enddo
 
     close(lunin)
diff --git a/src/gsi/correlated_obsmod.F90 b/src/gsi/correlated_obsmod.F90
index 17cd94efe1..9c8791b0f3 100644
--- a/src/gsi/correlated_obsmod.F90
+++ b/src/gsi/correlated_obsmod.F90
@@ -422,7 +422,7 @@ subroutine set_(instrument,fname,mask,method,kreq,kmut,ErrorCov)
       ErrorCov%nch_active = coun
       if (.not.GMAO_ObsErrorCov) ErrorCov%nctot = nctot
       call create_(coun,ErrorCov)
-      allocate(indxRf(nch_active),indxR(nch_active),Rcov(nctot,nctot))
+      allocate(indxRf(coun),indxR(nch_active),Rcov(nctot,nctot))
 
 !     Read GSI-like channel numbers used in estimating R for this instrument
       read(lu,IOSTAT=ioflag) indxR
diff --git a/src/gsi/crtm_interface.f90 b/src/gsi/crtm_interface.f90
index 4bb1191001..c30a2e06c8 100644
--- a/src/gsi/crtm_interface.f90
+++ b/src/gsi/crtm_interface.f90
@@ -169,7 +169,8 @@ module crtm_interface
   real(r_kind)   , save ,allocatable,dimension(:,:) :: cloudefr     ! effective radius of cloud type in CRTM
   real(r_kind)   , save ,allocatable,dimension(:,:) :: cloud_cont   ! cloud content fed into CRTM 
   real(r_kind)   , save ,allocatable,dimension(:,:) :: cloud_efr    ! effective radius of cloud type in CRTM
-  real(r_kind)   , save ,allocatable,dimension(:)   :: cf           ! effective radius of cloud type in CRTM 
+  real(r_kind)   , save ,allocatable,dimension(:)   :: cf           ! cloud fraction
+  real(r_kind)   , save ,allocatable,dimension(:)   :: ni,nr        ! ice and rain number concentration  
   real(r_kind)   , save ,allocatable,dimension(:)   :: hwp_guess    ! column total for each hydrometeor  
   real(r_kind)   , save ,allocatable,dimension(:)  :: table,table2,tablew ! GFDL saturation water vapor pressure tables
   real(r_kind)   , save ,allocatable,dimension(:)  :: des2,desw           ! GFDL saturation water vapor presure
@@ -435,13 +436,17 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo
     allocate(cloud(nsig,n_clouds_fwd))
     allocate(cloudefr(nsig,n_clouds_fwd))
     allocate(icloud(n_actual_clouds))
-    allocate(cf(nsig))  
+    allocate(cf(nsig))
+    allocate(ni(nsig))   
+    allocate(nr(nsig))     
     allocate(hwp_guess(n_clouds_fwd))   
     cloud_cont=zero
     cloud_efr =zero
     cloud     =zero
     cloudefr  =zero
-    cf        =zero  
+    cf        =zero
+    ni        =zero
+    nr        =zero  
     hwp_guess =zero    
 
     call gsi_bundlegetpointer(gsi_metguess_bundle(1),cloud_names,icloud,ier)
@@ -859,10 +864,9 @@ subroutine init_crtm(init_pass,mype_diaghdr,mype,nchanl,nreal,isis,obstype,radmo
     endif ! nvege_type
  endif ! regional or IGBP
     
-! Initial GFDL saturation water vapor pressure tables
+! Initial GFDL saturation water vapor pressure tables 
   if (use_gfdl_qsat) then
      if (n_actual_aerosols_wk>0 .or. n_clouds_fwd_wk>0 .and. imp_physics==11) then
-
         if (mype==0) write(6,*)myname_,':initial and load GFDL saturation water vapor pressure tables'
   
         allocate(table (length))
@@ -965,7 +969,9 @@ subroutine destroy_crtm
   if(allocated(jcloud)) deallocate(jcloud)
   if(allocated(cloud_cont)) deallocate(cloud_cont)
   if(allocated(cloud_efr)) deallocate(cloud_efr)
-  if(allocated(cf)) deallocate(cf)  
+  if(allocated(cf)) deallocate(cf)
+  if(allocated(ni)) deallocate(ni)   
+  if(allocated(nr)) deallocate(nr)    
   if(allocated(hwp_guess)) deallocate(hwp_guess)   
   if(allocated(icw)) deallocate(icw)
   if(allocated(lcloud4crtm_wk)) deallocate(lcloud4crtm_wk)
@@ -1076,7 +1082,8 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
   use crtm_module, only: limit_exp,o3_id,toa_pressure
   use obsmod, only: iadate
   use aeroinfo, only: nsigaerojac
-  use chemmod, only: lread_ext_aerosol !for separate aerosol input file
+  use chemmod, only: lread_ext_aerosol !for separate aerosol input file 
+
   use crtm_cloudcover_define, only: cloudcover_maximum_overlap, &
                                     cloudcover_random_overlap,  &
                                     cloudcover_maxran_overlap,  &
@@ -1123,7 +1130,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
       reshape((/0.0_r_kind, 1.0_r_kind, 1.0_r_kind, 2.0_r_kind, 1.0_r_kind, &
                -1.0_r_kind, 1.0_r_kind, -1.0_r_kind/), (/4, 2/))
   real(r_kind),parameter:: jac_pert = 1.0_r_kind
-
+  real(r_kind),parameter:: xrc3 = 200.0_r_kind
 ! Declare local variables  
   integer(i_kind):: iquadrant  
   integer(i_kind):: ier,ii,kk,kk2,i,itype,leap_day,day_of_year
@@ -1132,12 +1139,12 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
   integer(i_kind):: i_minus, i_plus, j_minus, j_plus
   integer(i_kind):: itsig,itsigp,itsfc,itsfcp,itaer,itaerp
   integer(i_kind):: istyp00,istyp01,istyp10,istyp11
-  integer(i_kind):: iqs,iozs,icfs 
+  integer(i_kind):: iqs,iozs,icfs
+  integer(i_kind):: inis,inrs   
   integer(i_kind):: error_status_clr
   integer(i_kind):: idx700,dprs,dprs_min  
   integer(i_kind),dimension(8)::obs_time,anal_time
   integer(i_kind),dimension(msig) :: klevel
-
 ! ****************************** 
 ! Constrained indexing for lai
 ! CRTM 2.1 implementation change
@@ -1169,6 +1176,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
   real(r_kind),dimension(nsig) :: rho_air   ! density of air (kg/m3)
   real(r_kind),dimension(nsig) :: cf_calc   ! GFDL cloud fraction calculation
   real(r_kind),dimension(nsig) :: qmix      ! water vapor mixing ratio
+  real(r_kind),dimension(nsig) :: qcond
   real(r_kind),allocatable,dimension(:,:) :: tgas1d
   real(r_kind),pointer,dimension(:,:  )::psges_itsig =>NULL()
   real(r_kind),pointer,dimension(:,:  )::psges_itsigp=>NULL()
@@ -1185,8 +1193,12 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
   real(r_kind),pointer,dimension(:,:,:)::aeroges_itsig =>NULL()
   real(r_kind),pointer,dimension(:,:,:)::aeroges_itsigp=>NULL()
   real(r_kind),pointer,dimension(:,:,:)::cfges_itsig =>NULL()  
-  real(r_kind),pointer,dimension(:,:,:)::cfges_itsigp=>NULL()  
-
+  real(r_kind),pointer,dimension(:,:,:)::cfges_itsigp=>NULL()
+  real(r_kind),pointer,dimension(:,:,:)::niges_itsig =>NULL()
+  real(r_kind),pointer,dimension(:,:,:)::niges_itsigp=>NULL()
+  real(r_kind),pointer,dimension(:,:,:)::nrges_itsig =>NULL()
+  real(r_kind),pointer,dimension(:,:,:)::nrges_itsigp=>NULL()  
+  logical  :: lmfdeep2
   logical :: sea,icmask
 
   integer(i_kind),parameter,dimension(12):: mday=(/0,31,59,90,&
@@ -1194,7 +1206,6 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
   real(r_kind) ::   lai
 
   m1=mype+1
-
   if (n_clouds_fwd_wk>0) hwp_guess=zero  
   hwp_total=zero  
   theta_700=zero
@@ -1210,6 +1221,8 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
      cloud_efr = zero
      cloudefr = zero
      cf = zero
+     ni = zero 
+     nr = zero 
   endif
 
   dx  = data_s(ilat)                 ! grid relative latitude
@@ -1352,10 +1365,19 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
   iqs=istatus
   call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'q',qges_itsigp,istatus)
   iqs=iqs+istatus
-  call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'cf',cfges_itsig ,icfs)
+  call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'cf',cfges_itsig ,istatus)
   icfs=istatus
-  call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'cf',cfges_itsigp,icfs)
+  call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'cf',cfges_itsigp,istatus)
   icfs=icfs+istatus
+  call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'ni',niges_itsig,istatus)
+  inis=istatus
+  call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'ni',niges_itsigp,istatus)
+  inis=inis+istatus
+  call gsi_bundlegetpointer(gsi_metguess_bundle(itsig ),'nr',nrges_itsig,istatus)
+  inrs=istatus
+  call gsi_bundlegetpointer(gsi_metguess_bundle(itsigp),'nr',nrges_itsigp,istatus)
+  inrs=inrs+istatus
+
 
 ! Space-time interpolation of temperature (h) and q fields from sigma files
 !$omp parallel do  schedule(dynamic,1) private(k,ii,iii)
@@ -1730,8 +1752,10 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
             ges_qsat (ixp,iy ,k, itsigp)*w10+ &
             ges_qsat (ix ,iyp,k, itsigp)*w01+ &
             ges_qsat (ixp,iyp,k, itsigp)*w11)*dtsigp
+    
      c3(k)=r1000/(one-q(k))
      qmix(k)=q(k)*c3(k)  !convert specific humidity to mixing ratio
+     rh(k) = q(k)/qs(k)  !calculate relative humidity; added for cloud fraction computation
 ! Space-time interpolation of ozone(poz)
      if (iozs==0) then
          poz(k)=((ozges_itsig (ix ,iy ,k)*w00+ &
@@ -1758,10 +1782,36 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
                  cfges_itsigp(ix ,iyp,k)*w01+ &
                  cfges_itsigp(ixp,iyp,k)*w11)*dtsigp)
 
-!        Ensure ozone is greater than ozsmall
+!        Ensure cloud fraction is between 0 and 1
 
          cf(k)=min(max(zero,cf(k)),one)
      endif ! cf
+
+! Space-time interpolation of ice number concentration (ni)
+     if (n_clouds_fwd_wk>0 .and. inis==0) then
+         ni(k)=((niges_itsig (ix ,iy ,k)*w00+ &
+                 niges_itsig (ixp,iy ,k)*w10+ &
+                 niges_itsig (ix ,iyp,k)*w01+ &
+                 niges_itsig (ixp,iyp,k)*w11)*dtsig + &
+                (niges_itsigp(ix ,iy ,k)*w00+ &
+                 niges_itsigp(ixp,iy ,k)*w10+ &
+                 niges_itsigp(ix ,iyp,k)*w01+ &
+                 niges_itsigp(ixp,iyp,k)*w11)*dtsigp)
+         ni(k)=max(zero,ni(k))
+     endif ! ni
+
+! Space-time interpolation of rain number concentration (nr)
+     if (n_clouds_fwd_wk>0 .and. inrs==0) then
+         nr(k)=((nrges_itsig (ix ,iy ,k)*w00+ &
+                 nrges_itsig (ixp,iy ,k)*w10+ &
+                 nrges_itsig (ix ,iyp,k)*w01+ &
+                 nrges_itsig (ixp,iyp,k)*w11)*dtsig + &
+                (nrges_itsigp(ix ,iy ,k)*w00+ &
+                 nrges_itsigp(ixp,iy ,k)*w10+ &
+                 nrges_itsigp(ix ,iyp,k)*w01+ &
+                 nrges_itsigp(ixp,iyp,k)*w11)*dtsigp)
+         nr(k)=max(zero,nr(k))
+     endif ! nr
 ! Quantities required for MW cloudy radiance calculations
 
      if (n_clouds_fwd_wk>0) then
@@ -1822,8 +1872,17 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
      end do
   endif
 
+
+  ! Calculate Thompson effective radius for each hydrometeor 
+  if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==08 .and. lprecip_wk) then
+     do ii = 1, n_clouds_fwd_wk
+        iii=jcloud(ii)
+        call calc_thompson_reff(rho_air,h,cloud(:,ii),cloud_names(iii),cloudefr(:,ii))
+     end do
+  endif
+
   ! Calculate GFDL cloud fraction (if no cf in metguess table) based on PDF scheme 
-! if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and.  lcalc_gfdl_cfrac ) then
+  ! if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and.  lcalc_gfdl_cfrac ) then
   if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==11 .and.  lprecip_wk ) then
      cf_calc  = zero
      call calc_gfdl_cloudfrac(rho_air,h,q,cloud,hs,garea,qs,cf_calc)
@@ -1831,6 +1890,16 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
      icfs = 0        ! load cloud fraction into CRTM 
   endif
 
+  ! Calculate Thompson cloud fraction 
+  if ( icmask .and. n_clouds_fwd_wk > 0 .and. imp_physics==08 .and.  lprecip_wk ) then
+     lmfdeep2 = .true.     
+     cf_calc  = zero
+         ! sum of the cloud condensate amount for all 5 hydrometeos (ql + qi + qs + qg + qh + qr) 
+      qcond(:) = cloud(:,1) + cloud(:,2) + cloud(:,3) + cloud(:,4) + cloud(:,5)
+      call calc_thompson_cloudfrac (nsig , lmfdeep2, xrc3, prsl, qcond(:), rh, qs, cf_calc)
+     cf   = cf_calc
+     icfs = 0        ! load cloud fraction into CRTM
+  endif
 ! Interpolate level pressure to observation point for top interface
   prsi(nsig+1)=(ges_prsi(ix ,iy ,nsig+1,itsig )*w00+ &
               ges_prsi(ixp,iy ,nsig+1,itsig )*w10+ &
@@ -2046,7 +2115,7 @@ subroutine call_crtm(obstype,obstime,data_s,nchanl,nreal,ich, &
               do ii=1,n_clouds_fwd_wk
                 !cloud_cont(k,ii)=cloud(kk2,ii)*kgkg_kgm2 
                  cloud_cont(k,ii)=cloud(kk2,ii)*c6(k)
-                 if (imp_physics==11 .and. lprecip_wk .and.  cloud_cont(k,ii) > 1.0e-6_r_kind) then
+                 if (lprecip_wk .and.  cloud_cont(k,ii) > 1.0e-6_r_kind) then
                     cloud_efr(k,ii)=cloudefr(kk2,ii)
                  else
                     cloud_efr(k,ii)=zero
@@ -2597,6 +2666,187 @@ subroutine calc_gfdl_reff(rho_air,tsen,qxmr,cloud_name,reff)
 
   end subroutine calc_gfdl_reff
 
+  subroutine calc_thompson_reff(rho_air,tsen,qxmr,cloud_name,reff)
+!$$$  subprogram documentation block
+!                .      .    .                                       .
+! subprogram:    calc_thompson_reff:  calculate effective radius based on Thompson scheme
+!   prgmmr:      Azadeh.Gholoubi
+!
+! program history log:
+!
+!   2023-01-24
+!
+!   language: f90
+!
+!$$$
+!--------
+  use constants, only: zero, pi
+  implicit none
+
+! Declare passed variables
+  character(10)                 ,intent(in   ) :: cloud_name
+  real(r_kind), dimension(nsig) ,intent(in   ) :: rho_air   ! [ kg/m3  ]
+  real(r_kind), dimension(nsig) ,intent(in   ) :: tsen      ! [ K      ]
+  real(r_kind), dimension(nsig) ,intent(in   ) :: qxmr      ! [ kg/kg  ]
+  real(r_kind), dimension(nsig) ,intent(inout) :: reff      ! [ micron ]
+
+! Declare local variables
+  character(len=*), parameter :: myname_ = 'calc_thompson_reff'
+  integer(i_kind) :: k
+  integer(i_kind) :: mu_w
+  real(r_kind)    :: qx
+  real(r_kind)    :: reff_min, reff_max
+  real(r_kind):: lam_i,lam_w,lam_r, lam_g,lam_exp, am_r,am_w,am_i,am_g
+  
+  ! Parameters
+  real(r_kind), parameter :: qmin = 1.0e-12_r_kind          ! [kg/kg ]
+  !.. droplet number concentration.
+  real, parameter :: Nt_c = 100.E6_r_kind  ![m-3]
+
+  ! Parameters for water cloud
+  real(r_kind), parameter :: ccn        =  1.0e8_r_kind
+  real(r_kind), parameter :: rho_w      = 1000.0_r_kind     ! [kg/m3 ]
+  real(r_kind), parameter :: reff_w_min =    2.0_r_kind     ! previous value was 5.0_r_kind
+  real(r_kind), parameter :: reff_w_max =   25.0_r_kind     ! previous value was 10.0_r_kind
+
+  ! Parameters for ice cloud (Hemisfield and mcFarquhar 1996)
+  real(r_kind), parameter :: rho_i      = 890.0_r_kind      ! [kg/m3 ]
+  real(r_kind), parameter :: reff_i_min =    2.5_r_kind     ! previous value was 10_r_kind
+  real(r_kind), parameter :: reff_i_max = 250.0_r_kind      ! previous value was 150_r_kind
+  real(r_kind), parameter:: mu_i = 0.0_r_kind
+  real(r_kind), parameter:: ni_min = 1.0e-6_r_kind          ! minimum number concentration (ccpp-physics)
+
+  ! Parameters for  rain (Lin 1983)
+  real(r_kind), parameter :: rho_r      =    1000.0_r_kind  ! [kg/m3 ]
+  real(r_kind), parameter :: reff_r_min =       50.0_r_kind ! [micron] ! previous value was 0.0_r_kind
+  real(r_kind), parameter :: reff_r_max =   1000.0_r_kind   ! [micron] ! previous value was 10000.0_r_kind
+  real(r_kind), parameter:: mu_r = 0.0_r_kind
+  ! Parameters for snow
+  real(r_kind), parameter :: reff_s_min =       5.0_r_kind  ! [micron] ! previous value was 0.0_r_kind
+  real(r_kind), parameter :: reff_s_max =   5000.0_r_kind   ! [micron] ! previous value was 10000.0_r_kind
+  real(r_kind), parameter:: nr_min = 1.0e-6_r_kind          ! minimum number concentration (ccpp-physics)
+
+!For snow moments conversions  (from Field et al. 2005)
+  real(r_kind), dimension(10), parameter:: &
+      sa = (/ 5.065339_r_kind, -0.062659_r_kind, -3.032362_r_kind, 0.029469_r_kind, -0.000285_r_kind,      &
+     &        0.31255_r_kind,   0.000204_r_kind,  0.003199_r_kind, 0.0_r_kind,      -0.015952_r_kind/)
+      real(r_kind), dimension(10), parameter:: &
+      sb = (/ 0.476221_r_kind, -0.015896_r_kind,  0.165977_r_kind, 0.007468_r_kind, -0.000141_r_kind,      &
+     &        0.060366_r_kind,  0.000079_r_kind,  0.000594_r_kind, 0.0_r_kind,      -0.003577_r_kind/)
+  real(r_kind), parameter :: am_s      =     0.069_r_kind
+  real(r_kind), parameter :: bm_s = 2.0_r_kind
+  real(r_kind), dimension(1), parameter :: cse = (/ bm_s + 1.0_r_kind /)
+  real(r_kind) :: tc0, smob, smoc, a_, b_, loga_
+
+  ! Parameters for graupel  (Lin 1983)
+  real(r_kind), parameter :: rho_g      =     500.0_r_kind  ! [kg/m3 ]
+  real(r_kind), parameter :: reff_g_min =       150.0_r_kind  ! [micron] ! previous value was 0.0_r_kind
+  real(r_kind), parameter :: reff_g_max =   5000.0_r_kind  ! [micron] ! previous value was 10000.0_r_kind
+  real(r_kind), parameter:: mu_g = 0.0_r_kind
+  real(r_kind), parameter :: no_exp     =  10.0_r_kind
+
+!cloud water
+  if (trim(cloud_name)=='ql') then
+     am_w = rho_w*pi/6.0_r_kind
+     reff_min = reff_w_min
+     reff_max = reff_w_max
+     do k = 1, nsig
+        qx = qxmr(k) * rho_air(k)  ! convert mixing ratio (kg/kg) to water content (kg/m3)
+        if (qx > qmin) then
+           mu_w = MAX(2, MIN((NINT(1000.E6_r_kind/Nt_c) + 2), 15))
+           lam_w=exp(1.0_r_kind / 3.0_r_kind * log ((am_w*Nt_c *gamma(mu_w + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_w+1.0_r_kind))))
+
+           reff(k) = 0.5_r_kind * ((3.0_r_kind+mu_w)/lam_w)*1.0e6_r_kind
+           reff(k) = max(reff_min, min(reff_max, reff(k)))
+        else
+           reff(k) = zero
+        endif
+     enddo
+
+  ! Cloud Ice
+  else if (trim(cloud_name)=='qi') then
+     am_i = rho_i*pi/6.0_r_kind
+     reff_min = reff_i_min
+     reff_max = reff_i_max
+     do k = 1, nsig
+        qx = qxmr(k) * rho_air(k)  ! convert mixing ratio (kg/kg) to water content (kg/m3)
+        if (qx > qmin .and. ni(k)>ni_min) then
+           lam_i=exp(1.0_r_kind / 3.0_r_kind * log((am_i*ni(k) *gamma(mu_i + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_i+1.0_r_kind))))
+           reff(k) = 0.5_r_kind * (3.0_r_kind /lam_i)*1.0e6_r_kind
+           reff(k) = max(reff_min, min(reff_max, reff(k)))
+        else
+           reff(k) = zero
+        endif
+     enddo
+  !Rain
+  else if (trim(cloud_name)=='qr') then
+     am_r = rho_r*pi/6.0_r_kind
+     reff_min = reff_r_min
+     reff_max = reff_r_max
+     do k = 1, nsig
+        qx = qxmr(k) * rho_air(k)  ! convert mixing ratio (kg/kg) to water content (kg/m3)
+        if (qx > qmin .and. nr(k)>nr_min) then
+           lam_r=exp(1.0_r_kind / 3.0_r_kind * log ((am_r*nr(k) *gamma(mu_r + 3.0_r_kind + 1.0_r_kind))/(qx*gamma(mu_r + 1.0_r_kind))))
+           reff(k) = 0.5_r_kind *(3.0_r_kind/lam_r)*1.0e6_r_kind
+           reff(k) = max(reff_min, min(reff_max, reff(k)))
+        else
+           reff(k) = zero
+        endif
+     enddo
+
+! Snow (Field et al. 2005)
+
+  else if (trim(cloud_name)=='qs') then
+     reff_min = reff_s_min
+     reff_max = reff_s_max
+     do k = 1, nsig
+!..Calculate bm_s+1 (th) moment, smoc.  Useful for diameter calcs.
+        tc0 = MIN(-0.1_r_kind, tsen(k)-273.15_r_kind)
+        qx = qxmr(k) * rho_air(k)  ! convert mixing ratio (kg/kg) to water content (kg/m3)
+        if (qx > qmin) then
+           smob =qx/am_s
+           loga_ = sa(1) + sa(2)*tc0 + sa(3)*cse(1) &
+                 + sa(4)*tc0*cse(1) + sa(5)*tc0*tc0 &
+                 + sa(6)*cse(1)*cse(1) + sa(7)*tc0*tc0*cse(1) &
+                 + sa(8)*tc0*cse(1)*cse(1) + sa(9)*tc0*tc0*tc0 &
+                 + sa(10)*cse(1)*cse(1)*cse(1)
+           a_ = 10.0_r_kind**loga_
+           b_ = sb(1)+ sb(2)*tc0 + sb(3)*cse(1) + sb(4)*tc0*cse(1) &
+              + sb(5)*tc0*tc0 + sb(6)*cse(1)*cse(1) &
+              + sb(7)*tc0*tc0*cse(1) + sb(8)*tc0*cse(1)*cse(1) &
+              + sb(9)*tc0*tc0*tc0 + sb(10)*cse(1)*cse(1)*cse(1)
+           smoc = a_ * smob**b_
+           reff(k) = MAX(2.51E-6_r_kind, MIN(0.5_r_kind*(smoc/smob), 1999.E-6_r_kind))
+           reff(k) = max(reff_min, min(reff_max, reff(k)*1.0e6_r_kind))
+        else
+           reff(k) = zero
+        endif
+     enddo
+
+  ! Graupel
+  else if (trim(cloud_name)=='qg') then
+     reff_min = reff_g_min
+     reff_max = reff_g_max
+     am_g = rho_g*pi/6.0_r_kind
+     do k = 1, nsig
+        qx = qxmr(k)*rho_air(k)    ! convert mixing ratio (kg/kg) to water content (kg/m3)
+        if (qx > qmin) then
+           lam_exp = no_exp* exp(1.0_r_kind/(3.0_r_kind+1.0_r_kind) * log ((am_g*gamma(3.0_r_kind +1.0_r_kind))/qx))
+           lam_g=lam_exp*exp(1/3.0_r_kind*log(gamma(3.0_r_kind+mu_g+1.0_r_kind)/((3.0_r_kind+mu_g+1.0_r_kind)*(mu_g+1.0_r_kind))))
+           reff(k) = 0.5_r_kind *(3.0_r_kind/ lam_g)*1.0e6_r_kind
+           reff(k) = max(reff_min, min(reff_max, reff(k)))
+        else
+           reff(k) = zero
+        endif
+     enddo
+
+  ! Mysterious
+  else
+     call die(myname_,"cannot recognize cloud name <"//trim(myname_)//">")
+  endif
+
+  end subroutine calc_thompson_reff
+
   subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,qsat,cfrac)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
@@ -2810,6 +3060,80 @@ subroutine calc_gfdl_cloudfrac(den,pt1,qv,cloud,hs,area,qsat,cfrac)
   enddo !k-loop
   end subroutine calc_gfdl_cloudfrac
 
+  subroutine calc_thompson_cloudfrac(nsig, lmfdeep2, xrc3, prsl, clwf, rhly, qstl,cldtot)
+!$$$  subprogram documentation block
+!                .      .    .                                       .
+! subprogram:    calc_thompson_cloudfrac  calculate cloud fractions from cloud using Thompson scheme
+!
+!   prgmmr:      Azadeh.Gholoubi
+!
+! abstract:
+!
+! program history log:
+!
+!   2023-01-24  
+!
+!   language: f90
+!   
+!  ====================  definition of variables  ====================  !
+! input variables:                                                      !
+!   prsl  (msig) : model layer mean pressure in cb                      !
+!   qstl  (nsig) : layer saturate humidity in gm/gm                     !
+!   rhly  (nsig) : layer relative humidity (=qlyr/qstl)                 !
+!   nsig         : vertical layer/level dimensions                      !
+!   lmfdeep2     : logical - true for mass flux deep convection         !
+!   clwf         : sum of the all 5 Hydrometeors                        !
+!   lmfdeep2     : scale-aware mass-flux deep conv scheme flag          !                                                        !
+! output variables:                                                     !
+!   cldtot (nsig) - layer total cloud fraction                          !
+!  ====================    end of description    =====================  !
+!$$$
+!--------
+  use constants
+  implicit none
+
+!  ---  inputs:
+  integer(i_kind), intent(in) :: nsig
+  real(r_kind) , dimension(nsig), intent(in) :: prsl, clwf, rhly, qstl  
+  logical, intent(in) :: lmfdeep2
+  real(r_kind) , intent(in):: xrc3
+
+!  ---  outputs
+  real(r_kind) , dimension(nsig), intent(inout) :: cldtot
+
+!  ---  local variables:
+  real(r_kind)  :: clwmin, clwm, clwt, onemrh, value, tem1, tem2
+  integer(i_kind) :: k
+ 
+!> - Compute layer cloud fraction.
+  clwmin = 0.0_r_kind
+  do k = 1, nsig-1
+! converting prsl from cb to bar     
+     clwt = 1.0e-10_r_kind * (prsl(k)*0.01_r_kind)
+     if (clwf(k) > clwt) then
+        if(rhly(k) > 0.99_r_kind) then
+           cldtot(k) = 1.0_r_kind
+        else
+           onemrh= max( 1.e-10_r_kind, 1.0_r_kind-rhly(k) )
+           clwm  = clwmin / max( 0.01_r_kind, prsl(k)*0.01_r_kind )
+           tem1  = min(max((onemrh*qstl(k))**0.49_r_kind,0.0001_r_kind),1.0_r_kind) 
+           if (lmfdeep2) then
+              tem1  = xrc3 / tem1
+           else
+              tem1  = 100.0_r_kind / tem1
+           endif 
+           value = max( min( tem1*(clwf(k)-clwm), 50.0_r_kind ), 0.0_r_kind )
+           tem2  = sqrt( sqrt(rhly(k)) )
+           cldtot(k) = max( tem2*(1.0_r_kind-exp(-value)), 0.0_r_kind )
+        endif
+     else
+        cldtot(k) = 0.0_r_kind
+     endif
+  enddo
+
+  end subroutine calc_thompson_cloudfrac
+!............................................ 
+
 subroutine qs_table(n)
 !$$$  subprogram documentation block
 !                .      .    .                                       .
diff --git a/src/gsi/deter_sfc_mod.f90 b/src/gsi/deter_sfc_mod.f90
index 271e81c5d2..2fa85f8ba5 100644
--- a/src/gsi/deter_sfc_mod.f90
+++ b/src/gsi/deter_sfc_mod.f90
@@ -1410,7 +1410,7 @@ subroutine deter_sfc_gmi(dlat_earth,dlon_earth,isflg)
      n_grid=int(40000 / grid_dist) + 1
      klatn = max(klat1 - n_grid, 1)
      klonn = klon1 - n_grid
-     if (klonn < 0)  klonn = nlon_sfc - klonn
+     if (klonn < 1)  klonn = nlon_sfc - klonn
      klatpn = min((klat1 + n_grid), nlat_sfc)
      klonpn = klon1 + n_grid
      if (klonpn > nlon_sfc)  klonpn = klonpn - nlon_sfc
diff --git a/src/gsi/general_read_gfsatm.f90 b/src/gsi/general_read_gfsatm.f90
index 2971f27148..e1a5406a13 100755
--- a/src/gsi/general_read_gfsatm.f90
+++ b/src/gsi/general_read_gfsatm.f90
@@ -189,13 +189,14 @@ subroutine general_reload(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz,g_cwmr,
 
 end subroutine general_reload
 subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-           g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vdflag,g_cf)
+           g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vdflag,g_ni,g_nr,g_cf)
 
 ! !USES:
 
   use kinds, only: r_kind,i_kind
   use mpimod, only: npe,mpi_comm_world,ierror,mpi_rtype
   use general_sub2grid_mod, only: sub2grid_info
+  use ncepnems_io, only: imp_physics
   implicit none
 
 ! !INPUT PARAMETERS:
@@ -211,7 +212,8 @@ subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
   real(r_kind),dimension(grd%lat2,grd%lon2),         intent(  out) :: g_ps
   real(r_kind),dimension(grd%lat2,grd%lon2),         intent(inout) :: g_z
   real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent(  out) :: g_u,g_v,&
-       g_vor,g_div,g_q,g_oz,g_tv,g_ql,g_qi,g_qr,g_qs,g_qg
+       g_vor,g_div,g_q,g_oz,g_tv,g_ql,g_qi,g_qr,g_qs,g_qg,g_ni,g_nr
+
   real(r_kind),dimension(grd%lat2,grd%lon2,grd%nsig),intent(  out),optional :: g_cf
 
 
@@ -392,7 +394,25 @@ subroutine general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
                g_qg(i,j,klev)=sub(ij,k)
             enddo
          enddo
-      elseif ( iflag(k) == 15 .and. present(g_cf) ) then  
+      elseif ( iflag(k) == 15 .and. imp_physics == 8) then
+         klev=ilev(k)
+         ij=0
+         do j=1,grd%lon2
+            do i=1,grd%lat2
+               ij=ij+1
+               g_ni(i,j,klev)=sub(ij,k)
+            enddo
+         enddo
+      elseif ( iflag(k) == 16 .and. imp_physics == 8) then
+         klev=ilev(k)
+         ij=0
+         do j=1,grd%lon2
+            do i=1,grd%lat2
+               ij=ij+1
+               g_nr(i,j,klev)=sub(ij,k)
+            enddo
+         enddo
+      elseif ( iflag(k) == 17 .and. present(g_cf) ) then 
          klev=ilev(k)
          ij=0
          do j=1,grd%lon2
@@ -2828,7 +2848,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
    real(r_kind),pointer,dimension(:,:)       :: g_t2m, g_q2m
    real(r_kind),pointer,dimension(:,:,:)     :: g_vor,g_div,&
                                                 g_q,g_oz,g_tv
-   real(r_kind),pointer,dimension(:,:,:)     :: g_ql,g_qi,g_qr,g_qs,g_qg
+   real(r_kind),pointer,dimension(:,:,:)     :: g_ql,g_qi,g_qr,g_qs,g_qg,g_ni,g_nr
 
    real(r_kind),allocatable,dimension(:,:)   :: g_z
    real(r_kind),allocatable,dimension(:,:,:) :: g_u,g_v
@@ -3038,6 +3058,8 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
        call gsi_bundlegetpointer(gfs_bundle,'qs',g_qs  ,ier);istatus1=istatus1+ier
        call gsi_bundlegetpointer(gfs_bundle,'qg',g_qg  ,ier);istatus1=istatus1+ier
     !  call gsi_bundlegetpointer(gfs_bundle,'cf',g_cf  ,ier);istatus1=istatus1+ier
+       call gsi_bundlegetpointer(gfs_bundle,'ni',g_ni  ,ier);istatus1=istatus1+ier
+       call gsi_bundlegetpointer(gfs_bundle,'nr',g_nr  ,ier);istatus1=istatus1+ier
        if ( istatus1 /= 0 ) then
           if ( mype == 0 ) then
              write(6,*) 'general_read_gfsatm_allhydro_nc: ERROR'
@@ -3100,7 +3122,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
       endif
       if ( icount == icm ) then
          call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-              g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+              g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
       endif
    endif
 
@@ -3135,7 +3157,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
        endif
        if ( icount == icm ) then
           call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-               g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+               g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
        endif
 
        !   Thermodynamic variable:  s-->g transform, communicate to all tasks
@@ -3172,7 +3194,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm ) then
              call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
           endif
        end do
 
@@ -3229,7 +3251,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
              endif
              if ( icount == icm ) then
                 call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-                     g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+                     g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
              endif
           end do
           do k=1,nlevs
@@ -3285,7 +3307,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
              endif
              if ( icount == icm ) then
                 call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-                     g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+                     g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
              endif
 
           end do
@@ -3321,7 +3343,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
              endif
              if ( icount == icm ) then
                 call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-                     g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+                     g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
              endif
 
              icount=icount+1
@@ -3352,7 +3374,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
              endif
              if ( icount == icm ) then
                 call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-                     g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+                     g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
              endif
           end do
        endif ! if ( uvflag )
@@ -3383,7 +3405,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm ) then
              call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
           endif
        end do
       
@@ -3414,7 +3436,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm ) then
              call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
           endif
        end do
 
@@ -3445,7 +3467,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
 
           if ( icount == icm ) then
              call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
           endif
        enddo ! do k=1,nlevs
 
@@ -3475,7 +3497,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm ) then
              call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
           endif
        enddo ! do k=1,nlevs
 
@@ -3505,7 +3527,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm ) then
              call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
           endif
        enddo ! do k=1,nlevs
 
@@ -3535,7 +3557,7 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm ) then
              call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
           endif
        enddo ! do k=1,nlevs
 
@@ -3565,40 +3587,101 @@ subroutine general_read_gfsatm_allhydro_nc(grd,sp_a,filename,uvflag,vordivflag,z
           endif
           if ( icount == icm .or. k==nlevs) then
              call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag)
+                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
           endif
        enddo ! do k=1,nlevs
 
-!   do k=1,nlevs
-!      icount=icount+1
-!      iflag(icount)=15
-!      ilev(icount)=k
-!      kr = levs+1-k ! netcdf is top to bottom, need to flip
-!
-!      if (mype==mype_use(icount)) then
-!         call read_vardata(filges, 'cld_amt', rwork3d0, nslice=kr, slicedim=3)
-!         ! Cloud amount (cloud fraction). 
-!         if ( diff_res ) then
-!            grid_b=rwork3d0(:,:,1)
-!            vector(1)=.false.
-!            call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
-!            call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
-!            do kk=1,grd%itotsub
-!               i=grd%ltosi_s(kk)
-!               j=grd%ltosj_s(kk)
-!               work(kk)=grid2(i,j,1)
-!            enddo
-!         else
-!            grid=rwork3d0(:,:,1)
-!            call general_fill_ns(grd,grid,work)
-!         endif
+       ! Read fields specific to Thompson microphysics
+       if (imp_physics == 8) then
+          do k=1,nlevs
+             icount=icount+1
+             iflag(icount)=15
+             ilev(icount)=k
+             kr = levs+1-k ! netcdf is top to bottom, need to flip
+
+             if (mype==mype_use(icount)) then
+                call read_vardata(filges, 'nccice', rwork3d0, nslice=kr, slicedim=3)
+                ! cloud ice water number concentration.
+                if ( diff_res ) then
+                   grid_b=rwork3d0(:,:,1)
+                   vector(1)=.false.
+                   call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
+                   call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
+                   do kk=1,grd%itotsub
+                      i=grd%ltosi_s(kk)
+                      j=grd%ltosj_s(kk)
+                      work(kk)=grid2(i,j,1)
+                   enddo
+                else
+                   grid=rwork3d0(:,:,1)
+                   call general_fill_ns(grd,grid,work)
+                endif
+             endif
+             if ( icount == icm .or. k==nlevs ) then
+                call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
+                     g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
+             endif
+          enddo ! do k=1,nlevs
+
+          do k=1,nlevs
+             icount=icount+1
+             iflag(icount)=16
+             ilev(icount)=k
+             kr = levs+1-k ! netcdf is top to bottom, need to flip
+
+             if (mype==mype_use(icount)) then
+                call read_vardata(filges, 'nconrd', rwork3d0, nslice=kr, slicedim=3)
+                ! rain number concentration.
+                if ( diff_res ) then
+                   grid_b=rwork3d0(:,:,1)
+                   vector(1)=.false.
+                   call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
+                   call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
+                   do kk=1,grd%itotsub
+                      i=grd%ltosi_s(kk)
+                      j=grd%ltosj_s(kk)
+                      work(kk)=grid2(i,j,1)
+                   enddo
+                else
+                   grid=rwork3d0(:,:,1)
+                   call general_fill_ns(grd,grid,work)
+                endif
+             endif
+             if ( icount == icm .or. k==nlevs) then
+                call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
+                     g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_ni,g_nr)
+                endif
+          enddo ! do k=1,nlevs
+       endif ! imp_physics
+!       do k=1,nlevs
+!          icount=icount+1
+!          iflag(icount)=17
+!          ilev(icount)=k
+!          kr = levs+1-k ! netcdf is top to bottom, need to flip
 !
-!      endif
-!      if ( icount == icm .or. k==nlevs ) then
-!         call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
-!              g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_cf)
-!      endif
-!   enddo ! do k=1,nlevs
+!          if (mype==mype_use(icount)) then
+!             call read_vardata(filges, 'cld_amt', rwork3d0, nslice=kr, slicedim=3)
+!             ! Cloud amount (cloud fraction). 
+!             if ( diff_res ) then
+!                grid_b=rwork3d0(:,:,1)
+!                vector(1)=.false.
+!                call fill2_ns(grid_b,grid_c(:,:,1),latb+2,lonb)
+!                call g_egrid2agrid(p_high,grid_c,grid2,1,1,vector)
+!                do kk=1,grd%itotsub
+!                   i=grd%ltosi_s(kk)
+!                   j=grd%ltosj_s(kk)
+!                   work(kk)=grid2(i,j,1)
+!                enddo
+!             else
+!                grid=rwork3d0(:,:,1)
+!                call general_fill_ns(grd,grid,work)
+!             endif
+!          endif
+!          if ( icount == icm .or. k==nlevs ) then
+!             call general_reload2(grd,g_z,g_ps,g_tv,g_vor,g_div,g_u,g_v,g_q,g_oz, &
+!                  g_ql,g_qi,g_qr,g_qs,g_qg,icount,iflag,ilev,work,uvflag,vordivflag,g_cf)
+!          endif
+!       enddo ! do k=1,nlevs
 
    else ! read_2m
 
diff --git a/src/gsi/gsi_files.cmake b/src/gsi/gsi_files.cmake
index 95d885e2ee..97dcc3101f 100644
--- a/src/gsi/gsi_files.cmake
+++ b/src/gsi/gsi_files.cmake
@@ -218,6 +218,7 @@ gsi_dbzOper.F90
 gsi_dwOper.F90
 gsi_enscouplermod.f90
 gsi_fedOper.F90
+gsi_gnssrspdOper.F90
 gsi_gpsbendOper.F90
 gsi_gpsrefOper.F90
 gsi_gustOper.F90
@@ -276,6 +277,7 @@ intco.f90
 intdbz.f90
 intfed.f90
 intdw.f90
+intgnssrspd.f90
 intgps.f90
 intgust.f90
 inthowv.f90
@@ -342,6 +344,7 @@ m_dwNode.F90
 m_extOzone.F90
 m_fedNode.F90
 m_find.f90
+m_gnssrspdNode.F90
 m_gpsNode.F90
 m_gpsrhs.F90
 m_gsiBiases.f90
@@ -486,6 +489,7 @@ read_files.f90
 read_fl_hdob.f90
 read_gfs_ozone_for_regional.f90
 read_gmi.f90
+read_gnssrspd.f90
 read_goesglm.f90
 read_goesimg.f90
 read_goesimgr_skycover.f90
@@ -537,6 +541,7 @@ setupdbz.f90
 setupdbz_lib.f90
 setupdw.f90
 setupfed.f90
+setupgnssrspd.f90
 setupgust.f90
 setuphowv.f90
 setuplag.f90
@@ -597,6 +602,7 @@ stpco.f90
 stpdbz.f90
 stpfed.f90
 stpdw.f90
+stpgnssrspd.f90
 stpgps.f90
 stpgust.f90
 stphowv.f90
diff --git a/src/gsi/gsi_gnssrspdOper.F90 b/src/gsi/gsi_gnssrspdOper.F90
new file mode 100644
index 0000000000..92f0d0a2a6
--- /dev/null
+++ b/src/gsi/gsi_gnssrspdOper.F90
@@ -0,0 +1,161 @@
+module gsi_gnssrspdOper
+!$$$  subprogram documentation block
+!                .      .    .                                       .
+! subprogram:	 module gsi_gnssrspdOper
+!                (based on gsi_*Oper.F90 routines by j guo <jguo@nasa.gov>)
+!   prgmmr:	 K. Apodaca <Karina.Apodaca@gcc.spire.com>
+!      org:      Spire Global, Inc. 
+!     date:      2023-04-28
+!
+! abstract: an obOper extension for gnssrspdNode type
+!
+! program history log:
+!   2023-04-21  k apodaca - initial version
+!   input argument list: see Fortran 90 style document below
+!
+!   output argument list: see Fortran 90 style document below
+!
+! attributes:
+!   language: Fortran 90 and/or above
+!   machine:
+!
+!$$$  end subprogram documentation block
+
+! module interface:
+
+  use gsi_obOper, only: obOper
+  use m_gnssrspdNode , only: gnssrspdNode
+  implicit none
+  public:: gnssrspdOper      ! data stracture
+
+  type,extends(obOper):: gnssrspdOper
+  contains
+    procedure,nopass:: mytype
+    procedure,nopass:: nodeMold
+    procedure:: setup_
+    procedure:: intjo1_
+    procedure:: stpjo1_
+  end type gnssrspdOper
+
+!~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+  character(len=*),parameter :: myname='gsi_gnssrspdOper'
+  type(gnssrspdNode),save,target:: myNodeMold_
+
+contains
+  function mytype(nodetype)
+    implicit none
+    character(len=:),allocatable:: mytype
+    logical,optional, intent(in):: nodetype
+    mytype="[gnssrspdOper]"
+    if(present(nodetype)) then
+      if(nodetype) mytype=myNodeMold_%mytype()
+    endif
+  end function mytype
+
+  function nodeMold()
+  !> %nodeMold() returns a mold of its corresponding obsNode
+    use m_obsNode, only: obsNode
+    implicit none
+    class(obsNode),pointer:: nodeMold
+    nodeMold => myNodeMold_
+  end function nodeMold
+
+  subroutine setup_(self, lunin, mype, is, nobs, init_pass,last_pass)
+    use gnssrspd_setup, only: setup
+    use kinds, only: i_kind
+    use gsi_obOper, only: len_obstype
+    use gsi_obOper, only: len_isis
+
+    use m_rhs  , only: awork => rhs_awork
+    use m_rhs  , only: bwork => rhs_bwork
+    use m_rhs  , only: iwork => i_gnssrspd
+
+    use obsmod  , only: write_diag
+    use convinfo, only: diag_conv
+    use jfunc   , only: jiter
+
+    use mpeu_util, only: die
+    implicit none
+    class(gnssrspdOper ), intent(inout):: self
+    integer(i_kind), intent(in):: lunin
+    integer(i_kind), intent(in):: mype
+    integer(i_kind), intent(in):: is
+    integer(i_kind), intent(in):: nobs
+    logical        , intent(in):: init_pass     ! supporting multi-pass setup()
+    logical        , intent(in):: last_pass     ! with incremental backgrounds.
+
+    !----------------------------------------
+    character(len=*),parameter:: myname_=myname//"::setup_"
+
+    character(len=len_obstype):: obstype
+    character(len=len_isis   ):: isis
+    integer(i_kind):: nreal,nchanl,ier,nele
+    logical:: diagsave
+
+    if(nobs == 0) return
+
+    read(lunin,iostat=ier) obstype,isis,nreal,nchanl
+    if(ier/=0) call die(myname_,'read(obstype,...), iostat =',ier)
+    nele = nreal+nchanl
+
+    diagsave  = write_diag(jiter) .and. diag_conv
+
+    call setup(self%obsLL(:), self%odiagLL(:), &
+        lunin,mype,bwork,awork(:,iwork),nele,nobs,is,diagsave)
+
+  end subroutine setup_
+
+  subroutine intjo1_(self, ibin, rval,sval, qpred,sbias)
+    use intgnssrspdmod, only: intjo => intgnssrspd
+    use gsi_bundlemod  , only: gsi_bundle
+    use bias_predictors, only: predictors
+    use m_obsNode , only: obsNode
+    use m_obsLList, only: obsLList_headNode
+    use kinds     , only: i_kind, r_quad
+    implicit none
+    class(gnssrspdOper  ),intent(in   ):: self
+    integer(i_kind ),intent(in   ):: ibin
+    type(gsi_bundle),intent(inout):: rval   ! (ibin)
+    type(gsi_bundle),intent(in   ):: sval   ! (ibin)
+    real(r_quad    ),target,dimension(:),intent(inout):: qpred  ! (ibin)
+    type(predictors),target,             intent(in   ):: 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 stpgnssrspdmod, only: stpjo => stpgnssrspd
+    use gsi_bundlemod, only: gsi_bundle
+    use bias_predictors, only: predictors
+    use m_obsNode , only: obsNode
+    use m_obsLList, only: obsLList_headNode
+    use kinds, only: r_quad,r_kind,i_kind
+    implicit none
+    class(gnssrspdOper  ),intent(in):: self
+    integer(i_kind ),intent(in):: ibin
+    type(gsi_bundle),intent(in):: dval
+    type(gsi_bundle),intent(in):: xval
+    real(r_quad    ),dimension(:),intent(inout):: pbcjo ! (1:4)
+    real(r_kind    ),dimension(:),intent(in   ):: sges
+    integer(i_kind),intent(in):: nstep
+
+    type(predictors),target, intent(in):: dbias
+    type(predictors),target, intent(in):: 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_gnssrspdOper
diff --git a/src/gsi/gsi_obOperTypeManager.F90 b/src/gsi/gsi_obOperTypeManager.F90
index 6db7921905..dfcde80a32 100644
--- a/src/gsi/gsi_obOperTypeManager.F90
+++ b/src/gsi/gsi_obOperTypeManager.F90
@@ -12,7 +12,8 @@ module gsi_obOperTypeManager
 !   2018-07-12  j guo   - a type-manager for all obOper extensions.
 !                       - an enum mapping of obsinput::dtype(:) to obOper type
 !                         extensions.
-!
+!   2022-03-15  k apodaca - add GNSS-R L2 ocean wind speed type-manager
+!   2023-03-20  k apodaca - add GNSS-R DDM type-manager !
 !   input argument list: see Fortran 90 style document below
 !
 !   output argument list: see Fortran 90 style document below
@@ -52,6 +53,7 @@ module gsi_obOperTypeManager
   use gsi_radOper     , only: radOper
   use gsi_rwOper      , only: rwOper
   use gsi_spdOper     , only: spdOper
+  use gsi_gnssrspdOper  , only: gnssrspdOper
   use gsi_sstOper     , only: sstOper
   use gsi_swcpOper    , only: swcpOper
   use gsi_tcamtOper   , only: tcamtOper
@@ -102,6 +104,7 @@ module gsi_obOperTypeManager
   public:: iobOper_w
   public:: iobOper_q
   public:: iobOper_spd
+  public:: iobOper_gnssrspd
   public:: iobOper_rw
   public:: iobOper_dw
   public:: iobOper_sst
@@ -148,6 +151,7 @@ module gsi_obOperTypeManager
     enumerator:: iobOper_w
     enumerator:: iobOper_q
     enumerator:: iobOper_spd
+    enumerator:: iobOper_gnssrspd
     enumerator:: iobOper_rw
     enumerator:: iobOper_dw
     enumerator:: iobOper_sst
@@ -210,6 +214,7 @@ module gsi_obOperTypeManager
   type(      wOper), target, save::        wOper_mold
   type(      qOper), target, save::        qOper_mold
   type(    spdOper), target, save::      spdOper_mold
+  type( gnssrspdOper), target, save::   gnssrspdOper_mold
   type(     rwOper), target, save::       rwOper_mold
   type(     dwOper), target, save::       dwOper_mold
   type(    sstOper), target, save::      sstOper_mold
@@ -264,6 +269,7 @@ function dtype2index_(dtype) result(index_)
 
   case("q"      ,"[qoper]"      ); index_= iobOper_q
   case("spd"    ,"[spdoper]"    ); index_= iobOper_spd
+  case("gnssrspd" ,"[gnssrspdoper]" ); index_= iobOper_gnssrspd
   case("rw"     ,"[rwoper]"     ); index_= iobOper_rw
   case("dw"     ,"[dwoper]"     ); index_= iobOper_dw
   case("sst"    ,"[sstoper]"    ); index_= iobOper_sst
@@ -293,7 +299,7 @@ function dtype2index_(dtype) result(index_)
     case("ompslp"   ); index_= iobOper_o3l
     case("ompslpuv" ); index_= iobOper_o3l
     case("ompslpvis"); index_= iobOper_o3l
-    case("ompslpnc" ); index_= iobOper_o3l
+    case("ompslpnc" ); index_= iobOper_o3l     
 
   case("gpsbend","[gpsbendoper]"); index_= iobOper_gpsbend
     case("gps_bnd"); index_= iobOper_gpsbend
@@ -457,6 +463,7 @@ function index2vmold_(iobOper) result(vmold_)
   case(iobOper_w        ); vmold_ =>       wOper_mold
   case(iobOper_q        ); vmold_ =>       qOper_mold
   case(iobOper_spd      ); vmold_ =>     spdOper_mold
+  case(iobOper_gnssrspd   ); vmold_ =>  gnssrspdOper_mold
   case(iobOper_rw       ); vmold_ =>      rwOper_mold
   case(iobOper_dw       ); vmold_ =>      dwOper_mold
   case(iobOper_sst      ); vmold_ =>     sstOper_mold
@@ -573,6 +580,7 @@ subroutine cobstype_config_()
     cobstype(iobOper_w          )  ="wind                " ! w_ob_type
     cobstype(iobOper_q          )  ="moisture            " ! q_ob_type
     cobstype(iobOper_spd        )  ="wind speed          " ! spd_ob_type
+    cobstype(iobOper_gnssrspd   )  ="gnss-r wind speed   " ! gnssrspd_ob_type
     cobstype(iobOper_rw         )  ="radial wind         " ! rw_ob_type
     cobstype(iobOper_dw         )  ="doppler wind        " ! dw_ob_type
     cobstype(iobOper_sst        )  ="sst                 " ! sst_ob_type
diff --git a/src/gsi/gsimain.f90 b/src/gsi/gsimain.f90
index 5eaa38286a..5c56259bc7 100644
--- a/src/gsi/gsimain.f90
+++ b/src/gsi/gsimain.f90
@@ -131,6 +131,8 @@ program gsi
 !                         related to the GOES/GLM lightnig assimilation
 !   2019-07-09  todling - add initialization of abstract layer defining use of GFS ensemble
 !   2019-08-04  guo     - moved ensemble object configuration into module gsi_fixture.
+!   2022-03-15  K Apodaca - add CYGNSS and Spire Ocean wind speed observations
+!   2023-03-15  K Apodaca - add GNSS-R Doppler Delay Map observations
 !
 ! usage:
 !   input files:
@@ -194,7 +196,7 @@ program gsi
 !         read_goesglm, read_goesndr, read_gps_ref, read_guess, read_ieeetovs, 
 !         read_lidar, read_obs, read_ozone, read_pcp, read_prepbufr, read_radar, 
 !         read_superwinds, read_wrf_mass_files, read_wrf_mass_guess, 
-!         read_wrf_nmm_files, read_wrf_nmm_guess, rfdpar, rsearch, satthin,
+!         read_wrf_nmm_files, read_wrf_nmm_guess, read_gnssrspd, rfdpar, rsearch, satthin,
 !         setupdw, setupoz, setuppcp, setupps, setuppw, setupq, setuprad,
 !         setupref, setupbend, setuplight, setuprhsall, setuprw, setupspd, setupsst,
 !         setupt, setupw, simpin1, simpin1_init, smooth121, smoothrf,
diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90
index d7f5667252..8c6977d78c 100644
--- a/src/gsi/gsimod.F90
+++ b/src/gsi/gsimod.F90
@@ -2320,7 +2320,7 @@ subroutine gsimain_initialize
      endif
      do i=1,ndat
         write(6,401)dfile(i),dtype(i),dplat(i),dsis(i),dval(i),dthin(i),dsfcalc(i),time_window(i)
- 401    format(1x,a20,1x,a10,1x,a10,1x,a20,1x,f10.2,1x,I3,1x,I3,1x,f10.2)
+ 401    format(1x,a20,1x,a10,1x,a11,1x,a20,1x,f10.2,1x,I3,1x,I3,1x,f10.2)
      end do
      write(6,superob_radar)
      write(6,lag_data)
diff --git a/src/gsi/intgnssrspd.f90 b/src/gsi/intgnssrspd.f90
new file mode 100644
index 0000000000..1310ec7e14
--- /dev/null
+++ b/src/gsi/intgnssrspd.f90
@@ -0,0 +1,239 @@
+module intgnssrspdmod
+!$$$ module documentation block
+!           .      .    .                                       .
+! module:   intgnssrspdmod    module for intgnssrspd and its tangent linear i ntgnssrspd_tl
+!           (Based on intspd.f90)
+!   prgmmr: Apodaca, K. - Add CYGNSS TL and AD modules
+!   email: Karina.Apodaca@Spire.com
+!
+! abstract: module for intgnssrspd and its tangent linear intgnssrspd_tl
+!
+! program history log:
+!   2023-04-21  k apodaca - initial version
+!
+! subroutines included:
+!   sub intgnssrspd_
+!
+! variable definitions:
+!
+! attributes:
+!   language: f90
+!   machine:
+!
+!$$$ end documentation block
+
+use m_obsNode, only: obsNode
+use m_gnssrspdNode, only: gnssrspdNode
+use m_gnssrspdNode, only: gnssrspdNode_typecast
+use m_gnssrspdNode, only: gnssrspdNode_nextcast
+use m_obsdiagNode, only: obsdiagNode_set
+implicit none
+
+PRIVATE
+PUBLIC intgnssrspd
+
+interface intgnssrspd; module procedure &
+          intgnssrspd_
+end interface
+
+contains
+
+subroutine intgnssrspd_(gnssrspdhead,rval,sval)
+!$$$  subprogram documentation block
+!                .      .    .                                       .
+! subprogram:    intgnssrspd      apply nonlin qc obs operator for wind speed
+!   prgmmr: K. Apodaca       org: Spire Global, Inc. date: 2023-03-15
+!   email: Karina.Apodaca@Spire.com
+!   Based on intspd.f90 by  prgmmr: derber     org: np23   date: 1991-02-26
+!
+! abstract: apply nonlinear observation operator and adjoint for winds
+!
+! program history log:
+!   2023-04-21  k apodaca - initial version`
+!
+!   input argument list:
+!     gnssrspdhead  - obs type pointer to obs structure
+!     su       - u increment in grid space
+!     sv       - v increment in grid space
+!     ru
+!     rv
+!
+!   output argument list:
+!     gnssrspdhead  - obs type pointer to obs structure
+!     ru       - u results from observation operator 
+!     rv       - v results from observation operator 
+!
+! attributes:
+!   language: f90
+!   machine:  ibm RS/6000 SP
+!
+!$$$
+  use kinds, only: r_kind,i_kind
+  use obsmod, only: lsaveobsens,l_do_adjoint,luse_obsdiag
+  use qcmod, only: nlnqc_iter,varqc_iter
+  use constants, only: zero, half, one, tiny_r_kind,cg_term,r3600
+  use gsi_4dvar, only: ltlint
+  use jfunc, only: jiter
+  use gsi_bundlemod, only: gsi_bundle
+  use gsi_bundlemod, only: gsi_bundlegetpointer
+  use gsi_4dvar, only: ladtest_obs
+  implicit none
+
+! Declare passed variables
+  class(obsNode),  pointer, intent(in   ) :: gnssrspdhead
+  type(gsi_bundle),         intent(in   ) :: sval
+  type(gsi_bundle),         intent(inout) :: rval
+
+! Declare local variables
+  integer(i_kind) ier,istatus
+  integer(i_kind) j1,j2,j3,j4
+  real(r_kind) w1,w2,w3,w4,term
+! real(r_kind) penalty
+  real(r_kind) uanl,vanl,gnssrspdanl,gnssrspd,valv,valu
+  real(r_kind) uatl,vatl,gnssrspdatl,gnssrspdtra,grad
+  real(r_kind) cg_gnssrspd,p0,wnotgross,wgross,pg_gnssrspd
+  real(r_kind),pointer,dimension(:) :: su,sv
+  real(r_kind),pointer,dimension(:) :: ru,rv
+  type(gnssrspdNode), pointer :: gnssrspdptr
+  logical :: ltlint_tmp
+
+!  If no gnssrspd data return
+  if(.not. associated(gnssrspdhead))return
+
+! Retrieve pointers
+! Simply return if any pointer not found
+  ier=0
+  call gsi_bundlegetpointer(sval,'u',su,istatus);ier=istatus+ier
+  call gsi_bundlegetpointer(sval,'v',sv,istatus);ier=istatus+ier
+  call gsi_bundlegetpointer(rval,'u',ru,istatus);ier=istatus+ier
+  call gsi_bundlegetpointer(rval,'v',rv,istatus);ier=istatus+ier
+  if(ier/=0)return
+
+  if( ladtest_obs) then
+     ltlint_tmp = ltlint   
+     ltlint = .true.
+  end if
+  !gnssrspdptr => gnssrspdhead
+  gnssrspdptr => gnssrspdNode_typecast(gnssrspdhead)
+  do while (associated(gnssrspdptr))
+
+     j1 = gnssrspdptr%ij(1)
+     j2 = gnssrspdptr%ij(2)
+     j3 = gnssrspdptr%ij(3)
+     j4 = gnssrspdptr%ij(4)
+     w1 = gnssrspdptr%wij(1)
+     w2 = gnssrspdptr%wij(2)
+     w3 = gnssrspdptr%wij(3)
+     w4 = gnssrspdptr%wij(4)
+
+
+     valu=zero
+     valv=zero
+     gnssrspdtra=sqrt(gnssrspdptr%uges*gnssrspdptr%uges+gnssrspdptr%vges*gnssrspdptr%vges)
+
+     if (ltlint) then
+
+        if (gnssrspdtra>EPSILON(gnssrspdtra)) then
+!          Forward model
+           uatl=w1*su(j1)+w2*su(j2)+w3*su(j3)+w4*su(j4)
+           vatl=w1*sv(j1)+w2*sv(j2)+w3*sv(j3)+w4*sv(j4)
+           gnssrspdatl=gnssrspdptr%uges*uatl+gnssrspdptr%vges*vatl
+           gnssrspdatl=gnssrspdatl/gnssrspdtra
+
+           if(luse_obsdiag)then
+              if (lsaveobsens) then
+                 grad=gnssrspdptr%raterr2*gnssrspdptr%err2*gnssrspdatl
+                 !-- gnssrspdptr%diags%obssen(jiter)=grad
+                 call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,obssen=grad)
+              else
+                 !-- if (gnssrspdptr%luse) gnssrspdptr%diags%tldepart(jiter)=gnssrspdatl
+                 if (gnssrspdptr%luse) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,tldepart=gnssrspdatl)
+              endif
+           endif
+
+           if (l_do_adjoint) then
+              if (.not. lsaveobsens) then
+                 if( ladtest_obs ) then
+                    grad=gnssrspdatl
+                 else
+                    gnssrspd=gnssrspdatl-gnssrspdptr%res+sqrt(gnssrspdptr%uges*gnssrspdptr%uges+gnssrspdptr%vges*gnssrspdptr%vges)
+                    grad=gnssrspdptr%raterr2*gnssrspdptr%err2*gnssrspd
+                 end if
+              endif
+
+!             Adjoint
+              valu=grad*gnssrspdptr%uges/gnssrspdtra
+              valv=grad*gnssrspdptr%vges/gnssrspdtra
+           endif
+        else
+           if(luse_obsdiag)then
+              !-- if (gnssrspdptr%luse) gnssrspdptr%diags%tldepart(jiter)=zero
+              if (gnssrspdptr%luse) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,tldepart=zero)
+              !-- if (lsaveobsens) gnssrspdptr%diags%obssen(jiter)=zero
+              if (lsaveobsens) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,obssen=zero)
+           end if
+        endif
+
+
+     else ! < ltlint >
+
+!       Forward model
+        uanl=gnssrspdptr%uges+w1* su(j1)+w2* su(j2)+w3* su(j3)+w4* su(j4)
+        vanl=gnssrspdptr%vges+w1* sv(j1)+w2* sv(j2)+w3* sv(j3)+w4* sv(j4)
+        gnssrspdanl=sqrt(uanl*uanl+vanl*vanl)
+        !-- if (luse_obsdiag .and. gnssrspdptr%luse) gnssrspdptr%diags%tldepart(jiter)=gnssrspdanl-gnssrspdtra
+        if (luse_obsdiag .and. gnssrspdptr%luse) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,tldepart=gnssrspdanl-gnssrspdtra)
+
+        if (l_do_adjoint) then
+           valu=zero
+           valv=zero
+           gnssrspd=gnssrspdanl-gnssrspdptr%res
+           grad=gnssrspdptr%raterr2*gnssrspdptr%err2*gnssrspd
+
+!          Adjoint
+!          if(gnssrspdanl > tiny_r_kind*100._r_kind) then
+           if (gnssrspdanl>EPSILON(gnssrspdanl)) then
+              !-- if (luse_obsdiag .and. lsaveobsens) gnssrspdptr%diags%obssen(jiter)=grad
+              if (luse_obsdiag .and. lsaveobsens) call obsdiagNode_set(gnssrspdptr%diags,jiter=jiter,obssen=grad)
+              valu=uanl/gnssrspdanl
+              valv=vanl/gnssrspdanl
+              if (nlnqc_iter .and. gnssrspdptr%pg > tiny_r_kind .and.  &
+                                   gnssrspdptr%b  > tiny_r_kind) then
+                 pg_gnssrspd=gnssrspdptr%pg*varqc_iter
+                 cg_gnssrspd=cg_term/gnssrspdptr%b
+                 wnotgross= one-pg_gnssrspd
+                 wgross = pg_gnssrspd*cg_gnssrspd/wnotgross
+                 p0 = wgross/(wgross+exp(-half*gnssrspdptr%err2*gnssrspd**2))
+                 term = (one-p0)
+                 grad = grad*term
+              endif
+           end if
+        endif ! < l_do_adjoint >
+
+        valu=valu*grad
+        valv=valv*grad
+
+     endif ! < ltlint >
+
+
+     if (l_do_adjoint) then
+        ru(j1)=ru(j1)+w1*valu
+        ru(j2)=ru(j2)+w2*valu
+        ru(j3)=ru(j3)+w3*valu
+        ru(j4)=ru(j4)+w4*valu
+        rv(j1)=rv(j1)+w1*valv
+        rv(j2)=rv(j2)+w2*valv
+        rv(j3)=rv(j3)+w3*valv
+        rv(j4)=rv(j4)+w4*valv
+
+     endif
+
+     !gnssrspdptr => gnssrspdptr%llpoint
+     gnssrspdptr => gnssrspdNode_nextcast(gnssrspdptr)
+
+  end do
+  if( ladtest_obs) ltlint = ltlint_tmp
+  return
+end subroutine intgnssrspd_
+
+end module intgnssrspdmod
diff --git a/src/gsi/intjo.f90 b/src/gsi/intjo.f90
index a68355471b..4ce2952062 100644
--- a/src/gsi/intjo.f90
+++ b/src/gsi/intjo.f90
@@ -14,6 +14,7 @@ module intjomod
 !   2016-08-29  J Guo   - Separated calls to intozlay() and intozlev()
 !   2018-08-10  guo     - a new generic intjo() implementation replaced type
 !                         specific intXYZ() calls with polymorphic %intjo().
+!   2022-03-15 K Apodaca - add CYGNSS and Spire ocean wind speed                        
 !
 ! subroutines included:
 !   sub intjo_
@@ -31,7 +32,7 @@ module intjomod
 use gsi_obOperTypeManager, only: &
   iobOper_t,          iobOper_pw,         iobOper_q,                                                    &
                                           iobOper_cldtot,     iobOper_w,          iobOper_dw,           &
-  iobOper_rw,         iobOper_dbz,        iobOper_fed,                                                  &
+  iobOper_rw,         iobOper_dbz,        iobOper_fed,        iobOper_gnssrspd,                         &
                       iobOper_spd,        iobOper_oz,         iobOper_o3l,        iobOper_colvk,        &
   iobOper_pm2_5,      iobOper_pm10,       iobOper_ps,         iobOper_tcp,        iobOper_sst,          &
   iobOper_gpsbend,    iobOper_gpsref,                                                                   &
@@ -60,7 +61,7 @@ module intjomod
 integer(i_kind),parameter,dimension(obOper_count):: ix_obtype = (/ &
   iobOper_t,          iobOper_pw,         iobOper_q,                                                    &
                                           iobOper_cldtot,     iobOper_w,          iobOper_dw,           &
-  iobOper_rw,         iobOper_dbz,        iobOper_fed,                                                  &
+  iobOper_rw,         iobOper_dbz,        iobOper_fed,        iobOper_gnssrspd,                         &
                       iobOper_spd,        iobOper_oz,         iobOper_o3l,        iobOper_colvk,        &
   iobOper_pm2_5,      iobOper_pm10,       iobOper_ps,         iobOper_tcp,        iobOper_sst,          &
   iobOper_gpsbend,    iobOper_gpsref,                                                                   &
diff --git a/src/gsi/m_gnssrspdNode.F90 b/src/gsi/m_gnssrspdNode.F90
new file mode 100644
index 0000000000..1a8400440f
--- /dev/null
+++ b/src/gsi/m_gnssrspdNode.F90
@@ -0,0 +1,257 @@
+module m_gnssrspdNode
+!$$$  subprogram documentation block
+!                .      .    .                                       .
+! subprogram:	 module m_gnssrspdNode
+!   prgmmr:      K. Apodaca <Karina.Apodaca@spire.com>
+!      org:      Spire Global, Inc.
+!     date:      2023-03-05
+!                (based on m_*Node.F90 modules by:
+!   prgmmr:	 j guo <jguo@nasa.gov>
+!      org:	 NASA/GSFC, Global Modeling and Assimilation Office, 610.3
+!     date:	 2016-05-18
+!
+! abstract: class-module of obs-type gnssrspdNode (wind speed)
+!
+! program history log:
+!   2016-05-18  j guo   - added this document block for the initial polymorphic
+!                         implementation.
+!
+!   input argument list: see Fortran 90 style document below
+!
+!   output argument list: see Fortran 90 style document below
+!
+! attributes:
+!   language: Fortran 90 and/or above
+!   machine:
+!
+!$$$  end subprogram documentation block
+
+! module interface:
+  use m_obsdiagNode, only: obs_diag
+  use m_obsdiagNode, only: obs_diags
+  use kinds , only: i_kind,r_kind
+  use mpeu_util, only: assert_,die,perr,warn,tell
+  use m_obsNode, only: obsNode
+  implicit none
+  private
+
+  public:: gnssrspdNode
+
+  type,extends(obsNode):: gnssrspdNode
+     !type(gnssrspd_ob_type),pointer :: llpoint => NULL()
+     type(obs_diag), pointer :: diags => NULL()
+     real(r_kind)    :: res           !  speed observation
+     real(r_kind)    :: err2          !  speed error squared
+     real(r_kind)    :: raterr2       !  square of ratio of final obs error 
+                                      !  to original obs error
+     !real(r_kind)    :: time          !  observation time in sec     
+     real(r_kind)    :: b             !  variational quality control parameter
+     real(r_kind)    :: pg            !  variational quality control parameter
+     real(r_kind)    :: wij(4)        !  horizontal interpolation weights
+     real(r_kind)    :: uges          !  guess u value        
+     real(r_kind)    :: vges          !  guess v value        
+     integer(i_kind) :: ij(4)         !  horizontal locations
+     !logical         :: luse          !  flag indicating if ob is used in pen.
+
+     !integer(i_kind) :: idv,iob	      ! device id and obs index for sorting
+     !real   (r_kind) :: elat, elon      ! earth lat-lon for redistribution
+     !real   (r_kind) :: dlat, dlon      ! earth lat-lon for redistribution
+     real   (r_kind) :: factw           ! factor of 10m wind
+  contains
+    procedure,nopass::  mytype
+    procedure::  setHop => obsNode_setHop_
+    procedure::   xread => obsNode_xread_
+    procedure::  xwrite => obsNode_xwrite_
+    procedure:: isvalid => obsNode_isvalid_
+    procedure::  gettlddp => gettlddp_
+
+    ! procedure, nopass:: headerRead  => obsHeader_read_
+    ! procedure, nopass:: headerWrite => obsHeader_write_
+    ! procedure:: init  => obsNode_init_
+    ! procedure:: clean => obsNode_clean_
+  end type gnssrspdNode
+
+  public:: gnssrspdNode_typecast
+  public:: gnssrspdNode_nextcast
+        interface gnssrspdNode_typecast; module procedure typecast_ ; end interface
+        interface gnssrspdNode_nextcast; module procedure nextcast_ ; end interface
+
+  public:: gnssrspdNode_appendto
+        interface gnssrspdNode_appendto; module procedure appendto_ ; end interface
+
+  character(len=*),parameter:: MYNAME="m_gnssrspdNode"
+
+#include "myassert.H"
+#include "mytrace.H"
+contains
+function typecast_(aNode) result(ptr_)
+!-- cast a class(obsNode) to a type(gnssrspdNode)
+  use m_obsNode, only: obsNode
+  implicit none
+  type(gnssrspdNode ),pointer:: ptr_
+  class(obsNode),pointer,intent(in):: aNode
+  ptr_ => null()
+  if(.not.associated(aNode)) return
+        ! logically, typecast of a null-reference is a null pointer.
+  select type(aNode)
+  type is(gnssrspdNode)
+    ptr_ => aNode
+  end select
+return
+end function typecast_
+
+function nextcast_(aNode) result(ptr_)
+!-- cast an obsNode_next(obsNode) to a type(gnssrspdNode)
+  use m_obsNode, only: obsNode,obsNode_next
+  implicit none
+  type(gnssrspdNode ),pointer:: ptr_
+  class(obsNode),target ,intent(in):: aNode
+
+  class(obsNode),pointer:: inode_
+  inode_ => obsNode_next(aNode)
+  ptr_ => typecast_(inode_)
+return
+end function nextcast_
+
+subroutine appendto_(aNode,oll)
+!-- append aNode to linked-list oLL
+  use m_obsNode , only: obsNode
+  use m_obsLList, only: obsLList,obsLList_appendNode
+  implicit none
+  type(gnssrspdNode),pointer,intent(in):: aNode
+  type(obsLList),intent(inout):: oLL
+
+  class(obsNode),pointer:: inode_
+  inode_ => aNode
+  call obsLList_appendNode(oLL,inode_)
+  inode_ => null()
+end subroutine appendto_
+
+! obsNode implementations
+
+function mytype()
+  implicit none
+  character(len=:),allocatable:: mytype
+  mytype="[gnssrspdNode]"
+end function mytype
+
+subroutine obsNode_xread_(aNode,iunit,istat,diagLookup,skip)
+  use m_obsdiagNode, only: obsdiagLookup_locate
+  implicit none
+  class(gnssrspdNode),intent(inout):: aNode
+  integer(i_kind),intent(in   ):: iunit
+  integer(i_kind),intent(  out):: istat
+  type(obs_diags),intent(in   ):: diagLookup
+  logical,optional,intent(in   ):: skip
+
+  character(len=*),parameter:: myname_=MYNAME//'.obsNode_xread_'
+  logical:: skip_
+_ENTRY_(myname_)
+  skip_=.false.
+  if(present(skip)) skip_=skip
+
+  istat=0
+  if(skip_) then
+    read(iunit,iostat=istat)
+                if(istat/=0) then
+                  call perr(myname_,'skipping read(%(res,err2,...)), iostat =',istat)
+                  _EXIT_(myname_)
+                  return
+                endif
+
+  else
+    read(iunit,iostat=istat)    aNode%res    , &
+                                aNode%err2   , &
+                                aNode%raterr2, &
+                                aNode%b      , &
+                                aNode%pg     , &
+                                aNode%uges   , &
+                                aNode%vges   , &
+                                aNode%factw  , &
+                                aNode%wij    , &
+                                aNode%ij
+                if (istat/=0) then
+                  call perr(myname_,'read(%(res,err2,...)), iostat =',istat)
+                  _EXIT_(myname_)
+                  return
+                end if
+
+    aNode%diags => obsdiagLookup_locate(diagLookup,aNode%idv,aNode%iob,1_i_kind)
+                if(.not.associated(aNode%diags)) then
+                  call perr(myname_,'obsdiagLookup_locate(), %idv =',aNode%idv)
+                  call perr(myname_,'                        %iob =',aNode%iob)
+                  call  die(myname_)
+                endif
+  endif
+_EXIT_(myname_)
+return
+end subroutine obsNode_xread_
+
+subroutine obsNode_xwrite_(aNode,junit,jstat)
+  implicit none
+  class(gnssrspdNode),intent(in):: aNode
+  integer(i_kind),intent(in   ):: junit
+  integer(i_kind),intent(  out):: jstat
+
+  character(len=*),parameter:: myname_=MYNAME//'.obsNode_xwrite_'
+_ENTRY_(myname_)
+
+  jstat=0
+  write(junit,iostat=jstat)     aNode%res    , &
+                                aNode%err2   , &
+                                aNode%raterr2, &
+                                aNode%b      , &
+                                aNode%pg     , &
+                                aNode%uges   , &
+                                aNode%vges   , &
+                                aNode%factw  , &
+                                aNode%wij    , &
+                                aNode%ij
+                if (jstat/=0) then
+                  call perr(myname_,'write(%(res,err2,...)), iostat =',jstat)
+                  _EXIT_(myname_)
+                  return
+                end if
+_EXIT_(myname_)
+return
+end subroutine obsNode_xwrite_
+
+subroutine obsNode_setHop_(aNode)
+  use m_cvgridLookup, only: cvgridLookup_getiw
+  implicit none
+  class(gnssrspdNode),intent(inout):: aNode
+
+  character(len=*),parameter:: myname_=MYNAME//'::obsNode_setHop_'
+_ENTRY_(myname_)
+  call cvgridLookup_getiw(aNode%elat,aNode%elon,aNode%ij,aNode%wij)
+  aNode%wij(1:4) = aNode%wij(1:4)*aNode%factw
+_EXIT_(myname_)
+return
+end subroutine obsNode_setHop_
+
+function obsNode_isvalid_(aNode) result(isvalid_)
+  implicit none
+  logical:: isvalid_
+  class(gnssrspdNode),intent(in):: aNode
+
+  character(len=*),parameter:: myname_=MYNAME//'::obsNode_isvalid_'
+_ENTRY_(myname_)
+  isvalid_=associated(aNode%diags)
+_EXIT_(myname_)
+return
+end function obsNode_isvalid_
+
+pure subroutine gettlddp_(aNode,jiter,tlddp,nob)
+  use kinds, only: r_kind
+  implicit none
+  class(gnssrspdNode), intent(in):: aNode
+  integer(kind=i_kind),intent(in):: jiter
+  real(kind=r_kind),intent(inout):: tlddp
+  integer(kind=i_kind),optional,intent(inout):: nob
+
+  tlddp = tlddp + aNode%diags%tldepart(jiter)*aNode%diags%tldepart(jiter)
+  if(present(nob)) nob=nob+1
+return
+end subroutine gettlddp_
+
+end module m_gnssrspdNode
diff --git a/src/gsi/m_obsNodeTypeManager.F90 b/src/gsi/m_obsNodeTypeManager.F90
index 43b42e4bf2..9da356c808 100644
--- a/src/gsi/m_obsNodeTypeManager.F90
+++ b/src/gsi/m_obsNodeTypeManager.F90
@@ -15,6 +15,8 @@ module m_obsNodeTypeManager
 !   2018-01-23  k apodaca - add a new observation type i.e. lightning (light)
 !                           suitable for the GOES/GLM instrument
 !
+!   2024-04-07 k apodaca - add gnssrspdnode
+!
 !   input argument list: see Fortran 90 style document below
 !
 !   output argument list: see Fortran 90 style document below
@@ -54,6 +56,8 @@ module m_obsNodeTypeManager
   use m_wspd10mNode, only: wspd10mNode
   use m_uwnd10mNode, only: uwnd10mNode
   use m_vwnd10mNode, only: vwnd10mNode
+  
+  use m_gnssrspdNode, only: gnssrspdNode
 
   use m_td2mNode , only:  td2mNode
   use m_mxtmNode , only:  mxtmNode
@@ -122,10 +126,10 @@ module m_obsNodeTypeManager
   public:: iobsNode_cldch
   public:: iobsNode_swcp
   public:: iobsNode_lwcp
-
   public:: iobsNode_light
   public:: iobsNode_dbz
   public:: iobsNode_fed
+  public:: iobsNode_gnssrspd
 
   public :: obsNode_typeMold
   public :: obsNode_typeIndex
@@ -166,6 +170,8 @@ module m_obsNodeTypeManager
   type(wspd10mNode), target, save:: wspd10m_mold
   type(uwnd10mNode), target, save:: uwnd10m_mold
   type(vwnd10mNode), target, save:: vwnd10m_mold
+  
+  type(gnssrspdNode), target, save:: gnssrspd_mold
 
   type(   td2mNode), target, save::    td2m_mold
   type(   mxtmNode), target, save::    mxtm_mold
@@ -214,6 +220,7 @@ module m_obsNodeTypeManager
     enumerator:: iobsNode_w
     enumerator:: iobsNode_q
     enumerator:: iobsNode_spd
+    enumerator:: iobsNode_gnssrspd
     enumerator:: iobsNode_rw
     enumerator:: iobsNode_dw
     enumerator:: iobsNode_sst
@@ -301,6 +308,8 @@ function vname2index_(vname) result(index_)
              "[uwnd10mnode]"); index_ = iobsNode_uwnd10m
   case("vwnd10m", &
              "[vwnd10mnode]"); index_ = iobsNode_vwnd10m
+  case("gnssrspd", &
+             "[gnssrspdnode]"); index_ = iobsNode_gnssrspd
 
   case("td2m" , "[td2mnode]"); index_ = iobsNode_td2m
   case("mxtm" , "[mxtmnode]"); index_ = iobsNode_mxtm
@@ -365,6 +374,8 @@ function vmold2index_select_(mold) result(index_)
   type is(wspd10mNode); index_ = iobsNode_wspd10m
   type is(uwnd10mNode); index_ = iobsNode_uwnd10m
   type is(vwnd10mNode); index_ = iobsNode_vwnd10m
+  
+  type is(gnssrspdNode); index_ = iobsNode_gnssrspd
 
   type is( td2mNode); index_ = iobsNode_td2m
   type is( mxtmNode); index_ = iobsNode_mxtm
@@ -423,6 +434,8 @@ function index2vmold_(i_obType) result(obsmold_)
   case(iobsNode_wspd10m); obsmold_ => wspd10m_mold
   case(iobsNode_uwnd10m); obsmold_ => uwnd10m_mold
   case(iobsNode_vwnd10m); obsmold_ => vwnd10m_mold
+  
+  case(iobsNode_gnssrspd  ); obsmold_ => gnssrspd_mold
 
   case(iobsNode_td2m ); obsmold_ =>    td2m_mold
   case(iobsNode_mxtm ); obsmold_ =>    mxtm_mold
diff --git a/src/gsi/m_rhs.F90 b/src/gsi/m_rhs.F90
index aea417fe27..723b0a6918 100644
--- a/src/gsi/m_rhs.F90
+++ b/src/gsi/m_rhs.F90
@@ -17,6 +17,7 @@ module m_rhs
 !                         through an enum block.
 !                       - removed external dimension argument aworkdim2 of
 !                         rhs_alloc().
+!   2022-03-15  K. Apodaca - add GNSS-R L2 ocean wind speed index (i_gnssrspd)
 !
 !   input argument list: see Fortran 90 style document below
 !
@@ -66,6 +67,7 @@ module m_rhs
   public:: i_vis
   public:: i_pblh
   public:: i_wspd10m
+  public:: i_gnssrspd
   public:: i_td2m
   public:: i_mxtm
   public:: i_mitm
@@ -133,6 +135,7 @@ module m_rhs
     enumerator:: i_vis
     enumerator:: i_pblh
     enumerator:: i_wspd10m
+    enumerator:: i_gnssrspd
     enumerator:: i_td2m
     enumerator:: i_mxtm
     enumerator:: i_mitm
@@ -207,7 +210,6 @@ subroutine rhs_alloc()
   allocate(rhs_stats_oz(9,jpch_oz))
 
   allocate(rhs_toss_gps(max(1,nprof_gps)))
-
   rhs_awork    =zero
   rhs_bwork    =zero
   rhs_aivals   =zero
diff --git a/src/gsi/mrmsmod.f90 b/src/gsi/mrmsmod.f90
index 713b974d7a..1cbad54904 100644
--- a/src/gsi/mrmsmod.f90
+++ b/src/gsi/mrmsmod.f90
@@ -27,7 +27,9 @@ subroutine load_mrms_data_info (mrms_listfile,nrows0,ntot_mrms,nrows_mrms,nrows,
    integer(i_kind),parameter::nobstype_mrms=2 ! first for vr, and the second for ref
    integer(i_kind),intent(in):: nrows0,ntot_mrms,nrows_mrms,nrows
    character(len=*),intent(in),optional ::mrms_listfile, rcname ! input filename
-   character(10),intent(inout),dimension(nrows):: dtype,ditype,dplat
+   character(len=*),intent(inout),dimension(nrows):: dtype
+   character(len=*),intent(inout),dimension(nrows):: ditype
+   character(len=*),intent(inout),dimension(nrows):: dplat
    character(20),intent(inout),dimension(nrows):: obsfile_all
    character(*),intent(inout),dimension(nrows):: dfile
    character(20),intent(inout),dimension(nrows):: dsis
@@ -42,7 +44,8 @@ subroutine load_mrms_data_info (mrms_listfile,nrows0,ntot_mrms,nrows_mrms,nrows,
 
 
    real(r_kind),allocatable,dimension(:):: dmesh_mrms
-   character(10),allocatable,dimension(:):: dtype_mrms,ditype_mrms,dplat_mrms
+   character(10),allocatable,dimension(:):: dtype_mrms,ditype_mrms
+   character(11),allocatable,dimension(:):: dplat_mrms
    character(120),allocatable,dimension(:):: dfile_mrms
    character(20),allocatable,dimension(:):: dsis_mrms
    real(r_kind) ,allocatable,dimension(:):: dval_mrms
diff --git a/src/gsi/netcdfgfs_io.f90 b/src/gsi/netcdfgfs_io.f90
index 8aa6c4cc6a..7db02636cd 100644
--- a/src/gsi/netcdfgfs_io.f90
+++ b/src/gsi/netcdfgfs_io.f90
@@ -159,17 +159,17 @@ subroutine read_
     real(r_kind),pointer,dimension(:,:,:):: ges_qs_it  => NULL()
     real(r_kind),pointer,dimension(:,:,:):: ges_qg_it  => NULL()
     real(r_kind),pointer,dimension(:,:,:):: ges_cf_it  => NULL()
-
+    real(r_kind),pointer,dimension(:,:,:):: ges_ni_it  => NULL()
+    real(r_kind),pointer,dimension(:,:,:):: ges_nr_it  => NULL()
     type(sub2grid_info) :: grd_t
     logical regional
     logical:: l_cld_derived,zflag,inithead
-
     type(gsi_bundle) :: atm_bundle
     type(gsi_grid)   :: atm_grid
     integer(i_kind),parameter :: n2d=2
     ! integer(i_kind),parameter :: n3d=8
     integer(i_kind),parameter :: n2d_2m=4
-    integer(i_kind),parameter :: n3d=14
+    integer(i_kind),parameter :: n3d=16
     character(len=4), parameter :: vars2d(n2d) = (/ 'z   ', 'ps  ' /)
     character(len=4), parameter :: vars2d_with2m(n2d_2m) = (/ 'z   ', 'ps  ','t2m ','q2m ' /)
     ! character(len=4), parameter :: vars3d(n3d) = (/ 'u   ', 'v   ', &
@@ -182,13 +182,16 @@ subroutine read_
                                                     'cw  ', 'oz  ', &
                                                     'ql  ', 'qi  ', &
                                                     'qr  ', 'qs  ', &
-                                                    'qg  ', 'cf  ' /)
+                                                    'qg  ', 'ni  ', &
+                                                    'nr  ', 'cf  ' /)
+
     real(r_kind),pointer,dimension(:,:):: ptr2d   =>NULL()
     real(r_kind),pointer,dimension(:,:,:):: ptr3d =>NULL()
 
     regional=.false.
     inner_vars=1
-    num_fields=min(14*grd_a%nsig+2,npe)
+
+    num_fields=min(n3d*grd_a%nsig+2,npe)
 !  Create temporary communication information fore read routines
     call general_sub2grid_create_info(grd_t,inner_vars,grd_a%nlat,grd_a%nlon, &
           grd_a%nsig,num_fields,regional)
@@ -200,6 +203,7 @@ subroutine read_
     else
         call gsi_bundlecreate(atm_bundle,atm_grid,'aux-atm-read',istatus,names2d=vars2d,names3d=vars3d)
     endif
+
     if(istatus/=0) then
       write(6,*) myname_,': trouble creating atm_bundle'
       call stop2(999)
@@ -248,6 +252,8 @@ subroutine read_
                    if (associated(ges_qr_it)) ges_qr_it(i,j,k) = max(qcmin,ges_qr_it(i,j,k))
                    if (associated(ges_qs_it)) ges_qs_it(i,j,k) = max(qcmin,ges_qs_it(i,j,k))
                    if (associated(ges_qg_it)) ges_qg_it(i,j,k) = max(qcmin,ges_qg_it(i,j,k))
+                   if (associated(ges_ni_it)) ges_ni_it(i,j,k) = max(qcmin,ges_ni_it(i,j,k))
+                   if (associated(ges_nr_it)) ges_nr_it(i,j,k) = max(qcmin,ges_nr_it(i,j,k))
                    if (associated(ges_cf_it)) ges_cf_it(i,j,k) = min(max(zero,ges_cf_it(i,j,k)),one)
                 enddo
              enddo
@@ -256,6 +262,8 @@ subroutine read_
           l_cld_derived = associated(ges_cwmr_it).and.&
                           associated(ges_q_it)   .and.&
                           associated(ges_ql_it)  .and.&
+                          associated(ges_ni_it)  .and.&
+                          associated(ges_nr_it)  .and.&
                           associated(ges_qi_it)  .and.& 
                           associated(ges_tv_it)
 !         call set_cloud_lower_bound(ges_cwmr_it)
@@ -343,6 +351,16 @@ subroutine set_guess_
        call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ql',ges_ql_it,istatus)
        if(istatus==0) ges_ql_it = ptr3d
     endif
+    call gsi_bundlegetpointer (atm_bundle,'ni',ptr3d,istatus)
+    if (istatus==0) then
+       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'ni',ges_ni_it,istatus)
+       if(istatus==0) ges_ni_it = ptr3d
+    endif
+    call gsi_bundlegetpointer (atm_bundle,'nr',ptr3d,istatus)
+    if (istatus==0) then
+       call gsi_bundlegetpointer (gsi_metguess_bundle(it),'nr',ges_nr_it,istatus)
+       if(istatus==0) ges_nr_it = ptr3d
+    endif
     call gsi_bundlegetpointer (atm_bundle,'qi',ptr3d,istatus)
     if (istatus==0) then
        call gsi_bundlegetpointer (gsi_metguess_bundle(it),'qi',ges_qi_it,istatus)
diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90
index 1c45c62bc8..eab9c828f7 100644
--- a/src/gsi/obsmod.F90
+++ b/src/gsi/obsmod.F90
@@ -160,6 +160,7 @@ module obsmod
 !   2021-11-16 Zhao      - add option l_obsprvdiag (if true) to trigger the output of
 !                          observation provider and sub-provider information into
 !                          obsdiags files (used for AutoObsQC)
+!  2022-03-15  K. Apodaca - add GNSS-R L2 ocean wind speed observations (CYGNSS, Spire) 
 !   2023-07-10  Y. Wang, D. Dowell - add variables for flash extent density
 !   2023-10-10  H. Wang (GSL) - add variables for flash extent density EnVar DA
 ! 
@@ -463,8 +464,8 @@ module obsmod
   public :: iout_lag,iout_uv,iout_gps,iout_ps,iout_light,mype_light
   public :: mype_gust,mype_vis,mype_pblh,iout_gust,iout_vis,iout_pblh
   public :: mype_tcamt,mype_lcbas,iout_tcamt,iout_lcbas
-  public :: mype_wspd10m,mype_td2m,iout_wspd10m,iout_td2m
-  public :: mype_uwnd10m,mype_vwnd10m,iout_uwnd10m,iout_vwnd10m
+  public :: mype_wspd10m,mype_gnssrspd,mype_td2m,iout_wspd10m,iout_gnssrspd,iout_td2m
+  public :: mype_uwnd10m,mype_vwnd10m,iout_uwnd10m,iout_vwnd10m 
   public :: mype_mxtm,mype_mitm,iout_mxtm,iout_mitm
   public :: mype_pmsl,mype_howv,iout_pmsl,iout_howv
   public :: mype_swcp,mype_lwcp,iout_swcp,iout_lwcp
@@ -600,12 +601,12 @@ module obsmod
   integer(i_kind) iout_dw,iout_gps,iout_sst,iout_tcp,iout_lag
   integer(i_kind) iout_co,iout_gust,iout_vis,iout_pblh,iout_tcamt,iout_lcbas
   integer(i_kind) iout_cldch
-  integer(i_kind) iout_wspd10m,iout_td2m,iout_mxtm,iout_mitm,iout_pmsl,iout_howv
+  integer(i_kind) iout_wspd10m,iout_gnssrspd,iout_td2m,iout_mxtm,iout_mitm,iout_pmsl,iout_howv
   integer(i_kind) iout_uwnd10m,iout_vwnd10m,iout_fed
   integer(i_kind) mype_t,mype_q,mype_uv,mype_ps,mype_pw, &
                   mype_rw,mype_dw,mype_gps,mype_sst, &
                   mype_tcp,mype_lag,mype_co,mype_gust,mype_vis,mype_pblh, &
-                  mype_wspd10m,mype_td2m,mype_mxtm,mype_mitm,mype_pmsl,mype_howv,&
+                  mype_wspd10m,mype_gnssrspd,mype_td2m,mype_mxtm,mype_mitm,mype_pmsl,mype_howv,&
                   mype_uwnd10m,mype_vwnd10m, mype_tcamt,mype_lcbas, mype_dbz, mype_fed
   integer(i_kind) mype_cldch
   integer(i_kind) iout_swcp, iout_lwcp
@@ -625,7 +626,8 @@ module obsmod
   character(128) dirname
   character(128) obs_input_common
   character(20),allocatable,dimension(:):: obsfile_all
-  character(10),allocatable,dimension(:):: dtype,ditype,dplat
+  character(10),allocatable,dimension(:):: dtype,ditype
+  character(11),allocatable,dimension(:):: dplat
   character(120),allocatable,dimension(:):: dfile
   character(20),allocatable,dimension(:):: dsis
   real(r_kind) ,allocatable,dimension(:):: dval
@@ -887,6 +889,7 @@ subroutine init_obsmod_dflts
     iout_light=237 ! lightning
     iout_dbz=238 ! radar reflectivity
     iout_fed=239   ! flash extent density
+    iout_gnssrspd=240 ! GNSS-R wind speed
 
     mype_ps = npe-1          ! surface pressure
     mype_t  = max(0,npe-2)   ! temperature
@@ -922,6 +925,7 @@ subroutine init_obsmod_dflts
     mype_light=max(0,npe-32)! GOES/GLM lightning
     mype_dbz=max(0,npe-33)   ! radar reflectivity
     mype_fed= max(0,npe-34)  ! flash extent density
+    mype_gnssrspd= max(0,npe-35) ! surface GNSS-R speed
 
 
 !   Initialize arrays used in namelist obs_input 
diff --git a/src/gsi/prt_guess.f90 b/src/gsi/prt_guess.f90
index 97c1e0a7c1..dbce7e0f79 100644
--- a/src/gsi/prt_guess.f90
+++ b/src/gsi/prt_guess.f90
@@ -287,7 +287,8 @@ subroutine prt_guess2(sgrep)
   integer(i_kind), parameter :: nvars3=5
   integer(i_kind) :: nvars,nvars2,nvarsc,nc
   integer(i_kind) ii,istatus,ier,ivar
-  integer(i_kind) iql,iqi,iqr,iqs,iqg,icf    
+  integer(i_kind) iql,iqi,iqr,iqs,iqg,icf
+  integer(i_kind) ini,inr    
   integer(i_kind) ntsig
   integer(i_kind) ntsfc
   integer(i_kind) n_actual_clouds
@@ -307,6 +308,8 @@ subroutine prt_guess2(sgrep)
   real(r_kind),pointer,dimension(:,:,:)::ges_qr_it => NULL()
   real(r_kind),pointer,dimension(:,:,:)::ges_qs_it => NULL()
   real(r_kind),pointer,dimension(:,:,:)::ges_qg_it => NULL()
+  real(r_kind),pointer,dimension(:,:,:)::ges_ni_it => NULL()
+  real(r_kind),pointer,dimension(:,:,:)::ges_nr_it => NULL()
   real(r_kind),pointer,dimension(:,:,:)::ges_cf_it => NULL()  
 ! character(len=4) :: cvar(nvars+3)
   character(len=4),allocatable,dimension(:) :: cvar
@@ -338,6 +341,7 @@ subroutine prt_guess2(sgrep)
 ! get pointer to cloud water condensate
   ier=0;nvarsc=0
   iql=0;iqi=0;iqr=0;iqs=0;iqg=0
+  ini=0;inr=0
   call gsi_metguess_get('clouds::3d',n_actual_clouds,ier)
   if (mype==0) write(6,*)'prt_guess2: n_actual_clouds = ', n_actual_clouds
   if (n_actual_clouds>0) then
@@ -345,36 +349,52 @@ subroutine prt_guess2(sgrep)
      if (ivar > 0) then 
         call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'ql',ges_ql_it,istatus)
         ier=ier+istatus
+        if (ier==0) nvarsc=nvarsc+1
      endif
      call gsi_metguess_get ( 'var::qi', ivar, ier ); iqi=ivar
      if (ivar > 0) then
         call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'qi',ges_qi_it,istatus)
         ier=ier+istatus
+        if (ier==0) nvarsc=nvarsc+1
      endif
      call gsi_metguess_get ( 'var::qr', ivar, ier ); iqr=ivar
      if (ivar > 0) then
         call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'qr',ges_qr_it,istatus)
         ier=ier+istatus
+        if (ier==0) nvarsc=nvarsc+1
      endif
      call gsi_metguess_get ( 'var::qs', ivar, ier ); iqs=ivar
      if (ivar > 0) then
         call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'qs',ges_qs_it,istatus)
         ier=ier+istatus
+        if (ier==0) nvarsc=nvarsc+1
      endif
      call gsi_metguess_get ( 'var::qg', ivar, ier ); iqg=ivar
      if (ivar > 0) then
         call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'qg',ges_qg_it,istatus)
         ier=ier+istatus
+        if (ier==0) nvarsc=nvarsc+1
+     endif
+     call gsi_metguess_get ( 'var::ni', ivar, ier ); ini=ivar
+     if (ivar > 0) then
+        call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'ni',ges_ni_it,istatus)
+        ier=ier+istatus
+        if (ier==0) nvarsc=nvarsc+1
+     endif
+     call gsi_metguess_get ( 'var::nr', ivar, ier ); inr=ivar
+     if (ivar > 0) then
+        call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'nr',ges_nr_it,istatus)
+        ier=ier+istatus
+        if (ier==0) nvarsc=nvarsc+1
      endif
   end if
-  nvarsc=n_actual_clouds
   call gsi_metguess_get ( 'var::cf', ivar, ier ); icf=ivar
   if (ivar > 0) then
      call gsi_bundlegetpointer (gsi_metguess_bundle(ntsig),'cf',ges_cf_it,istatus)
      ier=ier+istatus
      nvarsc=nvarsc+1
   endif
-! if (ier/=0) return ! this is a fundamental routine, when some not found just ! emily
+! if (ier/=0) return ! this is a fundamental routine, when some not found just !
 
   nvars  = nvars1+nvarsc+nvars3
   nvars2 = nvars1+nvarsc
@@ -415,6 +435,14 @@ subroutine prt_guess2(sgrep)
      nc=nc+1
      cvar(nvars1+nc)='QG  '
   endif
+  if(ini>0) then
+     nc=nc+1
+     cvar(nvars1+nc)='NI  '
+  endif
+  if(inr>0) then
+     nc=nc+1
+     cvar(nvars1+nc)='NR  '
+  endif
   if(icf>0) then
      nc=nc+1
      cvar(nvars1+nc)='CF  '
@@ -455,6 +483,17 @@ subroutine prt_guess2(sgrep)
      nc=nc+1
      zloc(nvars1+nc) = sum   (ges_qg_it (2:lat1+1,2:lon1+1,1:nsig))
   endif
+
+  if(ini>0) then
+     nc=nc+1
+     zloc(nvars1+nc) = sum   (ges_ni_it (2:lat1+1,2:lon1+1,1:nsig))
+  endif
+
+  if(inr>0) then
+     nc=nc+1
+     zloc(nvars1+nc) = sum   (ges_nr_it (2:lat1+1,2:lon1+1,1:nsig))
+  endif
+
   if(icf>0) then 
      nc=nc+1
      zloc(nvars1+nc) = sum   (ges_cf_it (2:lat1+1,2:lon1+1,1:nsig))
@@ -491,6 +530,16 @@ subroutine prt_guess2(sgrep)
      nc=nc+1
      zloc(nvars+nvars1+nc)   = minval(ges_qg_it (2:lat1+1,2:lon1+1,1:nsig))
   endif
+
+  if(ini>0) then
+     nc=nc+1
+     zloc(nvars+nvars1+nc)   = minval(ges_ni_it (2:lat1+1,2:lon1+1,1:nsig))
+  endif
+  if(inr>0) then
+     nc=nc+1
+     zloc(nvars+nvars1+nc)   = minval(ges_nr_it (2:lat1+1,2:lon1+1,1:nsig))
+  endif
+
   if(icf>0) then 
      nc=nc+1
      zloc(nvars+nvars1+nc)   = minval(ges_cf_it (2:lat1+1,2:lon1+1,1:nsig))
@@ -527,6 +576,15 @@ subroutine prt_guess2(sgrep)
      nc=nc+1
      zloc(2*nvars+nvars1+nc) = maxval(ges_qg_it (2:lat1+1,2:lon1+1,1:nsig))
   endif
+  if(ini>0) then
+     nc=nc+1
+     zloc(2*nvars+nvars1+nc) = maxval(ges_ni_it (2:lat1+1,2:lon1+1,1:nsig))
+  endif
+  if(inr>0) then
+     nc=nc+1
+     zloc(2*nvars+nvars1+nc) = maxval(ges_nr_it (2:lat1+1,2:lon1+1,1:nsig))
+  endif
+
   if(icf>0) then 
      nc=nc+1
      zloc(2*nvars+nvars1+nc) = maxval(ges_cf_it (2:lat1+1,2:lon1+1,1:nsig))
diff --git a/src/gsi/qcmod.f90 b/src/gsi/qcmod.f90
index 9804965573..2c4cbf6317 100644
--- a/src/gsi/qcmod.f90
+++ b/src/gsi/qcmod.f90
@@ -4240,7 +4240,7 @@ subroutine qc_goesimg(nchanl,is,ndat,nsig,ich,dplat,sea,land,ice,snow,luse,   &
   real(r_kind),dimension(nsig,nchanl),intent(in ) :: temp,wmix
   real(r_kind),dimension(nchanl),   intent(in   ) :: tb_obs,tb_obs_sdv,tbc,tnoise,emissivity_k,ts
   real(r_kind),dimension(nchanl),   intent(inout) :: errf,varinv
-  character(10),                    intent(in   ) :: dplat
+  character(len=*),                 intent(in   ) :: dplat
 
 ! Declare local parameters
 
diff --git a/src/gsi/radinfo.f90 b/src/gsi/radinfo.f90
index 4ad17626e6..8bb496c645 100644
--- a/src/gsi/radinfo.f90
+++ b/src/gsi/radinfo.f90
@@ -805,10 +805,10 @@ subroutine radinfo_read
     end do
     close(lunin)
 100 format(a1,a120)
-110 format(i4,1x,a20,' chan= ',i4,  &
+110 format(i4,1x,a20,' chan= ',i5,  &
           ' var= ',f7.3,' varch_cld=',f7.3,' use= ',i2,' ermax= ',F7.3, &
           ' b_rad= ',F7.2,' pg_rad=',F7.2,' icld_det=',I2,' icloud=',I2,' iaeros=',I2)
-111 format(i4,1x,a20,' chan= ',i4,  &
+111 format(i4,1x,a20,' chan= ',i5,  &
           ' var= ',f7.3,' varch_cld=',f7.3,' use= ',i2,' ermax= ',F7.3, &
           ' b_rad= ',F7.2,' pg_rad=',F7.2,' iextra_det=',I2, 'icloud=',I2,'iaeros=', I2)
 
diff --git a/src/gsi/read_cris.f90 b/src/gsi/read_cris.f90
index d5668f8864..84288a7f04 100644
--- a/src/gsi/read_cris.f90
+++ b/src/gsi/read_cris.f90
@@ -112,9 +112,9 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,&
   integer(i_kind)  ,intent(in   ) :: mype_sub
   integer(i_kind)  ,intent(in   ) :: npe_sub
   integer(i_kind)  ,intent(in   ) :: mpi_comm_sub  
-  character(len=*), intent(in   ) :: infile
-  character(len=10),intent(in   ) :: jsatid
-  character(len=*), intent(in   ) :: obstype
+  character(len=*) ,intent(in   ) :: infile
+  character(len=*) ,intent(in   ) :: jsatid
+  character(len=*) ,intent(in   ) :: obstype
   character(len=20),intent(in   ) :: sis
   real(r_kind)     ,intent(in   ) :: twind
   real(r_kind)     ,intent(inout) :: val_cris
@@ -200,7 +200,6 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,&
   real(r_kind),dimension(83,7) :: imager_info
   real(r_kind),dimension(7)    :: imager_cluster_size
   real(r_kind),dimension(2)    :: imager_mean, imager_std_dev, imager_conversion
-  real(r_kind)                 :: imager_cluster_tot
 
 ! bufr error codes
 !   real(r_kind),dimension(7,3)  :: error_codes
@@ -848,10 +847,9 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,&
            if ( cris_cads ) then
              call ufbseq(lnbufr,imager_info,83,7,iret,'CRISCS')
              if ( iret == 7 .and. imager_info(3,1) <= 100.0_r_kind .and. &
-                  imager_info(3,1) >= zero .and. imager_coeff ) then   ! if imager cluster info exists
+                  sum(imager_info(3,:)) > zero .and. imager_coeff ) then   ! if imager cluster info exists
                imager_mean = zero
                imager_std_dev = zero
-               imager_cluster_tot = zero
                imager_cluster_flag = .TRUE.
                imager_cluster_size = imager_info(3,1:7)
                imager_cluster_size(:) = imager_cluster_size(:) / sum(imager_cluster_size(:))
@@ -871,51 +869,62 @@ subroutine read_cris(mype,val_cris,ithin,isfcalc,rmesh,jsatid,gstime,&
                imager_cluster_info: do j=1,7
                  i = imager_cluster_index(j)
 
-                 data_all(maxinfo+j,itx) =  imager_cluster_size(i)                ! Imager cluster fraction
-                 imager_cluster_tot = imager_cluster_tot + imager_info(3,i)
+!   If the cluster size, or radiance values of channel 4 or 5 are zero, do not compute statistics for that cluster
 
-                 iexponent = -(nint(imager_info(75,i)) -11)                        ! channel 15 radiance for each cluster
-                 imager_info(76,i) =  imager_info(76,i) * imager_conversion(1) * (ten ** iexponent) 
+                 if ( imager_cluster_size(i) > zero .and. imager_info(76,i) > zero .and. imager_info(81,i) > zero ) then
+                   data_all(maxinfo+j,itx) =  imager_cluster_size(i)                  ! Imager cluster fraction
 
-                 iexponent = -(nint(imager_info(77,i)) -11)                        ! channel 15 radiance std dev for each cluster.
-                 imager_info(78,i) =  imager_info(78,i) * imager_conversion(1) * (ten ** iexponent) 
+                   iexponent = -(nint(imager_info(75,i)) -11)                         ! channel 15 radiance for each cluster
+                   imager_info(76,i) =  imager_info(76,i) * imager_conversion(1) * (ten ** iexponent) 
 
-                 call crtm_planck_temperature(sensorindex_imager,4,imager_info(76,i),data_all(maxinfo+7+j,itx))
-                 data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero)
+                   iexponent = -(nint(imager_info(77,i)) -11)                         ! channel 15 radiance std dev for each cluster.
+                   imager_info(78,i) =  imager_info(78,i) * imager_conversion(1) * (ten ** iexponent) 
 
-                 iexponent = -(nint(imager_info(80,i)) -11)                        ! channel 16 radiance for each cluster
-                 imager_info(81,i) =  imager_info(81,i) * imager_conversion(2) * (ten ** iexponent)
+                   iexponent = -(nint(imager_info(80,i)) -11)                         ! channel 16 radiance for each cluster
+                   imager_info(81,i) =  imager_info(81,i) * imager_conversion(2) * (ten ** iexponent)
 
-                 iexponent = -(nint(imager_info(82,i))-5 )                        ! channel 16 radiance std dev for each cluster.
-                 iexponent = -(nint(imager_info(82,i)) -11)                        ! channel 16 radiance std dev for each cluster.
-                 imager_info(83,i) =  imager_info(83,i) * imager_conversion(2) * (ten ** iexponent)
+                   iexponent = -(nint(imager_info(82,i)) -11)                         ! channel 16 radiance std dev for each cluster.
+                   imager_info(83,i) =  imager_info(83,i) * imager_conversion(2) * (ten ** iexponent)
 
-                 call crtm_planck_temperature(sensorindex_imager,5,imager_info(81,i),data_all(maxinfo+14+j,itx))
-                 data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero)
-
-
-                end do imager_cluster_info
+                   call crtm_planck_temperature(sensorindex_imager,4,imager_info(76,i),data_all(maxinfo+7+j,itx))
+                   data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero)
+                   call crtm_planck_temperature(sensorindex_imager,5,imager_info(81,i),data_all(maxinfo+14+j,itx))
+                   data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero)
+                 else
+                   data_all(maxinfo+j,itx) = zero                                     ! something is wrong 
+                   data_all(maxinfo+7+j,itx) = zero                                   ! set everything to zero
+                   data_all(maxinfo+14+j,itx) = zero
+                 endif
+               end do imager_cluster_info
 
 !  Compute cluster averages for each channel
 
-                imager_mean(1) = sum(imager_cluster_size(:) * imager_info(76,:))     ! Channel 15 radiance cluster average
-                imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(76,:)**2 + imager_info(78,:)**2)) - imager_mean(1)**2
-                imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero))                ! Channel 15 radiance RMSE
-                call crtm_planck_temperature(sensorindex_imager,4,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1))
-                call crtm_planck_temperature(sensorindex_imager,4,imager_mean(1),imager_mean(1))    ! Channel 15 average BT
-                imager_std_dev(1) = imager_std_dev(1) - imager_mean(1)               ! Channel 15 BT std dev
-                data_all(maxinfo+22,itx) = imager_std_dev(1)
-
-                imager_mean(2) = sum(imager_cluster_size(:) * imager_info(81,:))     ! Channel 16 radiance cluster average
-                imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(81,:)**2 + imager_info(83,:)**2)) - imager_mean(1)**2
-                imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero))                ! Channel 16 radiance RMSE
-                call crtm_planck_temperature(sensorindex_imager,5,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2))
-                call crtm_planck_temperature(sensorindex_imager,5,imager_mean(2),imager_mean(2))    ! Channel 16 average BT
-                imager_std_dev(2) = imager_std_dev(2) - imager_mean(2)               ! Channel 16 BT std dev
-                data_all(maxinfo+23,itx) = imager_std_dev(2)
+               imager_mean(1) = sum(imager_cluster_size(:) * imager_info(76,:))       ! Channel 15 radiance cluster average
+               imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(76,:)**2 + imager_info(78,:)**2)) - imager_mean(1)**2
+               imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero))                  ! Channel 15 radiance RMSE
+               if ( imager_mean(1) > zero .and. imager_std_dev(1) > zero ) then
+                 call crtm_planck_temperature(sensorindex_imager,4,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1))
+                 call crtm_planck_temperature(sensorindex_imager,4,imager_mean(1),imager_mean(1))    ! Channel 15 average BT
+                 imager_std_dev(1) = imager_std_dev(1) - imager_mean(1)               ! Channel 15 BT std dev
+                 data_all(maxinfo+22,itx) = imager_std_dev(1)
+               else
+                 data_all(maxinfo+22,itx) = zero
+               endif
+
+               imager_mean(2) = sum(imager_cluster_size(:) * imager_info(81,:))       ! Channel 16 radiance cluster average
+               imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(81,:)**2 + imager_info(83,:)**2)) - imager_mean(1)**2
+               imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero))                  ! Channel 16 radiance RMSE
+               if ( imager_mean(2) > zero .and. imager_std_dev(2) > zero ) then
+                 call crtm_planck_temperature(sensorindex_imager,5,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2))
+                 call crtm_planck_temperature(sensorindex_imager,5,imager_mean(2),imager_mean(2))    ! Channel 16 average BT
+                 imager_std_dev(2) = imager_std_dev(2) - imager_mean(2)               ! Channel 16 BT std dev
+                 data_all(maxinfo+23,itx) = imager_std_dev(2)
+               else
+                 data_all(maxinfo+23,itx) = zero
+               endif
 
              else    !  Imager cluster information is missing.  Set everything to zero
-               data_all(maxinfo+1 : maxinfo+25,itx) = zero
+               data_all(maxinfo+1 : maxinfo+cads_info,itx) = zero
              endif
           endif    ! cris_cads
 
diff --git a/src/gsi/read_diag.f90 b/src/gsi/read_diag.f90
index e389a708d5..f69f2f5133 100644
--- a/src/gsi/read_diag.f90
+++ b/src/gsi/read_diag.f90
@@ -88,7 +88,7 @@ module read_diag
 ! Declare structures for radiance diagnostic file information
   type diag_header_fix_list
      character(len=20) :: isis           ! sat and sensor type
-     character(len=10) :: id             ! sat type
+     character(len=11) :: id             ! sat type
      character(len=10) :: obstype        ! observation type
      integer(i_kind) :: jiter            ! outer loop counter
      integer(i_kind) :: nchan            ! number of channels in the sensor
@@ -414,7 +414,8 @@ subroutine read_radiag_header_nc(ftin,header_fix,header_chan,iflag)
   real(r_kind),allocatable,dimension(:)               :: r_var_stor
   integer(i_kind),allocatable,dimension(:)            :: i_var_stor
   character(20)                          :: isis
-  character(10)                          :: id, obstype
+  character(11)                          :: id
+  character(10)                          :: obstype
 !  integer(i_kind),dimension(:),allocatable           :: jiter, nchan_diag, npred, idate, &
   integer(i_kind)                        :: jiter, nchan_diag, npred, idate, &
                                             ireal, ipchan, iextra, jextra,   &
@@ -515,7 +516,8 @@ subroutine read_radiag_header_bin(ftin,npred_radiag,retrieval,header_fix,header_
 
 !  Declare local variables
   character(len=2):: string
-  character(len=10):: satid,sentype
+  character(len=11):: satid
+  character(len=10):: sentype
   character(len=20):: sensat
   integer(i_kind) :: i,ich
   integer(i_kind):: jiter,nchanl,npred,ianldate,ireal,ipchan,iextra,jextra
diff --git a/src/gsi/read_gnssrspd.f90 b/src/gsi/read_gnssrspd.f90
new file mode 100644
index 0000000000..859345a9cc
--- /dev/null
+++ b/src/gsi/read_gnssrspd.f90
@@ -0,0 +1,463 @@
+subroutine read_gnssrspd(nread,ndata,nodata,infile,obstype,lunout,twind,sis,&
+                       nobs)
+
+!$$$  subprogram documentation block
+!                .      .    .                                       .
+! subprogram:  read_gnssrspd            read obs from gnssrspd bufr file
+!   prgmmr: kapodaca          org: Spire Global, Inc.                date: 2022-03-12
+!   Largely based on other read_* routines. For the  !Winds --- surface wind speed section
+!   ~L438, we used the exixting read_fl_hdob.f90 routine
+!
+! abstract:  This routine reads GNSSRSPD L2 wind speed observations
+!
+!            When running the gsi in regional mode, the code only
+!            retains those observations that fall within the regional
+!            domain
+
+! program history log:
+!   2022-03-12  k apodaca- initial coding
+!
+!   input argument list:
+!     infile    - unit from which to read BUFR data
+!     obstype   - observation type to process
+!     lunout    - unit to which to write data for further processing
+!     twind     - input group time window (hours)
+!     sis       - satellite/instrument/sensor indicator
+!
+!   output argument list:
+!     nread     - number of type "obstype" observations read
+!     nodata    - number of individual "obstype" observations read
+!     ndata     - number of type "obstype" observations retained for further processing
+!     nobs     - array of observations on each subdomain for each processor
+!
+! attributes:
+!   language: f90
+!   machine:  
+!
+!$$$
+     use kinds, only: r_single,r_kind,r_double,i_kind
+     use constants, only: zero,one_tenth,one,two,ten,deg2rad,t0c,half,&
+         three,four,rad2deg,r0_01,&
+         r60inv,r10,r100,r2000,hvap
+     use gridmod, only: diagnostic_reg,regional,nlon,nlat,&
+         tll2xy,txy2ll,rotate_wind_ll2xy,rotate_wind_xy2ll,&
+         rlats,rlons,twodvar_regional
+     use convinfo, only: nconvtype, &
+         icuse,ictype,ioctype
+     use obsmod, only: ran01dom
+     use obsmod, only: iadate,bmiss,offtime_data
+     use gsi_4dvar, only: l4dvar,l4densvar,time_4dvar,winlen
+     use qcmod, only: errormod
+     use convthin, only: make3grids,del3grids
+     use ndfdgrids,only: init_ndfdgrid,destroy_ndfdgrid,relocsfcob,adjust_error
+     use deter_sfc_mod, only: deter_sfc_type,deter_sfc2
+     use mpimod, only: npe
+                                                                                                      
+     implicit none
+
+!    Declare passed variables
+     character(len=*), intent(in   ) :: infile,obstype
+     character(len=20),intent(in   ) :: sis
+     integer(i_kind) , intent(in   ) :: lunout
+     integer(i_kind) , dimension(npe), intent(inout) :: nobs
+     integer(i_kind) , intent(inout) :: nread,ndata,nodata
+     real(r_kind)    , intent(in   ) :: twind
+   
+!    Declare local variables
+!    Logical variables
+     logical :: outside 
+     logical :: inflate_error
+     logical :: lspdob
+
+!    Character variables
+     character(40) :: timestr,locstr,wndstr,oestr  
+     character( 8) :: subset
+     character( 8) :: c_prvstg,c_sprvstg
+     character( 8) :: c_station_id
+     character(10) date
+!    Integer variables
+     integer(i_kind), parameter :: mxib  = 31
+
+     integer(i_kind) :: i,k
+     integer(i_kind) :: ihh,idd,idate,im,iy
+     integer(i_kind) :: lunin 
+     integer(i_kind) :: ireadmg,ireadsb
+     integer(i_kind) :: ilat,ilon 
+     integer(i_kind) :: nlv
+     integer(i_kind) :: nreal,nchanl
+     integer(i_kind) :: idomsfc,isflg
+     integer(i_kind) :: iout 
+     integer(i_kind) :: nc,ncsave
+     integer(i_kind) :: ntmatch,ntb
+     integer(i_kind) :: nmsg   
+     integer(i_kind) :: maxobs 
+     integer(i_kind) :: itype,iecol
+     integer(i_kind) :: qcm,lim_qm
+     integer(i_kind) :: wspd_qm
+     integer(i_kind) :: ntest,nvtest
+     integer(i_kind) :: igood
+     integer(i_kind) :: iuse
+ 
+     integer(i_kind) :: idate5(5)
+     integer(i_kind) :: minobs,minan
+
+     integer(i_kind), allocatable,dimension(:) :: isort
+
+!    Real variables
+     real(r_kind), parameter :: r0_001  =  0.001_r_kind
+     real(r_kind), parameter :: r1_2    =    1.2_r_kind
+     real(r_kind), parameter :: r3_0    =    3.0_r_kind
+     real(r_kind), parameter :: r0_7    =    0.7_r_kind
+     real(r_kind), parameter :: r6      =    6.0_r_kind
+     real(r_kind), parameter :: r50     =   50.0_r_kind
+     real(r_kind), parameter :: r1200   = 1200.0_r_kind
+     real(r_kind), parameter :: emerr   =    0.2_r_kind ! RH
+     real(r_kind), parameter :: missing = 1.0e+11_r_kind
+     real(r_kind), parameter :: r180    = 180.0_r_kind
+     real(r_kind), parameter :: r360    = 360.0_r_kind
+
+     real(r_kind) :: toff,t4dv
+     real(r_kind) :: usage
+     real(r_kind) :: woe,obserr
+     real(r_kind) :: dlat,dlon,dlat_earth,dlon_earth
+     real(r_kind) :: dlat_earth_deg,dlon_earth_deg
+     real(r_kind) :: cdist,disterr,disterrmax,rlon00,rlat00
+     real(r_kind) :: vdisterrmax
+     real(r_kind) :: spdob
+     real(r_kind) :: gob
+     real(r_kind) :: dlnpsob
+     real(r_kind) :: tdiff 
+     real(r_kind) :: tsavg,ff10,sfcr,zz
+     real(r_kind) :: errmin
+     real(r_kind) :: log100  
+
+     real(r_kind) :: obstime(6,1)
+     real(r_kind) :: obsloc(2,1)
+     real(r_kind) :: gnssrw(2,1)
+      
+     real(r_double) :: rstation_id
+     real(r_double) :: r_prvstg(1,1),r_sprvstg(1,1)
+
+     real(r_kind), allocatable,dimension(:,:) :: cdata_all,cdata_out
+
+!    Equivalence to handle character names
+     equivalence(r_prvstg(1,1),c_prvstg)
+     equivalence(r_sprvstg(1,1),c_sprvstg)
+     equivalence(rstation_id,c_station_id)
+
+!    Data 
+     !data timestr  / 'YEAR MNTH DAYS HOUR MINU SECO' /
+     data timestr / 'DHR RPT' /
+     !data locstr   / 'CLAT CLON' /
+     data locstr   / 'XOB YOB' /
+     !data wndstr   / 'WSPD' / !GNSSRSPD Wind speed
+     data wndstr   / 'SOB' / !GNSSRSPD Wind speed
+     data oestr   / 'WSU' / !GNSSRSPD Wind speed uncertainty/error 
+     data lunin    / 13 /
+ 
+!------------------------------------------------------------------------------------------------
+
+     write(6,*)'READ_GNSSRSPD: begin to read gnssrspd satellite data ...'
+
+!    Initialize parameters
+
+!    Set common variables
+     lspdob = obstype == 'gnssrspd'
+
+     nreal  = 0
+     iecol  = 0
+     log100=log(100._r_kind)
+ 
+ 
+     lim_qm = 4
+     iecol=0
+     if (lspdob) then
+        nreal  = 23
+        iecol  =  4  
+        errmin = one
+     else 
+        write(6,*) ' illegal obs type in read_gnssrspd '
+        call stop2(94)
+     end if
+
+     inflate_error = .true.
+
+! Check if the obs type specified in the convinfo is in the fl hdob bufr file 
+! If found, get the index (nc) from the convinfo for the specified type
+ntmatch =  0
+ncsave  =  0
+do nc = 1, nconvtype
+    if (trim(ioctype(nc)) == trim(obstype)) then 
+        if (trim(ioctype(nc)) == 'gnssrspd' .and. ictype(nc) == 600 ) then
+            ntmatch = ntmatch + 1
+            ncsave  = nc
+            itype   = ictype(nc)
+        end if
+    end if
+end do
+if (ntmatch == 0) then  ! Return if not specified in convinfo 
+    write(6,*) ' READ_GNSSRSPD: No matching obstype found in obsinfo ', obstype
+    return
+else 
+    nc = ncsave
+    write(6,*) ' READ_GNSSRSPD: Processing GNSSRSPD data : ', ntmatch, nc, ioctype(nc), ictype(nc), itype 
+end if
+
+
+!------------------------------------------------------------------------------------------------
+
+!    Go through the bufr file to find out how mant subsets to process
+     nmsg   = 0
+     maxobs = 0
+     call closbf(lunin) 
+     open(lunin,file=trim(infile),form='unformatted')
+     call openbf(lunin,'IN',lunin)
+     call datelen(10)
+     
+     loop_msg1: do while(ireadmg(lunin,subset,idate) >= 0)
+        if(nmsg == 0) call time_4dvar(idate,toff)   ! time offset (hour)
+
+        nmsg = nmsg+1
+        loop_readsb1: do while(ireadsb(lunin) == 0)
+           maxobs = maxobs+1     
+        end do loop_readsb1
+     end do loop_msg1
+     call closbf(lunin)
+     write(6,*) 'READ_GNSSRSPD: total number of data found in the bufr file ',maxobs,obstype      
+     write(6,*) 'READ_GNSSRSPD: time offset is ',toff,' hours'
+
+!---------------------------------------------------------------------------------------------------
+
+!    Allocate array to hold data
+     allocate(cdata_all(nreal,maxobs))
+     allocate(isort(maxobs))
+
+!    Initialize
+     cdata_all = zero 
+     isort     = 0
+     nread     = 0
+     nchanl    = 0
+     ntest     = 0
+     nvtest    = 0
+     ilon      = 2 
+     ilat      = 3 
+
+!    Open bufr file again for reading
+     call closbf(lunin)
+     open(lunin,file=trim(infile),form='unformatted')
+     call openbf(lunin,'IN',lunin)
+     call datelen(10)
+     ntb   = 0     
+     igood = 0
+!    Loop through BUFR file
+     loop_msg2: do while(ireadmg(lunin,subset,idate) >= 0)
+        loop_readsb2: do while(ireadsb(lunin) == 0)
+
+           ntb = ntb+1
+
+           c_station_id = subset
+
+!          QC mark 9: will be monitored but not assimilated
+!          QC mark 4: reject - will not be monitored nor assimilated 
+!          QC mark 3: suspect
+!          QC mark 2: neutral or not checked 
+!          QC mark 1: good
+!          QC mark 0: keep - will be always assimilated
+           qcm     = 0 
+           wspd_qm = 0 
+
+
+!          Read observation time 
+           call ufbint(lunin,obstime,2,1,nlv,timestr) 
+
+! If date in gnssrspd file does not agree with analysis date, 
+! print message and stop program execution.
+              write(date,'( i10)') idate
+              read (date,'(i4,3i2)') iy,im,idd,ihh
+           if(offtime_data) then
+
+!             in time correction for observations to account for analysis
+              idate5(1)=iy
+              idate5(2)=im
+              idate5(3)=idd
+              idate5(4)=ihh
+              idate5(5)=0
+              call w3fs21(idate5,minobs)    !  obs ref time in minutes relative to historic date
+              idate5(1)=iadate(1)
+              idate5(2)=iadate(2)
+              idate5(3)=iadate(3)
+              idate5(4)=iadate(4)
+              idate5(5)=0
+              call w3fs21(idate5,minan)    !  analysis ref time in minutes relative to historic date
+!             Add obs reference time, then subtract analysis time to get obs time relative to analysis
+
+              tdiff=float(minobs-minan)*r60inv
+
+           else
+              tdiff=zero
+           end if
+
+           t4dv = toff+obstime(1,1)
+
+           if (l4dvar.or.l4densvar) then
+              if (t4dv < zero .OR. t4dv > winlen) cycle loop_readsb2
+           else
+              if (abs(tdiff)>twind) cycle loop_readsb2
+           endif
+           nread = nread+1
+
+           usage = zero                ! will be considered for assimilation
+                                       ! subject to further QC in setupt subroutine
+           iuse  = icuse(nc)           ! assimilation flag 
+           if (iuse <=0) usage = r100  ! will be monitored but not assimilated
+
+!          Read observation location (lat/lon degree) 
+           call ufbint(lunin,obsloc,2,1,nlv,locstr)
+
+           if (obsloc(1,1) == missing .or. abs(obsloc(1,1)) < -180.0_r_kind .or. &     
+               obsloc(1,1) == missing .or. abs(obsloc(1,1)) >  180.0_r_kind .or. &
+               obsloc(2,1) == missing .or. abs(obsloc(2,1)) <  -90.0_r_kind .or. &
+               obsloc(2,1) == missing .or. abs(obsloc(2,1)) >   90.0_r_kind ) then 
+               write(6,*) 'READ_GNSSRSPD: bad lon/lat values: ', obsloc(1,1),obsloc(2,1)              
+               cycle loop_readsb2     
+           endif
+! GNSSRSPD BUFR longitudes are in the +/- 180 range, need to convert to 0 to 360
+! deg
+           if (obsloc(1,1) < 0.0_r_kind) obsloc(1,1) = obsloc(1,1) + 360.0_r_kind
+           !if (obsloc(1,1) >= 0.0_r_kind .or. obsloc(1,1) <= 180.0_r_kind) obsloc(1,1) = obsloc(1,1) + 180.0_r_kind
+ 
+           dlon_earth_deg = obsloc(1,1)
+           dlat_earth_deg = obsloc(2,1)
+           dlon_earth = obsloc(1,1)*deg2rad ! degree to radian
+           dlat_earth = obsloc(2,1)*deg2rad ! degree to radian
+
+!          Convert obs lat/lon to rotated coordinate and check 
+!          if the obs is outside of the regional domain
+           if (regional) then
+              call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)
+              if (diagnostic_reg) then
+                 call txy2ll(dlon,dlat,rlon00,rlat00)
+                 ntest      = ntest+1
+                 cdist      = sin(dlat_earth)*sin(rlat00)+cos(dlat_earth)*cos(rlat00)* &
+                             (sin(dlon_earth)*sin(rlon00)+cos(dlon_earth)*cos(rlon00))
+                 cdist      = max(-one,min(cdist,one))
+                 disterr    = acos(cdist)*rad2deg
+                 disterrmax = max(disterrmax,disterr)
+              end if
+              if(outside) cycle loop_readsb2
+           else
+              dlon = dlon_earth
+              dlat = dlat_earth
+              call grdcrd1(dlat,rlats,nlat,1)
+              call grdcrd1(dlon,rlons,nlon,1)
+           endif
+
+
+!          Read surface wind speed [m/s] from GNSSRSPD 
+           if (lspdob) then
+              !usage=r100 
+              !usage = zero  ! will be considered for assimilation
+
+              ! Get Wind Speed observations from bufr file
+              call ufbint(lunin,gnssrw,1,1,nlv,wndstr)
+              spdob = gnssrw(1,1) ! surface wind speed 
+
+              ! Don't permit observations with ws <= 1 m/s
+              if (spdob <= one) cycle loop_readsb2
+              if (spdob >= missing) cycle loop_readsb2
+
+              ! Get observation error from bufr file
+              call ufbint(lunin,gnssrw,1,1,nlv,oestr)
+              obserr = max(gnssrw(1,1),1.5_r_kind) ! surface wind speed error
+           endif
+
+
+           if ( .not. twodvar_regional) then
+              call deter_sfc_type(dlat_earth,dlon_earth,t4dv,isflg,tsavg)
+           endif
+
+           ! Get information from surface file necessary for conventional data
+           call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,tsavg,ff10,sfcr,zz)                                                                      
+           ! Process data passed quality control 
+           igood = igood + 1
+           ndata = ndata + 1
+           nodata = nodata + 1
+           iout = ndata
+
+!         Read extrapolated surface pressure [pa] and convert to [cb]
+           dlnpsob = log100         ! default (1000mb)           
+
+!-------------------------------------------------------------------------------------------------          
+
+           ! Winds --- surface wind speed 
+           if (lspdob) then
+              woe = obserr
+              !if (inflate_error) woe = woe*r3_0
+              !if (inflate_error) woe = woe*r1_2
+              !if (qcm > lim_qm ) woe = woe*1.0e6_r_kind
+              cdata_all( 1,iout)=woe                    ! wind error
+              cdata_all( 2,iout)=dlon                   ! grid relative longitude             
+              cdata_all( 3,iout)=dlat                   ! grid relative latitude                  
+              cdata_all( 4,iout)=dlnpsob                ! ln(surface pressure in cb)
+              cdata_all( 5,iout)=spdob*sqrt(two)*half   ! u obs
+              cdata_all( 6,iout)=spdob*sqrt(two)*half   ! v obs
+              cdata_all( 7,iout)=rstation_id            ! station id
+              cdata_all( 8,iout)=t4dv                   ! time
+              cdata_all( 9,iout)=nc                     ! type
+              cdata_all(10,iout)=r10                    !  elevation of observation ! 10-m wind       
+              cdata_all(11,iout)=qcm                    !  quality mark 
+              cdata_all(12,iout)=obserr                 !  original obs error 
+              cdata_all(13,iout)=usage                  ! usage parameter 
+              cdata_all(14,iout)=idomsfc                !  dominate surface type        
+              cdata_all(15,iout)=tsavg                  ! skin temperature 
+              cdata_all(16,iout)=ff10                   ! 10 meter wind factor     
+              cdata_all(17,iout)=sfcr                   ! surface roughness 
+              cdata_all(18,iout)=dlon_earth_deg         ! earth relative longitude (degree)                
+              cdata_all(19,iout)=dlat_earth_deg         ! earth relative latitude (degree)                  
+              cdata_all(20,iout)=gob                    !  station elevation (m)    
+              cdata_all(21,iout)=zz                     !  terrain height at ob location        
+              cdata_all(22,iout)=r_prvstg(1,1)          !  provider name 
+              cdata_all(23,iout)=r_sprvstg(1,1)         !  subprovider name 
+           endif 
+
+        end do loop_readsb2
+     end do loop_msg2
+
+!    Close unit to bufr file
+     call closbf(lunin)
+ 
+!    Write header record and data to output file for further processing
+     write(6,*) "READ_GNSSRSPD: nreal=",nreal," ndata=",ndata
+     allocate(cdata_out(nreal,ndata))
+     do i=1,ndata
+        do k=1,nreal
+           cdata_out(k,i)=cdata_all(k,i)
+        end do
+     end do
+     deallocate(cdata_all)
+
+     call count_obs(ndata,nreal,ilat,ilon,cdata_out,nobs)
+     write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
+     write(lunout) cdata_out
+     deallocate(cdata_out)
+900  continue
+     if(diagnostic_reg .and. ntest>0)  write(6,*)'READ_GNSSRSPD:  ',&
+        'ntest,  disterrmax=', ntest,disterrmax
+     if(diagnostic_reg .and. nvtest>0) write(6,*)'READ_GNSSRSPD:  ',&
+        'nvtest,vdisterrmax=',ntest,vdisterrmax
+
+     if (ndata == 0) then
+        call closbf(lunin)
+        write(6,*)'READ_GNSSRSPD: no data to process'
+     endif
+     write(6,*)'READ_GNSSRSPD: nreal=',nreal
+     write(6,*)'READ_GNSSRSPD: ntb,nread,ndata,nodata=',ntb,nread,ndata,nodata
+
+     call closbf(lunin)
+     close(lunin)
+
+!    End of routine
+     return
+
+end subroutine read_gnssrspd
+
diff --git a/src/gsi/read_iasi.f90 b/src/gsi/read_iasi.f90
index edd7a9b50e..75dc50bb76 100644
--- a/src/gsi/read_iasi.f90
+++ b/src/gsi/read_iasi.f90
@@ -142,9 +142,9 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,&
   integer(i_kind)  ,intent(in   ) :: mype_sub
   integer(i_kind)  ,intent(in   ) :: npe_sub
   integer(i_kind)  ,intent(in   ) :: mpi_comm_sub  
-  character(len=*), intent(in   ) :: infile
-  character(len=10),intent(in   ) :: jsatid
-  character(len=*), intent(in   ) :: obstype
+  character(len=*) ,intent(in   ) :: infile
+  character(len=*) ,intent(in   ) :: jsatid
+  character(len=*) ,intent(in   ) :: obstype
   character(len=20),intent(in   ) :: sis
   real(r_kind)     ,intent(in   ) :: twind
   real(r_kind)     ,intent(inout) :: val_iasi
@@ -233,7 +233,6 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,&
   real(r_kind),dimension(33,7) :: imager_info
   real(r_kind),dimension(7)    :: imager_cluster_size
   real(r_kind),dimension(2)    :: imager_mean, imager_std_dev
-  real(r_kind)                 :: imager_cluster_tot
 
 ! Set standard parameters
   character(8),parameter:: fov_flag="crosstrk"
@@ -813,10 +812,9 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,&
            if ( iasi_cads ) then
              call ufbseq(lnbufr,imager_info,33,7,iret,'IASIL1CS')
              if (iret == 7 .and. imager_info(3,1) <= 100.0_r_kind .and. &
-                  imager_info(3,1) >= zero .and. imager_coeff ) then   ! if imager cluster info exists
+                  sum(imager_info(3,:)) > zero .and. imager_coeff ) then   ! if imager cluster info exists
                imager_mean = zero
                imager_std_dev = zero
-               imager_cluster_tot = zero
                imager_cluster_flag = .TRUE.
                imager_cluster_size = imager_info(3,1:7)
                imager_cluster_size(:) = imager_cluster_size(:) / sum(imager_cluster_size(:))
@@ -834,49 +832,62 @@ subroutine read_iasi(mype,val_iasi,ithin,isfcalc,rmesh,jsatid,gstime,&
                imager_cluster_info: do j=1,7
                  i = imager_cluster_index(j)
 
-                 data_all(maxinfo+j,itx) =  imager_cluster_size(i)                ! Imager cluster fraction
-                 imager_cluster_tot = imager_cluster_tot + imager_info(3,i)
+!   If the cluster size, or radiance values of channel 4 or 5 are zero, do not compute statistics for that cluster
+                 if ( imager_cluster_size(i) > zero .and. imager_info(26,i) > zero .and. imager_info(31,i) > zero ) then
+                   data_all(maxinfo+j,itx) =  imager_cluster_size(i)                   ! Imager cluster fraction
 
-                 iexponent = -(nint(imager_info(25,i))-5 )                        ! channel 4 radiance for each cluster.
-                 imager_info(26,i) =  imager_info(26,i) * (ten ** iexponent)
+                   iexponent = -(nint(imager_info(25,i))-5 )                           ! channel 4 radiance for each cluster.
+                   imager_info(26,i) =  imager_info(26,i) * (ten ** iexponent)
 
-                 iexponent = -(nint(imager_info(27,i))-5 )                        ! channel 4 radiance std dev for each cluster.
-                 imager_info(28,i) =  imager_info(28,i) * (ten ** iexponent)
+                   iexponent = -(nint(imager_info(27,i))-5 )                           ! channel 4 radiance std dev for each cluster.
+                   imager_info(28,i) =  imager_info(28,i) * (ten ** iexponent)
 
-                 call crtm_planck_temperature(sensorindex_imager,2,imager_info(26,i),data_all(maxinfo+7+j,itx))
-                 data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero)
+                   iexponent = -(nint(imager_info(30,i))-5 )                           ! channel 5 radiance for each cluster
+                   imager_info(31,i) =  imager_info(31,i) * (ten ** iexponent)
 
-                 iexponent = -(nint(imager_info(30,i))-5 )                        ! channel 5 radiance for each cluster
-                 imager_info(31,i) =  imager_info(31,i) * (ten ** iexponent)
+                   iexponent = -(nint(imager_info(32,i))-5 )                           ! channel 5 radiance std dev for each cluser.
+                   imager_info(33,i) =  imager_info(33,i) * (ten ** iexponent)
 
-                 iexponent = -(nint(imager_info(32,i))-5 )                        ! channel 5 radiance std dev for each cluser.
-                 imager_info(33,i) =  imager_info(33,i) * (ten ** iexponent)
-
-                 call crtm_planck_temperature(sensorindex_imager,3,imager_info(31,i),data_all(maxinfo+14+j,itx))
-                 data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero)
+                   call crtm_planck_temperature(sensorindex_imager,2,imager_info(26,i),data_all(maxinfo+7+j,itx))
+                   data_all(maxinfo+7+j,itx) = max(data_all(maxinfo+7+j,itx),zero)
+                   call crtm_planck_temperature(sensorindex_imager,3,imager_info(31,i),data_all(maxinfo+14+j,itx))
+                   data_all(maxinfo+14+j,itx) = max(data_all(maxinfo+14+j,itx),zero)
+                 else                                                                  ! something is wrong
+                   data_all(maxinfo+j,itx) = zero                                      ! set everything to zero
+                   data_all(maxinfo+7+j,itx) = zero
+                   data_all(maxinfo+14+j,itx) = zero
+                 endif
 
                end do imager_cluster_info
 
 ! Compute cluster averages for each channel
 
-               imager_mean(1) = sum(imager_cluster_size(:) * imager_info(26,:))      ! Channel 4 radiance cluster average
+               imager_mean(1) = sum(imager_cluster_size(:) * imager_info(26,:))        ! Channel 4 radiance cluster average
                imager_std_dev(1) = sum(imager_cluster_size(:) * (imager_info(26,:)**2 + imager_info(28,:)**2)) - imager_mean(1)**2
-               imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero))                 ! Channel 4 radiance RMSE
-               call crtm_planck_temperature(sensorindex_imager,2,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1))
-               call crtm_planck_temperature(sensorindex_imager,2,imager_mean(1),imager_mean(1))    ! Channel 4 average BT
-               imager_std_dev(1) = imager_std_dev(1) - imager_mean(1)                ! Channel 4 BT std dev
-               data_all(maxinfo+22,itx) = imager_std_dev(1)
-
-               imager_mean(2) = sum(imager_cluster_size(:) * imager_info(31,:))      ! Channel 5 radiance cluster average
+               imager_std_dev(1) = sqrt(max(imager_std_dev(1),zero))                   ! Channel 4 radiance RMSE
+               if ( imager_mean(1) > zero .and. imager_std_dev(1) > zero ) then
+                 call crtm_planck_temperature(sensorindex_imager,2,(imager_std_dev(1) + imager_mean(1)),imager_std_dev(1))
+                 call crtm_planck_temperature(sensorindex_imager,2,imager_mean(1),imager_mean(1))    ! Channel 4 average BT
+                 imager_std_dev(1) = imager_std_dev(1) - imager_mean(1)                ! Channel 4 BT std dev
+                 data_all(maxinfo+22,itx) = imager_std_dev(1)
+               else
+                 data_all(maxinfo+22,itx) = zero
+               endif
+
+               imager_mean(2) = sum(imager_cluster_size(:) * imager_info(31,:))        ! Channel 5 radiance cluster average
                imager_std_dev(2) = sum(imager_cluster_size(:) * (imager_info(31,:)**2 + imager_info(33,:)**2)) - imager_mean(1)**2
-               imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero))                 ! Channel 5 radiance RMSE
-               call crtm_planck_temperature(sensorindex_imager,3,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2))
-               call crtm_planck_temperature(sensorindex_imager,3,imager_mean(2),imager_mean(2))     ! Channel 5 average BT
-               imager_std_dev(2) = imager_std_dev(2) - imager_mean(2)                ! Channel 5 BT std dev
-               data_all(maxinfo+23,itx) = imager_std_dev(2)
+               imager_std_dev(2) = sqrt(max(imager_std_dev(1),zero))                   ! Channel 5 radiance RMSE
+               if ( imager_mean(2) > zero .and. imager_std_dev(2) > zero ) then 
+                 call crtm_planck_temperature(sensorindex_imager,3,(imager_std_dev(2) + imager_mean(2)),imager_std_dev(2))
+                 call crtm_planck_temperature(sensorindex_imager,3,imager_mean(2),imager_mean(2))     ! Channel 5 average BT
+                 imager_std_dev(2) = imager_std_dev(2) - imager_mean(2)                ! Channel 5 BT std dev
+                 data_all(maxinfo+23,itx) = imager_std_dev(2)
+               else
+                 data_all(maxinfo+23,itx) = zero
+               endif
 
              else  ! Imager cluster information is missing.  Set everything to zero
-               data_all(maxinfo+1 : maxinfo+25,itx) = zero
+               data_all(maxinfo+1 : maxinfo+cads_info,itx) = zero
              endif
            endif ! iasi_cads = .true.
 !
diff --git a/src/gsi/read_obs.F90 b/src/gsi/read_obs.F90
index aa0f11b4e3..b81bdb5752 100644
--- a/src/gsi/read_obs.F90
+++ b/src/gsi/read_obs.F90
@@ -141,6 +141,9 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread)
 !   2019-08-21  H. Shao  - add METOPC-C, COSMIC-2 and PAZ to the GPS check list                           
 !   2020-05-21  H. Shao  - add commercial GNSSRO (Spire, PlanetIQ, GeoOptics) and other existing missions to the check list                           
 !   2021-02-20  X.Li      - add viirs-m and get_hsst
+!   2022-03-12  K. Apodaca - Enable CYGNSS ocean wind speed reading capability 
+!   2023-03-12  K. Apodaca - Enable GNSS-R L2 ocean wind speed reading (CYGNSS, Spire)
+!   2023-03-12  K. Apodaca - Enable GNSS-R DDM reading capability 
 !
 !   input argument list:
 !    lexist    - file status
@@ -487,6 +490,16 @@ subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread)
             endif
             nread = nread + 1
          end do loop_hdob
+       else if(trim(filename) == 'gnssrwndbufr')then
+         lexist = .false.
+         gnssrwndloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
+           if (trim(dtype)=='gnssrspd') then
+               lexist = .true.
+               exit gnssrwndloop
+           endif
+           nread = nread + 1
+         end do gnssrwndloop
+
        else if(trim(dtype) == 'pm2_5')then
           if (oneobtest_chem .and. oneob_type_chem=='pm2_5') then
              lexist=.true.
@@ -778,7 +791,7 @@ subroutine read_obs(ndata,mype)
     logical,dimension(ndat):: belong,parallel_read,ears_possible,db_possible
     logical :: modis,use_sfc_any
     logical :: acft_profl_file
-    character(10):: obstype,platid
+    character(10):: obstype
     character(22):: string
     character(120):: infile
     character(20):: sis
@@ -901,7 +914,7 @@ subroutine read_obs(ndata,mype)
        if (obstype == 't'  .or. obstype == 'uv' .or. &
            obstype == 'q'  .or. obstype == 'ps' .or. &
            obstype == 'pw' .or. obstype == 'spd'.or. &
-           obstype == 'sst'.or. &
+           obstype == 'sst'.or. obstype == 'gnssrspd'.or. &
            obstype == 'tcp'.or. obstype == "lag".or. &
            obstype == 'dw' .or. obstype == 'rw' .or. &
            obstype == 'mta_cld' .or. obstype == 'gos_ctp' .or. &
@@ -1402,7 +1415,6 @@ subroutine read_obs(ndata,mype)
        task_belongs: &
        if (i > 0 .and. belong(i)) then
 
-          platid=dplat(i)                    !     platid   - satellites to read
           obstype=dtype(i)                   !     obstype  - observation types to process
           infile=trim(dfile(i))              !     infile   - units from which to read data
           sis=dsis(i)                        !     sensor/instrument/satellite indicator
@@ -1449,8 +1461,13 @@ subroutine read_obs(ndata,mype)
                    call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
                         prsl_full,nobs_sub1(1,i),read_rec(i))
                    string='READ_PREPBUFR'
-
                 endif
+
+             else if (obstype == 'gnssrspd' .and. index(infile,'gnssrwndbufr') /=0 ) then
+                call read_gnssrspd(nread,npuse,nouse,infile,obstype,lunout,twind,sis, &
+                      nobs_sub1(1,i))
+                string='READ_GNSSRSPD'
+             
              else if(obstype == 'howv') then
                  if ( index(infile,'satmar') /=0) then
                    
@@ -1535,12 +1552,12 @@ subroutine read_obs(ndata,mype)
 
 !            Process conventional SST (nsstbufr, at this moment) data
              elseif ( obstype == 'sst' ) then
-                string="--"//trim(ditype(i))//":sst:"//trim(platid)
-                if ( platid == 'nsst') then
+                string="--"//trim(ditype(i))//":sst:"//trim(dplat(i))
+                if ( dplat(i) == 'nsst') then
                    call read_nsstbufr(nread,npuse,nouse,gstime,infile,obstype, &
                         lunout,twind,sis,nobs_sub1(1,i))
                    string='READ_NSSTBUFR'
-                elseif ( platid == 'mods') then
+                elseif ( dplat(i) == 'mods') then
                    call read_modsbufr(nread,npuse,nouse,gstime,infile,obstype, &
                         lunout,twind,sis,nobs_sub1(1,i))
                    string='READ_MODSBUFR'
@@ -1698,12 +1715,12 @@ subroutine read_obs(ndata,mype)
 
 !            Process TOVS 1b data
              rad_obstype_select: &
-             if (platid /= 'aqua' .and. (obstype == 'amsua' .or. &
+             if (dplat(i) /= 'aqua' .and. (obstype == 'amsua' .or. &
                   obstype == 'amsub' .or. obstype == 'msu'   .or.  &
                   obstype == 'mhs'   .or. obstype == 'hirs4' .or.  &
                   obstype == 'hirs3' .or. obstype == 'hirs2' .or.  &
                   obstype == 'ssu' )) then
-                call read_bufrtovs(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
+                call read_bufrtovs(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), nobs_sub1(1,i), &
                      read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use,radmod)
@@ -1711,7 +1728,7 @@ subroutine read_obs(ndata,mype)
 
 !            Process atms data
              else if (obstype == 'atms') then
-                call read_atms(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
+                call read_atms(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i),&
                      read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use,radmod)
@@ -1719,7 +1736,7 @@ subroutine read_obs(ndata,mype)
 
 !            Process saphir data
              else if (obstype == 'saphir') then
-                call read_saphir(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
+                call read_saphir(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),  &
                      nobs_sub1(1,i),dval_use)
@@ -1727,9 +1744,9 @@ subroutine read_obs(ndata,mype)
 
 
 !            Process airs data        
-             else if(platid == 'aqua' .and. (obstype == 'airs' .or.   &
+             else if(dplat(i) == 'aqua' .and. (obstype == 'airs' .or.   &
                   obstype == 'amsua'  .or.  obstype == 'hsb' ))then
-                call read_airs(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
+                call read_airs(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                      nobs_sub1(1,i),read_rec(i),dval_use)
@@ -1737,7 +1754,7 @@ subroutine read_obs(ndata,mype)
 
 !            Process iasi data
              else if(obstype == 'iasi')then
-                call read_iasi(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
+                call read_iasi(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i), &
                      read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use)
@@ -1745,7 +1762,7 @@ subroutine read_obs(ndata,mype)
 
 !            Process cris data
              else if(obstype == 'cris' .or. obstype =='cris-fsr' )then
-                call read_cris(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
+                call read_cris(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i), &
                      read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use)
@@ -1756,7 +1773,7 @@ subroutine read_obs(ndata,mype)
              else if (obstype == 'sndr' .or.                            &
                       obstype == 'sndrd1' .or. obstype == 'sndrd2' .or. &
                       obstype == 'sndrd3' .or. obstype == 'sndrd4') then
-                call read_goesndr(mype,val_dat,ithin,rmesh,platid,&
+                call read_goesndr(mype,val_dat,ithin,rmesh,dplat(i),&
                      infile,lunout,obstype,nread,npuse,nouse,twind,gstime,sis,&
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                      nobs_sub1(1,i),read_rec(i),dval_use)
@@ -1764,7 +1781,7 @@ subroutine read_obs(ndata,mype)
 
 !            Process ssmi data
              else if (obstype == 'ssmi' ) then
-                call read_ssmi(mype,val_dat,ithin,rmesh,platid,gstime,&
+                call read_ssmi(mype,val_dat,ithin,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                      nobs_sub1(1,i),read_rec(i),dval_use)
@@ -1773,7 +1790,7 @@ subroutine read_obs(ndata,mype)
 !            Process amsre data
              else if ( obstype == 'amsre_low' .or. obstype == 'amsre_mid' .or. &
                        obstype == 'amsre_hig' ) then
-                call read_amsre(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
+                call read_amsre(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                      nobs_sub1(1,i),read_rec(i),dval_use)
@@ -1783,7 +1800,7 @@ subroutine read_obs(ndata,mype)
              else if (obstype == 'ssmis'     .or. &
                       obstype == 'ssmis_las' .or. obstype == 'ssmis_uas' .or. &
                       obstype == 'ssmis_img' .or. obstype == 'ssmis_env' ) then
-                call read_ssmis(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
+                call read_ssmis(mype,val_dat,ithin,isfcalc,rmesh,dplat(i),gstime,&
                       infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                       mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                       nobs_sub1(1,i),read_rec(i),dval_use)
@@ -1791,7 +1808,7 @@ subroutine read_obs(ndata,mype)
 
 !            Process AMSR2 data
              else if(obstype == 'amsr2')then
-                call read_amsr2(mype,val_dat,ithin,rmesh,platid,gstime,&
+                call read_amsr2(mype,val_dat,ithin,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),  &
                      nobs_sub1(1,i))
@@ -1799,7 +1816,7 @@ subroutine read_obs(ndata,mype)
 
 !            Process GOES IMAGER RADIANCE  data
              else if(obstype == 'goes_img') then
-                call read_goesimg(mype,val_dat,ithin,rmesh,platid,gstime,&
+                call read_goesimg(mype,val_dat,ithin,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                      nobs_sub1(1,i),read_rec(i),dval_use)
@@ -1807,7 +1824,7 @@ subroutine read_obs(ndata,mype)
 
 !            Process GMI data
              else if (obstype == 'gmi') then
-                call read_gmi(mype,val_dat,ithin,rmesh,platid,gstime,&
+                call read_gmi(mype,val_dat,ithin,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),  &
                      nobs_sub1(1,i),dval_use)
@@ -1815,14 +1832,14 @@ subroutine read_obs(ndata,mype)
 
 !            Process Meteosat SEVIRI RADIANCE  data
              else if(obstype == 'seviri') then
-                 call read_seviri(mype,val_dat,ithin,rmesh,platid,gstime,&
+                 call read_seviri(mype,val_dat,ithin,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                      nobs_sub1(1,i),read_rec(i),dval_use)
                 string='READ_SEVIRI'
 !            Process GOES-R ABI RADIANCE  data
              else if(obstype == 'abi') then
-                call read_abi(mype,val_dat,ithin,rmesh,platid,gstime,&
+                call read_abi(mype,val_dat,ithin,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                      nobs_sub1(1,i),read_rec(i),dval_use)
@@ -1830,7 +1847,7 @@ subroutine read_obs(ndata,mype)
 
         !    Process Himawari-8 AHI RADIANCE  data
              else if(obstype == 'ahi') then
-                call read_ahi(mype,val_dat,ithin,rmesh,platid,gstime,&
+                call read_ahi(mype,val_dat,ithin,rmesh,dplat(i),gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),  &
                      nobs_sub1(1,i),read_rec(i),dval_use)
@@ -1839,7 +1856,7 @@ subroutine read_obs(ndata,mype)
 
 !            Process NAVY AVHRR RADIANCE  data
                else if(obstype == 'avhrr_navy') then
-                  call read_avhrr_navy(mype,val_dat,ithin,rmesh,platid,gstime,&
+                  call read_avhrr_navy(mype,val_dat,ithin,rmesh,dplat(i),gstime,&
                        infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                        mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                        nobs_sub1(1,i),read_rec(i),dval_use)
@@ -1847,14 +1864,14 @@ subroutine read_obs(ndata,mype)
 
   !            Process NESDIS AVHRR RADIANCE  data
                else if(obstype == 'avhrr') then
-                  call read_avhrr(mype,val_dat,ithin,rmesh,platid,gstime,&
+                  call read_avhrr(mype,val_dat,ithin,rmesh,dplat(i),gstime,&
                        infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                        mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                        nobs_sub1(1,i),read_rec(i),dval_use)
                   string='READ_AVHRR'
   !            Process SST VIIRS RADIANCE  data
                else if(obstype == 'viirs-m') then
-                  call read_sst_viirs(mype,val_dat,ithin,rmesh,platid,gstime,&
+                  call read_sst_viirs(mype,val_dat,ithin,rmesh,dplat(i),gstime,&
                        infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                        mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                        nobs_sub1(1,i),read_rec(i),dval_use)
@@ -1867,7 +1884,7 @@ subroutine read_obs(ndata,mype)
              if (is_extOzone(infile,obstype,dplat(i))) then
 
                 call extOzone_read(infile,obstype,dplat(i),dsis(i), &
-                   iread,ipuse,iouse, platid,gstime,lunout,twind,ithin,rmesh, &
+                   iread,ipuse,iouse, dplat(i),gstime,lunout,twind,ithin,rmesh, &
                    nobs_sub1(:,i))
                 string='extOzone_read'
 
@@ -1877,7 +1894,7 @@ subroutine read_obs(ndata,mype)
 
              else
                 call read_ozone(nread,npuse,nouse,&
-                   platid,infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, &
+                   dplat(i),infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, &
                    nobs_sub1(1,i))
                 string='READ_OZONE'
              endif ozone_obstype_select
@@ -1903,7 +1920,7 @@ subroutine read_obs(ndata,mype)
 !         Process aerosol data
           else if (ditype(i) == 'aero' )then
              call read_aerosol(nread,npuse,nouse,&
-                  platid,infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, &
+                  dplat(i),infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, &
                   mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                   nobs_sub1(1,i))
              string='READ_AEROSOL'
@@ -1929,7 +1946,7 @@ subroutine read_obs(ndata,mype)
                call warn('read_obs','                infile =',trim(infile))
                call warn('read_obs','             ditype(i) =',trim(ditype(i)))
                call warn('read_obs','               obstype =',trim(obstype))
-               call warn('read_obs','                platid =',trim(platid))
+               call warn('read_obs','                 dplat =',trim(dplat(i)))
                call warn('read_obs','                string =',trim(string))
              endif
 
diff --git a/src/gsi/read_viirs.f90 b/src/gsi/read_viirs.f90
index ab80642f29..3ea50352ec 100644
--- a/src/gsi/read_viirs.f90
+++ b/src/gsi/read_viirs.f90
@@ -475,7 +475,7 @@ subroutine read_sst_viirs(mype,val_viirs,ithin,rmesh,jsatid,&
           nele,itxmax,nread,ndata_mesh,data_mesh,score_crit,nrec)
 
      if ( nread > 0 ) then
-        write(*,'(a,a10,I3,F6.1,3I10)') 'read_viirs,jsatid,imesh,amesh,itxmax,nread,ndata_mesh :',jsatid,imesh,amesh(imesh),itxmax,nread,ndata_mesh
+        write(*,'(a,a11,I3,F6.1,3I10)') 'read_viirs,jsatid,imesh,amesh,itxmax,nread,ndata_mesh :',jsatid,imesh,amesh(imesh),itxmax,nread,ndata_mesh
      endif
 !
 !    get data_all by combining data from all thinning box sizes
diff --git a/src/gsi/set_crtm_cloudmod.f90 b/src/gsi/set_crtm_cloudmod.f90
index 9e54bd73ce..b52acad2bd 100644
--- a/src/gsi/set_crtm_cloudmod.f90
+++ b/src/gsi/set_crtm_cloudmod.f90
@@ -187,7 +187,7 @@ subroutine setCloud (cloud_name, icmask, cloud_cont, cloud_efr,jcloud, dp, tp, p
               cloud(n)%Effective_Radius(:) = cloud_efr(:,n)
            else
              !cloud(n)%Effective_Radius(:) = EftSize_(cloud_name(jcloud(n))) 
-              if ( imp_physics==11 .and. lprecip ) then
+              if (lprecip) then
                  cloud(n)%Effective_Radius(:) = cloud_efr(:,n)
               else
                  do k = 1, km
diff --git a/src/gsi/setupgnssrspd.f90 b/src/gsi/setupgnssrspd.f90
new file mode 100644
index 0000000000..e90b94a1bd
--- /dev/null
+++ b/src/gsi/setupgnssrspd.f90
@@ -0,0 +1,940 @@
+module gnssrspd_setup
+  implicit none
+  private
+  public:: setup
+        interface setup; module procedure setupgnssrspd; end interface
+
+contains
+subroutine setupgnssrspd(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsave)
+!$$$  subprogram documentation block
+!                .      .    .                                       .
+! subprogram:    setupgnssrspd    compute rhs of oi for wind speed obs
+!   prgmmr: Karina Apodaca org: Spire Global, Inc.    date: 2023-03-15
+!   Based on: setup10mspd.f90 by:  
+!   prgmmr: parrish          org: np22                date: 1990-10-06
+!
+! abstract:  For wind speed observations, this routine
+!              a) reads obs assigned to given mpi task (geographic region),
+!              b) simulates obs from guess,
+!              c) apply some quality control to obs,
+!              d) load weight and innovation arrays used in minimization
+!              e) collects statistics for runtime diagnostic output
+!              f) writes additional diagnostic information to output file
+!
+! program history log:
+!   2023-03-15  K. Apodaca - add GNSS-R L2 ocean wind speed  
+!
+!   input argument list:
+!     lunin    - unit from which to read observations
+!     mype     - mpi task id
+!     nele     - number of data elements per observation
+!     nobs     - number of observations
+!
+!   output argument list:
+!     bwork    - array containing information about obs-ges statistics
+!     awork    - array containing information for data counts and gross checks
+!
+! attributes:
+!   language: f90
+!   machine:  ibm RS/6000 SP
+!
+!$$$
+  use mpeu_util, only: die,perr,getindex
+  use kinds, only: r_kind,r_single,r_double,i_kind
+  use m_obsdiagNode, only: obs_diag
+  use m_obsdiagNode, only: obs_diags
+  use m_obsdiagNode, only: obsdiagLList_nextNode
+  use m_obsdiagNode, only: obsdiagNode_set
+  use m_obsdiagNode, only: obsdiagNode_get
+  use m_obsdiagNode, only: obsdiagNode_assert
+
+  use obsmod, only: rmiss_single,&
+                    lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset,&
+                    lobsdiag_forenkf
+  use obsmod, only: netcdf_diag, binary_diag, dirname, ianldate
+  use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, &
+       nc_diag_write, nc_diag_data2d
+  use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_get_dim, nc_diag_read_close
+  use m_obsNode, only: obsNode
+  use m_gnssrspdNode, only: gnssrspdNode
+  use m_gnssrspdNode, only: gnssrspdNode_appendto
+  use m_obsLList, only: obsLList
+  use obsmod, only: luse_obsdiag
+  use gsi_4dvar, only: nobs_bins,hr_obsbin
+  use guess_grids, only: nfldsig,hrdifsig,ges_lnprsl, &
+           comp_fact10,sfcmod_gfs,sfcmod_mm5
+  use guess_grids, only: geop_hgtl
+  use gridmod, only: nsig,get_ij,twodvar_regional
+  use qcmod, only: npres_print,ptop,pbot
+  use constants, only: one,grav,rd,zero,four,tiny_r_kind, &
+       half,two,cg_term,huge_single,r1000,wgtlim,ten
+  use jfunc, only: jiter,last,miter,jiterstart
+  use state_vectors, only: svars3d, levels
+  use qcmod, only: dfact,dfact1
+  use convinfo, only: nconvtype,cermin,cermax,cgross,cvar_b,cvar_pg,ictype
+  use convinfo, only: icsubtype
+  use gsi_bundlemod, only : gsi_bundlegetpointer
+  use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle
+  use m_dtime, only: dtime_setup, dtime_check
+  use sparsearr, only: sparr2, new, size, writearray, fullarray
+
+  ! The following variables are the coefficients that describe the
+  ! linear regression fits that are used to define the dynamic
+  ! observation error (DOE) specifications for all reconnissance
+  ! observations collected within hurricanes/tropical cyclones; these
+  ! apply only to the regional forecast models (e.g., HWRF); Henry
+  ! R. Winterbottom (henry.winterbottom@noaa.gov).
+
+  
+  implicit none
+
+! Declare passed variables
+  type(obsLList ),target,dimension(:),intent(in):: obsLL
+  type(obs_diags),target,dimension(:),intent(in):: odiagLL
+
+  logical                                          ,intent(in   ) :: conv_diagsave
+  integer(i_kind)                                  ,intent(in   ) :: lunin,mype,nele,nobs
+  real(r_kind),dimension(100+7*nsig)               ,intent(inout) :: awork
+  real(r_kind),dimension(npres_print,nconvtype,5,3),intent(inout) :: bwork
+  integer(i_kind)                                  ,intent(in   ) :: is ! ndat index
+
+! Declare local variables
+  character(len=*),parameter:: myname='setupgnssrspd'
+
+! Declare external calls for code analysis
+  external:: tintrp2a1,tintrp2a11
+  external:: tintrp31
+  external:: grdcrd1
+  external:: stop2
+
+! Declare local variables
+  
+  real(r_double) rstation_id
+  
+  real(r_kind) uob,vob,gnssrspdges,gnssrspdob,gnssrspdob0,goverrd,ratio_errors
+  real(r_kind) presw,factw,dpres,ugesin,vgesin,sfcr,skint
+  real(r_kind) scale
+  real(r_kind) val2,ressw,ress,error,ddiff,dx10,rhgh,prsfc,r0_001
+  real(r_kind) sfcchk,prsln2,rwgt,tfact                        
+  real(r_kind) thirty,rsig,ratio,residual,obserrlm,obserror
+  real(r_kind) val,valqc,psges,drpx,dlat,dlon,dtime,rlow
+  real(r_kind) cg_gnssrspd,wgross,wnotgross,wgt,arg,exp_arg,term,rat_err2
+  real(r_kind) errinv_input,errinv_adjst,errinv_final
+  real(r_kind) err_input,err_adjst,err_final
+  real(r_kind),dimension(nele,nobs):: data
+  real(r_kind),dimension(nobs):: dup
+  real(r_kind),dimension(nsig)::prsltmp,tges
+  real(r_single),allocatable,dimension(:,:)::rdiagbuf
+
+  integer(i_kind) mm1,ibin,ioff,ioff0
+  integer(i_kind) ii,jj,i,nchar,nreal,k,j,l,nty,nn,ikxx
+  integer(i_kind) ier,ilon,ilat,ipres,iuob,ivob,id,itime,ikx
+  integer(i_kind) ihgt,iqc,ier2,iuse,ilate,ilone,istnelv,izz,iprvd,isprvd
+  integer(i_kind) idomsfc,iskint,iff10,isfcr,isli
+
+  type(sparr2) :: dhx_dx
+  integer(i_kind) :: iz, u_ind, v_ind, nnz, nind
+  real(r_kind) :: delz
+  
+  logical,dimension(nobs):: luse,muse
+  integer(i_kind),dimension(nobs):: ioid ! initial (pre-distribution) obs ID
+  logical proceed
+
+  character(8) station_id
+  character(8),allocatable,dimension(:):: cdiagbuf
+  character(8),allocatable,dimension(:):: cprvstg,csprvstg
+  character(8) c_prvstg,c_sprvstg
+  real(r_double) r_prvstg,r_sprvstg
+
+  logical:: in_curbin, in_anybin, save_jacobian
+  type(gnssrspdNode),pointer:: my_head
+  type(obs_diag),pointer:: my_diag
+  type(obs_diags),pointer:: my_diagLL
+
+  logical z_height
+  real(r_kind) zsges,dstn
+  real(r_kind),dimension(nsig):: zges
+  real(r_kind) dz,zob,z1,z2,p1,p2,dz21,dlnp21,pobl
+  integer(i_kind) k1,k2
+
+  equivalence(rstation_id,station_id)
+  equivalence(r_prvstg,c_prvstg)
+  equivalence(r_sprvstg,c_sprvstg)
+
+  real(r_kind),allocatable,dimension(:,:,:  ) :: ges_ps
+  real(r_kind),allocatable,dimension(:,:,:  ) :: ges_z
+  real(r_kind),allocatable,dimension(:,:,:,:) :: ges_u
+  real(r_kind),allocatable,dimension(:,:,:,:) :: ges_v
+  real(r_kind),allocatable,dimension(:,:,:,:) :: ges_tv
+
+  type(obsLList),pointer,dimension(:):: gnssrspdhead
+  gnssrspdhead => obsLL(:)
+
+!******************************************************************************
+! Read and reformat observations in work arrays.
+  read(lunin)data,luse,ioid
+
+!        index information for data array (see reading routine)
+  ier=1       ! index of obs error
+  ilon=2      ! index of grid relative obs location (x)
+  ilat=3      ! index of grid relative obs location (y)
+  ipres=4     ! index of pressure
+  iuob=5      ! index of u observation
+  ivob=6      ! index of v observation
+  id=7        ! index of station id
+  itime=8     ! index of observation time in data array
+  ikxx=9      ! index of ob type
+  ihgt=10     ! index of observation elevation
+  iqc=11      ! index of quality mark
+  ier2=12     ! index of original-original obs error ratio
+  iuse=13     ! index of use parameter
+  idomsfc=14  ! index of dominant surface type
+  iskint=15   ! index of surface skin temperature
+  iff10=16    ! index of 10 meter wind factor
+  isfcr=17    ! index of surface roughness
+  ilone=18    ! index of longitude (degrees)
+  ilate=19    ! index of latitude (degrees)
+  istnelv=20  ! index of station elevation (m)
+  izz=21      ! index of surface height
+  iprvd=22    ! index of observation provider
+  isprvd=23   ! index of  observation subprovider
+
+  mm1=mype+1
+  scale=one
+  rsig=nsig
+  thirty = 30.0_r_kind
+  r0_001=0.001_r_kind
+  goverrd=grav/rd
+  
+  save_jacobian = conv_diagsave .and. jiter==jiterstart .and. lobsdiag_forenkf
+
+! Check to see if required guess fields are available
+  call check_vars_(proceed)
+  if(.not.proceed) return  ! not all vars available, simply return
+
+! If require guess vars available, extract from bundle ...
+  call init_vars_
+
+! If requested, save select data for output to diagnostic file
+  if(conv_diagsave)then
+     ii=0
+     nchar=1
+     ioff0=21
+     nreal=ioff0
+     if (lobsdiagsave) nreal=nreal+4*miter+1
+     if (twodvar_regional) then; nreal=nreal+2; allocate(cprvstg(nobs),csprvstg(nobs)); endif
+     if (save_jacobian) then
+       nnz    = 4                   ! number of non-zero elements in dH(x)/dx profile
+       nind   = 2
+       call new(dhx_dx, nnz, nind)
+       nreal = nreal + size(dhx_dx)
+     endif
+     allocate(cdiagbuf(nobs),rdiagbuf(nreal,nobs))
+     if(netcdf_diag) call init_netcdf_diag_
+  end if
+
+  !!awork(:) = zero
+
+  do i=1,nobs
+     muse(i)=nint(data(iuse,i)) <= jiter
+  end do
+
+  dup=one
+  do k=1,nobs
+     do l=k+1,nobs
+        if(data(ilat,k) == data(ilat,l) .and.  &
+           data(ilon,k) == data(ilon,l) .and.  &
+           data(ipres,k)== data(ipres,l) .and. &
+           data(ier,k) < r1000 .and. data(ier,l) < r1000 .and. &
+           muse(l) .and. muse(k))then
+
+           tfact=min(one,abs(data(itime,k)-data(itime,l))/dfact1)
+           dup(k)=dup(k)+one-tfact*tfact*(one-dfact)
+           dup(l)=dup(l)+one-tfact*tfact*(one-dfact)
+        end if
+     end do
+  end do
+
+  call dtime_setup()
+  do i=1,nobs
+     dtime=data(itime,i)
+     call dtime_check(dtime, in_curbin, in_anybin)
+     if(.not.in_anybin) cycle
+
+     if(in_curbin) then
+        dlat=data(ilat,i)
+        dlon=data(ilon,i)
+ 
+        dpres=data(ipres,i)
+        error=data(ier2,i)
+        ikx=nint(data(ikxx,i))
+     endif
+
+!    Link observation to appropriate observation bin
+     if (nobs_bins>1) then
+        ibin = NINT( dtime/hr_obsbin ) + 1
+     else
+        ibin = 1
+     endif
+     IF (ibin<1.OR.ibin>nobs_bins) write(6,*)mype,'Error nobs_bins,ibin= ',nobs_bins,ibin
+
+     if (luse_obsdiag) my_diagLL => odiagLL(ibin)
+
+!    Link obs to diagnostics structure
+     if (luse_obsdiag) then
+        my_diag => obsdiagLList_nextNode(my_diagLL      ,&
+                create = .not.lobsdiag_allocated        ,&
+                   idv = is             ,&
+                   iob = ioid(i)        ,&
+                   ich = 1              ,&
+                  elat = data(ilate,i)  ,&
+                  elon = data(ilone,i)  ,&
+                  luse = luse(i)        ,&
+                 miter = miter          )
+
+        if(.not.associated(my_diag)) call die(myname, &
+                'obsdiagLList_nextNode(), create =', .not.lobsdiag_allocated)
+     endif
+
+     if(.not.in_curbin) cycle
+
+!    Load obs error and u,v obs
+     obserror = max(cermin(ikx),min(cermax(ikx),data(ier,i)))
+     uob = data(iuob,i)
+     vob = data(ivob,i)
+ 
+ 
+     gnssrspdob=sqrt(uob*uob+vob*vob)
+     call tintrp2a1(ges_tv,tges,dlat,dlon,dtime,hrdifsig,&
+          nsig,mype,nfldsig)
+     call tintrp2a11(ges_ps,psges,dlat,dlon,dtime,hrdifsig,&
+          mype,nfldsig)
+     call tintrp2a1(ges_lnprsl,prsltmp,dlat,dlon,dtime,hrdifsig,&
+          nsig,mype,nfldsig)
+
+     factw = data(iff10,i)
+     if(sfcmod_gfs .or. sfcmod_mm5)then
+        sfcr = data(isfcr,i)
+        skint = data(iskint,i)
+        isli=data(idomsfc,i)
+        call comp_fact10(dlat,dlon,dtime,skint,sfcr,isli,mype,factw)
+     end if
+
+     nty=ictype(ikx)
+
+!    Initialize logical     
+     z_height = .false.
+
+!    Process observations reported with height differently than those
+!    reported with pressure.  Type 260=nacelle 261=tower wind spd are
+!    encoded in NCEP prepbufr files with geometric height above
+!    sea level. 
+
+!    Process cygnss observations at mslp
+     if ( nty == 600 ) then
+          z_height = .true.
+          data(ihgt,i) = ten
+     endif
+     
+     if (z_height) then
+        
+        drpx = zero
+        dpres = data(ihgt,i)
+        dstn = data(istnelv,i)
+        call tintrp2a11(ges_z,zsges,dlat,dlon,dtime,hrdifsig,&
+             mype,nfldsig)
+
+!       Get guess surface elevation and geopotential height profile
+!       at observation location.
+        call tintrp2a1(geop_hgtl,zges,dlat,dlon,dtime,hrdifsig,&
+             nsig,mype,nfldsig)
+
+!       Convert observation height (in dpres) from meters to grid relative
+!       units.  Save the observation height in zob for later use.
+        zob = dpres
+        call grdcrd1(dpres,zges,nsig,1)
+        factw=one
+        rlow=zero
+        rhgh=zero
+
+!       Compute observation pressure (only used for diagnostics)
+!       Set indices of model levels below (k1) and above (k2) observation.
+        if (dpres<one) then
+           z1=zero;    p1=log(psges)
+           z2=zges(1); p2=prsltmp(1)
+        elseif (dpres>nsig) then
+           z1=zges(nsig-1); p1=prsltmp(nsig-1)
+           z2=zges(nsig);   p2=prsltmp(nsig)
+           drpx = 1.e6_r_kind
+        else
+           k=dpres
+           k1=min(max(1,k),nsig)
+           k2=max(1,min(k+1,nsig))
+           z1=zges(k1); p1=prsltmp(k1)
+           z2=zges(k2); p2=prsltmp(k2)
+        endif
+
+        dz21     = z2-z1
+        dlnp21   = p2-p1
+        dz       = zob-z1
+        pobl     = p1 + (dlnp21/dz21)*dz
+        presw    = ten*exp(pobl)
+
+        prsfc=psges
+        prsln2=log(exp(prsltmp(1))/prsfc)
+        sfcchk=log(psges)
+        if(dpres <= prsln2)then
+           factw=one
+        else
+           dx10=-goverrd*ten/tges(1)
+           if (dpres < dx10)then
+              factw=(dpres-dx10+factw*(prsln2-dpres))/(prsln2-dx10)
+           end if
+        end if
+
+!    Put obs pressure in correct units to get grid coord. number
+        dpres=log(exp(dpres)*prsfc)
+        call grdcrd1(dpres,prsltmp(1),nsig,-1)
+
+!    Get approx k value of sfc by using surface pressure of 1st ob
+        call grdcrd1(sfcchk,prsltmp(1),nsig,-1)
+
+
+!    Check to see if observations is below what is seen to be the surface
+        rlow=max(sfcchk-dpres,zero)
+
+        rhgh=max(dpres-r0_001-rsig,zero)
+
+     endif  ! end of process observations with reported pressure
+
+     if(luse(i))then
+        awork(1) = awork(1) + one
+        if(rlow/=zero) awork(2) = awork(2) + one
+        if(rhgh/=zero) awork(3) = awork(3) + one
+     end if
+
+     ratio_errors=error/(data(ier,i)+drpx+1.0e6_r_kind*rhgh+four*rlow)
+     
+
+! Interpolate guess u and v to observation location and time.
+     call tintrp31(ges_u,ugesin,dlat,dlon,dpres,dtime, &
+        hrdifsig,mype,nfldsig)
+     call tintrp31(ges_v,vgesin,dlat,dlon,dpres,dtime, &
+        hrdifsig,mype,nfldsig)
+ 
+
+! Apply 10-meter wind reduction factor to guess winds.  Compute
+! guess wind speed.
+     ugesin=factw*ugesin
+     vgesin=factw*vgesin
+     gnssrspdges=sqrt(ugesin*ugesin+vgesin*vgesin)
+
+     iz = max(1, min( int(dpres), nsig))
+     delz = max(zero, min(dpres - float(iz), one))
+
+     if (save_jacobian) then
+
+        u_ind = getindex(svars3d, 'u')
+        if (u_ind < 0) then
+           print *, 'Error: no variable u in state vector. Exiting.'
+           call stop2(1300)
+        endif
+        v_ind = getindex(svars3d, 'v')
+        if (v_ind < 0) then
+           print *, 'Error: no variable v in state vector. Exiting.'
+           call stop2(1300)
+        endif
+
+        dhx_dx%st_ind(1)  = iz               + sum(levels(1:u_ind-1))
+        dhx_dx%end_ind(1) = min(iz + 1,nsig) + sum(levels(1:u_ind-1))
+
+        dhx_dx%val(1) = (one - delz) * two * ugesin
+        dhx_dx%val(2) = delz * two * ugesin
+
+        dhx_dx%st_ind(2)  = iz               + sum(levels(1:v_ind-1))
+        dhx_dx%end_ind(2) = min(iz + 1,nsig) + sum(levels(1:v_ind-1))
+
+        dhx_dx%val(3) = (one - delz) * two * vgesin
+        dhx_dx%val(4) = delz * two * ugesin
+     endif
+
+
+     ddiff = gnssrspdob-gnssrspdges
+
+!    write(6,*) "cygnss gnssrspdob=",gnssrspdob," gnssrspdges=",gnssrspdges," ddiff=",ddiff
+     
+     error=one/error
+
+!    Check to see if observations is above the top of the model (regional mode)
+     if (dpres>rsig) ratio_errors=zero
+
+
+!    Gross error checks
+     obserror = one/max(ratio_errors*error,tiny_r_kind)
+     obserrlm = max(cermin(ikx),min(cermax(ikx),obserror))
+     residual = abs(ddiff)
+     ratio    = residual/obserrlm
+     if (ratio>cgross(ikx) .or. ratio_errors < tiny_r_kind) then
+        if (luse(i)) awork(4) = awork(4)+one
+        error = zero
+        ratio_errors = zero
+        muse(i)=.false.
+     else
+        ratio_errors=ratio_errors/sqrt(dup(i))
+     end if
+
+     if (ratio_errors*error <=tiny_r_kind) muse(i)=.false.
+     if (nobskeep>0.and.luse_obsdiag) call obsdiagNode_get(my_diag, jiter=nobskeep, muse=muse(i))
+
+!    Compute penalty terms (linear & nonlinear qc).
+     val      = error*ddiff
+     val2     = val*val
+     exp_arg  = -half*val2
+     rat_err2 = ratio_errors**2
+     if (cvar_pg(ikx) > tiny_r_kind .and. error > tiny_r_kind) then
+        arg  = exp(exp_arg)
+        wnotgross= one-cvar_pg(ikx)
+        cg_gnssrspd=cvar_b(ikx)
+        wgross = cg_term*cvar_pg(ikx)/(cg_gnssrspd*wnotgross)
+        term = log((arg+wgross)/(one+wgross))
+        wgt  = one-wgross/(arg+wgross)
+        wgt  = wgt/wgtlim
+     else
+        term = exp_arg
+        wgt  = wgtlim
+        rwgt = wgt/wgtlim
+     endif
+     valqc = -two*rat_err2*term
+
+
+!       Accumulate statistics for obs belonging to this task
+     if (luse(i) .and. muse(i)) then
+        if(rwgt < one) awork(61) = awork(61)+one
+        awork(5)=awork(5) + val2*rat_err2
+        awork(6)=awork(6) + one
+        awork(22)=awork(22) + valqc
+     end if
+
+! Loop over pressure level groupings and obs to accumulate statistics
+! as a function of observation type.
+     do k = 1,npres_print
+        if(luse(i) .and.presw >ptop(k) .and. presw<=pbot(k))then
+           ress  = scale*ddiff
+           ressw = ress*ress
+           nn=1
+           if (.not. muse(i)) then
+              nn=2
+              if(ratio_errors*error >=tiny_r_kind)nn=3
+           end if
+           bwork(k,ikx,1,nn) = bwork(k,ikx,1,nn)+one            ! count
+           bwork(k,ikx,2,nn) = bwork(k,ikx,2,nn)+ddiff          ! bias
+           bwork(k,ikx,3,nn) = bwork(k,ikx,3,nn)+ressw          ! (o-g)**2
+           bwork(k,ikx,4,nn) = bwork(k,ikx,4,nn)+val2*rat_err2  ! penalty
+           bwork(k,ikx,5,nn) = bwork(k,ikx,5,nn)+valqc          ! nonlin qc penalty
+
+        end if
+     end do
+
+     !write(6,*) "cygnss: i,luse(i),muse(i)=",i,luse(i),muse(i)
+
+     if (luse_obsdiag) then
+        call obsdiagNode_set(my_diag, wgtjo=(error*ratio_errors)**2, &
+           jiter=jiter, muse=muse(i), nldepart=gnssrspdob-sqrt(ugesin*ugesin+vgesin*vgesin))
+     endif
+
+!    If obs is "acceptable", load array with obs info for use
+!    in inner loop minimization (int* and stp* routines)
+     if (.not. last .and. muse(i)) then
+
+        allocate(my_head)
+        call gnssrspdNode_appendto(my_head,gnssrspdhead(ibin))
+
+        my_head%idv = is
+        my_head%iob = ioid(i)
+        my_head%elat= data(ilate,i)
+        my_head%elon= data(ilone,i)
+
+!       Set (i,j) indices of guess gridpoint that bound obs location
+        call get_ij(mm1,dlat,dlon,my_head%ij,my_head%wij)
+
+        my_head%factw=factw
+        do j=1,4
+           my_head%wij(j)=factw*my_head%wij(j)
+        end do
+        my_head%raterr2= ratio_errors**2     
+        my_head%res    = gnssrspdob
+        my_head%uges   = ugesin
+        my_head%vges   = vgesin
+        my_head%err2   = error**2
+        my_head%time   = dtime
+        my_head%luse   = luse(i)
+        my_head%b      = cvar_b(ikx)
+        my_head%pg     = cvar_pg(ikx)
+
+        if (luse_obsdiag) then
+           call obsdiagNode_assert(my_diag, my_head%idv,my_head%iob,1,myname,'my_diag:my_head')
+           my_head%diags => my_diag
+        endif
+
+        my_head => null()
+     end if
+! Save select output for diagnostic file
+     if(conv_diagsave .and. luse(i))then
+        ii=ii+1
+        gnssrspdob0      = sqrt(data(iuob,i)*data(iuob,i)+data(ivob,i)*data(ivob,i))
+        rstation_id = data(id,i)
+        err_input   = data(ier2,i)
+        err_adjst   = data(ier,i)
+        if (ratio_errors*error>tiny_r_kind) then
+           err_final = one/(ratio_errors*error)
+        else
+           err_final = huge_single
+        endif
+
+        errinv_input = huge_single
+        errinv_adjst = huge_single
+        errinv_final = huge_single
+        if (err_input>tiny_r_kind) errinv_input = one/err_input
+        if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst
+        if (err_final>tiny_r_kind) errinv_final = one/err_final
+
+        if(binary_diag) call contents_binary_diag_(my_diag)
+        if(netcdf_diag) call contents_netcdf_diag_(my_diag)
+
+     end if
+
+! End of loop over observations
+  end do
+
+! Release memory of local guess arrays
+  call final_vars_
+
+! Write information to diagnostic file
+  if(conv_diagsave) then
+     if(netcdf_diag) call nc_diag_write
+     if(binary_diag .and. ii>0)then
+        write(7)'gnssrspd',nchar,nreal,ii,mype,ioff0
+        write(7)cdiagbuf(1:ii),rdiagbuf(:,1:ii)
+        deallocate(cdiagbuf,rdiagbuf)
+
+        if (twodvar_regional) then
+           write(7)cprvstg(1:ii),csprvstg(1:ii)
+           deallocate(cprvstg,csprvstg)
+        endif
+     end if
+  end if
+
+! End of routine
+
+  return
+  contains
+
+  subroutine check_vars_ (proceed)
+  logical,intent(inout) :: proceed
+  integer(i_kind) ivar, istatus
+! Check to see if required guess fields are available
+  call gsi_metguess_get ('var::ps', ivar, istatus )
+  proceed=ivar>0
+  call gsi_metguess_get ('var::z' , ivar, istatus )
+  proceed=proceed.and.ivar>0
+  call gsi_metguess_get ('var::u' , ivar, istatus )
+  proceed=proceed.and.ivar>0
+  call gsi_metguess_get ('var::v' , ivar, istatus )
+  proceed=proceed.and.ivar>0
+  call gsi_metguess_get ('var::tv', ivar, istatus )
+  proceed=proceed.and.ivar>0
+  end subroutine check_vars_ 
+
+  subroutine init_vars_
+
+  real(r_kind),dimension(:,:  ),pointer:: rank2=>NULL()
+  real(r_kind),dimension(:,:,:),pointer:: rank3=>NULL()
+  character(len=5) :: varname
+  integer(i_kind) ifld, istatus
+
+! If require guess vars available, extract from bundle ...
+  if(size(gsi_metguess_bundle)==nfldsig) then
+!    get ps ...
+     varname='ps'
+     call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus)
+     if (istatus==0) then
+         if(allocated(ges_ps))then
+            write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
+            call stop2(999)
+         endif
+         allocate(ges_ps(size(rank2,1),size(rank2,2),nfldsig))
+         ges_ps(:,:,1)=rank2
+         do ifld=2,nfldsig
+            call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus)
+            ges_ps(:,:,ifld)=rank2
+         enddo
+     else
+         write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
+         call stop2(999)
+     endif
+!    get z ...
+     varname='z'
+     call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank2,istatus)
+     if (istatus==0) then
+         if(allocated(ges_z))then
+            write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
+            call stop2(999)
+         endif
+         allocate(ges_z(size(rank2,1),size(rank2,2),nfldsig))
+         ges_z(:,:,1)=rank2
+         do ifld=2,nfldsig
+            call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank2,istatus)
+            ges_z(:,:,ifld)=rank2
+         enddo
+     else
+         write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
+         call stop2(999)
+     endif
+!    get u ...
+     varname='u'
+     call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus)
+     if (istatus==0) then
+         if(allocated(ges_u))then
+            write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
+            call stop2(999)
+         endif
+         allocate(ges_u(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig))
+         ges_u(:,:,:,1)=rank3
+         do ifld=2,nfldsig
+            call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus)
+            ges_u(:,:,:,ifld)=rank3
+         enddo
+     else
+         write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
+         call stop2(999)
+     endif
+!    get v ...
+     varname='v'
+     call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus)
+     if (istatus==0) then
+         if(allocated(ges_v))then
+            write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
+            call stop2(999)
+         endif
+         allocate(ges_v(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig))
+         ges_v(:,:,:,1)=rank3
+         do ifld=2,nfldsig
+            call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus)
+            ges_v(:,:,:,ifld)=rank3
+         enddo
+     else
+         write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
+         call stop2(999)
+     endif
+!    get tv ...
+     varname='tv'
+     call gsi_bundlegetpointer(gsi_metguess_bundle(1),trim(varname),rank3,istatus)
+     if (istatus==0) then
+         if(allocated(ges_tv))then
+            write(6,*) trim(myname), ': ', trim(varname), ' already incorrectly alloc '
+            call stop2(999)
+         endif
+         allocate(ges_tv(size(rank3,1),size(rank3,2),size(rank3,3),nfldsig))
+         ges_tv(:,:,:,1)=rank3
+         do ifld=2,nfldsig
+            call gsi_bundlegetpointer(gsi_metguess_bundle(ifld),trim(varname),rank3,istatus)
+            ges_tv(:,:,:,ifld)=rank3
+         enddo
+     else
+         write(6,*) trim(myname),': ', trim(varname), ' not found in met bundle, ier= ',istatus
+         call stop2(999)
+     endif
+  else
+     write(6,*) trim(myname), ': inconsistent vector sizes (nfldsig,size(metguess_bundle) ',&
+                 nfldsig,size(gsi_metguess_bundle)
+     call stop2(999)
+  endif
+  end subroutine init_vars_
+
+  subroutine init_netcdf_diag_
+  character(len=80) string
+  character(len=128) diag_conv_file
+  integer(i_kind) ncd_fileid,ncd_nobs
+  logical append_diag
+  logical,parameter::verbose=.false. 
+     write(string,900) jiter
+900  format('conv_gnssrspd_',i2.2,'.nc4')
+     diag_conv_file=trim(dirname) // trim(string)
+
+     inquire(file=diag_conv_file, exist=append_diag)
+
+     if (append_diag) then
+        call nc_diag_read_init(diag_conv_file,ncd_fileid)
+        ncd_nobs = nc_diag_read_get_dim(ncd_fileid,'nobs')
+        call nc_diag_read_close(diag_conv_file)
+
+        if (ncd_nobs > 0) then
+           if(verbose) print *,'file ' // trim(diag_conv_file) // ' exists.  Appending.  nobs,mype=',ncd_nobs,mype
+        else
+           if(verbose) print *,'file ' // trim(diag_conv_file) // ' exists but contains no obs.  Not appending. nobs,mype=',ncd_nobs,mype
+           append_diag = .false. ! if there are no obs in existing file, then do not try to append
+        endif
+     end if
+
+     call nc_diag_init(diag_conv_file, append=append_diag)
+
+     if (.not. append_diag) then ! don't write headers on append - the module will break?
+        call nc_diag_header("date_time",ianldate )
+        if (save_jacobian) then
+          call nc_diag_header("jac_nnz", nnz)
+          call nc_diag_header("jac_nind", nind)
+        endif
+     endif
+  end subroutine init_netcdf_diag_
+  subroutine contents_binary_diag_(odiag)
+  type(obs_diag),pointer,intent(in):: odiag
+        cdiagbuf(ii)    = station_id         ! station id
+
+        rdiagbuf(1,ii)  = ictype(ikx)        ! observation type
+        rdiagbuf(2,ii)  = icsubtype(ikx)     ! observation subtype
+    
+        rdiagbuf(3,ii)  = data(ilate,i)      ! observation latitude (degrees)
+        rdiagbuf(4,ii)  = data(ilone,i)      ! observation longitude (degrees)
+        rdiagbuf(5,ii)  = data(istnelv,i)    ! station elevation (meters)
+        rdiagbuf(6,ii)  = presw              ! observation pressure (hPa)
+        rdiagbuf(7,ii)  = data(ihgt,i)       ! observation height (meters)
+        rdiagbuf(8,ii)  = dtime-time_offset  ! obs time (hours relative to analysis time)
+
+        rdiagbuf(9,ii)  = data(iqc,i)        ! input prepbufr qc or event mark
+        rdiagbuf(10,ii) = rmiss_single       ! setup qc or event mark
+        rdiagbuf(11,ii) = data(iuse,i)       ! read_prepbufr data usage flag
+        if(muse(i)) then
+           rdiagbuf(12,ii) = one             ! analysis usage flag (1=use, -1=not used)
+        else
+           rdiagbuf(12,ii) = -one
+        endif
+
+        rdiagbuf(13,ii) = rwgt               ! nonlinear qc relative weight
+        rdiagbuf(14,ii) = errinv_input       ! prepbufr inverse obs error (m/s)**-1
+        rdiagbuf(15,ii) = errinv_adjst       ! read_prepbufr inverse obs error (m/s)**-1
+        rdiagbuf(16,ii) = errinv_final       ! final inverse observation error (m/s)**-1
+
+        rdiagbuf(17,ii) = gnssrspdob              ! wind speed observation (m/s)
+        rdiagbuf(18,ii) = ddiff              ! obs-ges used in analysis (m/s)
+        rdiagbuf(19,ii) = gnssrspdob0-gnssrspdges      ! obs-ges w/o bias correction (m/s) (future slot)
+
+        rdiagbuf(20,ii) = factw              ! 10m wind reduction factor
+
+        rdiagbuf(21,ii) = 1.e+10_r_single    ! ges ensemble spread (filled in by EnKF)
+
+        ioff=ioff0
+        if (lobsdiagsave) then
+           do jj=1,miter 
+              ioff=ioff+1 
+              if (odiag%muse(jj)) then
+                 rdiagbuf(ioff,ii) = one
+              else
+                 rdiagbuf(ioff,ii) = -one
+              endif
+           enddo
+           do jj=1,miter+1
+              ioff=ioff+1
+              rdiagbuf(ioff,ii) = odiag%nldepart(jj)
+           enddo
+           do jj=1,miter
+              ioff=ioff+1
+              rdiagbuf(ioff,ii) = odiag%tldepart(jj)
+           enddo
+           do jj=1,miter
+              ioff=ioff+1
+              rdiagbuf(ioff,ii) = odiag%obssen(jj)
+           enddo
+        endif
+
+        if (twodvar_regional) then
+           ioff = ioff + 1
+           rdiagbuf(ioff,ii) = data(idomsfc,i) ! dominate surface type
+           ioff = ioff + 1
+           rdiagbuf(ioff,ii) = data(izz,i)     ! model terrain at ob location
+           r_prvstg          = data(iprvd,i)
+           cprvstg(ii)       = c_prvstg        ! provider name
+           r_sprvstg         = data(isprvd,i)
+           csprvstg(ii)      = c_sprvstg       ! subprovider name
+        endif
+
+        if (save_jacobian) then
+           call writearray(dhx_dx, rdiagbuf(ioff+1:nreal,ii))
+           ioff = ioff + size(dhx_dx)
+        endif
+
+  end subroutine contents_binary_diag_
+  subroutine contents_netcdf_diag_(odiag)
+  type(obs_diag),pointer,intent(in):: odiag
+! Observation class
+  character(8),parameter     :: obsclass = 'gnssrspd'
+  real(r_kind),dimension(miter) :: obsdiag_iuse
+           call nc_diag_metadata("Station_ID",              station_id             )
+           call nc_diag_metadata("Observation_Class",       obsclass               )
+           call nc_diag_metadata("Observation_Type",        ictype(ikx)            )
+           call nc_diag_metadata("Observation_Subtype",     icsubtype(ikx)         )
+           call nc_diag_metadata("Latitude",                sngl(data(ilate,i))    )
+           call nc_diag_metadata("Longitude",               sngl(data(ilone,i))    )
+           call nc_diag_metadata("Station_Elevation",       sngl(data(istnelv,i))  )
+           call nc_diag_metadata("Pressure",                sngl(presw)            )
+           call nc_diag_metadata("Height",                  sngl(data(ihgt,i))     )
+           call nc_diag_metadata("Time",                    sngl(dtime-time_offset))
+           call nc_diag_metadata("Prep_QC_Mark",            sngl(data(iqc,i))      )
+           call nc_diag_metadata("Prep_Use_Flag",           sngl(data(iuse,i))     )
+!          call nc_diag_metadata("Nonlinear_QC_Var_Jb",     var_jb                 )
+           call nc_diag_metadata("Nonlinear_QC_Rel_Wgt",    sngl(rwgt)             )                 
+           if(muse(i)) then
+              call nc_diag_metadata("Analysis_Use_Flag",    sngl(one)              )
+           else
+              call nc_diag_metadata("Analysis_Use_Flag",    sngl(-one)             )              
+           endif
+
+           call nc_diag_metadata("Errinv_Input",            sngl(errinv_input)     )
+           call nc_diag_metadata("Errinv_Adjust",           sngl(errinv_adjst)     )
+           call nc_diag_metadata("Errinv_Final",            sngl(errinv_final)     )
+
+           call nc_diag_metadata("Observation",                   sngl(gnssrspdob)      )
+           call nc_diag_metadata("Obs_Minus_Forecast_adjusted",   sngl(ddiff)      )
+           call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(gnssrspdob0-gnssrspdges) )
+ 
+           if (lobsdiagsave) then
+              do jj=1,miter
+                 if (odiag%muse(jj)) then
+                       obsdiag_iuse(jj) =  one
+                 else
+                       obsdiag_iuse(jj) = -one
+                 endif
+              enddo
+   
+              call nc_diag_data2d("ObsDiagSave_iuse",     obsdiag_iuse                             )
+              call nc_diag_data2d("ObsDiagSave_nldepart", odiag%nldepart )
+              call nc_diag_data2d("ObsDiagSave_tldepart", odiag%tldepart )
+              call nc_diag_data2d("ObsDiagSave_obssen",   odiag%obssen   )             
+           endif
+   
+           if (twodvar_regional) then
+              call nc_diag_metadata("Dominant_Sfc_Type", data(idomsfc,i)              )
+              call nc_diag_metadata("Model_Terrain",     data(izz,i)                  )
+              r_prvstg            = data(iprvd,i)
+              call nc_diag_metadata("Provider_Name",     c_prvstg                     )    
+              r_sprvstg           = data(isprvd,i)
+              call nc_diag_metadata("Subprovider_Name",  c_sprvstg                    )
+           endif
+           if (save_jacobian) then
+             call nc_diag_data2d("Observation_Operator_Jacobian_stind", dhx_dx%st_ind)
+             call nc_diag_data2d("Observation_Operator_Jacobian_endind", dhx_dx%end_ind)
+             call nc_diag_data2d("Observation_Operator_Jacobian_val", real(dhx_dx%val,r_single))
+           endif
+
+
+  end subroutine contents_netcdf_diag_
+
+  subroutine final_vars_
+    if(allocated(ges_tv)) deallocate(ges_tv)
+    if(allocated(ges_v )) deallocate(ges_v )
+    if(allocated(ges_u )) deallocate(ges_u )
+    if(allocated(ges_z )) deallocate(ges_z )
+    if(allocated(ges_ps)) deallocate(ges_ps)
+  end subroutine final_vars_
+
+end subroutine setupgnssrspd
+end module gnssrspd_setup
diff --git a/src/gsi/setuprhsall.f90 b/src/gsi/setuprhsall.f90
index 8075956431..6c00de5803 100644
--- a/src/gsi/setuprhsall.f90
+++ b/src/gsi/setuprhsall.f90
@@ -101,6 +101,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass)
 !   2019-03-15  Ladwig  - add option for cloud analysis in observer
 !   2019-03-28  Ladwig  - add metar cloud obs as pseudo water vapor in var analysis
 !   2020-09-08  CAPS(G. Zhao) - add 'l_use_dbz_directDA' flag not to sort obsdiag
+!   2023-03-20  K Apodaca - add GNSS-R L2 Ocean Wind Speed
 !
 !   input argument list:
 !     ndata(*,1)- number of prefiles retained for further processing
@@ -165,7 +166,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass)
   use m_rhs, only: toss_gps_sub => rhs_toss_gps
 
   use m_rhs, only: i_ps,i_uv,i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag, &
-                   i_gust,i_vis,i_pblh,i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, &
+                   i_gust,i_vis,i_pblh,i_wspd10m,i_gnssrspd,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, &
                    i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp
   use m_rhs, only: i_dbz
   use m_rhs, only: i_fed
@@ -625,7 +626,7 @@ subroutine setuprhsall(ndata,mype,init_pass,last_pass)
 !    Compute and print statistics for "conventional" data
      call statsconv(mype,&
           i_ps,i_uv,i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag, &
-          i_gust,i_vis,i_pblh,i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, &
+          i_gust,i_vis,i_pblh,i_wspd10m,i_gnssrspd,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, &
           i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_fed,i_dbz, &
           size(awork1,2),bwork1,awork1,ndata)
 
diff --git a/src/gsi/statsconv.f90 b/src/gsi/statsconv.f90
index fc105515ff..b6bc849a37 100644
--- a/src/gsi/statsconv.f90
+++ b/src/gsi/statsconv.f90
@@ -1,6 +1,6 @@
 subroutine statsconv(mype,&
      i_ps,i_uv,i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag, &
-     i_gust,i_vis,i_pblh,i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & 
+     i_gust,i_vis,i_pblh,i_wspd10m,i_gnssrspd,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv, & 
      i_tcamt,i_lcbas,i_cldch,i_uwnd10m,i_vwnd10m,&
      i_swcp,i_lwcp,i_fed,i_dbz,i_ref,bwork,awork,ndata)
 !$$$  subprogram documentation block
@@ -44,6 +44,7 @@ subroutine statsconv(mype,&
 !   2015-07-10  pondeca - add cldch
 !   2016-05-05  pondeca - add uwnd10m, vwnd10m
 !   2017-05-12  Y. Wang and X. Wang - add dbz, POC: xuguang.wang@ou.edu
+!   2022-03-15  K. Apodaca - add GNSS-R L2 ocean wind speed 
 !
 !   input argument list:
 !     mype     - mpi task number
@@ -62,6 +63,7 @@ subroutine statsconv(mype,&
 !     i_vis    - index in awork array holding vis info
 !     i_pblh   - index in awork array holding pblh info
 !     i_wspd10m- index in awork array holding wspd10m info
+!     i_gnssrspd - index in awork array holding gnssrspd info
 !     i_td2m   - index in awork array holding td2m info
 !     i_mxtm   - index in awork array holding mxtm info
 !     i_mitm   - index in awork array holding mitm info
@@ -94,13 +96,13 @@ subroutine statsconv(mype,&
   use constants, only: zero,three,five
   use obsmod, only: iout_sst,iout_pw,iout_t,iout_rw,iout_dw,&
        iout_uv,iout_gps,iout_ps,iout_q,iout_tcp,iout_lag,&
-       iout_gust,iout_vis,iout_pblh,iout_wspd10m,iout_td2m,& 
+       iout_gust,iout_vis,iout_pblh,iout_wspd10m,iout_gnssrspd,iout_td2m,& 
        iout_mxtm,iout_mitm,iout_pmsl,iout_howv,iout_tcamt,iout_lcbas,iout_cldch,&
        iout_uwnd10m,iout_vwnd10m,&
        iout_fed,iout_dbz,iout_swcp,iout_lwcp,&
        mype_dw,mype_rw,mype_sst,mype_gps,mype_uv,mype_ps,&
        mype_t,mype_pw,mype_q,mype_tcp,ndat,dtype,mype_lag,mype_gust,&
-       mype_vis,mype_pblh,mype_wspd10m,mype_td2m,mype_mxtm,mype_mitm,&
+       mype_vis,mype_pblh,mype_wspd10m,mype_gnssrspd,mype_td2m,mype_mxtm,mype_mitm,&
        mype_pmsl,mype_howv,mype_tcamt,mype_lcbas,mype_cldch,mype_uwnd10m,mype_vwnd10m,&
        mype_fed,mype_dbz,mype_swcp,mype_lwcp
   use qcmod, only: npres_print,ptop,pbot,ptopq,pbotq
@@ -112,7 +114,7 @@ subroutine statsconv(mype,&
 ! Declare passed variables
   integer(i_kind)                                  ,intent(in   ) :: mype,i_ps,i_uv,&
        i_t,i_q,i_pw,i_rw,i_dw,i_gps,i_sst,i_tcp,i_lag,i_gust,i_vis,i_pblh,&
-       i_wspd10m,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv,i_tcamt,i_lcbas,&
+       i_wspd10m,i_gnssrspd,i_td2m,i_mxtm,i_mitm,i_pmsl,i_howv,i_tcamt,i_lcbas,&
        i_cldch,i_uwnd10m,i_vwnd10m,i_swcp,i_lwcp,i_fed,i_dbz,i_ref
   real(r_kind),dimension(7*nsig+100,i_ref)     ,intent(in   ) :: awork
   real(r_kind),dimension(npres_print,nconvtype,5,3),intent(in   ) :: bwork
@@ -122,7 +124,7 @@ subroutine statsconv(mype,&
   character(100) mesage
 
   integer(i_kind) numgrspw,numsst,nsuperp,nump,nhitopo,ntoodif
-  integer(i_kind) numgrsq,numhgh,numgust,numvis,numpblh,numwspd10m,numuwnd10m,numvwnd10m
+  integer(i_kind) numgrsq,numhgh,numgust,numvis,numpblh,numwspd10m,numgnssrspd,numuwnd10m,numvwnd10m
   integer(i_kind) numtd2m,nummxtm,nummitm,numpmsl,numhowv,numtcamt,numlcbas,numcldch
   integer(i_kind) numgrsswcp,numgrslwcp
   integer(i_kind) ntot,numlow,k,numssm,i,j
@@ -677,6 +679,49 @@ subroutine statsconv(mype,&
      end if
   end if
 
+! Summary report for conventional gnssrspd
+if (mype == mype_gnssrspd) then
+    nread = 0
+    nkeep = 0
+    do i = 1, ndat
+        if (dtype(i) == 'gnssrspd') then
+            nread = nread + ndata(i, 2)
+            nkeep = nkeep + ndata(i, 3)
+        end if
+    end do
+    if (nread > 0) then
+        if (first) then
+            open(iout_gnssrspd)
+        else
+            open(iout_gnssrspd, position = 'append')
+        end if
+
+        numgnssrspd = nint(awork(5, i_gnssrspd))
+        pw = zero
+        pw3 = zero
+        if (nkeep > 0) then
+            mesage = 'current fit of conventional gnssrspd data, ranges in m/s'
+            do j = 1, nconvtype
+                pflag(j) = trim(ioctype(j)) == 'gnssrspd'
+            end do
+            call dtast(bwork, 1, pbotall, ptopall, mesage, jiter, iout_gnssrspd, pflag)
+
+            numgross = nint(awork(6, i_gnssrspd))
+            numfailqc = nint(awork(21, i_gnssrspd))
+            if (numgnssrspd > 0) then
+                pw = awork(4, i_gnssrspd) / numgnssrspd
+                pw3 = awork(22, i_gnssrspd) / numgnssrspd
+            end if
+            write(iout_gnssrspd, 925) 'gnssrspd', numgross, numfailqc
+        end if
+        write(iout_gnssrspd, 950) 'gnssrspd', jiter, nread, nkeep, numgnssrspd
+        write(iout_gnssrspd, 951) 'gnssrspd', awork(4, i_gnssrspd), awork(22, i_gnssrspd), pw, pw3
+
+        close(iout_gnssrspd)
+    end if
+end if
+
+
 ! Summary report for conventional td2m
   if(mype==mype_td2m) then
      nread=0
diff --git a/src/gsi/statsrad.f90 b/src/gsi/statsrad.f90
index d42a53f7d6..dfedcc92ab 100644
--- a/src/gsi/statsrad.f90
+++ b/src/gsi/statsrad.f90
@@ -100,7 +100,7 @@ subroutine statsrad(aivals,stats,ndata)
            write(iout_rad,4002) 'qcpenalty','qc1','qc2','qc3','qc4','qc5','qc6','qc7'    
            write(iout_rad,4003) aivals(39,is),(nint(aivals(i,is)),i=8,14)
 4000       format(2x,a3,7x,a4,14x,a7,1x,8(a7,1x))
-4001       format(1x,a10,1x,a10,f16.8,8(i7,1x))
+4001       format(1x,a11,1x,a10,f16.8,8(i7,1x))
 4002       format(28x,a9,1x,8(a7,1x))
 4003       format(22x,f16.8,8(i7,1x))
            
@@ -161,11 +161,11 @@ subroutine statsrad(aivals,stats,ndata)
 2011 format(8x,f16.8,8(i7,1x))
 2012 format(12x,A7,5x,8(a7,1x))
 2999 format(' Illegal satellite type ')
-1102 format(1x,i4,i5,1x,a16,2i7,1x,f10.3,1x,6(f11.7,1x))
+1102 format(1x,i4,i6,1x,a20,2i7,1x,f10.3,1x,6(f11.7,1x))
 1109 format(t5,'it',t13,'satellite',t23,'instrument',t40, &
           '# read',t53,'# keep',t65,'# assim',&
           t75,'penalty',t88,'qcpnlty',t104,'cpen',t115,'qccpen')
-1115 format('o-g',1x,i2.2,1x,'rad',2x,2A10,2x,3(i11,2x),4(g12.5,1x))
+1115 format('o-g',1x,i2.2,1x,'rad',2x,2A12,2x,3(i11,2x),4(g12.5,1x))
 
 ! Close output unit
   close(iout_rad)
diff --git a/src/gsi/stpgnssrspd.f90 b/src/gsi/stpgnssrspd.f90
new file mode 100644
index 0000000000..3909c40534
--- /dev/null
+++ b/src/gsi/stpgnssrspd.f90
@@ -0,0 +1,175 @@
+module stpgnssrspdmod
+
+
+!$$$ module documentation block
+!           .      .    .                                       .
+! module:   stpgnssrspdmod    module for stpgnssrspd and its tangent linear stpgnssrspd_tl
+!  prgmmr: K. Apodaca         org: Spire Global, Inc.      date: 2022-03-12
+!   Largely based on other stp_* routines 
+!
+! abstract: module for stpgnssrspd and its tangent linear stpgnssrspd_tl
+!
+! program history log:
+!   2023-09-21  K. Apodaca  - add documentation
+! subroutine included:
+!   sub stpgnssrspd
+!
+! attributes:
+!   language: f90
+!   machine:
+!
+!$$$ end documentation block
+
+implicit none
+
+PRIVATE
+PUBLIC stpgnssrspd
+
+contains
+
+subroutine stpgnssrspd(gnssrspdhead,rval,sval,out,sges,nstep)
+!$$$  subprogram documentation block
+!                .      .    .                                       .
+! subprogram:    stpgnssrspd  calculate penalty and stepsize terms
+!                for wind speed, with nonlinear qc.
+!   2023-03-15  K. Apodaca - add GNSS-R L2 ocean wind speed 
+!
+! abstract: calculate penalty and stepsize terms for wind speed
+!
+! program history log:
+
+!
+!   input argument list:
+!     gnssrspdhead
+!     ru       - search direction for u
+!     rv       - search direction for v
+!     su       - analysis increment for u
+!     sv       - analysis increment for v
+!     sges     - step size estimates (nstep)
+!     nstep    - number of stepsizes  (==0 means use outer iteration values)
+!
+!   output argument list 
+!     out(1:nstep)   - contribution to penalty from wind speed sges(1:nstep)
+!
+! attributes:
+!   language: f90
+!   machine:  ibm RS/6000 SP
+!
+!$$$
+  use kinds, only: r_kind,i_kind,r_quad
+  use qcmod, only: nlnqc_iter,varqc_iter
+  use constants, only: zero,half,one,two,tiny_r_kind,cg_term,zero_quad,r3600
+  use gsi_4dvar, only: ltlint
+  use gsi_bundlemod, only: gsi_bundle
+  use gsi_bundlemod, only: gsi_bundlegetpointer
+  use m_obsNode, only: obsNode
+  use m_gnssrspdNode, only: gnssrspdNode
+  use m_gnssrspdNode, only: gnssrspdNode_typecast
+  use m_gnssrspdNode, only: gnssrspdNode_nextcast
+  implicit none
+
+! Declare passed variables
+  class(obsNode), pointer             ,intent(in   ) :: gnssrspdhead
+  integer(i_kind)                     ,intent(in   ) :: nstep
+  real(r_kind),dimension(max(1,nstep)),intent(in   ) :: sges
+  real(r_quad),dimension(max(1,nstep)),intent(inout) :: out
+  type(gsi_bundle)                    ,intent(in   ) :: rval,sval
+
+! Declare local variables
+  integer(i_kind) j1,j2,j3,j4,kk,ier,istatus
+  real(r_kind) w1,w2,w3,w4,time_gnssrspd
+  real(r_kind) valu,valv,ucur,vcur,gnssrspdnl,gnssrspdtl,uu,vv,gnssrspd
+  real(r_kind),dimension(max(1,nstep)):: pen
+  real(r_kind) cg_gnssrspd,pencur,wgross,wnotgross
+  real(r_kind) pg_gnssrspd,pentl
+  real(r_kind),pointer,dimension(:) :: su,sv
+  real(r_kind),pointer,dimension(:) :: ru,rv
+  type(gnssrspdNode), pointer :: gnssrspdptr
+
+  out=zero_quad
+
+!  If no gnssrspd data return
+  if(.not. associated(gnssrspdhead))return
+
+  time_gnssrspd=zero
+! Retrieve pointers
+! Simply return if any pointer not found
+  ier=0
+  call gsi_bundlegetpointer(sval,'u',su,istatus);ier=istatus+ier
+  call gsi_bundlegetpointer(sval,'v',sv,istatus);ier=istatus+ier
+  call gsi_bundlegetpointer(rval,'u',ru,istatus);ier=istatus+ier
+  call gsi_bundlegetpointer(rval,'v',rv,istatus);ier=istatus+ier
+  if(ier/=0)return
+
+  gnssrspdptr => gnssrspdNode_typecast(gnssrspdhead)
+  do while (associated(gnssrspdptr))
+
+     if(gnssrspdptr%luse)then
+        if(nstep > 0)then
+           j1 = gnssrspdptr%ij(1)
+           j2 = gnssrspdptr%ij(2)
+           j3 = gnssrspdptr%ij(3)
+           j4 = gnssrspdptr%ij(4)
+           w1 = gnssrspdptr%wij(1)
+           w2 = gnssrspdptr%wij(2)
+           w3 = gnssrspdptr%wij(3)
+           w4 = gnssrspdptr%wij(4)
+ 
+           valu=w1* ru(j1)+w2* ru(j2)+w3* ru(j3)+w4* ru(j4)
+           valv=w1* rv(j1)+w2* rv(j2)+w3* rv(j3)+w4* rv(j4)
+           ucur=w1* su(j1)+w2* su(j2)+w3* su(j3)+w4* su(j4)+gnssrspdptr%uges
+           vcur=w1* sv(j1)+w2* sv(j2)+w3* sv(j3)+w4* sv(j4)+gnssrspdptr%vges
+
+           if (ltlint) then
+              gnssrspd=sqrt(ucur*ucur+vcur*vcur)-gnssrspdptr%res
+              pencur=gnssrspd*gnssrspd*gnssrspdptr%err2
+              do kk=1,nstep
+                 gnssrspdnl=sqrt(ucur*ucur+vcur*vcur)
+                 gnssrspdtl=ucur*valu+vcur*valv
+                 if (gnssrspdnl>tiny_r_kind*100._r_kind) then
+                    gnssrspdtl=gnssrspdtl/gnssrspdnl
+                 else
+                    gnssrspdtl=zero
+                 endif
+                 pentl  =two*gnssrspdtl*gnssrspd*gnssrspdptr%err2
+                 pen(kk)=pencur+sges(kk)*pentl
+              end do
+           else
+              do kk=1,nstep
+                 uu=ucur+sges(kk)*valu
+                 vv=vcur+sges(kk)*valv
+                 gnssrspd=sqrt(uu*uu+vv*vv)-gnssrspdptr%res
+                 pen(kk)=gnssrspd*gnssrspd*gnssrspdptr%err2
+              end do
+           end if
+        else
+           pen(1)=gnssrspdptr%res*gnssrspdptr%res*gnssrspdptr%err2
+        end if
+
+!  Modify penalty term if nonlinear QC
+        if (nlnqc_iter .and. gnssrspdptr%pg > tiny_r_kind .and. &
+                             gnssrspdptr%b  > tiny_r_kind) then
+           pg_gnssrspd=gnssrspdptr%pg*varqc_iter
+           cg_gnssrspd=cg_term/gnssrspdptr%b
+           wnotgross= one-pg_gnssrspd
+           wgross = pg_gnssrspd*cg_gnssrspd/wnotgross
+           pencur = -two*log((exp(-half*pencur) + wgross)/(one+wgross))
+           do kk=1,max(1,nstep)
+              pen(kk) = -two*log((exp(-half*pen(kk)  ) + wgross)/(one+wgross))
+           enddo
+        endif
+
+        out(1) = out(1)+pen(1)*gnssrspdptr%raterr2
+        do kk=2,nstep
+           out(kk) = out(kk)+(pen(kk)-pen(1))*gnssrspdptr%raterr2
+        end do
+
+     end if
+    
+     gnssrspdptr => gnssrspdNode_nextcast(gnssrspdptr)
+
+  end do
+  return
+end subroutine stpgnssrspd
+
+end module stpgnssrspdmod
diff --git a/src/mgbf/mg_input.f90 b/src/mgbf/mg_input.f90
index 80b0772c12..ad2a47034d 100644
--- a/src/mgbf/mg_input.f90
+++ b/src/mgbf/mg_input.f90
@@ -57,7 +57,7 @@ subroutine input_2d &
 integer(i_kind),intent(in):: imax0
 integer(i_kind),intent(in):: ampl
 real(r_kind),dimension(imin:imax,jmin:jmax),intent(out):: V
-real(i_kind):: ng,mg,L,m,n
+integer(i_kind):: ng,mg,L,m,n
 !-----------------------------------------------------------------------
 
      do m=imin,jmax
@@ -134,7 +134,7 @@ subroutine input_3d &
 integer(i_kind),intent(in):: imax0
 integer(i_kind),intent(in):: ampl,incrm
 real(r_kind),dimension(lmin:lmax,imin:imax,jmin:jmax),intent(out):: V
-real(i_kind):: ng,mg,L,m,n
+integer(i_kind):: ng,mg,L,m,n
 !-----------------------------------------------------------------------
 
    do l=lmin,lmax