Skip to content
Snippets Groups Projects
Commit bfe92558 authored by Lucie Bakels's avatar Lucie Bakels
Browse files

Merge branch 'fix_GFS' into development

Former-commit-id: b938e8be
parents 7a8200bf 7e51d823
No related branches found
No related tags found
No related merge requests found
Showing
with 1104 additions and 442 deletions
......@@ -24,11 +24,30 @@ alma8-build:
- ./src/FLEXPART_ETA
expire_in: never
rocky9-build:
image: harbor.wolke.img.univie.ac.at/flexpart/rockylinux9:latest
stage: build
script:
- export FC=gfortran
- export LIBRARY_PATH=/usr/lib64:/usr/local/lib64
- export CPATH=/usr/include:/usr/local/include:/usr/lib64/gfortran/modules
- make -j -C src -f makefile_gfortran eta=no
- make -j -C src -f makefile_gfortran
artifacts:
when: on_success
paths:
- ./src/FLEXPART
- ./src/FLEXPART_ETA
expire_in: never
options-test:
image: harbor.wolke.img.univie.ac.at/flexpart/almalinux8:latest
# image: harbor.wolke.img.univie.ac.at/flexpart/almalinux8:latest
image: harbor.wolke.img.univie.ac.at/flexpart/rockylinux9:latest
stage: test
needs:
- alma8-build
- rocky9-build
script:
- ulimit -s unlimited
......@@ -42,12 +61,30 @@ options-test:
- ./tests/output_bkw_eta
expire_in: 5 mins
nests-test:
# image: harbor.wolke.img.univie.ac.at/flexpart/almalinux8:latest
image: harbor.wolke.img.univie.ac.at/flexpart/rockylinux9:latest
stage: test
when: manual
needs:
- alma8-build
- rocky9-build
script:
- ulimit -s unlimited
- bash ./tests/run_nests_test.sh
artifacts:
when: on_success
paths:
openmp-simulation:
image: harbor.wolke.img.univie.ac.at/flexpart/almalinux8:latest
# image: harbor.wolke.img.univie.ac.at/flexpart/almalinux8:latest
image: harbor.wolke.img.univie.ac.at/flexpart/rockylinux9:latest
stage: test
needs:
- alma8-build
- rocky9-build
script:
- ulimit -s unlimited
......@@ -63,11 +100,13 @@ openmp-simulation:
expire_in: 5 mins
etex-simulation:
image: harbor.wolke.img.univie.ac.at/flexpart/almalinux8:latest
# image: harbor.wolke.img.univie.ac.at/flexpart/almalinux8:latest
image: harbor.wolke.img.univie.ac.at/flexpart/rockylinux9:latest
stage: test
when: manual
needs:
- alma8-build
- rocky9-build
script:
- ulimit -s unlimited
......@@ -133,3 +172,10 @@ etex-test:
paths:
- ./tests/etex_test.txt
expire_in: never
documentation:
image: harbor.wolke.img.univie.ac.at/podman/mkdocs-computer:latest
stage: build
script:
- cd ./documentation && mkdocs build -c --verbose
- sshpass -p "$WOLKE_PASSWORD" rsync -autv --delete -e "ssh -o StrictHostKeyChecking=no" /tmp/cr-site/* "$WOLKE_USER@wolke.img.univie.ac.at:/var/www/html/documentation/flexpart"
\ No newline at end of file
......@@ -2,26 +2,31 @@
# Dockerfile for CI and Flexpart container images
# Build Examples:
# - podman build -t harbor.wolke.img.univie.ac.at/flexpart/almalinux8 -f Dockerfile
# - podman build -t harbor.wolke.img.univie.ac.at/flexpart/rockylinux9 -f Dockerfile
#
FROM almalinux:8-minimal
# FROM rockylinux:8-minimal
FROM rockylinux:9-minimal
#
# Build development image (with/without jasper)
# jasper was used in FP 10.4 (eccodes, emoslib)
#
# 8: --enablerepo=powertools
# since 9, there is 2.31 eccodes in epel
RUN microdnf install -y epel-release && \
microdnf install -y --enablerepo=powertools make netcdf-fortran-devel.x86_64 netcdf.x86_64 cmake tar gcc-c++ perl && \
microdnf clean all -y
microdnf install -y --enablerepo=crb make netcdf-fortran-devel.x86_64 netcdf.x86_64 eccodes eccodes-devel cmake tar gcc-c++ perl git && \
microdnf clean all -y && \
rm -rf /var/cache/yum
#
# Download ECCODES Version
# note 2.30.0 has an issue!!!
#
RUN curl https://confluence.ecmwf.int/download/attachments/45757960/eccodes-2.31.0-Source.tar.gz | tar xz
RUN mkdir build && \
cd build && \
cmake -DENABLE_ECCODES_OMP_THREADS=ON ../eccodes-*/ && \
make -j8 && \
make install
# RUN curl https://confluence.ecmwf.int/download/attachments/45757960/eccodes-2.31.0-Source.tar.gz | tar xz
# RUN mkdir build && \
# cd build && \
# cmake -DENABLE_ECCODES_OMP_THREADS=ON ../eccodes-*/ && \
# make -j8 && \
# make install
#
# set environment variables
#
......
......@@ -4,12 +4,13 @@
# - COMMIT = ... (use git log --pretty=format:'%h %D (%s %at)' -n 1)
# with this option Flexpart is build inside the container
# Build Examples:
# - podman build -t harbor.wolke.img.univie.ac.at/flexpart/flexpartv11:$(git rev-parse --abbrev-ref HEAD) -f Dockerfile_flexpart --build-arg COMMIT=$(git log --pretty=format:'%h %D %ad' -n 1)
# - podman build -t harbor.wolke.img.univie.ac.at/flexpart/flexpartv11:$(git rev-parse --abbrev-ref HEAD) -f Dockerfile_flexpart --build-arg COMMIT="$(git log --pretty=format:'%h %D %ad' -n 1)"
#
#
# Build image with Flexpart inside
#
FROM harbor.wolke.img.univie.ac.at/flexpart/almalinux8:latest
#FROM harbor.wolke.img.univie.ac.at/flexpart/almalinux8:latest
FROM harbor.wolke.img.univie.ac.at/flexpart/rockylinux9:latest
ARG COMMIT=0
ENV COMMIT=$COMMIT
COPY ./src /src
......@@ -17,7 +18,16 @@ COPY ./tests/default_options /options
COPY ./tests/default_winds /inputs
COPY ./entrypoint.sh /entrypoint.sh
WORKDIR /src
RUN make -f makefile_gfortran && \
mkdir -p /output && \
echo -e "/options/\n/output/\n/inputs/\n/inputs/AVAILABLE" > /pathnames
#
# compile using a standardized
# need to fix march to x86_64_v3 , remove mtune
# here generic means core_avx2
RUN make -f makefile_gfortran eta=no arch=generic \
&& make -f makefile_gfortran arch=generic \
&& mkdir -p /output \
&& echo -e "/options/\n/output/\n/inputs/\n/inputs/AVAILABLE" > /pathnames
# In the Dockerfile, the ENTRYPOINT command defines the executable,
# while CMD sets the default parameter.
# so here our entrypoint run script and default FLEXPART_ETA or FLEXPART
ENTRYPOINT ["/entrypoint.sh"]
CMD ["FLEXPART_ETA"]
......@@ -63,11 +63,354 @@ Now you are almost ready to run.
In the winds are available in flex_ecmwf/work it should suffice to execute
`./src/FLEXPART` in the main directory
### Container version
since version 11, FLEXPART is also available as a container. There is a Dockerfile and some instructions on how to build the container from scratch or you can download the container from our [registry](https://harbor.wolke.img.univie.ac.at) using for example podman/docker or singularity/apptainer:
Specifications:
- compiled using `march=core-avx2`, compatible with CPUs above AVX2 (since Haswell or Zen)
- ecCodes
```sh
# registry download, replace podman to docker or vice versa
podman pull harbor.wolke.img.univie.ac.at/flexpart/flexpartv11:master
# using singularity/apptainer for running FLEXPART
# there might be some warnings about EPERM (can be ignored)
apptainer pull flexpart.sif docker://harbor.wolke.img.univie.ac.at/flexpart/flexpartv11:master
```
running it using podman/docker or singularity/apptainer:
```sh
# simple run the container with default settings
podman run harbor.wolke.img.univie.ac.at/flexpart/flexpartv11:master
# running the container requires a writable output directory
# mounting the local directory to /output inside the container
apptainer run -B .:/output flexpart.sif
```
<details>
<summary>podman/docker flexpart run log</summary>
<pre><code class="shell">
Welcome, running FLEXPART
Using defaults (/pathnames)
/options/
/output/
/inputs/
/inputs/AVAILABLE
Mount volumes to change inputs
Git: c5bdd94 HEAD -> master, origin/master, origin/HEAD Tue Nov 21 16:15:27 2023 +0100
EXECUTING FLEXPART
trying to execute: /src/FLEXPART_ETA
Executing: /src/FLEXPART_ETA
Welcome to FLEXPART Version 11
Git: undefined
FLEXPART is free software released under the GNU General Public License.
FLEXPART is running with ETA coordinates.
----------------
INFORMATION: SUBGRIDSCALE TERRAIN EFFECT IS
NOT PARAMETERIZED DURING THIS SIMULATION.
----------------
*********** WARNING **********************************
* FLEXPART running in parallel mode *
* Number of uncertainty classes in *
* set to number of threads: 1 *
* All other computations are done with 8 threads. *
*******************************************************
FLEXPART WARNING: TIME DIFFERENCE BETWEEN TWO
WIND FIELDS IS BIG. THIS MAY CAUSE A DEGRADATION
OF SIMULATION QUALITY.
ECMWF metdata detected
NXSHIFT is set to 0
grid dim: 181 91 92 92 91 92
Vertical levels in ECMWF data: 92 92
Mother domain:
Longitude range: -178.00000 to 182.00000 Grid distance: 2.00000
Latitude range : -90.00000 to 90.00000 Grid distance: 2.00000
Number of receptors: 2
Releasepoints : 1
reading SPECIES 24
Particle shape SPHERE for particle 24
SPECIES: 24 AIRTRACER (GAS)
Wet removal for gases is turned: OFF
Dry deposition for gases is turned: OFF
Below-cloud scavenging: OFF
In-cloud scavenging: OFF
Particles released (numpartmax): 10000
Total mass released: 1.0000000E+00
Allocating fields for global output (x,y): 85 65
Concentrations are calculated using kernel
WARNING: turbulence switched off.
Simulated 0.0 hours ( 0 s), 0 particles
Time: 0 seconds. Total spawned: 0 alive: 0 terminated: 0
Allocating 10000 particles 0 0 0
Finished allocation
Time: 900 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 1800 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 2700 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 3600 seconds. Total spawned: 10000 alive: 10000 terminated: 0
3600 Seconds simulated: 10000 Particles: Uncertainty: 0.000 0.000 0.000
Time: 4500 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 5400 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 6300 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 7200 seconds. Total spawned: 10000 alive: 10000 terminated: 0
7200 Seconds simulated: 10000 Particles: Uncertainty: 0.000 0.000 0.000
Time: 8100 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 9000 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 9900 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 10800 seconds. Total spawned: 10000 alive: 10000 terminated: 0
10800 Seconds simulated: 10000 Particles: Uncertainty: 0.000 0.000 0.000
Read wind fields: 0.609375000 seconds
Timemanager: 1.17187500 seconds,first timestep: 1.09375000 seconds
Write particle files: 4.68750000E-02 seconds
Total running time: 1.81250000 seconds
tps,io,tot: 1.95312500E-02 0.131249994 1.81250000
CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!
FINISHED
</code></pre>
</details>
<details>
<summary>singularity/apptainer flexpart run log</summary>
<pre>
<code class="language-shell">
INFO: gocryptfs not found, will not be able to use gocryptfs
Welcome, running FLEXPART
Using defaults (/pathnames)
/options/
/output/
/inputs/
/inputs/AVAILABLE
Mount volumes to change inputs
Git: c5bdd94 HEAD -> master, origin/master, origin/HEAD Tue Nov 21 16:15:27 2023 +0100
EXECUTING FLEXPART
trying to execute: /src/FLEXPART_ETA
Executing: /src/FLEXPART_ETA
Welcome to FLEXPART Version 11
Git: undefined
FLEXPART is free software released under the GNU General Public License.
FLEXPART is running with ETA coordinates.
----------------
INFORMATION: SUBGRIDSCALE TERRAIN EFFECT IS
NOT PARAMETERIZED DURING THIS SIMULATION.
----------------
*********** WARNING **********************************
* FLEXPART running in parallel mode *
* Number of uncertainty classes in *
* set to number of threads: 1 *
* All other computations are done with 8 threads. *
*******************************************************
FLEXPART WARNING: TIME DIFFERENCE BETWEEN TWO
WIND FIELDS IS BIG. THIS MAY CAUSE A DEGRADATION
OF SIMULATION QUALITY.
ECMWF metdata detected
NXSHIFT is set to 0
grid dim: 181 91 92 92 91 92
Vertical levels in ECMWF data: 92 92
Mother domain:
Longitude range: -178.00000 to 182.00000 Grid distance: 2.00000
Latitude range : -90.00000 to 90.00000 Grid distance: 2.00000
Number of receptors: 2
Releasepoints : 1
reading SPECIES 24
Particle shape SPHERE for particle 24
SPECIES: 24 AIRTRACER (GAS)
Wet removal for gases is turned: OFF
Dry deposition for gases is turned: OFF
Below-cloud scavenging: OFF
In-cloud scavenging: OFF
Particles released (numpartmax): 10000
Total mass released: 1.0000000E+00
Allocating fields for global output (x,y): 85 65
Concentrations are calculated using kernel
WARNING: turbulence switched off.
Simulated 0.0 hours ( 0 s), 0 particles
Time: 0 seconds. Total spawned: 0 alive: 0 terminated: 0
Allocating 10000 particles 0 0 0
Finished allocation
Time: 900 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 1800 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 2700 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 3600 seconds. Total spawned: 10000 alive: 10000 terminated: 0
3600 Seconds simulated: 10000 Particles: Uncertainty: 0.000 0.000 0.000
Time: 4500 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 5400 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 6300 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 7200 seconds. Total spawned: 10000 alive: 10000 terminated: 0
7200 Seconds simulated: 10000 Particles: Uncertainty: 0.000 0.000 0.000
Time: 8100 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 9000 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 9900 seconds. Total spawned: 10000 alive: 10000 terminated: 0
Time: 10800 seconds. Total spawned: 10000 alive: 10000 terminated: 0
10800 Seconds simulated: 10000 Particles: Uncertainty: 0.000 0.000 0.000
Read wind fields: 0.703125000 seconds
Timemanager: 1.25000000 seconds,first timestep: 1.15625000 seconds
Write particle files: 1.56250000E-02 seconds
Total running time: 2.00000000 seconds
tps,io,tot: 2.34375000E-02 0.143749997 2.00000000
CONGRATULATIONS: YOU HAVE SUCCESSFULLY COMPLETED A FLEXPART MODEL RUN!
FINISHED
</code></pre>
</details>
<details>
<summary>ncdump of default run output</summary>
<pre><code class="language-shell">
# check the output
ncdump -h grid_conc_20090101000000.nc
netcdf grid_conc_20090101000000 {
dimensions:
time = UNLIMITED ; // (3 currently)
longitude = 85 ;
latitude = 65 ;
height = 4 ;
numspec = 1 ;
pointspec = 1 ;
nageclass = 1 ;
nchar = 45 ;
ncharrec = 16 ;
numpoint = 1 ;
receptor = UNLIMITED ; // (2 currently)
variables:
int time(time) ;
time:units = "seconds since 2009-01-01 00:00" ;
time:calendar = "proleptic_gregorian" ;
float longitude(longitude) ;
longitude:long_name = "longitude in degree east" ;
longitude:axis = "Lon" ;
longitude:units = "degrees_east" ;
longitude:standard_name = "grid_longitude" ;
longitude:description = "grid cell centers" ;
float latitude(latitude) ;
latitude:long_name = "latitude in degree north" ;
latitude:axis = "Lat" ;
latitude:units = "degrees_north" ;
latitude:standard_name = "grid_latitude" ;
latitude:description = "grid cell centers" ;
float height(height) ;
height:units = "meters" ;
height:positive = "up" ;
height:standard_name = "height" ;
height:long_name = "height above ground" ;
char RELCOM(numpoint, nchar) ;
RELCOM:long_name = "release point name" ;
float RELLNG1(numpoint) ;
RELLNG1:units = "degrees_east" ;
RELLNG1:long_name = "release longitude lower left corner" ;
float RELLNG2(numpoint) ;
RELLNG2:units = "degrees_east" ;
RELLNG2:long_name = "release longitude upper right corner" ;
float RELLAT1(numpoint) ;
RELLAT1:units = "degrees_north" ;
RELLAT1:long_name = "release latitude lower left corner" ;
float RELLAT2(numpoint) ;
RELLAT2:units = "degrees_north" ;
RELLAT2:long_name = "release latitude upper right corner" ;
float RELZZ1(numpoint) ;
RELZZ1:units = "meters" ;
RELZZ1:long_name = "release height bottom" ;
float RELZZ2(numpoint) ;
RELZZ2:units = "meters" ;
RELZZ2:long_name = "release height top" ;
int RELKINDZ(numpoint) ;
RELKINDZ:long_name = "release kind" ;
int RELSTART(numpoint) ;
RELSTART:units = "seconds" ;
RELSTART:long_name = "release start relative to simulation start" ;
int RELEND(numpoint) ;
RELEND:units = "seconds" ;
RELEND:long_name = "release end relative to simulation start" ;
int RELPART(numpoint) ;
RELPART:long_name = "number of release particles" ;
float RELXMASS(numspec, numpoint) ;
RELXMASS:long_name = "total release particle mass" ;
int LAGE(nageclass) ;
LAGE:units = "seconds" ;
LAGE:long_name = "age class" ;
int ORO(latitude, longitude) ;
ORO:standard_name = "surface altitude" ;
ORO:long_name = "outgrid surface altitude" ;
ORO:units = "m" ;
char receptor(receptor, ncharrec) ;
receptor:long_name = "receptor name" ;
float spec001_mr(nageclass, pointspec, time, height, latitude, longitude) ;
spec001_mr:units = "ng m-3" ;
spec001_mr:long_name = "AIRTRACER" ;
spec001_mr:decay = -0.07001485f ;
spec001_mr:weightmolar = 29.f ;
spec001_mr:ohcconst = -9.e-10f ;
spec001_mr:ohdconst = -9.9f ;
spec001_mr:vsetaver = 0.f ;
float receptor_conc001(receptor, time) ;
receptor_conc001:units = "ng m-3" ;
receptor_conc001:_FillValue = -1.f ;
receptor_conc001:positive = "up" ;
receptor_conc001:standard_name = "receptor_conc" ;
receptor_conc001:long_name = "receptor_concentration" ;
// global attributes:
:Conventions = "CF-1.6" ;
:title = "FLEXPART model output" ;
:git = "undefined" ;
:source = "Version 11 model output" ;
:history = "2023-11-21 16:22 +0100 created by mblaschek on NB513" ;
:references = "Stohl et al., Atmos. Chem. Phys., 2005, doi:10.5194/acp-5-2461-200" ;
:outlon0 = -25.f ;
:outlat0 = 10.f ;
:dxout = 1.f ;
:dyout = 1.f ;
:ldirect = 1 ;
:ibdate = "20090101" ;
:ibtime = "000000" ;
:iedate = "20090101" ;
:ietime = "030000" ;
:loutstep = 3600 ;
:loutaver = 3600 ;
:loutsample = 900 ;
:loutrestart = -1 ;
:lsynctime = 900 ;
:ctl = -0.2f ;
:ifine = 1 ;
:iout = 1 ;
:ipout = 0 ;
:lsubgrid = 0 ;
:lconvection = 0 ;
:lagespectra = 0 ;
:ipin = 0 ;
:ioutputforeachrelease = 0 ;
:iflux = 0 ;
:mdomainfill = 0 ;
:ind_source = 1 ;
:ind_receptor = 1 ;
:mquasilag = 0 ;
:nested_output = 0 ;
:sfc_only = 0 ;
:linit_cond = 0 ;
}
</code>
</pre>
</details>
### Contribution guidelines
* The version contributed should compile on a reference version of the system and compiler.
- `FLEXPART 10.4` used as reference gfortran 5.4 on Ubuntu 16.04
- `FLEXPART 11` uses as reference gfortran 8.5.0 on AlmaLinux 8
- `FLEXPART 11` uses as reference gfortran 8.5.0 on AlmaLinux 8/RockyLinux 8 or gfortran 11.4.1 on RockyLinux 9
* Code contribution including new features and bug fixes should be complemented with appropriate tests
An essential test consists of a set of input files and directories that allow FLEXPART to run.
......
# Installation
# Building
## Download FLEXPART
There are two options to download _FLEXPART_:
......@@ -21,7 +21,7 @@ _FLEXPART_ 11 is written in Fortran 2018. The following compilers can be used to
- GNU Fortran compiler version 8+ (`gfortran`)
- Intel Fortran compiler (`ifort`)
For running _FLEXPART_ in parallel mode, a compiler supporting [OpenMP](https://www.openmp.org/) is required.
For running _FLEXPART_ in parallel mode, a compiler supporting [OpenMP](https://www.openmp.org/) is required. In addition, libraries (in particular hdf5) should be compiled threadsafe.
## Libraries
_FLEXPART_ uses the following libraries:
......
This diff is collapsed.
# Welcome to the FLEXPART 11 documentation
Information on how to download and install FLEXPART can be found [here](installation.md).
How to set up a simulation, with an explanation of all input files can be found [here](running.md), and some examples can be found [here](examples.md).
Information on how to download and install FLEXPART can be found [here](building.md).
How to set up a simulation, with an explanation of all input files can be found [here](configuration.md#config), and some examples can be found [here](examples.md).
A list of all possible output options and files can be found [here](output.md).
An overview of all processes relation to the direct transport of particles can be found [here](transport.md), and internal particle processes [here](evolution.md).
......
This diff is collapsed.
# Trouble Shooting
Here we provide a list of common problems and their solutions.
<span style="color:seagreen;">
If you have a problem that is not listed or clearly explained by an error message, please create a ticket on our gitlab page.
</span>
#### **My application crashes with a segmentation fault shortly after its launch**
This could be due to having compiled with OpenMP, but not setting the following `ulimit` before launching your application, resulting in too little memory per core:
~~~
ulimit -s unlimited
~~~
#### **My application is unexpectedly slow**
There could be many reasons for this to happen, here are some tips to make your application faster:
- Make sure you compiled using appropriate optimisation flags (see [Optimisation](building.md#paths)).
- Check if you need the resolution of output grid, number of particles, and timestep you are currently using and reduce these where possible. For an explanation of all options, see [Configuration](configuration.md).
- Use more OpenMP threads when running your application.
- Make sure to set the following options in your command line (or submit script if you use one) before launching your application:
~~~
export OMP_PLACES=cores
export OMP_PROC_BIND=true
~~~
......@@ -30,10 +30,12 @@ markdown_extensions:
# - pymdownx.tilde
nav:
- 'index.md'
- 'installation.md'
- 'building.md'
- 'configuration.md'
- 'running.md'
- 'output.md'
- 'transport.md'
- 'evolution.md'
- 'examples.md'
- 'troubleshooting.md'
#!/bin/bash
# By MB
# run FLEXPART, but provide a bit of information
echo "Welcome, running FLEXPART "
echo "Using defaults:"
echo "Using defaults (/pathnames)"
cat /pathnames
echo "Mount volumes to change inputs"
echo "Git: $COMMIT"
echo "EXECUTING FLEXPART"
/src/FLEXPART /pathnames
if [ $# -eq 1 ]; then
echo "trying to execute: /src/$1"
if [ -e /src/"$1" ]; then
echo "Executing: /src/$1"
/src/"$1" /pathnames
else
echo "Falling back to default, executing: /src/FLEXPART_ETA"
/src/FLEXPART_ETA /pathnames
fi
else
echo "Executing: /src/FLEXPART_ETA"
/src/FLEXPART_ETA /pathnames
fi
echo "FINISHED"
......@@ -41,4 +41,5 @@
MAXTHREADGRID= 1, ! Set maximum number of threads for doing grid computations. Recommended to set this no higher than 16. High numbers create more overhead and a larger memory footprint, 1=no parallelisation on grid.
MAXFILESIZE= 10000, ! Maximum output of each partoutput NetCDF-4 file in Mb before a new one is created
LOGVERTINTERP= 0, ! Flag to set all vertical interpolation to logarithmic instead of linear
BCSCHEME= 1, ! Below-cloud scheme; [1] Laakso & Kyro [2] Wang et al [3] modified Wang et al
/
......@@ -40,6 +40,7 @@ program flexpart
implicit none
integer :: i
real :: s_timemanager
character(len=256) :: &
inline_options ! pathfile, flexversion, arg2
......@@ -50,10 +51,10 @@ program flexpart
CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
s_total = (count_clock - count_clock0)/real(count_rate)
! FLEXPART version string
flexversion_major = '11' ! Major version number, also used for species file names
flexversion='Version '//trim(flexversion_major)
flexversion='Version '//trim(flexversion_major)//'.0 (2023-07-11)'
verbosity=0
call update_gitversion(gitversion_tmp)
......@@ -119,6 +120,21 @@ program flexpart
CALL SYSTEM_CLOCK(count_clock, count_rate, count_max)
s_total = (count_clock - count_clock0)/real(count_rate) - s_total
if (verbosity.gt.0) then
! NIK 16.02.2005
do i=1,nspec
if (icnt_incld(i).gt.0) then
write(*,*) '**********************************************'
write(*,*) 'Scavenging statistics for species ', species(i), ':'
write(*,*) 'Total number of occurences of below-cloud scavenging', &
& icnt_belowcld(i)
write(*,*) 'Total number of occurences of in-cloud scavenging', &
& icnt_incld(i)
write(*,*) '**********************************************'
endif
end do
endif
write(*,*) 'Read wind fields: ', s_readwind, ' seconds'
write(*,*) 'Timemanager: ', s_timemanager, ' seconds,', 'first timestep: ',s_firstt, 'seconds'
write(*,*) 'Write particle files: ', s_writepart, ' seconds'
......@@ -215,7 +231,6 @@ subroutine read_options_and_initialise_flexpart
!********************************
call readageclasses
! ! Allocate memory for windfields
! !*******************************
! call alloc_windfields
......@@ -393,6 +408,8 @@ subroutine initialise_particles
if (ipin.le.2) then
call readreleases
! needs to be called after maxspec is defined in readreleases or readinitconditions
if (ipout.ne.0) call readpartoptions
else
#ifdef USE_NCF
call readinitconditions_netcdf
......@@ -401,7 +418,7 @@ subroutine initialise_particles
#endif
endif
if (ipout.ne.0) call readpartoptions
if (iout.ne.0) then
call alloc_grid
......
......@@ -410,10 +410,16 @@ subroutine adv_above_pbl(itime,itimec,dxsave,dysave,ux,vy,tropop,nrand,ipart)
call get_settling(xts,yts,zts,nsp,part(ipart)%settling)
#ifdef ETA
call update_zeta_to_z(itime,ipart)
if ((ldirect.eq.1).and.(part(ipart)%z+part(ipart)%settling*dt.lt.0.)) then
part(ipart)%settling=-part(ipart)%z/dt
endif
call w_to_weta(itime,dt,part(ipart)%xlon,part(ipart)%ylat, &
part(ipart)%z,part(ipart)%zeta,part(ipart)%settling,weta_settling)
weta=weta+weta_settling
#else
if ((ldirect.eq.1).and.(part(ipart)%z+part(ipart)%settling*dt.lt.0.)) then
part(ipart)%settling=-part(ipart)%z/dt
endif
w=w+part(ipart)%settling
#endif
end if
......@@ -586,6 +592,9 @@ subroutine adv_in_pbl(itime,itimec, dxsave,dysave,dawsave,dcwsave, abovePBL, &
endif
if (density(nsp).gt.0.) then
call get_settling(xts,yts,zts,nsp,part(ipart)%settling) !bugfix
if ((ldirect.eq.1).and.(part(ipart)%z+part(ipart)%settling*dt.lt.0.)) then
part(ipart)%settling=-part(ipart)%z/dt
endif
w=w+part(ipart)%settling
end if
end if
......@@ -712,6 +721,10 @@ subroutine petterssen_corr(itime,ipart)
call update_z_to_zeta(itime+part(ipart)%idt,ipart)
zts=real(part(ipart)%z)
call get_settling(xts,yts,zts,nsp,part(ipart)%settling) !bugfix
if ((ldirect.eq.1).and. &
(part(ipart)%z+part(ipart)%settling*part(ipart)%idt.lt.0)) then
part(ipart)%settling=-part(ipart)%z/part(ipart)%idt
endif
call w_to_weta( &
itime+part(ipart)%idt, real(part(ipart)%idt), part(ipart)%xlon, &
part(ipart)%ylat, part(ipart)%z, part(ipart)%zeta, &
......@@ -721,6 +734,10 @@ subroutine petterssen_corr(itime,ipart)
!real(part(ipart)%zeta-part(ipart)%zeta_prev)/real(part(ipart)%idt*ldirect)
#else
call get_settling(xts,yts,zts,nsp,part(ipart)%settling)
if ((ldirect.eq.1).and. &
(part(ipart)%z+part(ipart)%settling*part(ipart)%idt.lt.0)) then
part(ipart)%settling=-part(ipart)%z/part(ipart)%idt
endif
w=w+part(ipart)%settling
#endif
end if
......
......@@ -6,14 +6,17 @@
! *
! June 1996 *
! *
! Last update:15 August 2013 IP *
! Update: 15 August 2013 IP *
! PS 19 Nov 2020: correct comment about lcw *
! Anne Tipka, Petra Seibert 2021-02: implement new interpolation *
! for precipitation according to #295 using 2 additional fields *
! *
!*******************************************************************************
module com_mod
use par_mod, only: dp, numpath, maxnests, &
numclass, maxcolumn, maxrand, numwfmem
numclass, maxcolumn, maxrand, numwfmem, numpf
implicit none
......@@ -21,10 +24,12 @@ module com_mod
!**********************************************************************************
type :: particleoptions
character(2) :: name
character(20) :: long_name
character(28) :: long_name
character(7) :: short_name
logical :: print
logical :: average=.false.
integer :: i_average=0
integer :: ncid
end type particleoptions
integer :: num_partopt=34
......@@ -159,24 +164,33 @@ module com_mod
! nageclass number of ageclasses for the age spectra calculation
! lage [s] ageclasses for the age spectra calculation
!ESO: Disable settling if more than 1 species per release point
logical :: lsettling=.true.
logical :: gdomainfill
! gdomainfill .T., if domain-filling is global, .F. if not
!ZHG SEP 2015 wheather or not to read clouds from GRIB
logical :: readclouds=.false.
!ESO DEC 2015 whether or not both clwc and ciwc are present (if so they are summed)
logical :: sumclouds=.false.
logical :: lcw=.false. ! ZHG Sep 2015 ! AT renamed
! whether or not cloud water data found in GRIB, overwritten if CW is found
!ESO: Disable settling if more than 1 species per release point
logical :: lsettling=.true.
logical :: lcwsum=.false. ! ESO Dec 2015 ! AT renamed
! whether or not both clwc and ciwc are present (if so they are summed)
logical,dimension(maxnests) :: readclouds_nest, sumclouds_nest
logical :: lprecint ! AT, PS 2021
! true if new interpolation using additional precip fields is used
logical,dimension(maxnests) :: lcw_nest=.false.
logical,dimension(maxnests) :: lcwsum_nest=.false.
logical,dimension(maxnests) :: lprecintn
!NIK 16.02.2015
integer(selected_int_kind(16)),allocatable,dimension(:) :: &
tot_blc_count,tot_inc_count
integer(selected_int_kind(16)),allocatable,dimension(:) :: icnt_belowcld, &
&icnt_incld
integer :: bcscheme ! AT 2021
! = 1 for below cloud scheme of Grythe et al 2017 (based on Laakso and Kyro)
! = 2 for below-cloud scheme of Tipka et al 2023 (based on Wang et al 2014)
! = 3 for below-cloud scheme of Tipka et al 2023 (based on modified Wang et al 2014)
!*********************************************************************
! Variables defining the release locations, released species and their
......@@ -527,8 +541,8 @@ contains
implicit none
integer :: stat
allocate( tot_blc_count(maxspec),tot_inc_count(maxspec),stat=stat)
if (stat.ne.0) error stop "Could not allocate tot_blc_count or tot_inc_count"
allocate( icnt_belowcld(maxspec),icnt_incld(maxspec),stat=stat)
if (stat.ne.0) error stop "Could not allocate cnt_belowcld or icnt_incld"
allocate( specnum(maxspec),decay(maxspec),weta_gas(maxspec), &
wetb_gas(maxspec),crain_aero(maxspec),csnow_aero(maxspec), &
ccn_aero(maxspec),in_aero(maxspec),ndia(maxspec), &
......@@ -552,8 +566,9 @@ contains
if (stat.ne.0) error stop "Could not allocate DRYDEPSPEC or WETDEPSPEC"
allocate( creceptor(numreceptor,maxspec),stat=stat)
if (stat.ne.0) error stop "Could not allocate creceptor"
tot_blc_count=0
tot_inc_count=0
icnt_belowcld=0
icnt_incld=0
end subroutine alloc_com
subroutine alloc_com_ndia
......@@ -565,7 +580,7 @@ contains
end subroutine alloc_com_ndia
subroutine dealloc_com
deallocate(tot_blc_count,tot_inc_count,specnum,decay,weta_gas,wetb_gas, &
deallocate(icnt_belowcld,icnt_incld,specnum,decay,weta_gas,wetb_gas, &
crain_aero,csnow_aero,ccn_aero,in_aero,reldiff,henry,f0,density,dquer, &
dsigma,ndia,vsetaver,cunningham,weightmolar,vset,schmi,fract,ri,rac,rcl, &
rgs,rlu,rm,dryvel,ohcconst,ohdconst,ohnconst,Fn,Fs,ks1,ks2,kn2,ishape, &
......
......@@ -853,11 +853,9 @@ subroutine redist(itime,ipart,ktop,ipconv,ithread)
loop1: do k = 1,nconvtop
! for backward runs use the transposed matrix
if (ldirect.eq.1) then
ffraction=ffraction+fmassfrac(levold,k,ithread) &
/totlevmass
ffraction=ffraction+fmassfrac(levold,k,ithread) / totlevmass
else
ffraction=ffraction+fmassfrac(k,levold,ithread) &
/totlevmass
ffraction=ffraction+fmassfrac(k,levold,ithread) / totlevmass
endif
if (rn.le.ffraction) then
levnew=k
......@@ -915,8 +913,10 @@ subroutine redist(itime,ipart,ktop,ipconv,ithread)
(tconv(levold,ithread)-tconv(levold-1,ithread)) &
*(pconv(levold-1,ithread)-phconv(levold,ithread))/ &
(pconv(levold-1,ithread)-pconv(levold,ithread))
! Bug fix: Added lsynctime to make units correct
sub_levold = sub(levold,ithread)/(1.-ga*sub(levold,ithread)*lsynctime/dpr(levold,ithread))
! LB: the units seem to not add up correctly, but adding lsynctime gives incorrect mixing
! in the lowest km and too many right above the ground
! sub_levold = sub(levold,ithread)/(1.-ga*sub(levold,ithread)*lsynctime/dpr(levold,ithread))
sub_levold = sub(levold,ithread)/(1.-ga*sub(levold,ithread)/dpr(levold,ithread))
wsub(levold,ithread)=-1.*sub_levold*r_air*temp_levold/(phconv(levold,ithread))
else
wsub(levold,ithread)=0.
......@@ -926,8 +926,8 @@ subroutine redist(itime,ipart,ktop,ipconv,ithread)
(tconv(levold+1,ithread)-tconv(levold,ithread)) &
*(pconv(levold,ithread)-phconv(levold+1,ithread))/ &
(pconv(levold,ithread)-pconv(levold+1,ithread))
! Bug fix: Added lsynctime to make units correct
sub_levold1 = sub(levold+1,ithread)/(1.-ga*sub(levold+1,ithread)*lsynctime/dpr(levold+1,ithread))
!sub_levold1 = sub(levold+1,ithread)/(1.-ga*sub(levold+1,ithread)*lsynctime/dpr(levold+1,ithread))
sub_levold1 = sub(levold+1,ithread)/(1.-ga*sub(levold+1,ithread)/dpr(levold+1,ithread))
wsub(levold+1,ithread)=-1.*sub_levold1*r_air*temp_levold1/ &
(phconv(levold+1,ithread))
......
......@@ -23,6 +23,9 @@ subroutine caldate(juliandate,yyyymmdd,hhmiss)
! *
! AUTHOR: Andreas Stohl (21 January 1994), adapted from Numerical Recipes*
! *
! PS 2020-07-27: add a check to avoid giving back 240000 for hhmiss *
! *
! *
! Variables: *
! dd Day *
! hh Hour *
......@@ -50,6 +53,7 @@ subroutine caldate(juliandate,yyyymmdd,hhmiss)
integer,parameter :: igreg=2299161
julday=int(juliandate)
! PS check to avoid 240000 as hhmiss:
if ((juliandate-julday)*86400._dp .ge. 86399.5_dp) then
juliandate = juliandate + juliandate-julday-86399.5_dp/86400._dp
julday=int(juliandate)
......
......@@ -85,16 +85,19 @@ subroutine assignland
implicit none
integer :: ix,jy,k,l,li,nrefine,iix,jjy
integer :: ix,jy,k,l,li,nrefine,iix,jjy,stat
integer,parameter :: lumaxx=1200,lumaxy=600
integer,parameter :: xlon0lu=-180,ylat0lu=-90
real,parameter :: dxlu=0.3
real :: xlon,ylat,sumperc,p,xi,yj
real :: xlandusep(lumaxx,lumaxy,numclass)
real,allocatable,dimension(:,:,:) :: xlandusep
! character*2 ck
if (.not.DRYDEP) return
allocate( xlandusep(lumaxx,lumaxy,numclass), stat=stat)
if (stat.ne.0) error stop "Could not allocate xlandusep in assignland"
do ix=1,lumaxx
do jy=1,lumaxy
do k=1,numclass
......@@ -1302,10 +1305,34 @@ subroutine partdep(nc,density,fract,schmi,vset,ra,ustar,nyl,rhoa,vdep_tmp)
if (ustar.gt.eps) then
if (ishape(ic).eq.0) then
reynolds=dquer(ic)/1.e6*vset(ic,j)/nyl
settling_old=-1.0*vset(ic,j)
do i=1,20
if (reynolds.le.0.02) then
c_d=(24.0/reynolds)
else ! Clif and Gauvin scheme is used
c_d=(24.0/reynolds)*(1+0.15*(reynolds**0.687))+ &
0.42/(1.0+42500.0/(reynolds**1.16))
endif
! Settling velocity of a particle is defined by the Newton's impact law:
settling=-1.* &
sqrt(4.*ga*dquer(ic)/1.e6*density(ic)*cunningham(ic)/ &
(3.*c_d*rhoa))
if (abs((settling-settling_old)/settling).lt.0.01) exit
reynolds=dquer(ic)/1.e6*abs(settling)/nyl
settling_old=settling
end do
! Stokes number for each diameter interval
!*****************************************
! Use this stokes number for different shapes
stokes=vset(ic,j)/ga*ustar*ustar/nyl
stokes=abs(settling)/ga*ustar*ustar/nyl
alpha=-3./stokes
! Deposition layer resistance
......@@ -1317,7 +1344,7 @@ subroutine partdep(nc,density,fract,schmi,vset,ra,ustar,nyl,rhoa,vdep_tmp)
rdp=1./((schmi(ic,j)+10.**alpha)*ustar)
endif
vdepj=vset(ic,j)+1./(ra+rdp+ra*rdp*vset(ic,j))
vdepj=abs(settling)+1./(ra+rdp+ra*rdp*abs(settling))
else ! Daria Tatsii: Drag coefficient scheme by Bagheri & Bonadonna 2016
! Settling velocities of other shapes
......
......@@ -155,6 +155,7 @@ subroutine getfields(itime,nstop)
return
endif
if ((ldirect*memtime(1).le.ldirect*itime).and. &
(ldirect*memtime(2).gt.ldirect*itime)) then
......@@ -166,7 +167,6 @@ subroutine getfields(itime,nstop)
else if ((ldirect*memtime(2).le.ldirect*itime).and. &
(memtime(2).ne.999999999)) then
! Current time is after 2nd wind field
! -> Resort wind field pointers, so that current time is between 1st and 2nd
!***************************************************************************
......@@ -256,11 +256,13 @@ subroutine getfields(itime,nstop)
call readwind_nest(indj,memind(1),uuhn,vvhn,wwhn)
call calcpar(memind(1))
call calcpar_nest(memind(1))
if (metdata_format.eq.GRIBFILE_CENTRE_ECMWF) then
call verttransform_ecmwf(memind(1),uuh,vvh,wwh,pvh)
else
call verttransform_gfs(memind(1),uuh,vvh,wwh,pvh)
end if
call verttransform_nest(memind(1),uuhn,vvhn,wwhn,pvhn)
memtime(1)=wftime(indj)
memind(2)=2
......@@ -886,7 +888,7 @@ subroutine calcpar(n)
! *
! 21 May 1995 *
! *
! ------------------------------------------------------------------ *
!*****************************************************************************
! Petra Seibert, Feb 2000: *
! convection scheme: *
! new variables in call to richardson *
......@@ -901,7 +903,9 @@ subroutine calcpar(n)
! - Merged calcpar and calcpar_gfs into one routine using if-then *
! for meteo-type dependent code *
!*****************************************************************************
! Changes Anne Tipka June 2023: *
! sum up precipitation fields over number of available fields in a single *
! time interval (newWetDepoScheme) *
!*****************************************************************************
! *
! Variables: *
......@@ -1059,7 +1063,7 @@ subroutine calcpar(n)
call getvdep(n,ix,jy,ustar(ix,jy,1,n), &
tt2(ix,jy,1,n),ps(ix,jy,1,n),1./oli(ix,jy,1,n), &
ssr(ix,jy,1,n),rh,lsprec(ix,jy,1,n)+convprec(ix,jy,1,n), &
ssr(ix,jy,1,n),rh,sum(lsprec(ix,jy,1,:,n))+sum(convprec(ix,jy,1,:,n)), &
sd(ix,jy,1,n),vd)
do i=1,nspec
......@@ -1159,7 +1163,7 @@ subroutine calcpar_nest(n)
! *
! 8 February 1999 *
! *
! ------------------------------------------------------------------ *
!*****************************************************************************
! Petra Seibert, Feb 2000: *
! convection scheme: *
! new variables in call to richardson *
......@@ -1171,6 +1175,10 @@ subroutine calcpar_nest(n)
! Unified ECMWF and GFS builds *
! Marian Harustak, 12.5.2017 *
!*****************************************************************************
! Changes Anne Tipka June 2023: *
! sum up precipitation fields over number of available fields in a single *
! time interval (newWetDepoScheme) *
!*****************************************************************************
! *
! Variables: *
! n temporal index for meteorological fields (1 to 3) *
......@@ -1301,8 +1309,8 @@ subroutine calcpar_nest(n)
call getvdep_nest(n,ix,jy,ustarn(ix,jy,1,n,l), &
tt2n(ix,jy,1,n,l),psn(ix,jy,1,n,l),1./olin(ix,jy,1,n,l), &
ssrn(ix,jy,1,n,l),rh,lsprecn(ix,jy,1,n,l)+ &
convprecn(ix,jy,1,n,l),sdn(ix,jy,1,n,l),vd,l)
ssrn(ix,jy,1,n,l),rh,sum(lsprecn(ix,jy,1,:,n,l))+ &
sum(convprecn(ix,jy,1,:,n,l)),sdn(ix,jy,1,n,l),vd,l)
do i=1,nspec
vdepn(ix,jy,i,n,l)=vd(i)
......
......@@ -98,14 +98,12 @@ subroutine releaseparticles(itime)
!real xaux,yaux,zaux,ran1,rfraction,xmasssave(maxpoint)
real :: xaux,yaux,zaux,rfraction
real :: topo,r,t
real :: dp1,dp2,xlonav,timecorrect(maxspec),press,pressold
real :: presspart,average_timecorrect
real :: xlonav,timecorrect(maxspec)
real :: average_timecorrect
integer :: itime,numrel,i,j,k,ipart,minpart
integer :: kz,istart,iend,totpart
integer :: istart,iend,totpart,iterm_index
integer :: nweeks,ndayofweek,nhour,jjjjmmdd,ihmmss,mm
real(kind=dp) :: julmonday,jul,jullocal,juldiff
real,parameter :: eps2=1.e-6
integer :: idummy = -7
!save idummy,xmasssave
......@@ -147,6 +145,7 @@ subroutine releaseparticles(itime)
endif
call get_totalpart_num(istart)
if (ipin.le.1 .and. ipout.eq.0) call rewrite_iterm()
minpart=1
do i=1,numpoint
if ((itime.ge.ireleasestart(i)).and. &! are we within release interval?
......@@ -214,11 +213,14 @@ subroutine releaseparticles(itime)
zaux=zpoint2(i)-zpoint1(i)
if (ipin.le.1 .and. ipout.eq.0) then
totpart = numrel-(count%allocated-count%alive)
call rewrite_iterm()
totpart = numrel-count%iterm_max
if (totpart.gt.0) call alloc_particles(totpart)
call rewrite_iterm()
iterm_index=1
endif
do j=1,numrel ! loop over particles to be released this time
call get_newpart_index(ipart)
call get_newpart_index(ipart,iterm_index)
call spawn_particle(itime, ipart)
! Particle coordinates are determined by using a random position within the release volume
......@@ -266,6 +268,64 @@ subroutine releaseparticles(itime)
! Interpolation of topography and density
!****************************************
! Transform the verticle particle position from pressure or sea level to above ground
! if necessary
!************************************************************************************
call kindz_to_z(ipart)
#ifdef ETA
call z_to_zeta(itime,part(ipart)%xlon,part(ipart)%ylat,part(ipart)%z,part(ipart)%zeta)
part(ipart)%etaupdate = .true. ! The z(meter) coordinate is up to date
part(ipart)%meterupdate = .true.
#endif
call init_mass_conversion(ipart,i)
call get_totalpart_num(numpart)
end do ! numrel
endif ! releasepoint
end do ! numpoint
if (ipin.le.1 .and. ipout.eq.0) call rewrite_iterm()
call get_totalpart_num(iend)
! NetCDF only: write initial positions of new particles
! #ifdef USE_NCF
! if ((iend-istart.gt.0).and.(ipout.ge.1)) then
! call wrt_part_initialpos(itime,istart,iend)
! call output_particles(itime,.true.)
! endif
! #endif
return
! 996 continue
! write(*,*) '#####################################################'
! write(*,*) '#### FLEXPART MODEL SUBROUTINE RELEASEPARTICLES: ####'
! write(*,*) '#### ####'
! write(*,*) '#### ERROR - TOTAL NUMBER OF PARTICLES REQUIRED ####'
! write(*,*) '#### EXCEEDS THE MAXIMUM ALLOWED NUMBER. REDUCE ####'
! write(*,*) '#### EITHER NUMBER OF PARTICLES PER RELEASE POINT####'
! write(*,*) '#### OR REDUCE NUMBER OF RELEASE POINTS. ####'
! write(*,*) '#####################################################'
! stop
end subroutine releaseparticles
subroutine kindz_to_z(ipart)
use point_mod
use xmass_mod
use output_mod
use interpol_mod
implicit none
integer,intent(in) :: ipart
integer :: kz
real :: dp1,dp2,press,presspart,pressold,topo,r,t
real,parameter :: eps2=1.e-6
! Determine the nest we are in
!*****************************
call find_ngrid(part(ipart)%xlon,part(ipart)%ylat)
......@@ -275,21 +335,14 @@ subroutine releaseparticles(itime)
call find_grid_indices(real(part(ipart)%xlon),real(part(ipart)%ylat))
call find_grid_distances(real(part(ipart)%xlon),real(part(ipart)%ylat))
if (ngrid.gt.0) then
topo=p1*oron(ix ,jy ,ngrid) &
+ p2*oron(ixp,jy ,ngrid) &
+ p3*oron(ix ,jyp,ngrid) &
+ p4*oron(ixp,jyp,ngrid)
else
topo=p1*oro(ix ,jy) &
+ p2*oro(ixp,jy) &
+ p3*oro(ix ,jyp) &
+ p4*oro(ixp,jyp)
endif
! If starting height is in pressure coordinates, retrieve pressure profile and convert zpart1 to meters
if (kindz(part(ipart)%npoint).eq.1) return ! Nothing needs to happen
! If starting height is in pressure coordinates, retrieve pressure profile and
! convert zpart1 to meters
!*****************************************************************************
if (kindz(i).eq.3) then
if (kindz(part(ipart)%npoint).eq.3) then
presspart=real(part(ipart)%z)
do kz=1,nz
if (ngrid.gt.0) then
......@@ -327,55 +380,37 @@ subroutine releaseparticles(itime)
endif
pressold=press
end do
endif
! If release positions are given in meters above sea level, subtract the
! topography from the starting height
!***********************************************************************
if (kindz(i).eq.2) call update_z(ipart,-topo)
else if (kindz(part(ipart)%npoint).eq.2) then
if (ngrid.gt.0) then
topo=p1*oron(ix ,jy ,ngrid) &
+ p2*oron(ixp,jy ,ngrid) &
+ p3*oron(ix ,jyp,ngrid) &
+ p4*oron(ixp,jyp,ngrid)
else
topo=p1*oro(ix ,jy) &
+ p2*oro(ixp,jy) &
+ p3*oro(ix ,jyp) &
+ p4*oro(ixp,jyp)
endif
call update_z(ipart,-topo)
endif
if (part(ipart)%z.lt.eps2) call set_z(ipart,eps2) ! Minimum starting height is eps2
if (part(ipart)%z.gt.height(nz)-0.5) &
call set_z(ipart,height(nz)-0.5) ! Maximum starting height is uppermost level - 0.5 meters
#ifdef ETA
call z_to_zeta(itime,part(ipart)%xlon,part(ipart)%ylat,part(ipart)%z,part(ipart)%zeta)
part(ipart)%etaupdate = .true. ! The z(meter) coordinate is up to date
part(ipart)%meterupdate = .true.
#endif
call init_mass_conversion(ipart,i)
call get_totalpart_num(numpart)
end do ! numrel
endif ! releasepoint
end do ! numpoint
call get_totalpart_num(iend)
! NetCDF only: write initial positions of new particles
#ifdef USE_NCF
if ((iend-istart.gt.0).and.(ipout.ge.1)) then
call wrt_part_initialpos(itime,istart,iend)
call output_particles(itime,.true.)
if (ipin.eq.3 .or. ipin.eq.4) then
if (part(ipart)%z.gt.zpoint2(part(ipart)%npoint)) &
zpoint2(part(ipart)%npoint)=real(part(ipart)%z)
if (part(ipart)%z.lt.zpoint1(part(ipart)%npoint)) &
zpoint1(part(ipart)%npoint)=real(part(ipart)%z)
endif
#endif
return
! 996 continue
! write(*,*) '#####################################################'
! write(*,*) '#### FLEXPART MODEL SUBROUTINE RELEASEPARTICLES: ####'
! write(*,*) '#### ####'
! write(*,*) '#### ERROR - TOTAL NUMBER OF PARTICLES REQUIRED ####'
! write(*,*) '#### EXCEEDS THE MAXIMUM ALLOWED NUMBER. REDUCE ####'
! write(*,*) '#### EITHER NUMBER OF PARTICLES PER RELEASE POINT####'
! write(*,*) '#### OR REDUCE NUMBER OF RELEASE POINTS. ####'
! write(*,*) '#####################################################'
! stop
end subroutine releaseparticles
end subroutine kindz_to_z
subroutine init_mass_conversion(ipart,ipoint)
! For special simulations, multiply particle concentration air density;
......@@ -411,8 +446,8 @@ subroutine init_mass_conversion(ipart,ipoint)
if ((ind_rel .eq. 1).or.(ind_rel .eq. 3).or.(ind_rel .eq. 4)) then
! Interpolate the air density
!****************************
! Interpolate the air density, horizontal values are computed in kindz_to_z
!**************************************************************************
call find_z_level_meters(real(part(ipart)%z))
call find_vert_vars_lin(height,real(part(ipart)%z),indz,dz1,dz2,lbounds)
......@@ -436,7 +471,7 @@ subroutine init_mass_conversion(ipart,ipoint)
! Multiply "mass" (i.e., mass mixing ratio in forward runs) with density
!********************************************************************
!***********************************************************************
mass(ipart,:)=mass(ipart,:)*rhoout
mass_init(ipart,:)=mass(ipart,:)
endif
......@@ -898,9 +933,9 @@ subroutine init_domainfill
call alloc_domainfill
nx_we(1)=max(int(xpoint1(1)),0)
nx_we(2)=min((int(xpoint2(1))+1),nxmin1)
nx_we(2)=min((ceiling(xpoint2(1))),nxmin1)
ny_sn(1)=max(int(ypoint1(1)),0)
ny_sn(2)=min((int(ypoint2(1))+1),nymin1)
ny_sn(2)=min((ceiling(ypoint2(1))),nymin1)
! For global simulations (both global wind data and global domain-filling),
! set a switch, such that no boundary conditions are used
......@@ -1358,7 +1393,7 @@ subroutine boundcond_domainfill(itime,loutend)
real :: windx,rhox
real :: deltaz,boundarea,fluxofmass
integer :: ixm,ixp,jym,jyp,indzm,mm
integer :: ixm,ixp,jym,jyp,indzm,mm,iterm_index
real :: pvpart,ddx,ddy,rddx,rddy,p1,p2,p3,p4,y1(2),yh1(2)
integer :: idummy = -11
......@@ -1382,6 +1417,8 @@ subroutine boundcond_domainfill(itime,loutend)
! Terminate trajectories that have left the domain, if domain-filling
! trajectory calculation domain is not global
!********************************************************************
if (ipin.le.1 .and. ipout.eq.0) call rewrite_iterm()
iterminate=0
do i=1,count%allocated
if (.not. part(i)%alive) cycle
......@@ -1417,6 +1454,7 @@ subroutine boundcond_domainfill(itime,loutend)
! ithread = 0
! #endif
ithread=0
iterm_index=1
! !$OMP DO
do jy=ny_sn(1),ny_sn(2)
......@@ -1530,7 +1568,7 @@ subroutine boundcond_domainfill(itime,loutend)
do m=1,mmass
!THIS WILL CAUSE PROBLEMS WITH OMP! because of dynamical allocation
call get_newpart_index(ipart)
call get_newpart_index(ipart,iterm_index)
call spawn_particle(itime, ipart)
! Assign particle positions
......@@ -1745,7 +1783,7 @@ subroutine boundcond_domainfill(itime,loutend)
endif
do m=1,mmass
call get_newpart_index(ipart)
call get_newpart_index(ipart,iterm_index)
call spawn_particle(itime, ipart)
! Assign particle positions
......@@ -1836,6 +1874,7 @@ subroutine boundcond_domainfill(itime,loutend)
end do ! north south
! !$OMP END DO
! !$OMP END PARALLEL
if (ipin.le.1 .and. ipout.eq.0) call rewrite_iterm()
numparticlecount = numparticlecount_tmp
! If particles shall be dumped, then accumulated masses at the domain boundaries
! must be dumped, too, to be used for later runs
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Please register or to comment