module read_obsmod
!$$$   module documentation block
!                .      .    .                                       .
! module:    read_obsmod extra inquire routine for reading obs
!   prgmmr: parrish          org: np22                date: 2005-06-06
!
! abstract:
!
! program history log:
!   2009-01-05  todling - add gsi_inquire
!   2015-05-01  Liu Ling - Add call to read_rapidscat 
!   2015-08-20  zhu  - add flexibility for enabling all-sky and using aerosol info in radiance 
!                      assimilation. Use radiance_obstype_search from radiance_mod.  
!   2017-05-12  Y. Wang and X. Wang - add dbz to be read in, POC: xuguang.wang@ou.edu
!
! subroutines included:
!   sub gsi_inquire   -  inquire statement supporting fortran earlier than 2003
!   sub read_obs      -  read, select, reformat obs data
!
! Variable Definitions:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

! set default to private
  private
! set subroutines to public
  public :: gsi_inquire
  public :: read_obs

contains

subroutine gsi_inquire (lbytes,lexist,filename,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    gsi_inquire        inquire file presence and size
!   prgmmr: todling      org: np22                date: 2009-01-05
!
! abstract:  Inquire file presence and size; to be used when fortran
!            2003 not available or non-compliant.
!
! program history log:
!   2009-01-05  todling
!   2013-05-14  guo     - changed the compiler #ifdef dependency on _size_
!                         inquire to a more portable way.
!
!   input argument list:
!     mype     - mpi task id
!    filename  - input filename
!
!   output argument list:
!    lexist     - file presence flag
!    lbytes     - file size (bytes)
!
! attributes:
!   language: f90
!   machine:  Linux-cluster
!
!$$$  end documentation block

  use kinds, only: i_kind,i_llong

  implicit none

  integer(i_llong),intent(  out) :: lbytes
  logical         ,intent(  out) :: lexist
  character(len=*),intent(in   ) :: filename
  integer(i_kind) ,intent(in   ) :: mype

  character(len=256) command, fname

#if defined(__INTEL_COMPILER) && (__INTEL_COMPILER < 1110)
#define __X_OR_Y_OR_Z_FORTRAN_COMPILERS__
#endif

#ifdef __X_OR_Y_OR_Z_FORTRAN_COMPILERS__
#define _DO_NOT_SUPPORT_SIZE_INQUIRE_
#endif

  lbytes=-1  ! in case that _size_ specifier is not supported.
#ifndef _DO_NOT_SUPPORT_SIZE_INQUIRE_
  inquire(file=trim(filename),exist=lexist,size=lbytes)
  ! Note that the value of _size_ is defined by Fortran in "file storage units",
  ! which is not neccesary in byte units.  It is not clear if this code had
  ! assumed the size value to be in byte units, or in whatever units.
#else
  inquire(file=trim(filename),exist=lexist)
#endif
  if(lexist) then
     ! Even with a compiler supporting 'size=' specifier, the size value may
     ! return -1, if a compiler considers that the size can not be determined.
     ! In that case, the size may be obtained through a user supported
     ! mechanism, such as reading the size from a system("wc -c") call result.
     if(lbytes<0) then
        write(fname,'(2a,i4.4)') 'fsize_',trim(filename),mype
        write(command,'(4a)') 'wc -c ', trim(filename),' > ', trim(fname)
        call system(command)
        open(unit=999,file=trim(fname),form='formatted')
        read(999,*) lbytes
        close(999)
        lexist = lbytes>0_i_llong ! skip this file if lbytes <=0
     endif
  endif

  return
end subroutine gsi_inquire

subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse,nread)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_obs_check     inquire file presence and size
!   prgmmr: todling      org: np22                date: 2010-03-05
!
! abstract:  Reset file status depending on whether observation time
!            matches analysis time and how offtime_date is set. This
!            also checks for consistency in satellite data files and 
!            known types.  
!            WARNING: some of it looks inconsistent with long-window 4dvar
!
! program history log:
!   2009-xx-xx  derber   - originally placed inside inquire
!   2009-01-05  todling  - move time/type-check out of inquire
!   2010-09-13  pagowski - add anow bufr and one obs chem
!   2013-01-26  parrish - WCOSS debug compile fails with satid not initialized.
!                         Set satid=1 at start of subroutine to allow debug compile.
!   2013-02-13  eliu     - add ssmis 
!   2013-07-01  todling/guo - allow user to bypass this check (old bufr support)
!   2014-10-01  ejones   - add gmi and amsr2
!   2015-01-16  ejones   - add saphir
!   2016-09-19  guo      - properly initialized nread, in case of for quick-return cases.
!   2017-11-16  dutta    - adding KOMPSAT5 bufr i.d for reading the data.
!   2019-03-27  h. liu   - add abi
!                           
!
!   input argument list:
!    lexist    - file status
!    filename  - input filename
!    jsatid    - satellite id
!    dtype     - satellite type
!
!   output argument list:
!    lexist    - file status
!
! attributes:
!   language: f90
!   machine:  Linux-cluster
!
!$$$  end documentation block

  use kinds, only: i_kind,i_llong,r_kind,r_double
  use gsi_4dvar, only: iadatebgn,iadateend
  use obsmod, only: offtime_data
  use convinfo, only: nconvtype,ictype,ioctype,icuse
  use chemmod, only : oneobtest_chem,oneob_type_chem,&
       code_pm25_ncbufr,code_pm25_anowbufr,code_pm10_ncbufr,code_pm10_anowbufr

  implicit none

  logical         ,intent(inout)  :: lexist
  character(len=*),intent(in)     :: filename
  character(len=*),intent(in)     :: jsatid
  character(len=*),intent(in)     :: dtype
  integer(i_kind) ,intent(in)     :: minuse
  integer(i_kind) ,intent(inout)  :: nread

  integer(i_kind) :: lnbufr,idate,idate2,iret,kidsat
  integer(i_kind) :: ireadsb,ireadmg,kx,nc,said
  real(r_double) :: satid,rtype
  character(8) subset

  satid=1      ! debug executable wants default value ???
  idate=0
  nread=0; if(lexist) nread=1   ! in case of a quick return
#ifdef _SKIP_READ_OBS_CHECK_
  return
#endif
  if(trim(dtype) == 'tcp' .or. trim(filename) == 'tldplrso')return
  if(trim(filename) == 'mitmdat' .or. trim(filename) == 'mxtmdat')return
  if(trim(filename) == 'satmar')return
  if(trim(dtype) == 'dbz' )return

! Use routine as usual

  if(lexist .and. trim(dtype) /= 'tcp' )then
      lnbufr = 15
      open(lnbufr,file=trim(filename),form='unformatted',status ='unknown')
      call openbf(lnbufr,'IN',lnbufr)
      call datelen(10)
      call readmg(lnbufr,subset,idate,iret)
      if(iret == 0)then

!        Extract date and check for consistency with analysis date
         if (idate<iadatebgn.or.idate>iadateend) then
            if(offtime_data) then
              write(6,*)'***read_obs_check analysis and data file date differ, but use anyway'
            else
               write(6,*)'***read_obs_check*** ',&
                 'incompatable analysis and observation date/time',trim(filename),trim(dtype)
               lexist=.false.
            end if
            write(6,*)'Analysis start  :',iadatebgn
            write(6,*)'Analysis end    :',iadateend
            write(6,*)'Observation time:',idate
        endif
      else
         write(6,*)'***read_obs_check*** iret/=0 for reading date for ',trim(filename),dtype,jsatid,iret
         lexist=.false.
      end if
      if(lexist)then
       if(jsatid == '')then
         kidsat=0
       else if(jsatid == 'metop-a')then
         kidsat=4
       else if(jsatid == 'metop-b')then
         kidsat=3
       else if(jsatid == 'metop-c')then
         kidsat=5
       else if(jsatid == 'm08')then
         kidsat = 55 
       else if(jsatid == 'm09')then
         kidsat = 56 
       else if(jsatid == 'm10')then
         kidsat = 57 
       else if(jsatid == 'm11')then
         kidsat = 70 
       else if(jsatid == 'n08')then
         kidsat=200
       else if(jsatid == 'n09')then
         kidsat=201
       else if(jsatid == 'n10')then
         kidsat=202
       else if(jsatid == 'n11')then
         kidsat=203
       else if(jsatid == 'n12')then
         kidsat=204
       else if(jsatid == 'n14')then
         kidsat=205
       else if(jsatid == 'n15')then
         kidsat=206
       else if(jsatid == 'n16')then
         kidsat=207
       else if(jsatid == 'n17')then
         kidsat=208
       else if(jsatid == 'n18')then
         kidsat=209
       else if(jsatid == 'n19')then
         kidsat=223
       else if(jsatid == 'npp')then
         kidsat=224
       else if(jsatid == 'n20')then
         kidsat=225
       else if(jsatid == 'n21')then
         kidsat=226
       else if(jsatid == 'f08')then
         kidsat=241
       else if(jsatid == 'f10')then
         kidsat=243
       else if(jsatid == 'f11')then
         kidsat=244
       else if(jsatid == 'f13')then
         kidsat=246
       else if(jsatid == 'f14')then
         kidsat=247
       else if(jsatid == 'f15')then
         kidsat=248
       else if(jsatid == 'f16')then
         kidsat=249    
       else if(jsatid == 'trmm')then
         kidsat=282    
       else if(jsatid == 'f17')then
         kidsat=285                  
       else if(jsatid == 'f18')then  
         kidsat=286                  
       else if(jsatid == 'f19')then  
         kidsat=287                  
       else if(jsatid == 'g08' .or. jsatid == 'g08_prep')then
         kidsat=252
       else if(jsatid == 'g09' .or. jsatid == 'g09_prep')then
         kidsat=253
       else if(jsatid == 'g10' .or. jsatid == 'g10_prep')then
         kidsat=254
       else if(jsatid == 'g11' .or. jsatid == 'g11_prep')then
         kidsat=255
       else if(jsatid == 'g12' .or. jsatid == 'g12_prep')then
         kidsat=256
       else if(jsatid == 'g13' .or. jsatid == 'g13_prep')then
         kidsat=257
       else if(jsatid == 'g14' .or. jsatid == 'g14_prep')then
         kidsat=258
       else if(jsatid == 'g15' .or. jsatid == 'g15_prep')then
         kidsat=259
       else if(jsatid == 'g16' .or. jsatid == 'g16_prep')then
         kidsat=270
       else if(jsatid == 'g17' .or. jsatid == 'g17_prep')then
         kidsat=271
       else if(jsatid == 'n05')then
         kidsat=705
       else if(jsatid == 'n06')then
         kidsat=706
       else if(jsatid == 'n07')then
         kidsat=707
       else if(jsatid == 'tirosn')then
         kidsat=708
       else if ( jsatid == 'terra' ) then
         kidsat = 783
       else if ( jsatid == 'aqua'  ) then
         kidsat = 784
       else if ( jsatid == 'aura'  ) then
         kidsat = 785
       else if ( jsatid == 'gcom-w1' ) then
         kidsat = 122
