!---------------------------------------------------------------------------------------
! 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