Merge branch 'NOAA-EMC:master' into feature/rm_sfcanl

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,ompslp_mult_fact
+  use obsmod, only: l_obsprvdiag
   use obsmod, only: aircraft_recon, &
 
        ! The following variables are the coefficients that describe
@@ -480,6 +481,9 @@ module gsimod
 !                       option for checking and adjusting the profile of Qr/Qs/Qg/Qnr
 !                       retrieved through cloud analysis to reduce the background
 !                       reflectivity ghost in analysis. (default is 0)
+!  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)
 !  01-07-2022 Hu        Add fv3_io_layout_y to let fv3lam interface read/write subdomain restart
 !                       files. The fv3_io_layout_y needs to match fv3lam model
 !                       option io_layout(2).
@@ -699,6 +703,8 @@ module gsimod
 !                     (.TRUE.: on; .FALSE.: off) / Inputfile: l2rwbufr_cltl (bufr format)
 !     l_use_dbz_directDA - option to assimilate radar reflectivity obs directly in GSI 
 !                     (.TRUE.: on; .FALSE.: off) / Inputfile: dbzbufr (bufr format)
+!     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.
@@ -743,7 +749,7 @@ module gsimod
        if_model_dbz,imp_physics,lupp,netcdf_diag,binary_diag,l_wcp_cwm,aircraft_recon,diag_version,&
        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,&
-       l_reg_update_hydro_delz, &
+       l_reg_update_hydro_delz, l_obsprvdiag,&
        l_use_dbz_directDA, l_use_rw_columntilt
 
 ! GRIDOPTS (grid setup variables,including regional specific variables):

src/gsi/obsmod.F90

@@ -157,6 +157,9 @@ module obsmod
 !                          GSI namelist level.  
 !  2020-09-15  Wu        - add option tcp_posmatch to mitigate possibility of erroneous TC initialization
 !  2020-09-19  CAPS(J. Park) - add 'vad_near_analtime' flag to assimilate newvad obs around analysis time only
+!   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
@@ -398,6 +401,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
@@ -484,6 +488,7 @@ module obsmod
   public :: nobs_sub
 
   public :: netcdf_diag, binary_diag
+  public :: l_obsprvdiag
 
   public :: l_wcp_cwm
   public :: aircraft_recon
@@ -553,6 +558,7 @@ module obsmod
 
   logical luse_obsdiag
   logical binary_diag, netcdf_diag 
+  logical l_obsprvdiag 
 
 ! Declare types
 
@@ -719,6 +725,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:
 !
@@ -911,6 +918,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 

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
@@ -904,6 +913,8 @@ subroutine contents_netcdf_diag_(odiag)
            call nc_diag_metadata("Observation",                   sngl(pob)        )
            call nc_diag_metadata("Obs_Minus_Forecast_adjusted",   sngl(pob-pges)   )
            call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(pob-pgesorig))
+           call nc_diag_metadata("Forecast_adjusted",             sngl(pges)       )
+           call nc_diag_metadata("Forecast_unadjusted",           sngl(pgesorig)   )
 
           if (lobsdiagsave) then
 
@@ -921,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)
@@ -936,6 +947,8 @@ subroutine contents_netcdf_diag_(odiag)
              call nc_diag_data2d("Observation_Operator_Jacobian_val", real(dhx_dx%val,r_single))
            endif
 
+           call nc_diag_data2d("atmosphere_pressure_coordinate", sngl(exp(prsltmp)*r1000))
+
   end subroutine contents_netcdf_diag_
 
   subroutine final_vars_