From 15f86f38bff70081d6c79a6d4f83308898fbc4b4 Mon Sep 17 00:00:00 2001
From: ronesy <rlt@nilu.no>
Date: Thu, 24 Oct 2024 15:46:36 +0200
Subject: [PATCH] changes for FPv11 with backwards compatibility with FPv10.4
and name change for aerosol species
---
grid_to_ncdf/grid_convert.f90 | 12 +-
grid_to_ncdf/init_convert.f90 | 6 +-
grid_to_ncdf/initialize.f90 | 19 +-
grid_to_ncdf/makefile | 61 ++--
grid_to_ncdf/mod_flexpart.f90 | 16 +-
grid_to_ncdf/run_grid2ncdf.sh | 16 +
prep_flexpart/SETTINGS | 22 +-
prep_flexpart/list_obsfiles.f90 | 28 +-
prep_flexpart/main.f90 | 56 +--
prep_flexpart/makefile | 8 +-
prep_flexpart/mod_prep_ageclass.f90 | 147 ++++++++
prep_flexpart/mod_prep_command.f90 | 359 ++++++++++++++++++++
prep_flexpart/mod_settings.f90 | 81 ++++-
prep_flexpart/mod_var.f90 | 1 +
prep_flexpart/prep_ageclass.f90 | 99 ------
prep_flexpart/prep_command.f90 | 241 -------------
prep_flexpart/prep_outgrid.f90 | 2 +-
prep_flexpart/prep_pathnames.f90 | 15 +-
prep_flexpart/prep_reagents.f90 | 73 ++++
prep_flexpart/prep_releases.f90 | 85 ++++-
prep_flexpart/prep_releases_reg.f90 | 80 ++++-
prep_flexpart/read_basic.f90 | 24 +-
prep_flexpart/read_icos.f90 | 26 +-
prep_flexpart/read_noaa.f90 | 24 +-
prep_flexpart/read_obspack.f90 | 23 +-
prep_flexpart/read_wdcgg.f90 | 25 +-
prep_flexpart/sbatch_prep_flexpart.sh | 2 +-
prep_fluxes/makefile | 62 ++--
prep_regions/get_surfinf.f90 | 6 +-
prep_regions/initialize.f90 | 35 +-
prep_regions/mod_flexpart.f90 | 16 +-
prep_regions/mod_settings.f90 | 2 +-
prep_regions/prep_regions.f90 | 9 -
prep_regions/run_prepregions.sh | 27 ++
prep_satellite/makefile | 70 ++--
prep_syndata/get_error_cov.f90 | 10 +-
prep_syndata/get_flux.f90 | 4 +-
prep_syndata/get_obs.f90 | 6 +-
prep_syndata/initialize.f90 | 40 ++-
prep_syndata/makefile | 72 ++--
prep_syndata/mod_flexpart.f90 | 16 +-
prep_syndata/mod_perturb.f90 | 6 +-
prep_syndata/mod_save.f90 | 48 +--
prep_syndata/mod_var.f90 | 2 +-
prep_syndata/prep_syndata.f90 | 3 -
prep_syndata/run_prepsyndata.sh | 24 ++
settings/SETTINGS_aero_config | 151 ++++++++
settings/SETTINGS_aero_files | 89 +++++
settings/SETTINGS_co2_config | 3 +
settings/SETTINGS_ghg_config | 3 +
source/average_fp.f90 | 6 +-
source/calc_conc.f90 | 15 +-
source/init_aero.f90 | 472 ++++++++++++++++++++++++++
source/init_cini.f90 | 6 +-
source/init_cini_month.f90 | 6 +-
source/init_ghg.f90 | 2 +-
source/initialize.f90 | 19 +-
source/job_flexinvert.sh | 4 +-
source/main.f90 | 6 +-
source/makefile | 1 +
source/mod_flexpart.f90 | 20 +-
source/mod_perturb.f90 | 3 +-
source/mod_save.f90 | 35 +-
source/mod_settings.f90 | 15 +-
source/read_obs.f90 | 76 ++---
source/sbatch_flexinvert.sh | 4 +-
source/simulate.f90 | 23 +-
67 files changed, 2221 insertions(+), 747 deletions(-)
create mode 100755 grid_to_ncdf/run_grid2ncdf.sh
create mode 100644 prep_flexpart/mod_prep_ageclass.f90
create mode 100644 prep_flexpart/mod_prep_command.f90
delete mode 100644 prep_flexpart/prep_ageclass.f90
delete mode 100644 prep_flexpart/prep_command.f90
create mode 100644 prep_flexpart/prep_reagents.f90
create mode 100755 prep_regions/run_prepregions.sh
create mode 100755 prep_syndata/run_prepsyndata.sh
create mode 100644 settings/SETTINGS_aero_config
create mode 100644 settings/SETTINGS_aero_files
create mode 100644 source/init_aero.f90
diff --git a/grid_to_ncdf/grid_convert.f90 b/grid_to_ncdf/grid_convert.f90
index 4c3adfd..6fb61a5 100644
--- a/grid_to_ncdf/grid_convert.f90
+++ b/grid_to_ncdf/grid_convert.f90
@@ -51,7 +51,7 @@ subroutine grid_convert(files, config, nr, timestamp)
integer :: ibdate, ibtime
integer, dimension(maxpoint) :: releases
real :: bndx, bndy, delx, dely
- integer :: numx, numy
+ integer :: numx, numy, numz
! list all grid files in directory
if ( trim(config%filefreq).eq.'month' ) then
@@ -90,17 +90,23 @@ subroutine grid_convert(files, config, nr, timestamp)
! read header for this month to get ngrid
if ( config%lnested ) then
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
filename = trim(path_flexpart)//'header'
+ inquire(file=trim(filename),exist=lexist)
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)
+ numx, numy, numz, bndx, bndy, delx, dely)
if ( config%linversionout ) then
ngrid = trajdays*ndgrid+2
else
diff --git a/grid_to_ncdf/init_convert.f90 b/grid_to_ncdf/init_convert.f90
index a57bfeb..1fb8f0f 100644
--- a/grid_to_ncdf/init_convert.f90
+++ b/grid_to_ncdf/init_convert.f90
@@ -52,7 +52,7 @@ subroutine init_convert(files, config, nr, timestamp)
integer :: ibdate, ibtime
integer, dimension(maxpoint) :: releases
real :: bndx, bndy, delx, dely
- integer :: numx, numy
+ integer :: numx, numy, numz
! list all grid files in directory
if ( trim(config%filefreq).eq.'month' ) then
@@ -92,8 +92,8 @@ subroutine init_convert(files, config, nr, timestamp)
endif
print*, 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely)
- print*, 'numx, numy, bndx, bndy, delx, dely = ',numx,numy,bndx,bndy,delx,dely
+ numx, numy, numz, 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( levels(nzgrid) )
! trick to use write_ncdf also for grid_initial
diff --git a/grid_to_ncdf/initialize.f90 b/grid_to_ncdf/initialize.f90
index ab0df23..451750a 100644
--- a/grid_to_ncdf/initialize.f90
+++ b/grid_to_ncdf/initialize.f90
@@ -43,7 +43,7 @@ subroutine initialize(files, config)
integer :: ibdate, ibtime
integer, dimension(maxpoint) :: releases
real :: bndx, bndy, delx, dely
- integer :: numx, numy
+ integer :: numx, numy, numz
integer :: i, ierr
! read list of receptors
@@ -74,15 +74,27 @@ subroutine initialize(files, config)
else
filename = trim(files%path_flexpart)//trim(areldate)//'/header'
endif
+ inquire(file=trim(filename),exist=lexist)
else if ( config%lnested ) then
if ( trim(config%filefreq).eq.'month' ) then
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
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
print*, filename
- inquire(file=trim(filename),exist=lexist)
i = i + 1
end do
if ( .not.lexist ) then
@@ -91,9 +103,10 @@ subroutine initialize(files, config)
endif
write(*,*) 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely)
+ numx, numy, numz, bndx, bndy, delx, dely)
nxgrid = numx
nygrid = numy
+ nzgrid = numz
llx = bndx
lly = bndy
dx = delx
diff --git a/grid_to_ncdf/makefile b/grid_to_ncdf/makefile
index 75feccc..55f38ba 100644
--- a/grid_to_ncdf/makefile
+++ b/grid_to_ncdf/makefile
@@ -1,34 +1,41 @@
-F90 = gfortran
-LIBPATH = /usr/lib/
-INCPATH = /usr/include/
-LNK = -o
-CMPL = -c
-LIBS = -lnetcdf -lnetcdff -llapack
-#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 \
- -fbacktrace
-LDFLAGS = $(FFLAGS) -L$(LIBPATH) -I$(INCPATH) $(LIBS)
+F90 = /apps/sw/ubuntu22.04/gcc-11.4.0/gcc-13.2.0-tkesyophy2o6rjlzknndu3b4oyasvuqm/bin/gfortran
+LIBPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/lib/
+INCPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/include/
+
+# OPTIMIZATION LEVEL
+O_LEV = 0 # [0,1,2,3,g,s,fast]
+
+LIBS = -lm -DUSE_NCF -lnetcdff
+
+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
-SRCS = mod_var.f90 \
- mod_settings.f90 \
- mod_dates.f90 \
- mod_flexpart.f90 \
- initialize.f90 \
- grid_convert.f90 \
- init_convert.f90 \
- factor_convert.f90 \
- write_ncdf.f90 \
- read_reclist.f90 \
- grid_to_ncdf.f90
-
-OBJECTS = $(SRCS:.f90=.o)
-$(MAIN): $(OBJECTS) $(MODULES)
- $(F90) $(LNK) $(MAIN) $(OBJECTS) $(LIBS)
+MODULES = mod_var.o \
+ mod_settings.o \
+ mod_dates.o \
+ mod_flexpart.o
+
+OBJECTS = initialize.o \
+ grid_convert.o \
+ init_convert.o \
+ factor_convert.o \
+ write_ncdf.o \
+ read_reclist.o \
+ grid_to_ncdf.o
+
+%.o: %.mod
+
+$(MAIN): $(MODULES) $(OBJECTS)
+ +$(F90) -o $@ $(MODULES) $(OBJECTS) $(LDFLAGS)
+
%.o : %.f90
- $(F90) $(LDFLAGS) $(CMPL) $< -o $@
+ +$(F90) -c $(FFLAGS) $<
clean:
- rm -f $(OBJECTS) $(MODULES)
+ rm -f *.o *.mod
+
+.SUFFIXES = $(SUFFIXES) .f90
diff --git a/grid_to_ncdf/mod_flexpart.f90 b/grid_to_ncdf/mod_flexpart.f90
index 0480e61..82810c3 100644
--- a/grid_to_ncdf/mod_flexpart.f90
+++ b/grid_to_ncdf/mod_flexpart.f90
@@ -49,7 +49,7 @@ module mod_flexpart
! --------------------------------------------------
subroutine read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely)
+ numx, numy, numz, bndx, bndy, delx, dely)
use mod_var
@@ -59,7 +59,7 @@ module mod_flexpart
integer, intent(in out) :: numpoint
integer, intent(in out) :: ibdate, ibtime
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
character(len=29) :: flexversion
@@ -75,7 +75,7 @@ module mod_flexpart
real, dimension(:), allocatable :: xpoint, ypoint, zpoint1, zpoint2
real, dimension(:,:), allocatable :: xmass
integer, dimension(10) :: lage
- integer :: n, i, j, ierr
+ integer :: n, i, j, ierr, dummy
inquire ( file=trim(filename),exist=lexist )
if ( .not.lexist ) then
@@ -91,14 +91,14 @@ module mod_flexpart
read(100) ibdate, ibtime, flexversion
read(100) loutstep, loutaver, loutsample
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) nspec, nzgrid
+ read(100) nspec, numz
nspec=nspec/3
do n=1,nspec
- read(100) nzgrid, species(n)
- read(100) nzgrid, species(n)
- read(100) nzgrid, species(n)
+ read(100) dummy, species(n)
+ read(100) dummy, species(n)
+ read(100) numz, species(n)
end do
read(100) numpoint
diff --git a/grid_to_ncdf/run_grid2ncdf.sh b/grid_to_ncdf/run_grid2ncdf.sh
new file mode 100755
index 0000000..3d5a737
--- /dev/null
+++ b/grid_to_ncdf/run_grid2ncdf.sh
@@ -0,0 +1,16 @@
+#!/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
diff --git a/prep_flexpart/SETTINGS b/prep_flexpart/SETTINGS
index ad393e6..c6634ba 100644
--- a/prep_flexpart/SETTINGS
+++ b/prep_flexpart/SETTINGS
@@ -8,16 +8,20 @@
# ==================================================
# 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
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_options: /mypath/flexpart_options/
# path where to write flexpart output (root only, recname and month are appended)
path_output: /mypath/flexpart_output/
-# path to OH fields (if use OH chemistry)
-path_ohfield: /mypath/ohfields/
# path to observations
path_obs: /mypath/obs/
+# path to OH fields (FPv10 only)
+path_ohfield: /mypath/oh_fields/
# path where to write observation output
path_obsout: /mypath/obs_out/
# observation file format (one of obspack, wdcgg, noaa, basic, icos)
@@ -29,6 +33,18 @@ file_recept: /mypath/input/reclist.txt
# AVAILABLE files
file_avail: /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
# make flexpart releases only when there are observations (1) or at regular intervals (0)
@@ -44,6 +60,8 @@ laverage: 0
intaverage: 0
# flexpart file format (namelist = 1, ascii = 0)
lnamelist: 1
+# if observation output already exists - append (1) or make new files (0)
+lappendfile: 0
# COMMAND settings
# start date of simulation (yyyymmdd)
diff --git a/prep_flexpart/list_obsfiles.f90 b/prep_flexpart/list_obsfiles.f90
index ce62464..fea8dc5 100644
--- a/prep_flexpart/list_obsfiles.f90
+++ b/prep_flexpart/list_obsfiles.f90
@@ -27,22 +27,40 @@ subroutine list_obsfiles(settings)
use mod_var
use mod_settings
+ use mod_dates
implicit none
type (settings_t), intent(in) :: settings
integer :: nf, ierr
+ character(len=6) :: adate
+ integer :: jjjjmmdd, hhmiss
! list observation files
! ----------------------
- 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')
+ if ( settings%lappendfile.eq.1 ) then
+ call system('rm -f '//trim(settings%path_obsout)//'obsfiles.dat')
+ 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
- write(*,*) 'ERROR: cannot open obsfiles.txt'
+ write(*,*) 'ERROR: cannot open obsfiles.dat'
stop
endif
read(100,*,iostat=ierr) nfiles
diff --git a/prep_flexpart/main.f90 b/prep_flexpart/main.f90
index f48a0c2..114e804 100644
--- a/prep_flexpart/main.f90
+++ b/prep_flexpart/main.f90
@@ -41,6 +41,8 @@ program main
use mod_dates
use mod_settings
use mod_obs
+ use mod_prep_ageclass
+ use mod_prep_command
implicit none
@@ -50,6 +52,7 @@ program main
character(max_path_len) :: settings_file
real(kind=8) :: jd
+ character(len=6) :: adate
logical :: lexist
integer :: nobs, nr
integer :: jjjjmmdd, hhmiss, eomday
@@ -70,37 +73,47 @@ program main
jreldatei = juldate(settings%reldatei, 0)
jreldatef = juldate(settings%reldatef, 0)
- print*, 'datei: ', settings%datei, jdatei
- print*, 'datef: ', settings%datef, jdatef
- print*, 'reldatei: ', settings%reldatei, jreldatei
- print*, 'reldatef: ', settings%reldatef, jreldatef
+ write(*,*) 'datei: ', settings%datei, jdatei
+ write(*,*) 'datef: ', settings%datef, jdatef
+ write(*,*) 'reldatei: ', settings%reldatei, jreldatei
+ write(*,*) 'reldatef: ', settings%reldatef, jreldatef
! checks
if( jdatei.ne.jreldatei.or.jdatef.ne.jreldatef ) then
- print*, 'WARNING: different run and release times'
+ write(*,*) 'WARNING: different run and release times'
endif
if( (jreldatei.lt.jdatei).or.(jreldatei.gt.jdatef) ) then
- print*, 'ERROR: release date outside run time'
+ write(*,*) 'ERROR: release date outside run time'
stop
endif
if( (jreldatef.lt.jdatei).or.(jreldatef.gt.jdatef) ) then
- print*, 'ERROR: release date outside run time'
+ write(*,*) 'ERROR: release date outside run time'
stop
endif
inquire(file=settings%file_avail,exist=lexist)
if( .not.lexist ) then
- print*, 'ERROR: AVAILABLE file not found'
+ write(*,*) 'ERROR: AVAILABLE file not found'
stop
endif
if( settings%laverage.eq.0 ) then
settings%intaverage = 0.
- print*, 'WARNING: laverage is 0 but intaverage > 0'
- print*, 'setting intaverage to 0'
+ write(*,*) 'WARNING: laverage is 0 but intaverage > 0'
+ write(*,*) 'setting intaverage to 0'
endif
if( settings%lrelease.eq.0 ) then
- print*, 'FLEXPART releases at regular intervals'
+ write(*,*) 'FLEXPART releases at regular intervals'
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
! read receptor list
@@ -121,18 +134,20 @@ program main
call prep_pathnames(settings, jd, nr)
! 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
call prep_outgrid(settings, jd, nr)
! write AGECLASS
- if ( settings%lFPversion.eq.0 ) then
- write(*,*) 'Preparing ageclass for FP10', settings%lFPversion
- call prep_ageclass_FP10(settings, jd, nr)
- else
- write(*,*) 'Preparing ageclass for FP11', settings%lFPversion
- call prep_ageclass(settings, jd, nr)
+ if ( settings%FPversion.eq.10 ) then
+ call prep_ageclass_FP10(settings, jd, nr)
+ else if ( settings%FPversion.eq.11 ) then
+ call prep_ageclass(settings, jd, nr)
endif
! read and process observations
@@ -147,7 +162,8 @@ program main
else if ( settings%obsformat.eq.'basic' ) then
call read_basic(settings, jd, nr, nobs, obs)
else
- print*, 'ERROR: unknown observation file format'
+ write(*,*) 'ERROR: unknown observation file format'
+ stop
endif
! write RELEASES
diff --git a/prep_flexpart/makefile b/prep_flexpart/makefile
index a3a0e4b..67c16ed 100644
--- a/prep_flexpart/makefile
+++ b/prep_flexpart/makefile
@@ -18,15 +18,15 @@ MODULES = mod_var.o \
mod_dates.o \
mod_tools.o \
mod_obs.o \
- mod_strings.o
+ mod_strings.o \
+ mod_prep_ageclass.o \
+ mod_prep_command.o
OBJECTS = read_reclist.o \
list_obsfiles.o \
prep_pathnames.o \
- prep_command.o \
prep_outgrid.o \
- prep_ageclass.o \
- prep_ageclass_FP10.o \
+ prep_reagents.o \
prep_releases.o \
prep_releases_reg.o \
process_obs.o \
diff --git a/prep_flexpart/mod_prep_ageclass.f90 b/prep_flexpart/mod_prep_ageclass.f90
new file mode 100644
index 0000000..d93df8a
--- /dev/null
+++ b/prep_flexpart/mod_prep_ageclass.f90
@@ -0,0 +1,147 @@
+!---------------------------------------------------------------------------------------
+! 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
+
diff --git a/prep_flexpart/mod_prep_command.f90 b/prep_flexpart/mod_prep_command.f90
new file mode 100644
index 0000000..7020fc5
--- /dev/null
+++ b/prep_flexpart/mod_prep_command.f90
@@ -0,0 +1,359 @@
+!---------------------------------------------------------------------------------------
+! PREP_FLEXPART: mod_prep_command
+!---------------------------------------------------------------------------------------
+! 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_command
+!!
+!! Purpose: Prepares the options file COMMAND.
+!!
+!! Interface:
+!!
+!! Inputs
+!! settings - settings data structure
+!! jd - julian day for start of month
+!! nr - index to the receptors list
+!!
+!! Externals
+!! caldate
+!! juldate
+!!
+!---------------------------------------------------------------------------------------
+
+module mod_prep_command
+
+ use mod_var
+ use mod_settings
+ use mod_dates
+ use mod_tools
+
+ implicit none
+ private
+
+ public :: prep_command, prep_command_FP10
+
+ character(len=max_path_len) :: filename
+ character(len=6) :: adate
+ real(kind=8) :: jdf
+ integer :: jjjjmm, jjjjmmdd, hhmiss
+ integer :: datei, datef, eomday
+ integer :: ierr
+ integer :: ldirect ! runtime mode (currently only backwards: ldirect = -1)
+ integer :: ibdate, ibtime
+ integer :: iedate, ietime
+ integer :: loutstep
+ integer :: loutaver
+ integer :: loutsample
+ integer :: loutrestart
+ integer :: itsplit ! time constant for particle splitting (secs)
+ integer :: lsynctime ! synchronisation interval of flexpart (secs)
+ 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 :: iout ! output type (1 = concentration/residence time, 2 = mixing ratio, 3 = both, 4 = plume trajectory)
+ integer :: ipout ! particle dump (0 = never, 1 = every output interval, 2 = only at end)
+ integer :: lsubgrid ! use subgrid terrain effect parameterization (0 = no, 1 = yes)
+ integer :: lconvection ! use convection (0 = no, 1 = yes)
+ integer :: lturbulence ! FP11
+ integer :: lturbulence_meso ! FP11
+ integer :: nxshift ! FP11
+ integer :: maxthreadgrid ! FP11
+ integer :: maxfilesize ! FP11
+ integer :: logvertinterp ! FP11
+ integer :: lagespectra ! calculate age spectra (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 :: mdomainfill ! domain filling option (0 = no, 1 = yes)
+ integer :: mquasilag ! quasilagrangian mode to track particles (0 = no, 1 = yes)
+ 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 :: cblflag !
+ integer :: linversionout ! one grid_time file per release = 1, or per timestep (original format) = 0
+ integer :: iflux ! calculate fluxes (0 = no, 1 = yes)
+ integer :: ind_source
+ integer :: ind_receptor
+ integer :: nested_output
+ integer :: linit_cond
+ character(len=max_path_len) :: ohfields_path
+
+ contains
+
+ ! --------------------------------------------------
+ ! prep_command
+ ! --------------------------------------------------
+ !
+ ! Prepares COMMAND for FLEXPARTv11
+ ! --------------------------------------------------
+
+ subroutine prep_command(settings, jd, nr)
+
+ implicit none
+
+ type (settings_t), intent(in) :: settings
+ real(kind=8), intent(in) :: jd
+ integer, intent(in) :: nr
+
+ namelist /command/ &
+ ldirect, &
+ ibdate,ibtime, &
+ iedate,ietime, &
+ loutstep, &
+ loutaver, &
+ loutsample, &
+ lsynctime, &
+ ctl, &
+ ifine, &
+ iout, &
+ ipout, &
+ lsubgrid, &
+ lconvection, &
+ lturbulence, &
+ lturbulence_meso, &
+ nxshift, &
+ maxthreadgrid, &
+ maxfilesize, &
+ logvertinterp, &
+ loutrestart, &
+ lagespectra, &
+ ipin, &
+ ioutputforeachrelease, &
+ iflux, &
+ mdomainfill, &
+ ind_source, &
+ ind_receptor, &
+ mquasilag, &
+ nested_output, &
+ linit_cond, &
+ lnetcdfout, &
+ sfc_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 = 0
+ lsubgrid = 1
+ lconvection = 1
+ lturbulence = 1
+ lturbulence_meso = 0
+ lnetcdfout = 0
+ nxshift = 1
+ maxthreadgrid = 1
+ maxfilesize = 10000
+ logvertinterp = 0
+ loutrestart = -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
+ sfc_only = 1
+ cblflag = 0
+ linversionout = 1
+ ohfields_path = trim(settings%path_ohfield)
+
+ open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
+ write(100,nml=command)
+ close(100)
+
+ end subroutine prep_command
+
+ ! --------------------------------------------------
+ ! prep_command_FP10
+ ! --------------------------------------------------
+ !
+ ! Prepares COMMAND for FLEXPARTv10
+ ! --------------------------------------------------
+
+ subroutine prep_command_FP10(settings, jd, nr)
+
+ implicit none
+
+ type (settings_t), intent(in) :: settings
+ real(kind=8), intent(in) :: jd
+ integer, intent(in) :: nr
+
+ namelist /command/ &
+ ldirect, &
+ ibdate,ibtime, &
+ iedate,ietime, &
+ loutstep, &
+ loutaver, &
+ loutsample, &
+ itsplit, &
+ lsynctime, &
+ ctl, &
+ 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
+
+
diff --git a/prep_flexpart/mod_settings.f90 b/prep_flexpart/mod_settings.f90
index beb460e..abc5061 100644
--- a/prep_flexpart/mod_settings.f90
+++ b/prep_flexpart/mod_settings.f90
@@ -48,10 +48,17 @@ module mod_settings
character(len=max_name_len) :: windfield
character(len=max_name_len) :: obsformat
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 :: lrelease
+ integer :: lappendfile
real :: relfreq
integer :: laverage
real :: intaverage
@@ -296,6 +303,9 @@ module mod_settings
stop
endif
+ ! default settings
+ settings%lscratch = 1 ! use scratch drive
+
! loop over lines
read_loop: do
read (100, fmt='(A)', iostat=ierr) line
@@ -313,6 +323,9 @@ module mod_settings
if ( len (trim (line)) .eq. 0 ) cycle read_loop
! Path and file settings
+ identifier = "lscratch:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) settings%lscratch = int(cn)
identifier = "path_flexpart:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%path_flexpart = cc
@@ -349,17 +362,64 @@ module mod_settings
identifier = "windfield:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%windfield = cc
+ identifier = "FPversion:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) settings%FPversion = int(cn)
- ! General settings
- identifier = "lFPversion:"
+ ! Species settings
+ 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)
- 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:"
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 = "lappendfile:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) settings%lappendfile = int(cn)
identifier = "relfreq:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) settings%relfreq = real(cn)
@@ -374,7 +434,7 @@ module mod_settings
if ( match ) settings%lnamelist = int(cn)
identifier = "measerr:"
call read_content (line, identifier, cc, cn, cl, match)
- if ( match ) settings%measerr = int(cn)
+ if ( match ) settings%measerr = real(cn)
! COMMAND settings
identifier = "datei:"
@@ -485,6 +545,17 @@ module mod_settings
print*, 'path_flexpart: ',trim(settings%path_flexpart)
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
diff --git a/prep_flexpart/mod_var.f90 b/prep_flexpart/mod_var.f90
index 05fbe5d..606809e 100644
--- a/prep_flexpart/mod_var.f90
+++ b/prep_flexpart/mod_var.f90
@@ -36,6 +36,7 @@ module mod_var
real(kind=8) :: jreldatei, jreldatef
integer :: nrec
integer :: nfiles
+ integer :: nreagent
character(len=recname_len), dimension(:), allocatable :: recname
character(len=256), dimension(:), allocatable :: filelist
real, dimension(:), allocatable :: reclon, reclat, recalt
diff --git a/prep_flexpart/prep_ageclass.f90 b/prep_flexpart/prep_ageclass.f90
deleted file mode 100644
index fbeac57..0000000
--- a/prep_flexpart/prep_ageclass.f90
+++ /dev/null
@@ -1,99 +0,0 @@
-!---------------------------------------------------------------------------------------
-! PREP_FLEXPART: 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
-!---------------------------------------------------------------------------------------
-!
-!> 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
-!!
-!---------------------------------------------------------------------------------------
-
-subroutine prep_ageclass(settings, jd, nr)
-
- use mod_var
- use mod_settings
- use mod_dates
-
- implicit none
-
- type (settings_t), intent(in) :: settings
- real(kind=8), intent(in) :: jd
- integer, intent(in) :: nr
-
- 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
- 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)
-
- if( settings%lnamelist.eq.1 ) then
- ! use namelist file format, FP11
- write(100,nml=nage)
- 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
diff --git a/prep_flexpart/prep_command.f90 b/prep_flexpart/prep_command.f90
deleted file mode 100644
index 9a310e9..0000000
--- a/prep_flexpart/prep_command.f90
+++ /dev/null
@@ -1,241 +0,0 @@
-!---------------------------------------------------------------------------------------
-! PREP_FLEXPART: prep_command
-!---------------------------------------------------------------------------------------
-! 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_command
-!!
-!! Purpose: Prepares the options file COMMAND.
-!!
-!! Interface:
-!!
-!! Inputs
-!! settings - settings data structure
-!! jd - julian day for start of month
-!! nr - index to the receptors list
-!!
-!! Externals
-!! caldate
-!! juldate
-!!
-!---------------------------------------------------------------------------------------
-
-subroutine prep_command(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
- character(len=6) :: adate
- real(kind=8) :: jdf
- integer :: jjjjmm, jjjjmmdd, hhmiss
- integer :: datei, datef, eomday
- integer :: ierr
-
- integer :: ldirect ! runtime mode (currently only backwards: ldirect = -1)
- integer :: ibdate, ibtime
- integer :: iedate, ietime
-! integer :: partsplit
- integer :: loutstep
- integer :: loutaver
- integer :: loutsample
- integer :: loutrestart
-! integer :: itsplit ! time constant for particle splitting (secs)
- integer :: lsynctime ! synchronisation interval of flexpart (secs)
- 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 :: iout ! output type (1 = concentration/residence time, 2 = mixing ratio, 3 = both, 4 = plume trajectory)
- integer :: ipout ! particle dump (0 = never, 1 = every output interval, 2 = only at end)
- integer :: lsubgrid ! use subgrid terrain effect parameterization (0 = no, 1 = yes)
- integer :: lconvection ! use convection (0 = no, 1 = yes)
- integer :: lturbulence ! FP11
- integer :: lturbulence_meso ! FP11
- integer :: nxshift ! FP11
- integer :: maxthreadgrid ! FP11
- integer :: maxfilesize ! FP11
- integer :: logvertinterp ! FP11
- integer :: bcscheme ! FP11
- integer :: lagespectra ! calculate age spectra (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 :: mdomainfill ! domain filling option (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 :: lnetcdfout ! netcdf output (0 = no, 1 = yes)
- integer :: cblflag !
- integer :: linversionout ! one grid_time file per release = 1, or per timestep (original format) = 0
- integer :: iflux ! calculate fluxes (0 = no, 1 = yes)
- integer :: ind_source
- integer :: ind_receptor
- integer :: nested_output
- integer :: linit_cond
- character(len=max_path_len) :: ohfields_path
-
- namelist /command/ &
- ldirect, &
- ibdate,ibtime, &
- iedate,ietime, &
- loutstep, &
- loutaver, &
- loutsample, &
-! itsplit, &
- lsynctime, &
- ctl, &
- ifine, &
- iout, &
- ipout, &
- lsubgrid, &
- lconvection, &
- lturbulence, &
- lturbulence_meso, &
- nxshift, &
- maxthreadgrid, &
- maxfilesize, &
- logvertinterp, &
- bcscheme, &
- loutrestart, &
- lagespectra, &
- ipin, &
- ioutputforeachrelease, &
- iflux, &
- mdomainfill, &
- ind_source, &
- ind_receptor, &
- mquasilag, &
- nested_output, &
- linit_cond, &
- lnetcdfout, &
- sfc_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 = 0
- lsubgrid = 1
- lconvection = 1
- lturbulence = 1
- lturbulence_meso = 0
- lnetcdfout = 0
- nxshift = 1
- maxthreadgrid = 1
- maxfilesize = 10000
- logvertinterp = 0
- bcscheme = 1
- loutrestart = -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
- sfc_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
-
diff --git a/prep_flexpart/prep_outgrid.f90 b/prep_flexpart/prep_outgrid.f90
index 4be6291..8d55f9f 100644
--- a/prep_flexpart/prep_outgrid.f90
+++ b/prep_flexpart/prep_outgrid.f90
@@ -84,7 +84,7 @@ subroutine prep_outgrid(settings, jd, nr)
if ( settings%windfield.eq.'oper' ) then
outlon0 = -179.
else if ( settings%windfield.eq.'era5' ) then
- outlon0 = -179.5
+ outlon0 = -179.
else
write(*,*) 'ERROR: prep_outgrid: unknown windfield type'
stop
diff --git a/prep_flexpart/prep_pathnames.f90 b/prep_flexpart/prep_pathnames.f90
index 7e16d98..28fa2a1 100644
--- a/prep_flexpart/prep_pathnames.f90
+++ b/prep_flexpart/prep_pathnames.f90
@@ -59,16 +59,21 @@ subroutine prep_pathnames(settings, jd, nr)
open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
- 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//'/'
+ if ( settings%lscratch.eq.0 ) then
+ ! 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='(A)') trim(settings%file_avail)
- if( settings%lnested.eq.1 ) then
+ if( settings%file_availnest.ne.'' ) then
write(100,fmt='(A1)') '/'
write(100,fmt='(A)') trim(settings%file_availnest)
- write(100,fmt='(A)') '============================================'
- write(100,fmt='(A1)') ' '
endif
close(100)
diff --git a/prep_flexpart/prep_reagents.f90 b/prep_flexpart/prep_reagents.f90
new file mode 100644
index 0000000..9673f7d
--- /dev/null
+++ b/prep_flexpart/prep_reagents.f90
@@ -0,0 +1,73 @@
+!---------------------------------------------------------------------------------------
+! PREP_FLEXPART: prep_reagents
+!---------------------------------------------------------------------------------------
+! 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_reagents
+!!
+!! Purpose: Prepares the options file REAGENTS
+!!
+!! Interface:
+!!
+!! Inputs
+!! settings - settings data structure
+!! reagent - character string of reagent name (see SPECIES files)
+!!
+!---------------------------------------------------------------------------------------
+
+subroutine prep_reagents(settings,reagents,nr,jd)
+
+ use mod_var
+ use mod_dates
+ use mod_settings
+
+ implicit none
+
+ type (settings_t), intent(in) :: settings
+ character(len=6), dimension(nreagent), intent(in) :: reagents
+ integer, intent(in) :: nr
+ real(kind=8), intent(in) :: jd
+
+ character(len=max_path_len) :: filename, preag_path
+ character(len=max_name_len) :: preagent, punit
+ character(len=6) :: adate
+ integer :: phourly
+ integer :: jjjjmmdd, hhmiss, ierr, n
+
+ namelist /reagent_params/ &
+ preagent, &
+ preag_path, &
+ phourly, &
+ punit
+
+ call caldate(jd, jjjjmmdd, hhmiss)
+ write(adate,fmt='(I6.6)') jjjjmmdd/100
+ filename = trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/REAGENTS'
+ open(100,file=trim(filename),status='replace',action='write',iostat=ierr)
+ do n = 1, nreagent
+ if (reagents(n).ne."") then
+ preagent = trim(reagents(n))
+ preag_path = trim(settings%path_reagents(n))
+ phourly = settings%reagent_interp(n)
+ punit = trim(settings%reagent_units(n))
+ write(100,nml=reagent_params)
+ endif
+ end do
+ close(100)
+
+end subroutine prep_reagents
+
diff --git a/prep_flexpart/prep_releases.f90 b/prep_flexpart/prep_releases.f90
index a015180..2dd83d1 100644
--- a/prep_flexpart/prep_releases.f90
+++ b/prep_flexpart/prep_releases.f90
@@ -43,6 +43,7 @@ subroutine prep_releases(settings, jd, nr, nobs, obs)
use mod_dates
use mod_obs
use mod_tools
+ use mod_strings
implicit none
@@ -70,6 +71,23 @@ subroutine prep_releases(settings, jd, nr, nobs, obs)
integer :: zkind, parts
real :: mass
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/ &
nspec, &
@@ -86,20 +104,65 @@ subroutine prep_releases(settings, jd, nr, nobs, obs)
parts, &
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)
write(adate,fmt='(I6.6)') jjjjmmdd/100
! 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
-! sec - reads the last char if ( line(len_trim(line)-nchar+1:len_trim(line)) == trim(settings%species) ) ierr = 1
- if ( line(5:5+nchar-1) == trim(settings%species) ) ierr = 1
- end do
-! sec - read(line(9:11),*) spec
- read(line(1:3),*) spec
- write(aspec,fmt='(I3.3)') spec
+ if ( settings%FPversion.eq.10 ) then
+ inquire(file=trim(settings%path_flexpart)//'options/SPECIES/spec_overview',exist=lexist)
+ if ( .not.lexist ) then
+ write(*,*) 'ERROR: cannot find file /options/SPECIES/spec_overview'
+ stop
+ endif
+ 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
+ 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
nspec = 1
@@ -108,7 +171,7 @@ subroutine prep_releases(settings, jd, nr, nobs, obs)
! copy standard input files to options folder
pathname = trim(settings%path_options)//trim(recname(nr))//'/'//adate//'/options/'
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(settings%path_flexpart)//'options/*.dat'//' '//trim(pathname))
call system('cp '//trim(settings%path_flexpart)//'options/*.t'//' '//trim(pathname))
diff --git a/prep_flexpart/prep_releases_reg.f90 b/prep_flexpart/prep_releases_reg.f90
index 9cc5ff3..1b2f182 100644
--- a/prep_flexpart/prep_releases_reg.f90
+++ b/prep_flexpart/prep_releases_reg.f90
@@ -45,6 +45,7 @@ subroutine prep_releases_reg(settings, jd, nr)
use mod_settings
use mod_dates
use mod_tools
+ use mod_strings
implicit none
@@ -72,6 +73,22 @@ subroutine prep_releases_reg(settings, jd, nr)
integer :: zkind, parts
real :: mass
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/ &
nspec, &
@@ -88,6 +105,16 @@ subroutine prep_releases_reg(settings, jd, nr)
parts, &
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
call caldate(jd, jjjjmmdd, hhmiss)
write(adate,fmt='(I6.6)') jjjjmmdd/100
@@ -100,14 +127,51 @@ subroutine prep_releases_reg(settings, jd, nr)
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
+ if ( settings%FPversion.eq.10 ) then
+ inquire(file=trim(settings%path_flexpart)//'options/SPECIES/spec_overview',exist=lexist)
+ if ( .not.lexist ) then
+ write(*,*) 'ERROR: cannot fine file /options/SPECIES/spec_overview'
+ stop
+ endif
+ 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
+ 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
nspec = 1
diff --git a/prep_flexpart/read_basic.f90 b/prep_flexpart/read_basic.f90
index 9f29efa..52f25ae 100644
--- a/prep_flexpart/read_basic.f90
+++ b/prep_flexpart/read_basic.f90
@@ -63,6 +63,7 @@ subroutine read_basic(settings, jd, nr, nobs, obs)
character(len=max_path_len) :: rowfmt
character(len=200) :: line
character(len=200), dimension(50) :: args
+ character(len=6) :: adate
character(len=3) :: flag
real, dimension(12) :: temp
logical :: lexist
@@ -166,7 +167,7 @@ subroutine read_basic(settings, jd, nr, nobs, obs)
! no observations go to end
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
go to 20
endif
@@ -181,7 +182,17 @@ subroutine read_basic(settings, jd, nr, nobs, obs)
if ( nobs.eq.0 ) go to 20
! 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)
if( lexist ) then
! append to existing
@@ -212,8 +223,13 @@ subroutine read_basic(settings, jd, nr, nobs, obs)
close(100)
! write detailed station list file
- print*, 'lat, lon , alt: ',obs%lat(1), obs%lon(1), obs%alt(1)
- file_out = trim(settings%path_obsout)//'stnlist_detail.txt'
+ if ( settings%lappendfile.eq.1 ) then
+ ! 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)
write(rowfmt,fmt='(A,I1,A)') '(A',recname_len,',1X,F7.2,1X,F7.2,1X,F7.1)'
if( lexist ) then
diff --git a/prep_flexpart/read_icos.f90 b/prep_flexpart/read_icos.f90
index 1368edd..ded6283 100644
--- a/prep_flexpart/read_icos.f90
+++ b/prep_flexpart/read_icos.f90
@@ -64,6 +64,7 @@ subroutine read_icos(settings, jd, nr, nobs, obs)
character(len=200) :: line
character(len=200), dimension(20) :: args
character(len=1) :: flag
+ character(len=6) :: adate
real, dimension(16) :: temp
logical :: lexist
integer :: ierr
@@ -153,7 +154,6 @@ subroutine read_icos(settings, jd, nr, nobs, obs)
do n = 2, 11
read(args(n),*) temp(n)
end do
- print*, temp
yyyy = int(temp(3))
mm = int(temp(4))
dd = int(temp(5))
@@ -162,7 +162,6 @@ subroutine read_icos(settings, jd, nr, nobs, obs)
conc = temp(9)
err = temp(10)
flag = args(12)
- print*, flag
! calculate julian date
jjjjmmdd = yyyy*10000+mm*100+dd
hhmiss = hh*10000+mi*100
@@ -190,12 +189,11 @@ subroutine read_icos(settings, jd, nr, nobs, obs)
! no observations go to end
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
go to 20
endif
-
! allocate obs
call alloc_obs(nobs, obs)
@@ -206,7 +204,17 @@ subroutine read_icos(settings, jd, nr, nobs, obs)
if ( nobs.eq.0 ) go to 20
! write formatted obs to file
- file_out = trim(settings%path_obsout)//trim(recname(nr))//'_'//trim(settings%species)//'.txt'
+ 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)
if( lexist ) then
! append to existing
@@ -237,7 +245,13 @@ subroutine read_icos(settings, jd, nr, nobs, obs)
close(100)
! write detailed station list file
- file_out = trim(settings%path_obsout)//'stnlist_detail.txt'
+ if ( settings%lappendfile.eq.1 ) then
+ ! 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)
write(rowfmt,fmt='(A,I1,A)') '(A',recname_len,',1X,F7.2,1X,F7.2,1X,F7.1)'
if( lexist ) then
diff --git a/prep_flexpart/read_noaa.f90 b/prep_flexpart/read_noaa.f90
index 06f5943..e11007e 100644
--- a/prep_flexpart/read_noaa.f90
+++ b/prep_flexpart/read_noaa.f90
@@ -64,6 +64,7 @@ subroutine read_noaa(settings, jd, nr, nobs, obs)
character(len=200) :: line
character(len=200), dimension(50) :: args
character(len=3) :: flag
+ character(len=6) :: adate
real, dimension(13) :: temp
logical :: lexist
integer :: ierr
@@ -180,7 +181,7 @@ subroutine read_noaa(settings, jd, nr, nobs, obs)
! no observations go to end
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
go to 20
endif
@@ -243,7 +244,17 @@ subroutine read_noaa(settings, jd, nr, nobs, obs)
if ( nobs.eq.0 ) go to 20
! write formatted obs to file
- file_out = trim(settings%path_obsout)//trim(recname(nr))//'_'//trim(settings%species)//'.txt'
+ 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)
if( lexist ) then
! append to existing
@@ -274,8 +285,13 @@ subroutine read_noaa(settings, jd, nr, nobs, obs)
close(100)
! write detailed station list file
- print*, 'lat, lon , alt: ',obs%lat(1), obs%lon(1), obs%alt(1)
- file_out = trim(settings%path_obsout)//'stnlist_detail.txt'
+ if ( settings%lappendfile.eq.1 ) then
+ ! 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)
write(rowfmt,fmt='(A,I1,A)') '(A',recname_len,',1X,F7.2,1X,F7.2,1X,F7.1)'
if( lexist ) then
diff --git a/prep_flexpart/read_obspack.f90 b/prep_flexpart/read_obspack.f90
index dc11b94..3dbb911 100644
--- a/prep_flexpart/read_obspack.f90
+++ b/prep_flexpart/read_obspack.f90
@@ -63,6 +63,7 @@ subroutine read_obspack(settings, jd, nr, nobs, obs)
character(len=max_path_len) :: rowfmt
character(len=200) :: line
character(len=200), dimension(20) :: args
+ character(len=6) :: adate
character(len=4) :: version
real, dimension(16) :: temp
logical :: lexist
@@ -193,7 +194,7 @@ subroutine read_obspack(settings, jd, nr, nobs, obs)
! no observations go to end
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
go to 20
endif
@@ -208,7 +209,17 @@ subroutine read_obspack(settings, jd, nr, nobs, obs)
if ( nobs.eq.0 ) go to 20
! write formatted obs to file
- file_out = trim(settings%path_obsout)//trim(recname(nr))//'_'//trim(settings%species)//'.txt'
+ 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)
if( lexist ) then
! append to existing
@@ -239,7 +250,13 @@ subroutine read_obspack(settings, jd, nr, nobs, obs)
close(100)
! write detailed station list file
- file_out = trim(settings%path_obsout)//'stnlist_detail.txt'
+ if ( settings%lappendfile.eq.1 ) then
+ ! 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)
write(rowfmt,fmt='(A,I1,A)') '(A',recname_len,',1X,F7.2,1X,F7.2,1X,F7.1)'
if( lexist ) then
diff --git a/prep_flexpart/read_wdcgg.f90 b/prep_flexpart/read_wdcgg.f90
index 48499d1..0783191 100644
--- a/prep_flexpart/read_wdcgg.f90
+++ b/prep_flexpart/read_wdcgg.f90
@@ -63,8 +63,7 @@ subroutine read_wdcgg(settings, jd, nr, nobs, obs)
character(len=max_path_len) :: rowfmt
character(len=200) :: line
character(len=200), dimension(30) :: args
- character(len=10) :: adate
- character(len=5) :: atime
+ character(len=6) :: adate
character(len=6) :: utcstr
real, dimension(22) :: temp
logical :: lexist
@@ -176,7 +175,7 @@ subroutine read_wdcgg(settings, jd, nr, nobs, obs)
! no observations go to end
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
go to 20
endif
@@ -193,7 +192,17 @@ subroutine read_wdcgg(settings, jd, nr, nobs, obs)
if ( nobs.eq.0 ) go to 20
! write formatted obs to file
- file_out = trim(settings%path_obsout)//trim(recname(nr))//'_'//trim(settings%species)//'.txt'
+ 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)
if( lexist ) then
! append to existing
@@ -224,7 +233,13 @@ subroutine read_wdcgg(settings, jd, nr, nobs, obs)
close(100)
! write detailed station list file
- file_out = trim(settings%path_obsout)//'stnlist_detail.txt'
+ if ( settings%lappendfile.eq.1 ) then
+ ! 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)
write(rowfmt,fmt='(A,I1,A)') '(A',recname_len,',1X,F7.2,1X,F7.2,1X,F7.1)'
if( lexist ) then
diff --git a/prep_flexpart/sbatch_prep_flexpart.sh b/prep_flexpart/sbatch_prep_flexpart.sh
index 6fc859b..919cc6d 100755
--- a/prep_flexpart/sbatch_prep_flexpart.sh
+++ b/prep_flexpart/sbatch_prep_flexpart.sh
@@ -1,6 +1,6 @@
#!/bin/bash
#---------------------------------------------------
-partition=nilu
+partition=main
settings_files='./SETTINGS'
#---------------------------------------------------
diff --git a/prep_fluxes/makefile b/prep_fluxes/makefile
index 97883d0..e09dae2 100644
--- a/prep_fluxes/makefile
+++ b/prep_fluxes/makefile
@@ -1,34 +1,42 @@
-F90 = gfortran
-LIBPATH = /usr/lib/
-INCPATH = /usr/include/
-LNK = -o
-CMPL = -c
-LIBS = -lnetcdf -lnetcdff -llapack
-#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 \
- -fbacktrace
-LDFLAGS = $(FFLAGS) -L$(LIBPATH) -I$(INCPATH) $(LIBS)
+F90 = /apps/sw/ubuntu22.04/gcc-11.4.0/gcc-13.2.0-tkesyophy2o6rjlzknndu3b4oyasvuqm/bin/gfortran
+LIBPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/lib/
+INCPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/include/
+
+# OPTIMIZATION LEVEL
+O_LEV = 0 # [0,1,2,3,g,s,fast]
+
+LIBS = -lm -DUSE_NCF -lnetcdff
+
+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 = prep_flux
-SRCS = mod_var.f90 \
- mod_settings.f90 \
- mod_dates.f90 \
- sort.f90 \
- read_ncdf.f90 \
- write_ncdf.f90 \
- gridarea.f90 \
- convertgrid.f90 \
- calc_monthly.f90 \
- adjust_time.f90 \
- prep_flux.f90
-
-OBJECTS = $(SRCS:.f90=.o)
-$(MAIN): $(OBJECTS) $(MODULES)
- $(F90) $(LNK) $(MAIN) $(OBJECTS) $(LIBS)
+MODULES = mod_var.o \
+ mod_settings.o \
+ mod_dates.o
+
+OBJECTS = sort.o \
+ read_ncdf.o \
+ write_ncdf.o \
+ gridarea.o \
+ convertgrid.o \
+ calc_monthly.o \
+ adjust_time.o \
+ prep_flux.o
+
+%.o: %.mod
+
+$(MAIN): $(MODULES) $(OBJECTS)
+ +$(F90) -o $@ $(MODULES) $(OBJECTS) $(LDFLAGS)
+
%.o : %.f90
- $(F90) $(LDFLAGS) $(CMPL) $< -o $@
+ +$(F90) -c $(FFLAGS) $<
clean:
- rm -f $(OBJECTS) $(MODULES)
+ rm -f *.o *.mod
+
+.SUFFIXES = $(SUFFIXES) .f90
+
diff --git a/prep_regions/get_surfinf.f90 b/prep_regions/get_surfinf.f90
index bd81909..7ad8cd7 100644
--- a/prep_regions/get_surfinf.f90
+++ b/prep_regions/get_surfinf.f90
@@ -70,7 +70,7 @@ subroutine get_surfinf(files, config, recs, cdryair, obstimes, avetimes, surfinf
real(kind=8), dimension(maxngrid) :: gtime_tmp ! flux sensitivity time stamp (julian days) of first footprint per observation
real :: bndx, bndy, delx, dely
- integer :: numx, numy, xshift, new_grid_time, standard_grid_time
+ integer :: numx, numy, numz, xshift, new_grid_time, standard_grid_time
integer :: numpoint, release_nr, nr_footprints, ngrid, nretr
integer :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2) :: releases
@@ -129,7 +129,7 @@ subroutine get_surfinf(files, config, recs, cdryair, obstimes, avetimes, surfinf
if ( lexist ) then
print*, 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
else
print*, 'WARNING: cannot find ',filename
go to 10
@@ -140,7 +140,7 @@ subroutine get_surfinf(files, config, recs, cdryair, obstimes, avetimes, surfinf
if ( lexist ) then
print*, 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
else
print*, 'WARNING: cannot find ',filename
go to 10
diff --git a/prep_regions/initialize.f90 b/prep_regions/initialize.f90
index 1f11671..fcbafce 100644
--- a/prep_regions/initialize.f90
+++ b/prep_regions/initialize.f90
@@ -56,7 +56,7 @@ subroutine initialize(files, config)
integer :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2) :: releases
real :: bndx, bndy, delx, dely
- integer :: numx, numy, xshift, ngrid
+ integer :: numx, numy, numz, xshift, ngrid
integer :: i, ierr
! read receptor list file
@@ -79,9 +79,18 @@ subroutine initialize(files, config)
if ( config%satellite.and.nrec.eq.0 ) then
! satellite only
write(areldate,'(I8)') config%datei
- if ( .not.config%nested ) filename = trim(files%path_flexsat)//areldate//'/header'
- if ( config%nested ) filename = trim(files%path_flexsat)//areldate//'/header_nest'
- inquire(file=trim(filename),exist=lexist)
+ if ( .not.config%nested ) then
+ filename = trim(files%path_flexsat)//areldate//'/header'
+ inquire(file=trim(filename),exist=lexist)
+ else
+ filename = trim(files%path_flexsat)//areldate//'/header_nest'
+ inquire(file=trim(filename),exist=lexist)
+ if ( .not.lexist ) then
+ ! test if is FPv11 file
+ filename = trim(files%path_flexsat)//areldate//'/header_nest_grid_time'
+ inquire(file=trim(filename),exist=lexist)
+ endif
+ endif
else
! ground-based output
yyyymm=config%datei/100
@@ -89,10 +98,18 @@ subroutine initialize(files, config)
lexist = .false.
i = 1
do while ( (.not.lexist).and.(i.le.nrec) )
- if ( .not.config%nested ) filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header'
- if ( config%nested ) filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header_nest'
- print*, filename
- inquire(file=trim(filename),exist=lexist)
+ if ( .not.config%nested ) then
+ filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header'
+ inquire(file=trim(filename),exist=lexist)
+ else
+ 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
+ endif
i = i + 1
end do
endif
@@ -103,7 +120,7 @@ subroutine initialize(files, config)
print*, 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
nxgrid = numx
nygrid = numy
nxshift = xshift
diff --git a/prep_regions/mod_flexpart.f90 b/prep_regions/mod_flexpart.f90
index 3316434..5c10a04 100644
--- a/prep_regions/mod_flexpart.f90
+++ b/prep_regions/mod_flexpart.f90
@@ -49,7 +49,7 @@ module mod_flexpart
! --------------------------------------------------
subroutine read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
use mod_var
@@ -59,7 +59,7 @@ module mod_flexpart
integer, intent(in out) :: numpoint
integer, intent(in out) :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2), intent(in out) :: releases
- integer, intent(in out) :: numx, numy, xshift
+ integer, intent(in out) :: numx, numy, numz, xshift
real, intent(in out) :: bndx, bndy, delx, dely
character(len=29) :: flexversion
@@ -75,7 +75,7 @@ module mod_flexpart
real, dimension(:), allocatable :: xpoint, ypoint, zpoint1, zpoint2
real, dimension(:,:), allocatable :: xmass
integer, dimension(10) :: lage
- integer :: n, i, j, ierr
+ integer :: n, i, j, ierr, dummy
inquire ( file=trim(filename),exist=lexist )
if ( .not.lexist ) then
@@ -91,14 +91,14 @@ module mod_flexpart
read(100) ibdate, ibtime, flexversion
read(100) loutstep, loutaver, loutsample
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) nspec, nzgrid
+ read(100) nspec, numz
nspec=nspec/3
do n=1,nspec
- read(100) nzgrid, species(n)
- read(100) nzgrid, species(n)
- read(100) nzgrid, species(n)
+ read(100) dummy, species(n)
+ read(100) dummy, species(n)
+ read(100) numz, species(n)
end do
read(100) numpoint
diff --git a/prep_regions/mod_settings.f90 b/prep_regions/mod_settings.f90
index ddbdb4a..6025ddd 100644
--- a/prep_regions/mod_settings.f90
+++ b/prep_regions/mod_settings.f90
@@ -82,7 +82,7 @@ module mod_settings
type :: config_t
- character(len=3) :: spec ! species ('co2' or 'ghg')
+ character(len=3) :: spec ! species ('co2' or 'ghg' or 'aero')
integer :: datei ! start date (yyyymmdd)
integer :: datef ! end date (yyyymmdd)
logical :: nested ! use nested flexpart output (true or false)
diff --git a/prep_regions/prep_regions.f90 b/prep_regions/prep_regions.f90
index 1aa8d19..ad66a6d 100644
--- a/prep_regions/prep_regions.f90
+++ b/prep_regions/prep_regions.f90
@@ -130,15 +130,6 @@ program prep_regions
varname = 'regions'
call write_ncdf(files, filename, varname, nbox_xy)
- ! write the average fooptrint for all receptors
- open(50,file=trim(files%path_output)//'all_footprint_average.txt')
- do ix=1,nxregrid
- do jy=1,nyregrid
- write(50,*) reg_lon(ix), reg_lat(jy), surfinf(ix,jy)
- end do
- end do
- close(50)
-
! write lsm
allocate( mask(nxregrid,nyregrid), stat=ierr)
mask=lsm
diff --git a/prep_regions/run_prepregions.sh b/prep_regions/run_prepregions.sh
new file mode 100755
index 0000000..53c0e5f
--- /dev/null
+++ b/prep_regions/run_prepregions.sh
@@ -0,0 +1,27 @@
+#!/bin/bash
+#SBATCH -J prep_region
+#SBATCH -p debug
+#SBATCH -q debug
+#SBATCH -o prep_regions.out
+#SBATCH --time=4:00:00
+#SBATCH --mem=10000
+
+#---------------------------------------------------
+
+settings_file='/mypath/settings/SETTINGS_files'
+settings_config='/mypath/settings/SETTINGS_config'
+settings_region='/mypath/prep_regions/SETTINGS_regions'
+
+#---------------------------------------------------
+
+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 ./prep_regions ${settings_file} ${settings_config} ${settings_region}
+
+
diff --git a/prep_satellite/makefile b/prep_satellite/makefile
index 9b0516c..e54f331 100644
--- a/prep_satellite/makefile
+++ b/prep_satellite/makefile
@@ -1,40 +1,44 @@
-F90 = gfortran
-LIBPATH = /usr/lib/
-INCPATH = /usr/include/
-LNK = -o
-CMPL = -c
-LIBS = -lnetcdf -lnetcdff
-#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 \
- -fbacktrace
-LDFLAGS = $(FFLAGS) -L$(LIBPATH) -I$(INCPATH) $(LIBS)
+F90 = /apps/sw/ubuntu22.04/gcc-11.4.0/gcc-13.2.0-tkesyophy2o6rjlzknndu3b4oyasvuqm/bin/gfortran
+LIBPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/lib/
+INCPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/include/
+
+# OPTIMIZATION LEVEL
+O_LEV = 0 # [0,1,2,3,g,s,fast]
+
+LIBS = -lm -DUSE_NCF -lnetcdff
+
+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 = prep_satellite
-SRCS = mod_var.f90 \
- mod_dates.f90 \
- mod_settings.f90 \
- mod_save.f90 \
- gridarea.f90 \
- geodarea.f90 \
- prep_command.f90 \
- prep_outgrid.f90 \
- prep_ageclass.f90 \
- prep_pathnames.f90 \
- average.f90 \
- get_tropomi.f90 \
- get_bremen.f90 \
- get_oco2.f90 \
- main.f90
-
-
-OBJECTS = $(SRCS:.f90=.o)
-$(MAIN): $(OBJECTS) $(MODULES)
- $(F90) $(LNK) $(MAIN) $(OBJECTS) $(LIBS)
+MODULES = mod_var.o \
+ mod_dates.o \
+ mod_settings.o \
+ mod_save.o
+
+OBJECTS = gridarea.o \
+ geodarea.o \
+ prep_command.o \
+ prep_outgrid.o \
+ prep_ageclass.o \
+ prep_pathnames.o \
+ average.o \
+ get_tropomi.o \
+ get_bremen.o \
+ get_oco2.o \
+ main.o
+
+%.o: %.mod
+
+$(MAIN): $(MODULES) $(OBJECTS)
+ +$(F90) -o $@ $(MODULES) $(OBJECTS) $(LDFLAGS)
+
%.o : %.f90
- $(F90) $(LDFLAGS) $(CMPL) $< -o $@
+ +$(F90) -c $(FFLAGS) $<
clean:
- rm -f $(OBJECTS) $(MODULES)
-
+ rm -f *.o *.mod
+.SUFFIXES = $(SUFFIXES) .f90
diff --git a/prep_syndata/get_error_cov.f90 b/prep_syndata/get_error_cov.f90
index 9c05c9a..4a45ca3 100644
--- a/prep_syndata/get_error_cov.f90
+++ b/prep_syndata/get_error_cov.f90
@@ -46,7 +46,7 @@ subroutine get_error_cov(files, config)
filename = trim(files%path_output)//'evals.txt'
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
- write(logid,*) 'ERROR: cannot find '//trim(filename)
+ write(logid,*) 'ERROR get_error_cov: cannot find '//trim(filename)
stop
endif
@@ -59,7 +59,7 @@ subroutine get_error_cov(files, config)
end do
close(100)
neig = neig - 1
- print*, 'neig = ',neig
+ write(logid,*) 'get_error_cov: neig = ',neig
allocate ( evals(neig) )
allocate ( evecs(ndvar,neig) )
@@ -79,7 +79,7 @@ subroutine get_error_cov(files, config)
filename = trim(files%path_output)//'evecs.txt'
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
- write(logid,*) 'ERROR: cannot find '//trim(filename)
+ write(logid,*) 'ERROR get_error_cov: cannot find '//trim(filename)
stop
endif
@@ -99,13 +99,13 @@ subroutine get_error_cov(files, config)
stat_time(i) = juldatei + (i-1)*statres
end do
- print*, 'stat_time = ',stat_time
+ write(logid,*) 'get_error_cov: stat_time = ',stat_time
allocate ( cort(ntstate,ntstate) )
filename = trim(files%path_output)//'cort.txt'
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
- write(logid,*) 'ERROR: cannot find '//trim(filename)
+ write(logid,*) 'ERROR get_error_cov: cannot find '//trim(filename)
stop
endif
diff --git a/prep_syndata/get_flux.f90 b/prep_syndata/get_flux.f90
index 8554d8a..78f25c9 100644
--- a/prep_syndata/get_flux.f90
+++ b/prep_syndata/get_flux.f90
@@ -153,7 +153,7 @@ subroutine get_flux(config, files, flx, state)
else
filename = str_replace(files%filename_nee, 'YYYY', adate)
endif
- print*, 'get_flux: filename = ',filename
+ write(logid,*) 'get_flux: filename = ',filename
filename = trim(files%path_prior)//trim(filename)
inquire(file=trim(filename),exist=lexist)
if (.not.lexist) then
@@ -208,7 +208,7 @@ subroutine get_flux(config, files, flx, state)
state(:) = 0.
! endif
- print*, 'flx_time = ',flx_time
+ write(logid,*) 'get_flux: flx_time = ',flx_time
end subroutine get_flux
diff --git a/prep_syndata/get_obs.f90 b/prep_syndata/get_obs.f90
index 03c82d4..63daf7f 100644
--- a/prep_syndata/get_obs.f90
+++ b/prep_syndata/get_obs.f90
@@ -66,7 +66,7 @@ subroutine get_obs(config, files, obs, obserr)
filename = trim(files%path_output)//'monitor.txt'
inquire(file=trim(filename), exist=lexist)
if ( .not.lexist ) then
- write(logid,*) 'ERROR: cannot find '//trim(filename)
+ write(logid,*) 'ERROR get_obs: cannot find '//trim(filename)
stop
endif
@@ -85,12 +85,14 @@ subroutine get_obs(config, files, obs, obserr)
do i = 1, maxobs
read(100,fmt=rowfmt,iostat=ierr) recs, jjjjmmdd, hhmiss, jdate, conc, cini(i), cinipos, bkg(i), &
ghg(i), prior(i), post, cpri(i), cakpri(i), delta, obserr(i)
+ print*, recs, jjjjmmdd, hhmiss, jdate, conc, cini(i), cinipos, bkg(i), &
+ ghg(i), prior(i), post, cpri(i), cakpri(i), delta, obserr(i)
if ( ierr.ne.0 ) exit
end do
endif
close(100)
nobs = i - 1
- print*, 'nobs: ',nobs
+ write(logid,*) 'get_obs: nobs = ',nobs
! check obserr
do i = 1, nobs
diff --git a/prep_syndata/initialize.f90 b/prep_syndata/initialize.f90
index a2dbb47..1b06a9e 100644
--- a/prep_syndata/initialize.f90
+++ b/prep_syndata/initialize.f90
@@ -49,7 +49,7 @@ subroutine initialize(files, config)
real :: area
integer :: i, ix, jy
real :: bndx, bndy, delx, dely
- integer :: numx, numy, xshift
+ integer :: numx, numy, numz, xshift
integer :: ierr
! open log file
@@ -84,9 +84,18 @@ subroutine initialize(files, config)
lexist = .false.
i = 1
do while ( (.not.lexist).and.(i.le.nrec) )
- if ( .not.config%nested ) filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header'
- if ( config%nested ) filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header_nest'
- inquire(file=trim(filename),exist=lexist)
+ if ( .not.config%nested ) then
+ filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header'
+ inquire(file=trim(filename),exist=lexist)
+ else
+ filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header_nest'
+ inquire(file=trim(filename),exist=lexist)
+ if ( .not.lexist ) then
+ ! try if is FPv11 file
+ filename = trim(files%path_flexpart)//trim(recname(i))//'/'//adate//'/header_nest_grid_time'
+ inquire(file=trim(filename),exist=lexist)
+ endif
+ endif
i = i + 1
end do
endif ! satellite
@@ -96,7 +105,7 @@ subroutine initialize(files, config)
endif
write(logid,*) 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
! global domain variables
nxgrid = numx
nygrid = numy
@@ -179,15 +188,15 @@ subroutine initialize(files, config)
stop
endif
- print*, 'statres: ',statres
- print*, 'flxtres: ',flxtres
- print*, 'ntstate: ',ntstate
- print*, 'nt_flx: ',nt_flx
- print*, 'nflx: ',nflx
- print*, 'nt_nee: ',nt_nee
- print*, 'nd_nee: ',nd_nee
- print*, 'nday: ',nday
- print*, 'nyr: ',nyr
+ write(logid,*) 'statres: ',statres
+ write(logid,*) 'flxtres: ',flxtres
+ write(logid,*) 'ntstate: ',ntstate
+ write(logid,*) 'nt_flx: ',nt_flx
+ write(logid,*) 'nflx: ', nflx
+ write(logid,*) 'nt_nee: ',nt_nee
+ write(logid,*) 'nd_nee: ',nd_nee
+ write(logid,*) 'nday: ',nday
+ write(logid,*) 'nyr: ',nyr
! initialize grid variables
! -------------------------
@@ -267,6 +276,9 @@ subroutine initialize(files, config)
ndvar = nbox
endif
+ write(logid,*) 'nvar = ',nvar
+ write(logid,*) 'ndvar = ',ndvar
+
end subroutine initialize
diff --git a/prep_syndata/makefile b/prep_syndata/makefile
index 1fb0149..15bddca 100644
--- a/prep_syndata/makefile
+++ b/prep_syndata/makefile
@@ -1,39 +1,47 @@
-F90 = gfortran
-LIBPATH = /usr/lib/
-INCPATH = /usr/include/
-LNK = -o
-CMPL = -c
-LIBS = -lnetcdf -lnetcdff -llapack
-#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 \
- -fbacktrace -ffree-line-length-0 -fcheck=all -Wall
-LDFLAGS = $(FFLAGS) -L$(LIBPATH) -I$(INCPATH) $(LIBS)
+F90 = /apps/sw/ubuntu22.04/gcc-11.4.0/gcc-13.2.0-tkesyophy2o6rjlzknndu3b4oyasvuqm/bin/gfortran
+LIBPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/lib/
+INCPATH1 = /apps/sw/ubuntu22.04/gcc-13.2.0/netcdf-fortran-4.6.1-ubcs4l6pjwpgoeyvz32wpzyzykib6bwm/include/
+
+# OPTIMIZATION LEVEL
+O_LEV = 0 # [0,1,2,3,g,s,fast]
+
+LIBS = -llapack -lm -DUSE_NCF -lnetcdff
+
+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 -ffree-form -ffree-line-length-0 # -march=native
+
+LDFLAGS = $(FFLAGS) -L$(LIBPATH1) -Wl,-rpath,$(LIBPATH1) $(LIBS)
MAIN = prep_syndata
-SRCS = mod_var.f90 \
- mod_settings.f90 \
- mod_strings.f90 \
- mod_dates.f90 \
- mod_tools.f90 \
- mod_flexpart.f90 \
- mod_perturb.f90 \
- mod_save.f90 \
- read_reclist.f90 \
- read_ncdf.f90 \
- initialize.f90 \
- get_flux.f90 \
- get_error_cov.f90 \
- get_obs.f90 \
- calc_error_cov.f90 \
- prep_syndata.f90
-
-OBJECTS = $(SRCS:.f90=.o)
-$(MAIN): $(OBJECTS) $(MODULES)
- $(F90) $(LNK) $(MAIN) $(OBJECTS) $(LIBS)
+MODULES = mod_var.o \
+ mod_settings.o \
+ mod_strings.o \
+ mod_dates.o \
+ mod_tools.o \
+ mod_flexpart.o \
+ mod_perturb.o \
+ mod_save.o
+
+OBJECTS = read_reclist.o \
+ read_ncdf.o \
+ initialize.o \
+ get_flux.o \
+ get_error_cov.o \
+ get_obs.o \
+ calc_error_cov.o \
+ prep_syndata.o
+
+%.o: %.mod
+
+$(MAIN): $(MODULES) $(OBJECTS)
+ +$(F90) -o $@ $(MODULES) $(OBJECTS) $(LDFLAGS)
+
%.o : %.f90
- $(F90) $(LDFLAGS) $(CMPL) $< -o $@
+ +$(F90) -c $(FFLAGS) $<
clean:
- rm -f $(OBJECTS) $(MODULES)
+ rm -f *.o *.mod
+
+.SUFFIXES = $(SUFFIXES) .f90
+
diff --git a/prep_syndata/mod_flexpart.f90 b/prep_syndata/mod_flexpart.f90
index 97abe47..40c99e0 100644
--- a/prep_syndata/mod_flexpart.f90
+++ b/prep_syndata/mod_flexpart.f90
@@ -47,7 +47,7 @@ module mod_flexpart
! --------------------------------------------------
subroutine read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
use mod_var
@@ -57,7 +57,7 @@ module mod_flexpart
integer, intent(in out) :: numpoint
integer, intent(in out) :: ibdate, ibtime
integer, dimension(maxpoint), intent(in out) :: releases
- integer, intent(in out) :: numx, numy, xshift
+ integer, intent(in out) :: numx, numy, numz, xshift
real, intent(in out) :: bndx, bndy, delx, dely
character(len=29) :: flexversion
@@ -73,7 +73,7 @@ module mod_flexpart
real, dimension(:), allocatable :: xpoint, ypoint, zpoint1, zpoint2
real, dimension(:,:), allocatable :: xmass
integer, dimension(10) :: lage
- integer :: n, i, j, ierr
+ integer :: n, i, j, ierr, dummy
inquire ( file=trim(filename),exist=lexist )
if ( .not.lexist ) then
@@ -89,14 +89,14 @@ module mod_flexpart
read(100) ibdate, ibtime, flexversion
read(100) loutstep, loutaver, loutsample
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) nspec, nzgrid
+ read(100) nspec, numz
nspec=nspec/3
do n=1,nspec
- read(100) nzgrid, species(n)
- read(100) nzgrid, species(n)
- read(100) nzgrid, species(n)
+ read(100) dummy, species(n)
+ read(100) dummy, species(n)
+ read(100) numz, species(n)
end do
read(100) numpoint
diff --git a/prep_syndata/mod_perturb.f90 b/prep_syndata/mod_perturb.f90
index 2bdd1b7..309456b 100644
--- a/prep_syndata/mod_perturb.f90
+++ b/prep_syndata/mod_perturb.f90
@@ -70,12 +70,12 @@ module mod_perturb
work = cort
call dpotrf('L',ntstate,work,ntstate,info)
if ( info.gt.0 ) then
- write(logid,*) 'ERROR: lower Cholesky triangle'
+ write(logid,*) 'ERROR perturb_flux: lower Cholesky triangle'
stop
endif
call dpotri('L',ntstate,work,ntstate,info)
if ( info.gt.0 ) then
- write(logid,*) 'ERROR: Cholesky inverse'
+ write(logid,*) 'ERROR perturb_flux: Cholesky inverse'
stop
endif
! returns lower ('L') triangle of inverse of work
@@ -106,7 +106,7 @@ module mod_perturb
end do
end do
end do
- print*, 'min/max perturbation = ',minval(perturb), maxval(perturb)
+ write(logid,*) 'perturb_flux: min/max perturbation = ',minval(perturb), maxval(perturb)
! add perturbation to state vector
state(:) = state(:) + real(perturb(:),kind=4)
diff --git a/prep_syndata/mod_save.f90 b/prep_syndata/mod_save.f90
index 9ad81e3..4e56592 100644
--- a/prep_syndata/mod_save.f90
+++ b/prep_syndata/mod_save.f90
@@ -89,9 +89,9 @@ module mod_save
jy1 = minloc(gbl_lat,dim=1,mask=abs(gbl_lat-rlly).eq.minval(abs(gbl_lat-rlly)))
jy2 = minloc(gbl_lat,dim=1,mask=abs(gbl_lat-(rury-rdy)).eq.minval(abs(gbl_lat-(rury-rdy))))
- print*, 'save_flux: gbl_lon = ',gbl_lon
- print*, 'save_flux: gbl_lat = ',gbl_lat
- print*, 'save_flux: ix1,ix2,jy1,jy2 = ',ix1,ix2,jy1,jy2
+ write(logid,*) 'save_flux: gbl_lon = ',gbl_lon
+ write(logid,*) 'save_flux: gbl_lat = ',gbl_lat
+ write(logid,*) 'save_flux: ix1,ix2,jy1,jy2 = ',ix1,ix2,jy1,jy2
! space and time conversions
! --------------------------
@@ -155,7 +155,7 @@ module mod_save
do n = 1, nflx
nt = minloc(stat_time,dim=1,mask=abs(stat_time-flx_time(n)).eq.minval(abs(stat_time-flx_time(n))))
nt = min(nt,ntstate)
- print*, 'stat_time(nt), flx_time(n) = ',stat_time(nt), flx_time(n)
+ write(logid,*) 'save_flux: stat_time(nt), flx_time(n) = ',stat_time(nt), flx_time(n)
do jy = 1, nyregrid
do ix = 1, nxregrid
! if inc_ocean is true state vector includes all fluxes, if false need to add ocean fluxes
@@ -213,7 +213,7 @@ module mod_save
! convert to kg/m2/h
flx = flx*3600.
- print*, 'mod_save: min/max flx = ',minval(flx), maxval(flx)
+ write(logid,*) 'save_flux: min/max flx = ',minval(flx), maxval(flx)
! adjust time to days since 1-Jan-1900
flx_time = flx_time - juldate(19000101, 0) !+ 1d0
@@ -333,7 +333,7 @@ module mod_save
filename = trim(files%path_output)//'monitor.txt'
inquire(file=trim(filename), exist=lexist)
if ( .not.lexist ) then
- write(logid,*) 'ERROR: cannot find '//trim(filename)
+ write(logid,*) 'ERROR save_obs: cannot find '//trim(filename)
stop
endif
@@ -402,7 +402,7 @@ module mod_save
! third instance of delim
call split(string,"_",species,sep)
satname = trim(prefix)//'_'//trim(species)
- write(logid,*) 'prefix satellite file = ',satname
+ write(logid,*) 'save_obs: prefix satellite file = ',satname
exit
endif
end do
@@ -432,23 +432,25 @@ module mod_save
vmr(cnt) = obs(i)
endif
end do ! nobs
- ! copy original ncdf file and replace vmr with obs
- call caldate(jd, currentdate, currenttime)
- write(adate,fmt='(I8.8)') currentdate
- filename = 'retrieval_'//trim(satname)//'_'//adate//'.nc'
- print*, 'mod_save: filename = ',filename
- call system('cp '//trim(files%path_satobs)//trim(filename)//' '//trim(files%path_output)//trim(filename))
- call check( nf90_open(trim(files%path_output)//trim(filename),nf90_WRITE,ncid) )
- call check( nf90_inq_dimid(ncid,'retrieval',dimid) )
- call check( nf90_inquire_dimension(ncid,dimid,len=numretr) )
- ! check consistent number of retrievals
- if ( cnt.ne.numretr ) then
- write(logid,*) 'ERROR: inconsistent number of retrievals for day ',currentdate,': numretr, cnt = ',numretr,cnt
- stop
+ if ( cnt>0 ) then
+ ! copy original ncdf file and replace vmr with obs
+ call caldate(jd, currentdate, currenttime)
+ write(adate,fmt='(I8.8)') currentdate
+ filename = 'retrieval_'//trim(satname)//'_'//adate//'.nc'
+ write(logid,*) 'save_obs: filename = ',filename
+ call system('cp '//trim(files%path_satobs)//trim(filename)//' '//trim(files%path_output)//trim(filename))
+ call check( nf90_open(trim(files%path_output)//trim(filename),nf90_WRITE,ncid) )
+ call check( nf90_inq_dimid(ncid,'retrieval',dimid) )
+ call check( nf90_inquire_dimension(ncid,dimid,len=numretr) )
+ ! check consistent number of retrievals
+ if ( cnt.ne.numretr ) then
+ write(logid,*) 'ERROR save_obs: inconsistent number of retrievals for day ',currentdate,': numretr, cnt = ',numretr,cnt
+ stop
+ endif
+ call check( nf90_inq_varid(ncid,'vmr',varid) )
+ call check( nf90_put_var(ncid, varid, vmr(1:numretr)) )
+ call check( nf90_close(ncid) )
endif
- call check( nf90_inq_varid(ncid,'vmr',varid) )
- call check( nf90_put_var(ncid, varid, vmr(1:numretr)) )
- call check( nf90_close(ncid) )
! next day
jd = jd + 1d0
vmr(:) = 0.
diff --git a/prep_syndata/mod_var.f90 b/prep_syndata/mod_var.f90
index 5d6e10b..4a9b190 100644
--- a/prep_syndata/mod_var.f90
+++ b/prep_syndata/mod_var.f90
@@ -30,7 +30,7 @@ module mod_var
integer, parameter :: max_path_len=200 ! max character length of paths and files
integer, parameter :: max_name_len=50 ! max character length of variable names
- integer, parameter :: recname_len=3 ! length of receptor names
+ integer, parameter :: recname_len=7 ! length of receptor names
real, parameter :: numscale=1.e12 ! numeric scaling factor
integer, parameter :: maxpoint=1000 ! max number of releases per receptor per month
integer, parameter :: maxspec=10 ! max number of species in a flexpart run
diff --git a/prep_syndata/prep_syndata.f90 b/prep_syndata/prep_syndata.f90
index 673a13d..43835e6 100644
--- a/prep_syndata/prep_syndata.f90
+++ b/prep_syndata/prep_syndata.f90
@@ -60,9 +60,6 @@ program prep_syndata
stop 'ERROR: need to specify SETTINGS_config file'
endif
- write(*,*) 'Using: '//trim(settings_files)
- write(*,*) 'Using: '//trim(settings_config)
-
! initialization
! --------------
diff --git a/prep_syndata/run_prepsyndata.sh b/prep_syndata/run_prepsyndata.sh
new file mode 100755
index 0000000..c9075af
--- /dev/null
+++ b/prep_syndata/run_prepsyndata.sh
@@ -0,0 +1,24 @@
+#!/bin/bash
+#SBATCH -J prep_syndata
+#SBATCH -p main
+#SBATCH -o prep_syndata.out
+#SBATCH --time=4:00:00
+#SBATCH --mem=8000
+
+#---------------------------------------------------
+
+settings_file='/mypath/settings/SETTINGS_files'
+settings_config='/mypath/settings/SETTINGS_config'
+
+#---------------------------------------------------
+
+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 ./prep_syndata ${settings_file} ${settings_config}
+
diff --git a/settings/SETTINGS_aero_config b/settings/SETTINGS_aero_config
new file mode 100644
index 0000000..511cd86
--- /dev/null
+++ b/settings/SETTINGS_aero_config
@@ -0,0 +1,151 @@
+# ==================================================
+#
+# FLEXINVERT CONFIGURATION SETTINGS
+#
+# comment lines start with '#'
+# each parameter line starts with 'parameter name:'
+#
+# ==================================================
+
+# Run mode:
+# 0 = run forward model
+# 1 = run optimization
+# 2 = randomly perturb for MC
+run_mode: 1
+
+# Random number seed
+# (only used if run_mode = 2)
+seed: 100
+
+# Run dates: format yyyymmdd
+# datei = start date
+# datef = end date
+datei: 20120101
+datef: 20120131
+
+# Use lognormal distribution (true or false) (only for non-CO2 species)
+# if true need to use m1qn3 method
+lognormal: .false.
+
+# Truncation of eigenvalues of B (cut at trunc x largest eigenvalue)
+trunc: 1.e-4
+
+# Inversion method ('analytic', 'congrad' or 'm1qn3')
+method: congrad
+
+# Average/select flexpart releases (true or false)
+average_fp: .true.
+
+# Number of iterations
+# only used if method is 'congrad' or 'm1qn3'
+maxiter: 10
+
+# Optimize ocean boxes (true or false)
+inc_ocean: .false.
+
+# Optimize initial mixing ratios (true or false)
+# (not relevant for aerosols)
+opt_cini: .false.
+
+# Use spatial correlation in error covariance matrix (true or false)
+# if use regions based ecosystems then should be false
+spa_corr: .true.
+
+# Use best guess estimate from previous inversion (congrad only)
+# 0 = no prior best guess
+# 1 = use prior best guess (file must be specified in SETTINGS_files)
+prior_bg: 0
+
+# Restart a crashed run
+# for congrad/m1qn3 will pick-up from last iteration
+# for analytic will only use pre-calculated covariance matrix and boundary conditions
+# 0 = new run
+# 1 = restart crashed run
+restart: 0
+
+# Verbose output
+# only use for debugging small runs
+verbose: .true.
+
+# Use satellite observations (true or false)
+satellite: .false.
+
+# Use ground-based observations (true or false)
+ground: .true.
+
+# Species ("co2" or "ghg" or "aero")
+spec: aero
+
+# Coefficient to convert from grid units of ppt to ppb
+coeff: 1.0
+
+# Use nested flexpart output (true or false)
+nested: .false.
+
+# Inversion domain:
+# if using nested output it must match the inversion domain
+# w_edge_lon = lon of western edge of inversion grid
+# s_edge_lat = lat of southern edge inversion grid
+# e_edge_lon = lon of eastern edge of inversion grid
+# n_edge_lat = lat of northern inversion grid
+# xres = longitudinal grid resolution
+# yres = latitudinal grid resolution
+w_edge_lon: -15
+s_edge_lat: 30
+e_edge_lon: 35
+n_edge_lat: 75
+xres: 1.0
+yres: 1.0
+
+# Spatial aggregation of grid (true or false)
+regions: .true.
+
+# Select ground observations for day/night for low altitude/mountain sites
+selectobs: .true.
+
+# State vector (days):
+# temporal resolution at which NEE fluxes are optimized
+# must be an integer multiple of 1 day
+statres: 10.
+
+# prior fluxes:
+# nstep_flx = time step of other fluxes (integer hours, for monthly data use 720)
+nstep_flx: 720
+
+# Measurement error: unit same as obs
+measerr: 100.0
+
+# Initial mixing ratio error: fraction
+# (not relevant for aerosols)
+cinierr: 0.005
+
+# Prior flux error: fraction
+flxerr: 0.5
+
+# Lower limit flux error: unit (kg/m2/h)
+flxerr_ll: 1.e-8
+
+# Spatial correlation length for land: unit (km)
+sigma_land: 250.
+
+# Spatial correlation length for ocean: unit (km)
+sigma_ocean: 1000.
+
+# Temporal correlation length: unit (days)
+sigmatime: 30.
+
+# Total error for inversion domain (Tg/y)
+# globerr <= 0: prior error covariance matrix not scaled
+globerr: -1.
+
+# Settings for optimization of initial mixing ratios
+# (not relevant for aerosols)
+# comma separated list of northern edges of latitude bands
+cini_lat: -30.,0.,30.,90.
+# comma separated list of upper level of vertical bands (upper most level > outheight(nzgrid))
+#cini_alt: 2000., 10000., 55000.
+cini_alt: 55000.
+# time resolution for initial mixing ratio scalars (days)
+cinires: 30.
+
+
diff --git a/settings/SETTINGS_aero_files b/settings/SETTINGS_aero_files
new file mode 100644
index 0000000..3274512
--- /dev/null
+++ b/settings/SETTINGS_aero_files
@@ -0,0 +1,89 @@
+# ==================================================
+#
+# FLEXINVERT FILE SETTINGS
+#
+# comment lines start with '#'
+# each parameter line starts with 'parameter name:'
+#
+# Prior fluxes file dimensions:
+# time : days since 1900
+# latitude : latitude degrees north midpoints
+# longitude : longitude degrees east midpoints
+#
+# Regional mask file dimensions:
+# latitude : latitude degrees north midpoints
+# longitude : longitude degrees east midpoints
+#
+# Land-sea mask file dimensions:
+# latitude : latitude degrees north midpoints
+# longitude : longitude degrees east midpoints
+#
+# Initial concentrations file dimensions:
+# latitude : latitude degrees north south boundary
+# longitude : longitude degrees east west boundary
+#
+# Receptor list file:
+# ascii file containing one column list of
+# receptor names
+#
+# ==================================================
+
+# Paths:
+path_obs: /mypath/obs/
+path_prior: /mypath/flux/
+path_output: /mypath/output/
+path_flexpart: /mypath/flexpart_output/
+path_flexncdf: /mypath/flexpart_ncdf_out/
+
+# Prior best guess file (if prior_bg = 1 in SETTINGS_config)
+file_bg:
+
+# Suffix by which to identify observation files
+suffix: .txt
+
+# Log file:
+file_log: flexinvert.log
+
+# Regions mask file:
+# needs to correspond to flexpart domain and resolution (if using nested output needs to correspond to nest)
+file_regions: /mypath/output/regions_aero.nc
+varname_regs: regions
+lonname_regs: longitude
+latname_regs: latitude
+
+# Orography file:
+# needs to correspond to global flexpart domain and resolution
+file_orog: /mypath/input/elev.1-deg.nc
+varname_orog: data
+lonname_orog: lon
+latname_orog: lat
+
+# Land-sea mask file:
+# needs to correspond to flexpart domain and resolution (if using nested output needs to correspond to nest)
+file_lsm: /mypath/input/lsm_1x1.nc
+varname_lsm: lsm
+lonname_lsm: lon
+latname_lsm: lat
+
+# Global prior flux file: flux kg/m2/h
+# needs to correspond to flexpart domain and resolution
+# generic name with year YYYY
+filename_flx: AERO_TOTAL_YYYY_10x10.nc
+varname_flx: emis
+lonname_flx: longitude
+latname_flx: latitude
+timename_flx: time
+
+# Regional prior flux file: flux kg/m2/h
+# only needed if using nested output (needs to correspond to nest)
+# generic name with year YYYY
+filenest_flx:
+varnest_flx:
+lonnest_flx:
+latnest_flx:
+timenest_flx:
+
+# Receptor list file:
+file_recept: /mypath/input/reclist_aero.txt
+
+
diff --git a/settings/SETTINGS_co2_config b/settings/SETTINGS_co2_config
index c35647a..a80f8fd 100644
--- a/settings/SETTINGS_co2_config
+++ b/settings/SETTINGS_co2_config
@@ -103,6 +103,9 @@ yres: 1.0
# Spatial aggregation of grid (true or false)
regions: .false.
+# Select ground observations for day/night for low altitude/mountain sites
+selectobs: .true.
+
# State vector time resolution:
# time resolution at which NEE fluxes are optimized
# statres = averaging interval over 1 or more days (given in days)
diff --git a/settings/SETTINGS_ghg_config b/settings/SETTINGS_ghg_config
index 225cc30..901327c 100644
--- a/settings/SETTINGS_ghg_config
+++ b/settings/SETTINGS_ghg_config
@@ -102,6 +102,9 @@ yres: 1.0
# Spatial aggregation of grid (true or false)
regions: .true.
+# Select ground observations for day/night for low altitude/mountain sites
+selectobs: .true.
+
# State vector (days):
# temporal resolution at which NEE fluxes are optimized
# must be an integer multiple of 1 day
diff --git a/source/average_fp.f90 b/source/average_fp.f90
index 4600328..1a52cc6 100755
--- a/source/average_fp.f90
+++ b/source/average_fp.f90
@@ -79,7 +79,7 @@ subroutine average_fp(files, config, obs)
integer :: ngrid, nretr, ierr
real :: bndx, bndy, delx, dely
- integer :: numx, numy, xshift, new_grid_time, standard_grid_time
+ integer :: numx, numy, numz, xshift, new_grid_time, standard_grid_time
integer :: numpoint, release_nr, nr_footprints
integer :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2) :: releases
@@ -127,12 +127,12 @@ subroutine average_fp(files, config, obs)
if ( i.eq.1 .or. month.ne.prevmonth ) then
write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
else
if(trim(obs%recs(i)).ne.trim(obs%recs(i-1)))then
write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
endif
endif
diff --git a/source/calc_conc.f90 b/source/calc_conc.f90
index c3e1a0f..b70c602 100644
--- a/source/calc_conc.f90
+++ b/source/calc_conc.f90
@@ -76,8 +76,8 @@ subroutine calc_conc(config, fluxes, obs, ngrid, gtime, hnest, hbkg, iobs, ix1,
flxbg = fluxes%flxnee(:,:,ind) + &
fluxes%flxff(:,:,ind) + &
fluxes%flxocn(:,:,ind)
- else if (( trim(config%spec).eq.'ghg' ) .or. ( trim(config%spec).eq.'bca' )) then
- ! GHG or BC aerosol species
+ else if (( trim(config%spec).eq.'ghg' ) .or. ( trim(config%spec).eq.'aero' )) then
+ ! GHG or aerosol species
flxbg = fluxes%flx(:,:,ind)
endif
@@ -129,7 +129,7 @@ subroutine calc_conc(config, fluxes, obs, ngrid, gtime, hnest, hbkg, iobs, ix1,
obs%ocn(iobs) = obs%ocn(iobs) + sum(hnest(:,:,n) * fluxes%flxocn(ix1:ix2,jy1:jy2,ilo))
endif
else
- ! GHG species
+ ! GHG and aerosol species
if ( config%nested ) then
! nested domain
obs%ghg(iobs) = obs%ghg(iobs) + sum(hnest(:,:,n) * fluxes%flx_nest(:,:,ilo))
@@ -138,7 +138,7 @@ subroutine calc_conc(config, fluxes, obs, ngrid, gtime, hnest, hbkg, iobs, ix1,
obs%ghg(iobs) = obs%ghg(iobs) + sum(hnest(:,:,n) * fluxes%flx(ix1:ix2,jy1:jy2,ilo))
endif
if ( config%lognormal.and..not.config%inc_ocean ) then
- ! since for lognormal fluxes are optimized (not offsets) if ocean fluxes are not optimized need
+ ! since for lognormal scalars of fluxes are optimized, if ocean fluxes are not optimized need
! to still account for influence on mixing ratios
obs%ocn(iobs) = obs%ocn(iobs) + sum(hnest(:,:,n) * fluxes%flxocn(:,:,ilo))
endif
@@ -151,9 +151,12 @@ subroutine calc_conc(config, fluxes, obs, ngrid, gtime, hnest, hbkg, iobs, ix1,
if ( trim(config%spec).eq.'co2' ) then
! CO2
obs%err(iobs) = sqrt(obs%measerr(iobs)**2 + (sum(obs%cini(iobs,:))*config%cinierr)**2 + bkgerr + ffferr)
- else
- ! GHG
+ else if ( trim(config%spec).eq.'ghg' ) then
+ ! GHG
obs%err(iobs) = sqrt(obs%measerr(iobs)**2 + (sum(obs%cini(iobs,:))*config%cinierr)**2 + bkgerr)
+ else if ( trim(config%spec).eq.'aero' ) then
+ ! aerosols
+ obs%err(iobs) = sqrt(obs%measerr(iobs)**2 + bkgerr)
endif
end subroutine calc_conc
diff --git a/source/init_aero.f90 b/source/init_aero.f90
new file mode 100644
index 0000000..7563d20
--- /dev/null
+++ b/source/init_aero.f90
@@ -0,0 +1,472 @@
+!---------------------------------------------------------------------------------------
+! FLEXINVERT: init_aero
+!---------------------------------------------------------------------------------------
+! 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
+!---------------------------------------------------------------------------------------
+!
+!> init_aero
+!! Initialize aerosol variables
+!!
+!! Purpose: Reads prior fluxes and observations of an aerosol and calculates the
+!! state and observation vectors as well as their uncertainties.
+!!
+!! Interface:
+!!
+!! Inputs
+!! files - file data structure
+!! config - configuration settings data structure
+!! fluxes - flux data structure
+!! obs - observation data structure
+!! covar - error covariance data structure
+!!
+!! Externals
+!! alloc_fluxes
+!! alloc_states
+!! caldate
+!! read_ncdf
+!! read_bg
+!! correl
+!! error_cov
+!! read_obs
+!! assign_obs
+!! perturb_state
+!! perturb_obs
+!!
+!---------------------------------------------------------------------------------------
+
+subroutine init_aero(files, config, fluxes, obs, states, covar)
+
+ use mod_var
+ use mod_settings
+ use mod_dates
+ use mod_tools
+ use mod_strings
+ use mod_fluxes
+ use mod_obs
+ use mod_states
+ use mod_covar
+ use mod_perturb
+ use mod_save, only : save_for_restart
+
+ implicit none
+
+ type (files_t), intent (in) :: files
+ type (config_t), intent (in) :: config
+ type (fluxes_t), intent (in out) :: fluxes
+ type (obs_t), intent (in out) :: obs
+ type (states_t), intent (in out) :: states
+ type (covar_t), intent (in out) :: covar
+
+ character(max_path_len) :: filename
+ character(max_path_len) :: rowfmt
+ character(len=4) :: adate
+ character(len=2) :: amonth
+ logical :: lexist
+ real, dimension(:,:,:), allocatable :: flx, flxghg
+ real, dimension(:,:), allocatable :: datain, dataout
+ real, dimension(nbox,ntstate) :: flx_box
+ real, dimension(nbox,ntstate) :: err_box
+ real, dimension(nbox) :: xerr
+ real, dimension(nbox,nbox) :: corr
+ integer :: yyyymmdd, hhmiss
+ real :: area
+ real(kind=8) :: jd
+ integer :: ix, jy, ix0, jy0, ix1, ix2, jy1, jy2
+ integer :: i, n, nb, nread, num, eomday
+ integer :: ierr
+ real(kind=8), dimension(:), allocatable :: jdate, timeghg
+ character(len=recname_len), dimension(:), allocatable:: recs
+ real :: conc
+ character(len=20) :: dump
+ real, parameter :: smallnum=1.e-15
+
+ ! global prior fluxes
+ ! -------------------
+
+ ! if nested run global fluxes are only used to calculate the contribution to mixing ratios
+ ! from fluxes outside the inversion domain, otherwise, global fluxes are used to
+ ! calculate the contribution to mixing ratios and as prior fluxes inside the domain.
+
+ call alloc_fluxes(config, fluxes)
+ fluxes%flx(:,:,:) = 0.
+ fluxes%flxocn(:,:,:) = 0.
+
+ ! loop over prior flux files
+
+ allocate( flxghg(nxgrid,nygrid,max(1,(nday*24/nt_flx))), stat = ierr )
+ if ( ierr.ne.0 ) stop 'ERROR init_aero: not enough memory'
+ allocate( timeghg(max(1,nday*24/nt_flx)) )
+
+ n = 0
+ jd = juldatei
+ do while ( jd.le.juldatef )
+ call caldate(jd, yyyymmdd, hhmiss)
+ write(adate,'(I4.4)') yyyymmdd/10000
+ write(amonth,'(I2.2)') (yyyymmdd-(yyyymmdd/10000)*10000)/100
+ eomday = calceomday(yyyymmdd/100)
+ nread = eomday*24/nt_flx
+ if ( nread.eq.0 ) nread = 1
+ print*, 'init_aero: nread = ',nread
+ allocate( flx(nxgrid,nygrid,nread), stat = ierr )
+ if ( ierr.ne.0 ) stop 'ERROR init_aero: not enough memory'
+ filename = str_replace(files%filename_flx, 'YYYY', adate)
+ if ( index(filename, 'MM').ne.0 ) then
+ filename = str_replace(filename, 'MM', amonth)
+ endif
+ filename = trim(files%path_prior)//trim(filename)
+ write(logid,*) 'Reading prior fluxes :'//trim(filename)
+ inquire(file=trim(filename),exist=lexist)
+ if (.not.lexist) then
+ write(logid,*) 'ERROR: cannot find '//trim(filename)
+ stop
+ endif
+ call read_ncdf(filename, &
+ files%varname_flx, &
+ files%lonname_flx, &
+ files%latname_flx, &
+ files%timename_flx,&
+ llx, lly, &
+ nxgrid, nygrid, &
+ jd, nread, num, &
+ flx)
+ do i = 1, nread
+ if ( ((n+i).gt.(nday*24/nt_flx)).and.(nday*24/nt_flx.ne.0) ) exit
+ flxghg(:,:,n+i) = flx(:,:,i)
+ timeghg(n+i) = jd + dble((i-1)*nt_flx)/24d0
+ end do
+ deallocate( flx )
+ n = n + nread
+ jd = jd + dble(eomday)
+ end do
+ n = min(n, nday*24/nt_flx)
+ n = max(n, 1)
+ print*, 'init_aero: n = ',n
+ print*, 'init_aero: timeghg = ',timeghg
+
+ ! convert to kg/m2/s
+ flxghg = flxghg/3600.
+ ! apply numerical scaling
+ flxghg = flxghg*numscale
+
+ ! average/interpolate to state vector timestamp
+ allocate( datain(nxgrid,n) )
+ allocate( dataout(nxgrid,ntstate) )
+ do i = 1, ntstate
+ fluxes%time(i) = juldatei + dble(i-1)*dble(statres)
+ enddo
+ if ( statres.gt.real(nt_flx/24) ) then
+ ! average flux
+ do jy = 1, nygrid
+ datain = flxghg(:,jy,:)
+ call average(nxgrid, n, timeghg, datain, ntstate, fluxes%time, dataout)
+ fluxes%flx(:,jy,:) = dataout(:,:)
+ end do
+ else if ( statres.lt.real(nt_flx/24) ) then
+ ! interpolate flux
+ do jy = 1, nygrid
+ datain = flxghg(:,jy,:)
+ call interp(nxgrid, n, timeghg, datain, ntstate, fluxes%time, dataout)
+ fluxes%flx(:,jy,:) = dataout(:,:)
+ end do
+ else if ( statres.eq.real(nt_flx/24) ) then
+ fluxes%flx(:,:,:) = flxghg
+ endif
+ deallocate( datain )
+ deallocate( dataout )
+ deallocate( flxghg )
+ deallocate( timeghg )
+
+ ! ocean fluxes for lognormal case when inc_ocean is false
+ ! cover only inversion domain
+ if ( config%lognormal.and..not.config%inc_ocean.and..not.config%nested ) then
+ ! use land-sea mask to set all land fluxes to zero
+ ix1 = int((rllx - llx + dx)/dx)
+ ix2 = int((rurx - llx)/dx)
+ jy1 = int((rlly - lly + dy)/dy)
+ jy2 = int((rury - lly)/dy)
+ do n = 1, ntstate
+ fluxes%flxocn(:,:,n) = abs(lsm(:,:) - 1.)*fluxes%flx(ix1:ix2,jy1:jy2,n)
+ end do
+ endif
+
+ ! regional prior fluxes
+ ! ---------------------
+
+ if ( config%nested ) then
+
+ allocate( flxghg(nxregrid,nyregrid,(max(1,nday*24/nt_flx))), stat = ierr )
+ if ( ierr.ne.0 ) stop 'ERROR init_aero: not enough memory'
+ allocate( timeghg(max(1,nday*24/nt_flx)) )
+
+ n = 0
+ jd = juldatei
+ do while ( jd.le.juldatef )
+ call caldate(jd, yyyymmdd, hhmiss)
+ write(adate,'(I4.4)') yyyymmdd/10000
+ write(amonth,'(I2.2)') (yyyymmdd-(yyyymmdd/10000)*10000)/100
+ eomday = calceomday(yyyymmdd/100)
+ nread = eomday*24/nt_flx
+ if ( nread.eq.0 ) nread = 1
+ print*, 'init_aero: reg nread = ',nread
+ allocate( flx(nxregrid,nyregrid,nread), stat = ierr )
+ if ( ierr.ne.0 ) stop 'ERROR init_aero: not enough memory'
+ filename = str_replace(files%filenest_flx, 'YYYY', adate)
+ if ( index(filename, 'MM').ne.0 ) then
+ filename = str_replace(filename, 'MM', amonth)
+ endif
+ filename = trim(files%path_prior)//trim(filename)
+ write(logid,*) 'Reading prior fluxes :'//trim(filename)
+ inquire(file=trim(filename),exist=lexist)
+ if (.not.lexist) then
+ write(logid,*) 'ERROR: cannot find '//trim(filename)
+ stop
+ endif
+ call read_ncdf(filename, &
+ files%varnest_flx, &
+ files%lonnest_flx, &
+ files%latnest_flx, &
+ files%timenest_flx,&
+ rllx, rlly, &
+ nxregrid, nyregrid, &
+ jd, nread, num, &
+ flx)
+ do i = 1, nread
+ if ( ((n+i).gt.(nday*24)/nt_flx).and.((nday*24)/nt_flx.ne.0) ) exit
+ timeghg(n+i) = jd + dble((i-1)*nt_flx)/24d0
+ flxghg(:,:,n+i) = flx(:,:,i)
+ end do
+ deallocate( flx )
+ n = n + nread
+ jd = jd + dble(eomday)
+ end do
+ n = min(n, (nday*24)/nt_flx)
+ n = max(n, 1)
+ print*, 'init_aero: reg n = ',n
+ print*, 'init_aero: reg timeghg = ',timeghg
+
+ ! convert to kg/m2/s
+ flxghg = flxghg/3600.
+ ! apply numerical scaling
+ flxghg = flxghg*numscale
+
+ ! average/interpolate to state vector timestamp
+ allocate( datain(nxregrid,n) )
+ allocate( dataout(nxregrid,ntstate) )
+ if ( statres.gt.real(nt_flx/24) ) then
+ ! average flux
+ do jy = 1, nyregrid
+ datain = flxghg(:,jy,:)
+ call average(nxregrid, n, timeghg, datain, ntstate, fluxes%time, dataout)
+ fluxes%flx_nest(:,jy,:) = dataout(:,:)
+ end do
+ else if ( statres.lt.real(nt_flx/24) ) then
+ ! interpolate flux
+ do jy = 1, nyregrid
+ datain = flxghg(:,jy,:)
+ call interp(nxregrid, n, timeghg, datain, ntstate, fluxes%time, dataout)
+ fluxes%flx_nest(:,jy,:) = dataout(:,:)
+ end do
+ else if ( statres.eq.real(nt_flx/24) ) then
+ fluxes%flx_nest(:,:,:) = flxghg
+ endif
+ deallocate( datain )
+ deallocate( dataout )
+ deallocate( flxghg )
+ deallocate( timeghg )
+
+ ! ocean fluxes for lognormal case when inc_ocean is false
+ if ( config%lognormal.and..not.config%inc_ocean ) then
+ ! use land-sea mask to set all land fluxes to zero
+ do n = 1, ntstate
+ fluxes%flxocn(:,:,n) = abs(lsm(:,:) - 1.)*fluxes%flx_nest(:,:,n)
+ end do
+ endif
+
+ ! nest -> use nested fluxes for state vector
+ allocate( flx(nxregrid,nyregrid,ntstate) )
+ flx(:,:,:) = fluxes%flx_nest
+
+ else
+
+ ! no nest -> use global fluxes for state vector
+ allocate( flx(nxgrid,nygrid,ntstate) )
+ flx(:,:,:) = fluxes%flx
+
+ endif
+
+ ! prior state vector
+ ! ------------------
+
+ ! calculate area-weighted mean flux and uncertainty for each box
+
+ if ( config%nested ) then
+ ix0 = 1
+ jy0 = 1
+ else
+ ix0 = minloc(gbl_lon,dim=1,mask=abs(gbl_lon-rllx).eq.minval(abs(gbl_lon-rllx)))
+ jy0 = minloc(gbl_lat,dim=1,mask=abs(gbl_lat-rlly).eq.minval(abs(gbl_lat-rlly)))
+ endif
+ write(logid,*) 'init_aero: ix0, jy0 = ',ix0,jy0
+
+ flx_box(:,:) = 0.
+ err_box(:,:) = 0.
+ do n = 1, ntstate
+ do ix = 1, nxregrid
+ do jy = 1, nyregrid
+ area = grid_area(reg_lat(jy)+0.5*rdy, rdy, rdx)
+ ! if inc_ocean flx_box contains all fluxes, if false only land fluxes
+ if ( nbox_xy(ix,jy).gt.0 ) then
+ flx_box(nbox_xy(ix,jy),n) = flx_box(nbox_xy(ix,jy),n) + &
+ flx(ix0+ix-1,jy0+jy-1,n) * area
+ err_box(nbox_xy(ix,jy),n) = err_box(nbox_xy(ix,jy),n) + &
+ flx(ix0+ix-1,jy0+jy-1,n) * config%flxerr * area
+ endif
+ end do
+ end do
+ flx_box(:,n) = flx_box(:,n)/area_box
+ err_box(:,n) = err_box(:,n)/area_box
+ end do
+ fluxes%flx_box(:,:,1) = flx_box(:,:)
+
+ ! set minimum error
+ do n = 1, ntstate
+ do nb = 1, nbox
+ err_box(nb,n) = max(config%flxerr_ll*numscale/3600.,err_box(nb,n))
+ end do
+ end do
+
+ ! verbose output
+ ! --------------
+
+ if ( config%verbose ) then
+ ! check prior fluxes read correctly (1-time dimension only)
+ write(rowfmt,'(A,I6,A)') '(',nxregrid,'(E11.4,1X))'
+ open(100,file=trim(files%path_output)//'flx_prior.txt',status='replace',action='write',iostat=ierr)
+ do n = 1, nyregrid
+ write(100,fmt=rowfmt) flx(ix0:(ix0+nxregrid-1),jy0+n-1,1)
+ end do
+ close(100)
+ ! check regridding to variable resolution
+ write(rowfmt,'(A,I6,A)') '(',nbox,'(E11.4,1X))'
+ open(100,file=trim(files%path_output)//'flx_box.txt',status='replace',action='write',iostat=ierr)
+ do n = 1, ntstate
+ write(100,fmt=rowfmt) flx_box(:,n)
+ end do
+ close(100)
+ endif
+
+ ! prior state vector
+ ! ------------------
+
+ ! number of state variables
+ npvar = ntstate*nbox
+ ndvar = nbox
+ nvar = npvar
+ write(logid,*) 'Number of flux state variables: ',npvar
+ write(logid,*) 'Number of state variables per flux time step: ',ndvar
+ write(logid,*) 'Total number of state variables: ',nvar
+
+ ! allocate state and error vectors
+ call alloc_states(states)
+
+ if ( config%lognormal ) then
+ ! lognormal distribution
+ ! optimize scalar z = ln(x/xb)
+ ! prior state vector z is zero
+ states%px0(1:npvar) = 0.
+ else
+ ! normal distribution
+ ! prior state vector zero since optimize offsets
+ states%px0(1:npvar) = 0.
+ endif
+
+ ! time stamp of flux variables
+ states%xtime(:) = fluxes%time
+
+ if ( config%lognormal ) then
+ ! lognormal
+ ! uncertainty of a scalar of fluxes -> use constant uncertainty
+ ! transform to error of z
+ xerr(:) = log(1.+config%flxerr)
+ else
+ ! normal
+ ! average error over time steps
+ xerr = abs(sum(err_box(:,:),dim=2)/real(ntstate) )
+ endif
+
+ ! assign best guess of state vector
+ if ( (config%prior_bg.eq.1).and.((config%method.eq.'congrad').or.(config%method.eq.'m1qn3')) ) then
+ ! from previous inversion
+ call read_bg(files, config, states)
+ else
+ ! from prior
+ states%px(:) = states%px0(:)
+ endif
+
+ ! testing
+ open(100,file=trim(files%path_output)//'prior.txt',status='replace',action='write',iostat=ierr)
+ do n = 1, nvar
+ write(100,fmt='(E14.8)') states%px(n)
+ end do
+ close(100)
+
+ ! correlation and error covariance matrices
+ ! -----------------------------------------
+
+ corr(:,:) = 0.
+ if ( config%spa_corr ) then
+ write(logid,*) 'Using spatial correlation in B'
+ ! check if file already exists (in case of rerun)
+ inquire(file=trim(files%path_output)//'correl.txt',exist=lexist)
+ if ( lexist ) then
+ write(logid,*) 'Reading spatial correlation from correl.txt'
+ open(100,file=trim(files%path_output)//'correl.txt',action='read',status='old',iostat=ierr)
+ write(rowfmt,'(A,I6,A)') '(',nbox,'(E11.4,1X))'
+ print*, rowfmt
+ do nb = 1, nbox
+ read(100,fmt=rowfmt) corr(nb,:)
+ end do
+ close(100)
+ else
+ call correl(files, config, corr)
+ endif
+ else
+ do n = 1, nbox
+ corr(n,n) = 1.
+ end do
+ endif
+ call error_cov(config, files, states, covar, corr, xerr)
+
+ ! observations
+ ! ------------
+
+ call read_obs(config, files)
+ call assign_obs(config, files, obs)
+
+ ! add random perturbation for MC
+ ! ------------------------------
+
+ if ( config%run_mode.eq.2 ) then
+ write(logid,*) 'Adding random perturbation to state vector and observations...'
+ ! perturb state vector
+ call perturb_state(config, states, covar)
+ ! perturb observations
+ call perturb_obs(config, obs)
+ endif
+
+end subroutine init_aero
+
+
diff --git a/source/init_cini.f90 b/source/init_cini.f90
index 9577c63..0447d07 100644
--- a/source/init_cini.f90
+++ b/source/init_cini.f90
@@ -72,7 +72,7 @@ subroutine init_cini(files, config, obs)
real, dimension(:,:,:,:), allocatable :: concini
real :: bndx, bndy, delx, dely
- integer :: numx, numy, xshift
+ integer :: numx, numy, numz, xshift
integer :: numpoint, release_nr, nr_init
integer :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2) :: releases
@@ -146,12 +146,12 @@ subroutine init_cini(files, config, obs)
if ( i.eq.1 .or. month.ne.prevmonth ) then
write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
else
if(trim(obs%recs(ind(i))).ne.trim(obs%recs(ind(i-1))))then
write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
endif
endif
! find the first release number corresponding to the current observation
diff --git a/source/init_cini_month.f90 b/source/init_cini_month.f90
index f2f56d9..ee6856b 100644
--- a/source/init_cini_month.f90
+++ b/source/init_cini_month.f90
@@ -70,7 +70,7 @@ subroutine init_cini_month(files, config, obs)
real, dimension(nxgrid,nygrid,nzgrid,1):: concini, preconcini, proconcini
real :: bndx, bndy, delx, dely
- integer :: numx, numy, xshift
+ integer :: numx, numy, numz, xshift
integer :: numpoint, release_nr, nr_init
integer :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2) :: releases
@@ -131,12 +131,12 @@ subroutine init_cini_month(files, config, obs)
if ( i.eq.1 .or. month.ne.prevmonth ) then
write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
else
if(trim(obs%recs(ind(i))).ne.trim(obs%recs(ind(i-1))))then
write(logid,*) 'Get header:', i, ' '//trim(path_flexrec)//'header'
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
endif
endif
! find the first release number corresponding to the current observation
diff --git a/source/init_ghg.f90 b/source/init_ghg.f90
index 15ab159..63cb925 100644
--- a/source/init_ghg.f90
+++ b/source/init_ghg.f90
@@ -481,7 +481,7 @@ subroutine init_ghg(files, config, fluxes, obs, states, covar)
! initial concentrations from termination of trajectories
! only do for a fresh run
- if ( ( config%restart.eq.0 ).and.( config%spec.ne.'bca') ) then
+ if ( config%restart.eq.0 ) then
if ( files%init_type.lt.5 ) then
! daily or higher frequency
call init_cini(files, config, obs)
diff --git a/source/initialize.f90 b/source/initialize.f90
index adbb03d..f0b7b1a 100644
--- a/source/initialize.f90
+++ b/source/initialize.f90
@@ -55,7 +55,7 @@ subroutine initialize(files, config)
integer :: ibdate, ibtime, jjjjmmdd, hhmiss
real(kind=8), dimension(maxpoint,2) :: releases
real :: bndx, bndy, delx, dely
- integer :: numx, numy, xshift
+ integer :: numx, numy, numz, xshift
integer :: i, ierr
! initialize log file
@@ -146,11 +146,14 @@ subroutine initialize(files, config)
endif
write(logid,*) 'Reading flexpart header: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
! global domain variables
nxgrid = numx
nygrid = numy
+ ! nzgrid is the number full vertical layers in grid_initial
+ ! read only from header for global grid_time files
+ nzgrid = numz
llx = bndx
lly = bndy
dx = delx
@@ -172,6 +175,11 @@ subroutine initialize(files, config)
write(areldate,'(I8)') config%datei
filename = trim(files%path_flexsat)//areldate//'/header_nest'
inquire(file=trim(filename),exist=lexist)
+ if ( .not.lexist ) then
+ ! test if is FPv11 file
+ filename = trim(files%path_flexsat)//areldate//'/header_nest_grid_time'
+ inquire(file=trim(filename),exist=lexist)
+ endif
else
! read from ground-based output
lexist = .false.
@@ -179,6 +187,11 @@ subroutine initialize(files, config)
do while ( (.not.lexist).and.(i.le.nrec) )
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
i = i + 1
end do
endif ! satellite
@@ -188,7 +201,7 @@ subroutine initialize(files, config)
endif
write(logid,*) 'Reading flexpart header_nest: '//trim(filename)
call read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
! nested domain variables
nnxgrid = numx
diff --git a/source/job_flexinvert.sh b/source/job_flexinvert.sh
index 479ee24..0f92bf6 100755
--- a/source/job_flexinvert.sh
+++ b/source/job_flexinvert.sh
@@ -1,7 +1,7 @@
#!/bin/bash
#---------------------------------------------------
-settings_files='/mypath/settings/SETTINGS_ghg_files'
-settings_config='/mypath/settings/SETTINGS_ghg_config'
+settings_files='../settings/SETTINGS_ghg_files'
+settings_config='../settings/SETTINGS_ghg_config'
./main ${settings_files} ${settings_config}
diff --git a/source/main.f90 b/source/main.f90
index d5e9971..f5ee070 100644
--- a/source/main.f90
+++ b/source/main.f90
@@ -89,10 +89,12 @@ program main
if ( config%spec.eq.'co2' ) then
call init_co2(files, config, fluxes, obs, states, covar)
- else if ((config%spec.eq.'ghg').or.(config%spec.eq.'bca')) then
+ else if ( config%spec.eq.'ghg' ) then
call init_ghg(files, config, fluxes, obs, states, covar)
+ else if ( config%spec.eq.'aero' ) then
+ call init_aero(files, config, fluxes, obs, states, covar)
else
- write(logid,*) 'ERROR: spec must be one of co2 or ghg'
+ write(logid,*) 'ERROR: spec must be one of co2, ghg, aero'
stop
endif
diff --git a/source/makefile b/source/makefile
index 4f2c631..38974c9 100644
--- a/source/makefile
+++ b/source/makefile
@@ -46,6 +46,7 @@ OBJECTS = read_reclist.o \
init_cini_month.o \
init_co2.o \
init_ghg.o \
+ init_aero.o \
calc_conc.o \
simulate.o \
average_fp.o \
diff --git a/source/mod_flexpart.f90 b/source/mod_flexpart.f90
index 08780c0..f4f54fe 100755
--- a/source/mod_flexpart.f90
+++ b/source/mod_flexpart.f90
@@ -49,7 +49,7 @@ module mod_flexpart
! --------------------------------------------------
subroutine read_header(filename, numpoint, ibdate, ibtime, releases, &
- numx, numy, bndx, bndy, delx, dely, xshift)
+ numx, numy, numz, bndx, bndy, delx, dely, xshift)
use mod_var
@@ -59,7 +59,7 @@ module mod_flexpart
integer, intent(in out) :: numpoint
integer, intent(in out) :: ibdate, ibtime
real(kind=8), dimension(maxpoint,2), intent(in out) :: releases
- integer, intent(in out) :: numx, numy, xshift
+ integer, intent(in out) :: numx, numy, numz, xshift
real, intent(in out) :: bndx, bndy, delx, dely
character(len=29) :: flexversion
@@ -75,7 +75,7 @@ module mod_flexpart
real, dimension(:), allocatable :: xpoint, ypoint, zpoint1, zpoint2
real, dimension(:,:), allocatable :: xmass
integer, dimension(10) :: lage
- integer :: n, i, j, ierr
+ integer :: n, i, j, ierr, dummy
inquire ( file=trim(filename),exist=lexist )
if ( .not.lexist ) then
@@ -91,14 +91,14 @@ module mod_flexpart
read(100) ibdate, ibtime, flexversion
read(100) loutstep, loutaver, loutsample
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) nspec, nzgrid
+ read(100) nspec, numz
nspec=nspec/3
do n=1,nspec
- read(100) nzgrid, species(n)
- read(100) nzgrid, species(n)
- read(100) nzgrid, species(n)
+ read(100) dummy, species(n)
+ read(100) dummy, species(n)
+ read(100) numz, species(n)
end do
read(100) numpoint
@@ -382,6 +382,8 @@ module mod_flexpart
integer, dimension(nxgrid*nygrid*nzgrid) :: sparse_dump_i
real, dimension(nxgrid*nygrid*nzgrid) :: sparse_dump_r
+ print*, 'read_init: filename = ',filename
+
gridinit(:,:,:)=0.
inquire ( file=trim(filename),exist=lexist )
@@ -398,6 +400,8 @@ module mod_flexpart
read(100) jjjjmmdd
read(100) hhmiss
+ print*, 'read_init: jjjjmmdd, hhmiss = ',jjjjmmdd, hhmiss
+
fact = 1.
read(100) sp_count_i
read(100) (sparse_dump_i(ix), ix=1, sp_count_i)
diff --git a/source/mod_perturb.f90 b/source/mod_perturb.f90
index d443d49..61e8834 100644
--- a/source/mod_perturb.f90
+++ b/source/mod_perturb.f90
@@ -98,7 +98,8 @@ module mod_perturb
! add perturbation to state vector
if ( config%lognormal ) then
! lognormal
- states%px(:) = real(perturb(:),kind=4)
+ states%px0(:) = real(perturb(:),kind=4)
+ states%px(:) = states%px0(:)
else
! normal
states%px0(:) = states%px0(:) + real(perturb(:),kind=4)
diff --git a/source/mod_save.f90 b/source/mod_save.f90
index aed5666..ec3c977 100644
--- a/source/mod_save.f90
+++ b/source/mod_save.f90
@@ -88,7 +88,7 @@ module mod_save
obs%model(i), obs%cpri(i), obs%cakpri(i), obs%delta(i), obs%err(i)
end do
else
- ! GHG species
+ ! GHG and aerosol species
write(100,*) 'rec yyyymmdd hhmmss juldate conc cini cinipos bkg ghg prior post cpri cakpri diff error'
write(rowfmt,'(A,I1,A,I6,A)') '(A',recname_len,',1X,I8,1X,I6,1X,F14.6,1X,',11,'(F10.4,1X))'
do i = 1, nobs
@@ -412,7 +412,7 @@ module mod_save
! --------------------------------------------------
!> save_ghg
!!
- !! Purpose: Saves the GHG fluxes to the file
+ !! Purpose: Saves the GHG and aerosol fluxes to
!! analysis.nc
!!
! --------------------------------------------------
@@ -448,9 +448,10 @@ module mod_save
integer, dimension(3) :: dimids
integer :: lon_varid, lat_varid, time_varid
integer :: pri_varid, pos_varid, epri_varid, epos_varid, xpos_varid
+ integer :: mcpri_varid, mcpos_varid
integer :: ix, jy, n, ix0, jy0, nb
real(kind=8), dimension(ntstate) :: timeout
- real, dimension(nxregrid,nyregrid,ntstate) :: fpri, fpos, epri, epos, xpos, focn_reg
+ real, dimension(nxregrid,nyregrid,ntstate) :: fpri, fpos, epri, epos, xpos, focn_reg, mcpri, mcpos
real :: area, totpos
! regrid fluxes and errors
@@ -528,6 +529,24 @@ module mod_save
epos(ix,jy,n) = states%pxerr((n-1)*nbox+nbox_xy(ix,jy))/numscale
xpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy))/numscale
endif
+ if ( config%run_mode.eq.2 ) then
+ ! monte carlo ensemble run for posterior uncertainty
+ ! save posterior and perturbed prior state vectors with no redistribution
+ if ( config%lognormal ) then
+ ! lognormal distribution
+ if ( config%nested ) then
+ mcpri(ix,jy,n) = exp(states%px0((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx_nest(ix,jy,n)/numscale
+ mcpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx_nest(ix,jy,n)/numscale
+ else
+ mcpri(ix,jy,n) = exp(states%px0((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx(ix0+ix-1,jy0+jy-1,n)/numscale
+ mcpos(ix,jy,n) = exp(states%px((n-1)*nbox+nbox_xy(ix,jy))) * fluxes%flx(ix0+ix-1,jy0+jy-1,n)/numscale
+ endif
+ else
+ ! normal distribution
+ mcpri(ix,jy,n) = states%px0((n-1)*nbox+nbox_xy(ix,jy))/numscale + fpri(ix,jy,n)
+ mcpos(ix,jy,n) = states%px((n-1)*nbox+nbox_xy(ix,jy))/numscale + fpri(ix,jy,n)
+ endif
+ endif
else
! include ocean fluxes if not optimized
if ( config%nested ) then
@@ -596,6 +615,12 @@ module mod_save
call check( nf90_put_att(ncid, epos_varid, units, flxunit) )
call check( nf90_def_var(ncid, xpos_name, nf90_real, dimids, xpos_varid) )
call check( nf90_put_att(ncid, xpos_varid, units, flxunit) )
+ if ( config%run_mode.eq.2 ) then
+ call check( nf90_def_var(ncid, 'mcpri', nf90_real, dimids, mcpri_varid) )
+ call check( nf90_put_att(ncid, mcpri_varid, units, flxunit) )
+ call check( nf90_def_var(ncid, 'mcpos', nf90_real, dimids, mcpos_varid) )
+ call check( nf90_put_att(ncid, mcpos_varid, units, flxunit) )
+ endif
call check( nf90_enddef(ncid) )
@@ -608,6 +633,10 @@ module mod_save
call check( nf90_put_var(ncid, pos_varid, fpos) )
call check( nf90_put_var(ncid, epos_varid, epos) )
call check( nf90_put_var(ncid, xpos_varid, xpos) )
+ if ( config%run_mode.eq.2 ) then
+ call check( nf90_put_var(ncid, mcpri_varid, mcpri) )
+ call check( nf90_put_var(ncid, mcpos_varid, mcpos) )
+ endif
call check( nf90_close(ncid) )
diff --git a/source/mod_settings.f90 b/source/mod_settings.f90
index 4b54e1c..2320915 100644
--- a/source/mod_settings.f90
+++ b/source/mod_settings.f90
@@ -117,7 +117,7 @@ module mod_settings
integer :: datei ! start date (yyyymmdd)
integer :: datef ! end date (yyyymmdd)
logical :: average_fp ! average footprints to match observation (true or false)
- character(len=3) :: spec ! species ('co2', 'bca' or 'ghg')
+ character(len=3) :: spec ! species ('co2' or 'ghg' or 'aero')
character(len=max_name_len) :: method ! inversion method ('analytic' or 'congrad' or 'm1qn3')
logical :: lognormal ! use log-normal distribution in state space
real :: trunc ! truncation factor for eigenvalues of B
@@ -130,8 +130,11 @@ module mod_settings
integer :: restart ! restart a run that crashed (0=false, 1=true)
real :: molarmass ! molar mass of species (corresponds to flux files, e.g. C=12, CH4=16)
real :: coeff ! conversion coefficient from ppt to ppm or ppbv
+ real :: saterrfactor ! factor by which to scale satellite observation error
+ real :: satlowlim ! exclude satellite observations below this value
logical :: satellite ! use satellite observations (true or false)
logical :: ground ! use ground-based observations (true or false)
+ logical :: selectobs ! select ground-based observations for day/night
real :: coeffsat ! conversion coefficient for satellite data to ppm or ppbv
logical :: nested ! use nested flexpart output (true or false)
real :: w_edge_lon ! lon of western edge of inversion grid
@@ -642,6 +645,7 @@ module mod_settings
config%lognormal = .false.
config%satellite = .false.
config%ground = .false.
+ config%selectobs = .false.
config%trunc = 1.e-3
! open file
@@ -741,9 +745,18 @@ module mod_settings
identifier = "satellite:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%satellite = cl
+ identifier = "saterrfactor:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%saterrfactor = real(cn,kind=4)
+ identifier = "satlowlim:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%satlowlim = real(cn,kind=4)
identifier = "ground:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%ground = cl
+ identifier = "selectobs:"
+ call read_content (line, identifier, cc, cn, cl, match)
+ if ( match ) config%selectobs = cl
identifier = "coeffsat:"
call read_content (line, identifier, cc, cn, cl, match)
if ( match ) config%coeffsat = real(cn,kind=4)
diff --git a/source/read_obs.f90 b/source/read_obs.f90
index d987322..7110f02 100755
--- a/source/read_obs.f90
+++ b/source/read_obs.f90
@@ -52,13 +52,15 @@ subroutine read_obs(config, files)
integer :: ncid, varid, dimid
integer, dimension(:), allocatable :: idate, itime
real, dimension(:), allocatable :: cpri, cakpri, cdryair, vmr, vmr_err
-! real, dimension(:,:), allocatable :: cak
real :: xcpri, xcakpri, scl_cpri, scl_cakpri, scl_cdryair !, cak_av
logical :: lexist
! list observation files
! ----------------------
+ print*, 'readobs: files%path_output = ',trim(files%path_output)
+ print*, 'readobs: files%path_satobs = ',trim(files%path_satobs)
+
if ( config%ground.and.config%satellite ) then
call system('ls '//trim(files%path_obs)//'*'//trim(files%suffix)//&
' '//trim(files%path_satobs)//'*.nc'//' | wc -l > '//trim(files%path_output)//'obsfiles.txt')
@@ -94,7 +96,6 @@ subroutine read_obs(config, files)
write(logid,*) 'ERROR: cannot open obsread.txt'
stop
endif
-! write(200,*) 'rec yyyymmdd hhmmss juldate avetime conc error cpri cakpri cdryair cakav'
write(200,*) 'rec yyyymmdd hhmmss juldate avetime conc error cpri cakpri cdryair'
cnt = 0
@@ -102,7 +103,9 @@ subroutine read_obs(config, files)
! check if ground-based or satellite observation
string = filelist(nf)
+ print*, 'readobs: string = ',string
call split(string,".",before,sep)
+ print*, 'readobs: before, string = ',before, string
if ( string.eq.'nc' ) then
@@ -111,7 +114,6 @@ subroutine read_obs(config, files)
! check if data is in inversion time interval
read(before(19:26),*) yyyymmdd
- print*, 'readobs: string before = ',before(19:26)
jdate = juldate(yyyymmdd, 0)
if ( jdate.lt.juldatei.or.jdate.gt.juldatef ) cycle
@@ -142,52 +144,34 @@ subroutine read_obs(config, files)
call check( nf90_get_var(ncid,varid,vmr_err) )
call check( nf90_close(ncid) )
- ! read ncdf file of releases (for AK info)
- ! this is not needed
-! write(adate,'(I8)') yyyymmdd
-! print*, trim(filelist(nf))
-! filerelease=str_replace(trim(filelist(nf)),'retrieval','releases')
-! inquire(file=trim(files%path_satrel)//adate//'/options/'//trim(filerelease),exist=lexist)
-! if (.not.lexist) then
-! write(logid,*) 'ERROR: file not found ',trim(files%path_satrel)//adate//'/options/'//trim(filerelease)
-! stop
-! endif
-! call check( nf90_open(trim(files%path_satrel)//adate//'/options/'//trim(filerelease),nf90_NOWRITE,ncid) )
-! call check( nf90_inq_dimid(ncid,'retrieval',dimid) )
-! call check( nf90_inquire_dimension(ncid,dimid,len=nrelease) )
-! call check( nf90_inq_dimid(ncid,'nlayer',dimid) )
-! call check( nf90_inquire_dimension(ncid,dimid,len=nlayer) )
-! if (nrelease.ne.nretr) then
-! write(logid,*) 'ERROR: number of releases does not much number of retrievals'
-! stop
-! endif
-! allocate( cak(nretr,nlayer), stat=ierr )
-! call check( nf90_inq_varid(ncid,'cak',varid) )
-! call check( nf90_get_var(ncid,varid,cak) )
-! call check( nf90_close(ncid) )
-
! loop over retrievals
avetime = 0d0
write(numfmt,fmt='(A,I1,A,I1,A)') '(I',recname_len,'.',recname_len,')'
-! write(rowfmt,fmt='(A,I1,A)') '(A',recname_len,',1X,I8,1X,I6,1X,F14.6,1X,F10.6,1X,F10.4,1X,F10.4,1X,F10.5,1X,F10.5,1X,F10.2,1X,F6.4)'
write(rowfmt,fmt='(A,I1,A)') '(A',recname_len,',1X,I8,1X,I6,1X,F14.6,1X,F10.6,1X,F10.4,1X,F10.4,1X,F10.5,1X,F10.5,1X,F10.2)'
do nr = 1, nretr
- cnt = cnt + 1
write(recs,fmt=numfmt) nr
+ ! exclude vmr values that are anomolously low
+ if ((vmr(nr).lt.config%satlowlim).and.(config%satlowlim.gt.0)) then
+ write(logid,*) 'WARNING: excluding low satellite value for = ',yyyymmdd, nr
+ cycle
+ endif
+ cnt = cnt + 1
! convert from input unit (e.g. mol/m2) to mole fraction
xcpri = cpri(nr)*config%coeffsat/cdryair(nr)
xcakpri = cakpri(nr)*config%coeffsat/cdryair(nr)
-! cak_av = sum(cak(nr,:))/real(nlayer)
print*, 'read_obs: xcpri, xcakpri = ',xcpri, xcakpri
print*, 'read_obs: cpri, cdryair = ',cpri(nr), cdryair(nr)
print*, 'read_obs: coeffsat = ',config%coeffsat
! julian time
jdate = juldate(idate(nr), itime(nr))
+ ! scale retrieval error by factor
+ if (config%saterrfactor.gt.0) then
+ vmr_err(nr) = vmr_err(nr)*config%saterrfactor
+ endif
write(200,fmt=rowfmt) recs, idate(nr), itime(nr), jdate, avetime, vmr(nr), vmr_err(nr), &
xcpri, xcakpri, cdryair(nr) !, cak_av
end do
-! deallocate( idate, itime, cpri, cakpri, cdryair, vmr, vmr_err, cak )
deallocate( idate, itime, cpri, cakpri, cdryair, vmr, vmr_err )
else
@@ -249,30 +233,30 @@ subroutine read_obs(config, files)
if ( ierr.ne.0 ) exit read_loop
if ( jdate.ge.(juldatef+1d0) ) exit read_loop
if ( jdate.lt.juldatei ) cycle read_loop
- ! select day/night for low/high alt sites
- hh = hhmiss/10000
- hl = hh + int(lon*24./360.)
- if ( hl.lt.0 ) then
- hl = hl + 24
- else if ( hl.ge.24 ) then
- hl = hl - 24
- endif
- if ( alt.le.1000. ) then
-! if ( (hl.lt.11).or.(hl.gt.15) ) cycle read_loop
- else
-! if ( (hl.gt.3).and.(hl.lt.23) ) cycle read_loop
+ if ( config%selectobs ) then
+ ! select day/night for low/high alt sites
+ hh = hhmiss/10000
+ hl = hh + int(lon*24./360.)
+ if ( hl.lt.0 ) then
+ hl = hl + 24
+ else if ( hl.ge.24 ) then
+ hl = hl - 24
+ endif
+ if ( alt.le.1000. ) then
+ if ( (hl.lt.11).or.(hl.gt.15) ) cycle read_loop
+ else
+ if ( (hl.gt.3).and.(hl.lt.23) ) cycle read_loop
+ endif
endif
if ( conc.le.-999. ) cycle read_loop
measerr = max(config%measerr, err)
scl_cpri = 0.
scl_cakpri = 0.
scl_cdryair = 0.
-! cak_av = 1.
cnt = cnt + 1
! write to output file
-! write(rowfmt,'(A,I1,A)') '(A',recname_len,',1X,I8,1X,I6,1X,F14.6,1X,F10.6,1X,F10.4,1X,F10.4,F10.5,1X,F10.5,1X,F10.2,1X,F6.4)'
write(rowfmt,'(A,I1,A)') '(A',recname_len,',1X,I8,1X,I6,1X,F14.6,1X,F10.6,1X,F10.4,1X,F10.4,F10.5,1X,F10.5,1X,F10.2)'
- write(200,fmt=rowfmt) recs, yyyymmdd, hhmiss, jdate, avetime, conc, measerr, scl_cpri, scl_cakpri, scl_cdryair !, cak_av
+ write(200,fmt=rowfmt) recs, yyyymmdd, hhmiss, jdate, avetime, conc, measerr, scl_cpri, scl_cakpri, scl_cdryair
end do read_loop
! close input file
diff --git a/source/sbatch_flexinvert.sh b/source/sbatch_flexinvert.sh
index b82c3b9..ce2779f 100755
--- a/source/sbatch_flexinvert.sh
+++ b/source/sbatch_flexinvert.sh
@@ -2,8 +2,8 @@
#---------------------------------------------------
partition=nilu
jobname=ghg
-settings_files='/mypath/settings/SETTINGS_ghg_files'
-settings_config='/mypath/settings/SETTINGS_ghg_config'
+settings_files='../settings/SETTINGS_ghg_files'
+settings_config='../settings/SETTINGS_ghg_config'
#---------------------------------------------------
cat <<EOF > run_job.sh
diff --git a/source/simulate.f90 b/source/simulate.f90
index 999c532..23bebb8 100644
--- a/source/simulate.f90
+++ b/source/simulate.f90
@@ -173,7 +173,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
ftime(n) = juldate(int(dates(n)/1d6),int(dates(n)-floor(dates(n)/1d6)*1d6))
end do
factor(:,:,:) = 1.
- if ( config%spec.ne.'bca' ) then
+ if ( config%spec.ne.'aero' ) then
filefactor = trim(path_flexrec)//'factor_drygrid'
inquire ( file=trim(filefactor),exist=lexist )
if ( lexist ) then
@@ -239,11 +239,11 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
end do
! convert s.m3/kg to s.m2/kg
grid = grid/outheight(1)
- ! convert from equivalent ppt to observation units (e.g. ppmv)
- ! sec not do it for bca
- if ( config%spec.eq.'bca' ) then
+ if ( config%spec.eq.'aero' ) then
+ ! convert to match concentration units
grid = grid*config%coeff
else
+ ! convert from equivalent ppt to observation units (e.g. ppmv)
grid = grid*config%coeff*mmair/config%molarmass
endif
@@ -262,7 +262,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
if ( allocated(factor_nest) ) deallocate(factor_nest)
allocate ( factor_nest(nnxgrid,nnygrid,nread) )
factor(:,:,:) = 1.
- if ( config%spec.ne.'bca' ) then
+ if ( config%spec.ne.'aero' ) then
filefactor = trim(path_flexrec)//'factor_drygrid_nest'
inquire ( file=trim(filefactor),exist=lexist )
if ( lexist ) then
@@ -316,8 +316,7 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
! convert s.m3/kg to s.m2/kg
gridnest = gridnest/outheight(1)
! convert from equivalent ppt to observation units (e.g. ppmv)
- ! sec not do it for bca
- if ( config%spec.eq.'bca' ) then
+ if ( config%spec.eq.'aero' ) then
gridnest = gridnest*config%coeff
else
gridnest = gridnest*config%coeff*mmair/config%molarmass
@@ -592,11 +591,11 @@ subroutine simulate(iter, files, config, fluxes, obs, states, grad_o, cost_o)
endif
endif
- ! if background missing then observation should have no influence!
-! if ( sum(obs%cini(i,:)).eq.0. ) then
-! write(logid,*) 'WARNING: background zero for observation ',i
-! obs%delta(i) = 0.
-! endif
+ ! check if background missing (not for aerosols)
+ if ( config%spec.ne.'aero'.and.sum(obs%cini(i,:)).eq.0. ) then
+ write(logid,*) 'ERROR: background zero for observation ',i
+ stop
+ endif
! calculate gradient in observation space
! ---------------------------------------
--
GitLab