! Temporary comment gpm out here; discrepancy between SAID in bufr file and
! kidsat.
!       else if ( jsatid == 'gpm' ) then
!         kidsat = 288
       else if ( jsatid == 'meghat' ) then
         kidsat = 440
       else
         kidsat = 0
       end if

       call closbf(lnbufr)
       open(lnbufr,file=trim(filename),form='unformatted',status ='unknown')
       call openbf(lnbufr,'IN',lnbufr)
       call datelen(10)

       if(kidsat /= 0)then
        lexist = .false.
        satloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
           if(ireadsb(lnbufr)==0)then
              call ufbint(lnbufr,satid,1,1,iret,'SAID')
           end if
           if(nint(satid) == kidsat) then
             lexist=.true.
             exit satloop
           end if
           nread = nread + 1
        end do satloop
       else if(trim(filename) == 'prepbufr')then  ! RTod: wired-in filename is not a good idea
         lexist = .false.
         fileloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
          do while(ireadsb(lnbufr)>=0)
           call ufbint(lnbufr,rtype,1,1,iret,'TYP')
           kx=nint(rtype)
           do nc=1,nconvtype
             if(trim(ioctype(nc)) == trim(dtype) .and. kx == ictype(nc) .and. icuse(nc) > minuse)then
               lexist = .true.
               exit fileloop
             end if
           end do
          end do 
          nread = nread + 1
         end do fileloop
        else if(trim(filename) == 'wcpbufr')then
         lexist = .false.
         file2loop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
          do while(ireadsb(lnbufr)>=0)
           call ufbint(lnbufr,rtype,1,1,iret,'TYP')
           kx=nint(rtype)
           do nc=1,nconvtype
             if(trim(ioctype(nc)) == trim(dtype) .and. kx == ictype(nc) .and. icuse(nc) > minuse)then
               lexist = .true.
               exit file2loop
             end if
           end do
          end do
          nread = nread + 1
         end do file2loop
       else if(trim(filename) == 'gps_ref' .or.  trim(filename) == 'gps_bnd')then
         lexist = .false.
         gpsloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
           if(ireadsb(lnbufr)==0)then 
              call ufbint(lnbufr,satid,1,1,iret,'SAID') 
           end if 
 
           said=nint(satid) 
           if(((said > 739) .and.(said < 746)).or.(said == 820) .or. (said == 825) .or. &
               (said == 786).or. (said == 4)  .or.(said == 3).or. &
               (said == 421).or. (said == 440).or.(said == 821)) then
             lexist=.true. 
             exit gpsloop 
           end if 
           nread = nread + 1
         end do gpsloop
       else if(trim(filename) == 'prepbufr_profl')then  
         lexist = .false.
         airploop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
          do while(ireadsb(lnbufr)>=0)
           call ufbint(lnbufr,rtype,1,1,iret,'TYP')
           kx=nint(rtype)
           if (trim(dtype)=='uv') then
              if (kx==330 .or. kx==430 .or. kx==530) kx=230
              if (kx==331 .or. kx==431 .or. kx==531) kx=231
              if (kx==332 .or. kx==432 .or. kx==532) kx=232
              if (kx==333 .or. kx==433 .or. kx==533) kx=233
              if (kx==334 .or. kx==434 .or. kx==534) kx=234
              if (kx==335 .or. kx==435 .or. kx==535) kx=235
           else
              if (kx==330 .or. kx==430 .or. kx==530) kx=130
              if (kx==331 .or. kx==431 .or. kx==531) kx=131
              if (kx==332 .or. kx==432 .or. kx==532) kx=132
              if (kx==333 .or. kx==433 .or. kx==533) kx=133
              if (kx==334 .or. kx==434 .or. kx==534) kx=134
              if (kx==335 .or. kx==435 .or. kx==535) kx=135
           end if
           do nc=1,nconvtype
             if(trim(ioctype(nc)) == trim(dtype) .and. kx == ictype(nc) .and. icuse(nc) > minuse)then
               lexist = .true.
               exit airploop
             end if
           end do
          end do
          nread = nread + 1
         end do airploop
       else if(trim(filename) == 'satwndbufr')then
         lexist = .false.
         loop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
!        5 GOES-R AMVs (NC005030, NC005031, NC005032, NC005034 and NC005039)
!        are added as the GOES-R bufr file provide do not contain other winds.
!        May not be necessary with the operational satwnd BUFR
            if(trim(subset) == 'NC005010' .or. trim(subset) == 'NC005011' .or.&
               trim(subset) == 'NC005070' .or. trim(subset) == 'NC005071' .or.&
               trim(subset) == 'NC005044' .or. trim(subset) == 'NC005045' .or.&
               trim(subset) == 'NC005046' .or. trim(subset) == 'NC005064' .or.&
               trim(subset) == 'NC005065' .or. trim(subset) == 'NC005066' .or.& 
               trim(subset) == 'NC005030' .or. trim(subset) == 'NC005031' .or.& 
               trim(subset) == 'NC005032' .or. trim(subset) == 'NC005034' .or.&
               trim(subset) == 'NC005039' .or. trim(subset) == 'NC005067' .or.&
               trim(subset) == 'NC005068' .or. trim(subset) == 'NC005069' ) then
               lexist = .true.
               exit loop
            endif
            nread = nread + 1
         end do loop
       else if(trim(filename) == 'oscatbufr')then
         lexist = .false.
         oscatloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
            if(trim(subset) == 'NC012255') then
               lexist = .true.
               exit oscatloop
            endif
         end do oscatloop
       else if(trim(filename) == 'rapidscatbufr')then
         lexist = .false.
         rapidscatloop: do while(ireadmg(lnbufr,subset,idate2) >= 0)
            if(trim(subset) == 'NC012255') then
               lexist = .true.
               exit rapidscatloop
            endif
            nread = nread + 1
         end do rapidscatloop
       else if(trim(filename) == 'hdobbufr')then
         lexist = .false.
         loop_hdob: do while(ireadmg(lnbufr,subset,idate2) >= 0)
            if(trim(subset) == 'NC004015') then
               lexist = .true.
               exit loop_hdob
            endif
            nread = nread + 1
         end do loop_hdob
       else if(trim(dtype) == 'pm2_5')then
          if (oneobtest_chem .and. oneob_type_chem=='pm2_5') then
             lexist=.true.
          else
             lexist = .false.
             fileloopanow_pm2_5:do while(ireadmg(lnbufr,subset,idate2) >= 0)
                do while(ireadsb(lnbufr)>=0)
                   if (subset == 'ANOWPM') then
                      call ufbint(lnbufr,rtype,1,1,iret,'TYP')
                      kx=nint(rtype)
                   else if ( (subset == 'NC008031') .or. &
                          (subset == 'NC008032' ) ) then
                      call ufbint(lnbufr,rtype,1,1,iret,'TYPO')
                      kx=nint(rtype)
                      if (kx/=code_pm25_ncbufr) then
                         cycle
                      else
                         kx=code_pm25_anowbufr
                      endif
                   else
                      cycle
                   endif
                   
                   do nc=1,nconvtype
                      if(trim(ioctype(nc)) == trim(dtype) .and. &
                           kx == ictype(nc) .and. icuse(nc) > minuse)then
                         lexist = .true.
                         exit fileloopanow_pm2_5
                      end if
                   end do
                end do
                nread = nread + 1
             enddo fileloopanow_pm2_5
          endif

          if (lexist) then
             write(6,*)'found pm2_5 in anow bufr'
          else
             write(6,*)'did not find pm2_5 in anow bufr'
          endif
           
       else if(trim(dtype) == 'pm10')then
          lexist = .false.
          fileloopanow_pm10:do while(ireadmg(lnbufr,subset,idate2) >= 0)
             do while(ireadsb(lnbufr)>=0)
                if (subset == 'NC008033') then
                   call ufbint(lnbufr,rtype,1,1,iret,'TYPO')
                   kx=nint(rtype)
                   IF (kx/=code_pm10_ncbufr) then
                      cycle
                   else
                      kx=code_pm10_anowbufr
                   endif
                else
                   cycle
                endif

                do nc=1,nconvtype
                   if(trim(ioctype(nc)) == trim(dtype) .and. &
                        kx == ictype(nc) .and. icuse(nc) > minuse)then
                      lexist = .true.
                      exit fileloopanow_pm10
                   end if
                end do
             end do
             nread = nread + 1
          enddo fileloopanow_pm10

          if (lexist) then
             write(6,*)'found pm10 in anow bufr'
          else
             write(6,*)'did not find pm10 in anow bufr'
          endif


       end if
      end if

      call closbf(lnbufr)
  end if
  if(lexist)then
      write(6,*)'read_obs_check: bufr file date is ',idate,trim(filename),' ',dtype,jsatid
  else
      write(6,*)'read_obs_check: bufr file ',dtype,jsatid,' not available ',trim(filename)
  end if
  return
