Skip to content

Commit

Permalink
Merge remote-tracking branch 'origin/develop' into feature/ss_160
Browse files Browse the repository at this point in the history
  • Loading branch information
DavidHuber-NOAA committed Feb 13, 2024
2 parents e56921c + 86ad20e commit 3951289
Show file tree
Hide file tree
Showing 3 changed files with 69 additions and 4 deletions.
2 changes: 1 addition & 1 deletion fix
46 changes: 43 additions & 3 deletions src/gsi/read_prepbufr.f90
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
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 3951289

Please sign in to comment.