Lorenz-63 model with PDAF

Note: The data assimilation with the Lorenz-63 model has been added in version 1.14 of PDAF. Here, the use is described for PDAF 1.15 where we moved the Lorenz-63 case into the directory models.

The implementation of the Lorenz-63 model coupled to PDAF is in the directory models/lorenz63/ of the PDAF package. Provided is a full implementation of PDAF with the nonlinear Lorenz-63 model (E. N. Lorenz (1963) Deterministic non-periodic flows. J. Atmos. Sci. 20, 130-141) providing various filter and smoother methods. The model has state dimenion 3. So it's not usable for localized filters, but it's a good test case for the particle filter.

Next to the implementation of the Lorenz-63 model with PDAF, the test case provides tool programs and scripts 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-63 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 files with observations and a covariance matrix for the initialization of the initial ensemble. In the second step, one compiles the Lorenz-63 model with activated coupling to PDAF and runs the data assimilation experiments.

1. Compile and run the forward model without assimilation

First change in the file make.arch/linux_gfortran.h the line

CPP_DEFS = -DUSE_PDAF

to

CPP_DEFS = #-DUSE_PDAF

to deactivate the coupling to PDAF in the model code. Now build the forward model program with

  cd models/lorenz63
  make lorenz_63 PDAF_ARCH=linux_gfortran

You have to ensure that in the machine-specific make include file linux_gfortran.h -DUSE_PDAF is not defined for CPP_DEFS (such that calls to PDAF are not active). You can replace linux_gfortran by any other make include file from make.arch/, e.g. specify macos_gfortran for compiling on MacOS. The executable is generated in the directory.

Note: The implementation uses the NetCDF library for file outputs. If the compilation above fails, please ensure the netcdf-library ist installed. On computers running Linux, it is usually available as a package of the operating system. On MacOS one can install the netcdf library e.g. using Fink 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_63 -total_steps 5000

This runs the Lorenz-63 model for 5000 time steps and the trajectory is written into a file state_l63.nc.

2. Generate observations

To build the executable to generate observations from the state trajectore use

  cd tools
  make generate_obs PDAF_ARCH=linux_gfortran

Now run

  ./generate_obs

to generate a file holding observations (obs_l63.nc in models/lorenz63/).

3. Build and run the assimilation program

Now do

  cd ../../../make.arch

and change in the file linux_gfortran.h 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-63 model with activated coupling to PDAF. First clean the directories for the main driver and the Lorenz-63 model using

  cd ../models/lorenz63
  make clean PDAF_ARCH=linux_gfortran

(This removes object files that were compiled without support for PDAF) Then build the executable using

  make pdaf_lorenz_63 PDAF_ARCH=linux_gfortran

The program pdaf_lorenz_63 is generated in the current directory.

To run the assimilation program, do

  ./tools/run_ESTKF.sh 

The script run_ESTKF.sh runs an experiment with the ESTKF filter method for ensmeble size 20 in which only the variable X os the model is observed at each 10th time step. The execution should only take seconds. Further you can run

  ./tools/run_PF.sh 

to run a similar experiment with the particle filter.

4. Plot output from the assimilation experiments

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

cd plotting/

The following plot functions are available:

Script	Description
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_obs_series	Plot a time series of the observations for a specified variable of the model. 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.
plot_state_series	Plot a time series of the state from the file holding a true state trajectory (state_l63.nc) or the estimated forecast or analysis states from an assimilation experiment for a selected state variable.

The scripts need the specification of the input file and the time step or variable to plot. The syntax of the functions is identical for Matlab and Python (except for plot_rms):

Matlab	Python	Comment
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_obs_series(filename, variable)	plot_obs_series(filename, variable)
plot_rms(filename [, plot_forecast=1, plot_analysis=1])	plot_rms(filename [, plot_forecast=True, plot_analysis=True])	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
plot_state_series(filename, variable [, choice='t'])	plot_state(filename, variable [, choice='t'])	choice is optional; options for choice: 't' true (use with truth file state_l63.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.

Matlab/Octave plotting examples:
plot_obs('../obs.nc',100) plots the observation at time step 100
plot_obs_series('../obs.nc',1) plots the time series of observations for variable X (likely set 2=Y, 3=Z)
plot_state('../ESTKF_N20.nc',100,'f') plots the forecast state estimate from the ESTKF run example 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_l63.nc',1101) plots the true state at model time step 1101 (= analysis step 100)
plot_state_series('../ESTKF_N20.nc',1,'f') plots the time series for the ESTKF run example of the forecast state for variable X

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

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

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

Run options

The implementation for the Lorenz-63 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 lorenz63. 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 experiment in run_ESTKF.sh is executed like

./pdaf_lorenz_63 -total_steps 5000 -step_null 1000 -dim_ens 20 \
        -forget 0.8 -filtertype 6 -file_asml ESTKF_N20.nc

and the experiment in run_PF.sh is executed like

./pdaf_lorenz_63 -total_steps 5000 -step_null 1000 -dim_ens 20 \
        -pf_res_type 2 -pf_noise_type 2 -pf_noise_amp 0.2 -filtertype 12 -file_asml PF_N20.nc

The meaning of the different options is the following:

Variable	Description	Default value
dim_ens	ensemble size	20
file_asml	the name of the output file	ESTKF_N20.nc
filtertype	Specifies the filter to use; 6 sets the ESTKF filter, 12 the particle filter	6
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
pf_res_type	type of resampling scheme in the particle filter (2 set stochastic unversal resampling)	1
pf_noise_type	type of noise to perturb resampled particles (2 sets noise with amplitude relative to ensemble standard deviation)	0
pf_noise_amp	amplitude of noise to perturb resampled particles	0.0

Some further options:

Variable	Description	Default value
delt_obs	time interval between observations, i.e. the forecast length	10
dim_lag	set the time lag for the smoother (the smoother is active for dim_lag>0)	0
dx_obs	distance between observations on the grid; let's you specify an incompletely observated state (only used when use_obs_mask=.true.)	3
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.)	1
use_obs_mask	whether to use incomplete observations (need to be .true. for any other settings on the observation density to be used)	.true.

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