GitHub Issue NOAA-EMC#604 Undefined values found in radar reflectivit…

…y direct DA (NOAA-EMC#605)

**Description**

To prevent some undefined values found in the radar reflectivity direct
DA (if_model_dbz=T and l_use_dbz_directDA=F), corresponding parts are
fixed. It doesn't change the result except for the case of the execution
with the debug option.

Fixes NOAA-EMC#604

**Type of change**

- [x] Bug fix (non-breaking change which fixes an issue)

**How Has This Been Tested?**

The radar reflectivity DA test with the RRFS setting was done on Orion.
After this modification, EnVar was completed even with the debug option.
This modification didn't change the result in the test without the debug
option.
  
**Checklist**

- [x] My code follows the style guidelines of this project
- [x] I have performed a self-review of my own code
- [x] I have commented my code, particularly in hard-to-understand areas
- [x] New and existing tests pass with my changes
- [x] Any dependent changes have been merged and published

**Due date for this PR is 9/15/2023.** If this PR is not merged into
`develop` by this date, the PR will be closed and returned to the
developer.

Co-authored-by: Sho Yokota <[email protected]>

src/gsi/read_dbz_nc.f90

@@ -69,11 +69,12 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
   use kinds, only: r_kind,r_double,i_kind,r_single
   use constants, only: zero,half,one,two,deg2rad,rad2deg, &
                        one_tenth,r1000,r60,r60inv,r100,r400,grav_equator, &
-                       eccentricity,somigliana,grav_ratio,grav,semi_major_axis,flattening 
+                       eccentricity,somigliana,grav_ratio,grav,semi_major_axis,flattening,r_missing
   use gridmod, only: tll2xy,nsig,nlat,nlon
   use obsmod, only: iadate,doradaroneob,oneoblat,oneoblon,oneobheight, &
                     mintiltdbz,maxtiltdbz,minobrangedbz,maxobrangedbz,&
                     static_gsi_nopcp_dbz,rmesh_dbz,zmesh_dbz
+  use gsi_4dvar, only: iwinbgn
   use hybrid_ensemble_parameters,only : l_hyb_ens
   use obsmod,only: radar_no_thinning,missing_to_nopcp
   use convinfo, only: nconvtype,ctwind,icuse,ioctype
@@ -147,7 +148,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
   integer(i_kind) :: maxobs,nchanl,ilat,ilon,scount
 
   real(r_kind) :: thistiltr,thisrange,this_stahgt,thishgt                           
-  real(r_kind) :: thisazimuthr,t4dv, &
+  real(r_kind) :: thisazimuthr, &
                   dlat,dlon,thiserr,thislon,thislat, &
                   timeb
   real(r_kind) :: radartwindow
@@ -337,6 +338,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
 
   call w3fs21(iadate,mins_an)  !mins_an -integer number of mins snce 01/01/1978
   rmins_an=mins_an             !convert to real number
+  timeb=real(mins_an-iwinbgn,r_kind)  !assume all observations are at the analysis time
 
   ivar = 1
 
@@ -453,7 +455,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
 
 
            ntmp=ndata  ! counting moved to map3gridS
-           timedif=abs(t4dv) !don't know about this
+           timedif=zero  ! assume all observations are at the analysis time
            crit1 = timedif/r6+half
 
            call map3grids(1,zflag,zl_thin,nlevz,thislat,thislon,&
@@ -481,7 +483,10 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
 
         !!end modified for thinning
 
-         thisazimuthr=0.0_r_kind
+         thisazimuthr=r_missing
+         thistiltr=r_missing
+         this_stahgt=r_missing
+         thisrange=r_missing
          this_staid=radid                !Via equivalence in declaration, value is propagated
                                                            !  to rstation_id used below.
          cdata_all(1,iout) = thiserr                       ! reflectivity obs error (dB) - inflated/adjusted
@@ -491,7 +496,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
          cdata_all(5,iout) = dbzQC(i,j,k)                      ! radar reflectivity factor 
          cdata_all(6,iout) = thisazimuthr                  ! 90deg-azimuth angle (radians)
 
-         cdata_all(7,iout) = timeb*r60inv                  ! obs time (analyis relative hour)
+         cdata_all(7,iout) = timeb*r60inv                  ! obs time (relative hour from beginning of the DA window)
          cdata_all(8,iout) = ikx                           ! type                  
          cdata_all(9,iout) = thistiltr                     ! tilt angle (radians)
          cdata_all(10,iout)= this_stahgt                   ! station elevation (m)
@@ -521,7 +526,6 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
   !---all looping done now print diagnostic output
 
   write(6,*)'READ_dBZ: Reached eof on radar reflectivity file'
-  write(6,*)'READ_dBZ: # volumes in input file             =',nvol
   write(6,*)'READ_dBZ: # read in obs. number               =',nread
   write(6,*)'READ_dBZ: # elevations outside time window    =',numbadtime
   write(6,*)'READ_dBZ: # of noise obs to no precip obs     =',num_nopcp

src/gsi/setupdbz.f90

@@ -590,14 +590,17 @@ subroutine setupdbz(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,radardbz_d
 !    Compute observation pressure (only used for diagnostics)
      dz     = zges(k2)-zges(k1)
      dlnp   = prsltmp(k2)-prsltmp(k1)
-     pobl   = prsltmp(k1) + (dlnp/dz)*(zob-zges(k1))
-
-     presw  = ten*exp(pobl)
-
      if ( l_use_dbz_directDA ) then
-        presq  = presw         
+        pobl   = prsltmp(k1) + (dlnp/dz)*(zob-zges(k1))
+        presw  = ten*exp(pobl)
+        presq  = presw
      else
-        if( (k1 == k2) .and. (k1 == 1) ) presw=ten*exp(prsltmp(k1)) 
+        if( (k1 == k2) .and. (k1 == 1) ) then
+           presw  = ten*exp(prsltmp(k1))
+        else
+           pobl   = prsltmp(k1) + (dlnp/dz)*(zob-zges(k1))
+           presw  = ten*exp(pobl)
+        end if
      end if
 
 !    solution to Nan in some members only for EnKF which causes problem?