Lorenz-96 model with PDAF

The implementation of the Lorenz-96 model coupled to PDAF is in the directory models/lorenz96/ of the PDAF package. Provided is a full implementation of PDAF with the nonlinear Lorenz96 model (E. N. Lorenz (1996) Predictability - a problem partly solved. Proceedings Seminar on Predictability, ECMWF, READING, UK) providing the possibility to assess various filter and smoother methods. We used this implementation for different publications in which we studied the behavior of different data assimilation methods.

Next to the implementation of the Lorenz-96 model with PDAF, the test case provides tool programs and scripts (in the subdirectory tools/) that allow to run a test case and to display the outputs. Note, that this implementation runs without parallelization.

Running the test case

Runnning a data assimilation experiment with the Lorenz-96 model is a two step process: First one runs the model without PDAF to generate a file holding the trajectory of a forward run. Then, one generates a file with observations and a file holding information on the covariance matrix (EOFs and singular values, see the page on ensemble generation). The second file will be used to generate the initial ensemble with second-order exact sampling. Afterwards, one compiles the Lorenz-96 model with activated coupling to PDAF and runs the data assimilation experiments.

1. Compile and run the forward model without assimilation

First, we compile the Lorenz-96 model without assimilation. For this we need to deactivate the code parts that couple to PDAF. For this, change in the file make.arch/linux_gfortran.h the line

CPP_DEFS = -DUSE_PDAF

to

CPP_DEFS = #-DUSE_PDAF

Now build the forward model program with

  cd models/lorenz96
  make lorenz_96 PDAF_ARCH=linux_gfortran

You can replace linux_gfortran by any other make include file from make.arch/, e.g. specify macos_gfortran_openmpi for compiling on MacOS. The executable is generated in the directory models/lorenz96.

Note: The implementation uses the NetCDF library for file outputs. If the compilation above fails, please ensure the netCDF-library is installed. On computers running Linux, netCDF is usually available as a package of the operating system. On MacOS one can install the netcdf library e.g. using Homebrew or MacPorts. NetCDF is a self-describing binary output format, but here it is not required that you know details about it. Anyway, if you like to look 'into' a NetCDF file, please try to use ncdump FILENAME | less.

To run the forward model use

  ./lorenz_96 -total_steps 10000

This runs the Lorenz-96 model for 10000 time steps and the trajectory is written into a file state.nc.

2. Generate observations and a covariance matrix

To compile the tool programs use

  cd tools
  make all PDAF_ARCH=linux_gfortran

Now run

  ./generate_obs

to generate a file holding observations (obs.nc in models/lorenz96/) and

  ./generate_covar

to generate a file holding the covariance matrix information (covar.nc in models/lorenz96/), which is used to generate an initial ensemble for the data assimilation experiments. The covariance information consists of singular vectors (EOFs), corresponding singular values and the mean state over the 10000 time steps.

3. Build and run the assimilation program

Now edit the file make.arch/linux_gfortran.h again and change the line

CPP_DEFS = #-DUSE_PDAF

back to

CPP_DEFS = -DUSE_PDAF

to activate the coupling to PDAF in the model code.

Now compile the Lorenz-96 model with activated coupling to PDAF. First remove the binary files from the previous compilation using

  make clean PDAF_ARCH=linux_gfortran

Then build the assimilative model using

  make pdaf_lorenz_96 PDAF_ARCH=linux_gfortran

The program pdaf_lorenz_96 is generated.

To run the assimilation program, do

  ./tools/runasml.sh 

The script runasml.sh runs 11 experiments with a fixed ensemble size, but different covariance inflations (forgetting factors). The execution can take a few minutes.

4. Plot output from the assimilation experiments

To display the output of the assimilation experiments we provide plotting scripts for Python or Matlab, Octave. To use them do

cd plotting

The following plot functions are available:

Script	Description
plot_example	Plot the mean true and estimated RMS errors for the example output files generated with the script tools/runasml.sh as a function of the forgetting factor
plot_eofs	Plot the mean state and the eigenvectors (EOFs) from the file holdig the covariance matrix that is generated with tools/generate_covar.F90.
plot_obs	Plot the observations for a specified time step. The observation file is generated using tools/generate_obs.F90.
plot_rms	Plot the true and estimated RMS errors over time from an assimilation experiment.
plot_sigma	Plot the singular values from the file holdig the covariance matrix that is generated with tools/generate_covar.F90.
plot_state	Plot the state from the file holding a true state trajectory (state.nc) or the estimated forecast or analysis states from an assimilation experiment.

Except for plot_example all scripts require the specification of the directory and name of the file to be read. Sometimes, there are additional arguments like the time step index. The syntax of the functions is identical for Matlab and Python (except for plot_rms):

