Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

GitHub Issue NOAA-EMC/GSI#604 Undefined values found in radar reflectivity direct DA #605

Merged
merged 1 commit into from
Sep 30, 2023
Merged

GitHub Issue NOAA-EMC/GSI#604 Undefined values found in radar reflectivity direct DA #605

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
16 changes: 10 additions & 6 deletions src/gsi/read_dbz_nc.f90
View file
Open in desktop
Copy link
Contributor

@TingLei-NOAA TingLei-NOAA Aug 9, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@shoyokota I think t4dv is also undefined as shown below. Is this an issue when you ran your case with debug mode gsi? That will be great if it is also taken care of in this PR. From my point of view, it think it could be defined as timeb+ (analysis time - the beginning time of the time window). You could see how it is set in other observation 's setup subroutine
image

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@TingLei-NOAA This t4dv is only used in thinning independently in each timeslot. However, it is assumed that all the inputted reflectivity observations were observed at the analysis time here. So, we don't have to input any specific times as t4dv here, and the result doesn't depend on t4dv as long as it is the same in all reflectivity observations. Here, I modified here to timedif=zero simply to prevent undefined values.

Copy link
Contributor

@TingLei-NOAA TingLei-NOAA Aug 16, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@shoyokota if you runs doesn't go through this part. We could let it be the current codes (with that comment) for future developers to take care . It should be analysis time - the begin of the time window. We don't need to enforce this part also work with that all obs are being at the analysis time . And, even with assumption/requirement that dbz obs being on the analysis time, the fgat/4dvar could still be invoked for other observations could be distributed on different time levels. So, I prefer leaving it to be updated in the future rather than making a change with unnecessary assumptions.

Copy link
Collaborator Author

@shoyokota shoyokota Aug 16, 2023 •
edited
Loading

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

This part is computed in the case of radar_no_thinning=F even now, so I think some kind of fix is needed. Is it better to add comments such as 'observation time (hour) should be input for 4D observations in 3D thinning, but 3D observations here' at the right of 'timedif=zero'?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@shoyokota Thanks for the clarification. I will add my thoughts on those specific codes respectively.

Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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
Expand Down Expand Up @@ -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

Expand Down Expand Up @@ -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,&
Expand Down Expand Up @@ -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
Expand All @@ -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)
Expand Down Expand Up @@ -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
Expand Down
15 changes: 9 additions & 6 deletions src/gsi/setupdbz.f90
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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?
Expand Down