From 6e9da0ed17183d7041741f2458b673be44dd1b47 Mon Sep 17 00:00:00 2001 From: Sho Yokota 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 | 9 +++++---- src/gsi/setupdbz.f90 | 15 +++++++++------ 2 files changed, 14 insertions(+), 10 deletions(-) diff --git a/src/gsi/read_dbz_nc.f90 b/src/gsi/read_dbz_nc.f90 index cddbd14de4..6c31caae75 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) !assume all observations are at the analysis time 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 ! assume all observations are at the analysis time crit1 = timedif/r6+half call map3grids(1,zflag,zl_thin,nlevz,thislat,thislon,& @@ -490,7 +492,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) @@ -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?