end subroutine read_obs_check

subroutine read_obs(ndata,mype)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_obs              read, select, reformat obs data
!   prgmmr: parrish          org: np22                date: 1990-10-07
!
! abstract:  This routine is a driver for routines which read different
!            types of observational data.
!
! program history log:
!   1990-10-07  parrish
!   1998-05-15  weiyu yang 
!   1999-08-24  derber, j., treadon, r., yang, w., first frozen mpp version
!   2004-06-16  treadon - update documentation
!   2004-07-23  derber - modify to include conventional sst
!   2004-07-29  treadon - add only to module use, add intent in/out
!   2005-01-20  okamoto - add calling read_ssmi
!   2005-06-14  wu      - add OMI oz
!   2005-07-06  derber - add mhs, hirs/4 and ears data 
!   2005-08-02  derber - modify to use convinfo file
!   2005-09-08  derber - modify to use input group time window
!   2005-09-20  xu & pawlak - modify calling read_ssmis and read_amsre
!   2005-09-28  derber - modify to simplify obs handling      
!   2005-10-17  derber - pass obs_load1 into read_amsre and read_ssmis 
!   2005-10-18  treadon - remove obs_load and obs_load1
!   2005-10-20  kazumori - modify to read real AMSR-E data
!   2005-11-28  derber move determination of which ob data sets to read inside read_obs
!   2005-11-14  li, xu - modify sst obs read and add avhrr gac 1b obs read
!   2006-01-25  treadon - remove read_ieeetovs
!   2006-02-01  parrish - add getsfc and destroy_sfc for full surface fields
!   2006-02-03  derber  - modify for new obs control and obs count
!   2005-02-03  treadon - gather guess 3d pressure to full grid array
!   2006-02-08  derber  - modify to use new convinfo module
!   2006-02-01  liu - add ssu
!   2006-02-24  derber  - modify to take advantage of convinfo module
!   2006-03-07  derber - consolidate processing of 1x1 and 5x5 goessndr
!   2006-04-20  kistler - moved conv_read to gsisub
!   2006-05-25  treadon - rename goesimg and goes_img and pcp_ssm/i as
!                         pcp_ssmi to make consistent with other obstype
!   2006-09-20  treadon - add mpi_io for select data file (obstype)s
!   2006-10-12  treadon - remove tendsflag check for pcp data (now in gsimain)
!   2007-02-21  sienkiewicz - bring in changes for MLS ozone
!   2007-03-15       su - add reading conventional error table option
!   2007-06-05  treadon - restructure mpi_querybf section to improve efficiency
!   2007-10-03  todling - skip most of this in 4dvar inner loop
!   2008-03-28       wu - move random seed for perturb_obs from setuprhsall
!   2008-04-18  safford - rm unused vars and uses
!   2008-05-01    h.liu - add gome ozone
!   2008-06-20   derber - move destroy_sfc to this routine 
!   2008-09-08   lueken - merged ed''s cahnges into q1fy09 code
!   2008-12-30  todling - handle inquire for diff versions of fortran
!   2009-01-05  todling - need tendency alloc in observer mode
!   2009-01-23  todling - echo surface state info 
!   2009-03-18  meunier - add a if statement to read lagrangian data
!   2009-08-19  guo     - moved destroy_sfc_grid() to observer_finalize().
!   2009-12-20  gayno - modify argument lists so that fov-based surface
!                       calculation may be used.
!   2010-03-29  hu    - add code to read in cloud observations  including:
!                            prepbufr (metar, nesdis cloud product)
!                            radar reflectivity, lightning, NASA LaRC cloud
!   2010-04-01  treadon - move strip and reorder to gridmod
!   2010-04-08  hliu - add seviri
!
!   2010-04-05  huang   - add aero and modis for reading modis aod from satellite
!                         currently read BUFR file only
!   2010-04-22  tangborn - read for carbon monoxide
!   2010-08-23  tong    - add calcuation of hgtl_full used in read_radar.f90
!   2011-04-02  li       - add nst_gsi, getnst and destroy_nst
!   2011-05-20  mccarty  - add cris/atms handling
!   2011-05-26  todling  - add call to create_nst
!   2012-03-07  akella   - set NST variables; init & destroy in guess_grids 
!                          [see create_sfc_grids & destroy_sfc_grids] 
!   2012-05-16  todling  - protect call to prebufr when NST is on (should be mods)
!   2012-09-08  wargan   - add OMI with efficiency factors
!   2013-01-26  parrish - WCOSS debug compile fails--extra arguments in call read_aerosol.
!                         Commented out extra line of arguments not used.
!   2013-02-13  eliu     - turn off parallel I/O for SSMIS (due to the need to
!                          do spatial averaging for noise reduction) 
!   2013-06-01  zhu     - add mype_airobst to handle aircraft temperature bias correction 
!   2013-08-08  s.liu     - add read NASA_LaRC_cloud product
!   2013-10-25  todling - reposition ltosi and others to commvars
!   2014-01-01  xli     - add option to read NSST marine BUFR data file nsstbufr (on the
!                         top of prepbufr and modsbufr)
!   2014-02-03  guo     - Hid processing (read) of non-EMC ozone obstypes
!                         through module m_extOzone, separated from read_ozone.
!                       - Added some -do- and -if- construct names, for easier
!                         understanding of the code.
!   2014-06-19  carley/zhu - Add tcamt and lcbas
!   2014-11-12  carley  - Add call to read goes imager sky cover data for tcamt
!   2014-12-03  derber - modify for no radiance cases and read processor for
!                        surface fields
!   2015-01-16  ejones  - added saphir, gmi, and amsr2 handling
!   2015-03-23  zaizhong ma - add Himawari-8 ahi
!   2015-05-30  li     - modify for no radiance cases but sst (nsstbufr) and read processor for
!                        surface fields (use_sfc = .true. for data type of sst),
!                        to use deter_sfc in read_nsstbufr.f90)
!   2015-07-10  pondeca - add cloud ceiling height (cldch)
!   2015-08-12  pondeca - add capability to read min/maxT obs from ascii file
!   2015-08-20  zhu     - use centralized radiance info from radiance_mod to specify
!                         all-sky and aerosol usages in radiance assimilation
!   2015-09-04  J. Jung - Added mods for CrIS full spectral resolution (FSR)
!   2015-10-19  lippi   - Added logic to read l2rw or rw radial winds.
!   2016-03-02  s.liu/carley - remove use_reflectivity and use i_gsdcldanal_type
!   2016-04-28  J. Jung - added logic for RARS and direct broadcast data from NESDIS/UW.
!   2016-05-05  pondeca - add 10-m u-wind and v-wind (uwnd10m, vwnd10m)
!   2016-09-19  Guo     - replaced open(obs_input_common) with "call unformatted_open(obs_input_common)"
!   2017-05-12  Y. Wang and X. Wang - add multi-interface to read in dBZ (nc) and radial velocity (ascii)
!   2016-12-14  lippi   - Fixed bug of using observations twice when both
!                         l2rwbufr and radarbufr are in the OBS_INPUT table.
!                         Changed the dsis entries for l2rwbufr and radarbufr to
!                         l2rw and l3rw respectively. Also make use of nml
!                         option vadwnd_l2rw_qc. 
!   2017-08-31  Li      - move gsi_nstcoupler_init & gsi_nstcoupler_read to getsfc in sathin.F90
!                       - move gsi_nstcoupler_final from create_sfc_grids to here
!   2018-01-23 Apodaca  - add GOES/GLM lightning data
!   2019-01-15  Li      - add to handle mbuoyb
!   2019-03-27  h. liu   - add abi
!   
!
!   input argument list:
!     mype     - mpi task id
!
!   output argument list:
!     ndata(*,1)- number of prefiles retained for further processing
!     ndata(*,2)- number of observations read
!     ndata(*,3)- number of observations keep after read
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$  end documentation block
    use kinds, only: r_kind,i_kind,i_llong,r_double
    use gridmod, only: nlon,nlat,nsig,iglobal,ijn,itotsub,lat1,lon1,&
         displs_g,strip,reorder
    use general_commvars_mod, only: ltosi,ltosj
    use obsmod, only: iadate,ndat,time_window,dplat,dsfcalc,dfile,dthin, &
           dtype,dval,dmesh,obsfile_all,ref_obs,nprof_gps,dsis,ditype,&
           perturb_obs,lobserver,lread_obs_save,obs_input_common, &
           reduce_diag,nobs_sub,dval_use
    use gsi_nstcouplermod, only: nst_gsi,gsi_nstcoupler_final
    use qcmod, only: njqc,vadwnd_l2rw_qc
    use gsi_4dvar, only: l4dvar
    use satthin, only: super_val,super_val1,superp,makegvals,getsfc,destroy_sfc
    use mpimod, only: ierror,mpi_comm_world,mpi_sum,mpi_rtype,mpi_integer,npe,&
         setcomm
    use constants, only: one,zero
    use converr, only: converr_read
    use converr_ps, only: converr_ps_read
    use converr_q, only: converr_q_read
    use converr_t, only: converr_t_read
    use converr_uv, only: converr_uv_read
    use converr_pw, only: converr_pw_read
    use convb_ps,only: convb_ps_read
    use convb_q,only:convb_q_read
    use convb_t,only:convb_t_read
    use convb_uv,only:convb_uv_read
    use guess_grids, only: ges_prsl,geop_hgtl,ntguessig
    use radinfo, only: nusis,iuse_rad,jpch_rad,diag_rad
    use insitu_info, only: mbuoy_info,mbuoyb_info,read_ship_info
    use aeroinfo, only: nusis_aero,iuse_aero,jpch_aero,diag_aero
    use ozinfo, only: nusis_oz,iuse_oz,jpch_oz,diag_ozone
    use pcpinfo, only: npcptype,nupcp,iusep,diag_pcp
    use convinfo, only: nconvtype,ioctype,icuse,diag_conv,ithin_conv
    use chemmod, only : oneobtest_chem,oneob_type_chem,oneobschem
    use lightinfo, only: nlighttype,iuse_light,diag_light
    use aircraftinfo, only: aircraft_t_bc,aircraft_t_bc_pof,aircraft_t_bc_ext,mype_airobst
    use gsi_io, only: mype_io
    use rapidrefresh_cldsurf_mod, only: i_gsdcldanal_type
    use radiance_mod, only: rad_obs_type,radiance_obstype_search

    use m_extOzone, only: is_extOzone
    use m_extOzone, only: extOzone_read
    use mpeu_util, only: warn
    use gsi_unformatted, only: unformatted_open

    use mrmsmod,only: l_mrms_sparse_netcdf

    implicit none

