Commit b895297b authored by ronesy's avatar ronesy
Browse files

Update to loop over receptors and months

parent c00264e4
......@@ -13,6 +13,15 @@ path_flexpart: /mypath/
# Path where to write ncdf files
path_output: /mypath/
# Receptor list file
file_recept: /myfile.txt
# Run dates: format yyyymmdd
# datei = start date
# datef = end date
datei: 20140701
datef: 20140831
# Nested output (true/false)
lnested: .true.
......
......@@ -23,7 +23,7 @@
!!
!---------------------------------------------------------------------------------------
subroutine grid_convert(files, config)
subroutine grid_convert(files, config, nr, amonth)
use mod_var
use mod_settings
......@@ -34,28 +34,40 @@ subroutine grid_convert(files, config)
type (files_t), intent(in out) :: files
type (config_t), intent(in) :: config
integer, intent(in) :: nr
character(max_path_len) :: filename, filedates
character(max_path_len) :: filename, filedates, path_flexpart, path_output
character(len=max_path_len), dimension(:), allocatable :: filelist
character(len=8) :: adate
character(len=6) :: atime
character(len=6) :: atime, amonth
logical :: lexist
integer :: nfiles, i, ierr, n
integer :: jjjjmmdd, hhmiss
real, dimension(nxgrid,nygrid,ngrid) :: grid
real(kind=8), dimension(ngrid) :: gtime, time
real, dimension(:,:,:), allocatable :: grid
real(kind=8), dimension(:), allocatable :: gtime, time
real(kind=8) :: jdate
integer :: numpoint
integer :: ibdate, ibtime
integer, dimension(maxpoint) :: releases
real :: bndx, bndy, delx, dely
integer :: numx, numy
! list all grid files in directory
path_flexpart = trim(files%path_flexpart)//trim(recname(nr))//'/'//trim(amonth)//'/'
path_output = trim(files%path_output)//trim(recname(nr))//'/'
print*, path_flexpart
print*, path_output
call system('mkdir '//trim(path_output))
if ( config%lnested ) then
call system('ls '//trim(files%path_flexpart)//' | grep grid_time_nest | wc -l > '//trim(files%path_output)//'files.txt')
call system('ls '//trim(files%path_flexpart)//' | grep grid_time_nest >> '//trim(files%path_output)//'files.txt')
call system('ls '//trim(path_flexpart)//' | grep grid_time_nest_ | wc -l > '//trim(path_output)//'files.txt')
call system('ls '//trim(path_flexpart)//' | grep grid_time_nest_ >> '//trim(path_output)//'files.txt')
else
call system('ls '//trim(files%path_flexpart)//' -I "*nest*" | grep grid_time | wc -l > '//trim(files%path_output)//'files.txt')
call system('ls '//trim(files%path_flexpart)//' -I "*nest*" | grep grid_time >> '//trim(files%path_output)//'files.txt')
call system('ls '//trim(path_flexpart)//' -I "*nest*" | grep grid_time_ | wc -l > '//trim(path_output)//'files.txt')
call system('ls '//trim(path_flexpart)//' -I "*nest*" | grep grid_time_ >> '//trim(path_output)//'files.txt')
endif
! read file list
open(100,file=trim(files%path_output)//'files.txt',status="old",action="read",iostat=ierr)
open(100,file=trim(path_output)//'files.txt',status="old",action="read",iostat=ierr)
if ( ierr.ne.0 ) then
print*, 'ERROR: problem read files.txt'
stop
......@@ -68,8 +80,32 @@ subroutine grid_convert(files, config)
end do
close(100)
! read header for this month to get ngrid
if ( config%lnested ) then
filename = trim(path_flexpart)//'header_nest'
else
filename = trim(path_flexpart)//'header'
endif
inquire(file=trim(filename),exist=lexist)
if ( .not.lexist ) then
print*, 'ERROR: cannot find flexpart header'
stop
endif
print*, 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely)
if ( config%linversionout ) then
ngrid = trajdays*ndgrid+2
else
ngrid = numpoint
endif
print*, 'ngrid = ',ngrid
allocate( grid(nxgrid,nygrid,ngrid) )
allocate( gtime(ngrid) )
allocate( time(ngrid) )
! loop over files
filedates = trim(files%path_flexpart)//'dates'
filedates = trim(path_flexpart)//'dates'
do i = 1, nfiles
! read binary grid files
if ( config%lnested ) then
......@@ -79,25 +115,25 @@ subroutine grid_convert(files, config)
adate = filelist(i)(11:18)
atime = filelist(i)(19:24)
endif
print*, 'adate = ',adate
print*, 'atime = ',atime
! print*, 'adate = ',adate
! print*, 'atime = ',atime
read(adate,*) jjjjmmdd
read(atime,*) hhmiss
print*, 'jjjjmmdd = ',jjjjmmdd
print*, 'hhmiss = ',hhmiss
! print*, 'jjjjmmdd = ',jjjjmmdd
! print*, 'hhmiss = ',hhmiss
if ( config%linversionout ) then
jdate = juldate(jjjjmmdd, hhmiss)
else
jdate = juldate(datei,timei)
endif
print*, 'jdate = ',jdate
filename = trim(files%path_flexpart)//trim(filelist(i))
! print*, 'jdate = ',jdate
filename = trim(path_flexpart)//trim(filelist(i))
if ( config%linversionout ) then
call read_grid(config, filename, filedates, jdate, nxgrid, nygrid, grid, gtime)
else
call read_grid_std(config, filename, jdate, nxgrid, nygrid, grid, gtime)
endif
print*, 'gtime = ',gtime
! print*, 'gtime = ',gtime
! write to ncdf
do n = 1, ngrid
call caldate(gtime(n), jjjjmmdd, hhmiss)
......@@ -105,13 +141,16 @@ subroutine grid_convert(files, config)
end do
print*, 'time = ',time
if ( config%lnested ) then
filename = trim(files%path_output)//'grid_time_nest_'//adate//atime//'_001.nc'
filename = trim(path_output)//'grid_time_nest_'//adate//atime//'_001.nc'
else
filename = trim(files%path_output)//'grid_time_'//adate//atime//'_001.nc'
filename = trim(path_output)//'grid_time_'//adate//atime//'_001.nc'
endif
call write_ncdf(filename, time, grid)
end do
deallocate( grid )
deallocate( gtime )
deallocate( time )
end subroutine grid_convert
......@@ -27,6 +27,7 @@ program grid_to_ncdf
use mod_flexpart
use mod_var
use mod_dates
use mod_settings
implicit none
......@@ -36,6 +37,10 @@ program grid_to_ncdf
type ( config_t ) :: config
character(max_path_len) :: settings_files
character(len=6) :: adate
integer :: nr
integer :: yyyymmdd, hhmiss, eomday
real(kind=8) :: juldatei, juldatef, jd
call getarg(1,settings_files)
if (settings_files == '') then
......@@ -47,9 +52,20 @@ program grid_to_ncdf
! initialize parameters
call initialize(files, config)
juldatei = juldate(config%datei, 0)
juldatef = juldate(config%datef, 0)
! convert grid files
call grid_convert(files, config)
do nr = 1, nrec
jd = juldatei
do while ( jd.le.juldatef )
call caldate(jd, yyyymmdd, hhmiss)
write(adate,'(I6.6)') yyyymmdd/100
call grid_convert(files, config, nr, adate)
eomday = calceomday(yyyymmdd/100)
jd = jd + dble(eomday)
end do
end do
end program
......
......@@ -36,6 +36,8 @@ subroutine initialize(files, config)
character(max_path_len) :: filename
logical :: lexist
integer :: yyyymm
character(len=6) :: adate
integer :: numpoint
integer :: ibdate, ibtime
integer, dimension(maxpoint) :: releases
......@@ -43,19 +45,30 @@ subroutine initialize(files, config)
integer :: numx, numy
integer :: i, ierr
! read list of receptors
! ----------------------
filename = trim(files%file_recept)
call read_reclist(filename)
! initialize flexpart variables
! -----------------------------
if ( .not.config%lnested ) filename = trim(files%path_flexpart)//'header'
if ( config%lnested ) filename = trim(files%path_flexpart)//'header_nest'
print*, filename
inquire(file=trim(filename),exist=lexist)
yyyymm=config%datei/100
write(adate,'(i6)') yyyymm
lexist = .false.
i = 1
do while ( (.not.lexist).and.(i.le.nrec) )
if ( .not.config%lnested ) filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header'
if ( config%lnested ) filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header_nest'
inquire(file=trim(filename),exist=lexist)
i = i + 1
end do
if ( .not.lexist ) then
print*, 'ERROR initialize: cannot find flexpart header'
write(*,*) 'ERROR initialize: cannot find flexpart header'
stop
endif
print*, 'Reading flexpart header: '//trim(filename)
write(*,*) 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely)
nxgrid = numx
......
......@@ -18,6 +18,7 @@ SRCS = mod_var.f90 \
initialize.f90 \
grid_convert.f90 \
write_ncdf.f90 \
read_reclist.f90 \
grid_to_ncdf.f90
OBJECTS = $(SRCS:.f90=.o)
......
......@@ -36,6 +36,7 @@ module mod_settings
character(len=max_path_len) :: path_output ! path to output
character(len=max_path_len) :: path_flexpart ! path to flexpart output
character(len=max_path_len) :: file_recept ! receptor list file
end type files_t
......@@ -43,6 +44,8 @@ module mod_settings
type :: config_t
integer :: datei ! start date (yyyymmdd)
integer :: datef ! end date (yyyymmdd)
logical :: lnested ! use nested flexpart output (true or false)
logical :: linversionout ! inversion formatted flexpart output (true or false)
......@@ -274,8 +277,17 @@ module mod_settings
identifier = "path_output:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) files%path_output = cc
identifier = "file_recept:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) files%file_recept = cc
! read config settings
identifier = "datei:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%datei = cn
identifier = "datef:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%datef = cn
identifier = "lnested:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%lnested = cl
......
......@@ -35,6 +35,7 @@ module mod_var
integer, parameter :: maxspec=10 ! max number of species in a flexpart run
integer, parameter :: maxlev=50 ! max number of vertical levels in flexpart output
integer, parameter :: maxobs=10000 ! max number of observations
integer, parameter :: recname_len=7 ! length of receptor names
! flexpart variables
! ------------------
......@@ -58,5 +59,7 @@ module mod_var
real(kind=8), dimension(:), allocatable :: reltime ! release times
integer :: datei ! start date of run
integer :: timei ! start time of run
integer :: nrec ! number of receptors
character(recname_len), dimension(:), allocatable :: recname ! names of receptors
end module mod_var
!---------------------------------------------------------------------------------------
! GRID_TO_NCDF: read_reclist
!---------------------------------------------------------------------------------------
! 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_reclist
!!
!---------------------------------------------------------------------------------------
subroutine read_reclist(filename)
use mod_var
implicit none
character(len=max_path_len), intent(in) :: filename
character(len=200) :: line
integer :: ierr
integer :: cnt
! count number of receptors
open(100,file=trim(filename),action='read',status='old',iostat=ierr)
if(ierr.gt.0) then
write(*,*) 'ERROR: cannot open: '//trim(filename)
stop
endif
write(*,*) 'Reading receptors file: '//trim(filename)
cnt = 0
do while (ierr.eq.0)
read(100,*,iostat=ierr,end=10) line
cnt = cnt + 1
end do
10 continue
close(100)
nrec = cnt
write(*,*) 'Number of receptors: ',nrec
! read receptors
allocate ( recname(nrec) )
open(100,file=trim(filename),action='read',status='old',iostat=ierr)
write(*,*) 'Receptors: '
do cnt = 1, nrec
read(100,*) recname(cnt)
write(*,*) recname(cnt)
end do
close(100)
end subroutine read_reclist
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment