read_diag.f90

Go to the documentation of this file.

!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_raddiag                       read rad diag file
!   prgmmr: tahara           org: np20                date: 2003-01-01
!
! abstract:  This module contains code to process radiance
!            diagnostic files.  The module defines structures
!            to contain information from the radiance
!            diagnostic files and then provides two routines
!            to access contents of the file.
!
! program history log:
!   2005-07-22 treadon - add this doc block
!   2010-10-05 treadon - refactor code to GSI standard
!   2010-10-08 zhu     - use data_tmp to handle various npred values
!   2011-02-22 kleist  - changes related to memory allocate/deallocate
!   2011-04-08 li      - add tref, dtw, dtc to diag_data_fix_list, add tb_tz to diag_data_chan_list
!                      - correspondingly, change ireal_radiag (26 -> 30) and ipchan_radiag (7 -> 8)
!   2011-07-24 safford - make structure size for reading data_fix data version dependent 
!   2013-11-21 todling - revisit how versions are set (add set/get_radiag)
!   2014-01-27 todling - add ob sensitivity index
!   2017-07-13 mccarty - incorporate hooks for nc4/binary diag reading
!
! contains
!   read_radiag_header - read radiance diagnostic file header
!   read_radiag_data   - read radiance diagnostic file data
!   set_netcdf_read    - call set_netcdf_read(.true.) to use nc4 hooks, otherwise read file as 
!                           traditional binary format
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

module read_diag

  use ncd_kinds, only:  i_kind,r_single,r_kind
  use nc_diag_read_mod, only: nc_diag_read_get_var,  nc_diag_read_get_global_attr
  use nc_diag_read_mod, only: nc_diag_read_get_dim
  use nc_diag_read_mod, only: nc_diag_read_init, nc_diag_read_close
  implicit none

! Declare public and private
  private

  public :: diag_header_fix_list
  public :: diag_header_chan_list
  public :: diag_data_name_list
  public :: diag_data_fix_list
  public :: diag_data_chan_list
  public :: diag_data_extra_list
  public :: open_radiag
  public :: close_radiag
  public :: read_radiag_header
  public :: read_radiag_data
  public :: set_netcdf_read
! public :: iversion_radiag
! public :: iversion_radiag_1
! public :: iversion_radiag_2
! public :: iversion_radiag_3
! public :: iversion_radiag_4
  public :: ireal_radiag
  public :: ipchan_radiag
  public :: set_radiag
  public :: get_radiag
  public :: read_all_radiag

  interface set_radiag
         module procedure set_radiag_int_ ! internal procedure for integers
  end interface
  interface get_radiag
         module procedure get_radiag_int_ ! internal procedure for integers
  end interface

  integer(i_kind),parameter :: ireal_radiag  = 30   ! number of real entries per spot in radiance diagnostic file
  integer(i_kind),parameter :: ireal_old_radiag  = 26   ! number of real entries per spot in versions older than iversion_radiag_2
  integer(i_kind),parameter :: ipchan_radiag = 8    ! number of entries per channel per spot in radiance diagnostic file

! Declare structures for radiance diagnostic file information
  type diag_header_fix_list
     character(len=20) :: isis           ! sat and sensor type
     character(len=10) :: id             ! sat type
     character(len=10) :: obstype        ! observation type
     integer(i_kind) :: jiter            ! outer loop counter
     integer(i_kind) :: nchan            ! number of channels in the sensor
     integer(i_kind) :: npred            ! number of updating bias correction predictors
     integer(i_kind) :: idate            ! time (yyyymmddhh)
     integer(i_kind) :: ireal            ! # of real elements in the fix part of a data record
     integer(i_kind) :: ipchan           ! # of elements for each channel except for bias correction terms
     integer(i_kind) :: iextra           ! # of extra elements for each channel
     integer(i_kind) :: jextra           ! # of extra elements
     integer(i_kind) :: idiag            ! first dimension of diag_data_chan
     integer(i_kind) :: angord           ! order of polynomial for adp_anglebc option
     integer(i_kind) :: iversion         ! radiance diagnostic file version number
     integer(i_kind) :: inewpc           ! indicator of newpc4pred (1 on, 0 off)
     integer(i_kind) :: isens            ! sensitivity index
  end type diag_header_fix_list

  type diag_data_name_list
     character(len=10),dimension(ireal_radiag) :: fix
     character(len=10),dimension(:),allocatable :: chn
  end type diag_data_name_list
  
  type diag_header_chan_list
     real(r_single) :: freq              ! frequency (Hz)
     real(r_single) :: polar             ! polarization
     real(r_single) :: wave              ! wave number (cm^-1)
     real(r_single) :: varch             ! error variance (or SD error?)
     real(r_single) :: tlapmean          ! mean lapse rate
     integer(i_kind):: iuse              ! use flag
     integer(i_kind):: nuchan            ! sensor relative channel number
     integer(i_kind):: iochan            ! satinfo relative channel number
  end type diag_header_chan_list

  type diag_data_fix_list
     real(r_single) :: lat               ! latitude (deg)
     real(r_single) :: lon               ! longitude (deg)
     real(r_single) :: zsges             ! guess elevation at obs location (m)
     real(r_single) :: obstime           ! observation time relative to analysis
     real(r_single) :: senscn_pos        ! sensor scan position (integer(i_kind))
     real(r_single) :: senscn_ang        ! sensor scan angle
     real(r_single) :: satzen_ang        ! satellite zenith angle (deg)
     real(r_single) :: satazm_ang        ! satellite azimuth angle (deg)
     real(r_single) :: solzen_ang        ! solar zenith angle (deg)
     real(r_single) :: solazm_ang        ! solar azimumth angle (deg)
     real(r_single) :: sungln_ang        ! sun glint angle (deg)
     real(r_single) :: water_frac        ! fractional coverage by water
     real(r_single) :: land_frac         ! fractional coverage by land
     real(r_single) :: ice_frac          ! fractional coverage by ice
     real(r_single) :: snow_frac         ! fractional coverage by snow
     real(r_single) :: water_temp        ! surface temperature over water (K)
     real(r_single) :: land_temp         ! surface temperature over land (K)
     real(r_single) :: ice_temp          ! surface temperature over ice (K)
     real(r_single) :: snow_temp         ! surface temperature over snow (K)
     real(r_single) :: soil_temp         ! soil temperature (K)
     real(r_single) :: soil_mois         ! soil moisture 
     real(r_single) :: land_type         ! land type (integer(i_kind))
     real(r_single) :: veg_frac          ! vegetation fraction
     real(r_single) :: snow_depth        ! snow depth
     real(r_single) :: sfc_wndspd        ! surface wind speed
     real(r_single) :: qcdiag1           ! ir=cloud fraction, mw=cloud liquid water
     real(r_single) :: qcdiag2           ! ir=cloud top pressure, mw=total column water
     real(r_single) :: tref              ! reference temperature (Tr) in NSST
     real(r_single) :: dtw               ! dt_warm at zob
     real(r_single) :: dtc               ! dt_cool at zob
     real(r_single) :: tz_tr             ! d(Tz)/d(Tr)
  end type diag_data_fix_list

  type diag_data_chan_list
     real(r_single) :: tbobs              ! Tb (obs) (K)
     real(r_single) :: omgbc              ! Tb_(obs) - Tb_(simulated w/ bc)  (K)
     real(r_single) :: omgnbc             ! Tb_(obs) - Tb_(simulated_w/o bc) (K)
     real(r_single) :: errinv             ! inverse error (K**(-1))
     real(r_single) :: qcmark             ! quality control mark
     real(r_single) :: emiss              ! surface emissivity
     real(r_single) :: tlap               ! temperature lapse rate
     real(r_single) :: tb_tz              ! d(Tb)/d(Tz)
     real(r_single) :: bicons             ! constant bias correction term
     real(r_single) :: biang              ! scan angle bias correction term
     real(r_single) :: biclw              ! CLW bias correction term
     real(r_single) :: bilap2             ! square lapse rate bias correction term
     real(r_single) :: bilap              ! lapse rate bias correction term
     real(r_single) :: bicos              ! node*cos(lat) bias correction term
     real(r_single) :: bisin              ! sin(lat) bias correction term
     real(r_single) :: biemis             ! emissivity sensitivity bias correction term
     real(r_single),dimension(:),allocatable :: bifix          ! angle dependent bias
     real(r_single) :: bisst              ! SST bias correction term
  end type diag_data_chan_list

  type diag_data_extra_list
     real(r_single) :: extra              ! extra information
  end type diag_data_extra_list

  integer(i_kind),save     :: iversion_radiag             ! Current version (see set routine)
  integer(i_kind),parameter:: iversion_radiag_1 = 11104   ! Version when bias-correction entries were modified 
  integer(i_kind),parameter:: iversion_radiag_2 = 13784   ! Version when NSST entries were added 
  integer(i_kind),parameter:: iversion_radiag_3 = 19180   ! Version when SSMIS added
  integer(i_kind),parameter:: iversion_radiag_4 = 30303   ! Version when emissivity predictor added

  real(r_single),parameter::  rmiss_radiag    = -9.9e11_r_single

  logical,save            ::  netcdf = .false. 

  type ncdiag_status
     logical :: nc_read
     integer(i_kind) :: cur_ob_idx
     integer(i_kind) :: num_records
     type(diag_data_fix_list),  allocatable  :: all_data_fix(:)
     type(diag_data_chan_list), allocatable  :: all_data_chan(:,:)
     type(diag_data_extra_list), allocatable :: all_data_extra(:,:,:)
  end type ncdiag_status

  integer(i_kind), parameter :: max_open_ncdiag = 2
  integer(i_kind), save :: nopen_ncdiag = 0
  integer(i_kind), dimension(MAX_OPEN_NCDIAG), save     :: ncdiag_open_id = (/-1, -1/)
  type(ncdiag_status), dimension(MAX_OPEN_NCDIAG), save :: ncdiag_open_status
      