!   Declare passed variables
    integer(i_kind)                  ,intent(in   ) :: mype
    integer(i_kind),dimension(ndat,3),intent(  out) :: ndata

!   Declare local parameters
    integer(i_llong),parameter:: lenbuf=8388608_i_llong  ! lenbuf=8*1024*1024

!   Declare local variables
    logical :: lexist,ssmis,amsre,sndr,hirs,avhrr,lexistears,lexistdb,use_prsl_full,use_hgtl_full
    logical :: use_sfc,nuse,use_prsl_full_proc,use_hgtl_full_proc,seviri,mls,abi
    logical,dimension(ndat):: belong,parallel_read,ears_possible,db_possible
    logical :: modis,use_sfc_any
    logical :: acft_profl_file
    character(10):: obstype,platid
    character(22):: string
    character(120):: infile
    character(20):: sis
    integer(i_kind) i,j,k,ii,nmind,lunout,isfcalc,ithinx,ithin,nread,npuse,nouse
    integer(i_kind) nprof_gps1,npem1,krsize,len4file,npemax,ilarge,nlarge,npestart
    integer(i_llong) :: lenbytes
    integer(i_kind):: npetot,npeextra,mmdat,nodata
    integer(i_kind):: iworld,iworld_group,next_mype,mm1,iix
    integer(i_kind):: mype_root
    integer(i_kind):: minuse,lunsave,maxproc,minproc
    integer(i_kind),dimension(ndat):: npe_sub,npe_sub3,mpi_comm_sub,mype_root_sub,npe_order
    integer(i_kind),dimension(ndat):: ntasks1,ntasks
    integer(i_kind),dimension(ndat):: read_rec1,read_rec
    integer(i_kind),dimension(ndat):: read_ears_rec1,read_ears_rec
    integer(i_kind),dimension(ndat):: read_db_rec1,read_db_rec
    integer(i_kind),dimension(ndat,3):: ndata1
    integer(i_kind),dimension(npe,ndat):: mype_work,nobs_sub1
    integer(i_kind),dimension(npe,ndat):: mype_sub
    integer(i_kind),allocatable,dimension(:):: nrnd
    integer(i_kind):: nmls_type,mype_io_sfc
    integer(i_kind):: iread,ipuse,iouse

    real(r_kind) gstime,val_dat,rmesh,twind,rseed
    real(r_kind),allocatable,dimension(:) :: prslsm,hgtlsm,work1
    real(r_kind),allocatable,dimension(:,:,:):: prsl_full,hgtl_full

    type(rad_obs_type) :: radmod

    data lunout / 81 /
    data lunsave  / 82 /

!*****************************************************************************

!   Set analysis time and allocate/initialize arrays and variables
    call w3fs21(iadate,nmind)
    gstime=real(nmind,r_kind)
    allocate(nobs_sub(npe,ndat))
    nobs_sub = 0
    nobs_sub1 = 0

    call makegvals
    do ii=1,ndat
       ndata1(ii,1)=0
       ndata1(ii,2)=0
       ndata1(ii,3)=0
       ntasks1(ii) =0
       parallel_read=.false.
    end do
    npem1=npe-1
    nprof_gps1=0

    if(njqc) then
       call converr_ps_read(mype)
       call converr_q_read(mype)
       call converr_t_read(mype)
       call converr_uv_read(mype)
       call converr_pw_read(mype)
       call convb_ps_read(mype)
       call convb_q_read(mype)
       call convb_t_read(mype)
       call convb_uv_read(mype)
    else
       call converr_read(mype)
    endif

!   Optionally set random seed to perturb observations
    if (perturb_obs) then
       rseed=iadate(4)+iadate(3)*100+iadate(2)*10000+iadate(1)*1000000+mype
       call random_seed(size=krsize)
       allocate(nrnd(krsize))
       do i=1,krsize
          nrnd(i)=rseed
       end do
       call random_seed(put=nrnd)
       deallocate(nrnd)
    endif



!   Set data class and number of reader tasks.  Set logical flag to indicate 
!   type type of GPS data (if present)
    ii=0
    ref_obs = .false.    !.false. = assimilate GPS bending angle
    ears_possible = .false.
    db_possible = .false.
    nmls_type=0
    read_rec1 = 0
    read_ears_rec1=0
    read_db_rec1=0
    do i=1,ndat
       obstype=dtype(i)                   !     obstype  - observation types to process
       amsre= index(obstype,'amsre') /= 0
       ssmis= index(obstype,'ssmis') /= 0
       sndr = index(obstype,'sndr') /= 0
       hirs = index(obstype,'hirs') /= 0
       avhrr = index(obstype,'avhrr') /= 0
       modis = index(obstype,'modis') /= 0
       seviri = index(obstype,'seviri') /= 0
       abi = index(obstype,'abi') /= 0
       mls = index(obstype,'mls') /= 0
       if(obstype == 'mls20' ) nmls_type=nmls_type+1
       if(obstype == 'mls22' ) nmls_type=nmls_type+1
       if(obstype == 'mls30' ) nmls_type=nmls_type+1
       if(nmls_type>1) then
          write(6,*) '******ERROR***********: there is more than one MLS data type, not allowed, please check'
          call stop2(339)
       end if
       if (obstype == 't'  .or. obstype == 'uv' .or. &
           obstype == 'q'  .or. obstype == 'ps' .or. &
           obstype == 'pw' .or. obstype == 'spd'.or. &
           obstype == 'sst'.or. &
           obstype == 'tcp'.or. obstype == "lag".or. &
           obstype == 'dw' .or. obstype == 'rw' .or. &
           obstype == 'mta_cld' .or. obstype == 'gos_ctp' .or. &
           obstype == 'rad_ref' .or. obstype=='lghtn' .or. &
           obstype == 'larccld' .or. obstype == 'pm2_5' .or. obstype == 'pm10' .or. &
           obstype == 'gust' .or. obstype=='vis' .or. &
           obstype == 'pblh' .or. obstype=='wspd10m' .or. &
           obstype == 'td2m' .or. obstype=='mxtm' .or. &
           obstype == 'mitm' .or. obstype=='pmsl' .or. &
           obstype == 'howv' .or. obstype=='tcamt' .or. &
           obstype=='lcbas' .or. obstype=='cldch' .or. obstype == 'larcglb' .or. &
           obstype=='uwnd10m' .or. obstype=='vwnd10m' .or. obstype=='dbz' ) then
          ditype(i) = 'conv'
       else if (obstype == 'swcp' .or. obstype == 'lwcp') then
          ditype(i) = 'wcp'
       else if( hirs   .or. sndr      .or.  seviri .or. abi .or.        &
               obstype == 'airs'      .or. obstype == 'amsua'     .or.  &
               obstype == 'msu'       .or. obstype == 'iasi'      .or.  &
               obstype == 'amsub'     .or. obstype == 'mhs'       .or.  &
               obstype == 'hsb'       .or. obstype == 'goes_img'  .or.  &
               obstype == 'ahi'       .or. avhrr                  .or.  &
               amsre  .or. ssmis      .or. obstype == 'ssmi'      .or.  &
               obstype == 'ssu'       .or. obstype == 'atms'      .or.  &
               obstype == 'cris'      .or. obstype == 'cris-fsr'  .or.  &
               obstype == 'amsr2'     .or.  &
               obstype == 'gmi'       .or. obstype == 'saphir'   ) then
          ditype(i) = 'rad'
       else if (is_extOzone(dfile(i),obstype,dplat(i))) then
          ditype(i) = 'ozone'
       else if (obstype == 'sbuv2' &
           .or. obstype == 'omi' &
           .or. obstype == 'ompstc8' &
           .or. obstype == 'ompsnp' &
           .or. obstype == 'gome' &
           .or. mls &
           ) then
          ditype(i) = 'ozone'
       else if (obstype == 'mopitt') then
          ditype(i) = 'co'
       else if (index(obstype,'pcp')/=0 )then
          ditype(i) = 'pcp'
       else if (obstype == 'gps_ref' .or. obstype == 'gps_bnd') then
          ditype(i) = 'gps'
       else if ( index(obstype,'aod') /= 0 ) then
          ditype(i) = 'aero'
       else if (obstype == 'goes_glm') then
          ditype(i) = 'light'
       else
          write(6,*)'READ_OBS:  ***ERROR*** - unknown ob type ',trim(obstype)
       end if

!   Set data class and number of reader tasks.  Set logical flag to indicate 
!   type type of GPS data (if present)
       if (index(dtype(i),'gps_ref') /= 0) ref_obs = .true.

