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
!
! 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
!
!   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,r_kind,r_double
  use gsi_4dvar, only: iadatebgn,iadateend
  use obsmod, only: offtime_data
  use convinfo, only: nconvtype,ictype,ioctype,icuse

  implicit none

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

  logical :: lhere
  integer(i_kind) :: lenb,iret
  character(len=256) command, fname
  
#ifdef _INTEL_11_0_083_
  lenb=0; lbytes = lenb
  inquire(file=trim(filename),exist=lhere)
  if(lhere)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,*) lenb
    close(999)
    lbytes=lenb
  endif
#else
  inquire(file=trim(filename),exist=lhere,size=lbytes)
#endif
  lexist=.false.
  if(lhere)then
     lexist=lbytes>0_i_llong
  end if
  return
end subroutine gsi_inquire

subroutine read_obs_check (lexist,filename,jsatid,dtype,minuse)
!$$$  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-??-??  derber   - originally placed inside inquire
!   2009-01-05  todling  - move time/type-check out of inquire
!
!   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

  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

  logical :: lhere
  integer(i_kind) :: lenb,lnbufr,idate,idate2,iret,kidsat
  integer(i_kind) :: ireadsb,ireadmg,kx,nc
  real(r_double) :: satid,rtype,bmixx
  character(len=256) command, fname
  character(8) subset

  idate=0
  if(trim(dtype) == 'tcp')return
! RTod: For some odd reason the block below does not work on the GMAO Linux Cluster
!       Anyone else on a Linux Cluster willing to try?
#ifdef _INTEL_11_0_083_
  return
#else
! Use routine as usual
  if(lexist)then
      lnbufr = 15
      open(lnbufr,file=trim(filename),form='unformatted',status ='unknown')
      call openbf(lnbufr,'IN',lnbufr)
      !!bmixx=10e10; call setbmiss(bmixx)
      call datelen(10)
      call readmg(lnbufr,subset,idate,iret)

!     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'
         end if
         write(6,*)'Analysis start  :',iadatebgn
         write(6,*)'Analysis end    :',iadateend
         write(6,*)'Observation time:',idate
         if(.not.offtime_data) lexist=.false.
      endif
      kidsat=0
      if(jsatid == 'metop-a')kidsat=4
      if(jsatid == 'metop-b')kidsat=5
      if(jsatid == 'metop-c')kidsat=6
      if(jsatid == 'n08')kidsat=200
      if(jsatid == 'n09')kidsat=201
      if(jsatid == 'n10')kidsat=202
      if(jsatid == 'n11')kidsat=203
      if(jsatid == 'n12')kidsat=204
      if(jsatid == 'n14')kidsat=205
      if(jsatid == 'n15')kidsat=206
      if(jsatid == 'n16')kidsat=207
      if(jsatid == 'n17')kidsat=208
      if(jsatid == 'n18')kidsat=209
      if(jsatid == 'n19')kidsat=223
      if(jsatid == 'f08')kidsat=241
      if(jsatid == 'f10')kidsat=243
      if(jsatid == 'f11')kidsat=244
      if(jsatid == 'f13')kidsat=246
      if(jsatid == 'f14')kidsat=247
      if(jsatid == 'f15')kidsat=248
      if(jsatid == 'f16')kidsat=249
      if(jsatid == 'f17')kidsat=250
      if(jsatid == 'g08' .or. jsatid == 'g08_prep')kidsat=252
      if(jsatid == 'g09' .or. jsatid == 'g09_prep')kidsat=253
      if(jsatid == 'g10' .or. jsatid == 'g10_prep')kidsat=254
      if(jsatid == 'g11' .or. jsatid == 'g11_prep')kidsat=255
      if(jsatid == 'g12' .or. jsatid == 'g12_prep')kidsat=256
      if(jsatid == 'g13' .or. jsatid == 'g13_prep')kidsat=257
      if(jsatid == 'n05')kidsat=705
      if(jsatid == 'n06')kidsat=706
      if(jsatid == 'n07')kidsat=707
      if(jsatid == 'tirosn')kidsat=708

      if(lexist)then
       if(kidsat /= 0)then
        lexist=.false.
        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
           end if
        end do
       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 
         end do fileloop
       end if
      end if

      call closbf(lnbufr)
      write(6,*)'read_obs_check: bufr file date is ',idate,trim(filename),' ',dtype
  else
      write(6,*)'read_obs_check: bufr file date is ',idate,trim(filename),' ',dtype,' not used '
  end if
