= 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.