Python	Matlab/Octave	Comment
plot_example()	plot_example()
plot_eofs(filename, index)	plot_eofs(filename, index)	index=0 for mean state, index>1 for EOFs
plot_obs(filename, timestep)	plot_obs(filename, timestep)
plot_rms(filename [, plot_forecast=True, plot_analysis=True])	plot_rms(filename [, plot_forecast=1, plot_analysis=1])	plot_forecast/plot_analysis are optional and switch on/off a plot
plot_sigma(filename)	plot_sigma(filename)
plot_state(filename, timestep [, choice='t'])	plot_state(filename, timestep [, choice='t'])	choice is optional; options for choice: 't' true (use with truth file state.nc); for use with assimilation experiment output files: 'f' forecast, 'a' analysis, 'i' initial

Here filename is the name of the file including its path. In Matlab use 'help' to display the information about required input.

Python plotting
For Python the scripts are provided by plot_l96.py.

The module can either be imported, e.g. to use its functions interactively:
>>> import plot_l96
>>> plot_l96.plot_obs('../obs.nc', 4)

Alternatively the script can be run from the command line, providing the function name and its argument as command line parameters:
./plot_l96.py plot_obs ../obs.nc 4
(If this fails you can also try to run the script as python plot_l96.py plot_obs ../obs.nc 4).
To plot only the analysis RMS error one can use
./plot_l96.py plot_rms ../t1_N30_f0.95.nc False True

Matlab/Octave plotting examples:
plot_obs('../obs.nc',100) plots the observation at time step 100
plot_state('../t1_N30_f0.97.nc',100,'f') plots the forecast state estimate at the 100th analysis step
plot_rms('../t1_N30_f0.97.nc') plots the true and estimated RMS errors over time for the chosen experiment
plot_state('../state.nc',1101) plots the true state at model time step 1101 (= analysis step 100)

Run options

The implementation for the Lorenz-96 model has a wide range of options. For the full set, we recommend to check the list in the file init_pdaf.F90 in the directory lorenz96. Essentially all of the options can be specified on the command line in the format '-VARIABLE VALUE', where VARIABLE is the name of the variable in the program and VALUE is its value. The experiments in runasml.sh are executed like

./pdaf_lorenz_96 -total_steps 5000 -step_null 1000 -dim_ens 30 \
     -filtertype 6 -forget 0.99 -file_asml t1_N30_f0.99.nc

The meaning of the different options is the following:

Variable	Description	Default value
dim_ens	ensemble size	30
file_asml	the name of the output file	assimilation.nc
filtertype	Specifies the filter to use; 6 sets the filter ESTKF	1
forget	value of the forgetting factor for multiplicative variance inflation	1.0
step_null	initial time step of experiment. This specifies the offset of the observations which are generated from a model forward run. A value 1000 avoids the spin-up phase of the model	0
total_steps	number of time step in the experiment	5000

Some further options:

Variable	Description	Default value
delt_obs	time interval between observations, i.e. the forecast length	1
dim_lag	set the time lag for the smoother (the smoother is active for dim_lag>0)	0
dim_state	set state dimension (if changed, the true state trajectory, observations, and covariance matrix file need to be re-generated)	40
dx_obs	distance between observations on the grid; let's you specify an incompletely observated state (only used when use_obs_mask=.true.)	1
cradius	localization cut-off radius in grid points	5
locweight	choose localization weight function, e.g. 4 is the 5th-order polynomial mimicking a Gaussian (see Gaspari and Cohn 1999)	0
model_error	a logical variable activating model error noise	.false.
model_err_amp	amplitude of the model error	0.1
numobs	number of observed grid points; the points 1 to numobs are observed (can be combined with dx_obs; only used when use_obs_mask=.true.)	dim_state
use_obs_mask	whether to use incomplete observations (need to be .true. for any other settings on the observation density to be used)	.false.
type_ensinit	Choose whether the initial ensemble is sampled using 2nd-order exact smapling ('eof') or using random draws from the foreward trajectory ('rnd'). Note that 'eof' can only be used up to dim_ens=41 (state dimension + 1).	'eof'

The localization is controlled by the variables locweight, cradius, and sradius:

locweight	0	1	2	3	4
function	unit weight	exponential	5-th order polynomial	5-th order polynomial	5-th order polynomial
regulation	-	-	-	regulation using mean variance	regulation using variance of single observation point
cradius	weight=0 if distance > cradius
sradius	no impact	weight = exp(-d / sradius)	weight = 0 if d >= sradius else weight = f(sradius, distance)

Here, 'regulation' refers to the regulated localization introduced in Nerger, L., Janjić, T., Schröter, J., Hiller, W. (2012). A regulated localization scheme for ensemble-based Kalman filters. Quarterly Journal of the Royal Meteorological Society, 138, 802-812. ​​doi:10.1002/qj.945.

For the full set of options, please see init_pdaf.F90 and mod_assimilation.F90 in the lorenz96 directory. There, also the possible settings e.g. for filtertype are described.