#endif /* ibm_sp */
  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
!
!   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
    use gridmod, only: nlon,nlat,nsig,iglobal,ijn,itotsub,lat1,lon1,&
         ltosi,ltosj,displs_g,strip,reorder
    use obsmod, only: iadate,ndat,time_window,dplat,dsfcalc,dfile,dthin, &
           dtype,dval,dmesh,obsfile_all,ref_obs,nprof_gps,dsis,ditype,&
           oberrflg,perturb_obs,lobserver,lread_obs_save,obs_input_common
    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 guess_grids, only: ges_prsl,ntguessig,destroy_sfc_grids
    use radinfo, only: nusis,iuse_rad,jpch_rad,diag_rad
    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

    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,use_prsl_full
    logical :: use_sfc,nuse
    logical,dimension(ndat):: belong,parallel_read
    character(10):: obstype,platid
    character(13):: string,infile
    character(16):: filesave
    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
    integer(i_kind):: iworld,iworld_group,next_mype,mm1,iix
    integer(i_kind):: mype_root,ntask_read,mpi_comm_sub_read,lll,llb
    integer(i_kind):: mype_sub_read,minuse,lunsave
    integer(i_kind):: iworld_group_r1,iworld_r1,iworld_group_r2,iworld_r2
    integer(i_kind),dimension(ndat):: npe_sub,npe_sub3,mpi_comm_sub,mype_root_sub,npe_order
    integer(i_kind),dimension(ndat):: mpi_comm_sub_r1,mpi_comm_sub_r2
    integer(i_kind),dimension(ndat):: ntasks1,ntasks
    integer(i_kind),dimension(ndat,3):: ndata1
    integer(i_kind),dimension(npe,ndat):: mype_work,mype_work_r1,mype_work_r2
    integer(i_kind),dimension(npe,ndat):: mype_sub,mype_sub_r1,mype_sub_r2
    integer(i_kind),allocatable,dimension(:):: nrnd

    real(r_kind) gstime,val_dat,rmesh,twind,rseed
    real(r_kind),dimension(lat1*lon1):: prslsm
    real(r_kind),dimension(max(iglobal,itotsub)):: work1
    real(r_kind),allocatable,dimension(:,:,:):: prsl_full

    data lunout / 81 /
    data lunsave  / 82 /

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

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

    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(oberrflg .or. perturb_obs) then
       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
    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
!  Control parallel read for each ob type (currently just rad obs).  
!  To remove parallel read comment out line.
       if(hirs .and. dthin(i) > 0)parallel_read(i)= .true.
       if(sndr .and. dthin(i) > 0)parallel_read(i)= .true.
       if(avhrr .and. dthin(i) > 0)parallel_read(i)= .true.
       if(amsre .and. dthin(i) > 0)parallel_read(i)= .true.
       if(ssmis .and. dthin(i) > 0)parallel_read(i)= .true.
       if(obstype == 'ssmi' .and. dthin(i) > 0)parallel_read(i)= .true.
       if(obstype == 'airs' .and. dthin(i) > 0)parallel_read(i)= .true.
       if(obstype == 'amsub' .and. dthin(i) > 0)parallel_read(i)= .true.
       if(obstype == 'hsb' .and. dthin(i) > 0)parallel_read(i)= .true.
       if(obstype == 'iasi' .and. dthin(i) > 0)parallel_read(i)= .true.
       if(obstype == 'amsua' .and. dthin(i) > 0)parallel_read(i)= .true.
       if(obstype == 'iasi' .and. dthin(i) > 0)parallel_read(i)= .true.
       if(obstype == 'mhs' .and. dthin(i) > 0)parallel_read(i)= .true.
       if(obstype == 'goes_img' .and. dthin(i) > 0)parallel_read(i)= .true.
       if(obstype == 'ssu' .and. dthin(i) > 0)parallel_read(i)= .true.
       if (obstype == 't'  .or. obstype == 'uv' .or. &
           obstype == 'q'  .or. obstype == 'ps' .or. &
           obstype == 'pw' .or. obstype == 'spd'.or. &
           obstype == 'sst'.or. obstype == 'srw'.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' )  then
          ditype(i) = 'conv'
       else if( hirs   .or. sndr      .or.  &
               obstype == 'seviri'    .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.  &
               avhrr .or.  &
               amsre  .or. ssmis      .or. obstype == 'ssmi'      .or.  &
               obstype == 'ssu' ) then
          ditype(i) = 'rad'
       else if (obstype == 'sbuv2' .or. obstype == 'omi' &
           .or. obstype == 'gome'  .or. obstype == 'o3lev') then
          ditype(i) = 'ozone'
       else if (index(obstype,'pcp')/=0 )then
          ditype(i) = 'pcp'
       else if (obstype == 'gps_ref' .or. obstype == 'gps_bnd') then
          ditype(i) = 'gps'
       else
          write(6,*)'READ_OBS:  ***ERROR*** - unknown ob type ',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 used.

       nuse=.false.
       minuse=-1
       if(ditype(i) == 'conv')then
          if(diag_conv)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)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)minuse=-2
          if (dtype(i) == 'o3lev') then
             do j=1,nconvtype
                if(trim(dtype(i)) == trim(ioctype(j)) .and. icuse(j) > minuse)nuse=.true.
             end do
          else
             do j=1,jpch_oz
                if(trim(dsis(i)) == trim(nusis_oz(j)) .and. iuse_oz(j) > minuse)nuse=.true.
             end do
          endif
       else if(ditype(i) == 'pcp')then
          if(diag_pcp)minuse=-2
             do j=1,npcptype
                if(trim(dsis(i)) == trim(nupcp(j)) .and. iusep(j) > minuse)nuse=.true.
             end do
          else
             nuse=.true.
          end if
       if(.not. nuse)then
          if(mype == 0)write(6,*) 'data type ',dsis(i), &
                'not used in info file -- do not read file',dfile(i)
       end if