!   Check info files to see if data is read.

       nuse=.false.
       minuse=-1
       if(ditype(i) == 'conv')then
          if(diag_conv .and. .not. reduce_diag)minuse=-2
          do j=1,nconvtype
             if(trim(dtype(i)) == trim(ioctype(j)) .and. icuse(j) > minuse)nuse=.true.
          end do
       else if(ditype(i) == 'rad')then
          if(diag_rad .and. .not. reduce_diag)minuse=-2
          do j=1,jpch_rad
             if(trim(dsis(i)) == trim(nusis(j)) .and. iuse_rad(j) > minuse)nuse=.true.
          end do
       else if(ditype(i) == 'ozone')then
          if(diag_ozone .and. .not. reduce_diag)minuse=-2
          do j=1,jpch_oz
             if(trim(dsis(i)) == trim(nusis_oz(j)) .and. iuse_oz(j) > minuse)nuse=.true.
          end do
       else if(ditype(i) == 'pcp')then
          if(diag_pcp .and. .not. reduce_diag)minuse=-2
          do j=1,npcptype
             if(trim(dsis(i)) == trim(nupcp(j)) .and. iusep(j) > minuse)nuse=.true.
          end do
       else if(ditype(i) == 'aero')then
          if(diag_aero .and. .not. reduce_diag)minuse=-2
          do j=1,jpch_aero
             if(trim(dsis(i)) == trim(nusis_aero(j)) .and. iuse_aero(j) > minuse)nuse=.true.
          end do
        else if(ditype(i) == 'light')then
          if(diag_light .and. .not. reduce_diag)minuse=-2
          do j=1,nlighttype
             if(iuse_light(j) > minuse) nuse=.true.
          end do
       else
          nuse=.true.
       end if

       if(nuse)then

!     Control parallel read for each ob type (currently just rad obs).  
!     To remove parallel read comment out line.
          ithin=dthin(i)
          if(ithin > 0 )then
            if(dmesh(ithin) > one)then
             if(hirs)then
                parallel_read(i)= .true.
             else if(obstype == 'amsua')then
                parallel_read(i)= .true.
             else if(obstype == 'airs' )then
                parallel_read(i)= .true.
             else if(obstype == 'iasi')then
                parallel_read(i)= .true.
             else if(obstype == 'amsub')then
                parallel_read(i)= .true.
             else if(obstype == 'mhs' )then
                parallel_read(i)= .true.
             else if(sndr )then
                parallel_read(i)= .true.
! N.B. ATMS must be run on one processor for the filtering code to work.
             else if(obstype == 'atms')then
!                 parallel_read(i)= .true.
             else if(ssmis)then
!               parallel_read(i)= .true.  
             else if(seviri)then
                parallel_read(i)= .true.
             else if(abi)then
                parallel_read(i)= .true.
             else if(obstype == 'cris' .or. obstype == 'cris-fsr')then
                parallel_read(i)= .true.
             else if(avhrr)then
                parallel_read(i)= .true.
             else if(amsre)then
                parallel_read(i)= .true.
             else if(obstype == 'goes_img' )then
                parallel_read(i)= .true.
             else if(obstype == 'ahi' )then
                parallel_read(i)= .true.
             else if(obstype == 'hsb' )then
                parallel_read(i)= .true.
             else if(obstype == 'ssmi' )then
                parallel_read(i)= .true.
             else if(obstype == 'ssu' )then
                parallel_read(i)= .true.
             else if(obstype == 'amsr2')then    
!                parallel_read(i)= .true.     ! turn parallel read off for spatial averaging
             else if(obstype == 'gmi')then
!                parallel_read(i)= .true.     ! turn parallel read off for spatial averaging
!   Parallel read for SAPHIR not currently working. Leave parallel read off.
!             else if(obstype == 'saphir')then
!                parallel_read(i)= .true.

             end if
           end if
          end if
! direct broadcast from EUMETSAT (EARS)
          ears_possible(i) = ditype(i) == 'rad'  .and. & 
                  (obstype == 'amsua' .or. obstype == 'amsub' .or. & 
                   obstype == 'mhs'   .or. obstype == 'hirs3' .or. &
                   obstype == 'cris'  .or. obstype == 'cris-fsr' .or. &
                   obstype == 'iasi'  .or. obstype == 'atms') .and. &
                  (dplat(i) == 'n17' .or. dplat(i) == 'n18' .or. & 
                   dplat(i) == 'n19' .or. dplat(i) == 'npp' .or. &
                   dplat(i) == 'n20' .or. &
                   dplat(i) == 'metop-a' .or. dplat(i) == 'metop-b' .or. &
                   dplat(i) == 'metop-c') 
! direct broadcast from NESDIS/UW
          db_possible(i) = ditype(i) == 'rad'  .and.       & 
                  (obstype == 'amsua' .or.  obstype == 'amsub' .or.  & 
                   obstype == 'mhs' .or. obstype == 'atms' .or. &
                   obstype == 'cris' .or. obstype == 'cris-fsr' .or. &
                   obstype == 'iasi') .and. &
                  (dplat(i) == 'n17' .or. dplat(i) == 'n18' .or. & 
                   dplat(i) == 'n19' .or. dplat(i) == 'npp' .or. &
                   dplat(i) == 'n20' .or. &
                   dplat(i) == 'metop-a' .or. dplat(i) == 'metop-b' .or. &
                   dplat(i) == 'metop-c') 

!   Inquire data set to deterimine if input data available and size of dataset
          ii=ii+1
          if (ii>npem1) ii=0
          if(mype==ii)then
             call gsi_inquire(lenbytes,lexist,trim(dfile(i)),mype)
             call read_obs_check (lexist,trim(dfile(i)),dplat(i),dtype(i),minuse,read_rec1(i))
             
