Fix HAFS GSI debug build and run issues (#679)

**DUE DATE for merger of this PR into `develop` is 2/19/2024 (six weeks
after PR creation).**
**DUE DATE for this PR is extended to 3/19/2024 because @XuLu-NOAA is on
leave.**
**Description**

Xu Lu (xu.lu@noaa.gov) and Biju Thomas (biju.thomas@noaa.gov) fixed bugs
regarding HAFS GSI debug build and run issues. This is in corresponding
to issue #661

Fixes #661

1. In read_radar.f90, uninitialized toff is making all the ground-based
radar observations be placed at -3h instead of 0h, which creates wrong
increments for FGAT and 4DEnVar.
2. In read_radar.f90, uninitialized zsges will crash the debug mode.
3. In read_radar.f90, t4dvo should be used instead of t4dv in the
read_radar_l2rw_novadqc subroutine.
4. In radinfo.90, maxscan should be increased to at least 252 to allow
more scans, or it will crash the debug mode.
5. In read_fl_hdob.f90, dlnpsob is replaced with 1000. since the SFMR
does not sample surface pressure, and the uninitialized dlnpsob creates
issues later in setupspd.f90 in the debug mode.
6. In mod_fv3_lola.f90, (i,j+1) should be used instead of (i+1,j) in
searching for V edges.
7. In stpcalc.f90, when tried to find the best stepsize from outpen
around L838-864, the minimum outstp(i) is stored in stp(ii), but the
istp_use is asigned with i instead of ii. Create inconsistency when
assigning stp(istp_use) to stpinout at L872. Should use istp_use=ii
instead.

**Type of change**
- [Yes] Bug fix (non-breaking change which fixes an issue)

**How Has This Been Tested?**
Regression test on Orion:
```
Test project /work/noaa/hwrf/save/xulu/mergeversions/GSI/build
CMake Warning (dev) at CTestTestfile.cmake:9 (subdirs):  Syntax Warning in cmake code at /work/noaa/hwrf/save/xulu/mergeversions/GSI/build/regression/CTestTestfile.cmake:7:10
1/7 Test #4: [=[netcdf_fv3_regional]=] ........   Passed  365.11 sec
2/7 Test #7: [=[global_enkf]=] ................   Passed  430.29 sec
3/7 Test #3: [=[rrfs_3denvar_glbens]=] ........   Passed  605.35 sec
4/7 Test #2: [=[rtma]=] .......................   Passed  969.78 sec
5/7 Test #6: [=[hafs_3denvar_hybens]=] ........***Failed  1455.47 sec
6/7 Test #1: [=[global_4denvar]=] .............   Passed  1682.40 sec
7/7 Test #5: [=[hafs_4denvar_glbens]=] ........***Failed  1758.90 sec
```

The failed hafs_3denvar and 4denvar are within expectation due to the
fix for toff. As demonstrated in the single observation tests in the
following figure, the uninitialized toff can result in increment
degradations due to wrongly assigned observation times:

![image](https://github.com/NOAA-EMC/GSI/assets/26603014/0de870e1-f8c8-4b6d-8039-57f417b76367)

src/gsi/mod_fv3_lola.f90

@@ -951,7 +951,6 @@ subroutine definecoef_regular_grids(nxen,nyen,grid_lon,grid_lont,grid_lat,grid_l
      do i=1,nxen
 !      center lat/lon of the edge 
         rlat=half*(grid_lat(i,j)+grid_lat(i+1,j))
-!        rlon=half*(grid_lon(i,j)+grid_lon(i+1,j))
         diff=(grid_lon(i,j)-grid_lon(i+1,j))**2
         if(diff < sq180)then
            rlon=half*(grid_lon(i,j)+grid_lon(i+1,j))
@@ -979,12 +978,11 @@ subroutine definecoef_regular_grids(nxen,nyen,grid_lon,grid_lont,grid_lat,grid_l
   do j=1,nyen
      do i=1,nxen+1
         rlat=half*(grid_lat(i,j)+grid_lat(i,j+1))
-!        rlon=half*(grid_lon(i,j)+grid_lon(i,j+1))
-        diff=(grid_lon(i,j)-grid_lon(i+1,j))**2
+        diff=(grid_lon(i,j)-grid_lon(i,j+1))**2
         if(diff < sq180)then
-           rlon=half*(grid_lon(i,j)+grid_lon(i+1,j))
+           rlon=half*(grid_lon(i,j)+grid_lon(i,j+1))
         else
-           rlon=half*(grid_lon(i,j)+grid_lon(i+1,j)-360._r_kind)
+           rlon=half*(grid_lon(i,j)+grid_lon(i,j+1)-360._r_kind)
         endif
         xr=cos(rlat*deg2rad)*cos(rlon*deg2rad)
         yr=cos(rlat*deg2rad)*sin(rlon*deg2rad)

src/gsi/radinfo.f90

@@ -897,7 +897,7 @@ subroutine radinfo_read
 !   Allocate arrays to receive angle dependent bias information.
 !   Open file to bias file (satang=satbias_angle).  Read data.
 
-    maxscan=250
+    maxscan=252
     if (.not.adp_anglebc) maxscan = 90 ! default value for old files
 
     if (adp_anglebc) then 

src/gsi/read_fl_hdob.f90

@@ -48,7 +48,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si
      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,tiny_r_kind,huge_r_kind,r0_01,&
-         r60inv,r10,r100,r2000,hvap,eps,omeps,rv,grav
+         r60inv,r10,r100,r2000,hvap,eps,omeps,rv,grav,r_missing
      use gridmod, only: diagnostic_reg,regional,nlon,nlat,nsig,&
          tll2xy,txy2ll,rotate_wind_ll2xy,rotate_wind_xy2ll,&
          rlats,rlons,twodvar_regional,fv3_regional
@@ -1133,7 +1133,7 @@ subroutine read_fl_hdob(nread,ndata,nodata,infile,obstype,lunout,gstime,twind,si
               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( 4,iout)=r_missing                ! ln(surface pressure in cb) !Since dlnpsob is not provided by SFMR, force it to be r_missing. Not used in setupspd.f90
               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

src/gsi/read_radar.f90

@@ -3282,7 +3282,7 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs)
   use convinfo, only: nconvtype,ncmiter,ncgroup,ncnumgrp,icuse,ioctype,pmot_conv
   use deter_sfc_mod, only: deter_sfc2,deter_zsfc_model
   use mpimod, only: npe
-  use obsmod, only: reduce_diag 
+  use obsmod, only: reduce_diag,time_offset 
 
   implicit none
 
@@ -3323,7 +3323,7 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs)
   integer(i_kind) iret,kx0
   integer(i_kind) nreal,nchanl,ilat,ilon,ikx
   integer(i_kind) idomsfc
-  real(r_kind) usage,ff10,sfcr,skint,t4dv,t4dvo,toff
+  real(r_kind) usage,ff10,sfcr,skint,t4dvo
   real(r_kind) eradkm,dlat_earth,dlon_earth
   real(r_kind) dlat,dlon,staheight,tiltangle,clon,slon,clat,slat
   real(r_kind) timeo,clonh,slonh,clath,slath,cdist,dist
@@ -3467,7 +3467,7 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs)
      staheight=this_stahgt    !station elevation
      tiltangle=corrected_tilt*deg2rad
 
-     t4dvo=toff+thistime
+     t4dvo=thistime+time_offset
      timemax=max(timemax,t4dvo)
      timemin=min(timemin,t4dvo)
 
@@ -3586,15 +3586,16 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs)
         end if
         if(usage >= 100._r_kind)rusage(ndata)=.true.
 
-        call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr)
+        call deter_sfc2(dlat_earth,dlon_earth,t4dvo,idomsfc,skint,ff10,sfcr)
+        call deter_zsfc_model(dlat,dlon,zsges)
 
         cdata(1) = error             ! wind obs error (m/s)
         cdata(2) = dlon              ! grid relative longitude
         cdata(3) = dlat              ! grid relative latitude
         cdata(4) = height            ! obs absolute height (m)
         cdata(5) = rwnd              ! wind obs (m/s)
         cdata(6) = azm*deg2rad       ! azimuth angle (radians)