contains

subroutine set_radiag_int_ (what,iv,ier)
character(len=*),intent(in) :: what
integer(i_kind),intent(in) :: iv
integer(i_kind),intent(out):: ier
ier=-1
if(trim(what)=='version') then
  iversion_radiag = iv
  ier=0
endif
end subroutine set_radiag_int_

subroutine get_radiag_int_ (what,iv,ier)
character(len=*),intent(in) :: what
integer(i_kind),intent(out):: iv
integer(i_kind),intent(out):: ier
ier=-1
if(trim(what)=='version') then
  iv = iversion_radiag
  ier=0
endif
end subroutine get_radiag_int_

subroutine set_netcdf_read(use_netcdf)
!                .      .    .                                       .
! subprogram:    read_diag_header_bin             read rad diag header
!   prgmmr: mccarty           org: gmao                date: 2015-08-06
!
! abstract:  This routine sets the routines to read from a netcdf file.
!            The default currently is to read binary files
!
! program history log:
!   2015-08-06 mccarty - created routine 
!
! input argument list:
!   use_netcdf    - logical .true. tells routine to read netcdf diag 
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  logical,intent(in)                     :: use_netcdf

  netcdf = use_netcdf
end subroutine set_netcdf_read


subroutine open_radiag(filename, ftin)
   character*500,   intent(in) :: filename
   integer(i_kind), intent(inout) :: ftin

   integer(i_kind) :: i

   if (netcdf) then
      if (nopen_ncdiag >= max_open_ncdiag) then
         write(6,*) 'OPEN_RADIAG:  ***ERROR*** Cannot open more than ', &
                    max_open_ncdiag, ' netcdf diag files.'
         call abort
      endif
      call nc_diag_read_init(filename,ftin)
      do i = 1, max_open_ncdiag
        if (ncdiag_open_id(i) < 0) then
           ncdiag_open_id(i) = ftin
           ncdiag_open_status(i)%nc_read = .false.
           ncdiag_open_status(i)%cur_ob_idx = -9999
           ncdiag_open_status(i)%num_records = -9999
           if (allocated(ncdiag_open_status(i)%all_data_fix)) then
              deallocate(ncdiag_open_status(i)%all_data_fix)
           endif
           if (allocated(ncdiag_open_status(i)%all_data_chan)) then
              deallocate(ncdiag_open_status(i)%all_data_chan)
           endif
           if (allocated(ncdiag_open_status(i)%all_data_extra)) then
              deallocate(ncdiag_open_status(i)%all_data_extra)
           endif
           nopen_ncdiag = nopen_ncdiag + 1
           exit
        endif
      enddo
   else
      open(ftin,form="unformatted",file=filename)
      rewind(ftin)
   endif

end subroutine open_radiag

subroutine close_radiag(filename, ftin)
   character*500,   intent(in) :: filename
   integer(i_kind), intent(inout) :: ftin

   integer(i_kind) :: id

   if (netcdf) then
      id = find_ncdiag_id(ftin)
      if (id < 0) then
         write(6,*) 'CLOSE_RADIAG:  ***ERROR*** ncdiag file ', filename,   &
                    ' was not opened'
         call abort
      endif
      call nc_diag_read_close(filename)
      ncdiag_open_id(id) = -1
      ncdiag_open_status(id)%nc_read = .false.
      ncdiag_open_status(id)%cur_ob_idx = -9999
      ncdiag_open_status(id)%num_records = -9999
      if (allocated(ncdiag_open_status(id)%all_data_fix)) then
         deallocate(ncdiag_open_status(id)%all_data_fix)
      endif
      if (allocated(ncdiag_open_status(id)%all_data_chan)) then
         deallocate(ncdiag_open_status(id)%all_data_chan)
      endif
      if (allocated(ncdiag_open_status(id)%all_data_extra)) then
         deallocate(ncdiag_open_status(id)%all_data_extra)
      endif
      nopen_ncdiag = nopen_ncdiag - 1
   else
      close(ftin)
   endif

end subroutine close_radiag

subroutine read_radiag_header(ftin,npred_radiag,retrieval,header_fix,header_chan,data_name,iflag,lverbose)
!                .      .    .                                       .
! subprogram:    read_diag_header_bin             read rad diag header
!   prgmmr: mccarty           org: gmao                date: 2015-08-06
!
! abstract:  This routine reads the header record from a radiance
!            diagnostic file
!
! program history log:
!   2015-08-06 mccarty - created routine w/ fork for ncdiag or binary
!
! input argument list:
!   ftin          - unit number connected to diagnostic file 
!   npred_radiag  - number of bias correction terms
!   retrieval     - .true. if sst retrieval
!
! output argument list:
!   header_fix    - header information structure
!   header_chan   - channel information structure
!   data_name     - diag file data names
!   iflag         - error code
!   lverbose      - optional flag to turn off default output to standard out 
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

! Declare passed arguments
  integer(i_kind),intent(in)             :: ftin
  integer(i_kind),intent(in)             :: npred_radiag
  logical,intent(in)                     :: retrieval
  type(diag_header_fix_list ),intent(out):: header_fix
  type(diag_header_chan_list),allocatable :: header_chan(:)
  type(diag_data_name_list)              :: data_name
  integer(i_kind),intent(out)            :: iflag
  logical,optional,intent(in)            :: lverbose

  iflag = 0
  if (netcdf) then
     print *,'netcdf slot'
     call read_radiag_header_nc(ftin,npred_radiag,retrieval,header_fix,header_chan,data_name,iflag,lverbose)
  else
     call read_radiag_header_bin(ftin,npred_radiag,retrieval,header_fix,header_chan,data_name,iflag,lverbose)
  endif

end subroutine read_radiag_header

subroutine read_radiag_header_nc(ftin,npred_radiag,retrieval,header_fix,header_chan,data_name,iflag,lverbose)
!                .      .    .                                       .
! subprogram:    read_diag_header_nc              read rad diag header
!   prgmmr: mccarty           org: gmao                date: 2003-01-01
!
! abstract:  This routine reads the header record from a radiance
!            diagnostic file
!
! program history log:
!   2015-08-06 mccarty - Created routine for ncdiag header reading
!
! input argument list:
!   ftin          - unit number connected to diagnostic file 
!   npred_radiag  - number of bias correction terms
!   retrieval     - .true. if sst retrieval
!
! output argument list:
!   header_fix    - header information structure
!   header_chan   - channel information structure
!   data_name     - diag file data names
!   iflag         - error code
!   lverbose      - optional flag to turn off default output to standard out 
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
! Declare passed arguments
  integer(i_kind),intent(in)             :: ftin
  integer(i_kind),intent(in)             :: npred_radiag
  logical,intent(in)                     :: retrieval
  type(diag_header_fix_list ),intent(out):: header_fix
  type(diag_header_chan_list),allocatable :: header_chan(:)
  type(diag_data_name_list)              :: data_name
  integer(i_kind),intent(out)            :: iflag
  logical,optional,intent(in)            :: lverbose