!   If no data set starting record to be 999999.  Note if this is not large
!   enough code should still work - just does a bit more work.

             if(.not. lexist)read_rec1(i) = 999999
             len4file=lenbytes/4
             if (ears_possible(i))then

                call gsi_inquire(lenbytes,lexistears,trim(dfile(i))//'ears',mype)
                call read_obs_check (lexistears,trim(dfile(i))//'ears',dplat(i),dtype(i),minuse, &
                    read_ears_rec1(i))

!   If no data set starting record to be 999999.  Note if this is not large
!   enough code should still work - just does a bit more work.

                if(.not. lexistears)read_ears_rec1(i) = 999999
                lexist=lexist .or. lexistears
                len4file=len4file+lenbytes/4
             else
               read_ears_rec1(i) = 999999
             end if
             if (db_possible(i))then

                call gsi_inquire(lenbytes,lexistdb,trim(dfile(i))//'_db',mype)
                call read_obs_check (lexistdb,trim(dfile(i))//'_db',dplat(i),dtype(i),minuse, &
                    read_db_rec1(i))

!   If no data set starting record to be 999999.  Note if this is not large
!   enough code should still work - just does a bit more work.

                if(.not. lexistdb)read_db_rec1(i) = 999999
                lexist=lexist .or. lexistdb
                len4file=len4file+lenbytes/4
             else
               read_db_rec1(i) = 999999
             end if

 
             if(lexist) then
!      Initialize number of reader tasks to 1.  For the time being
!      only allow number of reader tasks >= 1 for select obstype.

                ntasks1(i)=1
                if(parallel_read(i)) then

!  Allow up to 16 processors/file increase loop bounds to increase number of processors allowed
                   do j=1,4
                      if(len4file < lenbuf)exit
                      ntasks1(i)=2*ntasks1(i)
                      len4file=len4file/2
                   end do
!               if(ntasks1(i)*lenbuf < len4file) ntasks1(i)=ntasks1(i)+1
                end if
             end if
          end if
       else
          if(mype == 0)write(6,*) 'data type ',dsis(i), &
                'not used in info file -- do not read file ',trim(dfile(i))
       end if
    end do


!   Distribute optimal number of reader tasks to all mpi tasks
    call mpi_allreduce(ntasks1,ntasks,ndat,mpi_integer,mpi_sum,mpi_comm_world,ierror)
    call mpi_allreduce(read_rec1,read_rec,ndat,mpi_integer,mpi_sum,mpi_comm_world,ierror) 
    call mpi_allreduce(read_ears_rec1,read_ears_rec,ndat,mpi_integer,mpi_sum,mpi_comm_world,ierror) 
    call mpi_allreduce(read_db_rec1,read_db_rec,ndat,mpi_integer,mpi_sum,mpi_comm_world,ierror) 

!   Limit number of requested tasks per type to be <= total available tasks
    npemax=0
    npetot=0
    do i=1,ndat
       if (ntasks(i)>npe) then
          write(6,*)'read_obs:  ***WARNING*** i=',i,' dtype=',dtype(i),' dsis=',dsis(i),&
               ' requested ntasks=',ntasks(i),' > npe=',npe,' reset ntasks=',npe
          ntasks(i)=npe
       endif
       npe_sub(i)=ntasks(i)
       npetot=npetot+npe_sub(i)
       npemax=max(npemax,npe_sub(i))
    end do

    if(l4dvar.and.(.not.lobserver)) then
!_RTod use_sfc=.false.
!_RTod call getsfc(mype,use_sfc)
       return
    endif
    
    npeextra=0
    if(mod(npetot,npe) > 0) npeextra=npe-mod(npetot,npe)
    maxproc=32
    if(npeextra > 0)then
       if(mype == 0)write(6,*) ' number of extra processors ',npeextra
       extraloop: do j=1,npeextra
          npe_sub3=ntasks
          minproc=999999
          do i=1,ndat
            if(ntasks(i) > 0 .and.  parallel_read(i))minproc=min(minproc,ntasks(i))
          end do
          if(npeextra < minproc) exit extraloop
          iix=minproc
          do ii=1,6
             do i=1,ndat
                if(iix == npe_sub3(i) .and. parallel_read(i))then
                   if(ntasks(i) <= npeextra .and. ntasks(i) < maxproc)then
                      npeextra=npeextra-ntasks(i)
                      ntasks(i)=2*ntasks(i)
                      if(npeextra < iix)cycle extraloop
                   end if
                end if
             end do
             iix=2*iix
             if(iix >= maxproc) cycle extraloop
          end do
       end do extraloop
    end if

!   Set up locations of first processor

    ilarge=0
    npestart=0
    npe_sub3=npe_sub
    mype_root_sub=0
    mmdat=0
    loopx: do j=1,ndat
       nlarge=0
       do i=1,ndat
          if(npe_sub3(i) > nlarge .and. npe_sub3(i)+npestart <= npe)then
             ilarge=i
             nlarge=npe_sub3(i)
          end if
       end do
       if(nlarge == 0)exit loopx
       npe_order(j)=ilarge
       mype_root_sub(ilarge)=npestart
       npestart=npestart+npe_sub3(ilarge)
       mmdat=mmdat+1
       npe_sub3(ilarge)=0
       if(npestart + minval(npe_sub3, mask=npe_sub3>0)>= npe) npestart=0
    end do loopx

!   Define sub-communicators for each data file
    mm1=mype+1
    belong=.false.
    mype_sub=-999
    mype_root=0
    next_mype=0
    do ii=1,mmdat
       i=npe_order(ii)
       if(npe_sub(i) > 0)then
          next_mype=mype_root_sub(i)
          do k=1,npe_sub(i)
             mype_work(k,i) = next_mype
             mype_sub(mype_work(k,i)+1,i)=k-1
             if(next_mype == mype)belong(i) = .true.
             next_mype = next_mype + 1
             if (next_mype>npem1) next_mype=0
          end do               

          call setcomm(iworld,iworld_group,npe_sub(i),mype_work(1,i),&
                 mpi_comm_sub(i),ierror)
       end if

    end do
    mype_airobst = mype_root
    do ii=1,mmdat
       i=npe_order(ii)
       if(mype == 0 .and. npe_sub(i) > 0) write(6,'(1x,a,i4,1x,a,1x,2a,2i4,1x,i6,1x,i6,1x,i6)') &
        'READ_OBS:  read ',i,dtype(i),dsis(i),' using ntasks=',ntasks(i),mype_root_sub(i), & 
               read_rec(i),read_ears_rec(i),read_db_rec(i)

       acft_profl_file = index(dfile(i),'_profl')/=0
       if ((aircraft_t_bc_pof .or. aircraft_t_bc_ext .or. &
          (aircraft_t_bc .and. acft_profl_file)) .and. dtype(i) == 't') &
                   mype_airobst = mype_root_sub(i)

    end do


    use_prsl_full=.false.
    use_hgtl_full=.false.
    use_sfc=.false.
    use_prsl_full_proc=.false.
    use_hgtl_full_proc=.false.
    mype_io_sfc=mype_io
    do ii=1,mmdat
       i=npe_order(ii)
       if(ditype(i) =='conv')then
          obstype=dtype(i)
          if (obstype == 't' .or. obstype == 'q'  .or. &
              obstype == 'uv' .or. obstype == 'wspd10m' .or. &
              obstype == 'uwnd10m' .or. obstype == 'vwnd10m') then
              use_prsl_full=.true.
              if(belong(i))use_prsl_full_proc=.true.
          else
            do j=1,nconvtype
               if(obstype == trim(ioctype(j)) .and. ithin_conv(j) > 0)then
                  use_prsl_full=.true.
                  if(belong(i))use_prsl_full_proc=.true.
               end if
            end do
          end if
          if(obstype == 'rw')then
             use_hgtl_full=.true.
             if(belong(i))use_hgtl_full_proc=.true.
           else if(obstype == 'dbz')then
             use_hgtl_full=.true.
             if(belong(i))use_hgtl_full_proc=.true.
          end if
          if(obstype == 'sst')then
            if(belong(i))use_sfc=.true.
          endif
       else if(ditype(i) == 'rad' )then
          if(belong(i)) use_sfc=.true.
       end if
    end do
    use_sfc_any=.false.
    loop: do ii=1,mmdat
       i=npe_order(ii)
       if(ditype(i) == 'rad' .or. ditype(i) == 'sst')then
          mype_io_sfc=mype_root_sub(i)
          use_sfc_any=.true.
          exit loop
       end if
    end do loop
    if(use_sfc_any .and. mype == mype_io)use_sfc=.true.

!   Create full horizontal surface fields from local fields in guess_grids
    call getsfc(mype,mype_io_sfc,use_sfc,use_sfc_any)

    if(mype == mype_io) call prt_guessfc2('sfcges2',use_sfc)

!   Get guess 3d pressure on full grid
    allocate(work1(max(iglobal,itotsub)),prslsm(lat1*lon1))
    if(use_prsl_full)then
       if(use_prsl_full_proc)then
          allocate(prsl_full(nlat,nlon,nsig))
       else
          allocate(prsl_full(1,1,1))
       end if
       do k=1,nsig
          call strip(ges_prsl(:,:,k,ntguessig),prslsm)
          call mpi_allgatherv(prslsm,ijn(mype+1),mpi_rtype,&
               work1,ijn,displs_g,mpi_rtype,mpi_comm_world,ierror)
          if(use_prsl_full_proc)then
             call reorder(work1,1,1)
             do ii=1,iglobal
                i=ltosi(ii)
                j=ltosj(ii)
                prsl_full(i,j,k)=work1(ii)
             end do
          end if
       end do
    else
       allocate(prsl_full(1,1,1))
    end if
!   Get guess 3d geopotential height on full grid
    if(use_hgtl_full)then
       allocate(hgtlsm(lat1*lon1))
       if(use_hgtl_full_proc)then
          allocate(hgtl_full(nlat,nlon,nsig))
       else
          allocate(hgtl_full(1,1,1))
       end if
       do k=1,nsig
          call strip(geop_hgtl(:,:,k,ntguessig),hgtlsm)
          call mpi_allgatherv(hgtlsm,ijn(mype+1),mpi_rtype,&
               work1,ijn,displs_g,mpi_rtype,mpi_comm_world,ierror)
          if(use_hgtl_full_proc)then
             call reorder(work1,1,1)
             do ii=1,iglobal
                i=ltosi(ii)
                j=ltosj(ii)
                hgtl_full(i,j,k)=work1(ii)
             end do
           end if
       end do
       deallocate(hgtlsm)
    else
      allocate(hgtl_full(1,1,1))
    end if
    deallocate(work1,prslsm)

!   Create moored buoy station ID
    call mbuoy_info(mype)

!   Create moored buoy station ID for mbuoyb with 7-digit station ID
    call mbuoyb_info(mype)

!   Create ships info(ID, Depth & Instrument)
    call read_ship_info(mype)

!   Loop over data files.  Each data file is read by a sub-communicator
    loop_of_obsdata_files: &
    do ii=1,mmdat

       i=npe_order(ii)
       task_belongs: &
       if (i > 0 .and. belong(i)) then

          platid=dplat(i)                    !     platid   - satellites to read
          obstype=dtype(i)                   !     obstype  - observation types to process
          infile=trim(dfile(i))              !     infile   - units from which to read data
          sis=dsis(i)                        !     sensor/instrument/satellite indicator
          val_dat=dval(i)                    !     weighting factors applied to super obs
          ithin=dthin(i)                     !     ithin    - flags to thin data
          ithinx=max(1,abs(ithin))
          rmesh=dmesh(ithinx)                !     rmesh    - thinning mesh sizes (km)
          twind=time_window(i)               !     time window (hours) for input group
          isfcalc=dsfcalc(i)                 !     method to calculate surface fields within fov
          nread=0
          nouse=0
          npuse=0

          if (mype_sub(mm1,i)==mype_root) then
             open(lunout,file=obsfile_all(i),form='unformatted')
             rewind(lunout)
          endif

!         Process conventional (prepbufr) data
          string="--"//trim(ditype(i))
          ditype_select: &
          if(ditype(i) == 'conv')then
             conv_obstype_select: &
             if (obstype == 't' .or. obstype == 'q'  .or. obstype == 'ps' .or. &
                 obstype == 'pw' .or. obstype == 'spd'.or. & 
                 obstype == 'gust' .or. obstype == 'vis'.or. &
                 obstype == 'td2m' .or. &
!                 obstype=='mxtm' .or. obstype == 'mitm' .or. &
!                 obstype=='howv' .or. obstype=='pmsl' .or. &    #howv #ww3 
                 obstype=='pmsl' .or. &
                 obstype == 'mta_cld' .or. obstype == 'gos_ctp' .or. &
                 obstype == 'lcbas' .or. obstype == 'cldch' ) then

!               Process flight-letel high-density data not included in prepbufr
                if ( index(infile,'hdobbufr') /=0 ) then
                  call read_fl_hdob(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,&
                                    prsl_full,nobs_sub1(1,i))
                  string='READ_FL_HDOB'
                else
                   call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
                        prsl_full,nobs_sub1(1,i),read_rec(i))
                   string='READ_PREPBUFR'

                endif
             else if(obstype == 'howv') then
                 if ( index(infile,'satmar') /=0) then
                   
                    call read_satmar (nread, npuse, nouse                                  &
                                     , infile, obstype, lunout, gstime, twind, sis         &
                                     , nobs_sub1(1,i))
                    string='READ_SATMAR'
                 else 
                   call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
                        prsl_full,nobs_sub1(1,i),read_rec(i))
                   string='READ_PREPBUFR'

                 endif
             else if(obstype == 'mitm') then
                if ( index(infile,'mitmdat') /=0) then
                   call read_mitm_mxtm(nread,npuse,nouse,infile,obstype,lunout,gstime,sis, & 
                                       nobs_sub1(1,i))
                   string='READ_ASCII_MITM'
                 else
                   call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
                        prsl_full,nobs_sub1(1,i),read_rec(i))
                   string='READ_PREPBUFR'
                 endif

             else if(obstype == 'mxtm') then
                if ( index(infile,'mxtmdat') /=0) then
                   call read_mitm_mxtm(nread,npuse,nouse,infile,obstype,lunout,gstime,sis, & 
                                       nobs_sub1(1,i))
                   string='READ_ASCII_MXTM'
                 else
                   call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
                        prsl_full,nobs_sub1(1,i),read_rec(i))
                   string='READ_PREPBUFR'
                 endif

!            Process total cloud amount (tcamt) in prepbufr -or- from goes imager sky cover products
             else if(obstype == 'tcamt') then
!             Process GOES Imager Sky Cover product separately from prepbufr-based sky cover obs
                if ( index(infile,'goessky') /=0 ) then
                   call read_goesimgr_skycover(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,&
                        prsl_full,nobs_sub1(1,i))
                   string='READ_GOESIMGR_SKYCOVER'
                else
!              else read from prepbufr
                   call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,prsl_full, &
                         nobs_sub1(1,i),read_rec(i))
                   string='READ_PREPBUFR'
                end if

!             Process winds in the prepbufr
            else if(obstype == 'uv' .or. obstype == 'wspd10m' .or. &
                    obstype == 'uwnd10m' .or. obstype == 'vwnd10m') then
!             Process satellite winds which seperate from prepbufr
                if ( index(infile,'satwnd') /=0 ) then
                  call read_satwnd(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,&
                     prsl_full,nobs_sub1(1,i))
                  string='READ_SATWND'
!             Process oscat winds which seperate from prepbufr
                elseif ( index(infile,'oscatbufr') /=0 ) then
                  call read_sfcwnd(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,&
                     prsl_full,nobs_sub1(1,i))
                  string='READ_SFCWND'
!             Process rapidscat winds which seperate from prepbufr
                elseif ( index(infile,'rapidscatbufr') /=0 ) then
                  call read_rapidscat(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,&
                     prsl_full,nobs_sub1(1,i))
                  string='READ_RAPIDSCAT'
                else if ( index(infile,'hdobbufr') /=0 ) then
                  call read_fl_hdob(nread,npuse,nouse,infile,obstype,lunout,gstime,twind,sis,&
                       prsl_full,nobs_sub1(1,i))
                  string='READ_FL_HDOB'
                else
                  call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
                     prsl_full,nobs_sub1(1,i),read_rec(i))
                  string='READ_PREPBUFR'
                endif

!            Process conventional SST (nsstbufr, at this moment) data
             elseif ( obstype == 'sst' ) then
                string="--"//trim(ditype(i))//":sst:"//trim(platid)
                if ( platid == 'nsst') then
                   call read_nsstbufr(nread,npuse,nouse,gstime,infile,obstype, &
                        lunout,twind,sis,nobs_sub1(1,i))
                   string='READ_NSSTBUFR'
                elseif ( platid == 'mods') then
                   call read_modsbufr(nread,npuse,nouse,gstime,infile,obstype, &
                        lunout,twind,sis,nobs_sub1(1,i))
                   string='READ_MODSBUFR'
                else
                   if(nst_gsi>0)then
                      write(6,*)'read_obs: should not handle SST via read_prepbufr when NSST on'
                      call stop2(999)
                   endif

                   call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
                        prsl_full,nobs_sub1(1,i),read_rec(i))
                   string='READ_PREPBUFR'
                endif

!            Process radar reflectivity Mosaic
             else if (obstype == 'rad_ref' ) then
                call read_radarref_mosaic(nread,npuse,infile,obstype,lunout,twind,sis, &
                      nobs_sub1(1,i))
                string='READ_RADARREF_MOSAIC'

!            Process  lightning
             else if (obstype == 'lghtn' ) then
                if(i_gsdcldanal_type==2) then
                   call read_lightning(nread,npuse,infile,obstype,lunout,twind,sis,nobs_sub1(1,i))
                else if(i_gsdcldanal_type==1 .or. i_gsdcldanal_type==6 &
                        .or. i_gsdcldanal_type==3 .or. i_gsdcldanal_type==7) then
                   call read_lightning_grid(nread,npuse,infile,obstype,lunout,twind,sis,nobs_sub1(1,i))
                endif
                string='READ_LIGHTNING'

!            Process  NASA LaRC 
             ! for regional obs that are already mapped to analysis grid
             else if (obstype == 'larccld' ) then
                if(i_gsdcldanal_type==2) then
                   call read_NASA_LaRC_cloud(nread,npuse,nouse,infile,obstype,lunout,sis,nobs_sub1(1,i))
                else if(i_gsdcldanal_type==1 .or. i_gsdcldanal_type==6 &
                        .or. i_gsdcldanal_type==3 .or. i_gsdcldanal_type==7 .or. i_gsdcldanal_type==99) then
                   call read_nasa_larc(nread,npuse,infile,obstype,lunout,twind,sis,nobs_sub1(1,i))
                end if
                string='READ_NASA_LaRC'
             ! for global NASA LaRC obs 
             else if (obstype == 'larcglb' ) then
                call read_NASA_LaRC_cloud(nread,npuse,nouse,infile,obstype,lunout,twind,sis,nobs_sub1(1,i))
                string='READ_NASA_LaRC'

!            Process radar winds
             else if (obstype == 'rw') then
                if( trim(infile) == 'vr_vol' )then
                  call read_radar_wind_ascii(nread,npuse,nouse,infile,lunout,obstype,sis,&
                                  hgtl_full,nobs_sub1(1,i))
                else
                 if (vadwnd_l2rw_qc) then
                      write(6,*)'READ_OBS: radial wind,read_radar,dfile=',infile,',dsis=',sis 
                     call read_radar(nread,npuse,nouse,infile,lunout,obstype,twind,sis,&
                                     hgtl_full,nobs_sub1(1,i))
                     string='READ_RADAR'
                  else if (sis == 'l2rw') then
                      write(6,*)'READ_OBS: radial wind,read_radar_l2rw_novadqc,dfile=',infile,',dsis=',sis
                     call read_radar_l2rw_novadqc(npuse,nouse,lunout,obstype,sis,nobs_sub1(1,i))
                     string='READ_RADAR_L2RW_NOVADQC'
                  end if
                end if
                string='READ_RADAR_WIND'

!            Process radar reflectivity from MRMS
             else if (obstype == 'dbz' ) then
                print *, "calling read_dbz"
                if(trim(infile)=='dbzobs.nc')then
                  call read_dbz_nc(nread,npuse,nouse,infile,lunout,obstype,sis,hgtl_full,nobs_sub1(1,i))
                  string='READ_dBZ'
                else
                  call read_dbz_mrms_detect_format(infile,l_mrms_sparse_netcdf)
                  if(l_mrms_sparse_netcdf) then
                     call read_dbz_mrms_sparse_netcdf(nread,npuse,nouse,infile,obstype,lunout,sis,nobs_sub1(1,i))
                     string='READ_dbz_mrms_sparse_netcdf'
                  else
                     call read_dbz_mrms_netcdf(nread,npuse,nouse,infile,obstype,lunout,sis,nobs_sub1(1,i))
                     string='READ_dbz_mrms_netcdf'
                  endif
                end if

!            Process lagrangian data
             else if (obstype == 'lag') then
                call read_lag(nread,npuse,nouse,infile,lunout,obstype,&
                 &twind,gstime,sis)
                string='READ_LAG'

!            Process lidar winds
             else if (obstype == 'dw') then
                call read_lidar(nread,npuse,nouse,infile,obstype,lunout,twind,sis, &
                     nobs_sub1(1,i))
                string='READ_LIDAR'

!            Process synthetic tc-mslp obs
             else if (obstype == 'tcp') then
                call read_tcps(nread,npuse,nouse,infile,obstype,lunout,sis, &
                     nobs_sub1(1,i))
                string='READ_TCPS'

             else if (obstype == 'pm2_5' .or. obstype == 'pm10') then

                if (oneobtest_chem .and. oneob_type_chem=='pm2_5') then
                   call oneobschem(nread,npuse,nouse,gstime,&
                        &infile,obstype,lunout,sis,nobs_sub1(1,i))
                   string='ONEOBSCHEM'
                else
                   call read_anowbufr(nread,npuse,nouse,gstime,&
                        &infile,obstype,lunout,twind,sis,nobs_sub1(1,i))
                   string='READ_ANOWBUFR'
                endif

!            Process pblh
             else if (obstype == 'pblh') then
                call read_pblh(nread,npuse,nouse,infile,obstype,lunout,twind,sis, &
                    nobs_sub1(1,i))
                string='READ_PBLH'
             end if conv_obstype_select
!            Process swcp and lwcp
          else if (ditype(i) == 'wcp') then
             if ( obstype == 'swcp' .or. obstype == 'lwcp' ) then
                call read_wcpbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis, &
                   prsl_full,nobs_sub1(1,i),read_rec(i))
                string='READ_WCPBUFR'
             end if
          else if (ditype(i) == 'rad')then

             call radiance_obstype_search(obstype,radmod)

!            Process TOVS 1b data
             rad_obstype_select: &
             if (platid /= 'aqua' .and. (obstype == 'amsua' .or. &
                  obstype == 'amsub' .or. obstype == 'msu'   .or.  &
                  obstype == 'mhs'   .or. obstype == 'hirs4' .or.  &
                  obstype == 'hirs3' .or. obstype == 'hirs2' .or.  &
                  obstype == 'ssu' )) then
                call read_bufrtovs(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), nobs_sub1(1,i), &
                     read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use,radmod)
                string='READ_BUFRTOVS'

!            Process atms data
             else if (obstype == 'atms') then
                call read_atms(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i),&
                     read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use,radmod)
                string='READ_ATMS'

!            Process saphir data
             else if (obstype == 'saphir') then
                call read_saphir(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),  &
                     nobs_sub1(1,i),dval_use)
                string='READ_SAPHIR'


!            Process airs data        
             else if(platid == 'aqua' .and. (obstype == 'airs' .or.   &
                  obstype == 'amsua'  .or.  obstype == 'hsb' ))then
                call read_airs(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                     nobs_sub1(1,i),read_rec(i),dval_use)
                string='READ_AIRS'

