diff --git a/fix b/fix
index 6720dce3c..384204259 160000
--- a/fix
+++ b/fix
@@ -1 +1 @@
-Subproject commit 6720dce3c168ab3a211c265d2a870c20bd220e35
+Subproject commit 384204259f72b55ae8dd1b3e6959b5bc735ce904
diff --git a/src/gsi/gsimod.F90 b/src/gsi/gsimod.F90
index 95ff7f79b..7177719fd 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,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):
diff --git a/src/gsi/obsmod.F90 b/src/gsi/obsmod.F90
index 2d9209e77..40dc3cc0d 100644
--- a/src/gsi/obsmod.F90
+++ b/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 
diff --git a/src/gsi/setupps.f90 b/src/gsi/setupps.f90
index 69379f15f..04344a2f7 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
@@ -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_
diff --git a/src/gsi/setupq.f90 b/src/gsi/setupq.f90
index 8189b0bcb..719b14f75 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)
@@ -920,11 +931,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
@@ -945,20 +956,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
@@ -1148,7 +1164,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
@@ -1165,12 +1181,13 @@ 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)  = -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)
@@ -1201,8 +1218,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._r_single !  data(idomsfc,i) ! dominate surface type
+           ioff = ioff + 1
+           rdiagbufp(ioff,iip) = -9999._r_single !  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
 
@@ -1240,6 +1292,8 @@ subroutine contents_netcdf_diag_(odiag)
            call nc_diag_metadata("Observation",                   sngl(data(iqob,i)))
            call nc_diag_metadata("Obs_Minus_Forecast_adjusted",   sngl(ddiff)       )
            call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(qob-qges)    )
+           call nc_diag_metadata("Forecast_adjusted",             sngl(data(iqob,i)-ddiff))
+           call nc_diag_metadata("Forecast_unadjusted",           sngl(qges))
            call nc_diag_metadata("Forecast_Saturation_Spec_Hum",  sngl(qsges)       )
            if (lobsdiagsave) then
               do jj=1,miter
@@ -1256,7 +1310,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)
@@ -1270,16 +1324,20 @@ 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",exp(prsltmp)*r1000)
+
   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",     -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))  )
@@ -1302,13 +1360,42 @@ subroutine contents_netcdf_diagp_
            call nc_diag_metadata("Observation",                   sngl(data(iqob,i)))
            call nc_diag_metadata("Obs_Minus_Forecast_adjusted",   sngl(ddiff)       )
            call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(ddiff)       )
+           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",     "88888888"                   )    
+              r_sprvstg           = data(isprvd,i)
+              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)
              call nc_diag_data2d("Observation_Operator_Jacobian_val", real(dhx_dx%val,r_single))
            endif
+
+           call nc_diag_data2d("atmosphere_pressure_coordinate",exp(prsltmp)*r1000)
+
   end subroutine contents_netcdf_diagp_
 
   subroutine final_vars_
diff --git a/src/gsi/setupspd.f90 b/src/gsi/setupspd.f90
index 1e740fdd2..91b2467bf 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 d94d1f26d..e7e73d82c 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,13 @@ 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)  = -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 +1620,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._r_single ! data(idomsfc,i) ! dominate surface type
+       idia = idia + 1
+       rdiagbufp(idia,iip) = -9999._r_single ! 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
 
@@ -1645,6 +1697,9 @@ subroutine contents_netcdf_diag_(odiag)
     call nc_diag_metadata("Observation",             sngl(data(itob,i))     )
     call nc_diag_metadata("Obs_Minus_Forecast_adjusted",   sngl(ddiff)      )
     call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(tob-tges)   )
+    call nc_diag_metadata("Forecast_unadjusted",     sngl(tges))
+    call nc_diag_metadata("Forecast_adjusted",       sngl(data(itob,i)-ddiff))
+
     if (aircraft_t_bc_pof .or. aircraft_t_bc .or. aircraft_t_bc_ext) then
        call nc_diag_metadata("Data_Pof",             sngl(data(ipof,i))     )
        call nc_diag_metadata("Data_Vertical_Velocity", sngl(data(ivvlc,i))  )
@@ -1681,7 +1736,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)
@@ -1696,17 +1751,23 @@ subroutine contents_netcdf_diag_(odiag)
        call nc_diag_data2d("Observation_Operator_Jacobian_val", real(dhx_dx%val,r_single))
     endif
 
+    ! need additional arrays for GeoVaLs for T2
+    call nc_diag_data2d("atmosphere_pressure_coordinate", sngl(exp(prsltmp)*r1000))
+
   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",     -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))  )
@@ -1729,13 +1790,42 @@ subroutine contents_netcdf_diagp_
     call nc_diag_metadata("Observation",             sngl(data(itob,i))     )
     call nc_diag_metadata("Obs_Minus_Forecast_adjusted",   sngl(ddiff)      )
     call nc_diag_metadata("Obs_Minus_Forecast_unadjusted", sngl(ddiff)      )
+    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)                  )
+       call nc_diag_metadata("Provider_Name",     "88888888"                   )
+       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)
        call nc_diag_data2d("Observation_Operator_Jacobian_val", real(dhx_dx%val,r_single))
     endif
 
+    ! need additional arrays for GeoVaLs for T2
+    call nc_diag_data2d("atmosphere_pressure_coordinate", sngl(exp(prsltmp)*r1000))
+
   end subroutine contents_netcdf_diagp_
 
   subroutine final_vars_
diff --git a/src/gsi/setupw.f90 b/src/gsi/setupw.f90
index 33f0a6295..55e498fe6 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
@@ -1485,13 +1493,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
 
 
@@ -1752,7 +1760,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
@@ -1826,10 +1834,14 @@ subroutine contents_netcdf_diag_(udiag,vdiag)
               call nc_diag_metadata("u_Observation",                              sngl(uob_e)           )
               call nc_diag_metadata("u_Obs_Minus_Forecast_adjusted",              sngl(dudiff_e)        )
               call nc_diag_metadata("u_Obs_Minus_Forecast_unadjusted",            sngl(uob_e-uges_e)    )
+              call nc_diag_metadata("u_Forecast_adjusted",                        sngl(uob_e-dudiff_e))
+              call nc_diag_metadata("u_Forecast_unadjusted",                      sngl(uges_e))
 
               call nc_diag_metadata("v_Observation",                              sngl(vob_e)           )
               call nc_diag_metadata("v_Obs_Minus_Forecast_adjusted",              sngl(dvdiff_e)        )
               call nc_diag_metadata("v_Obs_Minus_Forecast_unadjusted",            sngl(vob_e-vges_e)    )
+              call nc_diag_metadata("v_Forecast_adjusted",                        sngl(vob_e-dvdiff_e))
+              call nc_diag_metadata("v_Forecast_unadjusted",                      sngl(vges_e))
            endif
 
            if (lobsdiagsave) then
@@ -1852,7 +1864,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)
@@ -1869,6 +1881,7 @@ subroutine contents_netcdf_diag_(udiag,vdiag)
               call nc_diag_data2d("v_Observation_Operator_Jacobian_val", real(dhx_dx_v%val,r_single))
            endif
 
+           call nc_diag_data2d("atmosphere_pressure_coordinate", exp(prsltmp)*r1000)
 
   end subroutine contents_netcdf_diag_
 
diff --git a/util/Analysis_Utilities/read_diag/read_diag_conv.f90 b/util/Analysis_Utilities/read_diag/read_diag_conv.f90
index b7c13f9b2..42618e256 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
diff --git a/util/Radiance_Monitor/data_extract/ush/nu_find_cycle.pl b/util/Radiance_Monitor/data_extract/ush/nu_find_cycle.pl
index 10a505adb..f304279a0 100755
--- a/util/Radiance_Monitor/data_extract/ush/nu_find_cycle.pl
+++ b/util/Radiance_Monitor/data_extract/ush/nu_find_cycle.pl
@@ -88,7 +88,7 @@
       $search_string = $run;
    }
 
-   my @mmdirs = grep { /$search_string/ } @alldirs;
+   my @mmdirs = grep {/$search_string/} @alldirs;
    
    #-----------------------------------------------------------------------   
    #  If there are no $run.yyyymmdd subdirectories, then exit without 
