Xu Lu ([email protected]) and Biju Thomas ([email protected]) fixed …

…bugs regarding uninitialized variables (e.g. toff) and dual_res in GSI.

        1. In read_radar.f90, uninitialized toff is making all the ground-based radar observations to be placed at -3h instead of 0h, which create wrong increments for FGAT and 4DEnVar.
        2. In read_radar.f90, unintialized 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
        5. In read_fl_hdob.f90, dlnpsob is replaced with 1000. since the SFMR does not sampling surface pressure, and the uninitiailzed dlnpsob create issues later in setupspd.f90.
        6. In mod_fv3_lola.f90, (i,j+1) should be used instead of (i+1,j) in searching for V edges

src/gsi/mod_fv3_lola.f90

@@ -955,7 +955,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))
@@ -983,12 +982,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

@@ -1121,7 +1121,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)=1000.                  ! ln(surface pressure in cb) !Since dlnpsob is not provided by SFMR, force it to be 1000hPa. 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

@@ -3226,8 +3226,9 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs)
   use gridmod, only: regional,nlat,nlon,tll2xy,rlats,rlons,rotate_wind_ll2xy,&
            fv3_regional
   use convinfo, only: nconvtype,ncmiter,ncgroup,ncnumgrp,icuse,ioctype
-  use deter_sfc_mod, only: deter_sfc2
+  use deter_sfc_mod, only: deter_sfc2,deter_zsfc_model
   use mpimod, only: npe
+  use obsmod, only: time_offset
 
   implicit none
 
@@ -3268,7 +3269,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
@@ -3407,7 +3408,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)
 
@@ -3526,15 +3527,16 @@ subroutine read_radar_l2rw_novadqc(ndata,nodata,lunout,obstype,sis,nobs)
            if(mod(ndata,ncnumgrp(ikx))== ncgroup(ikx)-1)usage=ncmiter(ikx)
         end if
 
-        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)
@@ -3598,10 +3600,10 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
   use oneobmod, only: oneobtest,learthrel_rw
   use gsi_4dvar, only: l4dvar,l4densvar,winlen,time_4dvar
   use gridmod, only: regional,nlat,nlon,tll2xy,rlats,rlons,rotate_wind_ll2xy,nsig 
-  use obsmod, only: doradaroneob,oneobradid,time_offset 
+  use obsmod, only: doradaroneob,oneobradid,time_offset
   use mpeu_util, only: gettablesize,gettable 
   use convinfo, only: nconvtype,icuse,ioctype
-  use deter_sfc_mod, only: deter_sfc2
+  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  
   use constants, only: eccentricity,somigliana,grav_ratio,grav,semi_major_axis,flattening,grav_equator 
@@ -3646,7 +3648,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,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
@@ -3973,7 +3975,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
@@ -4070,7 +4072,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(1,zflag,zl_thin,nlevz,dlat_earth,dlon_earth,&
@@ -4140,6 +4142,7 @@ subroutine read_radar_l2rw(ndata,nodata,lunout,obstype,sis,nobs,hgtl_full)
 
           usage = zero
           if(icuse(ikx) < 0)usage=r100
+          call deter_zsfc_model(dlat,dlon,zsges)
 
           nsuper2_kept=nsuper2_kept+1
           cdata(1) = error             ! wind obs error (m/s)
@@ -4148,7 +4151,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)