Skip to content

Commit

Permalink
Adding the dry-bulb air temperature and tv flag in the observation di…
Browse files Browse the repository at this point in the history 
…agnostic file of conventional Q obs (for 3D RTMA run only) (NOAA-EMC#688)

- Motivation:
In the Auto-QC utility of 3D-RTMA, the dry-bulb air temperature (hereafter as T_dry) which is accompanied with moisture observation is required in the procedure of the quality control for conventional moisture observations (Q obs). Since this auto-qc of 3DRTMA retrieves many information from the observation diagnostic files (output of GSI run), it should be convenient for the sake of auto-qc of Q obs, if the required T_dry could be output into the observation diagnostics file of Q obs in the GSI run.

- Modifications to the code:
**_read_prepbufr.f90_**:
for qob and rtma run, add lines to use tpc to identify the sensible temp and virtual temp, and mark with tvflag, then get tdry with the temp obs in different ways w.r.t the associated tvflag;
**_setupq.f90_**:
add line to use l_rtma3d from module rapidrefresh_cldsurf_mod variable nreal (the length of obs diag information record for each obs) needs to be increased by 2 (for tdry and tvflag, when running GSI for 3DRTMA only)
in subroutine contents_binary_diag_, add line to put tdry into the array rdiagbuff (for binary format obsdaig file)
in subroutine contents_netcdf_diag_, add line to put tdry into metadata (for netcdf format obsdiag file)

This PR is to address the issue NOAA-EMC#666 : Adding the (calculated) dry-bulb temperature in the observation diagnostic file for conventional Q obs (only for 3D RTMA)
  • Loading branch information
@GangZhao-NOAA
GangZhao-NOAA committed Feb 13, 2024
1 parent bae0342 commit cd62003
Show file tree
Hide file tree
Showing 2 changed files with 68 additions and 3 deletions.
46 changes: 43 additions & 3 deletions src/gsi/read_prepbufr.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -149,8 +149,10 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
! 2020-05-04 wu - no rotate_wind for fv3_regional
! 2020-09-05 CAPS(C. Tong) - add flag for new vadwind obs to assimilate around the analysis time only
! 2023-03-23 draper - add code for processing T2m and q2m for global system
! 2023-07-30 Zhao - added code to extract obs of significant wave height (howvob) from bufr record
! 2023-07-30 zhao - added code to extract obs of significant wave height (howvob) from bufr record
! in prepbufr file for 3D analysis
! 2024-01-11 zhao - added code to extract sensible temp (tdry) and tv flag
! for moisture obs(qob) when running (2D/3D)RTMA

! input argument list:
! infile - unit from which to read BUFR data
Expand Down Expand Up @@ -225,6 +227,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
use adjust_cloudobs_mod, only: adjust_convcldobs,adjust_goescldobs
use mpimod, only: npe
use rapidrefresh_cldsurf_mod, only: i_gsdsfc_uselist,i_gsdqc,i_ens_mean
use rapidrefresh_cldsurf_mod, only: l_rtma3d
use gsi_io, only: verbose
use phil2, only: denest ! hilbert curve

Expand Down Expand Up @@ -390,6 +393,10 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
logical, allocatable,dimension(:) :: rusage,rthin
! end of block

! for extracting sensible-vs-virtual temp obs
integer(i_kind),dimension(1,255):: tqm4q
real(r_kind),dimension(1,255):: tvflg4q
real(r_double),dimension(1,255):: tobs4q

! equivalence to handle character names
equivalence(r_prvstg(1,1),c_prvstg)
Expand Down Expand Up @@ -876,6 +883,9 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
glcd=-999._r_double
endif

if(print_verbose) write(6,'(1x,A,A,A,2(A,1x,F8.3))') 'read_prepbufr:', &
trim(adjustl(obstype)),':', ' vtcd= ',vtcd,' glcd= ',glcd

call init_rjlists
call init_aircraft_rjlists
if(i_gsdsfc_uselist==1) call init_gsd_sfcuselist
Expand Down Expand Up @@ -1503,6 +1513,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
if(driftl)call ufbint(lunin,drfdat,8,255,iret,drift)

! raob level enhancement on temp and q obs
! (note: levs is increased by sonde_ext, and not same as original value read from prepbufr)
if(ext_sonde .and. kx==120) call sonde_ext(obsdat,tpc,qcmark,obserr,drfdat,levs,kx,vtcd)

nread=nread+levs
Expand Down Expand Up @@ -1713,6 +1724,25 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
end if
end if

! If moisture ob (qob) and (2D/3D)RTMA, set tv flag information (based on tpc)
! regarding virtual vs. sensible temperaure, to get tdry (if virtual temp
! then compute tdry; if sensible temp, then tdry= tsen), then save tdry
! in q-obsdaig file for RTMA offline Auto-QC.
if (qob .and. (l_rtma3d .or. twodvar_regional)) then
tobs4q(1,:) = bmiss
tqm4q(1,:) = bmiss
tvflg4q(1,:)= -one
do k=1,levs
tvflg4q(1,k)=one ! initialize as sensible
tobs4q(1,k)=obsdat(3,k) ! temp obs read in prepbufr
tqm4q(1,k)=tqm(k)
do j=1,20
if (tpc(k,j)==vtcd) tvflg4q(1,k)=zero ! reset flag if virtual
if (tpc(k,j)>=bmiss) exit ! end of stack
end do
end do
end if ! if qob & rtma

if(i_gsdqc==2) then
! AMV acceptance for all obs (E. James)
if (kx >= 240 .and. kx <= 260) then
Expand Down Expand Up @@ -2398,7 +2428,15 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
end if
qobcon=obsdat(2,k)*convert
tdry=r999
if (tqm(k)<lim_tqm) tdry=(obsdat(3,k)+t0c)/(one+fv*qobcon)
if (l_rtma3d .or. twodvar_regional) then
if (tvflg4q(1,k) == one) then ! tobs is sensible temp (tdry=tsen)
if (tqm4q(1,k)<lim_tqm) tdry= tobs4q(1,k)+t0c
else if (tvflg4q(1,k) == zero) then ! tobs is virtual temp (computing tdry with tv)
if (tqm4q(1,k)<lim_tqm) tdry= (tobs4q(1,k)+t0c)/(one+fv*qobcon)
end if
else
if (tqm(k)<lim_tqm) tdry=(obsdat(3,k)+t0c)/(one+fv*qobcon)
end if
cdata_all(1,iout)=qoe ! q error
cdata_all(2,iout)=dlon ! grid relative longitude
cdata_all(3,iout)=dlat ! grid relative latitude
Expand All @@ -2408,7 +2446,7 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
cdata_all(7,iout)=t4dv ! time
cdata_all(8,iout)=nc ! type
cdata_all(9,iout)=qmaxerr ! q max error
cdata_all(10,iout)=tdry ! dry temperature (obs is tv)
cdata_all(10,iout)=tdry ! dry temperature (obs is tv? No, depending on tvflg)
cdata_all(11,iout)=qqm(k) ! quality mark
cdata_all(12,iout)=obserr(2,k)*one_tenth ! original obs error
cdata_all(13,iout)=usage ! usage parameter
Expand All @@ -2426,6 +2464,8 @@ subroutine read_prepbufr(nread,ndata,nodata,infile,obstype,lunout,twindin,sis,&
if(perturb_obs)cdata_all(24,iout)=ran01dom()*perturb_fact ! q perturbation
if (twodvar_regional) &
call adjust_error(cdata_all(15,iout),cdata_all(16,iout),cdata_all(12,iout),cdata_all(1,iout))
if (l_rtma3d .or. twodvar_regional) &
cdata_all(25,iout)=tvflg4q(1,k) ! saving tv flag for q-obsdiag

! Total precipitable water (ssm/i)
else if(pwob) then
Expand Down
25 changes: 25 additions & 0 deletions src/gsi/setupq.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -113,6 +113,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
! for 3D-RTMA (if l_obsprvdiag is true).
! 2023-03-09 Draper added option to interpolate screen-level q from model 2m output.
! (hofx_2m_sfcfile)
! 2024-01-11 zhao - added tdry/tvflg in obs diagnostic files for (2D/3D)RTMA
!
!
! input argument list:
Expand Down Expand Up @@ -172,6 +173,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
use rapidrefresh_cldsurf_mod, only: l_sfcobserror_ramp_q
use rapidrefresh_cldsurf_mod, only: l_pbl_pseudo_surfobsq,pblh_ration,pps_press_incr, &
i_use_2mq4b,l_closeobs,i_coastline
use rapidrefresh_cldsurf_mod, only: l_rtma3d
use gsi_bundlemod, only : gsi_bundlegetpointer
use gsi_metguess_mod, only : gsi_metguess_get,gsi_metguess_bundle
use sparsearr, only: sparr2, new, size, writearray, fullarray
Expand Down Expand Up @@ -245,6 +247,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
integer(i_kind) ier,ilon,ilat,ipres,iqob,id,itime,ikx,iqmax,iqc
integer(i_kind) ier2,iuse,ilate,ilone,istnelv,iobshgt,izz,iprvd,isprvd
integer(i_kind) idomsfc,iderivative
integer(i_kind) iqt
integer(i_kind) ibb,ikk,idddd
real(r_kind) :: delz
type(sparr2) :: dhx_dx
Expand Down Expand Up @@ -335,6 +338,9 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
icat =22 ! index of data level category
ijb =23 ! index of non linear qc parameter
iptrb=24 ! index of q perturbation
if (l_rtma3d .or. twodvar_regional) then
iqt =25 ! index of flag indicating if virtual temp is associated to this moisture obs
end if

do i=1,nobs
muse(i)=nint(data(iuse,i)) <= jiter .and. nint(data(iqc,i)) < 8
Expand Down Expand Up @@ -416,6 +422,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
iip=0
nchar=1
ioff0=21
if (l_rtma3d .or. twodvar_regional) ioff0 = ioff0 + 2 ! 22:tdry; 23:tvflag (in binary obsdiag for 2D/3DRTMA)
nreal=ioff0
if (lobsdiagsave) nreal=nreal+4*miter+1
if (twodvar_regional .or. l_obsprvdiag) then
Expand Down Expand Up @@ -1232,6 +1239,11 @@ subroutine contents_binary_diag_(odiag)
rdiagbuf(20,ii) = qsges ! guess saturation specific humidity
rdiagbuf(21,ii) = 1e+10_r_single ! spread (filled in by EnKF)

if (l_rtma3d .or. twodvar_regional) then ! in binary obsdiag for 2D/3DRTMA
rdiagbuf(22,ii) = data(itemp,i) ! dry temperature associated to qob
rdiagbuf(23,ii) = data(iqt, i) ! tv flag (0: virtual temp; 1: sensible temp)
end if

ioff=ioff0
if (lobsdiagsave) then
do jj=1,miter
Expand Down Expand Up @@ -1309,6 +1321,11 @@ subroutine contents_binary_diagp_(odiag)
rdiagbufp(20,iip) = qsges ! guess saturation specific humidity
rdiagbufp(21,iip) = 1e+10_r_single ! spread (filled in by EnKF)

if (l_rtma3d .or. twodvar_regional) then ! in binary obsdiag for 2D/3DRTMA
rdiagbufp(22,ii) = data(itemp,i) ! dry temperature associated to qob
rdiagbufp(23,ii) = data(iqt, i) ! tv flag (0: virtual temp; 1: sensible temp)
end if

ioff=ioff0
!----
if (lobsdiagsave) then
Expand Down Expand Up @@ -1387,6 +1404,10 @@ subroutine contents_netcdf_diag_(odiag)
call nc_diag_metadata_to_single("Obs_Minus_Forecast_adjusted",ddiff )
call nc_diag_metadata_to_single("Obs_Minus_Forecast_unadjusted",qob,qges,'-')
call nc_diag_metadata_to_single("Forecast_Saturation_Spec_Hum",qsges )
if (l_rtma3d .or. twodvar_regional) then
call nc_diag_metadata_to_single("Observation_Tdry", data(itemp,i) )
call nc_diag_metadata_to_single("Setup_QC_Mark", data(iqt, i) )
endif
if (lobsdiagsave) then
do jj=1,miter
if (odiag%muse(jj)) then
Expand Down Expand Up @@ -1451,6 +1472,10 @@ subroutine contents_netcdf_diagp_(odiag)
call nc_diag_metadata_to_single("Obs_Minus_Forecast_adjusted",ddiff )
call nc_diag_metadata_to_single("Obs_Minus_Forecast_unadjusted",ddiff )
call nc_diag_metadata_to_single("Forecast_Saturation_Spec_Hum",qsges )
if (l_rtma3d .or. twodvar_regional) then
call nc_diag_metadata_to_single("Observation_Tdry", data(itemp,i) )
call nc_diag_metadata_to_single("Setup_QC_Mark", data(iqt, i) )
endif
!----
if (lobsdiagsave) then
do jj=1,miter
Expand Down

0 comments on commit cd62003

Please sign in to comment.