-        cdata(7) = t4dv              ! obs time (hour)
+        cdata(7) = t4dvo             ! obs time (hour)
         cdata(8) = ikx               ! type
         cdata(9) = tiltangle         ! tilt angle (radians)
         cdata(10)= staheight         ! station elevation (m)
@@ -3699,6 +3700,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
   use obsmod, only: doradaroneob,oneobradid,time_offset,reduce_diag 
   use mpeu_util, only: gettablesize,gettable 
   use convinfo, only: nconvtype,icuse,ioctype
+  use deter_sfc_mod, only: deter_sfc2,deter_zsfc_model
   use mpimod, only: npe
   use read_l2bufr_mod, only: radar_sites,radar_rmesh,radar_zmesh,elev_angle_max,del_time,range_max,radar_pmot  
   use constants, only: eccentricity,somigliana,grav_ratio,grav,semi_major_axis,flattening,grav_equator 
@@ -3744,7 +3746,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
   integer(i_kind) iret,kx0
   integer(i_kind) nreal,nchanl,ilat,ilon,ikx
   integer(i_kind) idomsfc
-  real(r_kind) usage,ff10,sfcr,skint,t4dv,t4dvo,toff
+  real(r_kind) usage,ff10,sfcr,skint,t4dvo
   real(r_kind) eradkm,dlat_earth,dlon_earth
   real(r_kind) dlat,dlon,staheight,tiltangle,clon,slon,clat,slat
   real(r_kind) timeo,clonh,slonh,clath,slath,cdist,dist
