From 5f3a9b9215abdf28c9cbdc607f8d3088d7c19b8f Mon Sep 17 00:00:00 2001 From: Gang <53267411+GangZhao-NOAA@users.noreply.github.com> Date: Tue, 7 Dec 2021 12:48:32 -0500 Subject: [PATCH] Add observation provider info into obs diagnostic file for ps/w/spd/q/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. --- src/gsi/gsimod.F90 | 9 +- src/gsi/obsmod.F90 | 10 ++ src/gsi/setupps.f90 | 19 +- src/gsi/setupq.f90 | 112 ++++++++++-- src/gsi/setupspd.f90 | 18 +- src/gsi/setupt.f90 | 117 +++++++++++-- src/gsi/setupw.f90 | 18 +- .../read_diag/read_diag_conv.f90 | 162 +++++++++++++++--- 8 files changed, 395 insertions(+), 70 deletions(-) diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90 index 4b74ac9149..d6b2efae33 100644 --- a/src/gsi/gsimod.F90 +++ b/src/gsi/gsimod.F90 @@ -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 @@ -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 !------------------------------------------------------------------------- @@ -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. @@ -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 diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90 index 8d46a51d2a..9f3c5f766e 100644 --- a/src/gsi/obsmod.F90 +++ b/src/gsi/obsmod.F90 @@ -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 @@ -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 @@ -483,6 +487,7 @@ module obsmod public :: nobs_sub public :: netcdf_diag, binary_diag + public :: l_obsprvdiag public :: l_wcp_cwm public :: aircraft_recon @@ -547,6 +552,7 @@ module obsmod logical luse_obsdiag logical binary_diag, netcdf_diag + logical l_obsprvdiag ! Declare types @@ -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: ! @@ -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 diff --git a/src/gsi/setupps.f90 b/src/gsi/setupps.f90 index 8a74d7d9a4..04344a2f7e 100644 --- a/src/gsi/setupps.f90 +++ b/src/gsi/setupps.f90 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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 @@ -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) diff --git a/src/gsi/setupq.f90 b/src/gsi/setupq.f90 index f1731fe848..c98a52da8e 100644 --- a/src/gsi/setupq.f90 +++ b/src/gsi/setupq.f90 @@ -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: @@ -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 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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 @@ -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 @@ -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) @@ -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 @@ -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) @@ -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)) ) @@ -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) diff --git a/src/gsi/setupspd.f90 b/src/gsi/setupspd.f90 index 1e740fdd25..91b2467bf3 100644 --- a/src/gsi/setupspd.f90 +++ b/src/gsi/setupspd.f90 @@ -79,6 +79,10 @@ subroutine setupspd(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diags ! 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: ! lunin - unit from which to read observations @@ -108,6 +112,7 @@ subroutine setupspd(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diags lobsdiagsave,nobskeep,lobsdiag_allocated,time_offset,& lobsdiag_forenkf,aircraft_recon use obsmod, only: netcdf_diag, binary_diag, dirname, ianldate + 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 @@ -280,7 +285,10 @@ subroutine setupspd(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diags 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 + endif if (save_jacobian) then nnz = 4 ! number of non-zero elements in dH(x)/dx profile nind = 2 @@ -684,13 +692,13 @@ subroutine setupspd(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diags if(binary_diag .and. ii>0)then write(7)'spd',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 @@ -915,7 +923,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 @@ -980,7 +988,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) diff --git a/src/gsi/setupt.f90 b/src/gsi/setupt.f90 index 8891fa5c12..19c6b6ec45 100644 --- a/src/gsi/setupt.f90 +++ b/src/gsi/setupt.f90 @@ -40,6 +40,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav use gsi_4dvar, only: nobs_bins,hr_obsbin,min_offset 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 @@ -69,7 +70,7 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav use converr, only: ptabl use rapidrefresh_cldsurf_mod, only: l_gsd_terrain_match_surftobs,l_sfcobserror_ramp_t use rapidrefresh_cldsurf_mod, only: l_pbl_pseudo_surfobst, pblh_ration,pps_press_incr - use rapidrefresh_cldsurf_mod, only: i_use_2mt4b,i_sfct_gross,l_closeobs,i_coastline + use rapidrefresh_cldsurf_mod, only: i_use_2mt4b,i_sfct_gross,l_closeobs,i_coastline use aircraftinfo, only: npredt,predt,aircraft_t_bc_pof,aircraft_t_bc, & aircraft_t_bc_ext,ostats_t,rstats_t,upd_pred_t @@ -221,6 +222,10 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! observation error (DOE) calculation to ! the namelist level; they are now ! loaded by obsmod. +! 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). ! ! !REMARKS: ! language: f90 @@ -296,6 +301,7 @@ subroutine setupt(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 @@ -469,7 +475,11 @@ subroutine setupt(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav nreal=nreal+npredt+2 idia0=nreal 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 used in diag for RTMA + allocate(cprvstg(nobs),csprvstg(nobs)) ! provider/subprovider info + if(l_pbl_pseudo_surfobst) allocate(cprvstgp(nobs*3),csprvstgp(nobs*3)) ! provider of pseudo obs + endif if (save_jacobian) then nnz = 2 ! number of non-zero elements in dH(x)/dx profile nind = 1 @@ -1233,12 +1243,12 @@ subroutine setupt(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=nobs endif - if(netcdf_diag) call contents_netcdf_diagp_ + if(netcdf_diag) call contents_netcdf_diagp_(my_diag_pbl) end if prest = prest - pps_press_incr @@ -1260,20 +1270,25 @@ subroutine setupt(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)' t',nchar,nreal,ii+iip,mype,idia0 + write(7)' t',nchar,nreal,ii+iip,mype,idia0,iip if(l_pbl_pseudo_surfobst .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_surfobst .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_surfobst) deallocate(cprvstgp,csprvstgp) endif end if + deallocate(cdiagbuf,rdiagbuf) + if(l_pbl_pseudo_surfobst) deallocate(cdiagbufp,rdiagbufp) end if @@ -1550,7 +1565,7 @@ subroutine contents_binary_diag_(odiag) enddo endif - if (twodvar_regional) then + if (twodvar_regional .or. l_obsprvdiag) then idia = idia + 1 rdiagbuf(idia,ii) = data(idomsfc,i) ! dominate surface type idia = idia + 1 @@ -1568,12 +1583,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) @@ -1604,8 +1621,44 @@ subroutine contents_binary_diagp_ rdiagbufp(20,iip) = 1.e10_r_single ! spread (filled in by EnKF) idia=idia0 +!---- + if (lobsdiagsave) then + do jj=1,miter + idia=idia+1 + if (odiag%muse(jj)) then + rdiagbufp(idia,iip) = one + else + rdiagbufp(idia,iip) = -one + endif + enddo + do jj=1,miter+1 + idia=idia+1 + rdiagbufp(idia,iip) = odiag%nldepart(jj) + enddo + do jj=1,miter + idia=idia+1 + rdiagbufp(idia,iip) = odiag%tldepart(jj) + enddo + do jj=1,miter + idia=idia+1 + rdiagbufp(idia,iip) = odiag%obssen(jj) + enddo + endif + + if (twodvar_regional .or. l_obsprvdiag) then + idia = idia + 1 + rdiagbufp(idia,iip) = -9999. ! data(idomsfc,i) ! dominate surface type + idia = idia + 1 + rdiagbufp(idia,iip) = -9999. ! data(izz,i) ! model terrain at observation 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(idia+1:nreal,ii)) + call writearray(dhx_dx, rdiagbufp(idia+1:nreal,iip)) idia = idia + size(dhx_dx) endif @@ -1684,7 +1737,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) @@ -1704,15 +1757,19 @@ 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 = ' t' real(r_single),parameter:: missing = -9.99e9_r_single + 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)) ) @@ -1738,6 +1795,34 @@ subroutine contents_netcdf_diagp_ call nc_diag_metadata("Forecast_unadjusted", sngl(data(itob,i)-ddiff)) call nc_diag_metadata("Forecast_adjusted", sngl(data(itob,i)-ddiff)) +!---- + 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) diff --git a/src/gsi/setupw.f90 b/src/gsi/setupw.f90 index 57f0111a90..1c6ce62062 100644 --- a/src/gsi/setupw.f90 +++ b/src/gsi/setupw.f90 @@ -41,6 +41,7 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav use obsmod, only: luse_obsdiag use obsmod, only: netcdf_diag, binary_diag, dirname + use obsmod, only: l_obsprvdiag use obsmod, only: neutral_stability_windfact_2dvar,use_similarity_2dvar use nc_diag_write_mod, only: nc_diag_init, nc_diag_header, nc_diag_metadata, & nc_diag_write, nc_diag_data2d @@ -218,6 +219,10 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ! level; they are now loaded by ! aircraftinfo. ! 2020-05-04 wu - no rotate_wind for fv3_regional +! 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). ! ! REMARKS: ! language: f90 @@ -398,7 +403,10 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav ioff0=25 nreal=ioff0 if (lobsdiagsave) nreal=nreal+7*miter+2 - 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 = 2 ! number of non-zero elements in dH(x)/dx profile nind = 1 @@ -1482,13 +1490,13 @@ subroutine setupw(obsLL,odiagLL,lunin,mype,bwork,awork,nele,nobs,is,conv_diagsav if(binary_diag .and. ii>0)then write(7)' uv',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 @@ -1749,7 +1757,7 @@ subroutine contents_binary_diag_(udiag,vdiag) 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 @@ -1853,7 +1861,7 @@ subroutine contents_netcdf_diag_(udiag,vdiag) !++ call nc_diag_data2d("ObsDiagSave_obssen", vdiag%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) diff --git a/util/Analysis_Utilities/read_diag/read_diag_conv.f90 b/util/Analysis_Utilities/read_diag/read_diag_conv.f90 index b7c13f9b2d..42618e2566 100644 --- a/util/Analysis_Utilities/read_diag/read_diag_conv.f90 +++ b/util/Analysis_Utilities/read_diag/read_diag_conv.f90 @@ -48,35 +48,55 @@ PROGRAM read_diag_conv ! read in variables ! character(8),allocatable,dimension(:):: cdiagbuf + character(8),allocatable,dimension(:):: cprvstg,csprvstg ! obs provider/sub-provider + character(8)::cprovider,csubprovider real(r_single),allocatable,dimension(:,:)::rdiagbuf - integer(i_kind) nchar,nreal,ii,mype + integer(i_kind) nchar,nreal,ii,mype,iip ! iip: number of pseudo-obs if existing integer(i_kind) idate ! ! namelist files ! character(180) :: infilename ! file from GSI running directory character(180) :: outfilename ! file name saving results - namelist/iosetup/ infilename, outfilename + logical :: l_obsprvdiag ! if true, read/write obs provider/sub-provider info + logical :: dump_pseudo_obs_too !if true write out pseudo obs as well + namelist/iosetup/ infilename, outfilename, l_obsprvdiag, dump_pseudo_obs_too ! ! output variables ! character(len=3) :: var - real :: rlat,rlon,rprs,robs1,rdpt1,robs2,rdpt2,ruse,rerr - real :: rdhr, ddiff + real(r_single) :: rlat,rlon,rprs,rhgt,robs1,rdpt1,robs2,rdpt2,ruse,rerr + real(r_single) :: rdhr,iusev,ddiff character(8) :: stationID - integer :: itype,iuse,iusev + integer(i_kind) :: itype,iuse + integer(i_kind) :: isubtype,isubtype0 ! ! misc. ! character :: ch - integer :: i,j,k,ios - integer :: ic, iflg + integer(i_kind) :: i,j,k,ios + integer(i_kind) :: ic, iflg + + logical fexist ! - outfilename='diag_results' - open(11,file='namelist.conv') - read(11,iosetup) - close(11) +! initialization of variables in namelist + data infilename /'diag_conv.dat'/ + data outfilename /'diag_results'/ + data l_obsprvdiag /.false./ ! no obs provider info (by default) + data dump_pseudo_obs_too /.false./ +! outfilename='diag_results' + + inquire(file='namelist.conv',exist=fexist) + if(fexist) then + open(11,file='namelist.conv') + read(11,iosetup) + close(11) + else + write(6,*) "no reading from namelist file." + endif + write(6,*) "checking the input/output setup info:" + write(6,iosetup) ! open(42, file=trim(outfilename),IOSTAT=ios) if(ios > 0 ) then @@ -100,25 +120,41 @@ PROGRAM read_diag_conv write(*,*) var, nchar,nreal,ii,mype if (ii > 0) then allocate(cdiagbuf(ii),rdiagbuf(nreal,ii)) + if (l_obsprvdiag) allocate(cprvstg(ii),csprvstg(ii)) read(17,ERR=999,end=110) cdiagbuf, rdiagbuf + if (l_obsprvdiag) then + cprvstg='XXXXXXXX' + csprvstg='XXXXXXXX' + if (var(1:3)==' t' .or. var(1:3)==' q' .or. var(2:3)=='ps' .or. & + var(2:3)=='uv' .or. var(1:3)=='spd') read(17)cprvstg,csprvstg + endif do i=1,ii + if (l_obsprvdiag) then + cprovider=cprvstg(i) + csubprovider=csprvstg(i) + endif itype=rdiagbuf(1,i) ! observation type + isubtype=rdiagbuf(2,i) ! observation subtype rlat=rdiagbuf(3,i) ! observation latitude (degrees) rlon=rdiagbuf(4,i) ! observation longitude (degrees) rprs=rdiagbuf(6,i) ! observation pressure (hPa) + rhgt=rdiagbuf(7,i) ! observation height (meters) rdhr=rdiagbuf(8,i) ! obs time (hours relative to analysis time) iuse=int(rdiagbuf(12,i)) ! analysis usage flag (1=use, -1=monitoring ) iusev=int(rdiagbuf(11,i)) ! analysis usage flag ( value ) ddiff=rdiagbuf(18,i) ! obs-ges used in analysis (K) - rerr = 0 - if (rdiagbuf(16,i) > 0) then ! final inverse observation error (K**-1) - rerr=1.0/rdiagbuf(16,i) - end if + rerr = 0._r_single + if (rdiagbuf(16,i) > 1.0E-12_r_single) then ! final inverse observation error (K**-1) + rerr = 1.0/rdiagbuf(16,i) + else + rerr = 0._r_single + end if robs1=rdiagbuf(17,i) ! observation (K) rdpt1=rdiagbuf(18,i) ! obs-ges used in analysis ! get station ID stationID = cdiagbuf(i) +! Remove odd spaces in the station ID iflg = 0 do ic=8,1,-1 ch = stationID(ic:ic) @@ -131,12 +167,42 @@ PROGRAM read_diag_conv stationID(ic:ic) = '_' endif enddo + +! Remove odd spaces in the obs provider, and subprovider names + if (l_obsprvdiag) then + iflg = 0 + do ic=8,1,-1 + ch = cprovider(ic:ic) + if (ch > ' ' .and. ch <= 'z') then + iflg = 1 + else + cprovider(ic:ic) = ' ' + end if + if (ch == ' ' .and. iflg == 1) then + cprovider(ic:ic) = '_' + endif + enddo + + iflg = 0 + do ic=8,1,-1 + ch = csubprovider(ic:ic) + if (ch > ' ' .and. ch <= 'z') then + iflg = 1 + else + csubprovider(ic:ic) = ' ' + end if + if (ch == ' ' .and. iflg == 1) then + csubprovider(ic:ic) = '_' + endif + enddo + endif ! ! When the data is q, unit convert kg/kg -> g/kg **/ if (var == " q") then robs1 = robs1 * 1000.0 rdpt1 = rdpt1 * 1000.0 rerr = rerr * 1000.0 + ddiff = ddiff * 1000.0 end if ! When the data is pw, replase the rprs to -999.0 **/ if (var == " pw") rprs=-999.0 @@ -145,24 +211,74 @@ PROGRAM read_diag_conv robs1=-99999.9 ddiff=-99999.9 endif + +! check up the information in the obs provider, and subprovider names + if (l_obsprvdiag) then + if (cprovider(1:4)=='B7Hv' .or. cprovider(5:8)=='vH7B') then !this provider name comes with strange characters + cprovider(1:4)='B7Hv' + cprovider(5:8)=' ' + endif + + if (csubprovider(1:4)=='B7Hv' .or. csubprovider(5:8)=='vH7B') then + csubprovider(1:4)='B7Hv' + csubprovider(5:8)=' ' + endif + + if (itype==154 .and. trim(adjustl(var))=='tca') then + stationID='GOESSKY' + cprovider='GOESSKY' + csubprovider='GOESSKY' + end if + + if (trim(cprovider)=='') cprovider='EMPTY' + if (trim(csubprovider)=='') csubprovider='EMPTY' + endif + +! special treatment to obs tca/cei/vis/gst +! If we have ceiling or total cloud amount obs less than 0, +! set to missing + if ((trim(adjustl(var))=="tca" .or. trim(adjustl(var))=="cei" .or. & + trim(adjustl(var))=="vis" .or. trim(adjustl(var))=="gst") .and. & + robs1 < 0.) then + robs1=0.10000E+10 + ddiff=0.10000E+10 + end if + + if (dump_pseudo_obs_too) then + isubtype0=0 ! reset subtype (esp.for pseudo-obs) to be zero in order to print it our + else + isubtype0=isubtype + endif ! ! write out result for one variable on one pitch - if (var .ne. " uv") then - write (42,'(A3," @ ",A8," : ",I3,F10.2,F8.2,F8.2,F8.2,I5,2F10.2)') & + if (l_obsprvdiag) then ! write out obs provider info with other info together + if ( var .ne. " uv" .and. isubtype0 >=0 ) then + write (42,'(A3,1x,A8,1x,A8,1x,A8,1x,I3,1x,F10.2,F8.2,F8.2,2F20.5,I5,2E15.5,1x,"NaN NaN ",E15.5,F10.3)') & + var,stationID,cprovider,csubprovider,itype,rdhr,rlat,rlon,rprs,rhgt,iuse,robs1,ddiff,rerr,iusev + else if ( var .eq. " uv" .and. isubtype0 >=0 ) then +! ** When the data is uv, additional output is needed **/ + robs2=rdiagbuf(20,i) + rdpt2=rdiagbuf(21,i) + write (42,'(A3,1x,A8,1x,A8,1x,A8,1x,I3,1x,F10.2,F8.2,F8.2,2F20.5,I5,4E15.5,1x,E15.5,F10.3)') & + var,stationID,cprovider,csubprovider,itype,rdhr,rlat,rlon,rprs,rhgt,iuse,robs1,ddiff,robs2,rdpt2,rerr,iusev + endif + else ! if no need to write out obs provider info + if (var .ne. " uv") then + write (42,'(A3," @ ",A8," : ",I3,F10.2,F8.2,F8.2,F8.2,I5,2F10.2)') & var,stationID,itype,rdhr,rlat,rlon,rprs,iuse,robs1,ddiff - else + else ! ** When the data is uv, additional output is needed **/ - robs2=rdiagbuf(20,i) - rdpt2=rdiagbuf(21,i) - write (42,'(A3," @ ",A8," : ",I3,F10.2,F8.2,F8.2,F8.2,I5,4F10.2)') & + robs2=rdiagbuf(20,i) + rdpt2=rdiagbuf(21,i) + write (42,'(A3," @ ",A8," : ",I3,F10.2,F8.2,F8.2,F8.2,I5,4F10.2)') & var,stationID,itype,rdhr,rlat,rlon,rprs,iuse,robs1,ddiff,robs2, rdpt2 + endif endif - - enddo ! i end for one station deallocate(cdiagbuf,rdiagbuf) + if (l_obsprvdiag) deallocate(cprvstg,csprvstg) else read(17) endif