    module readobs
!  READOBS - routines for reading and dispatching the observations, used by
!  preanalysis.
!
!  06/29/05 - scale the observations to energy units.
!  01/13/05 - added filter by observation type.  For now, the first letter
!     of the FIELD component is used in each RAW_OBS type to determine whether
!     the obeervation is surfrace pressure, temperature, or u- or v- component
!     of the wind.  (These are the only observation types that we can handle
!     at the moment.)  Only the observations matching one of the selected types
!     are distributed to the processors.
!  Eventually we'll simply dispatch each observation to the appropriate
!  processors with MPI, but for now will buffer them and write them
!  out in files for distribution.
!
!  Required packages:
!
    use grid_obs
    use ferror
    use io
    use fsplit
    use resolution,only:LATG
    implicit none
!
    private::flush
    contains
!-----------------------------
    subroutine read_observations(obsfn,outfn,maxproc,ywin,list)
!  READ_OBSERVATIONS - read the observations in RAW_OBS format,
!  then create observations in the format that we need for the data
!  assimilation.
!  OBSFN : name of the file containing the observations.
!  OUTFN : name of the output file. If NPROC > 1, then OUTFN is a template name.
!  STARTLAT : Kth entry is lowest latitude that Kth processor will handle.
!  ENDLAT : Kth entry is highest latitude that Kth processor will handle.
!  MAXPROC : max number of processors to use.  The observations for the
!     data assimilation are broken into as many as MAXPROC chunks.
!  YWIN : Y (latitude) radius to allow for each strip.
!  LIST : some subset of PTUV, which determines the obs types to retain.
!  No other preprocessing or scaling of the observations is done, other
!  than to divide the observations into latitude strips that encircle
!  the globe and to retain only those observations of a given type
!  (see FILTER_BY_TYPE).
!
    implicit none
    character(*),intent(in)::obsfn,outfn,list
    integer,intent(in)::maxproc,ywin
!
!  Local variables
!
    type(raw_obs)::yobs
    integer,parameter::INITBUFSIZE=100000  ! must be > 1
    integer::nproc,nobs,idate,istat
    type(raw_obs),allocatable::obsbuf(:)
    integer,allocatable::iwidth(:),startlat(:),endlat(:),psouth(:),pnorth(:)
!
    if(verbosity().gt.VERB_SILENT) write(LOGFILE,100) trim(obsfn)
100 format('preanalysis: read_observations: input file: ',a)
    open(INNOSWAP,file=obsfn,form='unformatted',status='old',iostat=istat)
    if(istat.ne.0) call fexit('read_observations: cannot open'//trim(obsfn),2)
!
    nproc=maxproc
    call imod(nproc,iwidth,LATG,ywin,startlat,endlat,psouth,pnorth)
    if(nproc.lt.maxproc) &
       call finfo('read_observations: fewer processors used than requested')
    allocate(obsbuf(INITBUFSIZE))
!
!  The first record is the date in yyyymmddhh format, followed by one
!  observation per record.
!
    read(INNOSWAP,end=10) idate
    if(verbosity().gt.VERB_SILENT) write(LOGFILE,110) idate
110 format('date (yyyymmddhh): ',i12)
    nobs=0
    do
       read(INNOSWAP,end=10) yobs
       nobs=nobs+1
       if(nobs.gt.size(obsbuf)) call enlarge_obsbuf
       obsbuf(nobs)=yobs
    enddo
10  continue
    if(verbosity().gt.VERB_SILENT) write(LOGFILE,111) nobs,'input'
111 format('preanalysis: ',i0,' observations ',a)
!
    call filter_by_type(list,obsbuf,nobs)
    if(verbosity().gt.VERB_SILENT) write(LOGFILE,111) nobs,'retained'
    call flush(obsbuf,nobs,outfn,startlat,endlat,nproc)
    deallocate(iwidth,startlat,endlat)
    close(INNOSWAP)
    return
    contains
    !-------
       subroutine enlarge_obsbuf  ! by 50% or one element, whichever is greater
       !
       type(raw_obs)::otmp(size(obsbuf))  ! can cause stack overflow
       integer::oldsize,newsize
       !
       oldsize=size(obsbuf)
       newsize=max((oldsize*3)/2,oldsize+1)
       otmp=obsbuf
       deallocate(obsbuf)
       allocate(obsbuf(newsize))
       obsbuf(1:oldsize)=otmp
       return
       end subroutine enlarge_obsbuf
    end subroutine read_observations    
!-----------------------------
    subroutine filter_by_type(list,obsbuf,nobs)
!  FILTER_BY_TYPE - omit observations whose FIELD variable does not match
!  one of those in LIST.
!  NOBS  :=: on input, the number of observations in OBSBUF.  On return, the
!     number retained.
!  OBSBUF :=: on return, only the selected observations are retained.
!
    character(*),intent(in)::list
    integer,intent(inout)::nobs
    type(raw_obs),intent(inout)::obsbuf(nobs)
!
    integer::j,n
    character(4)::field
!
    n=0
    do j=1,nobs
       field=adjustl(obsbuf(j)%field)
       if(index(list,field(1:1)).gt.0) then
          n=n+1
          if(n.lt.j) obsbuf(n)=obsbuf(j)
       endif
    enddo
    if(n.lt.nobs) write(LOGFILE,100) n,nobs
100 format('filter_by_type: retaining ',i0,' obs of ',i0/ &
       ' input based on obs type')
    nobs=n
    return
    end subroutine filter_by_type
!-----------------------------
    subroutine flush(obsbuf,nobs,outfn,startlat,endlat,nproc)
!  FLUSH - flush the observation buffer by distributing all the observations
!  to each processor in latitude strips.
!  OBSBUF : buffer of NOBS model_obs to be distributed
!  OUTFN : template name for output observation arrays
!  STARTLAT : Kth entry is lowest latitude handled by Kth processor.
!  ENDLAT : Kth entry is highest latitude handled by Kth processor.
!  NPROC : number of processors.
!
!  Note: we create an output data file even if there are no observaitons
!  for a given latitude strip.
!
    use command_line_support,only:numeric_filename
    use grid_io
    implicit none
    character(*),intent(in)::outfn
    integer,intent(in)::nobs,nproc,startlat(nproc),endlat(nproc)
    type(raw_obs),intent(in)::obsbuf(nobs)
!
!  Local variables
!
    integer::k
    character(MAX_FILE_NAME_LEN)::fn
    character(4)::label='obs '
!
!  Distribute observations among processors in overlapping latitude strips.
!
    if(nproc.eq.1) then  ! no filename substitution
       fn=outfn
       call write_obs_slice(fn,'none',startlat(1),endlat(1),obsbuf,nobs)
       close(OUTNOSWAP)
    else
       do k=1,nproc
          call numeric_filename(outfn,'%n',k,fn)
          call write_obs_slice(fn,label,startlat(k),endlat(k),obsbuf,nobs)
          close(OUTNOSWAP)
       enddo
    endif
    return
    end subroutine flush
!-----------------------------
    end module readobs
