= Lorenz-63 model with PDAF = 

[[PageOutline(2-3,Contents of this page)]]

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 10000
}}}
This runs the Lorenz-63 model for 10000 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 cleandriver PDAF_ARCH=linux_gfortran
  make cleanlorenz_63 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._l63nc) 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:'''[[BR]]
`plot_obs('../obs.nc',100)` plots the observation at time step 100[[BR]]
`plot_obs_series('../obs.nc',1)` plots the time series of observations for variable X (likely set 2=Y, 3=Z) [[BR]]
`plot_state('../ESTKF_N20.nc',100,'f')` plots the forecast state estimate from the ESTKF run example at the 100th analysis step[[BR]]
`plot_rms('../t1_N30_f0.97.nc')` plots the true and estimated RMS errors over time for the chosen experiment[[BR]]
`plot_state('../state_l63.nc',1101)` plots the true state at model time step 1101 (= analysis step 100)[[BR]] 
`plot_state_series('../ESTKF_N20.nc',1,'f')` plots the time series for the ESTKF run example of the forecast state for variable X[[BR]] 

'''Python plotting'''[[BR]]
For Python the scripts are provided by `plot_l63.py`.

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

Alternatively the script can be run from the command line, providing the function
name and its argument as command line parameters:[[BR]]
`./plot_l63.py plot_obs ../obs_l63.nc 4`[[BR]]
(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[[BR]]
`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 resmpling 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.