Skip to content

Commit

Permalink
Add observation provider info into obs diagnostic file for ps/w/spd/q…
Browse files Browse the repository at this point in the history 
…/t (#8)

* Adding modifications to obs setup code(setup*.f90) for ps/q/t/w/spd, to include the observation provider and
sub-provider information into the diag files, whichwill be used in read_diag_conv_3drtma in Auto-ObsQC 

* Also adding a new spcific code read_diag_conv_3drtma.f90 under read_diag to read the osb-diag file which includes
the obs provider information for auto-qc.
  • Loading branch information
@GangZhao-NOAA
GangZhao-NOAA committed Dec 7, 2021
1 parent 2c8edc5 commit 5f3a9b9
Show file tree
Hide file tree
Showing 8 changed files with 395 additions and 70 deletions.
9 changes: 8 additions & 1 deletion src/gsi/gsimod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -35,6 +35,7 @@ module gsimod
use obsmod, only: luse_obsdiag
use obsmod, only: netcdf_diag, binary_diag
use obsmod, only: l_wcp_cwm
use obsmod, only: l_obsprvdiag
use obsmod, only: aircraft_recon, &

! The following variables are the coefficients that describe
Expand Down Expand Up @@ -457,6 +458,9 @@ module gsimod
! reflectivity ghost in analysis. (default is 0)
! 2021-09-08 Guoqing - Add npePgrp_rfv3,rfv3_pe_T,rfv3_pe_v,rfv3_pe_q,rfv3_pe_ps,rfv3_pe_dz
! To speed up writing out regional FV3 final analysis.
! 2021-11-16 Zhao - add option l_obsprvdiag (if true) to trigger the output of
! observation provider and sub-provider information into
! obsdiags files (used for AutoObsQC)
!
!EOP
!-------------------------------------------------------------------------
Expand Down Expand Up @@ -668,6 +672,9 @@ module gsimod
! npePgrp_rfv3 - The number of PEs in each write group for regional fv3
! rfv3_pe_T,rfv3_pe_v,rfv3_pe_q,rfv3_pe_ps,rfv3_pe_dz - specify starting pe for T,U/V,Q,PS,DZ write groups
!
! l_obsprvdiag - trigger (if true) writing out observation provider and sub-provider
! information into obsdiags files (used for AutoObsQC)
!
! NOTE: for now, if in regional mode, then iguess=-1 is forced internally.
! add use of guess file later for regional mode.

Expand Down Expand Up @@ -712,7 +719,7 @@ module gsimod
write_fv3_incr,incvars_to_zero,incvars_zero_strat,incvars_efold,diag_version,&
cao_check,lcalc_gfdl_cfrac,tau_fcst,efsoi_order,lupdqc,lqcoef,cnvw_option,l2rwthin,hurricane_radar,&
npePgrp_rfv3,rfv3_pe_T,rfv3_pe_v,rfv3_pe_q,rfv3_pe_ps,rfv3_pe_dz, &
l_reg_update_hydro_delz
l_reg_update_hydro_delz, l_obsprvdiag

! GRIDOPTS (grid setup variables,including regional specific variables):
! jcap - spectral resolution
Expand Down
10 changes: 10 additions & 0 deletions src/gsi/obsmod.F90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -156,6 +156,9 @@ module obsmod
! observation error (DOE) specification to
! GSI namelist level.
! 2020-09-15 Wu - add option tcp_posmatch to mitigate possibility of erroneous TC initialization
! 2021-11-16 Zhao - add option l_obsprvdiag (if true) to trigger the output of
! observation provider and sub-provider information into
! obsdiags files (used for AutoObsQC)
!
! Subroutines Included:
! sub init_obsmod_dflts - initialize obs related variables to default values
Expand Down Expand Up @@ -397,6 +400,7 @@ module obsmod
! (nobs_type,npe)
! def binary_diag - trigger binary diag-file output (being phased out)
! def netcdf_diag - trigger netcdf diag-file output
! def l_obsprvdiag - trigger obs provider info output into obsdiags files
! def l_wcp_cwm - namelist logical whether to use operator that
! includes cwm for both swcp and lwcp or not
! def neutral_stability_windfact_2dvar - logical, if .true., then use simple formula representing
Expand Down Expand Up @@ -483,6 +487,7 @@ module obsmod
public :: nobs_sub

public :: netcdf_diag, binary_diag
public :: l_obsprvdiag

public :: l_wcp_cwm
public :: aircraft_recon
Expand Down Expand Up @@ -547,6 +552,7 @@ module obsmod

logical luse_obsdiag
logical binary_diag, netcdf_diag
logical l_obsprvdiag

! Declare types

Expand Down Expand Up @@ -707,6 +713,7 @@ subroutine init_obsmod_dflts
! 2015-07-10 pondeca - add cldch
! 2015-10-27 todling - default to luse_obsdiag is true now
! 2016-03-07 pondeca - add uwnd10m,vwnd10m
! 2021-11-16 zhao - add initialization of l_obsprvdiag (.FALSE. as default)
!
! input argument list:
!
Expand Down Expand Up @@ -899,6 +906,9 @@ subroutine init_obsmod_dflts
netcdf_diag = .false. ! by default, do not write netcdf_diag
binary_diag = .true. ! by default, do write binary diag

! set default on triggering the output of obs provider info into obsdiags file
l_obsprvdiag = .false. ! by default, do not write obs provider info

l_wcp_cwm = .false. ! .true. = use operator that involves cwm
aircraft_recon = .false. ! .true. = use DOE for aircraft data
hurricane_radar = .false. ! .true. = use radar data for hurricane application
Expand Down
19 changes: 14 additions & 5 deletions src/gsi/setupps.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -81,6 +81,11 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa
! 2017-03-31 Hu - addd option l_closeobs to use closest obs to analysis
! time in analysis
! 2019-09-20 Su - remove current VQC part and add subroutine call on VQC and add new VQC option
! 2021-10-xx pondeca/morris/zhao - added observation provider/subprovider
! information in diagonostic file, which is used
! in offline observation quality control program (AutoObsQC)
! for 3D-RTMA (if l_obsprvdiag is true).
!
!
! input argument list:
! lunin - unit from which to read observations
Expand Down Expand Up @@ -118,6 +123,7 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa
use gsi_4dvar, only: nobs_bins,hr_obsbin,min_offset
use oneobmod, only: magoberr,maginnov,oneobtest
use obsmod, only: netcdf_diag, binary_diag, dirname
use obsmod, only: l_obsprvdiag
use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, &
nc_diag_write, nc_diag_data2d
use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_get_dim, nc_diag_read_close
Expand Down Expand Up @@ -321,7 +327,10 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa
ioff0=20
nreal=ioff0
if (lobsdiagsave) nreal=nreal+4*miter+1
if (twodvar_regional) then; nreal=nreal+2; allocate(cprvstg(nobs),csprvstg(nobs)); endif
if (twodvar_regional .or. l_obsprvdiag) then
nreal=nreal+2 !account for idomsfc,izz
allocate(cprvstg(nobs),csprvstg(nobs)) !obs provider info
endif
if (save_jacobian) then
nnz = 1 ! number of non-zero elements in dH(x)/dx profile
nind = 1
Expand Down Expand Up @@ -667,13 +676,13 @@ subroutine setupps(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsa
if(binary_diag .and. ii>0)then
write(7)' ps',nchar,nreal,ii,mype,ioff0
write(7)cdiagbuf(1:ii),rdiagbuf(:,1:ii)
deallocate(cdiagbuf,rdiagbuf)

if (twodvar_regional) then
if (twodvar_regional .or. l_obsprvdiag) then
write(7)cprvstg(1:ii),csprvstg(1:ii)
deallocate(cprvstg,csprvstg)
endif
end if
deallocate(cdiagbuf,rdiagbuf)
end if

! End of routine
Expand Down Expand Up @@ -856,7 +865,7 @@ subroutine contents_binary_diag_(odiag)
enddo
endif

if (twodvar_regional) then
if (twodvar_regional .or. l_obsprvdiag) then
ioff = ioff + 1
rdiagbuf(ioff,ii) = data(idomsfc,i) ! dominate surface type
ioff = ioff + 1
Expand Down Expand Up @@ -923,7 +932,7 @@ subroutine contents_netcdf_diag_(odiag)
call nc_diag_data2d("ObsDiagSave_obssen", odiag%obssen )
endif

if (twodvar_regional) then
if (twodvar_regional .or. l_obsprvdiag) then
call nc_diag_metadata("Dominant_Sfc_Type", data(idomsfc,i) )
call nc_diag_metadata("Model_Terrain", data(izz,i) )
r_prvstg = data(iprvd,i)
Expand Down
112 changes: 97 additions & 15 deletions src/gsi/setupq.f90
View file
Original file line number Diff line number Diff line change
Expand Up @@ -107,6 +107,10 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
! error (DOE) calculation to the namelist
! level; they are now loaded by
! aircraftinfo.
! 2021-10-xx pondeca/morris/zhao - added observation provider/subprovider
! information in diagonostic file, which is used
! in offline observation quality control program (AutoObsQC)
! for 3D-RTMA (if l_obsprvdiag is true).
!
!
! input argument list:
Expand Down Expand Up @@ -144,6 +148,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
use m_obsLList, only: obsLList
use obsmod, only: luse_obsdiag,ianldate
use obsmod, only: netcdf_diag, binary_diag, dirname
use obsmod, only: l_obsprvdiag
use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, &
nc_diag_write, nc_diag_data2d
use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_get_dim, nc_diag_read_close
Expand Down Expand Up @@ -245,6 +250,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
character(8) station_id
character(8),allocatable,dimension(:):: cdiagbuf,cdiagbufp
character(8),allocatable,dimension(:):: cprvstg,csprvstg
character(8),allocatable,dimension(:):: cprvstgp,csprvstgp ! <-- provider info array for pseudo obs
character(8) c_prvstg,c_sprvstg
real(r_double) r_prvstg,r_sprvstg

Expand Down Expand Up @@ -394,7 +400,11 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
ioff0=21
nreal=ioff0
if (lobsdiagsave) nreal=nreal+4*miter+1
if (twodvar_regional) then; nreal=nreal+2; allocate(cprvstg(nobs),csprvstg(nobs)); endif
if (twodvar_regional .or. l_obsprvdiag) then
nreal=nreal+2 ! account for idomsfc,izz
allocate(cprvstg(nobs),csprvstg(nobs)) ! obs provider info
if(l_pbl_pseudo_surfobsq) allocate(cprvstgp(nobs*3),csprvstgp(nobs*3)) ! obs provider info for pseudo obs
endif
if (save_jacobian) then
nnz = 2 ! number of non-zero elements in dH(x)/dx profile
nind = 1
Expand Down Expand Up @@ -442,6 +452,7 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
dpres=data(ipres,i)

rmaxerr=data(iqmax,i)
rstation_id = data(id,i)
error=data(ier2,i)
prest=r10*exp(dpres) ! in mb
var_jb=data(ijb,i)
Expand Down Expand Up @@ -915,11 +926,11 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
if (err_adjst>tiny_r_kind) errinv_adjst = one/err_adjst
if (err_final>tiny_r_kind) errinv_final = one/err_final

if(binary_diag) call contents_binary_diagp_()
if(binary_diag) call contents_binary_diagp_(my_diag_pbl)
else
iip=3*nobs
endif
if(netcdf_diag) call contents_netcdf_diagp_()
if(netcdf_diag) call contents_netcdf_diagp_(my_diag_pbl)
endif !conv_diagsave .and. luse(i))

prest = prest - pps_press_incr
Expand All @@ -940,20 +951,25 @@ subroutine setupq(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav
if(conv_diagsave) then
if(netcdf_diag) call nc_diag_write
if(binary_diag .and. ii>0)then
write(7)' q',nchar,nreal,ii+iip,mype,ioff0
write(7)' q',nchar,nreal,ii+iip,mype,ioff0,iip
if(l_pbl_pseudo_surfobsq .and. iip>0) then
write(7)cdiagbuf(1:ii),cdiagbufp(1:iip),rdiagbuf(:,1:ii),rdiagbufp(:,1:iip)
deallocate(cdiagbufp,rdiagbufp)
else
write(7)cdiagbuf(1:ii),rdiagbuf(:,1:ii)
endif
deallocate(cdiagbuf,rdiagbuf)

if (twodvar_regional) then
write(7)cprvstg(1:ii),csprvstg(1:ii)
if (twodvar_regional .or. l_obsprvdiag) then
if(l_pbl_pseudo_surfobsq .and. iip>0) then
write(7)cprvstg(1:ii),cprvstgp(1:iip),csprvstg(1:ii),csprvstgp(1:iip)
else
write(7)cprvstg(1:ii),csprvstg(1:ii)
endif
deallocate(cprvstg,csprvstg)
if(l_pbl_pseudo_surfobsq) deallocate(cprvstgp,csprvstgp)
endif
endif
deallocate(cdiagbuf,rdiagbuf)
if(l_pbl_pseudo_surfobsq) deallocate(cdiagbufp,rdiagbufp)
end if

! End of routine
Expand Down Expand Up @@ -1143,7 +1159,7 @@ subroutine contents_binary_diag_(odiag)
enddo
endif

if (twodvar_regional) then
if (twodvar_regional .or. l_obsprvdiag) then
ioff = ioff + 1
rdiagbuf(ioff,ii) = data(idomsfc,i) ! dominate surface type
ioff = ioff + 1
Expand All @@ -1160,12 +1176,14 @@ subroutine contents_binary_diag_(odiag)

end subroutine contents_binary_diag_

subroutine contents_binary_diagp_
subroutine contents_binary_diagp_(odiag)
type(obs_diag),pointer,intent(in):: odiag

cdiagbufp(iip) = station_id ! station id

rdiagbufp(1,iip) = ictype(ikx) ! observation type
rdiagbufp(2,iip) = icsubtype(ikx) ! observation subtype
! rdiagbufp(2,iip) = icsubtype(ikx) ! observation subtype
rdiagbufp(2,iip) = -1 ! observation subtype (-1 for pseudo obs sub-type)

rdiagbufp(3,iip) = data(ilate,i) ! observation latitude (degrees)
rdiagbufp(4,iip) = data(ilone,i) ! observation longitude (degrees)
Expand Down Expand Up @@ -1196,8 +1214,43 @@ subroutine contents_binary_diagp_
rdiagbufp(21,iip) = 1e+10_r_single ! spread (filled in by EnKF)

ioff=ioff0
!----
if (lobsdiagsave) then
do jj=1,miter
ioff=ioff+1
if (odiag%muse(jj)) then
rdiagbufp(ioff,iip) = one
else
rdiagbufp(ioff,iip) = -one
endif
enddo
do jj=1,miter+1
ioff=ioff+1
rdiagbufp(ioff,iip) = odiag%nldepart(jj)
enddo
do jj=1,miter
ioff=ioff+1
rdiagbufp(ioff,iip) = odiag%tldepart(jj)
enddo
do jj=1,miter
ioff=ioff+1
rdiagbufp(ioff,iip) = odiag%obssen(jj)
enddo
endif

if (twodvar_regional .or. l_obsprvdiag) then
ioff = ioff + 1
rdiagbufp(ioff,iip) = -9999. ! data(idomsfc,i) ! dominate surface type
ioff = ioff + 1
rdiagbufp(ioff,iip) = -9999. ! data(izz,i) ! model terrain at ob location
r_prvstg = data(iprvd,i)
cprvstgp(iip) = '88888888' !c_prvstg ! provider name
r_sprvstg = data(isprvd,i)
csprvstgp(iip) = '88888888' !c_sprvstg ! subprovider name
endif
!----
if (save_jacobian) then
call writearray(dhx_dx, rdiagbuf(ioff+1:nreal,ii))
call writearray(dhx_dx, rdiagbufp(ioff+1:nreal,iip))
ioff = ioff + size(dhx_dx)
endif

Expand Down Expand Up @@ -1253,7 +1306,7 @@ subroutine contents_netcdf_diag_(odiag)
call nc_diag_data2d("ObsDiagSave_obssen", odiag%obssen )
endif

if (twodvar_regional) then
if (twodvar_regional .or. l_obsprvdiag) then
call nc_diag_metadata("Dominant_Sfc_Type", data(idomsfc,i) )
call nc_diag_metadata("Model_Terrain", data(izz,i) )
r_prvstg = data(iprvd,i)
Expand All @@ -1271,14 +1324,17 @@ subroutine contents_netcdf_diag_(odiag)

end subroutine contents_netcdf_diag_

subroutine contents_netcdf_diagp_
subroutine contents_netcdf_diagp_(odiag)
type(obs_diag),pointer,intent(in):: odiag
! Observation class
character(7),parameter :: obsclass = ' q'
real(r_kind),dimension(miter) :: obsdiag_iuse

call nc_diag_metadata("Station_ID", station_id )
call nc_diag_metadata("Observation_Class", obsclass )
call nc_diag_metadata("Observation_Type", ictype(ikx) )
call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) )
! call nc_diag_metadata("Observation_Subtype", icsubtype(ikx) )
call nc_diag_metadata("Observation_Subtype", -1 ) ! (-1 for pseudo obs sub-type)
call nc_diag_metadata("Latitude", sngl(data(ilate,i)) )
call nc_diag_metadata("Longitude", sngl(data(ilone,i)) )
call nc_diag_metadata("Station_Elevation", sngl(data(istnelv,i)) )
Expand All @@ -1304,7 +1360,33 @@ subroutine contents_netcdf_diagp_
call nc_diag_metadata("Forecast_adjusted", sngl(data(iqob,i)-ddiff))
call nc_diag_metadata("Forecast_unadjusted", sngl(data(iqob,i)-ddiff))
call nc_diag_metadata("Forecast_Saturation_Spec_Hum", sngl(qsges) )
!----
if (lobsdiagsave) then
do jj=1,miter
if (odiag%muse(jj)) then
obsdiag_iuse(jj) = one
else
obsdiag_iuse(jj) = -one
endif
enddo

call nc_diag_data2d("ObsDiagSave_iuse", obsdiag_iuse )
call nc_diag_data2d("ObsDiagSave_nldepart", odiag%nldepart )
call nc_diag_data2d("ObsDiagSave_tldepart", odiag%tldepart )
call nc_diag_data2d("ObsDiagSave_obssen", odiag%obssen )
endif

if (twodvar_regional .or. l_obsprvdiag) then
call nc_diag_metadata("Dominant_Sfc_Type", data(idomsfc,i) )
call nc_diag_metadata("Model_Terrain", data(izz,i) )
r_prvstg = data(iprvd,i)
! call nc_diag_metadata("Provider_Name", c_prvstg )
call nc_diag_metadata("Provider_Name", "88888888" )
r_sprvstg = data(isprvd,i)
! call nc_diag_metadata("Subprovider_Name", c_sprvstg )
call nc_diag_metadata("Subprovider_Name", "88888888" )
endif
!----
if (save_jacobian) then
call nc_diag_data2d("Observation_Operator_Jacobian_stind", dhx_dx%st_ind)
call nc_diag_data2d("Observation_Operator_Jacobian_endind", dhx_dx%end_ind)
Expand Down
Loading

0 comments on commit 5f3a9b9

Please sign in to comment.