@@ -4071,7 +4073,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
         slat=sin(dlat_earth)
         staheight=this_stahgt    !station elevation
         tiltangle=corrected_tilt*deg2rad
-        t4dvo=toff+thistime
+        t4dvo=time_offset+thistime
         timemax=max(timemax,t4dvo)
         timemin=min(timemin,t4dvo)
 !    Exclude data if it does not fall within time window
@@ -4166,7 +4168,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
           if (l4dvar) then
             timedif = zero
           else
-            timedif=abs(t4dvo-toff) 
+            timedif=abs(t4dvo-time_offset)
           endif
           crit1 = timedif/r6+half
           call map3grids_m(1,save_all,zflag,zl_thin,nlevz, &
@@ -4233,8 +4235,8 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
           end if
 
 !    Get information from surface file necessary for conventional data here
-!         call deter_zsfc_model(dlat,dlon,zsges)
-!         call deter_sfc2(dlat_earth,dlon_earth,t4dv,idomsfc,skint,ff10,sfcr)
+          call deter_zsfc_model(dlat,dlon,zsges)
+          call deter_sfc2(dlat_earth,dlon_earth,t4dvo,idomsfc,skint,ff10,sfcr)
 
           nsuper2_kept=nsuper2_kept+1
           cdata(1) = error             ! wind obs error (m/s)
@@ -4243,7 +4245,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
           cdata(4) = height            ! obs absolute height (m)
           cdata(5) = rwnd              ! wind obs (m/s)
           cdata(6) = azm*deg2rad       ! azimuth angle (radians)
-          cdata(7) = t4dvo+time_offset ! obs time (hour)
+          cdata(7) = t4dvo             ! obs time (hour)
           cdata(8) = ikx               ! type
           cdata(9) = tiltangle         ! tilt angle (radians)
           cdata(10)= staheight         ! station elevation (m)

src/gsi/stpcalc.f90

@@ -844,10 +844,10 @@ subroutine stpcalc(stpinout,sval,sbias,dirx,dval,dbias, &
            if(outpen(i) < outpensave)then
               stp(ii)=outstp(i)
               outpensave=outpen(i)
-              istp_use=i
+              istp_use=ii
            end if
         end do
-        if(istp_use /= istp_iter) then
+        if(istp_use /= nsteptot) then
            final_ii=ii
            exit stepsize
         end if