Add Fortran version of MAOOAM that computes the Lyapunov Exponents.

By Sebastian Schubert and Jonathan Demaeyer.

README.md

@@ -1,4 +1,5 @@
 # Modular arbitrary-order ocean-atmosphere model: MAOOAM
+## Lyapunov exponents computation version
 
 ------------------------------------------------------------------------
 
@@ -40,6 +41,27 @@ exchanges.
 
 ------------------------------------------------------------------------
 
+## Note on the Lyapunov exponents computation
+
+This version of the code (Fortran version) allows for the computation of the
+backward Lyapunov exponents. The method used is described in
+
+* Benettin, G., Galgani, L., Giorgilli, A., and Strelcyn, J. M. : Lyapunov
+characteristic exponents for smooth dynamical systems; a method for computing
+all of them. Part 2: Numerical application. Meccanica, 15, 21-30,
+[doi:10.1007/BF02128237](http://dx.doi.org/10.1007/BF02128237), 1980.
+
+This version has been used to compute the Lyapunov exponents of MAOOAM in the
+following manuscript:
+
+* De Cruz, L., Schubert, S., Demaeyer, J., Lucarini, V., and Vannitsem, S.:
+Exploring the Lyapunov instability properties of high-dimensional atmospheric
+and climate models, Nonlin. Processes Geophys. Discuss.,
+[doi:10.5194/npg-2017-76](https://doi.org/10.5194/npg-2017-76), in
+review, 2018.
+
+------------------------------------------------------------------------
+
 ## Implementation notes ##
 
 As the system of differential equations is at most bilinear in y[j] (j=1..n), y

fortran/LICENSE.txt

@@ -1,7 +1,7 @@
 
 The MIT License (MIT)
 
-Copyright (c) 2015-2016 Lesley De Cruz and Jonathan Demaeyer
+Copyright (c) 2015-2016 Lesley De Cruz, Jonathan Demaeyer and Sebastian Schubert
 
 Permission is hereby granted, free of charge, to any person obtaining a copy
 of this software and associated documentation files (the "Software"), to deal

fortran/Makefile

@@ -68,9 +68,9 @@ endif
 .SUFFIXES:
 .SUFFIXES: .f90 .o .mod .out .test
 
-PROGRAMS = maooam test_tl_ad
+PROGRAMS = maooam_lyap maooam test_tl_ad maooam_lyap maooam_lyap_div tr_jacob_mat
 
-MODULES = util.mod params.mod inprod_analytic$(ST).mod aotensor_def.mod ic_def.mod $(INTEGRATOR).mod stat.mod tensor.mod tl_ad_tensor.mod $(TL_AD_INTEGRATOR).mod $(IFPORT)
+MODULES = util.mod params.mod inprod_analytic$(ST).mod aotensor_def.mod ic_def.mod $(INTEGRATOR).mod stat.mod tensor.mod tl_ad_tensor.mod $(TL_AD_INTEGRATOR).mod icdelta_def.mod lyap_vectors.mod lyap_stat.mod 
 MODULE_OBJECTS = $(MODULES:.mod=.o)
 
 TEST_PROGRAMS = test_inprod_analytic test_aotensor
@@ -109,6 +109,7 @@ $(INTEGRATOR).mod: aotensor_def.mod
 ic_def.mod: util.mod aotensor_def.mod
 stat.mod: params.mod
 tl_ad_tensor.mod: aotensor_def.mod
+lyap_vectors.mod: params.mod
 $(TL_AD_INTEGRATOR).mod: tl_ad_tensor.mod
 util.mod: $(IFPORT)
 

fortran/README.md

@@ -1,8 +1,11 @@
-# Modular arbitrary-order ocean-atmosphere model: MAOOAM -- Fortran implementation #
+# Modular arbitrary-order ocean-atmosphere model: MAOOAM -- Fortran implementation
+## Lyapunov exponents computation version
+
+------------------------------------------------------------------------
 
 ## About ##
 
-(c) 2013-2016 Lesley De Cruz and Jonathan Demaeyer
+(c) 2013-2016 Lesley De Cruz, Jonathan Demaeyer and Sebastian Schubert
 
 See [LICENSE.txt](../LICENSE.txt) for license information.
 
@@ -25,6 +28,27 @@ resources.
 
 ------------------------------------------------------------------------
 
+## Note on the Lyapunov exponents computation
+
+This version of the code allows for the computation of the backward Lyapunov
+exponents. The method used is described in
+
+* Benettin, G., Galgani, L., Giorgilli, A., and Strelcyn, J. M. : Lyapunov
+characteristic exponents for smooth dynamical systems; a method for computing
+all of them. Part 2: Numerical application. Meccanica, 15, 21-30,
+[doi:10.1007/BF02128237](http://dx.doi.org/10.1007/BF02128237), 1980.
+
+This version has been used to compute the Lyapunov exponents of MAOOAM in the
+following manuscript:
+
+* De Cruz, L., Schubert, S., Demaeyer, J., Lucarini, V., and Vannitsem, S.:
+Exploring the Lyapunov instability properties of high-dimensional atmospheric
+and climate models, Nonlin. Processes Geophys. Discuss.,
+[doi:10.5194/npg-2017-76](https://doi.org/10.5194/npg-2017-76), in
+review, 2018.
+
+------------------------------------------------------------------------
+
 ## Installation ##
 
 The program can be installed with Makefile. We provide configuration options for 
@@ -75,7 +99,10 @@ includes all the coefficients. This tensor is computed once at the program
 initialization.
 
 * maooam.f90 : Main program.
+* maooam_lyap.f90 : Version of the main program with the computation of the Lyapunov exponents.
+* maooam_lyap_div.f90 : Version of the main program with the computation of the first Lyapunov exponents with the divergence method.
 * aotensor_def.f90 : Tensor aotensor computation module.
+* lyap_vectors.f90 : Module for the computation of Lyapunov exponents and vectors
 * IC_def.f90 : A module which loads the user specified initial condition.
 * inprod_analytic.f90 : Inner products computation module.
 * rk2_integrator.f90 : A module which contains the Heun integrator for the model equations.
@@ -86,6 +113,7 @@ initialization.
 * rk2_tl_ad_integrator.f90 : Heun Tangent Linear (TL) and Adjoint (AD) model integrators module
 * rk4_tl_ad_integrator.f90 : RK4 Tangent Linear (TL) and Adjoint (AD) model integrators module
 * test_tl_ad.f90 : Tests for the Tangent Linear (TL) and Adjoint (AD) model versions
+* tr_jacob_mat.f90 : Test program to compute the trace of the Jacobian matrix. 
 * README.md : The present file.
 * LICENSE.txt : The license text of the program.
 * util.f90 : A module with various useful functions.
@@ -104,6 +132,8 @@ A documentation is available [here](./doc/html/index.html) (html) and [here](./d
 The user first has to fill the params.nml and int_params.nml namelist files according to their needs.
 Indeed, model and integration parameters can be specified respectively in the params.nml and int_params.nml namelist files. Some examples related to already published article are available in the params folder.
 
+The parameters related to the Lyapunov exponents computation are located in the int_params.nml namelist file.
+
 The modeselection.nml namelist can then be filled : 
 * NBOC and NBATM specify the number of blocks that will be used in respectively the ocean and
   the atmosphere. Each block corresponds to a given x and y wavenumber.
@@ -124,7 +154,7 @@ By default a RK2 scheme is selected.
 Finally, the IC.nml file specifying the initial condition should be defined. To
 obtain an example of this configuration file corresponding to the model you
 have previously defined, simply delete the current IC.nml file (if it exists)
-and run the program :
+and run the main program :
 
     ./maooam
 
@@ -141,6 +171,25 @@ The tangent linear and adjoint models of MAOOAM are provided in the
 tl_ad_tensor, rk2_tl_ad_integrator and rk4_tl_ad_integrator modules. It is
 documented [here](./doc/html/md_doc_tl_ad_doc.html).
 
+The computation of the Lyapunov exponents can be done by running the program :
+
+    ./maooam_lyap
+
+It will generate four files:
+ * evol_field.dat : the recorded time evolution of the variables.
+ * mean_field.dat : the mean field (the climatology)
+ * lyapunov_exponents.dat : the recorded time evolution of the Lyapunov exponents
+ * mean_lyapunov.dat : the mean Lyapunov exponents as well as their variance.
+
+Alternatively, the first Lyapunov exponent can be computed by the divergence method, 
+by running :
+
+   ./maooam_lyap_div
+
+It will generate three files :
+ * evol_field.dat : the recorded time evolution of the variables.
+ * mean_field.dat : the mean field (the climatology)
+ * lyapunov_exponents_div.dat : the recorded time evolution of the first Lyapunov exponent
 
 ------------------------------------------------------------------------