From 664d0726034b34df03ce4c54c79ba7f9f3cb6c41 Mon Sep 17 00:00:00 2001
From: Sho Yokota <syokota@Orion-login-2.HPC.MsState.Edu>
Date: Fri, 4 Aug 2023 16:00:16 -0500
Subject: [PATCH] GitHub Issue NOAA-EMC/GSI#604 Undefined values found in radar
 reflectivity direct DA

---
 src/gsi/read_dbz_nc.f90 |  7 ++++---
 src/gsi/setupdbz.f90    | 15 +++++++++------
 2 files changed, 13 insertions(+), 9 deletions(-)

diff --git a/src/gsi/read_dbz_nc.f90 b/src/gsi/read_dbz_nc.f90
index cddbd14de4..14fa322af7 100644
--- a/src/gsi/read_dbz_nc.f90
+++ b/src/gsi/read_dbz_nc.f90
@@ -74,6 +74,7 @@ subroutine read_dbz_nc(nread,ndata,nodata,infile,lunout,obstype,sis,hgtl_full,no
   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)
  
   ivar = 1
   
@@ -452,7 +454,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
            crit1 = timedif/r6+half
      
            call map3grids(1,zflag,zl_thin,nlevz,thislat,thislon,&
@@ -520,7 +522,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
diff --git a/src/gsi/setupdbz.f90 b/src/gsi/setupdbz.f90
index 96f0378c52..043adc112a 100644
--- a/src/gsi/setupdbz.f90
+++ b/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?