= Lorenz-96 model with PDAF = 

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

The implementation of the Lorenz-96 model coupled to PDAF is in the directory `testsuite/src/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 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 that
allow to run a test case and to display the outputs. Note, that this implementation runs without parallelization.

This description focuses on PDAF version 1.13.2 and later as we have revised the Lorenz-96 model case for version 1.13.2.


== 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 files with observations and a covariance matrix for the initialization of the initial ensemble. In the second step, 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 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 testsuite/src
  make lorenz_96 PDAF_ARCH=linux_gfortran
}}}
in the directory `testsuite/src/` of the PDAF package. 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 `testsuite/bin/`.

'''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
{{{
  cd ../bin
  ./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 build the executables for the tool programs use
{{{
  cd ../src/lorenz96/tools
  make all PDAF_ARCH=linux_gfortran
}}}

Now run
{{{
  ./generate_obs
}}}
and
{{{
  ./generate_covar
}}}
to generate a file holding observations (`obs.nc` in testsuite/bin/) and
a file holding the covariance matrix information (`covar.nc` in
testsuite/bin/), which is used to generate an initial ensemble for
the data assimilation experiments.


=== 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-96 model with activated
coupling to PDAF. First clean the directories for the main driver and the Lorenz-96 model using
{{{
  cd ../testsuite/src
  make cleandriver PDAF_ARCH=linux_gfortran
  make cleanlorenz_96 PDAF_ARCH=linux_gfortran
}}}
(This removes object files that were compiled without support for PDAF)
Then build the executable using
{{{
  make pdaf_lorenz_96 PDAF_ARCH=linux_gfortran
}}} 
The program `pdaf_lorenz_96` is generated in testsuite/bin.


To run the assimilation program, do
{{{
  cd ../bin
  ../src/lorenz96/tools/runasml.sh 
}}}
The script runsasml.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 Matlab, Octave and (from PDAF 1.13.2) Python. To use them do
{{{
cd ../src/lorenz96/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`):

||= Matlab =||= Python =||= 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=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 ||

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('../../../bin/obs.nc',100)` plots the observation at time step 100[[BR]]
`plot_state('../../../bin/t1_N30_f0.97.nc',100,'f')` plots the forecast state estimate at the 100th analysis step[[BR]]
`plot_rms('../../../bin/t1_N30_f0.97.nc')` plots the true and estimated RMS errors over time for the chosen experiment[[BR]]
`plot_state('../../../bin/state.nc',1101)` plots the true state at model time step 1101 (= analysis step 100)[[BR]] 

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

The module can either be imported, e.g. to use its functions interactively:[[BR]]
`>>> import plot_l96`[[BR]]
`>>> plot_l96.plot_obs('../../../bin/obs.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_l96.py plot_obs ../../../bin/obs.nc 4`[[BR]]
(If this fails you can also try to run the script as
`python plot_l96.py plot_obs ../../../bin/obs.nc 4`). To plot only the Analysis RMS error one can use[[BR]]
`plot_rms ../../../bin/t1_N30_f0.95.nc False True`


== 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 1 -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; 1 sets the SEIK filter || 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 ||
|| local_range || localization 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. ||

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