!---------------------------------------------------------------------------------------
! PREP_REGIONS: read_obs
!---------------------------------------------------------------------------------------
!  FLEXINVERT is free software: you can redistribute it and/or modify
!  it under the terms of the GNU General Public License as published by
!  the Free Software Foundation, either version 3 of the License, or
!  (at your option) any later version.
!
!  FLEXINVERT is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
!  GNU General Public License for more details.
!
!  You should have received a copy of the GNU General Public License
!  along with FLEXINVERT.  If not, see <http://www.gnu.org/licenses/>.
!
!  Copyright 2017, Rona Thompson
!---------------------------------------------------------------------------------------
!
!> read_obs
!!
!! Purpose:    Reads observations to match with FLEXPART grid_time files.
!!
!! Interface:
!!
!!    Inputs
!!             files    -  file data structure
!!             reclist  -  list of receptor IDs
!!             obstimes -  list of observation time stamps as julian days
!!             avetimes -  list of observation averaging times as julian days
!!
!---------------------------------------------------------------------------------------

subroutine read_obs(files, reclist, obstimes, avetimes)

  use mod_settings
  use mod_strings
  use mod_var
  use mod_dates

  implicit none

  type (files_t),                  intent(in)            :: files
  real(kind=8), dimension(maxobs), intent(in out)        :: obstimes
  real(kind=8), dimension(maxobs), intent(in out)        :: avetimes
  character(4), dimension(maxobs), intent(in out)        :: reclist

  character(len=max_path_len), dimension(:), allocatable :: filelist
  character(len=7)                                       :: recs
  character(len=200)                                     :: header
  character(len=200), dimension(18)                      :: args
  real, dimension(16)                                    :: temp
  integer                                                :: ierr
  integer                                                :: cnt
  integer                                                :: narg, n, reclen
  integer                                                :: nf, nfiles, nr
  integer                                                :: yyyymmdd, hhmiss, yyyy, mm, dd, hh, mi, ss, hl
  real(kind=8)                                           :: jdate, avetime
  real                                                   :: conc, err, lon, alt
  integer                                                :: num

  ! list observation files
  ! ----------------------

  print*, trim(files%path_obs)
  print*, trim(files%suffix)
  print*, trim(files%path_output)

  call system('ls '//trim(files%path_obs)//'*'//trim(files%suffix)//' | wc -l > '//trim(files%path_output)//'obsfiles.txt')
  call system('ls '//trim(files%path_obs)//' | grep '//trim(files%suffix)//' >> '//trim(files%path_output)//'obsfiles.txt')

  open(100,file=trim(files%path_output)//'obsfiles.txt',action='read',status='old',iostat=ierr)
  if ( ierr.ne. 0 ) then
    write(*,*) 'ERROR: cannot open obsfiles.txt'
    stop
  endif
  read(100,*,iostat=ierr) nfiles
  allocate ( filelist(nfiles), stat = ierr )
  do nf = 1, nfiles
    read(100,fmt='(A)',iostat=ierr) filelist(nf)
    if (ierr.ne.0) exit
  end do
  close(100)

  ! loop over files
  ! ---------------

  cnt = 0
  do nf = 1, nfiles

    ! check this file belongs to receptor in reclist
    ! added option for 4 character length names
    if (filelist(nf)(4:4).eq.'_' )then
      recs = filelist(nf)(1:3)
      reclen=3
    else
      recs = filelist(nf)(1:4)
      reclen=4
    endif
    print*, 'read_obs: rec = ',recs
    if ( .not.( any(recname(:)(1:4).eq.to_upper(recs)) ) .and. .not.( any(recname(:)(1:3).eq.to_upper(recs)) )) go to 10

    ! open input file
    open(100,file=trim(files%path_obs)//trim(filelist(nf)),action='read',status='old',iostat=ierr)
    if ( ierr.ne.0 ) then
      write(*,*) 'WARNING: cannot open: '//trim(files%path_obs)//trim(filelist(nf))
      go to 10
    endif 
    write(*,*) 'Reading file: '//trim(files%path_obs)//trim(filelist(nf))

    ! station coordinates for data selection
    do nr = 1, nrec
      if ( to_upper(recs).eq.recname(nr) ) then
        lon = reclon(nr)
        alt = recalt(nr)
        go to 20
      endif
    end do
20 continue

    ! read header
    read (100, fmt='(A)', iostat=ierr) header
    print*, 'read_obs: header = ',header

    ! read data
    read_loop: do 
      ! if obs file contains avetime
      if ( reclen.eq.4 ) then
        read(100,fmt='(A4,1X,I8,1X,I6,1X,F14.6,1X,F10.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
                recs, yyyymmdd, hhmiss, jdate, avetime, conc, err, num
      else
        read(100,fmt='(A3,1X,I8,1X,I6,1X,F14.6,1X,F10.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
                recs, yyyymmdd, hhmiss, jdate, avetime, conc, err, num
      endif
     ! if obs file doesn't contain avetime (backwards compatability)
!      avetime = 0d0
!      if ( reclen.eq.4 ) then
!        read(100,fmt='(A4,1X,I8,1X,I6,1X,F14.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
!                recs, yyyymmdd, hhmiss, jdate, conc, err, num
!      else
!        read(100,fmt='(A3,1X,I8,1X,I6,1X,F14.6,1X,F10.4,1X,F10.4,1X,I4)',iostat=ierr) &
!                recs, yyyymmdd, hhmiss, jdate, conc, err, num
!      endif
      if ( ierr.gt.0 ) exit read_loop
      if ( jdate.ge.(juldatef+1d0) ) exit read_loop
      if ( jdate.lt.juldatei ) cycle read_loop
      if ( conc.le.-999. ) cycle read_loop
      ! select day/night for low/high alt sites
      hh = hhmiss/10000
      hl = hh + int(lon*24./360.)
      if ( hl.lt.0 ) then
        hl = hl + 24
      else if ( hl.ge.24 ) then
        hl = hl - 24
      endif
      if ( alt.le.1000. ) then
        if ( (hl.lt.11).or.(hl.gt.15) ) cycle read_loop
      else
        if ( (hl.gt.3).and.(hl.lt.23) ) cycle read_loop
      endif
      if ( conc.le.-999. ) cycle read_loop

      cnt = cnt + 1
      obstimes(cnt) = jdate
      avetimes(cnt) = avetime
      reclist(cnt) = recs
    end do read_loop

    ! close input file
    close(100)

10  continue

  end do

  print*, 'read_obs: cnt = ',cnt

  ! count observations
  nobs = nobs + cnt

end subroutine read_obs