!   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,dfile(i),mype)
          call read_obs_check (lexist,dfile(i),dplat(i),dtype(i),minuse)

          len4file=lenbytes/4
          if (ditype(i) == 'rad' .and. nuse .and.           &
                    dplat(i) /= 'aqua' .and. dplat(i) /= 'metop-a' .and. &
                   (obstype == 'amsua' .or.  obstype == 'amsub' .or.     &
                    obstype == 'mhs' )) then
!                   obstype == 'mhs'   .or. hirs )) then

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

             lexist=lexist .or. lexistears
             len4file=len4file+lenbytes/4
          end if
!      Initialize number of reader tasks to 1.  For the time being
!      only allow number of reader tasks >= 1 for select obstype.

          if(lexist .and. nuse) then
             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
    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)

    npemax=0
    npetot=0
    do i=1,ndat
       npe_sub(i)=ntasks(i)
       npetot=npetot+npe_sub(i)
       npemax=max(npemax,npe_sub(i))
    end do

    if(l4dvar.and.(.not.lobserver)) return
    
    npeextra=npe-mod(npetot,npe)
    if(npeextra > 0)then
       if(mype == 0)write(6,*) ' number of extra processors ',npeextra
       npe_sub3=npe_sub
       extraloop: do j=1,npeextra
          iix=1
          do ii=1,5
             do i=1,ndat
                if(iix == npe_sub3(i) .and. parallel_read(i))then
                   if(ntasks(i) > 0 .and. ntasks(i) <= npeextra)then
                      npeextra=npeextra-ntasks(i)
                      npe_sub(i)=npe_sub(i)+ntasks(i)
                      ntasks(i)=2*ntasks(i)
                      if(npeextra < iix)cycle extraloop
                   end if
                end if
             end do
          end do
          iix=max(min(2*iix,8),npeextra)
       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
       if(npestart == npe)npestart=0
       npe_sub3(ilarge)=0
    end do loopx

        
!   Define sub-communicators for each data file
    mm1=mype+1
    belong=.false.
    mype_sub=-999
    mype_sub_r1=-999
    mype_sub_r2=-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
    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,3i4)') &
        'READ_OBS:  read ',i,dtype(i),dsis(i),' using ntasks=',ntasks(i),mype_root_sub(i),npe_sub(i) 
    end do


    use_prsl_full=.false.
    use_sfc=.false.
    do i=1,ndat
       if(belong(i) .and. ditype(i) =='conv')then
          obstype=dtype(i)                  
          if(obstype /= 'dw' .and. obstype /= 'rw' .and. obstype /= 'srw')then
             use_prsl_full=.true.
          end if
       else if(belong(i) .and. (ditype(i) == 'rad' .or. ditype(i)=='pcp'))then
          use_sfc=.true.
       end if
    end do
!   Get guess 3d pressure on full grid
    if(use_prsl_full)allocate(prsl_full(nlat,nlon,nsig))
    do k=1,nsig
       call strip(ges_prsl(1,1,k,ntguessig),prslsm,1)
       call mpi_allgatherv(prslsm,ijn(mype+1),mpi_rtype,&
            work1,ijn,displs_g,mpi_rtype,mpi_comm_world,ierror)
       if(use_prsl_full)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