! local variables
  integer(i_kind)                        :: nchan_dim
  real(r_kind),allocatable,dimension(:)               :: r_var_stor
  integer(i_kind),allocatable,dimension(:)            :: i_var_stor
  character(20)                          :: isis
  character(10)                          :: id, obstype
!  integer(i_kind),dimension(:),allocatable           :: jiter, nchan_diag, npred, idate, &
  integer(i_kind)                        :: jiter, nchan_diag, npred, idate, &
                                            ireal, ipchan, iextra, jextra,   &
                                            idiag, angord, iversion, inewpc, &
                                            isens

  iflag = 0
!  allocate(nchan_diag(1) )
  nchan_dim = nc_diag_read_get_dim(ftin,'nchans')
  write(*,*)'Number of channels=',nchan_dim
  header_fix%nchan = nchan_dim
  
  call nc_diag_read_get_global_attr(ftin, "Number_of_channels", nchan_diag)

  if (nchan_dim .ne. nchan_diag) then
     write(*,*)'ERROR: Number of channels from dimension do not match those from header, aborting.'
     call abort
  endif  

  call nc_diag_read_get_global_attr(ftin, "Satellite_Sensor", isis)      ; header_fix%isis = isis
  call nc_diag_read_get_global_attr(ftin, "Satellite", id)               ; header_fix%id = id
  call nc_diag_read_get_global_attr(ftin, "Observation_type", obstype)   ; header_fix%obstype = obstype
  call nc_diag_read_get_global_attr(ftin, "Outer_Loop_Iteration", jiter) ; header_fix%jiter = jiter
  call nc_diag_read_get_global_attr(ftin, "Number_of_Predictors", npred) ; header_fix%npred = npred
  call nc_diag_read_get_global_attr(ftin, "date_time", idate)            ; header_fix%idate = idate
  call nc_diag_read_get_global_attr(ftin, "ireal_radiag", ireal)         ; header_fix%ireal = ireal
  call nc_diag_read_get_global_attr(ftin, "ipchan_radiag", ipchan)       ; header_fix%ipchan = ipchan
  call nc_diag_read_get_global_attr(ftin, "iextra", iextra)              ; header_fix%iextra = iextra
  call nc_diag_read_get_global_attr(ftin, "jextra", jextra)              ; header_fix%jextra = jextra
  call nc_diag_read_get_global_attr(ftin, "idiag", idiag)                ; header_fix%idiag = idiag
  call nc_diag_read_get_global_attr(ftin, "angord", angord)              ; header_fix%angord = angord
  call nc_diag_read_get_global_attr(ftin, "iversion_radiag", iversion)   ; header_fix%iversion = iversion
  call nc_diag_read_get_global_attr(ftin, "New_pc4pred", inewpc)         ; header_fix%inewpc = inewpc
  call nc_diag_read_get_global_attr(ftin, "ioff0", isens)                ; header_fix%isens = isens


  allocate(header_chan(nchan_dim)   )

  allocate(r_var_stor(nchan_dim), &
           i_var_stor(nchan_dim)  )
!  call nc_diag_read_get_var('Var', var_stor)
  call nc_diag_read_get_var(ftin, 'frequency',r_var_stor)      ; header_chan%freq     = r_var_stor
  call nc_diag_read_get_var(ftin, 'polarization',i_var_stor)   ; header_chan%polar    = i_var_stor
  call nc_diag_read_get_var(ftin, 'wavenumber',r_var_stor)     ; header_chan%wave     = r_var_stor
  call nc_diag_read_get_var(ftin, 'error_variance',r_var_stor) ; header_chan%varch    = r_var_stor
  call nc_diag_read_get_var(ftin, 'mean_lapse_rate',r_var_stor); header_chan%tlapmean = r_var_stor
  call nc_diag_read_get_var(ftin, 'use_flag',i_var_stor)       ; header_chan%iuse     = i_var_stor
  call nc_diag_read_get_var(ftin, 'sensor_chan',i_var_stor)    ; header_chan%nuchan   = i_var_stor
  call nc_diag_read_get_var(ftin, 'satinfo_chan',i_var_stor)   ; header_chan%iochan   = i_var_stor


end subroutine read_radiag_header_nc

subroutine read_radiag_header_bin(ftin,npred_radiag,retrieval,header_fix,header_chan,data_name,iflag,lverbose)
!                .      .    .                                       .
! subprogram:    read_diag_header_bin             read rad diag header
!   prgmmr: tahara           org: np20                date: 2003-01-01
!
! abstract:  This routine reads the header record from a radiance
!            diagnostic file
!
! program history log:
!   2010-10-05 treadon - add this doc block
!   2011-02-22 kleist  - changes related to memory allocation and standard output
!   2014-07-25 sienkiewicz - supress warning if npred_radiag == 0
!   2017-07-17 mccarty - renamed routine to _bin suffix for ncdiag
!
! input argument list:
!   ftin          - unit number connected to diagnostic file 
!   npred_radiag  - number of bias correction terms
!   retrieval     - .true. if sst retrieval
!
! output argument list:
!   header_fix    - header information structure
!   header_chan   - channel information structure
!   data_name     - diag file data names
!   iflag         - error code
!   lverbose      - optional flag to turn off default output to standard out 
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

! Declare passed arguments
  integer(i_kind),intent(in)             :: ftin
  integer(i_kind),intent(in)             :: npred_radiag
  logical,intent(in)                     :: retrieval
  type(diag_header_fix_list ),intent(out):: header_fix
  type(diag_header_chan_list),allocatable :: header_chan(:)
  type(diag_data_name_list)              :: data_name
  integer(i_kind),intent(out)            :: iflag
  logical,optional,intent(in)            :: lverbose    

!  Declare local variables
  character(len=2):: string
  character(len=10):: satid,sentype
  character(len=20):: sensat
  integer(i_kind) :: i,ich
  integer(i_kind):: jiter,nchanl,npred,ianldate,ireal,ipchan,iextra,jextra
  integer(i_kind):: idiag,angord,iversion,inewpc,isens
  integer(i_kind):: iuse_tmp,nuchan_tmp,iochan_tmp
  real(r_single) :: freq_tmp,polar_tmp,wave_tmp,varch_tmp,tlapmean_tmp
  logical loutall

  loutall=.true.
  if(present(lverbose)) loutall=lverbose

! Read header (fixed_part).
  read(ftin,iostat=iflag)  sensat,satid,sentype,jiter,nchanl,npred,ianldate,&
       ireal,ipchan,iextra,jextra,idiag,angord,iversion,inewpc,isens
  if (iflag/=0) then
     rewind(ftin)
     read(ftin,iostat=iflag) sensat,satid,sentype,jiter,nchanl,npred,ianldate,&
          ireal,ipchan,iextra,jextra,idiag,angord,iversion,inewpc
     isens=0
  end if

  if (iflag/=0) then
     rewind(ftin)
     read(ftin,iostat=iflag) sensat,satid,sentype,jiter,nchanl,npred,ianldate,&
          ireal,ipchan,iextra,jextra
     idiag=ipchan+npred+1
     angord=0
     iversion=0
     inewpc=0
     isens=0
     if (iflag/=0) then
        write(6,*)'READ_RADIAG_HEADER:  ***ERROR*** Unknown file format.  Cannot read'
        return
     endif
  endif

  header_fix%isis    = sensat
  header_fix%id      = satid
  header_fix%obstype = sentype
  header_fix%jiter   = jiter
  header_fix%nchan   = nchanl
  header_fix%npred   = npred
  header_fix%idate   = ianldate
  header_fix%ireal   = ireal
  header_fix%ipchan  = ipchan
  header_fix%iextra  = iextra
  header_fix%jextra  = jextra
  header_fix%idiag   = idiag
  header_fix%angord  = angord
  header_fix%iversion= iversion
  header_fix%inewpc  = inewpc
  header_fix%isens   = isens

  if (loutall) then
     write(6,*)'READ_RADIAG_HEADER:  isis=',header_fix%isis,&
          ' nchan=',header_fix%nchan,&
          ' npred=',header_fix%npred,&
          ' angord=',header_fix%angord,&
          ' idiag=',header_fix%idiag,&
          ' iversion=',header_fix%iversion,&
          ' inewpc=',header_fix%inewpc,&
          ' isens=',header_fix%isens
     
     if ( header_fix%iextra /= 0) &
          write(6,*)'READ_RADIAG_HEADER:  extra diagnostic information available, ',&
          'iextra=',header_fix%iextra
  end if

  if (header_fix%npred  /= npred_radiag .and. npred_radiag /= 0) &
       write(6,*) 'READ_RADIAG_HEADER:  **WARNING** header_fix%npred,npred=',&
       header_fix%npred,npred_radiag
  
