Commit 6b9570df authored by ronesy's avatar ronesy
Browse files

Updates to prep_flexpart to set-up FP runs with regular releases

parent 7d702dfd
......@@ -31,6 +31,11 @@ file_avail: /mypath_to_available_file/AVAILABLE
file_availnest: /mypath_to_available_file/AVAILABLE
# General settings
# make flexpart releases only when there are observations (1) or at regular intervals (0)
# important: if set lrelease = 0 observations have to be from fixed locations (not aircraft, ship, satellite)
lrelease: 0
# if regular flexpart releases (lrelease = 0) set frequency (hours)
relfreq: 1
# select afternoon observations (night for mountains) (select = 1, no selection = 0)
lselect: 1
# average observations (average = 1, no average = 0)
......
......@@ -31,6 +31,11 @@ file_avail: /mypath_to_available_file/AVAILABLE
file_availnest: /mypath_to_available_file/AVAILABLE
# General settings
# make flexpart releases only when there are observations (1) or at regular intervals (0)
# important: if set lrelease = 0 observations have to be from fixed locations (not aircraft, ship, satellite)
lrelease: 0
# if regular flexpart releases (lrelease = 0) set frequency (hours)
relfreq: 1
# select afternoon observations (night for mountains) (select = 1, no selection = 0)
lselect: 1
# average observations (average = 1, no average = 0)
......
......@@ -31,6 +31,11 @@ file_avail: /mypath_to_available_file/AVAILABLE
file_availnest: /mypath_to_available_file/AVAILABLE
# General settings
# make flexpart releases only when there are observations (1) or at regular intervals (0)
# important: if set lrelease = 0 observations have to be from fixed locations (not aircraft, ship, satellite)
lrelease: 0
# if regular flexpart releases (lrelease = 0) set frequency (hours)
relfreq: 1
# select afternoon observations (night for mountains) (select = 1, no selection = 0)
lselect: 0
# average observations (average = 1, no average = 0)
......
......@@ -20,7 +20,7 @@ path_ohfield: /mypath_to_ohfields/
path_obs: /mypath/TEST_INPUT/OBS/GHG/
# path where to write observation output
path_obsout: /mypath/TEST_OUTPUT/OBS/GHG/
# observation file format (one of obspack, wdcgg, noaa, basic, icos)
# observation file format (one of obspack, wdcgg, noaa, icos)
obsformat: noaa
# observation file suffix
suffix: event.txt
......@@ -31,6 +31,11 @@ file_avail: /mypath_to_available_file/AVAILABLE
file_availnest: /mypath_to_available_file/AVAILABLE
# General settings
# make flexpart releases only when there are observations (1) or at regular intervals (0)
# important: if set lrelease = 0 observations have to be from fixed locations (not aircraft, ship, satellite)
lrelease: 0
# if regular flexpart releases (lrelease = 0) set frequency (hours)
relfreq: 1
# select afternoon observations (night for mountains) (select = 1, no selection = 0)
lselect: 0
# average observations (average = 1, no average = 0)
......
......@@ -36,6 +36,14 @@ subroutine list_obsfiles(settings)
! list observation files
! ----------------------
call system('rm -f '//trim(settings%path_obsout)//'obsfiles.txt')
! use find to descend into subdirectories (filelist then contains file path and requires editting read_xxx.f90)
! call system ('find ' //trim(settings%path_obs)// ' -not -path "*/\.*" -type f \( -iname "' //'*'//trim(settings%suffix)// &
! '" \) | wc -l > '//trim(settings%path_obsout)//'obsfiles.txt')
! call system ('find ' //trim(settings%path_obs)// ' -not -path "*/\.*" -type f \( -iname "' //'*'//trim(settings%suffix)// &
! '" \) >> '//trim(settings%path_obsout)//'obsfiles.txt')
call system('rm -f '//trim(settings%path_obsout)//'obsfiles.txt')
call system('ls '//trim(settings%path_obs)//'*'//trim(settings%suffix)//' | wc -l > '//trim(settings%path_obsout)//'obsfiles.txt')
call system('ls '//trim(settings%path_obs)//' | grep '//trim(settings%suffix)//' >> '//trim(settings%path_obsout)//'obsfiles.txt')
......
......@@ -97,6 +97,11 @@ program main
print*, 'WARNING: laverage is 0 but intaverage > 0'
print*, 'setting intaverage to 0'
endif
if( settings%lrelease.eq.0 ) then
print*, 'FLEXPART releases at regular intervals'
else
print*, 'FLEXPART releases at observation timestamps'
endif
! read receptor list
call read_reclist(settings%file_recept)
......@@ -139,11 +144,16 @@ program main
print*, 'ERROR: unknown observation file format'
endif
! no observations
if ( nobs.eq.0 ) go to 10
! write RELEASES
call prep_releases(settings, jd, nr, nobs, obs)
if( settings%lrelease.eq.0 ) then
! flexpart releases at regular intervals
call prep_releases_reg(settings, jd, nr)
else
! flexpart releases for observation timestamps
! if no observations don't create releases file
if ( nobs.eq.0 ) go to 10
call prep_releases(settings, jd, nr, nobs, obs)
endif
10 continue
......
......@@ -24,12 +24,13 @@ SRCS = mod_var.f90 \
prep_outgrid.f90 \
prep_ageclass.f90 \
prep_releases.f90 \
prep_releases_reg.f90 \
process_obs.f90 \
read_obspack.f90 \
read_wdcgg.f90 \
read_noaa.f90 \
read_basic.f90 \
read_icos.f90 \
read_basic.f90 \
main.f90
......
......@@ -49,6 +49,8 @@ module mod_settings
character(len=max_name_len) :: suffix
integer :: lselect
integer :: lrelease
real :: relfreq
integer :: laverage
real :: intaverage
integer :: lnamelist
......@@ -347,6 +349,12 @@ module mod_settings
identifier = "lselect:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%lselect = int(cn)
identifier = "lrelease:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%lrelease = int(cn)
identifier = "relfreq:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%relfreq = real(cn)
identifier = "laverage:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%laverage = int(cn)
......
......@@ -44,30 +44,38 @@ module mod_tools
!!
! --------------------------------------------------
subroutine sort(n, arr)
subroutine sort(n, arr, ind)
implicit none
integer, intent(in) :: n
real, dimension(n), intent(in out) :: arr
integer, intent(in) :: n
real(kind=8), dimension(n), intent(in out) :: arr
integer, dimension(n), intent(in out) :: ind
integer, parameter :: m=7, nstack=50
integer :: i, ir, j, k, l, jstack
integer, dimension(nstack) :: istack
real :: a, temp
integer, parameter :: m=7, nstack=50
integer :: i, ir, j, k, l, jstack
integer, dimension(nstack) :: istack
real(kind=8) :: a, temp
integer :: ia, itemp
do i=1,n
ind(i)=i
end do
jstack=0
l=1
ir=n
1 if(ir-l.lt.m) then
do j=l+1,ir
a=arr(j)
ia=ind(j)
do i=j-1,1,-1
if(arr(i).le.a) go to 2
arr(i+1)=arr(i)
ind(i+1)=ind(i)
end do
i=0
2 arr(i+1)=a
ind(i+1)=ia
end do
if(jstack.eq.0) return
ir=istack(jstack)
......@@ -76,26 +84,39 @@ module mod_tools
else
k=(l+ir)/2
temp=arr(k)
itemp=ind(k)
arr(k)=arr(l+1)
arr(l+1)=temp
ind(k)=ind(l+1)
ind(l+1)=itemp
if(arr(l+1).gt.arr(ir)) then
temp=arr(l+1)
itemp=ind(l+1)
arr(l+1)=arr(ir)
arr(ir)=temp
ind(l+1)=ind(ir)
ind(ir)=itemp
endif
if(arr(l).gt.arr(ir)) then
temp=arr(l)
itemp=ind(l)
arr(l)=arr(ir)
arr(ir)=temp
ind(l)=ind(ir)
ind(ir)=itemp
endif
if(arr(l+1).gt.arr(l)) then
temp=arr(l+1)
itemp=ind(l+1)
arr(l+1)=arr(l)
arr(l)=temp
ind(l+1)=ind(l)
ind(l)=itemp
endif
i=l+1
j=ir
a=arr(l)
ia=ind(l)
3 continue
i=i+1
if(arr(i).lt.a) go to 3
......@@ -104,11 +125,16 @@ module mod_tools
if(arr(j).gt.a) go to 4
if(j.lt.i) go to 5
temp=arr(i)
itemp=ind(i)
arr(i)=arr(j)
arr(j)=temp
ind(i)=ind(j)
ind(j)=itemp
go to 3
5 arr(l)=arr(j)
arr(j)=a
ind(l)=ind(j)
ind(j)=ia
jstack=jstack+2
if(jstack.gt.nstack) print*, 'ERROR sort: nstack too small'
if(ir-i+1.ge.j-l) then
......
......@@ -36,7 +36,8 @@ module mod_var
integer :: nrec
integer :: nfiles
character(len=3), dimension(:), allocatable :: recname
character(len=100), dimension(:), allocatable :: filelist
character(len=256), dimension(:), allocatable :: filelist
real, dimension(:), allocatable :: reclon, reclat, recalt
! flexpart variables
! ------------------
......
!---------------------------------------------------------------------------------------
! PREP_FLEXPART: prep_releases_reg
!---------------------------------------------------------------------------------------
! 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
!---------------------------------------------------------------------------------------
!
!> prep_releases_reg
!!
!! Purpose: Prepares the options file RELEASES for regular release times
!!
!! Interface:
!!
!! Inputs
!! settings - settings data structure
!! jd - julian day for start of month
!! nr - index to the receptors list
!!
!! Externals
!! caldate
!! skiplines
!!
!! Note: For releases made at regular intervals (and not matching the
!! observation timestamps) the lat,lon, alt coordinates are stationary.
!! For data from a moving platform (aircraft, ship, or satellite) need to
!! use prep_releases.f90 (so lrelease = 1)
!!
!---------------------------------------------------------------------------------------
subroutine prep_releases_reg(settings, jd, nr)
use mod_var
use mod_settings
use mod_dates
use mod_tools
implicit none
type (settings_t), intent(in) :: settings
real(kind=8), intent(in) :: jd
integer, intent(in) :: nr
character(len=max_path_len) :: filename, pathname
character(len=100) :: line
character(len=6) :: adate
character(len=3) :: acnt, aspec
integer :: jjjjmmdd, hhmiss, ss
integer :: ierr
integer :: spec
integer :: nchar, n
integer :: eomday
real(kind=8) :: jdi, jdf, jdi_start, jdi_end, relfreq
integer :: nspec, specnum_rel
integer :: idate1, itime1
integer :: idate2, itime2
real :: lon1, lon2
real :: lat1, lat2
real :: z1, z2
integer :: zkind, parts
real :: mass
character(len=40) :: comment
namelist /releases_ctrl/ &
nspec, &
specnum_rel
namelist /release/ &
idate1, itime1, &
idate2, itime2, &
lon1, lon2, &
lat1, lat2, &
z1, z2, &
zkind, &
mass, &
parts, &
comment
! initialize date/time variables
call caldate(jd, jjjjmmdd, hhmiss)
write(adate,fmt='(I6.6)') jjjjmmdd/100
eomday = calceomday(jjjjmmdd/100)
jdf = jd + real(eomday,kind=8)
jdi = jd
relfreq = dble(settings%relfreq)/24d0
print*, 'prep_releases_reg: relfreq = ',relfreq
print*, 'prep_releases_reg: jdi = ',jdi
print*, 'prep_releases_reg: jdf = ',jdf
! read species file
open(100,file=trim(settings%path_flexpart)//'options/SPECIES/spec_overview',status='old',action='read',iostat=ierr)
nchar = len_trim(settings%species)
do while ( ierr.eq.0 )
read (100, fmt='(A)', iostat=ierr) line
if ( line(len_trim(line)-nchar+1:len_trim(line)) == trim(settings%species) ) ierr = 1
end do
read(line(9:11),*) spec
write(aspec,fmt='(I3.3)') spec
! preset namelist variables releases_ctrl
nspec = 1
specnum_rel = spec
! copy standard input files to options folder
pathname = trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/'
print*, 'prep_releases_reg: pathname = ',pathname
print*, 'prep_releases_reg: pathflexpart = ',settings%path_flexpart
call system('mkdir -p '//trim(pathname)//'SPECIES/')
filename = trim(settings%path_flexpart)//'options/SPECIES/SPECIES_'//aspec
print*, 'prep_releases_reg: filename = ',filename
call system('cp '//trim(filename)//' '//trim(pathname)//'SPECIES/')
call system('cp '//trim(settings%path_flexpart)//'options/*.dat'//' '//trim(pathname))
call system('cp '//trim(settings%path_flexpart)//'options/*.t'//' '//trim(pathname))
! open RELEASES file
filename = trim(pathname)//'RELEASES'
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
! write RELEASES file
if( settings%lnamelist.eq.1 ) then
! use namelist format
write(100,nml=releases_ctrl)
n = 0
do while ( jdi.le.(jdf-relfreq) )
n = n + 1
! set namelist variables releases
jdi_start = dnint(jdi*1.e6)
jdi_start = jdi_start/1.e6
call caldate(jdi_start, jjjjmmdd, hhmiss)
idate1 = jjjjmmdd
itime1 = hhmiss
jdi_end = dnint((jdi+relfreq)*1.e6)
jdi_end = jdi_end/1.e6
call caldate(jdi_end, jjjjmmdd, hhmiss)
idate2 = jjjjmmdd
itime2 = hhmiss
lon1 = reclon(nr)
lon2 = reclon(nr)
lat1 = reclat(nr)
lat2 = reclat(nr)
z1 = recalt(nr)
z2 = recalt(nr)
zkind = settings%zref
mass = 1000.
parts = settings%npart
write(acnt,fmt='(I3.3)') n
comment = trim(recname(nr))//'_'//adate//'_'//acnt
write(100,nml=release)
jdi = jdi + relfreq
end do
else
! use old file format
write(100,fmt='(A)') '*************************************************************************'
write(100,fmt='(A)') '* *'
write(100,fmt='(A)') '* *'
write(100,fmt='(A)') '* *'
write(100,fmt='(A)') '* *'
write(100,fmt='(A)') '* Input file for the Lagrangian particle dispersion model FLEXPART *'
write(100,fmt='(A)') '* *'
write(100,fmt='(A)') '* *'
write(100,fmt='(A)') '* *'
write(100,fmt='(A)') '* *'
write(100,fmt='(A)') '*************************************************************************'
call skiplines(100,1)
write(100,fmt='(I3)') 1
write(100,fmt='(I3)') spec
write(100,fmt='(A)') '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
n = 0
do while ( jdi.le.(jdf-relfreq) )
n = n + 1
jdi_start = dnint(jdi*1.e6)
jdi_start = jdi_start/1.e6
call caldate(jdi_start, jjjjmmdd, hhmiss)
write(100,fmt='(I8,1X,I6)') jjjjmmdd,hhmiss
jdi_end = dnint((jdi+relfreq)*1.e6)
jdi_end = jdi_end/1.e6
call caldate(jdi_end, jjjjmmdd, hhmiss)
write(100,fmt='(I8,1X,I6)') jjjjmmdd,hhmiss
write(100,fmt='(F9.4)') reclon(nr)
write(100,fmt='(F9.4)') reclat(nr)
write(100,fmt='(F9.4)') reclon(nr)
write(100,fmt='(F9.4)') reclat(nr)
write(100,fmt='(I9)') settings%zref
write(100,fmt='(F10.3)') recalt(nr)
write(100,fmt='(F10.3)') recalt(nr)
write(100,fmt='(I9)') settings%npart
write(100,fmt='(F9.4)') mass
write(acnt,fmt='(I3.3)') n
write(100,fmt='(A)') trim(recname(nr))//'_'//adate//'_'//acnt
write(100,fmt='(A)') '+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++'
jdi = jdi + relfreq
end do
endif
close(100)
end subroutine prep_releases_reg
......@@ -61,22 +61,45 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
integer, intent(in out) :: nobs
real, intent(in) :: freq
real(kind=8), dimension(maxobs), intent(in) :: jdobs
real, dimension(maxobs), intent(in) :: latobs, lonobs, altobs, concobs, errobs
real, dimension(maxobs), intent(in out) :: latobs, lonobs, altobs, concobs, errobs
type (obs_t), intent(in out) :: obs
integer :: cnt
integer :: c, i, m, n
integer :: c, i, m, n, nprev
integer :: jjjjmmdd, hhmiss, yyyy, mm, dd, hh, mi, ss
integer :: eomday, hloc
real :: conc, err
real(kind=8) :: jdate, jdatestart, jdateend, jdt, jdsom, jdeom
real(kind=8), dimension(nobs,6) :: buffer
real, dimension(nobs) :: work
real(kind=8), dimension(nobs) :: jdobs_out
real(kind=8), dimension(nobs) :: work
real, dimension(nobs) :: concobs_out, latobs_out, lonobs_out, altobs_out, errobs_out
integer, dimension(nobs) :: ind
! check that time series are increasing
! if two different time series were patched together they need
! to be sorted before averaging can be done
jdobs_out = jdobs(1:nobs)
call sort(nobs, jdobs_out, ind)
do i = 1, nobs
concobs_out(ind(i)) = concobs(i)
latobs_out(ind(i)) = latobs(i)
lonobs_out(ind(i)) = lonobs(i)
altobs_out(ind(i)) = altobs(i)
errobs_out(ind(i)) = errobs(i)
end do
concobs(1:nobs) = concobs_out
latobs(1:nobs) = latobs_out
lonobs(1:nobs) = lonobs_out
altobs(1:nobs) = altobs_out
errobs(1:nobs) = errobs_out
print*, 'process_obs: ind = ',ind
! if average and/or select obs
if( settings%laverage.eq.1.or.settings%lselect.eq.1 ) then
n = 1
nprev = 0
i = 0
jdate = jdobs(n)
call caldate(jdate, jjjjmmdd, hhmiss)
......@@ -87,24 +110,32 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
print*, 'jdeom: ',jdeom
! loop over all observations in current month
do while( jdate.le.jdobs(nobs) )
jdt = dint(jdate*1.e6)/1.e6
if ( settings%intaverage.ge.6. ) then
if ( n.eq.1 ) then
! if long averaging interval start at 00:00
jdatestart = dint(jdate)
jdt = dnint(jdate*1.e6)/1.e6
if ( n.ne.nprev ) then
if ( settings%intaverage.ge.6. ) then
if ( n.eq.1 ) then
! if long averaging interval start at 00:00
jdatestart = dint(jdate)
else
! continue regular intervals (start at end of last interval)
jdatestart = jdateend
endif
else
! continue regular intervals (start at end of last interval)
jdatestart = jdateend
! otherwise start at current timestamp
jdatestart = jdt
endif
else
! otherwise start at current timestamp
jdatestart = jdt
endif
if( n.eq.nprev .and. settings%intaverage.ge.6. ) then
! if no observation in last averaging interval (i.e. did not enter while loop below)
! increment start date
print*, 'incrementing start date'
jdatestart = dnint((jdatestart+real(settings%intaverage/24.,kind=8))*1.e6)/1.e6
endif
if ( settings%intaverage.gt.0. ) then
! averaging interval extends up to jdateend
jdateend = dint((jdatestart+real(settings%intaverage/24.,kind=8))*1.e6)/1.e6
jdateend = dnint((jdatestart+real(settings%intaverage/24.,kind=8))*1.e6)/1.e6
else
jdateend = dint((jdatestart+real(freq/24.,kind=8))*1.e6)/1.e6
jdateend = dnint((jdatestart+real(freq/24.,kind=8))*1.e6)/1.e6
endif
! if interval crosses month set to min of current month
jdateend = min(jdateend,jdeom)
......@@ -160,12 +191,15 @@ subroutine process_obs(settings, nobs, freq, jdobs, latobs, lonobs, altobs, conc
buffer(cnt,5) = concobs(n)
buffer(cnt,6) = errobs(n)
endif
nprev = n
n = n + 1
if ( n.gt.nobs ) go to 10
print*, 'jdt: ',jdt
jdt = jdobs(n)
! jdt = jdobs(n)
jdt = dnint(jdobs(n)*1.e6)/1.e6
print*, 'next jdt: ',jdt
end do