diff --git a/util/Radiance_Monitor/image_gen/html/Install_html.sh b/util/Radiance_Monitor/image_gen/html/Install_html.sh
index ec81991e6..8ecb71d1d 100755
--- a/util/Radiance_Monitor/image_gen/html/Install_html.sh
+++ b/util/Radiance_Monitor/image_gen/html/Install_html.sh
@@ -51,12 +51,6 @@ else
 fi
 
 
-#--------------------------------------------------------------
-#  source plot_rad_conf to get WEB_SVR, WEB_USER, WEBDIR
-#
-. ${RADMON_IMAGE_GEN}/parm/plot_rad_conf
-
-
 #--------------------------------------------------------------
 #  call the appropriate child script for glb or rgn
 #
diff --git a/util/Radiance_Monitor/image_gen/html/install_glb.sh b/util/Radiance_Monitor/image_gen/html/install_glb.sh
index eb55e7c20..4c70951df 100755
--- a/util/Radiance_Monitor/image_gen/html/install_glb.sh
+++ b/util/Radiance_Monitor/image_gen/html/install_glb.sh
@@ -33,35 +33,8 @@ RAD_AREA="glb"
 this_file=`basename $0`
 this_dir=`dirname $0`
 
-#top_parm=${this_dir}/../../parm
-#
-#if [[ -s ${top_parm}/RadMon_config ]]; then
-#   . ${top_parm}/RadMon_config
-#else
-#   echo "ERROR:  Unable to source ${top_parm}/RadMon_config"
-#   exit
-#fi
-#
-#if [[ -s ${top_parm}/RadMon_user_settings ]]; then
-#   . ${top_parm}/RadMon_user_settings
-#else
-#   echo "ERROR:  Unable to source ${top_parm}/RadMon_user_settings"
-#   exit
-#fi
-
-
-#--------------------------------------------------------------
-#  source plot_rad_conf to get WEB_SVR, WEB_USER, WEBDIR
-#
-. ${RADMON_IMAGE_GEN}/parm/plot_rad_conf
-
-
-#--------------------------------------------------------------
-#  Get the area for this SUFFIX from the data_map file
-#
 
 new_webdir=${WEBDIR}/${SUFFIX}
-. ${RADMON_IMAGE_GEN}/parm/glbl_conf
 
 echo RAD_AREA = $RAD_AREA
 echo TANKverf = $TANKverf
@@ -85,7 +58,7 @@ cd $workdir
 #  Find the first date with data.  Start at today and work
 #  backwards.  Stop after 90 days and exit.
 #
-PDATE=`${IG_SCRIPTS}/find_cycle.pl --dir ${TANKverf} --cyc 1`
+PDATE=`${IG_SCRIPTS}/nu_find_cycle.pl --dir ${TANKverf} --cyc 1`
 echo PDATE= $PDATE
 
 limit=`$NDATE -2160 $PDATE`		# 90 days
@@ -98,35 +71,30 @@ echo limit, PDATE = $limit, $PDATE
 data_found=0
 while [[ data_found -eq 0 && $PDATE -ge $limit ]]; do
    PDY=`echo $PDATE|cut -c1-8`
+   CYC=`echo $PDATE|cut -c9-10`
 
-   test_dir=${TANKverf}/${RUN}.${PDY}/${MONITOR}
+   test_dir=${TANKverf}/${RUN}.${PDY}/${CYC}/${MONITOR}
    if [[ ! -d ${test_dir} ]]; then
-      test_dir=${TANKverf}/${RUN}.${PDY}
+      test_dir=${TANKverf}/${RUN}.${PDY}/${MONITOR}
    fi
    if [[ ! -d ${test_dir} ]]; then
-      test_dir=${TANKverf}/${MONITOR}.${PDY}
+      test_dir=${TANKverf}/${RUN}.${PDY}
    fi
    
    echo "test_dir = ${test_dir}"
 
    if [[ -d ${test_dir} ]]; then
       echo " test_dir is GO "