! Allocate and initialize as needed
     if (allocated(header_chan)) deallocate(header_chan)
     if (allocated(data_name%chn))  deallocate(data_name%chn)

     allocate(header_chan( header_fix%nchan))
     allocate(data_name%chn(header_fix%idiag))

     data_name%fix(1) ='lat       '
     data_name%fix(2) ='lon       '
     data_name%fix(3) ='zsges     '
     data_name%fix(4) ='obstim    '
     data_name%fix(5) ='scanpos   '
     data_name%fix(6) ='satzen    '
     data_name%fix(7) ='satazm    '
     data_name%fix(8) ='solzen    '
     data_name%fix(9) ='solazm    '
     data_name%fix(10)='sungln    '
     data_name%fix(11)='fwater    '
     data_name%fix(12)='fland     '
     data_name%fix(13)='fice      '
     data_name%fix(14)='fsnow     '
     data_name%fix(15)='twater    '
     data_name%fix(16)='tland     '
     data_name%fix(17)='tice      '
     data_name%fix(18)='tsnow     '
     data_name%fix(19)='tsoil     '
     data_name%fix(20)='soilmoi   '
     data_name%fix(21)='landtyp   '
     data_name%fix(22)='vegfrac   '
     data_name%fix(23)='snowdep   '
     data_name%fix(24)='wndspd    '
     data_name%fix(25)='qc1       '
     data_name%fix(26)='qc2       '
     data_name%fix(27)='tref      '
     data_name%fix(28)='dtw       '
     data_name%fix(29)='dtc       '
     data_name%fix(30)='tz_tr     '

     data_name%chn(1)='obs       '
     data_name%chn(2)='omgbc     '
     data_name%chn(3)='omgnbc    '
     data_name%chn(4)='errinv    '
     data_name%chn(5)='qcmark    '
     data_name%chn(6)='emiss     '
     data_name%chn(7)='tlap      '
     data_name%chn(8)='tb_tz     '

     if (header_fix%iversion < iversion_radiag_1) then
        data_name%chn( 8)= 'bifix     '
        data_name%chn( 9)= 'bilap     '
        data_name%chn(10)= 'bilap2    '
        data_name%chn(11)= 'bicons    '
        data_name%chn(12)= 'biang     '
        data_name%chn(13)= 'biclw     '
        if (retrieval) data_name%chn(13)= 'bisst     '
     elseif ( header_fix%iversion < iversion_radiag_2 .and. header_fix%iversion >= iversion_radiag_1 ) then
        data_name%chn( 8)= 'bicons    '
        data_name%chn( 9)= 'biang     '
        data_name%chn(10)= 'biclw     '
        data_name%chn(11)= 'bilap2    '
        data_name%chn(12)= 'bilap     '
        do i=1,header_fix%angord
           write(string,'(i2.2)') header_fix%angord-i+1
           data_name%chn(12+i)= 'bifix' // string
        end do
        data_name%chn(12+header_fix%angord+1)= 'bifix     '
        data_name%chn(12+header_fix%angord+2)= 'bisst     '
     elseif ( header_fix%iversion < iversion_radiag_3 .and. header_fix%iversion >= iversion_radiag_2 ) then
        data_name%chn( 9)= 'bicons    '
        data_name%chn(10)= 'biang     '
        data_name%chn(11)= 'biclw     '
        data_name%chn(12)= 'bilap2    '
        data_name%chn(13)= 'bilap     '
        do i=1,header_fix%angord
           write(string,'(i2.2)') header_fix%angord-i+1
           data_name%chn(13+i)= 'bifix' // string
        end do
        data_name%chn(13+header_fix%angord+1)= 'bifix     '
        data_name%chn(13+header_fix%angord+2)= 'bisst     '
     elseif ( header_fix%iversion < iversion_radiag_4 .and. header_fix%iversion >= iversion_radiag_3 ) then
        data_name%chn( 9)= 'bicons    '
        data_name%chn(10)= 'biang     '
        data_name%chn(11)= 'biclw     '
        data_name%chn(12)= 'bilap2    '
        data_name%chn(13)= 'bilap     '
        data_name%chn(14)= 'bicos     '
        data_name%chn(15)= 'bisin     '
        do i=1,header_fix%angord
           write(string,'(i2.2)') header_fix%angord-i+1
           data_name%chn(15+i)= 'bifix' // string
        end do
        data_name%chn(15+header_fix%angord+1)= 'bifix     '
        data_name%chn(15+header_fix%angord+2)= 'bisst     '
     else
        data_name%chn( 9)= 'bicons    '
        data_name%chn(10)= 'biang     '
        data_name%chn(11)= 'biclw     '
        data_name%chn(12)= 'bilap2    '
        data_name%chn(13)= 'bilap     '
        data_name%chn(14)= 'bicos     '
        data_name%chn(15)= 'bisin     '
        data_name%chn(16)= 'biemis    '
        do i=1,header_fix%angord
           write(string,'(i2.2)') header_fix%angord-i+1
           data_name%chn(16+i)= 'bifix' // string
        end do
        data_name%chn(16+header_fix%angord+1)= 'bifix     '
        data_name%chn(16+header_fix%angord+2)= 'bisst     '
     endif

! Read header (channel part)
  do ich=1, header_fix%nchan
      read(ftin,iostat=iflag) freq_tmp,polar_tmp,wave_tmp,varch_tmp,tlapmean_tmp,iuse_tmp,nuchan_tmp,iochan_tmp
      header_chan(ich)%freq     = freq_tmp
      header_chan(ich)%polar    = polar_tmp
      header_chan(ich)%wave     = wave_tmp
      header_chan(ich)%varch    = varch_tmp
      header_chan(ich)%tlapmean = tlapmean_tmp
      header_chan(ich)%iuse     = iuse_tmp
      header_chan(ich)%nuchan   = nuchan_tmp
      header_chan(ich)%iochan   = iochan_tmp
     if (iflag/=0) return
  end do

! Construct array containing menonics for data record entries
  
     
end subroutine read_radiag_header_bin

integer(i_kind) function find_ncdiag_id(ftin)
  integer, intent(in) :: ftin

  integer :: i

  find_ncdiag_id = -1
  do i = 1, max_open_ncdiag
     if (ncdiag_open_id(i) == ftin) then
        find_ncdiag_id = i
        return
     endif
  enddo
  return

end function find_ncdiag_id

subroutine read_radiag_data(ftin,header_fix,retrieval,data_fix,data_chan,data_extra,iflag )
!                .      .    .                                       .
! subprogram:    read_radiag_dat                    read rad diag data
!   prgmmr: tahara           org: np20                date: 2003-01-01
!
! abstract:  This routine reads the data record from a radiance
!            diagnostic file
!
! program history log:
!   2010-10-05 treadon - add this doc block
!   2011-02-22 kleist  - changes related to memory allocation
!   2017-07-17 mccarty - change routine to be generalized for bin/nc4 files
!
! input argument list:
!   ftin - unit number connected to diagnostic file
!   header_fix - header information structure
!
! output argument list:
!   data_fix   - spot header information structure
!   data_chan  - spot channel information structure
!   data_extra - spot extra information
!   iflag      - error code
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$


