Skip to content
Snippets Groups Projects
Commit 15f86f38 authored by ronesy's avatar ronesy
Browse files

changes for FPv11 with backwards compatibility with FPv10.4 and name change for aerosol species

parent 062adde1
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Showing
with 1002 additions and 109 deletions
grid_to_ncdf/grid_convert.f90
+ 9
− 3
View file @ 15f86f38
...@@ -51,7 +51,7 @@ subroutine grid_convert(files, config, nr, timestamp) ...@@ -51,7 +51,7 @@ subroutine grid_convert(files, config, nr, timestamp)
integer :: ibdate, ibtime integer :: ibdate, ibtime
integer, dimension(maxpoint) :: releases integer, dimension(maxpoint) :: releases
real :: bndx, bndy, delx, dely real :: bndx, bndy, delx, dely
integer :: numx, numy integer :: numx, numy, numz
! list all grid files in directory ! list all grid files in directory
if ( trim(config%filefreq).eq.'month' ) then if ( trim(config%filefreq).eq.'month' ) then
...@@ -90,17 +90,23 @@ subroutine grid_convert(files, config, nr, timestamp) ...@@ -90,17 +90,23 @@ subroutine grid_convert(files, config, nr, timestamp)
! read header for this month to get ngrid ! read header for this month to get ngrid
if ( config%lnested ) then if ( config%lnested ) then
filename = trim(path_flexpart)//'header_nest' filename = trim(path_flexpart)//'header_nest'
inquire(file=trim(filename),exist=lexist)
if ( .not.lexist ) then
! test if is FPv11 file
filename = trim(path_flexpart)//'header_nest_grid_time'
inquire(file=trim(filename),exist=lexist)
endif
else else
filename = trim(path_flexpart)//'header' filename = trim(path_flexpart)//'header'
inquire(file=trim(filename),exist=lexist)
endif endif
inquire(file=trim(filename),exist=lexist)
if ( .not.lexist ) then if ( .not.lexist ) then
print*, 'ERROR: cannot find flexpart header' print*, 'ERROR: cannot find flexpart header'
stop stop
endif endif
print*, 'Reading flexpart header: '//trim(filename) print*, 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, & call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely) numx, numy, numz, bndx, bndy, delx, dely)
if ( config%linversionout ) then if ( config%linversionout ) then
ngrid = trajdays*ndgrid+2 ngrid = trajdays*ndgrid+2
else else
......
grid_to_ncdf/init_convert.f90
+ 3
− 3
View file @ 15f86f38
...@@ -52,7 +52,7 @@ subroutine init_convert(files, config, nr, timestamp) ...@@ -52,7 +52,7 @@ subroutine init_convert(files, config, nr, timestamp)
integer :: ibdate, ibtime integer :: ibdate, ibtime
integer, dimension(maxpoint) :: releases integer, dimension(maxpoint) :: releases
real :: bndx, bndy, delx, dely real :: bndx, bndy, delx, dely
integer :: numx, numy integer :: numx, numy, numz
! list all grid files in directory ! list all grid files in directory
if ( trim(config%filefreq).eq.'month' ) then if ( trim(config%filefreq).eq.'month' ) then
...@@ -92,8 +92,8 @@ subroutine init_convert(files, config, nr, timestamp) ...@@ -92,8 +92,8 @@ subroutine init_convert(files, config, nr, timestamp)
endif endif
print*, 'Reading flexpart header: '//trim(filename) print*, 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, & call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely) numx, numy, numz, bndx, bndy, delx, dely)
print*, 'numx, numy, bndx, bndy, delx, dely = ',numx,numy,bndx,bndy,delx,dely print*, 'numx, numy, numz, bndx, bndy, delx, dely = ',numx,numy,numz,bndx,bndy,delx,dely
allocate( grid(nxgrid,nygrid,nzgrid) ) allocate( grid(nxgrid,nygrid,nzgrid) )
allocate( levels(nzgrid) ) allocate( levels(nzgrid) )
! trick to use write_ncdf also for grid_initial ! trick to use write_ncdf also for grid_initial
......
grid_to_ncdf/initialize.f90
+ 16
− 3
View file @ 15f86f38
...@@ -43,7 +43,7 @@ subroutine initialize(files, config) ...@@ -43,7 +43,7 @@ subroutine initialize(files, config)
integer :: ibdate, ibtime integer :: ibdate, ibtime
integer, dimension(maxpoint) :: releases integer, dimension(maxpoint) :: releases
real :: bndx, bndy, delx, dely real :: bndx, bndy, delx, dely
integer :: numx, numy integer :: numx, numy, numz
integer :: i, ierr integer :: i, ierr
! read list of receptors ! read list of receptors
...@@ -74,15 +74,27 @@ subroutine initialize(files, config) ...@@ -74,15 +74,27 @@ subroutine initialize(files, config)
else else
filename = trim(files%path_flexpart)//trim(areldate)//'/header' filename = trim(files%path_flexpart)//trim(areldate)//'/header'
endif endif
inquire(file=trim(filename),exist=lexist)
else if ( config%lnested ) then else if ( config%lnested ) then
if ( trim(config%filefreq).eq.'month' ) then if ( trim(config%filefreq).eq.'month' ) then
filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header_nest' filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header_nest'
inquire(file=trim(filename),exist=lexist)
if ( .not.lexist ) then
! test if is FPv11 file
filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header_nest_grid_time'
inquire(file=trim(filename),exist=lexist)
endif
else else
filename = trim(files%path_flexpart)//trim(areldate)//'/header_nest' filename = trim(files%path_flexpart)//trim(areldate)//'/header_nest'
inquire(file=trim(filename),exist=lexist)
if ( .not.lexist ) then
! test if is FPv11 file
filename = trim(files%path_flexpart)//trim(areldate)//'/header_nest_grid_time'
inquire(file=trim(filename),exist=lexist)
endif
endif endif
endif endif
print*, filename print*, filename
inquire(file=trim(filename),exist=lexist)
i = i + 1 i = i + 1
end do end do
if ( .not.lexist ) then if ( .not.lexist ) then
...@@ -91,9 +103,10 @@ subroutine initialize(files, config) ...@@ -91,9 +103,10 @@ subroutine initialize(files, config)
endif endif
write(*,*) 'Reading flexpart header: '//trim(filename) write(*,*) 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, & call read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely) numx, numy, numz, bndx, bndy, delx, dely)
nxgrid = numx nxgrid = numx
nygrid = numy nygrid = numy
nzgrid = numz
llx = bndx llx = bndx
lly = bndy lly = bndy
dx = delx dx = delx
......
grid_to_ncdf/makefile
+ 34
− 27
View file @ 15f86f38
F90 = gfortran F90 = /apps/sw/ubuntu22.04/gcc-11.4.0/gcc-13.2.0-tkesyophy2o6rjlzknndu3b4oyasvuqm/bin/gfortran
LIBPATH = /usr/lib/ LIBPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/lib/
INCPATH = /usr/include/ INCPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/include/
LNK = -o
CMPL = -c # OPTIMIZATION LEVEL
LIBS = -lnetcdf -lnetcdff -llapack O_LEV = 0 # [0,1,2,3,g,s,fast]
#FFLAGS = -O2 -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -ffree-form
FFLAGS = -O0 -g -m64 -fbounds-check -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -ffree-form \ LIBS = -lm -DUSE_NCF -lnetcdff
-fbacktrace
LDFLAGS = $(FFLAGS) -L$(LIBPATH) -I$(INCPATH) $(LIBS) FFLAGS = -I$(INCPATH1) -O$(O_LEV) -g -cpp -m64 -mcmodel=medium -fconvert=little-endian -frecord-marker=4 -fmessage-length=0 -flto=jobserver -O$(O_LEV) -DUSE_NCF -lnetcdff -fbacktrace -Warray-bounds -fcheck=all # -march=native
LDFLAGS = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS)
MAIN = grid_to_ncdf MAIN = grid_to_ncdf
SRCS = mod_var.f90 \ MODULES = mod_var.o \
mod_settings.f90 \ mod_settings.o \
mod_dates.f90 \ mod_dates.o \
mod_flexpart.f90 \ mod_flexpart.o
initialize.f90 \
grid_convert.f90 \ OBJECTS = initialize.o \
init_convert.f90 \ grid_convert.o \
factor_convert.f90 \ init_convert.o \
write_ncdf.f90 \ factor_convert.o \
read_reclist.f90 \ write_ncdf.o \
grid_to_ncdf.f90 read_reclist.o \
grid_to_ncdf.o
OBJECTS = $(SRCS:.f90=.o)
$(MAIN): $(OBJECTS) $(MODULES) %.o: %.mod
$(F90) $(LNK) $(MAIN) $(OBJECTS) $(LIBS)
$(MAIN): $(MODULES) $(OBJECTS)
+$(F90) -o $@ $(MODULES) $(OBJECTS) $(LDFLAGS)
%.o : %.f90 %.o : %.f90
$(F90) $(LDFLAGS) $(CMPL) $< -o $@ +$(F90) -c $(FFLAGS) $<
clean: clean:
rm -f $(OBJECTS) $(MODULES) rm -f *.o *.mod
.SUFFIXES = $(SUFFIXES) .f90
grid_to_ncdf/mod_flexpart.f90
+ 8
− 8
View file @ 15f86f38
...@@ -49,7 +49,7 @@ module mod_flexpart ...@@ -49,7 +49,7 @@ module mod_flexpart
! -------------------------------------------------- ! --------------------------------------------------
subroutine read_header(filename, numpoint, ibdate, ibtime, releases, & subroutine read_header(filename, numpoint, ibdate, ibtime, releases, &
numx, numy, bndx, bndy, delx, dely) numx, numy, numz, bndx, bndy, delx, dely)
use mod_var use mod_var
...@@ -59,7 +59,7 @@ module mod_flexpart ...@@ -59,7 +59,7 @@ module mod_flexpart
integer, intent(in out) :: numpoint integer, intent(in out) :: numpoint
integer, intent(in out) :: ibdate, ibtime integer, intent(in out) :: ibdate, ibtime
integer, dimension(maxpoint), intent(in out) :: releases integer, dimension(maxpoint), intent(in out) :: releases
integer, intent(in out) :: numx, numy integer, intent(in out) :: numx, numy, numz
real, intent(in out) :: bndx, bndy, delx, dely real, intent(in out) :: bndx, bndy, delx, dely
character(len=29) :: flexversion character(len=29) :: flexversion
...@@ -75,7 +75,7 @@ module mod_flexpart ...@@ -75,7 +75,7 @@ module mod_flexpart
real, dimension(:), allocatable :: xpoint, ypoint, zpoint1, zpoint2 real, dimension(:), allocatable :: xpoint, ypoint, zpoint1, zpoint2
real, dimension(:,:), allocatable :: xmass real, dimension(:,:), allocatable :: xmass
integer, dimension(10) :: lage integer, dimension(10) :: lage
integer :: n, i, j, ierr integer :: n, i, j, ierr, dummy
inquire ( file=trim(filename),exist=lexist ) inquire ( file=trim(filename),exist=lexist )
if ( .not.lexist ) then if ( .not.lexist ) then
...@@ -91,14 +91,14 @@ module mod_flexpart ...@@ -91,14 +91,14 @@ module mod_flexpart
read(100) ibdate, ibtime, flexversion read(100) ibdate, ibtime, flexversion
read(100) loutstep, loutaver, loutsample read(100) loutstep, loutaver, loutsample
read(100) bndx, bndy, numx, numy, delx, dely read(100) bndx, bndy, numx, numy, delx, dely
read(100) nzgrid, (outheight(i), i = 1, nzgrid) read(100) numz, (outheight(i), i = 1, numz)
read(100) jjjjmmdd, ihmmss read(100) jjjjmmdd, ihmmss
read(100) nspec, nzgrid read(100) nspec, numz
nspec=nspec/3 nspec=nspec/3
do n=1,nspec do n=1,nspec
read(100) nzgrid, species(n) read(100) dummy, species(n)
read(100) nzgrid, species(n) read(100) dummy, species(n)
read(100) nzgrid, species(n) read(100) numz, species(n)
end do end do
read(100) numpoint read(100) numpoint
......
grid_to_ncdf/run_grid2ncdf.sh 0 → 100755
+ 16
− 0
View file @ 15f86f38
#!/bin/bash
#SBATCH -J grid2ncdf
#SBATCH -p main
#SBATCH -o grid2ncdf.out
#SBATCH --time=2:00:00
#SBATCH --mem=8000
module purge
module load gcc/13.2.0-tkesy
module load openmpi/5.0.3-l3mqa
module load eccodes/2.34.0-wrusa
export LIBRARY_PATH=$LIBRARY_PATH:/apps/sw/ubuntu22.04/gcc-13.2.0/eccodes-2.34.0-wrusaoev7mlquymqrelfbxrhvnagtz6d/lib:/apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/lib
export LD_LIBRARY_PATH=$LD_LIBRARY_PATH:/apps/sw/ubuntu22.04/gcc-13.2.0/eccodes-2.34.0-wrusaoev7mlquymqrelfbxrhvnagtz6d/lib:/apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/lib
srun ./grid_to_ncdf SETTINGS
prep_flexpart/SETTINGS
+ 20
− 2
View file @ 15f86f38
...@@ -8,16 +8,20 @@ ...@@ -8,16 +8,20 @@
# ================================================== # ==================================================
# Path and file settings # Path and file settings
# use scratch drive, if yes writes relative paths for options and output (1 = yes, 0 = no)
lscratch : 1
# flexpart source path # flexpart source path
path_flexpart: /mypath/flexpart/ path_flexpart: /mypath/flexpart/
# flexpart version (either 10 or 11)
FPversion: 10
# path where to write options folder (root only, recname and month are appended) # path where to write options folder (root only, recname and month are appended)
path_options: /mypath/flexpart_options/ path_options: /mypath/flexpart_options/
# path where to write flexpart output (root only, recname and month are appended) # path where to write flexpart output (root only, recname and month are appended)
path_output: /mypath/flexpart_output/ path_output: /mypath/flexpart_output/
# path to OH fields (if use OH chemistry)
path_ohfield: /mypath/ohfields/
# path to observations # path to observations
path_obs: /mypath/obs/ path_obs: /mypath/obs/
# path to OH fields (FPv10 only)
path_ohfield: /mypath/oh_fields/
# path where to write observation output # path where to write observation output
path_obsout: /mypath/obs_out/ path_obsout: /mypath/obs_out/
# observation file format (one of obspack, wdcgg, noaa, basic, icos) # observation file format (one of obspack, wdcgg, noaa, basic, icos)
...@@ -29,6 +33,18 @@ file_recept: /mypath/input/reclist.txt ...@@ -29,6 +33,18 @@ file_recept: /mypath/input/reclist.txt
# AVAILABLE files # AVAILABLE files
file_avail: /mypath/available_files/AVAILABLE file_avail: /mypath/available_files/AVAILABLE
file_availnest: /mypath/available_files/AVAILABLE file_availnest: /mypath/available_files/AVAILABLE
# windfields to use (determines lon of domain for flexpart run) (oper, era5)
windfield: era5
# If species has chemical reactions (for FPv11 only)
# number of reactions
nreagent: 1
# comma-separated list of paths to reagents files (length of nreagent)
path_reagents: /mypath/oh_fields/
# comma-separated list of units in reagents files ("molecule/cm3" or "mol/mol") (length of nreagent)
reagent_units: molecule/cm3
# comma-sepated list indicating if reagent fields are to be interpolated to hourly based on SZA (0 = no, 1 = yes)
reagent_interp: 1
# General settings # General settings
# make flexpart releases only when there are observations (1) or at regular intervals (0) # make flexpart releases only when there are observations (1) or at regular intervals (0)
...@@ -44,6 +60,8 @@ laverage: 0 ...@@ -44,6 +60,8 @@ laverage: 0
intaverage: 0 intaverage: 0
# flexpart file format (namelist = 1, ascii = 0) # flexpart file format (namelist = 1, ascii = 0)
lnamelist: 1 lnamelist: 1
# if observation output already exists - append (1) or make new files (0)
lappendfile: 0
# COMMAND settings # COMMAND settings
# start date of simulation (yyyymmdd) # start date of simulation (yyyymmdd)
......
prep_flexpart/list_obsfiles.f90
+ 23
− 5
View file @ 15f86f38
...@@ -27,22 +27,40 @@ subroutine list_obsfiles(settings) ...@@ -27,22 +27,40 @@ subroutine list_obsfiles(settings)
use mod_var use mod_var
use mod_settings use mod_settings
use mod_dates
implicit none implicit none
type (settings_t), intent(in) :: settings type (settings_t), intent(in) :: settings
integer :: nf, ierr integer :: nf, ierr
character(len=6) :: adate
integer :: jjjjmmdd, hhmiss
! list observation files ! list observation files
! ---------------------- ! ----------------------
call system('rm -f '//trim(settings%path_obsout)//'obsfiles.txt') if ( settings%lappendfile.eq.1 ) then
call system('ls '//trim(settings%path_obs)//'*'//trim(settings%suffix)//' | wc -l > '//trim(settings%path_obsout)//'obsfiles.txt') call system('rm -f '//trim(settings%path_obsout)//'obsfiles.dat')
call system('ls '//trim(settings%path_obs)//' | grep '//trim(settings%suffix)//' >> '//trim(settings%path_obsout)//'obsfiles.txt') call system('ls '//trim(settings%path_obs)//'*'//trim(settings%suffix)//&
' | wc -l > '//trim(settings%path_obsout)//'obsfiles.dat')
call system('ls '//trim(settings%path_obs)//' | grep '//trim(settings%suffix)//&
' >> '//trim(settings%path_obsout)//'obsfiles.dat')
else
call caldate(jdatei, jjjjmmdd, hhmiss)
write(adate,'(I6)') jjjjmmdd/100
call system('ls '//trim(settings%path_obs)//'*'//trim(settings%suffix)//&
' | wc -l > '//trim(settings%path_obsout)//adate//'/'//'obsfiles.dat')
call system('ls '//trim(settings%path_obs)//' | grep '//trim(settings%suffix)//&
' >> '//trim(settings%path_obsout)//adate//'/'//'obsfiles.dat')
endif
open(100,file=trim(settings%path_obsout)//'obsfiles.txt',action='read',status='old',iostat=ierr) if ( settings%lappendfile.eq.1 ) then
open(100,file=trim(settings%path_obsout)//'obsfiles.dat',action='read',status='old',iostat=ierr)
else
open(100,file=trim(settings%path_obsout)//adate//'/'//'obsfiles.dat',action='read',status='old',iostat=ierr)
endif
if ( ierr.ne. 0 ) then if ( ierr.ne. 0 ) then
write(*,*) 'ERROR: cannot open obsfiles.txt' write(*,*) 'ERROR: cannot open obsfiles.dat'
stop stop
endif endif
read(100,*,iostat=ierr) nfiles read(100,*,iostat=ierr) nfiles
......
prep_flexpart/main.f90
+ 36
− 20
View file @ 15f86f38
...@@ -41,6 +41,8 @@ program main ...@@ -41,6 +41,8 @@ program main
use mod_dates use mod_dates
use mod_settings use mod_settings
use mod_obs use mod_obs
use mod_prep_ageclass
use mod_prep_command
implicit none implicit none
...@@ -50,6 +52,7 @@ program main ...@@ -50,6 +52,7 @@ program main
character(max_path_len) :: settings_file character(max_path_len) :: settings_file
real(kind=8) :: jd real(kind=8) :: jd
character(len=6) :: adate
logical :: lexist logical :: lexist
integer :: nobs, nr integer :: nobs, nr
integer :: jjjjmmdd, hhmiss, eomday integer :: jjjjmmdd, hhmiss, eomday
...@@ -70,37 +73,47 @@ program main ...@@ -70,37 +73,47 @@ program main
jreldatei = juldate(settings%reldatei, 0) jreldatei = juldate(settings%reldatei, 0)
jreldatef = juldate(settings%reldatef, 0) jreldatef = juldate(settings%reldatef, 0)
print*, 'datei: ', settings%datei, jdatei write(*,*) 'datei: ', settings%datei, jdatei
print*, 'datef: ', settings%datef, jdatef write(*,*) 'datef: ', settings%datef, jdatef
print*, 'reldatei: ', settings%reldatei, jreldatei write(*,*) 'reldatei: ', settings%reldatei, jreldatei
print*, 'reldatef: ', settings%reldatef, jreldatef write(*,*) 'reldatef: ', settings%reldatef, jreldatef
! checks ! checks
if( jdatei.ne.jreldatei.or.jdatef.ne.jreldatef ) then if( jdatei.ne.jreldatei.or.jdatef.ne.jreldatef ) then
print*, 'WARNING: different run and release times' write(*,*) 'WARNING: different run and release times'
endif endif
if( (jreldatei.lt.jdatei).or.(jreldatei.gt.jdatef) ) then if( (jreldatei.lt.jdatei).or.(jreldatei.gt.jdatef) ) then
print*, 'ERROR: release date outside run time' write(*,*) 'ERROR: release date outside run time'
stop stop
endif endif
if( (jreldatef.lt.jdatei).or.(jreldatef.gt.jdatef) ) then if( (jreldatef.lt.jdatei).or.(jreldatef.gt.jdatef) ) then
print*, 'ERROR: release date outside run time' write(*,*) 'ERROR: release date outside run time'
stop stop
endif endif
inquire(file=settings%file_avail,exist=lexist) inquire(file=settings%file_avail,exist=lexist)
if( .not.lexist ) then if( .not.lexist ) then
print*, 'ERROR: AVAILABLE file not found' write(*,*) 'ERROR: AVAILABLE file not found'
stop stop
endif endif
if( settings%laverage.eq.0 ) then if( settings%laverage.eq.0 ) then
settings%intaverage = 0. settings%intaverage = 0.
print*, 'WARNING: laverage is 0 but intaverage > 0' write(*,*) 'WARNING: laverage is 0 but intaverage > 0'
print*, 'setting intaverage to 0' write(*,*) 'setting intaverage to 0'
endif endif
if( settings%lrelease.eq.0 ) then if( settings%lrelease.eq.0 ) then
print*, 'FLEXPART releases at regular intervals' write(*,*) 'FLEXPART releases at regular intervals'
else else
print*, 'FLEXPART releases at observation timestamps' write(*,*) 'FLEXPART releases at observation timestamps'
endif
! clean-up previous run
if ( settings%lappendfile.eq.0 ) then
write(*,*) 'WARNING: cleaning-up existing directory'
call caldate(jdatei, jjjjmmdd, hhmiss)
write(adate,'(I6)') jjjjmmdd/100
print*, trim(settings%path_obsout)//adate//'/'
call system('rm -rf '//trim(settings%path_obsout)//adate//'/')
call system('mkdir -p '//trim(settings%path_obsout)//adate//'/')
endif endif
! read receptor list ! read receptor list
...@@ -121,18 +134,20 @@ program main ...@@ -121,18 +134,20 @@ program main
call prep_pathnames(settings, jd, nr) call prep_pathnames(settings, jd, nr)
! write COMMAND ! write COMMAND
call prep_command(settings, jd, nr) if ( settings%FPversion.eq.10 ) then
call prep_command_FP10(settings, jd, nr)
else if ( settings%FPversion.eq.11 ) then
call prep_command(settings, jd, nr)
endif
! write OUTGRID ! write OUTGRID
call prep_outgrid(settings, jd, nr) call prep_outgrid(settings, jd, nr)
! write AGECLASS ! write AGECLASS
if ( settings%lFPversion.eq.0 ) then if ( settings%FPversion.eq.10 ) then
write(*,*) 'Preparing ageclass for FP10', settings%lFPversion call prep_ageclass_FP10(settings, jd, nr)
call prep_ageclass_FP10(settings, jd, nr) else if ( settings%FPversion.eq.11 ) then
else call prep_ageclass(settings, jd, nr)
write(*,*) 'Preparing ageclass for FP11', settings%lFPversion
call prep_ageclass(settings, jd, nr)
endif endif
! read and process observations ! read and process observations
...@@ -147,7 +162,8 @@ program main ...@@ -147,7 +162,8 @@ program main
else if ( settings%obsformat.eq.'basic' ) then else if ( settings%obsformat.eq.'basic' ) then
call read_basic(settings, jd, nr, nobs, obs) call read_basic(settings, jd, nr, nobs, obs)
else else
print*, 'ERROR: unknown observation file format' write(*,*) 'ERROR: unknown observation file format'
stop
endif endif
! write RELEASES ! write RELEASES
......
prep_flexpart/makefile
+ 4
− 4
View file @ 15f86f38
...@@ -18,15 +18,15 @@ MODULES = mod_var.o \ ...@@ -18,15 +18,15 @@ MODULES = mod_var.o \
mod_dates.o \ mod_dates.o \
mod_tools.o \ mod_tools.o \
mod_obs.o \ mod_obs.o \
mod_strings.o mod_strings.o \
mod_prep_ageclass.o \
mod_prep_command.o
OBJECTS = read_reclist.o \ OBJECTS = read_reclist.o \
list_obsfiles.o \ list_obsfiles.o \
prep_pathnames.o \ prep_pathnames.o \
prep_command.o \
prep_outgrid.o \ prep_outgrid.o \
prep_ageclass.o \ prep_reagents.o \
prep_ageclass_FP10.o \
prep_releases.o \ prep_releases.o \
prep_releases_reg.o \ prep_releases_reg.o \
process_obs.o \ process_obs.o \
......
prep_flexpart/mod_prep_ageclass.f90 0 → 100644
+ 147
− 0
View file @ 15f86f38
!---------------------------------------------------------------------------------------
! PREP_FLEXPART: mod_prep_ageclass
!---------------------------------------------------------------------------------------
! 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
!---------------------------------------------------------------------------------------
!
!> mod_prep_ageclass
!!
!! Purpose: Prepares the options file AGECLASS.
!!
!! Interface:
!!
!! Inputs
!! settings - settings data structure
!! jd - julian day of start of month
!! nr - index to receptor list
!!
!! Externals
!! caldate
!!
!---------------------------------------------------------------------------------------
module mod_prep_ageclass
use mod_var
use mod_settings
use mod_dates
implicit none
private
public :: prep_ageclass_FP10, prep_ageclass
character(len=max_path_len) :: filename
character(len=6) :: adate
integer :: jjjjmm, jjjjmmdd, hhmiss
integer :: datei, datef, eomday
integer :: ierr
integer :: nageclass, lage
contains
! --------------------------------------------------
! prep_ageclass
! --------------------------------------------------
!
! Prepares the AGECLASS file for FLEXPART v11
! --------------------------------------------------
subroutine prep_ageclass(settings, jd, nr)
implicit none
type (settings_t), intent(in) :: settings
real(kind=8), intent(in) :: jd
integer, intent(in) :: nr
namelist /nage/ &
nageclass
namelist /ageclass/ &
lage
! preset namelist variables
nageclass = 1
lage = settings%ageclass
call caldate(jd, jjjjmmdd, hhmiss)
write(adate,fmt='(I6.6)') jjjjmmdd/100
filename = trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/AGECLASSES'
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
write(100,nml=nage)
write(100,nml=ageclass)
close(100)
end subroutine prep_ageclass
! --------------------------------------------------
! prep_ageclass_FP10
! --------------------------------------------------
!
! Prepares the AGECLASS file for FLEXPART v10.4
! --------------------------------------------------
subroutine prep_ageclass_FP10(settings, jd, nr)
implicit none
type (settings_t), intent(in) :: settings
real(kind=8), intent(in) :: jd
integer, intent(in) :: nr
namelist /ageclass/ &
nageclass, &
lage
! preset namelist variables
nageclass = 1
lage = settings%ageclass
call caldate(jd, jjjjmmdd, hhmiss)
write(adate,fmt='(I6.6)') jjjjmmdd/100
filename = trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/AGECLASSES'
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
if( settings%lnamelist.eq.1 ) then
! use namelist file format
write(100,nml=ageclass)
else
! use old file format
write(100,fmt='(A48)') '************************************************'
write(100,fmt='(A48)') '* *'
write(100,fmt='(A48)') '*Lagrangian particle dispersion model FLEXPART *'
write(100,fmt='(A48)') '* Please select your options *'
write(100,fmt='(A48)') '* *'
write(100,fmt='(A48)') '*This file determines the ageclasses to be used*'
write(100,fmt='(A48)') '* *'
write(100,fmt='(A48)') '*Ages are given in seconds. The first class *'
write(100,fmt='(A48)') '*starts at age zero and goes up to the first *'
write(100,fmt='(A48)') '*age specified. The last age gives the maximum *'
write(100,fmt='(A48)') '*time a particle is carried in the simulation. *'
write(100,fmt='(A48)') '* *'
write(100,fmt='(A48)') '************************************************'
write(100,fmt='(A48)') '1 Integer Number of age classes '
write(100,fmt=*) settings%ageclass
endif
close(100)
end subroutine prep_ageclass_FP10
end module mod_prep_ageclass
prep_flexpart/prep_command.f90 → prep_flexpart/mod_prep_command.f90
+ 359
− 0
View file @ 15f86f38
!--------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------
! PREP_FLEXPART: prep_command ! PREP_FLEXPART: mod_prep_command
!--------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------
! FLEXINVERT is free software: you can redistribute it and/or modify ! FLEXINVERT is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by ! it under the terms of the GNU General Public License as published by
...@@ -17,7 +17,7 @@ ...@@ -17,7 +17,7 @@
! Copyright 2017, Rona Thompson ! Copyright 2017, Rona Thompson
!--------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------
! !
!> prep_command !> mod_prep_command
!! !!
!! Purpose: Prepares the options file COMMAND. !! Purpose: Prepares the options file COMMAND.
!! !!
...@@ -34,7 +34,7 @@ ...@@ -34,7 +34,7 @@
!! !!
!--------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------
subroutine prep_command(settings, jd, nr) module mod_prep_command
use mod_var use mod_var
use mod_settings use mod_settings
...@@ -42,10 +42,9 @@ subroutine prep_command(settings, jd, nr) ...@@ -42,10 +42,9 @@ subroutine prep_command(settings, jd, nr)
use mod_tools use mod_tools
implicit none implicit none
private
type (settings_t), intent(in) :: settings public :: prep_command, prep_command_FP10
real(kind=8), intent(in) :: jd
integer, intent(in) :: nr
character(len=max_path_len) :: filename character(len=max_path_len) :: filename
character(len=6) :: adate character(len=6) :: adate
...@@ -53,16 +52,14 @@ subroutine prep_command(settings, jd, nr) ...@@ -53,16 +52,14 @@ subroutine prep_command(settings, jd, nr)
integer :: jjjjmm, jjjjmmdd, hhmiss integer :: jjjjmm, jjjjmmdd, hhmiss
integer :: datei, datef, eomday integer :: datei, datef, eomday
integer :: ierr integer :: ierr
integer :: ldirect ! runtime mode (currently only backwards: ldirect = -1) integer :: ldirect ! runtime mode (currently only backwards: ldirect = -1)
integer :: ibdate, ibtime integer :: ibdate, ibtime
integer :: iedate, ietime integer :: iedate, ietime
! integer :: partsplit
integer :: loutstep integer :: loutstep
integer :: loutaver integer :: loutaver
integer :: loutsample integer :: loutsample
integer :: loutrestart integer :: loutrestart
! integer :: itsplit ! time constant for particle splitting (secs) integer :: itsplit ! time constant for particle splitting (secs)
integer :: lsynctime ! synchronisation interval of flexpart (secs) integer :: lsynctime ! synchronisation interval of flexpart (secs)
real :: ctl ! factor by which time step must be smaller than tl real :: ctl ! factor by which time step must be smaller than tl
integer :: ifine ! factor by which to decrease time step for vertical motion integer :: ifine ! factor by which to decrease time step for vertical motion
...@@ -76,13 +73,13 @@ subroutine prep_command(settings, jd, nr) ...@@ -76,13 +73,13 @@ subroutine prep_command(settings, jd, nr)
integer :: maxthreadgrid ! FP11 integer :: maxthreadgrid ! FP11
integer :: maxfilesize ! FP11 integer :: maxfilesize ! FP11
integer :: logvertinterp ! FP11 integer :: logvertinterp ! FP11
integer :: bcscheme ! FP11
integer :: lagespectra ! calculate age spectra (0 = no, 1 = yes) integer :: lagespectra ! calculate age spectra (0 = no, 1 = yes)
integer :: ipin ! continue simulation with dumped particle data (0 = no, 1 = yes) integer :: ipin ! continue simulation with dumped particle data (0 = no, 1 = yes)
integer :: ioutputforeachrelease ! create output file for each release location (0 = no, 1 = yes) integer :: ioutputforeachrelease ! create output file for each release location (0 = no, 1 = yes)
integer :: mdomainfill ! domain filling option (0 = no, 1 = yes) integer :: mdomainfill ! domain filling option (0 = no, 1 = yes)
integer :: mquasilag ! quasilagrangian mode to track particles (0 = no, 1 = yes) integer :: mquasilag ! quasilagrangian mode to track particles (0 = no, 1 = yes)
integer :: sfc_only ! save only surface layer in grid_time files = 1, full resolution = 0 integer :: sfc_only ! FP11 save only surface layer in grid_time files = 1, full resolution = 0
integer :: surf_only ! save only surface layer in grid_time files = 1, full resolution = 0
integer :: lnetcdfout ! netcdf output (0 = no, 1 = yes) integer :: lnetcdfout ! netcdf output (0 = no, 1 = yes)
integer :: cblflag ! integer :: cblflag !
integer :: linversionout ! one grid_time file per release = 1, or per timestep (original format) = 0 integer :: linversionout ! one grid_time file per release = 1, or per timestep (original format) = 0
...@@ -93,149 +90,270 @@ subroutine prep_command(settings, jd, nr) ...@@ -93,149 +90,270 @@ subroutine prep_command(settings, jd, nr)
integer :: linit_cond integer :: linit_cond
character(len=max_path_len) :: ohfields_path character(len=max_path_len) :: ohfields_path
namelist /command/ & contains
ldirect, &
ibdate,ibtime, & ! --------------------------------------------------
iedate,ietime, & ! prep_command
loutstep, & ! --------------------------------------------------
loutaver, & !
loutsample, & ! Prepares COMMAND for FLEXPARTv11
! itsplit, & ! --------------------------------------------------
lsynctime, &
ctl, & subroutine prep_command(settings, jd, nr)
ifine, &
iout, & implicit none
ipout, &
lsubgrid, & type (settings_t), intent(in) :: settings
lconvection, & real(kind=8), intent(in) :: jd
lturbulence, & integer, intent(in) :: nr
lturbulence_meso, &
nxshift, & namelist /command/ &
maxthreadgrid, & ldirect, &
maxfilesize, & ibdate,ibtime, &
logvertinterp, & iedate,ietime, &
bcscheme, & loutstep, &
loutrestart, & loutaver, &
lagespectra, & loutsample, &
ipin, & lsynctime, &
ioutputforeachrelease, & ctl, &
iflux, & ifine, &
mdomainfill, & iout, &
ind_source, & ipout, &
ind_receptor, & lsubgrid, &
mquasilag, & lconvection, &
nested_output, & lturbulence, &
linit_cond, & lturbulence_meso, &
lnetcdfout, & nxshift, &
sfc_only, & maxthreadgrid, &
cblflag, & maxfilesize, &
linversionout, & logvertinterp, &
ohfields_path loutrestart, &
lagespectra, &
! end of month day ipin, &
call caldate(jd, jjjjmmdd, hhmiss) ioutputforeachrelease, &
eomday = calceomday(jjjjmmdd/100) iflux, &
mdomainfill, &
! initialize command file path and name ind_source, &
write(adate,fmt='(I6.6)') jjjjmmdd/100 ind_receptor, &
filename = trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/COMMAND' mquasilag, &
nested_output, &
call system('mkdir -p '//trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/') linit_cond, &
lnetcdfout, &
! get dates sfc_only, &
! correct start date for length of back trajectory cblflag, &
call caldate((jd-settings%ageclass/3600d0/24d0), jjjjmmdd, hhmiss) linversionout, &
datei = jjjjmmdd ohfields_path
call caldate(jd, jjjjmmdd, hhmiss)
jdf = juldate(jjjjmmdd+eomday,0) ! end of month day
call caldate(min(jdf,jdatef),jjjjmmdd,hhmiss) call caldate(jd, jjjjmmdd, hhmiss)
datef = jjjjmmdd eomday = calceomday(jjjjmmdd/100)
! preset namelist variables ! initialize command file path and name
ldirect = -1 write(adate,fmt='(I6.6)') jjjjmmdd/100
ibdate = datei filename = trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/COMMAND'
ibtime = 0
iedate = datef call system('mkdir -p '//trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/')
ietime = 0
loutstep = settings%outrate ! get dates
loutaver = settings%outaverage ! correct start date for length of back trajectory
loutsample = settings%outsample call caldate((jd-settings%ageclass/3600d0/24d0), jjjjmmdd, hhmiss)
! itsplit = 999999999 datei = jjjjmmdd
lsynctime = lsync call caldate(jd, jjjjmmdd, hhmiss)
ctl = -5 jdf = juldate(jjjjmmdd+eomday,0)
ifine = 4 call caldate(min(jdf,jdatef),jjjjmmdd,hhmiss)
iout = 1 datef = jjjjmmdd
ipout = 0
lsubgrid = 1 ! preset namelist variables
lconvection = 1 ldirect = -1
lturbulence = 1 ibdate = datei
lturbulence_meso = 0 ibtime = 0
lnetcdfout = 0 iedate = datef
nxshift = 1 ietime = 0
maxthreadgrid = 1 loutstep = settings%outrate
maxfilesize = 10000 loutaver = settings%outaverage
logvertinterp = 0 loutsample = settings%outsample
bcscheme = 1 itsplit = 999999999
loutrestart = -1 lsynctime = lsync
lagespectra = 1 ctl = -5
ipin = 0 ifine = 4
ioutputforeachrelease = 1 iout = 1
iflux = 0 ipout = 0
mdomainfill = 0 lsubgrid = 1
ind_source = settings%ind_source lconvection = 1
ind_receptor = settings%ind_receptor lturbulence = 1
mquasilag = 0 lturbulence_meso = 0
nested_output = settings%lnested lnetcdfout = 0
linit_cond = settings%linit_cond nxshift = 1
sfc_only = 1 maxthreadgrid = 1
cblflag = 0 maxfilesize = 10000
linversionout = 1 logvertinterp = 0
ohfields_path = trim(settings%path_ohfield) loutrestart = -1
lagespectra = 1
! open command file ipin = 0
open(100,file=trim(filename),status='replace',action='write',iostat=ierr) ioutputforeachrelease = 1
iflux = 0
if( settings%lnamelist.eq.1 ) then mdomainfill = 0
! use namelist file format ind_source = settings%ind_source
write(100,nml=command) ind_receptor = settings%ind_receptor
else mquasilag = 0
! use old file format nested_output = settings%lnested
write(100,fmt='(A)') '********************************************************************************' linit_cond = settings%linit_cond
write(100,fmt='(A)') '* *' sfc_only = 1
write(100,fmt='(A)') '* Input file for the Lagrangian particle dispersion model FLEXPART *' cblflag = 0
write(100,fmt='(A)') '* Please select your options *' linversionout = 1
write(100,fmt='(A)') '* *' ohfields_path = trim(settings%path_ohfield)
write(100,fmt='(A)') '* *'
write(100,fmt='(A)') '********************************************************************************' open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
write(100,fmt='(I2)') ldirect write(100,nml=command)
write(100,fmt='(I8,1X,I6.6)') datei,0 close(100)
write(100,fmt='(I8,1X,I6.6)') datef,0
write(100,fmt='(I5)') settings%outrate end subroutine prep_command
write(100,fmt='(I5)') settings%outaverage
write(100,fmt='(I5)') settings%outsample ! --------------------------------------------------
! write(100,fmt='(I9)') itsplit ! prep_command_FP10
write(100,fmt='(I3,3X,A14)') lsynctime,'SYNC' ! --------------------------------------------------
write(100,fmt='(F5.2,1X,A10)') ctl,'CTL' !
write(100,fmt='(I3,3X,A14)') ifine,'IFINE' ! Prepares COMMAND for FLEXPARTv10
write(100,fmt='(I1,5X,A14)') iout,'IOUT' ! --------------------------------------------------
write(100,fmt='(I1,5X,A14)') ipout,'IPOUT'
write(100,fmt='(I1,5X,A14)') lsubgrid,'LSUBGRID' subroutine prep_command_FP10(settings, jd, nr)
write(100,fmt='(I1,5X,A14)') lconvection,'LCONVECTION'
write(100,fmt='(I1,5X,A14)') lagespectra,'LAGESPECTRA' implicit none
write(100,fmt='(I1,5X,A14)') ipin,'IPIN'
write(100,fmt='(I1,5X,A14)') ioutputforeachrelease,'IOUTPUTFOREACHRELEASE' type (settings_t), intent(in) :: settings
write(100,fmt='(I1,5X,A14)') iflux,'IFLUX' real(kind=8), intent(in) :: jd
write(100,fmt='(I1,5X,A14)') mdomainfill,'MDOMAINFILL' integer, intent(in) :: nr
write(100,fmt='(I1,5X,A14)') settings%ind_source,'IND_SOURCE'
write(100,fmt='(I1,5X,A14)') settings%ind_receptor,'IND_RECEPTOR' namelist /command/ &
write(100,fmt='(I1,5X,A14)') mquasilag,'MQUASILAG' ldirect, &
write(100,fmt='(I1,5X,A14)') settings%lnested,'NESTED_OUTPUT' ibdate,ibtime, &
write(100,fmt='(I1,5X,A14)') settings%linit_cond,'LINIT_COND' iedate,ietime, &
write(100,fmt='(I1,5X,A14)') sfc_only,'SURF_ONLY' loutstep, &
write(100,fmt='(I1,5X,A14)') linversionout,'LINVERSIONOUT' loutaver, &
endif loutsample, &
itsplit, &
close(100) lsynctime, &
ctl, &
end subroutine prep_command ifine, &
iout, &
ipout, &
lsubgrid, &
lconvection, &
lagespectra, &
ipin, &
ioutputforeachrelease, &
iflux, &
mdomainfill, &
ind_source, &
ind_receptor, &
mquasilag, &
nested_output, &
linit_cond, &
lnetcdfout, &
surf_only, &
cblflag, &
linversionout, &
ohfields_path
! end of month day
call caldate(jd, jjjjmmdd, hhmiss)
eomday = calceomday(jjjjmmdd/100)
! initialize command file path and name
write(adate,fmt='(I6.6)') jjjjmmdd/100
filename = trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/COMMAND'
call system('mkdir -p '//trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/')
! get dates
! correct start date for length of back trajectory
call caldate((jd-settings%ageclass/3600d0/24d0), jjjjmmdd, hhmiss)
datei = jjjjmmdd
call caldate(jd, jjjjmmdd, hhmiss)
jdf = juldate(jjjjmmdd+eomday,0)
call caldate(min(jdf,jdatef),jjjjmmdd,hhmiss)
datef = jjjjmmdd
! preset namelist variables
ldirect = -1
ibdate = datei
ibtime = 0
iedate = datef
ietime = 0
loutstep = settings%outrate
loutaver = settings%outaverage
loutsample = settings%outsample
itsplit = 999999999
lsynctime = lsync
ctl = -5
ifine = 4
iout = 1
ipout = 2
lsubgrid = 1
lconvection = 1
lagespectra = 1
ipin = 0
ioutputforeachrelease = 1
iflux = 0
mdomainfill = 0
ind_source = settings%ind_source
ind_receptor = settings%ind_receptor
mquasilag = 0
nested_output = settings%lnested
linit_cond = settings%linit_cond
lnetcdfout = 0
surf_only = 1
cblflag = 0
linversionout = 1
ohfields_path = trim(settings%path_ohfield)
! open command file
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
if( settings%lnamelist.eq.1 ) then
! use namelist file format
write(100,nml=command)
else
! use old file format
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)') '* Please select your options *'
write(100,fmt='(A)') '* *'
write(100,fmt='(A)') '* *'
write(100,fmt='(A)') '********************************************************************************'
write(100,fmt='(I2)') ldirect
write(100,fmt='(I8,1X,I6.6)') datei,0
write(100,fmt='(I8,1X,I6.6)') datef,0
write(100,fmt='(I5)') settings%outrate
write(100,fmt='(I5)') settings%outaverage
write(100,fmt='(I5)') settings%outsample
write(100,fmt='(I9)') itsplit
write(100,fmt='(I3,3X,A14)') lsynctime,'SYNC'
write(100,fmt='(F5.2,1X,A10)') ctl,'CTL'
write(100,fmt='(I3,3X,A14)') ifine,'IFINE'
write(100,fmt='(I1,5X,A14)') iout,'IOUT'
write(100,fmt='(I1,5X,A14)') ipout,'IPOUT'
write(100,fmt='(I1,5X,A14)') lsubgrid,'LSUBGRID'
write(100,fmt='(I1,5X,A14)') lconvection,'LCONVECTION'
write(100,fmt='(I1,5X,A14)') lagespectra,'LAGESPECTRA'
write(100,fmt='(I1,5X,A14)') ipin,'IPIN'
write(100,fmt='(I1,5X,A14)') ioutputforeachrelease,'IOUTPUTFOREACHRELEASE'
write(100,fmt='(I1,5X,A14)') iflux,'IFLUX'
write(100,fmt='(I1,5X,A14)') mdomainfill,'MDOMAINFILL'
write(100,fmt='(I1,5X,A14)') settings%ind_source,'IND_SOURCE'
write(100,fmt='(I1,5X,A14)') settings%ind_receptor,'IND_RECEPTOR'
write(100,fmt='(I1,5X,A14)') mquasilag,'MQUASILAG'
write(100,fmt='(I1,5X,A14)') settings%lnested,'NESTED_OUTPUT'
write(100,fmt='(I1,5X,A14)') settings%linit_cond,'LINIT_COND'
write(100,fmt='(I1,5X,A14)') sfc_only,'SURF_ONLY'
write(100,fmt='(I1,5X,A14)') linversionout,'LINVERSIONOUT'
endif
close(100)
end subroutine prep_command_FP10
end module mod_prep_command
prep_flexpart/mod_settings.f90
+ 76
− 5
View file @ 15f86f38
...@@ -48,10 +48,17 @@ module mod_settings ...@@ -48,10 +48,17 @@ module mod_settings
character(len=max_name_len) :: windfield character(len=max_name_len) :: windfield
character(len=max_name_len) :: obsformat character(len=max_name_len) :: obsformat
character(len=max_name_len) :: suffix character(len=max_name_len) :: suffix
integer :: lscratch
integer :: FPversion
integer :: nreagent
character(len=max_path_len), dimension(:), allocatable :: path_reagents
character(len=max_name_len), dimension(:), allocatable :: reagent_units
integer, dimension(:), allocatable :: reagent_interp
integer :: lFPversion
integer :: lselect integer :: lselect
integer :: lrelease integer :: lrelease
integer :: lappendfile
real :: relfreq real :: relfreq
integer :: laverage integer :: laverage
real :: intaverage real :: intaverage
...@@ -296,6 +303,9 @@ module mod_settings ...@@ -296,6 +303,9 @@ module mod_settings
stop stop
endif endif
! default settings
settings%lscratch = 1 ! use scratch drive
! loop over lines ! loop over lines
read_loop: do read_loop: do
read (100, fmt='(A)', iostat=ierr) line read (100, fmt='(A)', iostat=ierr) line
...@@ -313,6 +323,9 @@ module mod_settings ...@@ -313,6 +323,9 @@ module mod_settings
if ( len (trim (line)) .eq. 0 ) cycle read_loop if ( len (trim (line)) .eq. 0 ) cycle read_loop
! Path and file settings ! Path and file settings
identifier = "lscratch:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%lscratch = int(cn)
identifier = "path_flexpart:" identifier = "path_flexpart:"
call read_content (line, identifier, cc, cn, cl, match) call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%path_flexpart = cc if ( match ) settings%path_flexpart = cc
...@@ -349,17 +362,64 @@ module mod_settings ...@@ -349,17 +362,64 @@ module mod_settings
identifier = "windfield:" identifier = "windfield:"
call read_content (line, identifier, cc, cn, cl, match) call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%windfield = cc if ( match ) settings%windfield = cc
identifier = "FPversion:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%FPversion = int(cn)
! General settings ! Species settings
identifier = "lFPversion:" identifier = "nreagent:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%nreagent = int(cn)
identifier = "path_reagents:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) then
call parse_string(cc, ',', temp, n)
if ( n.ne.settings%nreagent ) then
write(*,*) 'ERROR reading SETTINGS: number of reagent paths not equal to nreagent'
stop
endif
allocate( settings%path_reagents(settings%nreagent) )
do i=1,settings%nreagent
settings%path_reagents(i) = trim(temp(i))
end do
endif
identifier = "reagent_units:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) then
call parse_string(cc, ',', temp, n)
if ( n.ne.settings%nreagent ) then
write(*,*) 'ERROR reading SETTINGS: number of reagent units not equal to nreagent'
stop
endif
allocate( settings%reagent_units(settings%nreagent) )
do i=1,settings%nreagent
settings%reagent_units(i) = trim(temp(i))
end do
endif
identifier = "reagent_interp:"
call read_content (line, identifier, cc, cn, cl, match) call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%lFPversion = int(cn) if ( match ) then
call parse_string(cc, ',', temp, n)
if ( n.ne.settings%nreagent ) then
write(*,*) 'ERROR reading SETTINGS: number of flags for interpolation not equal to nreagent'
stop
endif
allocate( settings%reagent_interp(settings%nreagent) )
do i=1,settings%nreagent
read(temp(i),*) settings%reagent_interp(i)
end do
endif
! General settings
identifier = "lselect:" identifier = "lselect:"
call read_content (line, identifier, cc, cn, cl, match) call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%lselect = int(cn) if ( match ) settings%lselect = int(cn)
identifier = "lrelease:" identifier = "lrelease:"
call read_content (line, identifier, cc, cn, cl, match) call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%lrelease = int(cn) if ( match ) settings%lrelease = int(cn)
identifier = "lappendfile:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%lappendfile = int(cn)
identifier = "relfreq:" identifier = "relfreq:"
call read_content (line, identifier, cc, cn, cl, match) call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%relfreq = real(cn) if ( match ) settings%relfreq = real(cn)
...@@ -374,7 +434,7 @@ module mod_settings ...@@ -374,7 +434,7 @@ module mod_settings
if ( match ) settings%lnamelist = int(cn) if ( match ) settings%lnamelist = int(cn)
identifier = "measerr:" identifier = "measerr:"
call read_content (line, identifier, cc, cn, cl, match) call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%measerr = int(cn) if ( match ) settings%measerr = real(cn)
! COMMAND settings ! COMMAND settings
identifier = "datei:" identifier = "datei:"
...@@ -485,6 +545,17 @@ module mod_settings ...@@ -485,6 +545,17 @@ module mod_settings
print*, 'path_flexpart: ',trim(settings%path_flexpart) print*, 'path_flexpart: ',trim(settings%path_flexpart)
print*, 'path_options: ',trim(settings%path_options) print*, 'path_options: ',trim(settings%path_options)
! Some checks
if ( settings%FPversion.lt.10.or.settings%FPversion.gt.11 ) then
write(*,*) 'ERROR: FPversion must be either 10 or 11'
stop
endif
if ( settings%lnamelist.ne.1.and.settings%FPversion.eq.11 ) then
write(*,*) 'WARNING: only namelist file format available for FPversion 11'
settings%lnamelist = 1
endif
nreagent = settings%nreagent
end subroutine read_settings end subroutine read_settings
......
prep_flexpart/mod_var.f90
+ 1
− 0
View file @ 15f86f38
...@@ -36,6 +36,7 @@ module mod_var ...@@ -36,6 +36,7 @@ module mod_var
real(kind=8) :: jreldatei, jreldatef real(kind=8) :: jreldatei, jreldatef
integer :: nrec integer :: nrec
integer :: nfiles integer :: nfiles
integer :: nreagent
character(len=recname_len), dimension(:), allocatable :: recname character(len=recname_len), dimension(:), allocatable :: recname
character(len=256), dimension(:), allocatable :: filelist character(len=256), dimension(:), allocatable :: filelist
real, dimension(:), allocatable :: reclon, reclat, recalt real, dimension(:), allocatable :: reclon, reclat, recalt
......
prep_flexpart/prep_outgrid.f90
+ 1
− 1
View file @ 15f86f38
...@@ -84,7 +84,7 @@ subroutine prep_outgrid(settings, jd, nr) ...@@ -84,7 +84,7 @@ subroutine prep_outgrid(settings, jd, nr)
if ( settings%windfield.eq.'oper' ) then if ( settings%windfield.eq.'oper' ) then
outlon0 = -179. outlon0 = -179.
else if ( settings%windfield.eq.'era5' ) then else if ( settings%windfield.eq.'era5' ) then
outlon0 = -179.5 outlon0 = -179.
else else
write(*,*) 'ERROR: prep_outgrid: unknown windfield type' write(*,*) 'ERROR: prep_outgrid: unknown windfield type'
stop stop
......
prep_flexpart/prep_pathnames.f90
+ 10
− 5
View file @ 15f86f38
...@@ -59,16 +59,21 @@ subroutine prep_pathnames(settings, jd, nr) ...@@ -59,16 +59,21 @@ subroutine prep_pathnames(settings, jd, nr)
open(100,file=trim(filename),status='replace',action='write',iostat=ierr) open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
write(100,fmt='(A)') trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/' if ( settings%lscratch.eq.0 ) then
write(100,fmt='(A)') trim(settings%path_output)//trim(recname(nr))//'/'//adate//'/' ! not using scratch
write(100,fmt='(A)') trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/'
write(100,fmt='(A)') trim(settings%path_output)//trim(recname(nr))//'/'//adate//'/'
else
! using scratch
write(100,fmt='(A)') './options/'
write(100,fmt='(A)') './output/'
endif
write(100,fmt='(A1)') '/' write(100,fmt='(A1)') '/'
write(100,fmt='(A)') trim(settings%file_avail) write(100,fmt='(A)') trim(settings%file_avail)
if( settings%lnested.eq.1 ) then if( settings%file_availnest.ne.'' ) then
write(100,fmt='(A1)') '/' write(100,fmt='(A1)') '/'
write(100,fmt='(A)') trim(settings%file_availnest) write(100,fmt='(A)') trim(settings%file_availnest)
write(100,fmt='(A)') '============================================'
write(100,fmt='(A1)') ' '
endif endif
close(100) close(100)
......
prep_flexpart/prep_ageclass.f90 → prep_flexpart/prep_reagents.f90
+ 73
− 0
View file @ 15f86f38
!--------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------
! PREP_FLEXPART: prep_ageclass ! PREP_FLEXPART: prep_reagents
!--------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------
! FLEXINVERT is free software: you can redistribute it and/or modify ! FLEXINVERT is free software: you can redistribute it and/or modify
! it under the terms of the GNU General Public License as published by ! it under the terms of the GNU General Public License as published by
...@@ -17,83 +17,57 @@ ...@@ -17,83 +17,57 @@
! Copyright 2017, Rona Thompson ! Copyright 2017, Rona Thompson
!--------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------
! !
!> prep_ageclass !> prep_reagents
!! !!
!! Purpose: Prepares the options file AGECLASS. !! Purpose: Prepares the options file REAGENTS
!! !!
!! Interface: !! Interface:
!! !!
!! Inputs !! Inputs
!! settings - settings data structure !! settings - settings data structure
!! jd - julian day of start of month !! reagent - character string of reagent name (see SPECIES files)
!! nr - index to receptor list
!!
!! Externals
!! caldate
!! !!
!--------------------------------------------------------------------------------------- !---------------------------------------------------------------------------------------
subroutine prep_ageclass(settings, jd, nr) subroutine prep_reagents(settings,reagents,nr,jd)
use mod_var use mod_var
use mod_settings
use mod_dates use mod_dates
use mod_settings
implicit none implicit none
type (settings_t), intent(in) :: settings type (settings_t), intent(in) :: settings
real(kind=8), intent(in) :: jd character(len=6), dimension(nreagent), intent(in) :: reagents
integer, intent(in) :: nr integer, intent(in) :: nr
real(kind=8), intent(in) :: jd
character(len=max_path_len) :: filename
character(len=6) :: adate
integer :: jjjjmm, jjjjmmdd, hhmiss
integer :: datei, datef, eomday
integer :: ierr
integer :: nageclass, lage
! FP 11 character(len=max_path_len) :: filename, preag_path
namelist /nage/ & character(len=max_name_len) :: preagent, punit
nageclass character(len=6) :: adate
namelist /ageclass/ & integer :: phourly
lage integer :: jjjjmmdd, hhmiss, ierr, n
namelist /reagent_params/ &
preagent, &
! preset namelist variables preag_path, &
nageclass = 1 phourly, &
lage = settings%ageclass punit
call caldate(jd, jjjjmmdd, hhmiss) call caldate(jd, jjjjmmdd, hhmiss)
write(adate,fmt='(I6.6)') jjjjmmdd/100 write(adate,fmt='(I6.6)') jjjjmmdd/100
filename = trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/AGECLASSES' filename = trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/REAGENTS'
open(100,file=trim(filename),status='replace',action='write',iostat=ierr) open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
do n = 1, nreagent
if( settings%lnamelist.eq.1 ) then if (reagents(n).ne."") then
! use namelist file format, FP11 preagent = trim(reagents(n))
write(100,nml=nage) preag_path = trim(settings%path_reagents(n))
write(100,nml=ageclass) phourly = settings%reagent_interp(n)
else punit = trim(settings%reagent_units(n))
! use old file format write(100,nml=reagent_params)
write(100,fmt='(A48)') '************************************************' endif
write(100,fmt='(A48)') '* *' end do
write(100,fmt='(A48)') '*Lagrangian particle dispersion model FLEXPART *'
write(100,fmt='(A48)') '* Please select your options *'
write(100,fmt='(A48)') '* *'
write(100,fmt='(A48)') '*This file determines the ageclasses to be used*'
write(100,fmt='(A48)') '* *'
write(100,fmt='(A48)') '*Ages are given in seconds. The first class *'
write(100,fmt='(A48)') '*starts at age zero and goes up to the first *'
write(100,fmt='(A48)') '*age specified. The last age gives the maximum *'
write(100,fmt='(A48)') '*time a particle is carried in the simulation. *'
write(100,fmt='(A48)') '* *'
write(100,fmt='(A48)') '************************************************'
write(100,fmt='(A48)') '1 Integer Number of age classes '
write(100,fmt=*) settings%ageclass
endif
close(100) close(100)
end subroutine prep_reagents
end subroutine prep_ageclass
prep_flexpart/prep_releases.f90
+ 74
− 11
View file @ 15f86f38
...@@ -43,6 +43,7 @@ subroutine prep_releases(settings, jd, nr, nobs, obs) ...@@ -43,6 +43,7 @@ subroutine prep_releases(settings, jd, nr, nobs, obs)
use mod_dates use mod_dates
use mod_obs use mod_obs
use mod_tools use mod_tools
use mod_strings
implicit none implicit none
...@@ -70,6 +71,23 @@ subroutine prep_releases(settings, jd, nr, nobs, obs) ...@@ -70,6 +71,23 @@ subroutine prep_releases(settings, jd, nr, nobs, obs)
integer :: zkind, parts integer :: zkind, parts
real :: mass real :: mass
character(len=40) :: comment character(len=40) :: comment
logical :: lexist
character(len=16) :: pspecies, pemis_name
character(len=256) :: pemis_path, pemis_file
integer :: pemis_unit
real :: pemis_coeff
real :: pdecay, pweta_gas, pwetb_gas, preldiff, phenry, pf0, pdensity, pdquer
real :: pdsigma, pdryvel, pweightmolar, pdia
character(len=10), allocatable, dimension(:) :: preactions
real, allocatable, dimension(:) :: pcconst, pdconst, pnconst
real, allocatable, dimension(:) :: pohnconst, pohcconst, pohdconst
real :: pcrain_aero, pcsnow_aero, pccn_aero, pin_aero
real :: parea_dow(7), parea_hour(24), ppoint_dow(7), ppoint_hour(24)
integer :: ios
integer :: pshape,porient
real ::pla,pia,psa,f,e,paspectratio
namelist /releases_ctrl/ & namelist /releases_ctrl/ &
nspec, & nspec, &
...@@ -86,20 +104,65 @@ subroutine prep_releases(settings, jd, nr, nobs, obs) ...@@ -86,20 +104,65 @@ subroutine prep_releases(settings, jd, nr, nobs, obs)
parts, & parts, &
comment comment
namelist /species_params/ &
pspecies, pdecay, pweta_gas, pwetb_gas, &
pcrain_aero, pcsnow_aero, pccn_aero, pin_aero, &
preldiff, phenry, pf0, pdensity, pdquer, pdia, &
pdsigma, pdryvel, pweightmolar, pohnconst, &
preactions, pcconst, pdconst, pnconst, pohcconst, pohdconst, &
pemis_path, pemis_file, pemis_name, pemis_unit, pemis_coeff, &
parea_dow, parea_hour, ppoint_dow, ppoint_hour, &
pshape, paspectratio, pla, pia, psa, porient
call caldate(jd, jjjjmmdd, hhmiss) call caldate(jd, jjjjmmdd, hhmiss)
write(adate,fmt='(I6.6)') jjjjmmdd/100 write(adate,fmt='(I6.6)') jjjjmmdd/100
! read species file ! read species file
open(100,file=trim(settings%path_flexpart)//'options/SPECIES/spec_overview',status='old',action='read',iostat=ierr) if ( settings%FPversion.eq.10 ) then
nchar = len_trim(settings%species) inquire(file=trim(settings%path_flexpart)//'options/SPECIES/spec_overview',exist=lexist)
do while ( ierr.eq.0 ) if ( .not.lexist ) then
read (100, fmt='(A)', iostat=ierr) line write(*,*) 'ERROR: cannot find file /options/SPECIES/spec_overview'
! sec - reads the last char if ( line(len_trim(line)-nchar+1:len_trim(line)) == trim(settings%species) ) ierr = 1 stop
if ( line(5:5+nchar-1) == trim(settings%species) ) ierr = 1 endif
end do open(100,file=trim(settings%path_flexpart)//'options/SPECIES/spec_overview',status='old',action='read',iostat=ierr)
! sec - read(line(9:11),*) spec nchar = len_trim(settings%species)
read(line(1:3),*) spec do while ( ierr.eq.0 )
write(aspec,fmt='(I3.3)') spec 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
else if ( settings%FPversion.eq.11 ) then
! FPv11 new version SPECIES files contain species name not number
spec = 1
inquire(file=trim(settings%path_flexpart)//'options/SPECIES/SPECIES_'//to_upper(trim(settings%species)),&
exist=lexist)
if ( .not.lexist ) then
write(*,*) 'ERROR: cannot find file /options/SPECIES/SPECIES_'//to_upper(trim(settings%species))
stop
endif
write(aspec,fmt='(I3.3)') spec
call system('cp '//trim(settings%path_flexpart)//'options/SPECIES/SPECIES_'//to_upper(trim(settings%species))//&
' '//trim(settings%path_flexpart)//'options/SPECIES/SPECIES_'//aspec)
write(aspec,fmt='(I3.3)') spec
! read SPECIES file for reagents info
allocate(preactions(nreagent))
allocate(pcconst(nreagent))
allocate(pdconst(nreagent))
allocate(pnconst(nreagent))
open(100,file=trim(settings%path_flexpart)//'options/SPECIES/SPECIES_'//to_upper(trim(settings%species)),&
action='read',iostat=ierr)
read(100,species_params,iostat=ierr)
if ( ierr.ne.0 ) then
write(*,*) 'ERROR: reading SPECIES_'//to_upper(trim(settings%species))
stop
endif
write(*,*) 'Chemical reactions: ',preactions
if (any(preactions.ne."")) then
! write reagents file
call prep_reagents(settings,preactions,nr,jd)
endif
endif
! preset namelist variables releases_ctrl ! preset namelist variables releases_ctrl
nspec = 1 nspec = 1
...@@ -108,7 +171,7 @@ subroutine prep_releases(settings, jd, nr, nobs, obs) ...@@ -108,7 +171,7 @@ subroutine prep_releases(settings, jd, nr, nobs, obs)
! copy standard input files to options folder ! copy standard input files to options folder
pathname = trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/' pathname = trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/'
call system('mkdir -p '//trim(pathname)//'SPECIES/') call system('mkdir -p '//trim(pathname)//'SPECIES/')
filename = trim(settings%path_flexpart)//'options/SPECIES/SPECIES_'//aspec filename = trim(settings%path_flexpart)//'options/SPECIES/SPECIES_'//trim(aspec)
call system('cp '//trim(filename)//' '//trim(pathname)//'SPECIES/') 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/*.dat'//' '//trim(pathname))
call system('cp '//trim(settings%path_flexpart)//'options/*.t'//' '//trim(pathname)) call system('cp '//trim(settings%path_flexpart)//'options/*.t'//' '//trim(pathname))
......
prep_flexpart/prep_releases_reg.f90
+ 72
− 8
View file @ 15f86f38
...@@ -45,6 +45,7 @@ subroutine prep_releases_reg(settings, jd, nr) ...@@ -45,6 +45,7 @@ subroutine prep_releases_reg(settings, jd, nr)
use mod_settings use mod_settings
use mod_dates use mod_dates
use mod_tools use mod_tools
use mod_strings
implicit none implicit none
...@@ -72,6 +73,22 @@ subroutine prep_releases_reg(settings, jd, nr) ...@@ -72,6 +73,22 @@ subroutine prep_releases_reg(settings, jd, nr)
integer :: zkind, parts integer :: zkind, parts
real :: mass real :: mass
character(len=40) :: comment character(len=40) :: comment
logical :: lexist
character(len=16) :: pspecies, pemis_name
character(len=256) :: pemis_path, pemis_file
integer :: pemis_unit
real :: pemis_coeff
real :: pdecay, pweta_gas, pwetb_gas, preldiff, phenry, pf0, pdensity, pdquer
real :: pdsigma, pdryvel, pweightmolar, pdia
character(len=10), allocatable, dimension(:) :: preactions
real, allocatable, dimension(:) :: pcconst, pdconst, pnconst
real, allocatable, dimension(:) :: pohnconst, pohcconst, pohdconst
real :: pcrain_aero, pcsnow_aero, pccn_aero, pin_aero
real :: parea_dow(7), parea_hour(24), ppoint_dow(7), ppoint_hour(24)
integer :: ios
integer :: pshape,porient
real ::pla,pia,psa,f,e,paspectratio
namelist /releases_ctrl/ & namelist /releases_ctrl/ &
nspec, & nspec, &
...@@ -88,6 +105,16 @@ subroutine prep_releases_reg(settings, jd, nr) ...@@ -88,6 +105,16 @@ subroutine prep_releases_reg(settings, jd, nr)
parts, & parts, &
comment comment
namelist /species_params/ &
pspecies, pdecay, pweta_gas, pwetb_gas, &
pcrain_aero, pcsnow_aero, pccn_aero, pin_aero, &
preldiff, phenry, pf0, pdensity, pdquer, pdia, &
pdsigma, pdryvel, pweightmolar, pohnconst, &
preactions, pcconst, pdconst, pnconst, pohcconst, pohdconst, &
pemis_path, pemis_file, pemis_name, pemis_unit, pemis_coeff, &
parea_dow, parea_hour, ppoint_dow, ppoint_hour, &
pshape, paspectratio, pla, pia, psa, porient
! initialize date/time variables ! initialize date/time variables
call caldate(jd, jjjjmmdd, hhmiss) call caldate(jd, jjjjmmdd, hhmiss)
write(adate,fmt='(I6.6)') jjjjmmdd/100 write(adate,fmt='(I6.6)') jjjjmmdd/100
...@@ -100,14 +127,51 @@ subroutine prep_releases_reg(settings, jd, nr) ...@@ -100,14 +127,51 @@ subroutine prep_releases_reg(settings, jd, nr)
print*, 'prep_releases_reg: jdf = ',jdf print*, 'prep_releases_reg: jdf = ',jdf
! read species file ! read species file
open(100,file=trim(settings%path_flexpart)//'options/SPECIES/spec_overview',status='old',action='read',iostat=ierr) if ( settings%FPversion.eq.10 ) then
nchar = len_trim(settings%species) inquire(file=trim(settings%path_flexpart)//'options/SPECIES/spec_overview',exist=lexist)
do while ( ierr.eq.0 ) if ( .not.lexist ) then
read (100, fmt='(A)', iostat=ierr) line write(*,*) 'ERROR: cannot fine file /options/SPECIES/spec_overview'
if ( line(len_trim(line)-nchar+1:len_trim(line)) == trim(settings%species) ) ierr = 1 stop
end do endif
read(line(9:11),*) spec open(100,file=trim(settings%path_flexpart)//'options/SPECIES/spec_overview',status='old',action='read',iostat=ierr)
write(aspec,fmt='(I3.3)') spec 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
else if ( settings%FPversion.eq.11 ) then
! FPv11 new version SPECIES files contain species name not number
spec = 1
inquire(file=trim(settings%path_flexpart)//'options/SPECIES/SPECIES_'//to_upper(trim(settings%species)),&
exist=lexist)
if ( .not.lexist ) then
write(*,*) 'ERROR: cannot find file /options/SPECIES/SPECIES_'//to_upper(trim(settings%species))
stop
endif
write(aspec,fmt='(I3.3)') spec
call system('cp '//trim(settings%path_flexpart)//'options/SPECIES/SPECIES_'//to_upper(trim(settings%species))//&
' '//trim(settings%path_flexpart)//'options/SPECIES/SPECIES_'//aspec)
write(aspec,fmt='(I3.3)') spec
! read SPECIES file for reagents info
allocate(preactions(nreagent))
allocate(pcconst(nreagent))
allocate(pdconst(nreagent))
allocate(pnconst(nreagent))
open(100,file=trim(settings%path_flexpart)//'options/SPECIES/SPECIES_'//to_upper(trim(settings%species)),&
action='read',iostat=ierr)
read(100,species_params,iostat=ierr)
if ( ierr.ne.0 ) then
write(*,*) 'ERROR: reading SPECIES_'//to_upper(trim(settings%species))
stop
endif
write(*,*) 'Chemical reactions: ',preactions
if (any(preactions.ne."")) then
! write reagents file
call prep_reagents(settings,preactions,nr,jd)
endif
endif
! preset namelist variables releases_ctrl ! preset namelist variables releases_ctrl
nspec = 1 nspec = 1
......
prep_flexpart/read_basic.f90
+ 20
− 4
View file @ 15f86f38
...@@ -63,6 +63,7 @@ subroutine read_basic(settings, jd, nr, nobs, obs) ...@@ -63,6 +63,7 @@ subroutine read_basic(settings, jd, nr, nobs, obs)
character(len=max_path_len) :: rowfmt character(len=max_path_len) :: rowfmt
character(len=200) :: line character(len=200) :: line
character(len=200), dimension(50) :: args character(len=200), dimension(50) :: args
character(len=6) :: adate
character(len=3) :: flag character(len=3) :: flag
real, dimension(12) :: temp real, dimension(12) :: temp
logical :: lexist logical :: lexist
...@@ -166,7 +167,7 @@ subroutine read_basic(settings, jd, nr, nobs, obs) ...@@ -166,7 +167,7 @@ subroutine read_basic(settings, jd, nr, nobs, obs)
! no observations go to end ! no observations go to end
if ( nobs.eq.0 ) then if ( nobs.eq.0 ) then
call caldate(jdate, jjjjmmdd, hhmiss) call caldate(jd, jjjjmmdd, hhmiss)
print*, 'WARNING: no data for ',trim(recname(nr)),' in year/month ',jjjjmmdd/100 print*, 'WARNING: no data for ',trim(recname(nr)),' in year/month ',jjjjmmdd/100
go to 20 go to 20
endif endif
...@@ -181,7 +182,17 @@ subroutine read_basic(settings, jd, nr, nobs, obs) ...@@ -181,7 +182,17 @@ subroutine read_basic(settings, jd, nr, nobs, obs)
if ( nobs.eq.0 ) go to 20 if ( nobs.eq.0 ) go to 20
! write formatted obs to file ! write formatted obs to file
file_out = trim(settings%path_obsout)//trim(recname(nr))//'_'//trim(settings%species)//'.gob' if ( settings%lappendfile.eq.1 ) then
! append file
! generic output directory
file_out = trim(settings%path_obsout)//trim(recname(nr))//'_'//trim(settings%species)//'.txt'
else
! do not append to file when this is a new run
! create subdirectory for ouput based on start year, month
call caldate(jdatei, jjjjmmdd, hhmiss)
write(adate,'(I6)') jjjjmmdd/100
file_out = trim(settings%path_obsout)//adate//'/'//trim(recname(nr))//'_'//trim(settings%species)//'.txt'
endif
inquire(file=trim(file_out),exist=lexist) inquire(file=trim(file_out),exist=lexist)
if( lexist ) then if( lexist ) then
! append to existing ! append to existing
...@@ -212,8 +223,13 @@ subroutine read_basic(settings, jd, nr, nobs, obs) ...@@ -212,8 +223,13 @@ subroutine read_basic(settings, jd, nr, nobs, obs)
close(100) close(100)
! write detailed station list file ! write detailed station list file
print*, 'lat, lon , alt: ',obs%lat(1), obs%lon(1), obs%alt(1) if ( settings%lappendfile.eq.1 ) then
file_out = trim(settings%path_obsout)//'stnlist_detail.txt' ! append file
file_out = trim(settings%path_obsout)//'stnlist_detail.dat'
else
! do not append
file_out = trim(settings%path_obsout)//adate//'/'//'stnlist_detail.dat'
endif
inquire(file=trim(file_out),exist=lexist) inquire(file=trim(file_out),exist=lexist)
write(rowfmt,fmt='(A,I1,A)') '(A',recname_len,',1X,F7.2,1X,F7.2,1X,F7.1)' write(rowfmt,fmt='(A,I1,A)') '(A',recname_len,',1X,F7.2,1X,F7.2,1X,F7.1)'
if( lexist ) then if( lexist ) then
......
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