!            Process iasi data
             else if(obstype == 'iasi')then
                call read_iasi(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i), &
                     read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use)
                string='READ_IASI'

!            Process cris data
             else if(obstype == 'cris' .or. obstype =='cris-fsr' )then
                call read_cris(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),nobs_sub1(1,i), &
                     read_rec(i),read_ears_rec(i),read_db_rec(i),dval_use)
                string='READ_CRIS'

!            Process GOES sounder data
!            Process raw or prepbufr files (1x1 or 5x5)
             else if (obstype == 'sndr' .or.                            &
                      obstype == 'sndrd1' .or. obstype == 'sndrd2' .or. &
                      obstype == 'sndrd3' .or. obstype == 'sndrd4') then
                call read_goesndr(mype,val_dat,ithin,rmesh,platid,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,gstime,sis,&
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                     nobs_sub1(1,i),read_rec(i),dval_use)
                string='READ_GOESNDR'

!            Process ssmi data
             else if (obstype == 'ssmi' ) then
                call read_ssmi(mype,val_dat,ithin,rmesh,platid,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                     nobs_sub1(1,i),read_rec(i),dval_use)
                string='READ_SSMI'

!            Process amsre data
             else if ( obstype == 'amsre_low' .or. obstype == 'amsre_mid' .or. &
                       obstype == 'amsre_hig' ) then
                call read_amsre(mype,val_dat,ithin,isfcalc,rmesh,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                     nobs_sub1(1,i),read_rec(i),dval_use)
                string='READ_AMSRE'
                
!            Process ssmis data
             else if (obstype == 'ssmis'     .or. &
                      obstype == 'ssmis_las' .or. obstype == 'ssmis_uas' .or. &
                      obstype == 'ssmis_img' .or. obstype == 'ssmis_env' ) then
                call read_ssmis(mype,val_dat,ithin,isfcalc,rmesh,platid,gstime,&
                      infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                      mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                      nobs_sub1(1,i),read_rec(i),dval_use)
                string='READ_SSMIS'

!            Process AMSR2 data
             else if(obstype == 'amsr2')then
                call read_amsr2(mype,val_dat,ithin,rmesh,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),  &
                     nobs_sub1(1,i))
                string='READ_AMSR2'

!            Process GOES IMAGER RADIANCE  data
             else if(obstype == 'goes_img') then
                call read_goesimg(mype,val_dat,ithin,rmesh,platid,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                     nobs_sub1(1,i),read_rec(i),dval_use)
                string='READ_GOESMIMG'

!            Process GMI data
             else if (obstype == 'gmi') then
                call read_gmi(mype,val_dat,ithin,rmesh,platid,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis,&
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),  &
                     nobs_sub1(1,i),dval_use)
                string='READ_GMI'

!            Process Meteosat SEVIRI RADIANCE  data
             else if(obstype == 'seviri') then
                 call read_seviri(mype,val_dat,ithin,rmesh,platid,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                     nobs_sub1(1,i),read_rec(i),dval_use)
                string='READ_SEVIRI'
!            Process GOES-R ABI RADIANCE  data
             else if(obstype == 'abi') then
                call read_abi(mype,val_dat,ithin,rmesh,platid,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                     nobs_sub1(1,i),read_rec(i),dval_use)
                string='READ_ABI'

        !    Process Himawari-8 AHI RADIANCE  data
             else if(obstype == 'ahi') then
                call read_ahi(mype,val_dat,ithin,rmesh,platid,gstime,&
                     infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                     mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i),  &
                     nobs_sub1(1,i))
                string='READ_AHI'


!            Process NAVY AVHRR RADIANCE  data
               else if(obstype == 'avhrr_navy') then
                  call read_avhrr_navy(mype,val_dat,ithin,rmesh,platid,gstime,&
                       infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                       mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                       nobs_sub1(1,i),read_rec(i),dval_use)
                  string='READ_AVH_NAVY'

  !            Process NESDIS AVHRR RADIANCE  data
               else if(obstype == 'avhrr') then
                  call read_avhrr(mype,val_dat,ithin,rmesh,platid,gstime,&
                       infile,lunout,obstype,nread,npuse,nouse,twind,sis, &
                       mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                       nobs_sub1(1,i),read_rec(i),dval_use)
                  string='READ_AVHRR'
               end if rad_obstype_select

!         Process ozone data
          else if (ditype(i) == 'ozone')then
             ozone_obstype_select: &
             if (is_extOzone(infile,obstype,dplat(i))) then

                call extOzone_read(infile,obstype,dplat(i),dsis(i), &
                   iread,ipuse,iouse, platid,gstime,lunout,twind,ithin,rmesh, &
                   nobs_sub1(:,i))
                string='extOzone_read'

                nread=nread+iread
                npuse=npuse+ipuse
                nouse=nouse+iouse

             else
                call read_ozone(nread,npuse,nouse,&
                   platid,infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, &
                   nobs_sub1(1,i))
                string='READ_OZONE'
             endif ozone_obstype_select

!         Process co data
          else if (ditype(i) =='co')then 
             call read_co(nread,npuse,nouse,infile,gstime,lunout,obstype,sis, &
                  nobs_sub1(1,i))
             string='READ_CO'

!         Process precipitation             
          else if (ditype(i) == 'pcp')then
             call read_pcp(nread,npuse,nouse,gstime,infile,lunout,obstype,twind,sis, &
                  nobs_sub1(1,i))
             string='READ_PCP'

!         Process gps observations
          else if (ditype(i) == 'gps')then
             call read_gps(nread,npuse,nouse,infile,lunout,obstype,twind, &
                  nprof_gps1,sis,nobs_sub1(1,i))
             string='READ_GPS'

!         Process aerosol data
          else if (ditype(i) == 'aero' )then
             call read_aerosol(nread,npuse,nouse,&
                  platid,infile,gstime,lunout,obstype,twind,sis,ithin,rmesh, &
                  mype_root,mype_sub(mm1,i),npe_sub(i),mpi_comm_sub(i), &
                  nobs_sub1(1,i))
             string='READ_AEROSOL'

!         Process satellite lightning observations (e.g. GOES/GLM)                     
          else if(ditype(i) == 'light')then
               if (obstype == 'goes_glm' ) then
             call read_goesglm(nread,ndata,nodata,infile,obstype,lunout,sis)
             string='READ_GOESGLM'
               endif

          end if ditype_select

!         Close unit to data file

!         Accumulate data counts on "root" task
          if (mype_sub(mm1,i)==mype_root) then
             close(lunout)
             ndata1(i,1)=ndata1(i,1)+npuse
             ndata1(i,2)=ndata1(i,2)+nread
             ndata1(i,3)=ndata1(i,3)+nouse

             if(string(1:2)=='--') then
               call warn('read_obs','task without reader, i =',i)
               call warn('read_obs','                infile =',trim(infile))
               call warn('read_obs','             ditype(i) =',trim(ditype(i)))
               call warn('read_obs','               obstype =',trim(obstype))
               call warn('read_obs','                platid =',trim(platid))
               call warn('read_obs','                string =',trim(string))
             endif

             write(6,8000) adjustl(string),infile,obstype,sis,nread,ithin,&
                  rmesh,isfcalc,nouse,npe_sub(i)
8000         format(1x,a22,': file=',a15,&
                  ' type=',a10,  ' sis=',a20,  ' nread=',i10,&
                  ' ithin=',i2, ' rmesh=',f11.6,' isfcalc=',i2,&
                  ' nkeep=',i10,' ntask=',i3)

          endif
       endif task_belongs

    end do loop_of_obsdata_files
    deallocate(prsl_full)
    deallocate(hgtl_full)

!   Broadcast aircraft new tail numbers for aircraft
!   temperature bias correction
!   if (aircraft_t_bc) then
!      call mpi_barrier(mpi_comm_world,ierror)
!      call mpi_bcast(ntail_update,1,mpi_itype,mype_airobst,mpi_comm_world,ierror)
!      call mpi_bcast(idx_tail,max_tail,mpi_itype,mype_airobst,mpi_comm_world,ierror)
!      call mpi_bcast(taillist,max_tail,MPI_CHARACTER,mype_airobst,mpi_comm_world,ierror)
!   end if

!   Deallocate arrays containing full horizontal surface fields
    call destroy_sfc
!   Deallocate arrays containing full horizontal nsst fields
    if (nst_gsi > 0) call gsi_nstcoupler_final()
!   Sum and distribute number of obs read and used for each input ob group
    call mpi_allreduce(ndata1,ndata,ndat*3,mpi_integer,mpi_sum,mpi_comm_world,&
       ierror)

!   Collect super obs factors
    if(mype == 0)write(6,*) ' dval_use = ',dval_use
    if(dval_use)then
       call mpi_allreduce(super_val,super_val1,superp+1,mpi_rtype,&
            mpi_sum,mpi_comm_world,ierror)
    else
       super_val1=zero
    end if
    super_val1(0)=one
    deallocate(super_val)

!   Collect number of gps profiles (needed later for qc)
    call mpi_allreduce(nprof_gps1,nprof_gps,1,mpi_integer,mpi_sum,mpi_comm_world,ierror)
    call mpi_allreduce(nobs_sub1,nobs_sub,npe*ndat,mpi_integer,mpi_sum,mpi_comm_world,& 
         ierror)

!   Write collective obs selection information to scratch file.
    if (lread_obs_save .and. mype==0) then
       write(6,*)'READ_OBS:  write collective obs selection info to ',trim(obs_input_common)
       call unformatted_open(lunsave,file=obs_input_common,class=".obs_input.")
       write(lunsave) ndata,ndat,npe,superp,nprof_gps,ditype
       write(lunsave) super_val1
       write(lunsave) nobs_sub
       close(lunsave)
    endif

!   End of routine
    return
end subroutine read_obs


end module read_obsmod