! Declare passed arguments
  integer(i_kind),intent(in)             :: ftin
  type(diag_header_fix_list ),intent(in) :: header_fix
  logical,intent(in)                     :: retrieval
  type(diag_data_fix_list)   ,intent(out):: data_fix
  type(diag_data_chan_list)  ,allocatable :: data_chan(:)
  type(diag_data_extra_list) ,allocatable :: data_extra(:,:)
  integer(i_kind),intent(out)            :: iflag

  integer(i_kind) :: id

  if (netcdf) then

     id = find_ncdiag_id(ftin)
     if (id < 0) then
        write(6,*) 'READ_RADIAG_DATA:  ***ERROR*** netcdf diag file ', ftin, ' has not been opened yet.'
        call abort
     endif

     if (.not. ncdiag_open_status(id)%nc_read) then
        call read_radiag_data_nc_init(ftin, ncdiag_open_status(id), header_fix, retrieval)
     endif

     if (ncdiag_open_status(id)%cur_ob_idx .eq. ncdiag_open_status(id)%num_records ) then
       iflag = 0
     else if (ncdiag_open_status(id)%cur_ob_idx .gt. ncdiag_open_status(id)%num_records) then
       iflag = -1
     else
       iflag = 1
     endif

     if (iflag .ge. 0) then
        call read_radiag_data_nc(ftin,ncdiag_open_status(id),header_fix,retrieval,data_fix,data_chan,data_extra,iflag)
     endif

  else
     call read_radiag_data_bin(ftin,header_fix,retrieval,data_fix,data_chan,data_extra,iflag )
  endif

end subroutine read_radiag_data

subroutine read_all_radiag(ftin, header_fix, retrieval, all_data_fix, &
                           all_data_chan, all_data_extra, &
                           nobs, iflag)

  integer(i_kind),intent(in)             :: ftin
  type(diag_header_fix_list ),intent(in) :: header_fix
  logical,intent(in)                     :: retrieval
  integer(i_kind),intent(out)            :: iflag
  type(diag_data_fix_list),  allocatable  :: all_data_fix(:)
  type(diag_data_chan_list), allocatable  :: all_data_chan(:,:)
  type(diag_data_extra_list), allocatable :: all_data_extra(:,:,:)
  integer(i_kind), intent(out) :: nobs

  integer(i_kind) :: id

  if (netcdf) then
     id = find_ncdiag_id(ftin)
     if (id < 0) then
        write(6,*) 'READ_RADIAG_DATA:  ***ERROR*** netcdf diag file ', ftin, ' has not been opened yet.'
        call abort
     endif

     if (.not. ncdiag_open_status(id)%nc_read) then
        call read_radiag_data_nc_init(ftin, ncdiag_open_status(id), header_fix, retrieval)
     endif
     nobs = ncdiag_open_status(id)%num_records

     iflag = 0
     if (.not. allocated(all_data_fix)) allocate(all_data_fix(nobs))
     if (.not. allocated(all_data_chan)) allocate(all_data_chan(nobs,header_fix%nchan) )
     if (.not. allocated(all_data_extra)) allocate(all_data_extra(nobs,header_fix%iextra, header_fix%nchan) )

     all_data_fix   = ncdiag_open_status(id)%all_data_fix
     all_data_chan  = ncdiag_open_status(id)%all_data_chan
     all_data_extra = ncdiag_open_status(id)%all_data_extra
  else
    iflag = -1
  endif
end subroutine read_all_radiag

subroutine read_radiag_data_nc_init(ftin, diag_status, header_fix, retrieval)
!                .      .    .                                       .
! subprogram:    read_radiag_data_nc_init           read rad diag data
!   prgmmr: mccarty          org: np20                date: 2015-08-10
!
! abstract:  This routine reads the data record from a netcdf radiance
!            diagnostic file
!
! program history log:
!   2015-06-10 mccarty  - create routine
!
! input argument list:
!   ftin - unit number connected to diagnostic file
!   header_fix - header information structure
!
! output argument list:
!   data_fix   - spot header information structure
!   data_chan  - spot channel information structure
!   data_extra - spot extra information
!   iflag      - error code
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

! Declare passed arguments
  integer(i_kind),intent(in)             :: ftin
  type(ncdiag_status), intent(inout)     :: diag_status
  type(diag_header_fix_list ),intent(in) :: header_fix
  logical,intent(in)                     :: retrieval

