Initial commit

LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2021 Christopher Polster, Volkmar Wirth
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

Makefile

@@ -0,0 +1,167 @@
+# Number of samples for the Monte-Carlo confidence estimation
+NSAMPLES := 1000
+
+# How many ensemble members in the sampled ensemble?
+SAMPLE_SIZE := 75
+
+# Threshold for the correlation distribution width below which datapoints are
+# considered to be statistically significant
+SIGNIFICANCE_LEVEL := 0.1
+
+# Determine filename suffix for C-extensions to installed Python
+PYEXT := $(shell python3-config --extension-suffix)
+
+# Names of all figures to be plotted
+FIGURES := \
+		figures/forecast-combination.pdf \
+		figures/event24-reanalysis+nh18fig5.pdf \
+		figures/event24-plume.pdf \
+		figures/event24-maps+hovmoeller.pdf \
+		figures/event24-cluster+scatter.pdf \
+		figures/event24-maps-separate.pdf \
+		figures/event24-cluster-f2.pdf \
+		figures/event24-nh18fig4.pdf
+
+
+.PHONY: all clean local
+
+all: $(FIGURES)
+
+clean:
+	rm -f $(FIGURES)
+	rm -f scripts/common/extensions/*.o
+	rm -f scripts/common/_*.c
+	rm -f scripts/common/_*.o
+	rm -f scripts/common/_*.so
+	rm -f scripts/hn2016_falwa/*.so
+
+
+# Python scipts for data processing and plotting, C-extensions
+
+scripts/calculate-lwa.py: scripts/common/regrid.py scripts/hn2016_falwa/oopinterface.py
+	touch $@
+
+scripts/plot.py: scripts/common/plotting.py scripts/common/esa.py scripts/common/metrics.py scripts/common/texify.py
+	touch $@
+
+scripts/common/metrics.py: scripts/common/texify.py
+	touch $@
+
+scripts/common/esa.py: scripts/common/_esa$(PYEXT)
+	touch $@
+
+scripts/common/regrid.py: scripts/common/_regrid$(PYEXT)
+	touch $@
+
+scripts/common/_%$(PYEXT): scripts/common/extensions/build_%.py scripts/common/extensions/%.c
+	python3 $<
+
+
+# hn2016_falwa package
+
+scripts/hn2016_falwa/oopinterface.py: \
+		scripts/hn2016_falwa/constant.py \
+		scripts/hn2016_falwa/utilities.py \
+		scripts/hn2016_falwa/interpolate_fields$(PYEXT) \
+		scripts/hn2016_falwa/compute_lwa_and_barotropic_fluxes$(PYEXT) \
+		scripts/hn2016_falwa/compute_reference_states$(PYEXT)
+	touch $@
+
+scripts/hn2016_falwa/%$(PYEXT): scripts/hn2016_falwa/%.f90
+	f2py -c -m $* $<
+	mv $(notdir $@) scripts/hn2016_falwa
+
+
+
+# Data Processing (Part 1):
+# Calculating LWA and related quantities with the hn2016_falwa package.
+
+# Write processed data to netcdf files:
+event24_ana := data/ERA-2016-12-10-to-2016-12-24-lwa-qg.nc
+event24_ens := data/ENS-2016-12-10T00Z-lwa-qg.nc data/ENS-2016-12-10T12Z-lwa-qg.nc data/ENS-2016-12-11T00Z-lwa-qg.nc
+# Don't remove these intermediate files
+.SECONDARY: $(event24_ana) $(event24_ens)
+
+# Reanalysis data
+data/ERA-%-lwa-qg.nc: data/ERA-%-uvt.nc scripts/calculate-lwa.py
+	python3 -m scripts.calculate-lwa --half-resolution $<
+
+# Ensemble data
+data/ENS-%-lwa-qg.nc: scripts/calculate-lwa.py data/IFS-ENS/ENS-%/ENS-*-[uvt].nc
+	python3 -m scripts.calculate-lwa --half-resolution $(filter-out $<,$^)
+
+
+
+# Data Processing (Part 2):
+# Analysis and plotting of figures.
+
+# Figure 1: Ensemble combination
+figures/forecast-combination.pdf: figures/forecast-combination/forecast-combination.tex
+	cd $(dir $<) && lualatex $(notdir $<)
+	mv $(<:.tex=.pdf) $@
+
+# Figure 2: Reanalysis overview
+# Figure 3: Ensemble target metric overview
+# Figure 4: ESA maps and Hovmöller with full flux
+# Figure 5: Cluster and scatter analysis of upstream region
+figures/event24-%.pdf: scripts/plot.py $(event24_ana) $(event24_ens)
+	python3 -m scripts.plot $* 2016-12-18T00,65,10,35,-30 $(event24_ens) \
+			--reanalysis $(event24_ana) \
+			--output $@ \
+			--nsamples $(NSAMPLES) \
+			--sample-size $(SAMPLE_SIZE) \
+			--significance-level $(SIGNIFICANCE_LEVEL) \
+			--delta-hours "0,-24,-48,-72,-96" \
+			--source "f1+f2+f3" \
+			--scatter 2016-12-16T00,65,-90,35,-15 \
+			--hovmoeller-extent "target" \
+			--end 2016-12-24T00
+
+# Figure 6: ESA maps with flux terms
+figures/event24-maps-separate.pdf: figures/event24-maps-separate/event24-maps-separate.tex \
+		figures/event24-maps-separate/f1f3.pdf figures/event24-maps-separate/f2.pdf
+	cd $(dir $<) && lualatex $(notdir $<)
+	mv $(<:.tex=.pdf) $@
+
+figures/event24-maps-separate/f1f3.pdf: scripts/plot.py $(event24_ana) $(event24_ens)
+	python3 -m scripts.plot maps4 2016-12-18T00,65,10,35,-30 $(event24_ens) \
+			--reanalysis $(event24_ana) \
+			--output $@ \
+			--nsamples $(NSAMPLES) \
+			--sample-size $(SAMPLE_SIZE) \
+			--significance-level $(SIGNIFICANCE_LEVEL) \
+			--delta-hours "24, 0,-24, -48" \
+			--source "f1+f3" \
+			--panel-letters "abcd"
+
+figures/event24-maps-separate/f2.pdf: scripts/plot.py $(event24_ana) $(event24_ens)
+	python3 -m scripts.plot maps4 2016-12-18T00,65,10,35,-30 $(event24_ens) \
+			--reanalysis $(event24_ana) \
+			--output $@ \
+			--nsamples $(NSAMPLES) \
+			--sample-size $(SAMPLE_SIZE) \
+			--significance-level $(SIGNIFICANCE_LEVEL) \
+			--delta-hours "24, 0,-24, -48" \
+			--source "f2" \
+			--panel-letters "efgh"
+
+# Figure 7: Cluster analysis of F2 in block
+figures/event24-cluster-f2.pdf: scripts/plot.py $(event24_ana) $(event24_ens)
+	python3 -m scripts.plot cluster 2016-12-18T00,65,10,35,-30 $(event24_ens) \
+			--reanalysis $(event24_ana) \
+			--output $@ \
+			--nsamples $(NSAMPLES) \
+			--sample-size $(SAMPLE_SIZE) \
+			--significance-level $(SIGNIFICANCE_LEVEL) \
+			--source "f2"
+
+# Figure 8: Flux-LWA relationship in block (similar to NH18's Fig. 4)
+figures/event24-nh18fig4.pdf: scripts/plot.py $(event24_ana) $(event24_ens)
+	python3 -m scripts.plot nh18fig4 2016-12-18T00,65,10,35,-30 $(event24_ens) \
+			--reanalysis $(event24_ana) \
+			--output $@ \
+			--nsamples $(NSAMPLES) \
+			--sample-size $(SAMPLE_SIZE) \
+			--significance-level $(SIGNIFICANCE_LEVEL) \
+			--scatter 2016-12-18T00,50,-5,40,-15
+

README.md

@@ -0,0 +1,78 @@
+# Block-Dec2016-LWA-ESA
+
+Code to produce all figures of the article **Evidence for a "Traffic Jam" Onset of Blocked Flow from Ensemble Sensitivity Analysis**.
+
+- Research by [Christopher Polster](https://dynmet.ipa.uni-mainz.de/christopher-polster/) and [Volkmar Wirth](https://dynmet.ipa.uni-mainz.de/volkmar-wirth/).
+- Software by Christopher Polster.
+- [hn2016_falwa](https://github.com/csyhuang/hn2016_falwa/) package by [Clare S. Y. Huang](https://csyhuang.github.io/).
+
+[LICENSE](LICENSE). The contents of `scripts/hn2016_falwa` are re-distributed under the terms of their [LICENSE](scripts/hn2016_falwa/LICENSE.txt).
+
+
+## How To Run
+
+Clone this repository:
+
+    $ git clone https://github.com/wavestoweather/Block-Dec2016-LWA-ESA.git
+
+### Data Preparation
+
+Use the scripts and instructions in the [`data`](data) directory to obtain the required input data. Before continuing, the `data` folder should look like this:
+
+    data
+    + IFS-ENS
+    | + ENS-2016-12-10T00Z
+    | | + ENS-2016-12-10T00Z-t.nc
+    | | + ENS-2016-12-10T00Z-u.nc
+    | | + ENS-2016-12-10T00Z-v.nc
+    | + ENS-2016-12-10T12Z
+    | | + ENS-2016-12-10T12Z-t.nc
+    | | + ENS-2016-12-10T12Z-u.nc
+    | | + ENS-2016-12-10T12Z-v.nc
+    | + ENS-2016-12-11T00Z
+    | | + ENS-2016-12-11T00Z-t.nc
+    | | + ENS-2016-12-11T00Z-u.nc
+    | | + ENS-2016-12-11T00Z-v.nc
+    | + job.sh
+    + download-ERA5.py
+    + ERA-2016-12-10-to-2016-12-24-uvt.nc
+    + README.md
+
+### Software Requirements
+
+The following software needs to be available to run the data processing and plotting scripts:
+
+- make
+- Fortran compiler (for building the `hn2016_falwa` Python extensions)
+- C compiler with OpenMP (for building Python extensions)
+- Python 3 with packages listed in [`requirements.txt`](requirements.txt)
+- lualatex with fontspec, amsmath, tikz, graphicx (Fig. 1 and to combine panels for Fig. 6)
+
+### Data Analysis and Plots
+
+Run
+
+    $ make
+
+from the top-level of the repository. This will compile the Python extensions, compute the intermediate PV, LWA and flux fields and create all plots. Use `--dry-run` to check the sequence of commands to be executed by `make` first. Use `-jN` to run `N` jobs in parallel if you have sufficient CPU and memory resources to do so.
+
+### Output
+
+Figures will appear in the `figures` directory:
+
+- Fig. 1: `figures/forecast-combination.pdf`
+- Fig. 2: `figures/event24-reanalysis+nh18fig5.pdf`
+- Fig. 3: `figures/event24-plume.pdf`
+- Fig. 4: `figures/event24-maps+hovmoeller.pdf`
+- Fig. 5: `figures/event24-cluster+scatter.pdf`
+- Fig. 6: `figures/event24-maps-separate.pdf`
+- Fig. 7: `figures/event24-cluster-f2.pdf`
+- Fig. 8: `figures/event24-nh18fig4.pdf`, does not contain the fits of [Nakamura and Huang (2018)](https://doi.org/10.1126/science.aat0721) due to licensing restrictions
+
+Because the statistical significance test is based on a randomized procedure, results may vary slightly each time the plots are created.
+
+
+## Acknowledgements
+
+This research project has been carried out within the Transregional Collaborative Research Center [SFB/TRR 165](https://www.dfg.de/en/funded_projects/current_projects_programmes/list/projectdetails/index.jsp?id=257899354) "Waves to Weather" funded by the German Science Foundation (DFG). https://www.wavestoweather.de/
+

data/IFS-ENS/job.sh

@@ -0,0 +1,118 @@
+#!/bin/bash
+
+# /!\ FIXME marks places where a custom path, etc. must be specified
+
+#SBATCH --qos=normal
+        # Specifies that your job will run in the queue (Quality Of
+        # Service) normal, which, if this option is not specified, is
+        # the default.
+
+#SBATCH --job-name=request-ens
+        # Assigns the specified name to the request
+
+#SBATCH --output=stdout/request-ens.%j.%N.stdout.log
+        # Specifies the name and location of STDOUT where %N is the
+        # node where the job runs and %j the job-id. The file will be 
+        # written in the workdir directory if it is a relative path.
+        # If not given, the default is slurm-%j.out in the workdir.
+
+#SBATCH --error=stderr/request-ens.%j.%N.stderr.out
+        # Specifies the name and location of STDERR where %N is the
+        # node where the job runs and %j the job-id. The file will be 
+        # written in the workdir directory if it is a relative path. 
+        # If not specified, STDERR will be written to the output file 
+        # defined above, or otherwise to slurm-%j.out in the workdir. 
+
+#SBATCH --workdir=FIXME
+        # Sets the working directory of the batch script before it is
+        # executed.
+
+#SBATCH --mail-type=END
+#SBATCH --mail-user=FIXME
+        # Specifies that an email should be sent in case the job fails.
+        # Other options include BEGIN, END, REQUEUE and ALL (any state
+        # change).
+
+# Take time and date from command line
+DATE=$1
+TIME=$2
+if [ -z "$DATE" ] || [ -z "$TIME" ]
+then
+	>&2 echo "Must specify arguments YYYY-MM-DD HH."
+	exit 1
+fi
+
+# Output directory FIXME
+TRGD="FIXME/IFS-ENS/ENS-${DATE}T${TIME}Z"
+# Output filename prefix
+TRGF="${TRGD}/ENS-${DATE}T${TIME}Z-"
+# Request file name
+RQST="request-ens-${DATE}T${TIME}Z"
+
+mkdir "${TRGD}"
+
+# Wind data on pressure level 400 is not available for older forecasts.
+# To avoid the MARS request failing, don't include it in the level list.
+# Forecasts earlier than middle of 2008 are missing even more levels
+# (I only know that February 2008 has less than September 2008).
+if (( "${DATE:0:4}" >= 2010 )); then
+	LVLS="10/50/100/200/250/300/400/500/700/850/925/1000"
+elif (( "${DATE:0:4}" == 2009 || ( "${DATE:0:4}" == 2008 && ( "${DATE:5:1}" == 1 || "${DATE:6:1}" >= 6 ) ) )); then
+	LVLS="10/50/100/200/250/300/500/700/850/925/1000"
+else
+    LVLS="200/250/300/500/700/850/925/1000"
+fi
+
+
+# 130.128  Temperature
+# 131      U-component of wind
+# 132      V-component of wind
+
+cat > "${RQST}" <<EOF
+retrieve,
+    class=od,
+    stream=enfo,
+    expver=1,
+    type=pf,
+    grid=1.0/1.0,
+    number=1/TO/50,
+    step=0/TO/288/BY/6,
+    area=90/-180/0/180,
+    date=${DATE},
+    time=${TIME}:00:00,
+    param=130.128,
+    levtype=pl,
+    levelist=${LVLS},
+    target="${TRGF}t.grb"
+retrieve,
+    param=131,
+    levtype=pl,
+    levelist=${LVLS},
+    target="${TRGF}u.grb"
+retrieve,
+    param=132,
+    levtype=pl,
+    levelist=${LVLS},
+    target="${TRGF}v.grb"
+EOF
+
+exitstatus=0
+
+mars "${RQST}"
+
+if [ $? != 0 ]
+then
+    echo " The MARS request failed."
+    echo
+    exitstatus=1
+fi
+
+# Convert files to NetCDF and remove original GRIB files
+GRIBS="${TRGD}/*.grb"
+for GRB in ${GRIBS}
+do
+	grib_to_netcdf -u time -o "${GRB%.grb}.nc" "${GRB}"  && rm -f "${GRB}"
+done
+
+exit "$exitstatus"
+