!   Create full horizontal surface fields from local fields in guess_grids
    call getsfc(mype,use_sfc)
    if(use_sfc) call prt_guessfc2('sfcges2')

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

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

          platid=dplat(i)                    !     platid   - satellites to read
          obstype=dtype(i)                   !     obstype  - observation types to process
          infile=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
          if(ditype(i) == 'conv')then
             if (obstype == 't'  .or. obstype == 'uv' .or. &
                 obstype == 'q'  .or. obstype == 'ps' .or. &
                 obstype == 'pw' .or. obstype == 'spd'.or. & 
                 obstype == 'mta_cld' .or. obstype == 'gos_ctp'  ) then
                call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
                     prsl_full)
                string='READ_PREPBUFR'

!            Process conventional SST (modsbufr, at this moment) data
             elseif ( obstype == 'sst' ) then
                if ( platid == 'mods') then
                   call read_modsbufr(nread,npuse,nouse,gstime,infile,obstype, &
                        lunout,twind,sis)
                   string='READ_MODSBUFR'
                else
                   call read_prepbufr(nread,npuse,nouse,infile,obstype,lunout,twind,sis,&
                        prsl_full)
                   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)
                string='READ_RADARREF_MOSAIC'

!            Process  lightning
             else if (obstype == 'lghtn' ) then
                call read_lightning(nread,npuse,infile,obstype,lunout,twind,sis)
                string='READ_LIGHTNING'

!            Process  NASA LaRC 
             else if (obstype == 'larccld' ) then
                call read_NASA_LaRC(nread,npuse,infile,obstype,lunout,twind,sis)
                string='READ_NASA_LaRC'

!            Process radar winds
             else if (obstype == 'rw') then
                call read_radar(nread,npuse,nouse,infile,lunout,obstype,twind,sis)
                string='READ_RADAR'

!            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)
                string='READ_LIDAR'

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

!            Process radar superob winds
             else if (obstype == 'srw') then
                call read_superwinds(nread,npuse,nouse,infile,obstype,lunout, &
                     twind,sis)
                string='READ_SUPRWNDS'
             end if

          else if (ditype(i) == 'rad')then


!            Process TOVS 1b data
             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
                llb=1
                lll=1
                if((obstype == 'amsua' .or. obstype == 'amsub' .or. obstype == 'mhs') .and. &
                   (platid /= 'metop-a' .or. platid /='metop-b' .or. platid /= 'metop-c'))lll=2
                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),llb,lll)
                string='READ_BUFRTOVS'

!            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))
                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))
                string='READ_IASI'

!            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))
                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))
                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))
                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))
                string='READ_SSMIS'

!            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))
                string='READ_GOESMIMG'

!            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))
                string='READ_SEVIRI'

!            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))
                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))
                string='READ_AVHRR'
             end if

!         Process ozone data
          else if (ditype(i) == 'ozone')then
             call read_ozone(nread,npuse,nouse,&
                  platid,infile,gstime,lunout,obstype,twind,sis,ithin,rmesh)
             string='READ_OZONE'

!         Process precipitation             
          else if (ditype(i) == 'pcp')then
             call read_pcp(nread,npuse,nouse,gstime,infile, &
                  lunout,obstype,twind,sis)
             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)
             string='READ_GPS'
             
          end if

!         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

             write(6,8000) adjustl(string),infile,obstype,sis,nread,ithin,&
                  rmesh,isfcalc,nouse,npe_sub(i)
8000         format(a13,': file=',a10,&
                  ' type=',a10,  ' sis=',a20,  ' nread=',i10,&
                  ' ithin=',i2, ' rmesh=',f7.3,' isfcalc=',i2,&
                  ' ndata=',i10,' ntask=',i3)

          endif
       endif

    end do
    if(use_prsl_full)deallocate(prsl_full)

!   Deallocate arrays containing full horizontal surface fields
    call destroy_sfc

!   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
    call mpi_allreduce(super_val,super_val1,superp+1,mpi_rtype,&
         mpi_sum,mpi_comm_world,ierror)
    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)

!   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)
       open(lunsave,file=obs_input_common,form='unformatted')
       write(lunsave) ndata,superp,nprof_gps,ditype
       write(lunsave) super_val1
       close(lunsave)
    endif

!   End of routine
    return
end subroutine read_obs


end module read_obsmod