-      test00=`ls ${test_dir}/angle.*${PDY}00*.ieee_d* | wc -l`
-      test06=`ls ${test_dir}/angle.*${PDY}06*.ieee_d* | wc -l`
-      test12=`ls ${test_dir}/angle.*${PDY}12*.ieee_d* | wc -l`
-      test18=`ls ${test_dir}/angle.*${PDY}18*.ieee_d* | wc -l`
-      if [[ $test00 -gt 0 ]]; then
-         test_list=`ls ${test_dir}/angle.*${PDY}00*.ieee_d*`
-         data_found=1
-      elif [[ $test06 -gt 0 ]]; then
-         test_list=`ls ${test_dir}/angle.*${PDY}06*.ieee_d*`
-         data_found=1
-      elif [[ $test12 -gt 0 ]]; then
-         test_list=`ls ${test_dir}/angle.*${PDY}12*.ieee_d*`
-         data_found=1
-      elif [[ $test18 -gt 0 ]]; then
-         test_list=`ls ${test_dir}/angle.*${PDY}18*.ieee_d*`
-         data_found=1
+
+      if [[ -e ${test_dir}/angle.tar ]]; then
+         test_list=`tar -tv ${test_dir}/angle.tar`
+	 data_found=1 
+      else
+         test=`ls ${test_dir}/angle.*${PDATE}*.ieee_d* | wc -l`
+         if [[ $test -gt 0 ]]; then
+            test_list=`ls ${test_dir}/angle.*${PDATE}*.ieee_d*`
+            data_found=1
+	 fi
       fi
    else
       echo "test_dir is NOGO"
diff --git a/util/Radiance_Monitor/image_gen/html/install_rgn.sh b/util/Radiance_Monitor/image_gen/html/install_rgn.sh
index ce2a4de31..6ceb1a4ca 100755
--- a/util/Radiance_Monitor/image_gen/html/install_rgn.sh
+++ b/util/Radiance_Monitor/image_gen/html/install_rgn.sh
@@ -50,18 +50,10 @@ else
 fi
 
 
-#--------------------------------------------------------------
-#  source plot_rad_conf to get WEB_SVR, WEB_USER, WEBDIR
-#
-. ${RADMON_IMAGE_GEN}/parm/plot_rad_conf
-
-
 #--------------------------------------------------------------
 #  Get the area for this SUFFIX from the data_map file
 #
-
 new_webdir=${WEBDIR}/${SUFFIX}
-. ${RADMON_IMAGE_GEN}/parm/rgnl_conf
 
 echo RAD_AREA = $RAD_AREA
 echo TANKverf = $TANKverf
diff --git a/util/Radiance_Monitor/image_gen/ush/RadMon_IG_glb.sh b/util/Radiance_Monitor/image_gen/ush/RadMon_IG_glb.sh
index 6ab9b4f70..644653df3 100755
--- a/util/Radiance_Monitor/image_gen/ush/RadMon_IG_glb.sh
+++ b/util/Radiance_Monitor/image_gen/ush/RadMon_IG_glb.sh
@@ -351,8 +351,9 @@ if [[ $RUN_TRANSFER -eq 1 ]]; then
          echo "${IG_SCRIPTS}/Transfer.sh --nosrc ${RADMON_SUFFIX}" >$cmdfile
          chmod 755 $cmdfile
 
+         run_time="$rhr$cmin"	# HHMM format for qsub
          $SUB -q $transfer_queue -A $ACCOUNT -o ${transfer_log} -e ${LOGdir}/Transfer_${RADMON_SUFFIX}.err \
-	      -V -l select=1:mem=500M -l walltime=45:00 -N ${jobname} ${cmdfile}
+	      -V -l select=1:mem=500M -l walltime=45:00 -N ${jobname} -a ${run_time} ${cmdfile}
       else
          $SUB -P $PROJECT -q $transfer_queue -o ${transfer_log} -M 80 -W 0:45 \
               -R affinity[core] -J ${jobname} -cwd ${PWD} -b $run_time ${job}
diff --git a/util/Radiance_Monitor/image_gen/ush/nu_find_cycle.pl b/util/Radiance_Monitor/image_gen/ush/nu_find_cycle.pl
index 10a505adb..f304279a0 100755
--- a/util/Radiance_Monitor/image_gen/ush/nu_find_cycle.pl
+++ b/util/Radiance_Monitor/image_gen/ush/nu_find_cycle.pl
@@ -88,7 +88,7 @@
       $search_string = $run;
    }
 
-   my @mmdirs = grep { /$search_string/ } @alldirs;
+   my @mmdirs = grep {/$search_string/} @alldirs;
    
    #-----------------------------------------------------------------------   
    #  If there are no $run.yyyymmdd subdirectories, then exit without