! Declare local variables
  integer(i_kind)                          :: nrecord, ndatum, nangord
  integer(i_kind)                          :: cch, ic, ir, cdatum
  real(r_kind), allocatable, dimension(:)  :: Latitude, Longitude, Elevation, Obs_Time, Scan_Position, Scan_Angle, &
                                              Sat_Zenith_Angle, Sat_Azimuth_Angle, Sol_Zenith_Angle, Sol_Azimuth_Angle,  &
                                              Sun_Glint_Angle, Water_Fraction, Land_Fraction, Ice_Fraction,  &
                                              Snow_Fraction, Water_Temperature, Land_Temperature, Ice_Temperature,  &
                                              Snow_Temperature, Soil_Temperature, Soil_Moisture,  &
                                              tsavg5, sstcu, sstph, sstnv, dta, dqa, dtp_avh, Vegetation_Fraction,  &
                                              Snow_Depth, tpwc_amsua, clw_guess_retrieval, Sfc_Wind_Speed,  &
                                              Cloud_Frac, CTP, CLW, TPWC, clw_obs, clw_guess, Foundation_Temperature, SST_Warm_layer_dt,  &
                                              SST_Cool_layer_tdrop, SST_dTz_dTfound, Observation, Obs_Minus_Forecast_adjusted,  &
                                              Obs_Minus_Forecast_unadjusted, Inverse_Observation_Error, QC_Flag, Emissivity,  &
                                              Weighted_Lapse_Rate, dTb_dTs, BC_Constant, BC_Scan_Angle,  &
                                              BC_Cloud_Liquid_Water, BC_Lapse_Rate_Squared, BC_Lapse_Rate, BC_Cosine_Latitude_times_Node,  &
                                              BC_Sine_Latitude,BC_Emissivity,BC_Fixed_Scan_Position
  integer(i_kind), allocatable, dimension(:)  :: Channel_Index, Land_Type_Index
  real(r_kind), allocatable, dimension(:,:)   :: BC_angord ! (nobs, BC_angord_arr_dim) ;

  real(r_kind)                                :: clat, clon

  ndatum = nc_diag_read_get_dim(ftin,'nobs')
  if (header_fix%angord > 0) then
     nangord = nc_diag_read_get_dim(ftin,'BC_angord_arr_dim')
  end if

  nrecord = ndatum / header_fix%nchan
  diag_status%num_records = nrecord

  allocate( channel_index(ndatum),                                                                                                     &
            latitude(ndatum),                 longitude(ndatum),                        elevation(ndatum),                             &
            obs_time(ndatum),                 scan_position(ndatum),                    sat_zenith_angle(ndatum),                      &
            sat_azimuth_angle(ndatum),        sol_zenith_angle(ndatum),                 sol_azimuth_angle(ndatum),                     &
            sun_glint_angle(ndatum),          water_fraction(ndatum),                   land_fraction(ndatum),                         &
            ice_fraction(ndatum),             snow_fraction(ndatum),                    water_temperature(ndatum),                     &
            land_temperature(ndatum),         ice_temperature(ndatum),                  snow_temperature(ndatum),                      &
            soil_temperature(ndatum),         soil_moisture(ndatum),                    tsavg5(ndatum),                                &
            sstcu(ndatum),                    sstph(ndatum),                            sstnv(ndatum),                                 &
            dta(ndatum),                      dqa(ndatum),                              dtp_avh(ndatum),                               &
            vegetation_fraction(ndatum),      snow_depth(ndatum),                       tpwc_amsua(ndatum),                            & 
            clw_guess_retrieval(ndatum),      sfc_wind_speed(ndatum),                   cloud_frac(ndatum),                            &
            ctp(ndatum),                      clw(ndatum),                              tpwc(ndatum),                                  &
            clw_obs(ndatum),                  clw_guess(ndatum),                        foundation_temperature(ndatum),                &
            sst_warm_layer_dt(ndatum),        sst_cool_layer_tdrop(ndatum),             sst_dtz_dtfound(ndatum),                       &
            observation(ndatum),              obs_minus_forecast_adjusted(ndatum),obs_minus_forecast_unadjusted(ndatum),               &
            inverse_observation_error(ndatum),qc_flag(ndatum),                          emissivity(ndatum),                            &
            weighted_lapse_rate(ndatum),      dtb_dts(ndatum),                          bc_constant(ndatum),                           &
            bc_scan_angle(ndatum),            bc_cloud_liquid_water(ndatum),            bc_lapse_rate_squared(ndatum),                 &
            bc_lapse_rate(ndatum),            bc_cosine_latitude_times_node(ndatum),    bc_sine_latitude(ndatum),                      &
            bc_emissivity(ndatum),            bc_fixed_scan_position(ndatum),           land_type_index(ndatum),                       &
            scan_angle(ndatum)                        )

  if (header_fix%angord > 0) then
     allocate( bc_angord(nangord, ndatum)    )
  end if

  if (allocated(diag_status%all_data_fix))   deallocate(diag_status%all_data_fix)
  if (allocated(diag_status%all_data_chan))  deallocate(diag_status%all_data_chan)
  if (allocated(diag_status%all_data_extra)) deallocate(diag_status%all_data_extra)
  allocate(  diag_status%all_data_fix(nrecord)        )
  allocate(  diag_status%all_data_chan(nrecord, header_fix%nchan))
  allocate(  diag_status%all_data_extra(nrecord, header_fix%iextra, header_fix%jextra) )

  call nc_diag_read_get_var(ftin, 'Channel_Index', channel_index)
  call nc_diag_read_get_var(ftin, 'Latitude', latitude)
  call nc_diag_read_get_var(ftin, 'Longitude', longitude)
  call nc_diag_read_get_var(ftin, 'Elevation', elevation)
  call nc_diag_read_get_var(ftin, 'Obs_Time', obs_time)
  call nc_diag_read_get_var(ftin, 'Scan_Position', scan_position)
  call nc_diag_read_get_var(ftin, 'Scan_Angle', scan_angle)
  call nc_diag_read_get_var(ftin, 'Sat_Zenith_Angle', sat_zenith_angle)
  call nc_diag_read_get_var(ftin, 'Sat_Azimuth_Angle', sat_azimuth_angle)
  call nc_diag_read_get_var(ftin, 'Sol_Zenith_Angle', sol_zenith_angle)
  call nc_diag_read_get_var(ftin, 'Sol_Azimuth_Angle', sol_azimuth_angle)
  call nc_diag_read_get_var(ftin, 'Sun_Glint_Angle', sun_glint_angle)
  call nc_diag_read_get_var(ftin, 'Water_Fraction', water_fraction)
  call nc_diag_read_get_var(ftin, 'Land_Fraction', land_fraction)
  call nc_diag_read_get_var(ftin, 'Ice_Fraction', ice_fraction)
  call nc_diag_read_get_var(ftin, 'Snow_Fraction', snow_fraction)
  call nc_diag_read_get_var(ftin, 'Water_Temperature', water_temperature)
  call nc_diag_read_get_var(ftin, 'Land_Temperature', land_temperature)
  call nc_diag_read_get_var(ftin, 'Ice_Temperature', ice_temperature)
  call nc_diag_read_get_var(ftin, 'Snow_Temperature', snow_temperature)
  call nc_diag_read_get_var(ftin, 'Soil_Temperature', soil_temperature)
  call nc_diag_read_get_var(ftin, 'Soil_Moisture', soil_moisture)
  call nc_diag_read_get_var(ftin, 'tsavg5', tsavg5)
  call nc_diag_read_get_var(ftin, 'sstcu', sstcu)
  call nc_diag_read_get_var(ftin, 'sstph', sstph)
  call nc_diag_read_get_var(ftin, 'sstnv', sstnv)
  call nc_diag_read_get_var(ftin, 'dta', dta)
  call nc_diag_read_get_var(ftin, 'dqa', dqa)
  call nc_diag_read_get_var(ftin, 'dtp_avh', dtp_avh)
  call nc_diag_read_get_var(ftin, 'Vegetation_Fraction', vegetation_fraction)
  call nc_diag_read_get_var(ftin, 'Snow_Depth', snow_depth)
  call nc_diag_read_get_var(ftin, 'tpwc_amsua', tpwc_amsua)
  call nc_diag_read_get_var(ftin, 'clw_guess_retrieval', clw_guess_retrieval)
  call nc_diag_read_get_var(ftin, 'Sfc_Wind_Speed', sfc_wind_speed)
  call nc_diag_read_get_var(ftin, 'Cloud_Frac', cloud_frac)
  call nc_diag_read_get_var(ftin,'CTP', ctp)
  call nc_diag_read_get_var(ftin, 'CLW', clw)
  call nc_diag_read_get_var(ftin, 'TPWC', tpwc)
  call nc_diag_read_get_var(ftin, 'clw_obs', clw_obs)
  call nc_diag_read_get_var(ftin, 'clw_guess', clw_guess)
  call nc_diag_read_get_var(ftin, 'Foundation_Temperature', foundation_temperature)
  call nc_diag_read_get_var(ftin, 'SST_Warm_layer_dt', sst_warm_layer_dt)
  call nc_diag_read_get_var(ftin, 'SST_Cool_layer_tdrop', sst_cool_layer_tdrop)
  call nc_diag_read_get_var(ftin, 'SST_dTz_dTfound', sst_dtz_dtfound)
  call nc_diag_read_get_var(ftin, 'Observation', observation)
  call nc_diag_read_get_var(ftin, 'Obs_Minus_Forecast_adjusted', obs_minus_forecast_adjusted)
  call nc_diag_read_get_var(ftin, 'Obs_Minus_Forecast_unadjusted', obs_minus_forecast_unadjusted)
  call nc_diag_read_get_var(ftin, 'Inverse_Observation_Error', inverse_observation_error)
  call nc_diag_read_get_var(ftin, 'QC_Flag', qc_flag)
  call nc_diag_read_get_var(ftin, 'Emissivity', emissivity)
  call nc_diag_read_get_var(ftin, 'Weighted_Lapse_Rate', weighted_lapse_rate)
  call nc_diag_read_get_var(ftin, 'dTb_dTs', dtb_dts)
  call nc_diag_read_get_var(ftin, 'BC_Constant', bc_constant)
  call nc_diag_read_get_var(ftin, 'BC_Scan_Angle', bc_scan_angle)
  call nc_diag_read_get_var(ftin, 'BC_Cloud_Liquid_Water', bc_cloud_liquid_water)
  call nc_diag_read_get_var(ftin, 'BC_Lapse_Rate_Squared', bc_lapse_rate_squared)
  call nc_diag_read_get_var(ftin, 'BC_Lapse_Rate', bc_lapse_rate)
  call nc_diag_read_get_var(ftin, 'BC_Cosine_Latitude_times_Node', bc_cosine_latitude_times_node)
  call nc_diag_read_get_var(ftin, 'BC_Sine_Latitude', bc_sine_latitude)
  call nc_diag_read_get_var(ftin, 'BC_Emissivity', bc_emissivity)
  call nc_diag_read_get_var(ftin, 'BC_Fixed_Scan_Position', bc_fixed_scan_position)
  call nc_diag_read_get_var(ftin, 'Land_Type_Index', land_type_index)
  if (header_fix%angord > 0) then
     call nc_diag_read_get_var(ftin, 'BC_angord ', bc_angord )
  end if
  cdatum = 1

!  allocate(  all_data_fix(nrecord)        )
!  allocate(  all_data_chan(nrecord, nchan))

  
  do ir=1,nrecord
    clat = latitude(cdatum)
    clon = longitude(cdatum)
    diag_status%all_data_fix(ir)%lat               = latitude(cdatum)
    diag_status%all_data_fix(ir)%lon               = longitude(cdatum)     
    diag_status%all_data_fix(ir)%zsges             = elevation(cdatum)
    diag_status%all_data_fix(ir)%obstime           = obs_time(cdatum)
    diag_status%all_data_fix(ir)%senscn_pos        = scan_position(cdatum)
    diag_status%all_data_fix(ir)%senscn_ang        = scan_angle(cdatum)
    diag_status%all_data_fix(ir)%satzen_ang        = sat_zenith_angle(cdatum)
    diag_status%all_data_fix(ir)%satazm_ang        = sat_azimuth_angle(cdatum)
    diag_status%all_data_fix(ir)%solzen_ang        = sol_zenith_angle(cdatum)
    diag_status%all_data_fix(ir)%solazm_ang        = sol_azimuth_angle(cdatum)
    diag_status%all_data_fix(ir)%sungln_ang        = sun_glint_angle(cdatum)
    diag_status%all_data_fix(ir)%water_frac        = water_fraction(cdatum)
    diag_status%all_data_fix(ir)%land_frac         = land_fraction(cdatum)
    diag_status%all_data_fix(ir)%ice_frac          = ice_fraction(cdatum)
    diag_status%all_data_fix(ir)%snow_frac         = snow_fraction(cdatum)
    diag_status%all_data_fix(ir)%water_temp        = water_temperature(cdatum)
    diag_status%all_data_fix(ir)%land_temp         = land_temperature(cdatum)
    diag_status%all_data_fix(ir)%ice_temp          = ice_temperature(cdatum)
    diag_status%all_data_fix(ir)%snow_temp         = snow_temperature(cdatum)
    diag_status%all_data_fix(ir)%soil_temp         = soil_temperature(cdatum)
    diag_status%all_data_fix(ir)%soil_mois         = soil_moisture(cdatum)
    diag_status%all_data_fix(ir)%land_type         = land_type_index(cdatum)
    diag_status%all_data_fix(ir)%veg_frac          = vegetation_fraction(cdatum)
    diag_status%all_data_fix(ir)%snow_depth        = snow_depth(cdatum)
    diag_status%all_data_fix(ir)%sfc_wndspd        = sfc_wind_speed(cdatum)
    diag_status%all_data_fix(ir)%qcdiag1           = cloud_frac(cdatum)
    diag_status%all_data_fix(ir)%qcdiag2           = ctp(cdatum)
    diag_status%all_data_fix(ir)%tref              = foundation_temperature(cdatum)
    diag_status%all_data_fix(ir)%dtw               = sst_warm_layer_dt(cdatum)
    diag_status%all_data_fix(ir)%dtc               = sst_cool_layer_tdrop(cdatum)
    diag_status%all_data_fix(ir)%tz_tr             = sst_dtz_dtfound(cdatum)

    if (retrieval) then
      diag_status%all_data_fix(ir)%water_temp        = tsavg5(cdatum)
      diag_status%all_data_fix(ir)%land_temp         = sstcu(cdatum)
      diag_status%all_data_fix(ir)%ice_temp          = sstph(cdatum)
      diag_status%all_data_fix(ir)%snow_temp         = sstnv(cdatum)
      diag_status%all_data_fix(ir)%soil_temp         = dta(cdatum)
      diag_status%all_data_fix(ir)%soil_mois         = dqa(cdatum)
      diag_status%all_data_fix(ir)%land_type         = dtp_avh(cdatum)
    endif

    do ic=1,header_fix%nchan
      if (clat .ne. latitude(cdatum) .or. clon .ne. longitude(cdatum)) then
        write(*,*) 'ERROR: Lats & Lons are mismatched.  This is bad'
        print *,'irecord=',ir
        print *,'clat,clon=',clat,clon
        print *,'lat/lon(datum)=',latitude(cdatum), longitude(cdatum)
        call abort
      endif
      cch = channel_index(cdatum)
      if (allocated(diag_status%all_data_chan(ir,cch)%bifix)) deallocate(diag_status%all_data_chan(ir,cch)%bifix )
      if (header_fix%angord > 0) then
         allocate(diag_status%all_data_chan(ir,cch)%bifix(nangord))
      else
         allocate(diag_status%all_data_chan(ir,cch)%bifix(1))
      end if

      diag_status%all_data_chan(ir,cch)%tbobs = observation(cdatum)
      diag_status%all_data_chan(ir,cch)%omgbc = obs_minus_forecast_adjusted(cdatum)
      diag_status%all_data_chan(ir,cch)%omgnbc= obs_minus_forecast_unadjusted(cdatum)
      diag_status%all_data_chan(ir,cch)%errinv= inverse_observation_error(cdatum)
      diag_status%all_data_chan(ir,cch)%qcmark= qc_flag(cdatum)
      diag_status%all_data_chan(ir,cch)%emiss = emissivity(cdatum)
      diag_status%all_data_chan(ir,cch)%tlap  = weighted_lapse_rate(cdatum)
      diag_status%all_data_chan(ir,cch)%tb_tz = dtb_dts(cdatum)
      diag_status%all_data_chan(ir,cch)%bicons= bc_constant(cdatum)
      diag_status%all_data_chan(ir,cch)%biang = bc_scan_angle(cdatum)
      diag_status%all_data_chan(ir,cch)%biclw = bc_cloud_liquid_water(cdatum)
      diag_status%all_data_chan(ir,cch)%bilap2= bc_lapse_rate_squared(cdatum)
      diag_status%all_data_chan(ir,cch)%bilap = bc_lapse_rate(cdatum)
      diag_status%all_data_chan(ir,cch)%bicos = bc_cosine_latitude_times_node(cdatum)
      diag_status%all_data_chan(ir,cch)%bisin = bc_sine_latitude(cdatum)
      diag_status%all_data_chan(ir,cch)%biemis= bc_emissivity(cdatum)
      if (header_fix%angord > 0) then
         diag_status%all_data_chan(ir,cch)%bifix = bc_angord(1:nangord,cdatum)
      else
         diag_status%all_data_chan(ir,cch)%bifix(1) = bc_fixed_scan_position(cdatum)
      endif
      ! placeholder for SST BC
      cdatum = cdatum + 1
    enddo
  enddo

  diag_status%nc_read = .true.
  diag_status%cur_ob_idx = 1
end subroutine read_radiag_data_nc_init

subroutine read_radiag_data_nc(ftin,diag_status,header_fix,retrieval,data_fix,data_chan,data_extra,iflag )
!                .      .    .                                       .
! subprogram:    read_radiag_dat                    read rad diag data
!   prgmmr: tahara           org: np20                date: 2015-08-10
!
! abstract:  This routine reads the data record from a netcdf radiance
!            diagnostic file
!
! program history log:
!   2015-08-10 mccarty  - create routine
!
! input argument list:
!   ftin - unit number connected to diagnostic file
!   header_fix - header information structure
!
! output argument list:
!   data_fix   - spot header information structure
!   data_chan  - spot channel information structure
!   data_extra - spot extra information
!   iflag      - error code
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$

! Declare passed arguments
  integer(i_kind),intent(in)             :: ftin
  type(ncdiag_status), intent(inout)     :: diag_status
  type(diag_header_fix_list ),intent(in) :: header_fix
  logical,intent(in)                     :: retrieval
  type(diag_data_fix_list)   ,intent(out):: data_fix
  type(diag_data_chan_list)  ,allocatable :: data_chan(:)
  type(diag_data_extra_list) ,allocatable :: data_extra(:,:)
  integer(i_kind),intent(out)            :: iflag

  iflag = 0
  if (.not. allocated(data_chan)) allocate(data_chan(header_fix%nchan) )
  if (.not. allocated(data_extra)) allocate(data_extra(header_fix%iextra, header_fix%nchan) )

  data_fix        = diag_status%all_data_fix(diag_status%cur_ob_idx)
  data_chan(:)    = diag_status%all_data_chan(diag_status%cur_ob_idx,:)
  data_extra(:,:) = diag_status%all_data_extra(diag_status%cur_ob_idx,:,:) 

  diag_status%cur_ob_idx = diag_status%cur_ob_idx + 1

end subroutine read_radiag_data_nc

subroutine read_radiag_data_bin(ftin,header_fix,retrieval,data_fix,data_chan,data_extra,iflag )
!                .      .    .                                       .
! subprogram:    read_radiag_dat                    read rad diag data
!   prgmmr: tahara           org: np20                date: 2003-01-01
!
! abstract:  This routine reads the data record from a radiance
!            diagnostic file
!
! program history log:
!   2010-10-05 treadon - add this doc block
!   2011-02-22 kleist  - changes related to memory allocation
!   2017-07-17 mccarty - rename binary-specific procedure
!
! input argument list:
!   ftin - unit number connected to diagnostic file
!   header_fix - header information structure
!
! output argument list:
!   data_fix   - spot header information structure
!   data_chan  - spot channel information structure
!   data_extra - spot extra information
!   iflag      - error code
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$


! Declare passed arguments
  integer(i_kind),intent(in)             :: ftin
  type(diag_header_fix_list ),intent(in) :: header_fix
  logical,intent(in)                     :: retrieval
  type(diag_data_fix_list)   ,intent(out):: data_fix
  type(diag_data_chan_list)  ,allocatable :: data_chan(:)
  type(diag_data_extra_list) ,allocatable :: data_extra(:,:)
  integer(i_kind),intent(out)            :: iflag
    
  integer(i_kind) :: ich,iang,i,j
  real(r_single),dimension(:,:),allocatable :: data_tmp
  real(r_single),dimension(:),allocatable   :: fix_tmp
  real(r_single),dimension(:,:),allocatable :: extra_tmp

! Allocate arrays as needed
  if (allocated(data_chan)) deallocate(data_chan)
  allocate(data_chan(header_fix%nchan))

  do ich=1,header_fix%nchan
     if (allocated(data_chan(ich)%bifix)) deallocate(data_chan(ich)%bifix)
     allocate(data_chan(ich)%bifix(header_fix%angord+1))
  end do

  if (header_fix%iextra > 0) then
     if (allocated(data_extra))   deallocate(data_extra)
     allocate(data_extra(header_fix%iextra,header_fix%jextra))
     allocate(extra_tmp(header_fix%iextra,header_fix%jextra))
  end if

! Allocate arrays to hold data record
  allocate(data_tmp(header_fix%idiag,header_fix%nchan))

  if (header_fix%iversion < iversion_radiag_2) then
     allocate( fix_tmp( ireal_old_radiag ) )
  else
     allocate( fix_tmp( ireal_radiag ) )
  end if

! Read data record

  if (header_fix%iextra == 0) then
     read(ftin,iostat=iflag) fix_tmp, data_tmp
  else
     read(ftin,iostat=iflag) fix_tmp, data_tmp, extra_tmp
  endif

  if (iflag /= 0) return

! Transfer fix_tmp record to output structure
  data_fix%lat = fix_tmp(1)
  data_fix%lon = fix_tmp(2)
  data_fix%zsges = fix_tmp(3)
  data_fix%obstime = fix_tmp(4) 
  data_fix%senscn_pos = fix_tmp(5)
  data_fix%satzen_ang = fix_tmp(6)
  data_fix%satazm_ang = fix_tmp(7)
  data_fix%solzen_ang = fix_tmp(8)
  data_fix%solazm_ang = fix_tmp(9)
  data_fix%sungln_ang = fix_tmp(10)
  data_fix%water_frac = fix_tmp(11)
  data_fix%land_frac = fix_tmp(12)
  data_fix%ice_frac = fix_tmp(13)
  data_fix%snow_frac = fix_tmp(14)
  data_fix%water_temp = fix_tmp(15)
  data_fix%land_temp = fix_tmp(16)
  data_fix%ice_temp = fix_tmp(17)
  data_fix%snow_temp = fix_tmp(18)
  data_fix%soil_temp = fix_tmp(19)
  data_fix%soil_mois = fix_tmp(20)
  data_fix%land_type = fix_tmp(21)
  data_fix%veg_frac = fix_tmp(22)
  data_fix%snow_depth = fix_tmp(23)
  data_fix%sfc_wndspd = fix_tmp(24)
  data_fix%qcdiag1 = fix_tmp(25)
  data_fix%qcdiag2 = fix_tmp(26)

  if ( header_fix%iversion <= iversion_radiag_1 ) then
     data_fix%tref = rmiss_radiag
     data_fix%dtw = rmiss_radiag
     data_fix%dtc = rmiss_radiag
     data_fix%tz_tr = rmiss_radiag
  else
     data_fix%tref = fix_tmp(27)
     data_fix%dtw = fix_tmp(28)
     data_fix%dtc = fix_tmp(29)
     data_fix%tz_tr = fix_tmp(30)
  end if


! Transfer data record to output structure
  do ich=1,header_fix%nchan
     data_chan(ich)%tbobs =data_tmp(1,ich)
     data_chan(ich)%omgbc =data_tmp(2,ich)
     data_chan(ich)%omgnbc=data_tmp(3,ich)
     data_chan(ich)%errinv=data_tmp(4,ich)
     data_chan(ich)%qcmark=data_tmp(5,ich)
     data_chan(ich)%emiss =data_tmp(6,ich)
     data_chan(ich)%tlap  =data_tmp(7,ich)
     data_chan(ich)%tb_tz =data_tmp(8,ich)
  end do
  if (header_fix%iversion < iversion_radiag_1) then
     do ich=1,header_fix%nchan
        data_chan(ich)%bifix(1)=data_tmp(8,ich)
        data_chan(ich)%bilap   =data_tmp(9,ich)
        data_chan(ich)%bilap2  =data_tmp(10,ich)
        data_chan(ich)%bicons  =data_tmp(11,ich)
        data_chan(ich)%biang   =data_tmp(12,ich)
        data_chan(ich)%biclw   =data_tmp(13,ich)
        data_chan(ich)%bisst   = rmiss_radiag
        if (retrieval) then
           data_chan(ich)%biclw   =rmiss_radiag
           data_chan(ich)%bisst   =data_tmp(13,ich) 
        endif
     end do
  elseif ( header_fix%iversion < iversion_radiag_2 .and. header_fix%iversion >= iversion_radiag_1 ) then
     do ich=1,header_fix%nchan
        data_chan(ich)%bicons=data_tmp(8,ich)
        data_chan(ich)%biang =data_tmp(9,ich)
        data_chan(ich)%biclw =data_tmp(10,ich)
        data_chan(ich)%bilap2=data_tmp(11,ich)
        data_chan(ich)%bilap =data_tmp(12,ich)
     end do
     do ich=1,header_fix%nchan
        do iang=1,header_fix%angord+1
           data_chan(ich)%bifix(iang)=data_tmp(12+iang,ich)
        end do
        data_chan(ich)%bisst = data_tmp(12+header_fix%angord+2,ich)  
     end do
  elseif ( header_fix%iversion < iversion_radiag_3 .and. header_fix%iversion >= iversion_radiag_2 ) then
     do ich=1,header_fix%nchan
        data_chan(ich)%bicons=data_tmp(9,ich)
        data_chan(ich)%biang =data_tmp(10,ich)
        data_chan(ich)%biclw =data_tmp(11,ich)
        data_chan(ich)%bilap2=data_tmp(12,ich)
        data_chan(ich)%bilap =data_tmp(13,ich)
     end do
     do ich=1,header_fix%nchan
        do iang=1,header_fix%angord+1
           data_chan(ich)%bifix(iang)=data_tmp(13+iang,ich)
        end do
        data_chan(ich)%bisst = data_tmp(13+header_fix%angord+2,ich)
     end do
  elseif ( header_fix%iversion < iversion_radiag_4 .and. header_fix%iversion >= iversion_radiag_3 ) then
     do ich=1,header_fix%nchan
        data_chan(ich)%bicons=data_tmp(9,ich)
        data_chan(ich)%biang =data_tmp(10,ich)
        data_chan(ich)%biclw =data_tmp(11,ich)
        data_chan(ich)%bilap2=data_tmp(12,ich)
        data_chan(ich)%bilap =data_tmp(13,ich)
        data_chan(ich)%bicos =data_tmp(14,ich) ! 1st bias correction term node*cos(lat) for SSMIS
        data_chan(ich)%bisin =data_tmp(15,ich) ! 2nd bias correction term sin(lat) for SSMI      
     end do
     do ich=1,header_fix%nchan
        do iang=1,header_fix%angord+1
           data_chan(ich)%bifix(iang)=data_tmp(15+iang,ich)
        end do
        data_chan(ich)%bisst = data_tmp(15+header_fix%angord+2,ich)
     end do
  else
     do ich=1,header_fix%nchan
        data_chan(ich)%bicons=data_tmp(9,ich)
        data_chan(ich)%biang =data_tmp(10,ich)
        data_chan(ich)%biclw =data_tmp(11,ich)
        data_chan(ich)%bilap2=data_tmp(12,ich)
        data_chan(ich)%bilap =data_tmp(13,ich)
        data_chan(ich)%bicos =data_tmp(14,ich)
        data_chan(ich)%bisin =data_tmp(15,ich)
        data_chan(ich)%biemis=data_tmp(16,ich)
     end do
     do ich=1,header_fix%nchan
        do iang=1,header_fix%angord+1
           data_chan(ich)%bifix(iang)=data_tmp(16+iang,ich)
        end do
        data_chan(ich)%bisst = data_tmp(16+header_fix%angord+2,ich)  
     end do

  endif

  if (header_fix%iextra > 0) then
     do j=1,header_fix%jextra
        do i=1,header_fix%iextra
           data_extra(i,j)%extra=extra_tmp(i,j)
        end do
     end do
  endif

  deallocate(data_tmp, fix_tmp)
  if (header_fix%iextra > 0) deallocate(extra_tmp)

end subroutine read_radiag_data_bin

end module read_diag