Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/CHANGES.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/CHANGES.rst	(revision 12165)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/CHANGES.rst	(revision 12165)
@@ -0,0 +1,7 @@
+*******
+Changes
+*******
+
+.. todo::
+
+   List the main additions of the the new release
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/CONTRIBUTING.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/CONTRIBUTING.rst	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/CONTRIBUTING.rst	(revision 12165)
@@ -3,6 +3,10 @@
 ************
 
+.. todo::
+
+
+
 .. contents::
-	:local:
+   :local:
 
 Sending feedbacks
@@ -17,7 +21,7 @@
 - You have a question: create a topic in the appropriate :forge:`discussion <forum>`
 - You would like to raise and issue: open a new ticket of the right type depending of its severity
-  
+
   - "Unavoidable" :forge:`newticket?type=Bug       <bug>`
-    
+
   - "Workable"    :forge:`newticket?type=Defect <defect>`
 
@@ -27,5 +31,5 @@
 ===============
 
-You have build a development relevant for NEMO shared reference: an addition of the source code,  
+You have build a development relevant for NEMO shared reference: an addition of the source code,
 a full fork of the reference, ...
 
@@ -34,5 +38,5 @@
 
 The proposals for developments to be included in the shared NEMO reference are first examined by NEMO Developers
-Committee / Scientific Advisory Board. 
+Committee / Scientific Advisory Board.
 The implementation of a new development requires some additionnal work from the intial developer.
 These tasks will need to be scheduled with NEMO System Team.
@@ -42,7 +46,7 @@
 ----
 
-You only would like to inform NEMO community about your developments. 
+You only would like to inform NEMO community about your developments.
 You can promote your work on NEMO forum gathering  the contributions fromof the community by creating
-a specific topic here :forge:`discussion/forum/5 <dedicated forum>` 
+a specific topic here :forge:`discussion/forum/5 <dedicated forum>`
 
 
@@ -55,5 +59,5 @@
   routines to the ticket, to highlight the proposed changes by adding to the ticket the output of ``svn diff``
   or ``svn patch`` from your working copy.
-  
+
 | Your development seems relevant for addition into the future release of NEMO shared reference.
   Implementing it into NEMO shared reference following the usual quality control will require some additionnal work
@@ -61,9 +65,9 @@
   your suggestion should be send as a proposed enhancement here :forge:`newticket?type=Enhancement <enhancement>`
   including description of the development, its implementation, and the existing validations.
-  
-  The proposed enhancement will be examined  by NEMO Developers Committee / Scientific Advisory Board. 
+
+  The proposed enhancement will be examined  by NEMO Developers Committee / Scientific Advisory Board.
   Once approved by the  Committee, the assicated development task can be scheduled in NEMO development work plan,
   and tasks distributed between you as initial developer and PI of this development action, and the NEMO System Team.
-  
+
   Once sucessful (meeting the usual quality control steps) this action will allow the merge of these developments with
   other developments of the year, building the future NEMO.
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/INSTALL.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/INSTALL.rst	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/INSTALL.rst	(revision 12165)
@@ -3,4 +3,8 @@
 *******************
 
+.. todo::
+
+
+
 .. contents::
    :local:
@@ -10,5 +14,5 @@
 
 | The NEMO source code is written in *Fortran 95* and
-  some of its prerequisite tools and libraries are already included in the ``./ext`` subdirectory.
+  some of its prerequisite tools and libraries are already included in the download.
 | It contains the AGRIF_ preprocessing program ``conv``; the FCM_ build system and
   the IOIPSL_ library for parts of the output.
@@ -23,5 +27,5 @@
 - *Fortran* compiler (``ifort``, ``gfortran``, ``pgfortran``, ...),
 - *Message Passing Interface (MPI)* implementation (e.g. |OpenMPI|_ or |MPICH|_).
-- |NetCDF|_ library with its underlying |HDF|_ 
+- |NetCDF|_ library with its underlying |HDF|_
 
 **NEMO, by default, takes advantage of some MPI features introduced into the MPI-3 standard.**
@@ -40,258 +44,4 @@
    This will limit MPI features to those defined within the MPI-2 standard
    (but will lose some performance benefits).
-
-Specifics for NetCDF and HDF
-----------------------------
-
-NetCDF and HDF versions from .
-However access to all the options available with the XIOS IO-server will require
-the parallel IO support of these libraries which can be unavailable.
-
-| **To satisfy these requirements, it is common to have to compile from source
-  in this order HDF (C library) then NetCDF (C and Fortran libraries)**
-| It is also necessary to compile these libraries with the same version of the MPI implementation that
-  both NEMO and XIOS (see below) are compiled and linked with.
-
-.. hint::
-
-   | It is difficult to define the options for the compilation as
-     they differ from one architecture to another according to
-     the hardware used and the software installed.
-   | The following is provided without any warranty
-
-   .. code-block:: console
-
-      $ ./configure [--{enable-fortran,disable-shared,enable-parallel}] ...
-
-   It is recommended to build the tests ``--enable-parallel-tests`` and run them with ``make check``
-
-Particular versions of these libraries may have their own restrictions.
-State the following requirements for netCDF-4 support:
-
-.. caution::
-
-   | When building NetCDF-C library versions older than 4.4.1, use only HDF5 1.8.x versions.
-   | Combining older NetCDF-C versions with newer HDF5 1.10 versions will create superblock 3 files
-     that are not readable by lots of older software.
-   
-Extract and install XIOS
-========================
-
-With the sole exception of running NEMO in mono-processor mode
-(in which case output options are limited to those supported by the ``IOIPSL`` library),
-diagnostic outputs from NEMO are handled by the third party ``XIOS`` library.
-This can be used in two different modes:
-
-- *attached* - Every NEMO process also acts as a XIOS server
-- *detached* - Every NEMO process runs as a XIOS client.
-  Output is collected and collated by external, stand-alone XIOS server processors.
-
-.. important::
-
-   In either case, XIOS needs to be compiled before NEMO,
-   since the libraries are needed to successfully create the NEMO executable.
-
-Instructions on how to obtain and install the software can be found on the :xios:`XIOS wiki<wiki>`.
-
-.. hint::
-
-   It is recommended to use XIOS version 2.5.
-   This version should be more stable (in terms of future code changes) than the XIOS trunk.
-   It is also the version used by the NEMO system team when testing all developments and new releases.
-   
-   This particular version has its own branch and can be checked out and downloaded with:
-
-   .. code:: console
-
-      $ svn co https://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/branchs/xios-2.5
-
-Download the NEMO source code
-=============================
-
-.. code:: console
-
-   $ svn co https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/trunk
-
-Description of directory tree
------------------------------
-
-+-----------+------------------------------------------------------------+
-| Folder    | Purpose                                                    |
-+===========+============================================================+
-| ``arch``  | Settings (per architecture-compiler pair)                  |
-+-----------+------------------------------------------------------------+
-| ``cfgs``  | :doc:`Reference configurations <configurations>`           |
-+-----------+------------------------------------------------------------+
-| ``doc``   | - ``latex``    : LaTex source code for ref. manuals        |
-|           | - ``namelists``: k start guide                             |
-|           | - ``rst``      : ReST files for quick start guide          |
-+-----------+------------------------------------------------------------+
-| ``ext``   | Dependencies included (``AGRIF``, ``FCM`` & ``IOIPSL``)    |
-+-----------+------------------------------------------------------------+
-| ``mk``    | Building  routines                                         |
-+-----------+------------------------------------------------------------+
-| ``src``   | Modelling routines                                         |
-|           |                                                            |
-|           | - ``ICE``: |SI3| for sea ice                               |
-|           | - ``NST``: AGRIF for embedded zooms                        |
-|           | - ``OCE``: |OPA| for ocean dynamics                        |
-|           | - ``TOP``: |TOP| for tracers                               |
-+-----------+------------------------------------------------------------+
-| ``tests`` | :doc:`Test cases <test_cases>` (unsupported)               |
-+-----------+------------------------------------------------------------+
-| ``tools`` | :doc:`Utilities <tools>` to [pre|post]process data         |
-+-----------+------------------------------------------------------------+
-
-Setup your architecture configuration file
-==========================================
-
-All compiler options in NEMO are controlled using files in
-``./arch/arch-'my_arch'.fcm`` where 'my_arch' is the name of the computing
-architecture.  It is recommended to copy and rename an configuration file from
-an architecture similar to your owns. You will need to set appropriate values
-for all of the variables in the file. In particular the FCM variables:
-``%NCDF_HOME``; ``%HDF5_HOME`` and ``%XIOS_HOME`` should be set to the
-installation directories used for XIOS installation.
-
-.. code-block:: sh
-
-        %NCDF_HOME           /opt/local
-        %HDF5_HOME           /opt/local
-        %XIOS_HOME           /Users/$( whoami )/xios-2.5
-        %OASIS_HOME          /not/defined
-
-Compile and create NEMO executable
-==================================
-
-The main script to compile and create executable is called makenemo and located in the CONFIG directory, it is used to identify the routines you need from the source code, to build the makefile and run it.
-As an example, compile GYRE with 'my_arch' to create a 'MY_GYRE' configuration:
-
-.. code-block:: sh
-
-	./makenemo –m 'my_arch' –r GYRE -n 'MY_GYRE'
-
-The image below shows the structure and some content of "MY_CONFIG" directory from the launching of the configuration creation (directories and fundamental files created by makenemo).
-
-+------------+----------------------------------------------------+
-| Folder     | Purpose                                            |
-+============+====================================================+
-| ``BLD``    |                                                    |
-+------------+----------------------------------------------------+
-| ``EXP00``  |                                                    |
-+------------+----------------------------------------------------+
-| ``EXPREF`` |                                                    |
-+------------+----------------------------------------------------+
-| ``MY_SRC`` |                                                    |
-+------------+----------------------------------------------------+
-| ``WORK``   |                                                    |
-+------------+----------------------------------------------------+
-
-Folder with the symbolic links to all unpreprocessed routines considered in the configuration
-Compilation folder (executables, headers files, libraries, preprocessed routines, flags, …)
-Computation folder for running the model (namelists, xml, executables and inputs-outputs)
-Folder intended to contain your customised routines (modified from initial ones or new entire routines)
-
-After successful execution of makenemo command, the executable called opa is created in the EXP00 directory (in the example above, the executable is created in CONFIG/MY_GYRE/EXP00).
-
-More makenemo options
----------------------
-
-``makenemo`` has several other options that can control which source files are selected and
-the operation of the build process itself.
-These are::
-
-   Optional:
-      -d  Set of new sub-components (space separated list from ./src directory)
-      -e  Path for alternative patch  location (default: 'MY_SRC' in configuration folder)
-      -h  Print this help
-      -j  Number of processes to compile (0: no build)
-      -n  Name for new configuration
-      -s  Path for alternative source location (default: 'src' root directory)
-      -t  Path for alternative build  location (default: 'BLD' in configuration folder)
-      -v  Level of verbosity ([0-3])
-
-These options can be useful for maintaining several code versions with only minor differences but
-they should be used sparingly.
-Note however the ``-j`` option which should be used more routinely to speed up the build process.
-For example:
-
-.. code-block:: sh
-
-        ./makenemo –m 'my_arch' –r GYRE -n 'MY_GYRE' -j 8
-
-which will compile up to 8 modules simultaneously.
-
-
-Default behaviour
------------------
-
-At the first use, you need the -m option to specify the architecture
-configuration file (compiler and its options, routines and libraries to
-include), then for next compilation, it is assumed you will be using the
-same compiler.  If the –n option is not specified the last compiled configuration
-will be used.
-
-Tools used during the process
------------------------------
-
-*   functions.sh : bash functions used by makenemo, for instance to create the WORK directory
-*   cfg.txt : text list of configurations and source directories
-*   bld.cfg : FCM rules to compile 
-
-Examples
---------
-
-.. code-block:: sh
-
-        echo "Example to install a new configuration MY_CONFIG";
-        echo "with OPA_SRC and LIM_SRC_2 ";
-        echo "makenemo -n MY_CONFIG -d \"OPA_SRC LIM_SRC_2\"";
-        echo "";
-        echo "Available configurations :"; cat ${CONFIG_DIR}/cfg.txt;
-        echo "";
-        echo "Available unsupported (external) configurations :"; cat ${CONFIG_DIR}/uspcfg.txt;
-        echo "";
-        echo "Example to remove bad configuration ";
-        echo "./makenemo -n MY_CONFIG clean_config";
-        echo "";
-        echo "Example to clean ";
-        echo "./makenemo clean";
-        echo "";
-        echo "Example to list the available keys of a CONFIG ";
-        echo "./makenemo list_key";
-        echo "";
-        echo "Example to add and remove keys";
-        echo "./makenemo add_key \"key_iomput key_mpp_mpi\" del_key \"key_agrif\" ";
-        echo "";
-        echo "Example to add and remove keys for a new configuration, and do not compile";
-        echo "./makenemo -n MY_CONFIG -j0 add_key \"key_iomput key_mpp_mpi\" del_key \"key_agrif\" ";
-
-Running the model
-=================
-
-Once makenemo has run successfully, the opa executable is available in ``CONFIG/MY_CONFIG/EXP00``
-For the reference configurations, the EXP00 folder also contains the initial input files (namelists, \*xml files for the IOs…). If the configuration also needs NetCDF input files, this should be downloaded here from the corresponding tar file, see Users/Reference Configurations
-
-.. code-block:: sh
-
-        cd 'MY_CONFIG'/EXP00
-        mpirun -n $NPROCS ./opa    # $NPROCS is the number of processes ; mpirun is your MPI wrapper
-
-
-Viewing and changing list of active CPP keys
-============================================
-
-For a given configuration (here called MY_CONFIG), the list of active CPP keys can be found in:
-
-.. code-block:: sh
-
-        ./cfgs/'MYCONFIG'/cpp_'MY_CONFIG'.fcm
-
-
-This text file can be edited to change the list of active CPP keys. Once changed, one needs to recompile opa executable using makenemo command in order for this change to be taken in account.
-Note that most NEMO configurations will need to specify the following CPP keys:
-``key_iomput`` and ``key_mpp_mpi``
-
-.. Links and substitutions
 
 .. |OpenMPI| replace:: *OpenMPI*
@@ -300,5 +50,221 @@
 .. _MPICH:   https://www.mpich.org
 .. |NetCDF|  replace:: *Network Common Data Form (NetCDF)*
-.. _NetCDF:  https://www.unidata.ucar.edu/downloads/netcdf
+.. _NetCDF:  https://www.unidata.ucar.edu
 .. |HDF|     replace:: *Hierarchical Data Form (HDF)*
-.. _HDF:     https://www.hdfgroup.org/downloads
+.. _HDF:     https://www.hdfgroup.org
+
+Specifics for NetCDF and HDF
+----------------------------
+
+NetCDF and HDF versions from official repositories may have not been compiled with MPI support.
+However access to all the options available with the XIOS IO-server will require
+the parallelism of these libraries.
+
+| **To satisfy these requirements, it is common to have to compile from source
+  in this order HDF (C library) then NetCDF (C and Fortran libraries)**
+| It is also necessary to compile these libraries with the same version of the MPI implementation that
+  both NEMO and XIOS (see below) have been compiled and linked with.
+
+.. hint::
+
+   | It is difficult to define the options for the compilation as
+     they differ from one architecture to another according to
+     the hardware used and the software installed.
+   | The following is provided without any warranty
+
+   .. code-block:: console
+
+      $ ./configure [--{enable-fortran,disable-shared,enable-parallel}] ...
+
+   It is recommended to build the tests ``--enable-parallel-tests`` and run them with ``make check``
+
+Particular versions of these libraries may have their own restrictions.
+State the following requirements for netCDF-4 support:
+
+.. caution::
+
+   | When building NetCDF-C library versions older than 4.4.1, use only HDF5 1.8.x versions.
+   | Combining older NetCDF-C versions with newer HDF5 1.10 versions will create superblock 3 files
+     that are not readable by lots of older software.
+
+Extract and install XIOS
+========================
+
+With the sole exception of running NEMO in mono-processor mode
+(in which case output options are limited to those supported by the ``IOIPSL`` library),
+diagnostic outputs from NEMO are handled by the third party ``XIOS`` library.
+It can be used in two different modes:
+
+:*attached*:  Every NEMO process also acts as a XIOS server
+:*detached*:  Every NEMO process runs as a XIOS client.
+  Output is collected and collated by external, stand-alone XIOS server processors.
+
+Instructions on how to install XIOS can be found on its :xios:`wiki<>`.
+
+.. hint::
+
+   It is recommended to use XIOS 2.5 release.
+   This version should be more stable (in terms of future code changes) than the XIOS trunk.
+   It is also the one used by the NEMO system team when testing all developments and new releases.
+
+   This particular version has its own branch and can be checked out with:
+
+   .. code:: console
+
+      $ svn co https://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/branchs/xios-2.5
+
+Download and install the NEMO code
+==================================
+
+Checkout the NEMO sources
+-------------------------
+
+.. code:: console
+
+   $ svn co https://forge.ipsl.jussieu.fr/nemo/svn/NEMO/trunk
+
+Description of 1\ :sup:`st` level tree structure
+------------------------------------------------
+
++---------------+----------------------------------------+
+| :file:`arch`  | Compilation settings                   |
++---------------+----------------------------------------+
+| :file:`cfgs`  | :doc:`Reference configurations <cfgs>` |
++---------------+----------------------------------------+
+| :file:`doc`   | :doc:`Documentation <doc>`             |
++---------------+----------------------------------------+
+| :file:`ext`   | Dependencies included                  |
+|               | (``AGRIF``, ``FCM`` & ``IOIPSL``)      |
++---------------+----------------------------------------+
+| :file:`mk`    | Compilation scripts                    |
++---------------+----------------------------------------+
+| :file:`src`   | :doc:`Modelling routines <src>`        |
++---------------+----------------------------------------+
+| :file:`tests` | :doc:`Test cases <tests>`              |
+|               | (unsupported)                          |
++---------------+----------------------------------------+
+| :file:`tools` | :doc:`Utilities <tools>`               |
+|               | to {pre,post}process data              |
++---------------+----------------------------------------+
+
+Setup your architecture configuration file
+------------------------------------------
+
+All compiler options in NEMO are controlled using files in :file:`./arch/arch-'my_arch'.fcm` where
+``my_arch`` is the name of the computing architecture
+(generally following the pattern ``HPCC-compiler`` or ``OS-compiler``).
+It is recommended to copy and rename an configuration file from an architecture similar to your owns.
+You will need to set appropriate values for all of the variables in the file.
+In particular the FCM variables:
+``%NCDF_HOME``; ``%HDF5_HOME`` and ``%XIOS_HOME`` should be set to
+the installation directories used for XIOS installation
+
+.. code-block:: sh
+
+   %NCDF_HOME    /usr/local/path/to/netcdf
+   %HDF5_HOME    /usr/local/path/to/hdf5
+   %XIOS_HOME    /home/$( whoami )/path/to/xios-2.5
+   %OASIS_HOME   /home/$( whoami )/path/to/oasis
+
+Create and compile a new configuration
+======================================
+
+The main script to {re}compile and create executable is called :file:`makenemo` located at
+the root of the working copy.
+It is used to identify the routines you need from the source code, to build the makefile and run it.
+As an example, compile a :file:`MY_GYRE` configuration from GYRE with 'my_arch':
+
+.. code-block:: sh
+
+   ./makenemo –m 'my_arch' –r GYRE -n 'MY_GYRE'
+
+Then at the end of the configuration compilation,
+:file:`MY_GYRE` directory will have the following structure.
+
++------------+----------------------------------------------------------------------------+
+| Directory  | Purpose                                                                    |
++============+============================================================================+
+| ``BLD``    | BuiLD folder: target executable, headers, libs, preprocessed routines, ... |
++------------+----------------------------------------------------------------------------+
+| ``EXP00``  | Run   folder: link to executable, namelists, ``*.xml`` and IOs             |
++------------+----------------------------------------------------------------------------+
+| ``EXPREF`` | Files under version control only for :doc:`official configurations <cfgs>` |
++------------+----------------------------------------------------------------------------+
+| ``MY_SRC`` | New routines or modified copies of NEMO sources                            |
++------------+----------------------------------------------------------------------------+
+| ``WORK``   | Links to all raw routines from :file:`./src` considered                    |
++------------+----------------------------------------------------------------------------+
+
+After successful execution of :file:`makenemo` command,
+the executable called `nemo` is available in the :file:`EXP00` directory
+
+More :file:`makenemo` options
+-----------------------------
+
+``makenemo`` has several other options that can control which source files are selected and
+the operation of the build process itself.
+
+.. literalinclude:: ../../../makenemo
+   :language: text
+   :lines: 119-143
+   :caption: Output of ``makenemo -h``
+
+These options can be useful for maintaining several code versions with only minor differences but
+they should be used sparingly.
+Note however the ``-j`` option which should be used more routinely to speed up the build process.
+For example:
+
+.. code-block:: sh
+
+        ./makenemo –m 'my_arch' –r GYRE -n 'MY_GYRE' -j 8
+
+will compile up to 8 processes simultaneously.
+
+Default behaviour
+-----------------
+
+At the first use,
+you need the ``-m`` option to specify the architecture configuration file
+(compiler and its options, routines and libraries to include),
+then for next compilation, it is assumed you will be using the same compiler.
+If the ``-n`` option is not specified the last compiled configuration will be used.
+
+Tools used during the process
+-----------------------------
+
+* :file:`functions.sh`: bash functions used by ``makenemo``, for instance to create the WORK directory
+* :file:`cfg.txt`     : text list of configurations and source directories
+* :file:`bld.cfg`     : FCM rules for compilation
+
+Examples
+--------
+
+.. literalinclude:: ../../../makenemo
+   :language: text
+   :lines: 146-153
+
+Running the model
+=================
+
+Once :file:`makenemo` has run successfully,
+the ``nemo`` executable is available in :file:`./cfgs/MY_CONFIG/EXP00`.
+For the reference configurations, the :file:`EXP00` folder also contains the initial input files
+(namelists, ``*.xml`` files for the IOs, ...).
+If the configuration needs other input files, they have to be placed here.
+
+.. code-block:: sh
+
+   cd 'MY_CONFIG'/EXP00
+   mpirun -n $NPROCS ./nemo   # $NPROCS is the number of processes
+                              # mpirun is your MPI wrapper
+
+Viewing and changing list of active CPP keys
+============================================
+
+For a given configuration (here called ``MY_CONFIG``),
+the list of active CPP keys can be found in :file:`./cfgs/'MYCONFIG'/cpp_MY_CONFIG.fcm`
+
+This text file can be edited by hand or with :file:`makenemo` to change the list of active CPP keys.
+Once changed, one needs to recompile ``nemo`` in order for this change to be taken in account.
+Note that most NEMO configurations will need to specify the following CPP keys:
+``key_iomput`` for IOs and ``key_mpp_mpi`` for parallelism.
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/README.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/README.rst	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/README.rst	(revision 12165)
@@ -1,7 +1,7 @@
-:Release:  |release|
-:Date:     |today|
-:SVN rev.: |revision|
+.. todo::
 
-NEMO_ for **Nucleus for European Modelling of the Ocean** is a state-of-the-art modelling framework for
+
+
+NEMO_ for *Nucleus for European Modelling of the Ocean* is a state-of-the-art modelling framework for
 research activities and forecasting services in ocean and climate sciences,
 developed in a sustainable way by a European consortium since 2008.
@@ -15,13 +15,14 @@
 The NEMO ocean model has 3 major components:
 
-- |OPA| models the ocean {thermo}dynamics and solves the primitive equations
-  (``./src/OCE``)
-- |SI3| simulates seaice {thermo}dynamics, brine inclusions and subgrid-scale thickness variations
-  (``./src/ICE``)
-- |TOP| models the {on,off}line oceanic tracers transport and biogeochemical processes
-  (``./src/TOP``)
+- |OCE| models the ocean {thermo}dynamics and solves the primitive equations
+  (:file:`./src/OCE`)
+- |ICE| simulates sea-ice {thermo}dynamics, brine inclusions and
+  subgrid-scale thickness variations (:file:`./src/ICE`)
+- |MBG| models the {on,off}line oceanic tracers transport and biogeochemical processes
+  (:file:`./src/TOP`)
 
-These physical core engines are described in their respective `references`_ that
-must be cited for any work related to their use.
+These physical core engines are described in
+their respective `reference publications <#project-documentation>`_ that
+must be cited for any work related to their use (see :doc:`cite`).
 
 Assets and solutions
@@ -33,16 +34,16 @@
 - Create :doc:`embedded zooms<zooms>` seamlessly thanks to 2-way nesting package AGRIF_.
 - Opportunity to integrate an :doc:`external biogeochemistry model<tracers>`
-- Versatile :doc:`data assimilation<data_assimilation>`
-- Generation of :doc:`diagnostics<diagnostics>` through effective XIOS_ system
-- Roll-out Earth system modeling with :doc:`coupling interface<coupling>` based on OASIS_
+- Versatile :doc:`data assimilation<da>`
+- Generation of :doc:`diagnostics<diags>` through effective XIOS_ system
+- Roll-out Earth system modeling with :doc:`coupling interface<cplg>` based on OASIS_
 
-Several :doc:`built-in configurations<configurations>` are provided to
+Several :doc:`built-in configurations<cfgs>` are provided to
 evaluate the skills and performances of the model which
-can be used as templates for setting up a new configurations (``./cfgs``).
+can be used as templates for setting up a new configurations (:file:`./cfgs`).
 
-The user can also checkout available :doc:`idealized test cases<test_cases>` that
-address specific physical processes(``./tests``).
+The user can also checkout available :doc:`idealized test cases<tests>` that
+address specific physical processes (:file:`./tests`).
 
-A set of :doc:`utilities <tools>` is also provided to {pre,post}process your data (``./tools``).
+A set of :doc:`utilities <tools>` is also provided to {pre,post}process your data (:file:`./tools`).
 
 Project documentation
@@ -50,23 +51,20 @@
 
 A walkthrough tutorial illustrates how to get code dependencies, compile and execute NEMO
-(``./INSTALL.rst``) . 
+(:file:`./INSTALL.rst`).
 
 Reference manuals and quick start guide can be build from source and
-exported to HTML or PDF formats (``./doc``) or
-downloaded directly from the :website:`website<bibliography/documentation>`.
+exported to HTML or PDF formats (:file:`./doc`) or
+downloaded directly from the :forge:`development platform<wiki/Documentations>`.
 
-=========== ===================== ===============
- Component   Reference Manual      Quick start
-=========== ===================== ===============
- |OPA|       |NEMO manual|_        |NEMO guide|
-             :cite:`NEMO_manual`
- |SI3|       |SI3 manual|
-             :cite:`SI3_manual`
- |TOP|       |TOP manual|
-             :cite:`TOP_manual`
-=========== ===================== ===============
+============ ================== ===================
+ Component    Reference Manual   Quick Start Guide
+============ ================== ===================
+ |NEMO-OCE|   |DOI man OCE|_     |DOI qsg|
+ |NEMO-ICE|   |DOI man ICE|
+ |NEMO-MBG|   |DOI man MBG|
+============ ================== ===================
 
 Since 2014 the project has a `Special Issue`_ in the open-access journal
-Geoscientific Model Development (GMD) from the European Geosciences Union (EGU).
+Geoscientific Model Development (GMD) from the European Geosciences Union (EGU_).
 The main scope is to collect relevant manuscripts covering various topics and
 to provide a single portal to assess the model potential and evolution.
@@ -79,22 +77,23 @@
 =================
 
-The NEMO Consortium pulling together 5 European institutes (CMCC_, CNRS_, MOI_, `Met Office`_ and NERC_)
-plans the sustainable development in order to keep a reliable evolving framework since 2008.
+The NEMO Consortium pulling together 5 European institutes
+(CMCC_, CNRS_, MOI_, `Met Office`_ and NERC_) plans the sustainable development in order to
+keep a reliable evolving framework since 2008.
 
-It defines the |NEMO strategy|_ that is implemented by the System Team on a yearly basis in order to
-release a new version almost every four years.
+It defines the |DOI dev stgy|_ that is implemented by the System Team on a yearly basis
+in order to release a new version almost every four years.
 
 When the need arises, :forge:`working groups<wiki/WorkingGroups>` are created or resumed to
 gather the community expertise for advising on the development activities.
 
+.. |DOI dev stgy| replace:: multi-year development strategy
 
-.. Substitutions / Links
+Disclaimer
+==========
 
-.. |NEMO manual| image:: https://zenodo.org/badge/DOI/10.5281/zenodo.1464816.svg
-.. |NEMO guide|  image:: https://zenodo.org/badge/DOI/10.5281/zenodo.1475325.svg
-.. |SI3 manual|  image:: https://zenodo.org/badge/DOI/10.5281/zenodo.1471689.svg
-.. |TOP manual|  image:: https://zenodo.org/badge/DOI/10.5281/zenodo.1471700.svg
+The NEMO source code is freely available and distributed under
+:download:`CeCILL v2.0 license <../../../LICENSE>` (GNU GPL compatible).
 
-.. |NEMO strategy| replace:: multi-year development strategy
-
-.. _Special Issue: https://www.geosci-model-dev.net/special_issue40.html
+You can use, modify and/or redistribute the software under its terms,
+but users are provided only with a limited warranty and the software's authors and
+the successive licensor's have only limited liability.
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/REFERENCES.bib
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/REFERENCES.bib	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/REFERENCES.bib	(revision 12165)
@@ -1,59 +1,46 @@
-@manual{NEMO_manual,
-   title={NEMO ocean engine},
-   author={Madec Gurvan and NEMO System Team},
-   organization={NEMO Consortium},
-   journal={Notes du Pôle de modélisation de l\'Institut Pierre-Simon Laplace (IPSL)},
-   number={27},
-   publisher={Zenodo},
-   abstract={The ocean engine of NEMO is a primitive equation model adapted to regional and
-   global ocean circulation problems.
-   It is intended to be a flexible tool for studying the ocean and its interactions with
-   the others components of the earth climate system over a wide range of space and time scales.},
-   doi={10.5281/zenodo.1464816},
-   edition={},
-   year={}
+@manual{NEMO_man,
+   title="NEMO ocean engine",
+   author="NEMO System Team",
+   series="Scientific Notes of Climate Modelling Center",
+   number="27",
+   institution="Institut Pierre-Simon Laplace (IPSL)",
+   publisher="Zenodo",
+   doi="10.5281/zenodo.1464816",
 }
+%   edition="",
+%   year=""
 
-@manual{SI3_manual,
-   title={SI³ – Sea Ice modelling Integrated Initiative – The NEMO Sea Ice engine},
-   author={NEMO Sea Ice Working Group},
-   organization={NEMO Consortium},
-   journal={Notes du Pôle de modélisation de l\'Institut Pierre-Simon Laplace (IPSL)},
-   number={31},
-   publisher={Zenodo},
-   abstract={SI³ (Sea Ice modelling Integrated Initiative) is the sea ice engine of NEMO
-   (Nucleus for European Modelling of the Ocean).
-   SI³ is based on the Arctic Ice Dynamics Joint EXperiment (AIDJEX) framework,
-   combining the ice thickness distribution framework, the conservation of horizontal momentum,
-   an elastic-viscous plastic rheology, and energy-conserving halo-thermodynamics.
-   SI³ is interfaced with the NEMO ocean engine, and, via the OASIS coupler, with
-   several atmospheric general circulation models.
-   It also supports two-way grid embedding via the AGRIF software.},
-   doi={10.5281/zenodo.1471689},
-   edition={},
-   year={}
+@manual{SI3_man,
+   title="Sea Ice modelling Integrated Initiative (SI$^3$) -- The NEMO Sea Ice engine",
+   author="NEMO Sea Ice Working Group",
+   series="Scientific Notes of Climate Modelling Center",
+   number="31",
+   institution="Institut Pierre-Simon Laplace (IPSL)",
+   publisher="Zenodo",
+   doi="10.5281/zenodo.1471689",
 }
+%   edition="",
+%   year=""
 
-@manual{TOP_manual,
-   title={TOP – Tracers in Ocean Paradigm – The NEMO Tracers engine},
-   author={NEMO TOP Working Group},
-   organization={NEMO Consortium},
-   journal={Notes du Pôle de modélisation de l\'Institut Pierre-Simon Laplace (IPSL)},
-   number={28},
-   publisher={Zenodo},
-   abstract={},
-   doi={10.5281/zenodo.1471700},
-   edition={},
-   year={}
+@manual{TOP_man,
+   title="Tracers in Ocean Paradigm (TOP) -- The NEMO Tracers engine",
+   author="NEMO TOP Working Group",
+   series="Scientific Notes of Climate Modelling Center",
+   number="28",
+   institution="Institut Pierre-Simon Laplace (IPSL)",
+   publisher="Zenodo",
+   doi="10.5281/zenodo.1471700",
 }
+%   edition="",
+%   year=""
 
-@Article{gmd-8-1245-2015,
-   author = {Vidard, A. and Bouttier, P.-A. and Vigilant, F.},
-   title = {{NEMOTAM}: {T}angent and {A}djoint {M}odels for the ocean modelling platform {NEMO}},
-   journal = {Geoscientific Model Development},
-   volume = {8},
-   year = {2015},
-   number = {4},
-   pages = {1245--1257},
-   doi = {10.5194/gmd-8-1245-2015}
+@article{TAM_pub,
+   author = "Vidard, A. and Bouttier, P.-A. and Vigilant, F.",
+   title = "NEMOTAM: Tangent and Adjoint Models for the ocean modelling platform NEMO",
+   journal = "Geoscientific Model Development",
+   volume = "8",
+   year = "2015",
+   number = "4",
+   pages = "1245--1257",
+   doi = "10.5194/gmd-8-1245-2015"
 }
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/RELEASE_NOTES.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/RELEASE_NOTES.rst	(revision 12164)
+++ 	(revision )
@@ -1,129 +1,0 @@
-*************
-Release Notes
-*************
-
-New sea-ice component SI3 (in place of LIMx)
-============================================
-
-- Physics
-
-  - Lateral melting
-  - Melt ponds: constant or :doi:`Holland et al. 2012<10.1175/JCLI-D-11-00078.1>` formulation (and soon topographic melt ponds)
-  - Ice-atm. drags from :doi:`Lupkes et al. 2012<10.1029/2012JD017630>` (depending on ice concentration) or :doi:`Lupkes et al. 2015<10.1002/2014JD022418>` (depending on sea ice concentration and atm. stability)
-  - Landfast ice (:doi:`Lemieux et al. 2016<10.1002/2016JC012006>`)
-
-- Numerics
-
-  - Advection: Ultimate-Macho scheme
-  - Rheology: adaptive EVP (:doi:`Kimmritz et al. 2017<10.1016/j.ocemod.2016.03.004>`)
-  - Coupling interface: conductivity as surface forcing instead of heat fluxes (Met Office requirement)
-
-- Performance: all thermodynamics in 1D and reduced mpp communications
-- Usability (users & developers friendly)
-
-  - Comprehensive set of outputs (universal units and understandable names + includes CMIP)
-  - New architecture and namelist
-  - All processes can be decoupled from each other (switch on/off)
-  - Ice categories bounds can be defined by the user or set automatically
-  - For open boundaries, the number of ice categories from the forcing model can be different from the number of categories in the regional simulation
-  - Fully compatible with AGRIF_
-  - Revised documentation 
-
-AGRIF for embedded zooms
-========================
-
-- Now compatible with new sea ice component
-- Now compatible with z* coordinate
-- Extended ghost cells area to properly handle scheme with spatial order >2
-- Added vertical refinement (beta)
-- Nesting tools for setup now up to date and working 
-
-Enhancements
-============
-
-- Fix for tracer conservation with split explicit free surface
-- Bulk formulae : move to aerobulk package (:doi:`Brodeau et al. 2016<10.1175/JPO-D-16-0169.1>`), i.e. NCAR, COARE and ECMWF bulk (remove Clio and MFS bulk)
-- Wetting and drying
-- Added tidal self attraction and loading either read from a file or from usual "scalar" approximation
-- Add a 4th order centered (CEN) and Flux Corrected Transport (FCT) tracer advection
-  (using a 4th compact in the vertical)
-- iso-neutral mixing (iso and triad operators):
-  add the Method of Stabilizing Correction (MSC) (more accurate calculation) + add a bilaplacian case
-- Lateral physics (``LDF``): scale aware setting of eddy viscosity and diffusivity
-- Vorticity: 2 new energy conserving scheme: ENT with Coriolis defined at T-point (better for Flux form) and EET a variant of EEN where e3t is used instead of e3f (solved the issue with e3f specification but is not enstrophy conserving) 
-
-Test Cases
-==========
-
-The first test cases available for now are:
-
-- ``CANAL``: east-west periodic canal of variable size with several initial states and associated geostrophic currents (zonal jets or vortex)
-- ``ICE_AGRIF``: east-west + north-south periodic channel. The common configuration includes an AGRIF zoom (1:3) in the middle of the basin to test how an ice patch is advected through it but one can also test the advection schemes (Prather and Ultimate-Macho) by removing the ``key_agrif`` in the cpp keys
-- ``ISOMIP``: simple box configuration with an ice shelf with simple geometry on top. The purpose of this test case is to evaluate the impact of various schemes and new development with iceshelf cavities. The exact original setup is described ​here
-- ``LOCK-EXCHANGE``: classical fluid dynamics experiment that has been adapted by Haidvogel and Beckmann (1999) for testing advection schemes in ocean circulation models. It has been used by several authors including Burchard and Bolding (2002) and :doi:`Ilıcak et al. (2012)<10.1016/j.ocemod.2011.10.003>`. The ``LOCK EXCHANGE`` experiment can in particular illustrate the impact of different choices of numerical schemes and/or subgrid closures on spurious interior mixing
-- ``OVERFLOW``: illustrates the impact of different choices of numerical schemes and/or subgrid closures on spurious interior mixing close to bottom topography. It is adapted from the non-rotating overflow configuration described in Haidvogel and Beckmann (1999) and further used by :doi:`Ilıcak et al. (2012)<10.1016/j.ocemod.2011.10.003>`
-- ``VORTEX``: illustrates the propagation of an anticyclonic eddy over a Beta plan and flat bottom. It is implemented here with an online refined subdomain (thanks to ``AGRIF`` library) out of which the vortex propagates. It serves as a benchmark to diagnose nesting errors as in :doi:`Debreu et al. (2012)<10.1016/j.ocemod.2012.03.003>`, :doi:`Penven et al. (2006)<10.1016/j.ocemod.2005.05.002>` and :doi:`Spall and Holland (1991)<10.1175/1520-0485(1991)021<0205:ANPEMF>2.0.CO;2>`
-- ``WAD``: a set of simple closed basin geometries for testing the wetting and drying capabilities. Examples range from a closed channel with EW linear bottom slope to a parabolic EW channel with a Gaussian ridge
-
-New Reference configurations
-============================
-
-``AGRIF_DEMO``: 2 interlocked zooms (1:4 & 1:3) in the Nordic Seas + 1 zoom (1:1) at the equator
-
-``SPITZ12``: regional configuration around the Svalbard archipelago 
-
-Wave coupling
-=============
-
-Coupled interface to external wave model
-
-Large scale wave interaction process added in momentum and tracer equations 
-
-Passive tracer TOP and biogeochemical PISCES components
-=======================================================
-
-- The passive tracers transport component was redesigned toward a modular structure and users can enable each module directly through logical flags in namelist_top (no more fortran macros!)
-- :doc:`TOP on-line user documentation<tracers>`
-- TOP currently accounts for the following 5 modules: ``CFC`` contains inorganic carbon tracers (CFC11/CFC12/SF6), ``MY_TRC`` is a template for new modules (or external couplings), ``AGE`` deals with water age tracking, ``C14`` as a radiocarbon passive tracer, and the companion ecosystem model ``PISCES``
-- A generalized infrastructure was developed to handle the prescription of either surface, coastal, or open boundaries conditions for each passive tracer
-- A new configuration, named ``ORCA2_OFF_TRC``, was created to provide a benchmark simulation environment to deal with inert carbon tracers dynamics by exploiting the offline coupling with NEMO
-- PISCES model contains new developments and modifications:
-
-  - Particulate Organic Carbon (POC) component comes with a new liability scheme, while the former Kriest parametrisation was superseded
-  - A complex iron chemistry scheme is now available, with an improved description of ligands for the marine iron cycle
-  - Carbonate chemistry is based on MOCSY 2.0 routines (see :doi:`10.5194/gmd-8-485-2015 <Orr and Epitalon, 2015>`), by complying also with CMIP6 standards
-  - Ecosystem components can be optionally modelled by means of explicit nutrient quotas (PISCES-QUOTA) 
-
-High Performance Computing (HPC): performances improvements
-===========================================================
-
-- Reduce number of MPI communications
-  (suppression of redundant communications, gather multiple communications into one)
-- Use of MPI-3 asynchronous routines for performance (use ``key_mpi2`` if MPI-3 not available)
-- Back to standard dynamical allocation (remove of wrk_alloc/dealloc statements)
-- :xios:`XIOS software<>` for IOs version 2.5 as default, and optionally available for restarts 
-
-Simplification and robustness
-=============================
-
-- Revised structure of ``namelist_ref`` / ``namelist_cfg`` and default reference values
-- Lateral physics (``LDF``): simplification of user interface and removal of CPP keys
-- Vertical physics (``ZDF``) (modularity, share shear production calculation between TKE and GKS, removal of all ZDF CPP keys, removal of ``avmu`` & ``avmv``, minimization of MPP comm.: ~15 removed)
-- Remove the split-explicit ZDF scheme for both ``TRA`` and ``DYN``
-- Remove the acceleration of convergence
-- Generalised ``lbc_lnk`` and ``lbc_nfd``
-- Unify mppini
-- Use non uniform jpi/jpj with dynamic allocation to avoid ghost rows/columns
-- MPI Message passing re coded
-- Configuration interface completely rewritten (``DOM`` module mainly suppressed , and in place: ``domain_cfg.nc`` file, or ``usr_def`` module) 
-
-Collaborative Development Environment
-=====================================
-
-- Access to information on wiki reorganised through portals for users/developers/System Team and complete refactoring of all pages and their layout
-- Reorganisation of SVN repository to be compliant with usual directory tree and facilitate building of the executable
-- Define and install a separate repository for test cases to all easy contributions from the community
-- :forge:`Forums<discussion>` created
-- :website:`Public web site<>` has been revamped and cleaned using Wordpress
-- New mailing lists have been set up
-- Improvements of reliability through automatic and regular testing of the changes made in repository 
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/AGRIF_DEMO/EXPREF/namelist_ice_cfg
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/AGRIF_DEMO/EXPREF/namelist_ice_cfg	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/AGRIF_DEMO/EXPREF/namelist_ice_cfg	(revision 12165)
@@ -38,4 +38,5 @@
 &namdyn_rhg     !   Ice rheology
 !------------------------------------------------------------------------------
+      ln_aEVP       = .false.          !     adaptive rheology (Kimmritz et al. 2016 & 2017)
 /
 !------------------------------------------------------------------------------
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/AGRIF_DEMO/README.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/AGRIF_DEMO/README.rst	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/AGRIF_DEMO/README.rst	(revision 12165)
@@ -2,4 +2,8 @@
 Embedded zooms
 **************
+
+.. todo::
+
+
 
 .. contents::
@@ -9,24 +13,26 @@
 ========
 
-AGRIF (Adaptive Grid Refinement In Fortran) is a library that allows the seamless space and time refinement over
-rectangular regions in NEMO.
+AGRIF (Adaptive Grid Refinement In Fortran) is a library that
+allows the seamless space and time refinement over rectangular regions in NEMO.
 Refinement factors can be odd or even (usually lower than 5 to maintain stability).
-Interaction between grid is "two-ways" in the sense that the parent grid feeds the child grid open boundaries and
-the child grid provides volume averages of prognostic variables once a given number of time step is completed.
+Interaction between grid is "two-ways" in the sense that
+the parent grid feeds the child grid open boundaries and
+the child grid provides volume averages of prognostic variables once
+a given number of time step is completed.
 These pages provide guidelines how to use AGRIF in NEMO.
-For a more technical description of the library itself, please refer to http://agrif.imag.fr.
+For a more technical description of the library itself, please refer to AGRIF_.
 
 Compilation
 ===========
 
-Activating AGRIF requires to append the cpp key ``key_agrif`` at compilation time: 
+Activating AGRIF requires to append the cpp key ``key_agrif`` at compilation time:
 
 .. code-block:: sh
 
-	./makenemo add_key 'key_agrif'
+   ./makenemo [...] add_key 'key_agrif'
 
-Although this is transparent to users, the way the code is processed during compilation is different from
-the standard case:
-a preprocessing stage (the so called "conv" program) translates the actual code so that
+Although this is transparent to users,
+the way the code is processed during compilation is different from the standard case:
+a preprocessing stage (the so called ``conv`` program) translates the actual code so that
 saved arrays may be switched in memory space from one domain to an other.
 
@@ -34,25 +40,26 @@
 ================================
 
-An additional text file ``AGRIF_FixedGrids.in`` is required at run time.
+An additional text file :file:`AGRIF_FixedGrids.in` is required at run time.
 This is where the grid hierarchy is defined.
-An example of such a file, here taken from the ``ICEDYN`` test case, is given below::
+An example of such a file, here taken from the ``ICEDYN`` test case, is given below
 
-	1
-	34 63 34 63 3 3 3
-	0
+.. literalinclude:: ../../../tests/ICE_AGRIF/EXPREF/AGRIF_FixedGrids.in
 
 The first line indicates the number of zooms (1).
 The second line contains the starting and ending indices in both directions on the root grid
-(imin=34 imax=63 jmin=34 jmax=63) followed by the space and time refinement factors (3 3 3).
+(``imin=34 imax=63 jmin=34 jmax=63``) followed by the space and time refinement factors (3 3 3).
 The last line is the number of child grid nested in the refined region (0).
 A more complex example with telescoping grids can be found below and
-in the ``AGRIF_DEMO`` reference configuration directory.
+in the :file:`AGRIF_DEMO` reference configuration directory.
 
-[Add some plots here with grid staggering and positioning ?]
+.. todo::
 
-When creating the nested domain, one must keep in mind that the child domain is shifted toward north-east and
-depends on the number of ghost cells as illustrated by the (attempted) drawing below for nbghostcells=1 and
-nbghostcells=3.
-The grid refinement is 3 and nxfin is the number of child grid points in i-direction.  
+   Add some plots here with grid staggering and positioning?
+
+When creating the nested domain, one must keep in mind that
+the child domain is shifted toward north-east and
+depends on the number of ghost cells as illustrated by
+the *attempted* drawing below for ``nbghostcells=1`` and ``nbghostcells=3``.
+The grid refinement is 3 and ``nxfin`` is the number of child grid points in i-direction.
 
 .. image:: _static/agrif_grid_position.jpg
@@ -62,14 +69,15 @@
 boundary data exchange and update being only performed between root and child grids.
 Use of east-west periodic or north-fold boundary conditions is not allowed in child grids either.
-Defining for instance a circumpolar zoom in a global model is therefore not possible. 
+Defining for instance a circumpolar zoom in a global model is therefore not possible.
 
 Preprocessing
 =============
 
-Knowing the refinement factors and area, a ``NESTING`` pre-processing tool may help to create needed input files
+Knowing the refinement factors and area,
+a ``NESTING`` pre-processing tool may help to create needed input files
 (mesh file, restart, climatological and forcing files).
 The key is to ensure volume matching near the child grid interface,
-a step done by invoking the ``Agrif_create_bathy.exe`` program.
-You may use the namelists provided in the ``NESTING`` directory as a guide.
+a step done by invoking the :file:`Agrif_create_bathy.exe` program.
+You may use the namelists provided in the :file:`NESTING` directory as a guide.
 These correspond to the namelists used to create ``AGRIF_DEMO`` inputs.
 
@@ -78,30 +86,21 @@
 
 Each child grid expects to read its own namelist so that different numerical choices can be made
-(these should be stored in the form ``1_namelist_cfg``, ``2_namelist_cfg``, etc... according to their rank in
-the grid hierarchy).
+(these should be stored in the form :file:`1_namelist_cfg`, :file:`2_namelist_cfg`, etc...
+according to their rank in the grid hierarchy).
 Consistent time steps and number of steps with the chosen time refinement have to be provided.
 Specific to AGRIF is the following block:
 
-.. code-block:: fortran
-
-	!-----------------------------------------------------------------------
-	&namagrif      !  AGRIF zoom                                            ("key_agrif")
-	!-----------------------------------------------------------------------
-	   ln_spc_dyn    = .true.  !  use 0 as special value for dynamics
-	   rn_sponge_tra = 2880.   !  coefficient for tracer   sponge layer [m2/s]
-	   rn_sponge_dyn = 2880.   !  coefficient for dynamics sponge layer [m2/s]
-	   ln_chk_bathy  = .false. !  =T  check the parent bathymetry
-	/             
+.. literalinclude:: ../../namelists/namagrif
+   :language: fortran
 
 where sponge layer coefficients have to be chosen according to the child grid mesh size.
 The sponge area is hard coded in NEMO and applies on the following grid points:
-2 x refinement factor (from i=1+nbghostcells+1 to i=1+nbghostcells+sponge_area) 
+2 x refinement factor (from ``i=1+nbghostcells+1`` to ``i=1+nbghostcells+sponge_area``)
 
-References
-==========
+.. rubric:: References
 
 .. bibliography:: zooms.bib
-	:all:
-	:style: unsrt
-	:labelprefix: A
-	:keyprefix: a-
+   :all:
+   :style: unsrt
+   :labelprefix: A
+   :keyprefix: a-
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_ice_cfg
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_ice_cfg	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_ice_cfg	(revision 12165)
@@ -38,4 +38,5 @@
 &namdyn_rhg     !   Ice rheology
 !------------------------------------------------------------------------------
+      ln_aEVP       = .false.          !     adaptive rheology (Kimmritz et al. 2016 & 2017)
 /
 !------------------------------------------------------------------------------
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/README.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/README.rst	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/README.rst	(revision 12165)
@@ -1,5 +1,9 @@
-************************
-Reference configurations
-************************
+********************************
+Run the Reference configurations
+********************************
+
+.. todo::
+
+   Lack of illustrations for ref. cfgs, and more generally in the guide.
 
 NEMO is distributed with a set of reference configurations allowing both
@@ -7,4 +11,8 @@
 the developer to test/validate his NEMO developments (using SETTE package).
 
+.. contents::
+   :local:
+   :depth: 1
+
 .. attention::
 
@@ -21,42 +29,49 @@
 ===========================================================
 
-A user who wants to compile the ORCA2_ICE_PISCES_ reference configuration using ``makenemo``
-should use the following, by selecting among available architecture file or providing a user defined one:
+To compile the ORCA2_ICE_PISCES_ reference configuration using :file:`makenemo`,
+one should use the following, by selecting among available architecture file or
+providing a user defined one:
 
 .. code-block:: console
-                
-   $ ./makenemo -r 'ORCA2_ICE_PISCES' -m 'my-fortran.fcm' -j '4'
+
+   $ ./makenemo -r 'ORCA2_ICE_PISCES' -m 'my_arch' -j '4'
 
 A new ``EXP00`` folder will be created within the selected reference configurations,
-namely ``./cfgs/ORCA2_ICE_PISCES/EXP00``,
-where it will be necessary to uncompress the Input & Forcing Files listed in the above table.
+namely ``./cfgs/ORCA2_ICE_PISCES/EXP00``.
+It will be necessary to uncompress the archives listed in the above table for
+the given reference configuration that includes input & forcing files.
 
 Then it will be possible to launch the execution of the model through a runscript
 (opportunely adapted to the user system).
-   
+
 List of Configurations
 ======================
 
-All forcing files listed below in the table are available from |NEMO archives URL|_
-
-.. |NEMO archives URL| image:: https://www.zenodo.org/badge/DOI/10.5281/zenodo.1472245.svg
-.. _NEMO archives URL: https://doi.org/10.5281/zenodo.1472245
-
-====================== ===== ===== ===== ======== ======= ================================================
- Configuration                     Component(s)                            Input & Forcing File(s)
----------------------- ---------------------------------- ------------------------------------------------
- Name                   OPA   SI3   TOP   PISCES   AGRIF
-====================== ===== ===== ===== ======== ======= ================================================
- AGRIF_DEMO_             X     X                     X     AGRIF_DEMO_v4.0.tar, ORCA2_ICE_v4.0.tar
- AMM12_                  X                                 AMM12_v4.0.tar
- C1D_PAPA_               X                                 INPUTS_C1D_PAPA_v4.0.tar
- GYRE_BFM_               X           X                     *none*
- GYRE_PISCES_            X           X      X              *none*
- ORCA2_ICE_PISCES_       X     X     X      X              ORCA2_ICE_v4.0.tar, INPUTS_PISCES_v4.0.tar
- ORCA2_OFF_PISCES_                   X      X              ORCA2_OFF_v4.0.tar, INPUTS_PISCES_v4.0.tar
- ORCA2_OFF_TRC_                      X                     ORCA2_OFF_v4.0.tar
- ORCA2_SAS_ICE_                X                           ORCA2_ICE_v4.0.tar, INPUTS_SAS_v4.0.tar
- SPITZ12_                X     X                           SPITZ12_v4.0.tar
-====================== ===== ===== ===== ======== ======= ================================================
+All forcing files listed below in the table are available from |DOI data|_
+
+=================== === === === === === ==================================
+ Configuration       Component(s)        Archives (input & forcing files)
+------------------- ------------------- ----------------------------------
+ Name                O   S   T   P   A
+=================== === === === === === ==================================
+ AGRIF_DEMO_         X   X           X   AGRIF_DEMO_v4.0.tar,
+                                         ORCA2_ICE_v4.0.tar
+ AMM12_              X                   AMM12_v4.0.tar
+ C1D_PAPA_           X                   INPUTS_C1D_PAPA_v4.0.tar
+ GYRE_BFM_           X       X           *none*
+ GYRE_PISCES_        X       X   X       *none*
+ ORCA2_ICE_PISCES_   X   X   X   X       ORCA2_ICE_v4.0.tar,
+                                         INPUTS_PISCES_v4.0.tar
+ ORCA2_OFF_PISCES_           X   X       ORCA2_OFF_v4.0.tar,
+                                         INPUTS_PISCES_v4.0.tar
+ ORCA2_OFF_TRC_              X           ORCA2_OFF_v4.0.tar
+ ORCA2_SAS_ICE_          X               ORCA2_ICE_v4.0.tar,
+                                         INPUTS_SAS_v4.0.tar
+ SPITZ12_            X   X               SPITZ12_v4.0.tar
+=================== === === === === === ==================================
+
+.. admonition:: Legend for component combination
+
+   O for OCE, S for SI\ :sup:`3`, T for TOP, P for PISCES and A for AGRIF
 
 AGRIF_DEMO
@@ -72,9 +87,14 @@
 particular interest to test sea ice coupling.
 
+.. image:: _static/AGRIF_DEMO_no_cap.jpg
+   :scale: 66%
+   :align: center
+
 The 1:1 grid can be used alone as a benchmark to check that
-the model solution is not corrupted by grid exchanges. 
+the model solution is not corrupted by grid exchanges.
 Note that since grids interact only at the baroclinic time level,
 numerically exact results can not be achieved in the 1:1 case.
-Perfect reproducibility is obtained only by switching to a fully explicit setup instead of a split explicit free surface scheme.
+Perfect reproducibility is obtained only by switching to a fully explicit setup instead of
+a split explicit free surface scheme.
 
 AMM12
@@ -85,5 +105,8 @@
 a regular horizontal grid of ~12 km of resolution (see :cite:`ODEA2012`).
 
-This configuration allows to tests several features of NEMO specifically addressed to the shelf seas. 
+.. image:: _static/AMM_domain.png
+   :align: center
+
+This configuration allows to tests several features of NEMO specifically addressed to the shelf seas.
 In particular, ``AMM12`` accounts for vertical s-coordinates system, GLS turbulence scheme,
 tidal lateral boundary conditions using a flather scheme (see more in ``BDY``).
@@ -99,21 +122,33 @@
 --------
 
-``C1D_PAPA`` is a 1D configuration for the `PAPA station <http://www.pmel.noaa.gov/OCS/Papa/index-Papa.shtml>`_ located in the northern-eastern Pacific Ocean at 50.1°N, 144.9°W.
-See `Reffray et al. (2015) <http://www.geosci-model-dev.net/8/69/2015>`_ for the description of its physical and numerical turbulent-mixing behaviour.
-
-The water column setup, called NEMO1D, is activated with the inclusion of the CPP key ``key_c1d`` and
-has a horizontal domain of 3x3 grid points.
-
-This reference configuration uses 75 vertical levels grid (1m at the surface), GLS turbulence scheme with K-epsilon closure and the NCAR bulk formulae.
+.. figure:: _static/Papa2015.jpg
+   :height: 225px
+   :align:  left
+
+``C1D_PAPA`` is a 1D configuration for the `PAPA station`_ located in
+the northern-eastern Pacific Ocean at 50.1°N, 144.9°W.
+See :gmd:`Reffray et al. (2015) <8/69/2015>` for the description of
+its physical and numerical turbulent-mixing behaviour.
+
+| The water column setup, called NEMO1D, is activated with
+  the inclusion of the CPP key ``key_c1d`` and
+  has a horizontal domain of 3x3 grid points.
+| This reference configuration uses 75 vertical levels grid (1m at the surface),
+  GLS turbulence scheme with K-epsilon closure and the NCAR bulk formulae.
+
 Data provided with ``INPUTS_C1D_PAPA_v4.0.tar`` file account for:
 
-- ``forcing_PAPASTATION_1h_y201[0-1].nc`` : ECMWF operational analysis atmospheric forcing rescaled to 1h (with long and short waves flux correction) for years 2010 and 2011
-- ``init_PAPASTATION_m06d15.nc`` : Initial Conditions from observed data and Levitus 2009 climatology
-- ``chlorophyll_PAPASTATION.nc`` : surface chlorophyll file from Seawifs data
+- :file:`forcing_PAPASTATION_1h_y201[0-1].nc`:
+  ECMWF operational analysis atmospheric forcing rescaled to 1h
+  (with long and short waves flux correction) for years 2010 and 2011
+- :file:`init_PAPASTATION_m06d15.nc`: Initial Conditions from
+  observed data and Levitus 2009 climatology
+- :file:`chlorophyll_PAPASTATION.nc`: surface chlorophyll file from Seawifs data
 
 GYRE_BFM
 --------
 
-``GYRE_BFM`` shares the same physical setup of GYRE_PISCES_, but NEMO is coupled with the `BFM <http://www.bfm-community.eu/>`_ biogeochemical model as described in ``./cfgs/GYRE_BFM/README``.
+``GYRE_BFM`` shares the same physical setup of GYRE_PISCES_,
+but NEMO is coupled with the `BFM`_ biogeochemical model as described in ``./cfgs/GYRE_BFM/README``.
 
 GYRE_PISCES
@@ -123,43 +158,42 @@
 in the Beta-plane approximation with a regular 1° horizontal resolution and 31 vertical levels,
 with PISCES BGC model :cite:`gmd-8-2465-2015`.
-Analytical forcing for heat, freshwater and wind-stress fields are applied.  
-
-This configuration acts also as demonstrator of the **user defined setup** (``ln_read_cfg = .false.``) and
-grid setting are handled through the ``&namusr_def`` controls in ``namelist_cfg``:
+Analytical forcing for heat, freshwater and wind-stress fields are applied.
+
+This configuration acts also as demonstrator of the **user defined setup**
+(``ln_read_cfg = .false.``) and grid setting are handled through
+the ``&namusr_def`` controls in :file:`namelist_cfg`:
 
 .. literalinclude:: ../../../cfgs/GYRE_PISCES/EXPREF/namelist_cfg
    :language: fortran
-   :lines: 34-42
+   :lines:    35-41
 
 Note that, the default grid size is 30x20 grid points (with ``nn_GYRE = 1``) and
 vertical levels are set by ``jpkglo``.
-The specific code changes can be inspected in ``./src/OCE/USR``.
-
-**Running GYRE as a benchmark** :
-this simple configuration can be used as a benchmark since it is easy to increase resolution,
-with the drawback of getting results that have a very limited physical meaning.
-
-GYRE grid resolution can be increased at runtime by setting a different value of ``nn_GYRE`` (integer multiplier scaling factor), as described in the following table: 
-
-=========== ========= ========== ============ ===================
-``nn_GYRE``  *jpiglo*  *jpjglo*   ``jpkglo``   **Equivalent to**
-=========== ========= ========== ============ ===================
- 1           30        20         31           GYRE 1°
- 25          750       500        101          ORCA 1/2°
- 50          1500      1000       101          ORCA 1/4°
- 150         4500      3000       101          ORCA 1/12°
- 200         6000      4000       101          ORCA 1/16°
-=========== ========= ========== ============ ===================
-
-Note that, it is necessary to set ``ln_bench = .true.`` in ``namusr_def`` to
-avoid problems in the physics computation and that
-the model timestep should be adequately rescaled. 
-
-For example if ``nn_GYRE = 150``, equivalent to an ORCA 1/12° grid,
-the timestep ``rn_rdt = 1200`` should be set to 1200 seconds
-
-Differently from previous versions of NEMO,
-the code uses by default the time-splitting scheme and
-internally computes the number of sub-steps. 
+The specific code changes can be inspected in :file:`./src/OCE/USR`.
+
+.. rubric:: Running GYRE as a benchmark
+
+| This simple configuration can be used as a benchmark since it is easy to increase resolution,
+  with the drawback of getting results that have a very limited physical meaning.
+| GYRE grid resolution can be increased at runtime by setting a different value of ``nn_GYRE``
+  (integer multiplier scaling factor), as described in the following table:
+
+=========== ============ ============ ============ ===============
+``nn_GYRE``  ``jpiglo``   ``jpjglo``   ``jpkglo``   Equivalent to
+=========== ============ ============ ============ ===============
+ 1           30           20           31           GYRE 1°
+ 25          750          500          101          ORCA 1/2°
+ 50          1500         1000         101          ORCA 1/4°
+ 150         4500         3000         101          ORCA 1/12°
+ 200         6000         4000         101          ORCA 1/16°
+=========== ============ ============ ============ ===============
+
+| Note that, it is necessary to set ``ln_bench = .true.`` in ``&namusr_def`` to
+  avoid problems in the physics computation and that
+  the model timestep should be adequately rescaled.
+| For example if ``nn_GYRE = 150``, equivalent to an ORCA 1/12° grid,
+  the timestep ``rn_rdt`` should be set to 1200 seconds
+  Differently from previous versions of NEMO, the code uses by default the time-splitting scheme and
+  internally computes the number of sub-steps.
 
 ORCA2_ICE_PISCES
@@ -174,41 +208,37 @@
 the ratio of anisotropy is nearly one everywhere
 
-this configuration uses the three components 
-
-- |OPA|, the ocean dynamical core 
-- |SI3|, the thermodynamic-dynamic sea ice model.
-- |TOP|, passive tracer transport module and PISCES BGC model :cite:`gmd-8-2465-2015`
+This configuration uses the three components
+
+- |OCE|, the ocean dynamical core
+- |ICE|, the thermodynamic-dynamic sea ice model.
+- |MBG|, passive tracer transport module and PISCES BGC model :cite:`gmd-8-2465-2015`
 
 All components share the same grid.
-
 The model is forced with CORE-II normal year atmospheric forcing and
 it uses the NCAR bulk formulae.
 
-In this ``ORCA2_ICE_PISCES`` configuration,
-AGRIF nesting can be activated that includes a nested grid in the Agulhas region.
-
-To set up this configuration, after extracting NEMO:
-
-Build your AGRIF configuration directory from ``ORCA2_ICE_PISCES``,
-with the ``key_agrif`` CPP key activated:
-
-.. code-block:: console
-                
-        $ ./makenemo -r 'ORCA2_ICE_PISCES' -n 'AGRIF' add_key 'key_agrif'
-
-By using the input files and namelists for ``ORCA2_ICE_PISCES``,
-the AGRIF test configuration is ready to run.
-
-**Ocean Physics**
-
-- *horizontal diffusion on momentum*: the eddy viscosity coefficient depends on the geographical position. It is taken as 40000 m^2/s, reduced in the equator regions (2000 m^2/s) excepted near the western boundaries.
-- *isopycnal diffusion on tracers*: the diffusion acts along the isopycnal surfaces (neutral surface) with an eddy diffusivity coefficient of 2000 m^2/s.
-- *Eddy induced velocity parametrization* with a coefficient that depends on the growth rate of baroclinic instabilities (it usually varies from 15 m^2/s to 3000 m^2/s).
-- *lateral boundary conditions* : zero fluxes of heat and salt and no-slip conditions are applied through lateral solid boundaries.
-- *bottom boundary condition* : zero fluxes of heat and salt are applied through the ocean bottom.
-  The Beckmann [19XX] simple bottom boundary layer parameterization is applied along continental slopes.
-  A linear friction is applied on momentum.
-- *convection*: the vertical eddy viscosity and diffusivity coefficients are increased to 1 m^2/s in case of static instability.
-- *time step* is 5760sec (1h36') so that there is 15 time steps in one day.
+.. rubric:: Ocean Physics
+
+:horizontal diffusion on momentum:
+   the eddy viscosity coefficient depends on the geographical position.
+   It is taken as 40000 m\ :sup:`2`/s, reduced in the equator regions (2000 m\ :sup:`2`/s)
+   excepted near the western boundaries.
+:isopycnal diffusion on tracers:
+   the diffusion acts along the isopycnal surfaces (neutral surface) with
+   an eddy diffusivity coefficient of 2000 m\ :sup:`2`/s.
+:Eddy induced velocity parametrization:
+   With a coefficient that depends on the growth rate of baroclinic instabilities
+   (it usually varies from 15 m\ :sup:`2`/s to 3000 m\ :sup:`2`/s).
+:lateral boundary conditions:
+   Zero fluxes of heat and salt and no-slip conditions are applied through lateral solid boundaries.
+:bottom boundary condition:
+   Zero fluxes of heat and salt are applied through the ocean bottom.
+   The Beckmann [19XX] simple bottom boundary layer parameterization is applied along
+   continental slopes.
+   A linear friction is applied on momentum.
+:convection:
+   The vertical eddy viscosity and diffusivity coefficients are increased to 1 m\ :sup:`2`/s in
+   case of static instability.
+:time step: is 5760sec (1h36') so that there is 15 time steps in one day.
 
 ORCA2_OFF_PISCES
@@ -218,34 +248,34 @@
 but only PISCES model is an active component of TOP.
 
-
 ORCA2_OFF_TRC
 -------------
 
-``ORCA2_OFF_TRC`` is based on the ORCA2 global ocean configuration
-(see ORCA2_ICE_PISCES_ for general description) along with the tracer passive transport module (TOP), but dynamical fields are pre-calculated and read with specific time frequency.
-
-This enables for an offline coupling of TOP components,
-here specifically inorganic carbon compounds (cfc11, cfc12, sf6, c14) and water age module (age).
-See ``namelist_top_cfg`` to inspect the selection of each component with the dedicated logical keys.
+| ``ORCA2_OFF_TRC`` is based on the ORCA2 global ocean configuration
+  (see ORCA2_ICE_PISCES_ for general description) along with
+  the tracer passive transport module (TOP),
+  but dynamical fields are pre-calculated and read with specific time frequency.
+| This enables for an offline coupling of TOP components,
+  here specifically inorganic carbon compounds (CFC11, CFC12, SF6, C14) and water age module (age).
+  See :file:`namelist_top_cfg` to inspect the selection of
+  each component with the dedicated logical keys.
 
 Pre-calculated dynamical fields are provided to NEMO using
-the namelist ``&namdta_dyn``  in ``namelist_cfg``,
+the namelist ``&namdta_dyn``  in :file:`namelist_cfg`,
 in this case with a 5 days frequency (120 hours):
 
-.. literalinclude:: ../../../cfgs/GYRE_PISCES/EXPREF/namelist_ref
+.. literalinclude:: ../../namelists/namdta_dyn
    :language: fortran
-   :lines: 935-960
-
-Input dynamical fields for this configuration (``ORCA2_OFF_v4.0.tar``) comes from
+
+Input dynamical fields for this configuration (:file:`ORCA2_OFF_v4.0.tar`) comes from
 a 2000 years long climatological simulation of ORCA2_ICE using ERA40 atmospheric forcing.
 
-Note that, this configuration default uses linear free surface (``ln_linssh = .true.``) assuming that
-model mesh is not varying in time and
-it includes the bottom boundary layer parameterization (``ln_trabbl = .true.``) that
-requires the provision of bbl coefficients through ``sn_ubl`` and ``sn_vbl`` fields.
-
-It is also possible to activate PISCES model (see ``ORCA2_OFF_PISCES``) or
-a user defined set of tracers and source-sink terms with ``ln_my_trc = .true.``
-(and adaptation of ``./src/TOP/MY_TRC`` routines).
+| Note that,
+  this configuration default uses linear free surface (``ln_linssh = .true.``) assuming that
+  model mesh is not varying in time and
+  it includes the bottom boundary layer parameterization (``ln_trabbl = .true.``) that
+  requires the provision of BBL coefficients through ``sn_ubl`` and ``sn_vbl`` fields.
+| It is also possible to activate PISCES model (see ``ORCA2_OFF_PISCES``) or
+  a user defined set of tracers and source-sink terms with ``ln_my_trc = .true.``
+  (and adaptation of ``./src/TOP/MY_TRC`` routines).
 
 In addition, the offline module (OFF) allows for the provision of further fields:
@@ -254,34 +284,32 @@
    by including an input datastream similarly to the following:
 
-.. code-block:: fortran
-
-   sn_rnf  = 'dyna_grid_T', 120, 'sorunoff' , .true., .true., 'yearly', '', '', ''
-
-2. **VVL dynamical fields**,
-   in the case input data were produced by a dyamical core using variable volume (``ln_linssh = .false.``)
-   it necessary to provide also diverce and E-P at before timestep by
+   .. code-block:: fortran
+
+      sn_rnf  = 'dyna_grid_T', 120, 'sorunoff' , .true., .true., 'yearly', '', '', ''
+
+2. **VVL dynamical fields**, in the case input data were produced by a dyamical core using
+   variable volume (``ln_linssh = .false.``)
+   it is necessary to provide also diverce and E-P at before timestep by
    including input datastreams similarly to the following
 
-.. code-block:: fortran
-
-   sn_div  = 'dyna_grid_T', 120, 'e3t'      , .true., .true., 'yearly', '', '', ''
-   sn_empb = 'dyna_grid_T', 120, 'sowaflupb', .true., .true., 'yearly', '', '', ''
-
+   .. code-block:: fortran
+
+      sn_div  = 'dyna_grid_T', 120, 'e3t'      , .true., .true., 'yearly', '', '', ''
+      sn_empb = 'dyna_grid_T', 120, 'sowaflupb', .true., .true., 'yearly', '', '', ''
 
 More details can be found by inspecting the offline data manager in
-the routine ``./src/OFF/dtadyn.F90``.
+the routine :file:`./src/OFF/dtadyn.F90`.
 
 ORCA2_SAS_ICE
 -------------
 
-ORCA2_SAS_ICE is a demonstrator of the Stand-Alone Surface (SAS) module and
-it relies on ORCA2 global ocean configuration (see ORCA2_ICE_PISCES_ for general description).
-
-The standalone surface module allows surface elements such as sea-ice, iceberg drift, and
-surface fluxes to be run using prescribed model state fields.
-It can profitably be used to compare different bulk formulae or
-adjust the parameters of a given bulk formula.
-
-More informations about SAS can be found in NEMO manual.
+| ORCA2_SAS_ICE is a demonstrator of the Stand-Alone Surface (SAS) module and
+  it relies on ORCA2 global ocean configuration (see ORCA2_ICE_PISCES_ for general description).
+| The standalone surface module allows surface elements such as sea-ice, iceberg drift, and
+  surface fluxes to be run using prescribed model state fields.
+  It can profitably be used to compare different bulk formulae or
+  adjust the parameters of a given bulk formula.
+
+More informations about SAS can be found in :doc:`NEMO manual <cite>`.
 
 SPITZ12
@@ -290,5 +318,5 @@
 ``SPITZ12`` is a regional configuration around the Svalbard archipelago
 at 1/12° of horizontal resolution and 75 vertical levels.
-See `Rousset et al. (2015) <https://www.geosci-model-dev.net/8/2991/2015/>`_ for more details.
+See :gmd:`Rousset et al. (2015) <8/2991/2015>` for more details.
 
 This configuration references to year 2002,
@@ -296,12 +324,8 @@
 while lateral boundary conditions for dynamical fields have 3 days time frequency.
 
-References
-==========
-
-.. bibliography:: configurations.bib
+.. rubric:: References
+
+.. bibliography:: cfgs.bib
    :all:
    :style: unsrt
    :labelprefix: C
-
-.. Links and substitutions
-
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/SHARED/README.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/SHARED/README.rst	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/SHARED/README.rst	(revision 12165)
@@ -3,39 +3,58 @@
 ***********
 
+.. todo::
+
+
+
 .. contents::
-           :local:
+   :local:
 
 Output of diagnostics in NEMO is usually done using XIOS.
-This is an efficient way of writing diagnostics because the time averaging, file writing and even some simple arithmetic or regridding is carried out in parallel to the NEMO model run.
+This is an efficient way of writing diagnostics because
+the time averaging, file writing and even some simple arithmetic or regridding is carried out in
+parallel to the NEMO model run.
 This page gives a basic introduction to using XIOS with NEMO.
-Much more information is available from the XIOS homepage above and from the NEMO manual.
+Much more information is available from the :xios:`XIOS homepage<>` above and from the NEMO manual.
 
-Use of XIOS for diagnostics is activated using the pre-compiler key ``key_iomput``. 
+Use of XIOS for diagnostics is activated using the pre-compiler key ``key_iomput``.
 
 Extracting and installing XIOS
-------------------------------
+==============================
 
 1. Install the NetCDF4 library.
-   If you want to use single file output you will need to compile the HDF & NetCDF libraries to allow parallel IO.
-2. Download the version of XIOS that you wish to use. The recommended version is now XIOS 2.5:
-   
-.. code-block:: console
+   If you want to use single file output you will need to compile the HDF & NetCDF libraries to
+   allow parallel IO.
+2. Download the version of XIOS that you wish to use.
+   The recommended version is now XIOS 2.5:
 
-	$ svn co http://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/branchs/xios-2.5 xios-2.5
+   .. code-block:: console
 
-and follow the instructions in `XIOS documentation <http://forge.ipsl.jussieu.fr/ioserver/wiki/documentation>`_ to compile it.
-   If you find problems at this stage, support can be found by subscribing to the `XIOS mailing list <http://forge.ipsl.jussieu.fr/mailman/listinfo.cgi/xios-users>`_ and sending a mail message to it. 
+      $ svn co http://forge.ipsl.jussieu.fr/ioserver/svn/XIOS/branchs/xios-2.5
+
+and follow the instructions in :xios:`XIOS documentation <wiki/documentation>` to compile it.
+If you find problems at this stage, support can be found by subscribing to
+the :xios:`XIOS mailing list <../mailman/listinfo.cgi/xios-users>` and sending a mail message to it.
 
 XIOS Configuration files
 ------------------------
 
-XIOS is controlled using xml input files that should be copied to your model run directory before running the model.
-Examples of these files can be found in the reference configurations (``cfgs``). The XIOS executable expects to find a file called ``iodef.xml`` in the model run directory.
-In NEMO we have made the decision to use include statements in the ``iodef.xml`` file to include ``field_def_nemo-oce.xml`` (for physics), ``field_def_nemo-ice.xml`` (for ice), ``field_def_nemo-pisces.xml`` (for biogeochemistry) and ``domain_def.xml`` from the /cfgs/SHARED directory.
-Most users will not need to modify ``domain_def.xml`` or ``field_def_nemo-???.xml`` unless they want to add new diagnostics to the NEMO code.
-The definition of the output files is organized into separate ``file_definition.xml`` files which are included in the ``iodef.xml`` file.
+XIOS is controlled using XML input files that should be copied to
+your model run directory before running the model.
+Examples of these files can be found in the reference configurations (:file:`./cfgs`).
+The XIOS executable expects to find a file called :file:`iodef.xml` in the model run directory.
+In NEMO we have made the decision to use include statements in the :file:`iodef.xml` file to include:
+
+- :file:`field_def_nemo-oce.xml` (for physics),
+- :file:`field_def_nemo-ice.xml` (for ice),
+- :file:`field_def_nemo-pisces.xml` (for biogeochemistry) and
+- :file:`domain_def.xml` from the :file:`./cfgs/SHARED` directory.
+
+Most users will not need to modify :file:`domain_def.xml` or :file:`field_def_nemo-???.xml` unless
+they want to add new diagnostics to the NEMO code.
+The definition of the output files is organized into separate :file:`file_definition.xml` files which
+are included in the :file:`iodef.xml` file.
 
 Modes
------
+=====
 
 Detached Mode
@@ -44,11 +63,12 @@
 In detached mode the XIOS executable is executed on separate cores from the NEMO model.
 This is the recommended method for using XIOS for realistic model runs.
-To use this mode set ``using_server`` to ``true`` at the bottom of the ``iodef.xml`` file:
+To use this mode set ``using_server`` to ``true`` at the bottom of the :file:`iodef.xml` file:
 
 .. code-block:: xml
 
-	<variable id="using_server" type="boolean">true</variable>
+   <variable id="using_server" type="boolean">true</variable>
 
-Make sure there is a copy (or link to) your XIOS executable in the working directory and in your job submission script allocate processors to XIOS.
+Make sure there is a copy (or link to) your XIOS executable in the working directory and
+in your job submission script allocate processors to XIOS.
 
 Attached Mode
@@ -56,33 +76,44 @@
 
 In attached mode XIOS runs on each of the cores used by NEMO.
-This method is less efficient than the detached mode but can be more convenient for testing or with small configurations.
-To activate this mode simply set ``using_server`` to false in the ``iodef.xml`` file
+This method is less efficient than the detached mode but can be more convenient for testing or
+with small configurations.
+To activate this mode simply set ``using_server`` to false in the :file:`iodef.xml` file
 
 .. code-block:: xml
 
-	<variable id="using_server" type="boolean">false</variable>
+   <variable id="using_server" type="boolean">false</variable>
 
 and don't allocate any cores to XIOS.
-Note that due to the different domain decompositions between XIOS and NEMO if the total number of cores is larger than the number of grid points in the j direction then the model run will fail.
+
+.. note::
+
+   Due to the different domain decompositions between XIOS and NEMO,
+   if the total number of cores is larger than the number of grid points in the ``j`` direction then
+   the model run will fail.
 
 Adding new diagnostics
-----------------------
+======================
 
 If you want to add a NEMO diagnostic to the NEMO code you will need to do the following:
 
 1. Add any necessary code to calculate you new diagnostic in NEMO
-2. Send the field to XIOS using ``CALL iom_put( 'field_id', variable )`` where ``field_id`` is a unique id for your new diagnostics and variable is the fortran variable containing the data.
-   This should be called at every model timestep regardless of how often you want to output the field. No time averaging should be done in the model code. 
-3. If it is computationally expensive to calculate your new diagnostic you should also use "iom_use" to determine if it is requested in the current model run. For example,
-   
-.. code-block:: fortran
+2. Send the field to XIOS using ``CALL iom_put( 'field_id', variable )`` where
+   ``field_id`` is a unique id for your new diagnostics and
+   variable is the fortran variable containing the data.
+   This should be called at every model timestep regardless of how often you want to output the field.
+   No time averaging should be done in the model code.
+3. If it is computationally expensive to calculate your new diagnostic
+   you should also use "iom_use" to determine if it is requested in the current model run.
+   For example,
 
-		IF iom_use('field_id') THEN
-		   !Some expensive computation
-	   	!...
-		   !...
-	   	iom_put('field_id', variable)
-		ENDIF
+   .. code-block:: fortran
 
-4. Add a variable definition to the ``field_def_nemo-???.xml`` file.
-5. Add the variable to the ``iodef.xml`` or ``file_definition.xml`` file.
+      IF iom_use('field_id') THEN
+         !Some expensive computation
+         !...
+         !...
+	 iom_put('field_id', variable)
+      ENDIF
+
+4. Add a variable definition to the :file:`field_def_nemo-???.xml` file.
+5. Add the variable to the :file:`iodef.xml` or :file:`file_definition.xml` file.
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/SHARED/namelist_ice_ref
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/SHARED/namelist_ice_ref	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/SHARED/namelist_ice_ref	(revision 12165)
@@ -57,10 +57,8 @@
    ln_landfast_L16  = .false.         !  landfast: parameterization from Lemieux 2016
       rn_depfra     =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast
-                                      !          recommended range: [0.1 ; 0.25] - L16=0.125 - home=0.15
-      rn_icebfr     =  15.            !        ln_landfast_L16:  maximum bottom stress per unit volume [N/m3]
-                                      !        ln_landfast_home: maximum bottom stress per unit area of contact [N/m2]
-                                      !          recommended range: ?? L16=15 - home=10
+                                      !          recommended range: [0.1 ; 0.25]
+      rn_icebfr     =  15.            !        maximum bottom stress per unit volume [N/m3]
       rn_lfrelax    =   1.e-5         !        relaxation time scale to reach static friction [s-1]
-      rn_tensile    =   0.2           !        ln_landfast_L16: isotropic tensile strength
+      rn_tensile    =   0.2           !        isotropic tensile strength [0-0.5??]
 /
 !------------------------------------------------------------------------------
@@ -103,6 +101,6 @@
 &namdyn_adv     !   Ice advection
 !------------------------------------------------------------------------------
-   ln_adv_Pra       = .false.         !  Advection scheme (Prather)
-   ln_adv_UMx       = .true.          !  Advection scheme (Ultimate-Macho)
+   ln_adv_Pra       = .true.         !  Advection scheme (Prather)
+   ln_adv_UMx       = .false.          !  Advection scheme (Ultimate-Macho)
       nn_UMx        =   5             !     order of the scheme for UMx (1-5 ; 20=centered 2nd order)
 /
@@ -234,8 +232,8 @@
 &namdia         !   Diagnostics
 !------------------------------------------------------------------------------
-   ln_icediachk     = .false.         !  check online the heat, mass & salt budgets at each time step
-      !                               !     rate of ice spuriously gained/lost. For ex., rn_icechk=1. <=> 1mm/year, rn_icechk=0.1 <=> 1mm/10years                                   
-      rn_icechk_cel =  1.             !     check at any gridcell           => stops the code if violated (and writes a file)
-      rn_icechk_glo =  0.1            !     check over the entire ice cover => only prints warnings
+   ln_icediachk     = .false.         !  check online heat, mass & salt budgets
+      !                               !   rate of ice spuriously gained/lost at each time step => rn_icechk=1 <=> 1.e-6 m/hour
+      rn_icechk_cel =  100.           !     check at each gridcell          (1.e-4m/h)=> stops the code if violated (and writes a file)
+      rn_icechk_glo =  1.             !     check over the entire ice cover (1.e-6m/h)=> only prints warnings
    ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F)
    ln_icectl        = .false.         !  ice points output for debug (T or F)
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/SPITZ12/EXPREF/namelist_ice_cfg
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/SPITZ12/EXPREF/namelist_ice_cfg	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/cfgs/SPITZ12/EXPREF/namelist_ice_cfg	(revision 12165)
@@ -44,10 +44,11 @@
 &namdyn_rhg     !   Ice rheology
 !------------------------------------------------------------------------------
-   ln_rhg_EVP       = .true.          !  EVP rheology
-      ln_aEVP       = .true.          !     adaptive rheology (Kimmritz et al. 2016 & 2017)
 /
 !------------------------------------------------------------------------------
 &namdyn_adv     !   Ice advection
 !------------------------------------------------------------------------------
+   ln_adv_Pra       = .false.         !  Advection scheme (Prather)
+   ln_adv_UMx       = .true.          !  Advection scheme (Ultimate-Macho)
+      nn_UMx        =   5             !     order of the scheme for UMx (1-5 ; 20=centered 2nd order)
 /
 !------------------------------------------------------------------------------
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/latex/global/ametsoc.bst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/latex/global/ametsoc.bst	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/latex/global/ametsoc.bst	(revision 12165)
@@ -9,5 +9,5 @@
 %% *** Bibliography style file for ALL AMS Journals...version 1.0  ***
 %% *** Brian Papa - American Meteorological Society ***
-%% 
+%%
 %% Copyright 1994-2004 Patrick W Daly
  % ===============================================================
@@ -519,5 +519,5 @@
   duplicate$ empty$ 'skip$
     {
-      "\href{http://dx.doi.org/" swap$ * "}{DOI}" *
+      "\href{http://dx.doi.org/" swap$ * "}{\aiDoi}" *
     }
   if$
@@ -1192,5 +1192,5 @@
   crossref missing$
     { format.in.ed.booktitle "booktitle" output.check
-      format.publisher.address output      
+      format.publisher.address output
       format.bvolume output
       format.number.series output
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/rst/source/conf.py
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/rst/source/conf.py	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/rst/source/conf.py	(revision 12165)
@@ -230,5 +230,5 @@
 texinfo_documents = [
   ('guide', 'NEMO', u'NEMO Documentation',
-   u'NEMO System Team', 'NEMO', 'One line description of project.',
+   u'NEMO System Team', 'NEMO', 'Community Ocean Model',
    'Miscellaneous'),
 ]
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/rst/source/docutils.conf
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/rst/source/docutils.conf	(revision 12165)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/rst/source/docutils.conf	(revision 12165)
@@ -0,0 +1,2 @@
+[restructuredtext parser]
+syntax_highlight = short
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/rst/source/global.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/rst/source/global.rst	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/doc/rst/source/global.rst	(revision 12165)
@@ -1,3 +1,5 @@
-.. Roles (custom styles related to CSS classes in 'source/_static/style.css')
+.. Roles
+
+.. custom styles related to CSS classes in './_static/style.css'
 
 .. role:: blue
@@ -5,4 +7,18 @@
 .. role:: grey
 .. role:: greysup(sup)
+
+.. inline code snippets
+
+.. role:: python(code)
+   :language: python
+   :class: highlight
+
+.. role:: fortran(code)
+   :language: fortran
+   :class: highlight
+
+.. role:: console(code)
+   :language: console
+   :class: highlight
 
 .. Substitutions
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/ice.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/ice.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/ice.F90	(revision 12165)
@@ -328,12 +328,12 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sz_i     !: ice salinity          [PSS]
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip       !: melt pond fraction per grid cell area
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip       !: melt pond concentration
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_ip       !: melt pond volume per grid cell area      [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_frac  !: melt pond volume per ice area
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_ip       !: melt pond thickness                      [m]
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond fraction
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_frac  !: melt pond fraction (a_ip/a_i)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_ip       !: melt pond depth                          [m]
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond concentration
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hm_ip      !: mean melt pond depth                     [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vt_ip      !: total melt pond volume per unit area     [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vt_ip      !: total melt pond volume per gridcell area [m]
 
    !!----------------------------------------------------------------------
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icectl.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icectl.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icectl.F90	(revision 12165)
@@ -44,9 +44,9 @@
    PUBLIC   ice_prt3D
 
-   ! thresold values for conservation
+   ! thresold rates for conservation
    !    these values are changed by the namelist parameter rn_icechk, so that threshold = zchk * rn_icechk
-   REAL(wp), PARAMETER ::   zchk_m   = 1.e-5   ! kg/m2/s <=> 1mm of ice per year  spuriously gained/lost
-   REAL(wp), PARAMETER ::   zchk_s   = 1.e-4   ! g/m2/s  <=> 1mm of ice per year  spuriously gained/lost (considering s=10g/kg)
-   REAL(wp), PARAMETER ::   zchk_t   = 3.      ! W/m2    <=> 1mm of ice per year  spuriously gained/lost (considering Lf=3e5J/kg)
+   REAL(wp), PARAMETER ::   zchk_m   = 2.5e-7   ! kg/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost
+   REAL(wp), PARAMETER ::   zchk_s   = 2.5e-6   ! g/m2/s  <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg)
+   REAL(wp), PARAMETER ::   zchk_t   = 7.5e-2   ! W/m2    <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg)
    
    !! * Substitutions
@@ -68,6 +68,5 @@
       !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true
       !!              It prints in ocean.output if there is a violation of conservation at each time-step
-      !!              The thresholds (zchk_m, zchk_s, zchk_t) which determine violations are set to
-      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
+      !!              The thresholds (zchk_m, zchk_s, zchk_t) determine violations
       !!              For salt and heat thresholds, ice is considered to have a salinity of 10 
       !!              and a heat content of 3e5 J/kg (=latent heat of fusion) 
@@ -133,6 +132,6 @@
 
          ! -- advection scheme is conservative? -- !
-         zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0
-         zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t ) ! must be close to 0
+         zvtrp = glob_sum( 'icectl', ( diag_trp_vi * rhoi + diag_trp_vs * rhos ) * e1e2t ) ! must be close to 0 (only for Prather)
+         zetrp = glob_sum( 'icectl', ( diag_trp_ei        + diag_trp_es        ) * e1e2t ) ! must be close to 0 (only for Prather)
 
          ! ice area (+epsi10 to set a threshold > 0 when there is no ice) 
@@ -157,6 +156,9 @@
                &                   WRITE(numout,*)   cd_routine,' : violation a_i > amax      = ',zdiag_amax
             ! check if advection scheme is conservative
-            IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) &
-               &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice
+            !    only check for Prather because Ultimate-Macho uses corrective fluxes (wfx etc)
+            !    so the formulation for conservation is different (and not coded) 
+            !    it does not mean UM is not conservative (it is checked with above prints) => update (09/2019): same for Prather now
+            !IF( ln_adv_Pra .AND. ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) &
+            !   &                   WRITE(numout,*)   cd_routine,' : violation adv scheme [kg] = ',zvtrp * rdt_ice
          ENDIF
          !
@@ -173,6 +175,5 @@
       !! ** Method  : This is an online diagnostics which can be activated with ln_icediachk=true
       !!              It prints in ocean.output if there is a violation of conservation at each time-step
-      !!              The thresholds (zchk_m, zchk_s, zchk_t) which determine the violation are set to
-      !!              a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
+      !!              The thresholds (zchk_m, zchk_s, zchk_t) determine the violations
       !!              For salt and heat thresholds, ice is considered to have a salinity of 10 
       !!              and a heat content of 3e5 J/kg (=latent heat of fusion) 
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_adv_pra.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_adv_pra.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_adv_pra.F90	(revision 12165)
@@ -16,7 +16,9 @@
    !!   adv_pra_rst     : read/write Prather field in ice restart file, or initialized to zero
    !!----------------------------------------------------------------------
+   USE phycst         ! physical constant
    USE dom_oce        ! ocean domain
    USE ice            ! sea-ice variables
    USE sbc_oce , ONLY : nn_fsbc   ! frequency of sea-ice call
+   USE icevar         ! sea-ice: operations
    !
    USE in_out_manager ! I/O manager
@@ -25,5 +27,4 @@
    USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
    USE lbclnk         ! lateral boundary conditions (or mpp links)
-   USE prtctl         ! Print control
 
    IMPLICIT NONE
@@ -36,8 +37,7 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxice, syice, sxxice, syyice, sxyice   ! ice thickness 
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsn , sysn , sxxsn , syysn , sxysn    ! snow thickness
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxa  , sya  , sxxa  , syya  , sxya     ! lead fraction
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxa  , sya  , sxxa  , syya  , sxya     ! ice concentration
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxsal, sysal, sxxsal, syysal, sxysal   ! ice salinity
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxage, syage, sxxage, syyage, sxyage   ! ice age
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   sxopw, syopw, sxxopw, syyopw, sxyopw   ! open water in sea ice
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxc0 , syc0 , sxxc0 , syyc0 , sxyc0    ! snow layers heat content
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sxe  , sye  , sxxe  , syye  , sxye     ! ice layers heat content
@@ -81,184 +81,167 @@
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
       !
-      INTEGER  ::   jk, jl, jt              ! dummy loop indices
-      INTEGER  ::   initad                  ! number of sub-timestep for the advection
-      REAL(wp) ::   zcfl , zusnit           !   -      -
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zarea
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0opw
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ice, z0snw, z0ai, z0smi, z0oi
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ap , z0vp
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0es
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0ei
+      INTEGER  ::   ji,jj, jk, jl, jt       ! dummy loop indices
+      INTEGER  ::   icycle                  ! number of sub-timestep for the advection
+      REAL(wp) ::   zdt                     !   -      -
+      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication
+      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2
+      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx
+      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea
+      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi
+      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp
+      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es
+      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei
       !!----------------------------------------------------------------------
       !
       IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme'
       !
-      ALLOCATE( zarea(jpi,jpj)    , z0opw(jpi,jpj, 1 ) , z0ice(jpi,jpj,jpl) , z0snw(jpi,jpj,jpl) ,                       &
-         &      z0ai(jpi,jpj,jpl) , z0smi(jpi,jpj,jpl) , z0oi (jpi,jpj,jpl) , z0ap (jpi,jpj,jpl) , z0vp(jpi,jpj,jpl) ,   &
-         &      z0es (jpi,jpj,nlay_s,jpl), z0ei(jpi,jpj,nlay_i,jpl) )
-      !
-      ! --- If ice drift field is too fast, use an appropriate time step for advection (CFL test for stability) --- !        
-      zcfl  =            MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
-      zcfl  = MAX( zcfl, MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
-      CALL mpp_max( 'icedyn_adv_pra', zcfl )
+      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- !
+      !        Note: the advection split is applied at the next time-step in order to avoid blocking global comm.
+      !              this should not affect too much the stability
+      zcflnow(1) =                  MAXVAL( ABS( pu_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )
+      zcflnow(1) = MAX( zcflnow(1), MAXVAL( ABS( pv_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
       
-      IF( zcfl > 0.5 ) THEN   ;   initad = 2   ;   zusnit = 0.5_wp
-      ELSE                    ;   initad = 1   ;   zusnit = 1.0_wp
+      ! non-blocking global communication send zcflnow and receive zcflprv
+      CALL mpp_delay_max( 'icedyn_adv_pra', 'cflice', zcflnow(:), zcflprv(:), kt == nitend - nn_fsbc + 1 )
+
+      IF( zcflprv(1) > .5 ) THEN   ;   icycle = 2
+      ELSE                         ;   icycle = 1
       ENDIF
+      zdt = rdt_ice / REAL(icycle)
       
-      zarea(:,:) = e1e2t(:,:)
-      !-------------------------
-      ! transported fields                                        
-      !-------------------------
-      z0opw(:,:,1) = pato_i(:,:) * e1e2t(:,:)              ! Open water area 
-      DO jl = 1, jpl
-         z0snw(:,:,jl) = pv_s (:,:,  jl) * e1e2t(:,:)     ! Snow volume
-         z0ice(:,:,jl) = pv_i (:,:,  jl) * e1e2t(:,:)     ! Ice  volume
-         z0ai (:,:,jl) = pa_i (:,:,  jl) * e1e2t(:,:)     ! Ice area
-         z0smi(:,:,jl) = psv_i(:,:,  jl) * e1e2t(:,:)     ! Salt content
-         z0oi (:,:,jl) = poa_i(:,:,  jl) * e1e2t(:,:)     ! Age content
-         DO jk = 1, nlay_s
-            z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content
-         END DO
-         DO jk = 1, nlay_i
-            z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
-         END DO
-         IF ( ln_pnd_H12 ) THEN
-            z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction
-            z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume
+      ! --- transport --- !
+      zudy(:,:) = pu_ice(:,:) * e2u(:,:)
+      zvdx(:,:) = pv_ice(:,:) * e1v(:,:)
+
+      DO jt = 1, icycle
+
+         ! record at_i before advection (for open water)
+         zati1(:,:) = SUM( pa_i(:,:,:), dim=3 )
+         
+         ! --- transported fields --- !                                        
+         DO jl = 1, jpl
+            zarea(:,:,jl) = e1e2t(:,:)
+            z0snw(:,:,jl) = pv_s (:,:,jl) * e1e2t(:,:)        ! Snow volume
+            z0ice(:,:,jl) = pv_i (:,:,jl) * e1e2t(:,:)        ! Ice  volume
+            z0ai (:,:,jl) = pa_i (:,:,jl) * e1e2t(:,:)        ! Ice area
+            z0smi(:,:,jl) = psv_i(:,:,jl) * e1e2t(:,:)        ! Salt content
+            z0oi (:,:,jl) = poa_i(:,:,jl) * e1e2t(:,:)        ! Age content
+            DO jk = 1, nlay_s
+               z0es(:,:,jk,jl) = pe_s(:,:,jk,jl) * e1e2t(:,:) ! Snow heat content
+            END DO
+            DO jk = 1, nlay_i
+               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
+            END DO
+            IF ( ln_pnd_H12 ) THEN
+               z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction
+               z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume
+            ENDIF
+         END DO
+         !
+         !                                                                  !--------------------------------------------!
+         IF( MOD( (kt - 1) / nn_fsbc , 2 ) ==  MOD( (jt - 1) , 2 ) ) THEN   !==  odd ice time step:  adv_x then adv_y  ==!
+            !                                                               !--------------------------------------------!
+            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
+            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
+            CALL adv_x( zdt , zudy , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume
+            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  )
+            CALL adv_x( zdt , zudy , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
+            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
+            CALL adv_x( zdt , zudy , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration
+            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   )
+            CALL adv_x( zdt , zudy , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
+            CALL adv_y( zdt , zvdx , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage )
+            !
+            DO jk = 1, nlay_s                                                                           !--- snow heat content
+               CALL adv_x( zdt, zudy, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
+                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
+               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
+                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
+            END DO
+            DO jk = 1, nlay_i                                                                           !--- ice heat content
+               CALL adv_x( zdt, zudy, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
+                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
+               CALL adv_y( zdt, zvdx, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
+                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
+            END DO
+            !
+            IF ( ln_pnd_H12 ) THEN
+               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
+               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
+               CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
+               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
+            ENDIF
+            !                                                               !--------------------------------------------!
+         ELSE                                                               !== even ice time step:  adv_y then adv_x  ==!
+            !                                                               !--------------------------------------------!
+            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice ) !--- ice volume
+            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ice , sxice , sxxice , syice , syyice , sxyice )
+            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  ) !--- snow volume
+            CALL adv_x( zdt , zudy , 0._wp , zarea , z0snw , sxsn  , sxxsn  , sysn  , syysn  , sxysn  )
+            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal ) !--- ice salinity
+            CALL adv_x( zdt , zudy , 0._wp , zarea , z0smi , sxsal , sxxsal , sysal , syysal , sxysal )
+            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   ) !--- ice concentration
+            CALL adv_x( zdt , zudy , 0._wp , zarea , z0ai  , sxa   , sxxa   , sya   , syya   , sxya   )
+            CALL adv_y( zdt , zvdx , 1._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage ) !--- ice age
+            CALL adv_x( zdt , zudy , 0._wp , zarea , z0oi  , sxage , sxxage , syage , syyage , sxyage )
+            DO jk = 1, nlay_s                                                                           !--- snow heat content
+               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
+                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
+               CALL adv_x( zdt, zudy, 0._wp, zarea, z0es (:,:,jk,:), sxc0(:,:,jk,:),   &
+                  &                                 sxxc0(:,:,jk,:), syc0(:,:,jk,:), syyc0(:,:,jk,:), sxyc0(:,:,jk,:) )
+            END DO
+            DO jk = 1, nlay_i                                                                           !--- ice heat content
+               CALL adv_y( zdt, zvdx, 1._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
+                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
+               CALL adv_x( zdt, zudy, 0._wp, zarea, z0ei(:,:,jk,:), sxe(:,:,jk,:),   & 
+                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) )
+            END DO
+            IF ( ln_pnd_H12 ) THEN
+               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction
+               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )
+               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
+               CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
+            ENDIF
+            !
          ENDIF
+
+         ! --- Recover the properties from their contents --- !
+         DO jl = 1, jpl
+            pv_i (:,:,jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            pv_s (:,:,jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            psv_i(:,:,jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            poa_i(:,:,jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            pa_i (:,:,jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            DO jk = 1, nlay_s
+               pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            END DO
+            DO jk = 1, nlay_i
+               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            END DO
+            IF ( ln_pnd_H12 ) THEN
+               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            ENDIF
+         END DO
+         !
+         ! derive open water from ice concentration
+         zati2(:,:) = SUM( pa_i(:,:,:), dim=3 )
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               pato_i(ji,jj) = pato_i(ji,jj) - ( zati2(ji,jj) - zati1(ji,jj) ) &                        !--- open water
+                  &                          - ( zudy(ji,jj) - zudy(ji-1,jj) + zvdx(ji,jj) - zvdx(ji,jj-1) ) * r1_e1e2t(ji,jj) * zdt
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_adv_pra', pato_i, 'T',  1. )
+         !
+         ! --- Ensure non-negative fields --- !
+         !     Remove negative values (conservation is ensured)
+         !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
+         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+         !
+         ! --- Ensure snow load is not too big --- !
+         CALL Hsnow( zdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
+         !
       END DO
-
-      !                                                    !--------------------------------------------!
-      IF( MOD( ( kt - 1) / nn_fsbc , 2 ) == 0 ) THEN       !==  odd ice time step:  adv_x then adv_y  ==!
-         !                                                 !--------------------------------------------!
-         DO jt = 1, initad
-            CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
-               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-            CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
-               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-            DO jl = 1, jpl
-               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
-                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
-                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
-                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
-                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
-                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
-                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &    !--- ice age      ---     
-                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
-                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &    !--- ice concentrations ---
-                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   & 
-                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               DO jk = 1, nlay_s                                                               !--- snow heat contents ---
-                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   &
-                     &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )
-                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   &
-                     &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )
-               END DO
-               DO jk = 1, nlay_i                                                               !--- ice heat contents ---
-                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   & 
-                     &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   & 
-                     &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-               END DO
-               IF ( ln_pnd_H12 ) THEN
-                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction --
-                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
-                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   & 
-                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
-                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &    !--- melt pond volume   --
-                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
-                  CALL adv_y( zusnit, pv_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   & 
-                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
-               ENDIF
-            END DO
-         END DO
-      !                                                    !--------------------------------------------!
-      ELSE                                                 !== even ice time step:  adv_y then adv_x  ==!
-         !                                                 !--------------------------------------------!
-         DO jt = 1, initad
-            CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &             !--- ice open water area
-               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-            CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0opw (:,:,1), sxopw(:,:),   &
-               &                                      sxxopw(:,:)  , syopw(:,:), syyopw(:,:), sxyopw(:,:)  )
-            DO jl = 1, jpl
-               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &    !--- ice volume  ---
-                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ice (:,:,jl), sxice(:,:,jl),   &
-                  &                                      sxxice(:,:,jl), syice(:,:,jl), syyice(:,:,jl), sxyice(:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &    !--- snow volume  ---
-                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0snw (:,:,jl), sxsn (:,:,jl),   &
-                  &                                      sxxsn (:,:,jl), sysn (:,:,jl), syysn (:,:,jl), sxysn (:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &    !--- ice salinity ---
-                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0smi (:,:,jl), sxsal(:,:,jl),   &
-                  &                                      sxxsal(:,:,jl), sysal(:,:,jl), syysal(:,:,jl), sxysal(:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &   !--- ice age      ---
-                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0oi  (:,:,jl), sxage(:,:,jl),   &
-                  &                                      sxxage(:,:,jl), syage(:,:,jl), syyage(:,:,jl), sxyage(:,:,jl)  )
-               CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &   !--- ice concentrations ---
-                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ai  (:,:,jl), sxa  (:,:,jl),   &
-                  &                                      sxxa  (:,:,jl), sya  (:,:,jl), syya  (:,:,jl), sxya  (:,:,jl)  )
-               DO jk = 1, nlay_s                                                             !--- snow heat contents ---
-                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   &
-                     &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )
-                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0es (:,:,jk,jl), sxc0(:,:,jk,jl),   &
-                     &                                      sxxc0(:,:,jk,jl), syc0(:,:,jk,jl), syyc0(:,:,jk,jl), sxyc0(:,:,jk,jl) )
-               END DO
-               DO jk = 1, nlay_i                                                             !--- ice heat contents ---
-                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   & 
-                     &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ei(:,:,jk,jl), sxe(:,:,jk,jl),   & 
-                     &                                      sxxe(:,:,jk,jl), sye(:,:,jk,jl), syye(:,:,jk,jl), sxye(:,:,jk,jl) )
-               END DO
-               IF ( ln_pnd_H12 ) THEN
-                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction ---
-                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
-                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &
-                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
-                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &   !--- melt pond volume   ---
-                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
-                  CALL adv_x( zusnit, pu_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &
-                     &                                      sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
-               ENDIF
-            END DO
-         END DO
-      ENDIF
-
-      !-------------------------------------------
-      ! Recover the properties from their contents
-      !-------------------------------------------
-      pato_i(:,:) = z0opw(:,:,1) * r1_e1e2t(:,:) * tmask(:,:,1)
-      DO jl = 1, jpl
-         pv_i (:,:,  jl) = z0ice(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         pv_s (:,:,  jl) = z0snw(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         psv_i(:,:,  jl) = z0smi(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         poa_i(:,:,  jl) = z0oi (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         pa_i (:,:,  jl) = z0ai (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         DO jk = 1, nlay_s
-            pe_s(:,:,jk,jl) = z0es(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         END DO
-         DO jk = 1, nlay_i
-            pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         END DO
-         IF ( ln_pnd_H12 ) THEN
-            pa_ip  (:,:,jl) = z0ap (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-            pv_ip  (:,:,jl) = z0vp (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
-         ENDIF
-      END DO
-      !
-      DEALLOCATE( zarea , z0opw , z0ice, z0snw , z0ai , z0smi , z0oi , z0ap , z0vp , z0es, z0ei )
       !
       IF( lrst_ice )   CALL adv_pra_rst( 'WRITE', kt )   !* write Prather fields in the restart file
@@ -267,5 +250,5 @@
    
    
-   SUBROUTINE adv_x( pdf, put , pcrh, psm , ps0 ,   &
+   SUBROUTINE adv_x( pdt, put , pcrh, psm , ps0 ,   &
       &              psx, psxx, psy , psyy, psxy )
       !!----------------------------------------------------------------------
@@ -275,14 +258,14 @@
       !!                variable on x axis
       !!----------------------------------------------------------------------
-      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
-      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
+      REAL(wp)                  , INTENT(in   ) ::   pdt                ! the time step
+      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
+      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   put                ! i-direction ice velocity at U-point [m/s]
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments 
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
       !! 
-      INTEGER  ::   ji, jj                               ! dummy loop indices
-      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! local scalars
+      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
+      REAL(wp) ::   zs1max, zslpmax, ztemp               ! local scalars
       REAL(wp) ::   zs1new, zalf , zalfq , zbt           !   -      -
       REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !   -      -
@@ -291,154 +274,150 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q         !  -      -
       !-----------------------------------------------------------------------
-
-      ! Limitation of moments.                                           
-
-      zrdt = rdt_ice * pdf      ! If ice drift field is too fast, use an appropriate time step for advection.
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zslpmax = MAX( 0._wp, ps0(ji,jj) )
-            zs1max  = 1.5 * zslpmax
-            zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj) ) )
-            zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
-               &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
-            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
-
-            ps0 (ji,jj) = zslpmax  
-            psx (ji,jj) = zs1new      * rswitch
-            psxx(ji,jj) = zs2new      * rswitch
-            psy (ji,jj) = psy (ji,jj) * rswitch
-            psyy(ji,jj) = psyy(ji,jj) * rswitch
-            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
-         END DO
+      !
+      jcat = SIZE( ps0 , 3 )   ! size of input arrays
+      !
+      DO jl = 1, jcat   ! loop on categories
+         !
+         ! Limitation of moments.                                           
+         DO jj = 2, jpjm1
+            DO ji = 1, jpi
+               !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
+               psm (ji,jj,jl) = MAX( pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20 )
+               !
+               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
+               zs1max  = 1.5 * zslpmax
+               zs1new  = MIN( zs1max, MAX( -zs1max, psx(ji,jj,jl) ) )
+               zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
+                  &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj,jl) )  )
+               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
+
+               ps0 (ji,jj,jl) = zslpmax  
+               psx (ji,jj,jl) = zs1new         * rswitch
+               psxx(ji,jj,jl) = zs2new         * rswitch
+               psy (ji,jj,jl) = psy (ji,jj,jl) * rswitch
+               psyy(ji,jj,jl) = psyy(ji,jj,jl) * rswitch
+               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
+            END DO
+         END DO
+
+         !  Calculate fluxes and moments between boxes i<-->i+1              
+         DO jj = 2, jpjm1                      !  Flux from i to i+1 WHEN u GT 0 
+            DO ji = 1, jpi
+               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
+               zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl)
+               zalfq        =  zalf * zalf
+               zalf1        =  1.0 - zalf
+               zalf1q       =  zalf1 * zalf1
+               !
+               zfm (ji,jj)  =  zalf  *   psm (ji,jj,jl)
+               zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj,jl) + zalf1 * ( psx(ji,jj,jl) + (zalf1 - zalf) * psxx(ji,jj,jl) ) )
+               zfx (ji,jj)  =  zalfq * ( psx (ji,jj,jl) + 3.0 * zalf1 * psxx(ji,jj,jl) )
+               zfxx(ji,jj)  =  zalf  *   psxx(ji,jj,jl) * zalfq
+               zfy (ji,jj)  =  zalf  * ( psy (ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
+               zfxy(ji,jj)  =  zalfq *   psxy(ji,jj,jl)
+               zfyy(ji,jj)  =  zalf  *   psyy(ji,jj,jl)
+
+               !  Readjust moments remaining in the box.
+               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
+               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
+               psx (ji,jj,jl)  =  zalf1q * ( psx(ji,jj,jl) - 3.0 * zalf * psxx(ji,jj,jl) )
+               psxx(ji,jj,jl)  =  zalf1  * zalf1q * psxx(ji,jj,jl)
+               psy (ji,jj,jl)  =  psy (ji,jj,jl) - zfy(ji,jj)
+               psyy(ji,jj,jl)  =  psyy(ji,jj,jl) - zfyy(ji,jj)
+               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
+            END DO
+         END DO
+
+         DO jj = 2, jpjm1                      !  Flux from i+1 to i when u LT 0.
+            DO ji = 1, fs_jpim1
+               zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl) 
+               zalg  (ji,jj) = zalf
+               zalfq         = zalf * zalf
+               zalf1         = 1.0 - zalf
+               zalg1 (ji,jj) = zalf1
+               zalf1q        = zalf1 * zalf1
+               zalg1q(ji,jj) = zalf1q
+               !
+               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji+1,jj,jl)
+               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji+1,jj,jl) &
+                  &                                   - zalf1 * ( psx(ji+1,jj,jl) - (zalf1 - zalf ) * psxx(ji+1,jj,jl) ) )
+               zfx   (ji,jj) = zfx (ji,jj) + zalfq * (  psx (ji+1,jj,jl) - 3.0 * zalf1 * psxx(ji+1,jj,jl) )
+               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji+1,jj,jl) * zalfq
+               zfy   (ji,jj) = zfy (ji,jj) + zalf  * (  psy (ji+1,jj,jl) - zalf1 * psxy(ji+1,jj,jl) )
+               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji+1,jj,jl)
+               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji+1,jj,jl)
+            END DO
+         END DO
+
+         DO jj = 2, jpjm1                     !  Readjust moments remaining in the box. 
+            DO ji = fs_2, fs_jpim1
+               zbt  =       zbet(ji-1,jj)
+               zbt1 = 1.0 - zbet(ji-1,jj)
+               !
+               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji-1,jj) )
+               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji-1,jj) )
+               psx (ji,jj,jl) = zalg1q(ji-1,jj) * ( psx(ji,jj,jl) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj,jl) )
+               psxx(ji,jj,jl) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj,jl)
+               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) - zfy (ji-1,jj) )
+               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) - zfyy(ji-1,jj) )
+               psxy(ji,jj,jl) = zalg1q(ji-1,jj) * psxy(ji,jj,jl)
+            END DO
+         END DO
+
+         !   Put the temporary moments into appropriate neighboring boxes.    
+         DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
+            DO ji = fs_2, fs_jpim1
+               zbt  =       zbet(ji-1,jj)
+               zbt1 = 1.0 - zbet(ji-1,jj)
+               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj,jl)
+               zalf          = zbt * zfm(ji-1,jj) / psm(ji,jj,jl)
+               zalf1         = 1.0 - zalf
+               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji-1,jj)
+               !
+               ps0 (ji,jj,jl) =  zbt  * ( ps0(ji,jj,jl) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj,jl)
+               psx (ji,jj,jl) =  zbt  * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp ) + zbt1 * psx(ji,jj,jl)
+               psxx(ji,jj,jl) =  zbt  * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj,jl)                             &
+                  &                     + 5.0 * ( zalf * zalf1 * ( psx (ji,jj,jl) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  ) &
+                  &            + zbt1 * psxx(ji,jj,jl)
+               psxy(ji,jj,jl) =  zbt  * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj,jl)             &
+                  &                     + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj,jl) ) )   &
+                  &            + zbt1 * psxy(ji,jj,jl)
+               psy (ji,jj,jl) =  zbt  * ( psy (ji,jj,jl) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj,jl)
+               psyy(ji,jj,jl) =  zbt  * ( psyy(ji,jj,jl) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj,jl)
+            END DO
+         END DO
+
+         DO jj = 2, jpjm1                      !  Flux from i+1 to i IF u LT 0.
+            DO ji = fs_2, fs_jpim1
+               zbt  =       zbet(ji,jj)
+               zbt1 = 1.0 - zbet(ji,jj)
+               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
+               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
+               zalf1         = 1.0 - zalf
+               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
+               !
+               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) + zf0(ji,jj) )
+               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj,jl) + 3.0 * ztemp )
+               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( zalf * zalf * zfxx(ji,jj) + zalf1 * zalf1 * psxx(ji,jj,jl) &
+                  &                                           + 5.0 * ( zalf * zalf1 * ( - psx(ji,jj,jl) + zfx(ji,jj) )    &
+                  &                                           + ( zalf1 - zalf ) * ztemp ) )
+               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * ( zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
+                  &                                           + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj,jl) ) )
+               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * ( psy (ji,jj,jl) + zfy (ji,jj) )
+               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * ( psyy(ji,jj,jl) + zfyy(ji,jj) )
+            END DO
+         END DO
+
       END DO
 
-      !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
-      psm (:,:)  = MAX( pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
-
-      !  Calculate fluxes and moments between boxes i<-->i+1              
-      DO jj = 1, jpj                      !  Flux from i to i+1 WHEN u GT 0 
-         DO ji = 1, jpi
-            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) )
-            zalf         =  MAX( 0._wp, put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji,jj)
-            zalfq        =  zalf * zalf
-            zalf1        =  1.0 - zalf
-            zalf1q       =  zalf1 * zalf1
-            !
-            zfm (ji,jj)  =  zalf  *   psm (ji,jj)
-            zf0 (ji,jj)  =  zalf  * ( ps0 (ji,jj) + zalf1 * ( psx(ji,jj) + (zalf1 - zalf) * psxx(ji,jj) )  )
-            zfx (ji,jj)  =  zalfq * ( psx (ji,jj) + 3.0 * zalf1 * psxx(ji,jj) )
-            zfxx(ji,jj)  =  zalf  *   psxx(ji,jj) * zalfq
-            zfy (ji,jj)  =  zalf  * ( psy (ji,jj) + zalf1 * psxy(ji,jj) )
-            zfxy(ji,jj)  =  zalfq *   psxy(ji,jj)
-            zfyy(ji,jj)  =  zalf  *   psyy(ji,jj)
-
-            !  Readjust moments remaining in the box.
-            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
-            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
-            psx (ji,jj)  =  zalf1q * ( psx(ji,jj) - 3.0 * zalf * psxx(ji,jj) )
-            psxx(ji,jj)  =  zalf1  * zalf1q * psxx(ji,jj)
-            psy (ji,jj)  =  psy (ji,jj) - zfy(ji,jj)
-            psyy(ji,jj)  =  psyy(ji,jj) - zfyy(ji,jj)
-            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
-         END DO
-      END DO
-
-      DO jj = 1, jpjm1                      !  Flux from i+1 to i when u LT 0.
-         DO ji = 1, fs_jpim1
-            zalf          = MAX( 0._wp, -put(ji,jj) ) * zrdt * e2u(ji,jj) / psm(ji+1,jj) 
-            zalg  (ji,jj) = zalf
-            zalfq         = zalf * zalf
-            zalf1         = 1.0 - zalf
-            zalg1 (ji,jj) = zalf1
-            zalf1q        = zalf1 * zalf1
-            zalg1q(ji,jj) = zalf1q
-            !
-            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji+1,jj)
-            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji+1,jj) - zalf1 * ( psx(ji+1,jj) - (zalf1 - zalf ) * psxx(ji+1,jj) ) )
-            zfx   (ji,jj) = zfx (ji,jj) + zalfq * ( psx (ji+1,jj) - 3.0 * zalf1 * psxx(ji+1,jj) )
-            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji+1,jj) * zalfq
-            zfy   (ji,jj) = zfy (ji,jj) + zalf  * ( psy (ji+1,jj) - zalf1 * psxy(ji+1,jj) )
-            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji+1,jj)
-            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji+1,jj)
-         END DO
-      END DO
-
-      DO jj = 2, jpjm1                     !  Readjust moments remaining in the box. 
-         DO ji = fs_2, fs_jpim1
-            zbt  =       zbet(ji-1,jj)
-            zbt1 = 1.0 - zbet(ji-1,jj)
-            !
-            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji-1,jj) )
-            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji-1,jj) )
-            psx (ji,jj) = zalg1q(ji-1,jj) * ( psx(ji,jj) + 3.0 * zalg(ji-1,jj) * psxx(ji,jj) )
-            psxx(ji,jj) = zalg1 (ji-1,jj) * zalg1q(ji-1,jj) * psxx(ji,jj)
-            psy (ji,jj) = zbt * psy (ji,jj) + zbt1 * ( psy (ji,jj) - zfy (ji-1,jj) )
-            psyy(ji,jj) = zbt * psyy(ji,jj) + zbt1 * ( psyy(ji,jj) - zfyy(ji-1,jj) )
-            psxy(ji,jj) = zalg1q(ji-1,jj) * psxy(ji,jj)
-         END DO
-      END DO
-
-      !   Put the temporary moments into appropriate neighboring boxes.    
-      DO jj = 2, jpjm1                     !   Flux from i to i+1 IF u GT 0.
-         DO ji = fs_2, fs_jpim1
-            zbt  =       zbet(ji-1,jj)
-            zbt1 = 1.0 - zbet(ji-1,jj)
-            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji-1,jj) ) + zbt1 * psm(ji,jj)
-            zalf        = zbt * zfm(ji-1,jj) / psm(ji,jj)
-            zalf1       = 1.0 - zalf
-            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji-1,jj)
-            !
-            ps0 (ji,jj) = zbt * ( ps0(ji,jj) + zf0(ji-1,jj) ) + zbt1 * ps0(ji,jj)
-            psx (ji,jj) = zbt * ( zalf * zfx(ji-1,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp ) + zbt1 * psx(ji,jj)
-            psxx(ji,jj) = zbt * ( zalf * zalf * zfxx(ji-1,jj) + zalf1 * zalf1 * psxx(ji,jj)                               &
-               &                + 5.0 * ( zalf * zalf1 * ( psx (ji,jj) - zfx(ji-1,jj) ) - ( zalf1 - zalf ) * ztemp )  )   &
-               &                                                + zbt1 * psxx(ji,jj)
-            psxy(ji,jj) = zbt * ( zalf * zfxy(ji-1,jj) + zalf1 * psxy(ji,jj)             &
-               &                + 3.0 * (- zalf1*zfy(ji-1,jj)  + zalf * psy(ji,jj) ) )   &
-               &                                                + zbt1 * psxy(ji,jj)
-            psy (ji,jj) = zbt * ( psy (ji,jj) + zfy (ji-1,jj) ) + zbt1 * psy (ji,jj)
-            psyy(ji,jj) = zbt * ( psyy(ji,jj) + zfyy(ji-1,jj) ) + zbt1 * psyy(ji,jj)
-         END DO
-      END DO
-
-      DO jj = 2, jpjm1                     !  Flux from i+1 to i IF u LT 0.
-         DO ji = fs_2, fs_jpim1
-            zbt  =       zbet(ji,jj)
-            zbt1 = 1.0 - zbet(ji,jj)
-            psm(ji,jj)  = zbt * psm(ji,jj)  + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
-            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
-            zalf1       = 1.0 - zalf
-            ztemp       = - zalf * ps0(ji,jj) + zalf1 * zf0(ji,jj)
-            !
-            ps0(ji,jj)  = zbt * ps0 (ji,jj) + zbt1 * ( ps0(ji,jj) + zf0(ji,jj) )
-            psx(ji,jj)  = zbt * psx (ji,jj) + zbt1 * ( zalf * zfx(ji,jj) + zalf1 * psx(ji,jj) + 3.0 * ztemp )
-            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( zalf * zalf * zfxx(ji,jj)  + zalf1 * zalf1 * psxx(ji,jj)  &
-               &                                      + 5.0 *( zalf * zalf1 * ( - psx(ji,jj) + zfx(ji,jj) )      &
-               &                                      + ( zalf1 - zalf ) * ztemp ) )
-            psxy(ji,jj) = zbt * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)  &
-               &                                      + 3.0 * ( zalf1 * zfy(ji,jj) - zalf * psy(ji,jj) )  )
-            psy(ji,jj)  = zbt * psy (ji,jj)  + zbt1 * ( psy (ji,jj) + zfy (ji,jj) )
-            psyy(ji,jj) = zbt * psyy(ji,jj)  + zbt1 * ( psyy(ji,jj) + zfyy(ji,jj) )
-         END DO
-      END DO
-
       !-- Lateral boundary conditions
-      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm , 'T',  1., ps0 , 'T',  1.   &
-         &              , psx , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
-         &              , psxx, 'T',  1., psyy, 'T',  1.   &
-         &              , psxy, 'T',  1. )
-
-      IF(ln_ctl) THEN
-         CALL prt_ctl(tab2d_1=psm  , clinfo1=' adv_x: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
-         CALL prt_ctl(tab2d_1=psx  , clinfo1=' adv_x: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
-         CALL prt_ctl(tab2d_1=psy  , clinfo1=' adv_x: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
-         CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_x: psxy :')
-      ENDIF
+      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   &
+         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
+         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. )
       !
    END SUBROUTINE adv_x
 
 
-   SUBROUTINE adv_y( pdf, pvt , pcrh, psm , ps0 ,   &
+   SUBROUTINE adv_y( pdt, pvt , pcrh, psm , ps0 ,   &
       &              psx, psxx, psy , psyy, psxy )
       !!---------------------------------------------------------------------
@@ -448,14 +427,14 @@
       !!                variable on y axis
       !!---------------------------------------------------------------------
-      REAL(wp)                    , INTENT(in   ) ::   pdf                ! reduction factor for the time step
-      REAL(wp)                    , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psm                ! area
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ps0                ! field to be advected
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psx , psy          ! 1st moments 
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
+      REAL(wp)                  , INTENT(in   ) ::   pdt                ! time step
+      REAL(wp)                  , INTENT(in   ) ::   pcrh               ! call adv_x then adv_y (=1) or the opposite (=0)
+      REAL(wp), DIMENSION(:,:)  , INTENT(in   ) ::   pvt                ! j-direction ice velocity at V-point [m/s]
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psm                ! area
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ps0                ! field to be advected
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psx , psy          ! 1st moments 
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   psxx, psyy, psxy   ! 2nd moments
       !!
-      INTEGER  ::   ji, jj                               ! dummy loop indices
-      REAL(wp) ::   zs1max, zrdt, zslpmax, ztemp         ! temporary scalars
+      INTEGER  ::   ji, jj, jl, jcat                     ! dummy loop indices
+      REAL(wp) ::   zs1max, zslpmax, ztemp               ! temporary scalars
       REAL(wp) ::   zs1new, zalf , zalfq , zbt           !    -         -
       REAL(wp) ::   zs2new, zalf1, zalf1q, zbt1          !    -         -
@@ -464,152 +443,203 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   zalg, zalg1, zalg1q     !  -      -
       !---------------------------------------------------------------------
-
-      ! Limitation of moments.
-
-      zrdt = rdt_ice * pdf ! If ice drift field is too fast, use an appropriate time step for advection.
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zslpmax = MAX( 0._wp, ps0(ji,jj) )
-            zs1max  = 1.5 * zslpmax
-            zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj) ) )
-            zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
-               &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
-            rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
-            !
-            ps0 (ji,jj) = zslpmax  
-            psx (ji,jj) = psx (ji,jj) * rswitch
-            psxx(ji,jj) = psxx(ji,jj) * rswitch
-            psy (ji,jj) = zs1new * rswitch
-            psyy(ji,jj) = zs2new * rswitch
-            psxy(ji,jj) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj) ) ) * rswitch
-         END DO
+      !
+      jcat = SIZE( ps0 , 3 )   ! size of input arrays
+      !      
+      DO jl = 1, jcat   ! loop on categories
+         !
+         ! Limitation of moments.
+         DO jj = 1, jpj
+            DO ji = fs_2, fs_jpim1
+               !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
+               psm(ji,jj,jl) = MAX(  pcrh * e1e2t(ji,jj) + ( 1.0 - pcrh ) * psm(ji,jj,jl) , epsi20  )
+               !
+               zslpmax = MAX( 0._wp, ps0(ji,jj,jl) )
+               zs1max  = 1.5 * zslpmax
+               zs1new  = MIN( zs1max, MAX( -zs1max, psy(ji,jj,jl) ) )
+               zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
+                  &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj,jl) )  )
+               rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zslpmax) ) ) * tmask(ji,jj,1)   ! Case of empty boxes & Apply mask
+               !
+               ps0 (ji,jj,jl) = zslpmax  
+               psx (ji,jj,jl) = psx (ji,jj,jl) * rswitch
+               psxx(ji,jj,jl) = psxx(ji,jj,jl) * rswitch
+               psy (ji,jj,jl) = zs1new         * rswitch
+               psyy(ji,jj,jl) = zs2new         * rswitch
+               psxy(ji,jj,jl) = MIN( zslpmax, MAX( -zslpmax, psxy(ji,jj,jl) ) ) * rswitch
+            END DO
+         END DO
+ 
+         !  Calculate fluxes and moments between boxes j<-->j+1              
+         DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
+            DO ji = fs_2, fs_jpim1
+               zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
+               zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl)
+               zalfq        =  zalf * zalf
+               zalf1        =  1.0 - zalf
+               zalf1q       =  zalf1 * zalf1
+               !
+               zfm (ji,jj)  =  zalf  * psm(ji,jj,jl)
+               zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj,jl) + zalf1 * ( psy(ji,jj,jl)  + (zalf1-zalf) * psyy(ji,jj,jl) ) ) 
+               zfy (ji,jj)  =  zalfq *( psy(ji,jj,jl) + 3.0*zalf1*psyy(ji,jj,jl) )
+               zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj,jl)
+               zfx (ji,jj)  =  zalf  * ( psx(ji,jj,jl) + zalf1 * psxy(ji,jj,jl) )
+               zfxy(ji,jj)  =  zalfq * psxy(ji,jj,jl)
+               zfxx(ji,jj)  =  zalf  * psxx(ji,jj,jl)
+               !
+               !  Readjust moments remaining in the box.
+               psm (ji,jj,jl)  =  psm (ji,jj,jl) - zfm(ji,jj)
+               ps0 (ji,jj,jl)  =  ps0 (ji,jj,jl) - zf0(ji,jj)
+               psy (ji,jj,jl)  =  zalf1q * ( psy(ji,jj,jl) -3.0 * zalf * psyy(ji,jj,jl) )
+               psyy(ji,jj,jl)  =  zalf1 * zalf1q * psyy(ji,jj,jl)
+               psx (ji,jj,jl)  =  psx (ji,jj,jl) - zfx(ji,jj)
+               psxx(ji,jj,jl)  =  psxx(ji,jj,jl) - zfxx(ji,jj)
+               psxy(ji,jj,jl)  =  zalf1q * psxy(ji,jj,jl)
+            END DO
+         END DO
+         !
+         DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
+            DO ji = fs_2, fs_jpim1
+               zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl) 
+               zalg  (ji,jj) = zalf
+               zalfq         = zalf * zalf
+               zalf1         = 1.0 - zalf
+               zalg1 (ji,jj) = zalf1
+               zalf1q        = zalf1 * zalf1
+               zalg1q(ji,jj) = zalf1q
+               !
+               zfm   (ji,jj) = zfm (ji,jj) + zalf  *    psm (ji,jj+1,jl)
+               zf0   (ji,jj) = zf0 (ji,jj) + zalf  * (  ps0 (ji,jj+1,jl) &
+                  &                                   - zalf1 * (psy(ji,jj+1,jl) - (zalf1 - zalf ) * psyy(ji,jj+1,jl) ) )
+               zfy   (ji,jj) = zfy (ji,jj) + zalfq * (  psy (ji,jj+1,jl) - 3.0 * zalf1 * psyy(ji,jj+1,jl) )
+               zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *    psyy(ji,jj+1,jl) * zalfq
+               zfx   (ji,jj) = zfx (ji,jj) + zalf  * (  psx (ji,jj+1,jl) - zalf1 * psxy(ji,jj+1,jl) )
+               zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *    psxy(ji,jj+1,jl)
+               zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *    psxx(ji,jj+1,jl)
+            END DO
+         END DO
+
+         !  Readjust moments remaining in the box. 
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               zbt  =         zbet(ji,jj-1)
+               zbt1 = ( 1.0 - zbet(ji,jj-1) )
+               !
+               psm (ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) - zfm(ji,jj-1) )
+               ps0 (ji,jj,jl) = zbt * ps0(ji,jj,jl) + zbt1 * ( ps0(ji,jj,jl) - zf0(ji,jj-1) )
+               psy (ji,jj,jl) = zalg1q(ji,jj-1) * ( psy(ji,jj,jl) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj,jl) )
+               psyy(ji,jj,jl) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj,jl)
+               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) - zfx (ji,jj-1) )
+               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) - zfxx(ji,jj-1) )
+               psxy(ji,jj,jl) = zalg1q(ji,jj-1) * psxy(ji,jj,jl)
+            END DO
+         END DO
+
+         !   Put the temporary moments into appropriate neighboring boxes.    
+         DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
+            DO ji = fs_2, fs_jpim1
+               zbt  =       zbet(ji,jj-1)
+               zbt1 = 1.0 - zbet(ji,jj-1)
+               psm(ji,jj,jl) = zbt * ( psm(ji,jj,jl) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj,jl) 
+               zalf          = zbt * zfm(ji,jj-1) / psm(ji,jj,jl) 
+               zalf1         = 1.0 - zalf
+               ztemp         = zalf * ps0(ji,jj,jl) - zalf1 * zf0(ji,jj-1)
+               !
+               ps0(ji,jj,jl)  =   zbt  * ( ps0(ji,jj,jl) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj,jl)
+               psy(ji,jj,jl)  =   zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )  &
+                  &             + zbt1 * psy(ji,jj,jl)  
+               psyy(ji,jj,jl) =   zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj,jl)                           &
+                  &                      + 5.0 * ( zalf * zalf1 * ( psy(ji,jj,jl) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) ) & 
+                  &             + zbt1 * psyy(ji,jj,jl)
+               psxy(ji,jj,jl) =   zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj,jl)            &
+                  &                      + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj,jl) ) )  &
+                  &             + zbt1 * psxy(ji,jj,jl)
+               psx (ji,jj,jl) =   zbt * ( psx (ji,jj,jl) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj,jl)
+               psxx(ji,jj,jl) =   zbt * ( psxx(ji,jj,jl) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj,jl)
+            END DO
+         END DO
+
+         DO jj = 2, jpjm1                      !  Flux from j+1 to j IF v LT 0.
+            DO ji = fs_2, fs_jpim1
+               zbt  =       zbet(ji,jj)
+               zbt1 = 1.0 - zbet(ji,jj)
+               psm(ji,jj,jl) = zbt * psm(ji,jj,jl) + zbt1 * ( psm(ji,jj,jl) + zfm(ji,jj) )
+               zalf          = zbt1 * zfm(ji,jj) / psm(ji,jj,jl)
+               zalf1         = 1.0 - zalf
+               ztemp         = - zalf * ps0(ji,jj,jl) + zalf1 * zf0(ji,jj)
+               !
+               ps0 (ji,jj,jl) = zbt * ps0 (ji,jj,jl) + zbt1 * (  ps0(ji,jj,jl) + zf0(ji,jj) )
+               psy (ji,jj,jl) = zbt * psy (ji,jj,jl) + zbt1 * (  zalf * zfy(ji,jj) + zalf1 * psy(ji,jj,jl) + 3.0 * ztemp )
+               psyy(ji,jj,jl) = zbt * psyy(ji,jj,jl) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj,jl) &
+                  &                                            + 5.0 * ( zalf * zalf1 * ( - psy(ji,jj,jl) + zfy(ji,jj) )    &
+                  &                                            + ( zalf1 - zalf ) * ztemp ) )
+               psxy(ji,jj,jl) = zbt * psxy(ji,jj,jl) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj,jl)  &
+                  &                                            + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj,jl) ) )
+               psx (ji,jj,jl) = zbt * psx (ji,jj,jl) + zbt1 * ( psx (ji,jj,jl) + zfx (ji,jj) )
+               psxx(ji,jj,jl) = zbt * psxx(ji,jj,jl) + zbt1 * ( psxx(ji,jj,jl) + zfxx(ji,jj) )
+            END DO
+         END DO
+
       END DO
 
-      !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
-      psm(:,:)  = MAX(  pcrh * e1e2t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
-
-      !  Calculate fluxes and moments between boxes j<-->j+1              
-      DO jj = 1, jpj                     !  Flux from j to j+1 WHEN v GT 0   
-         DO ji = 1, jpi
-            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) )
-            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * zrdt * e1v(ji,jj) / psm(ji,jj)
-            zalfq        =  zalf * zalf
-            zalf1        =  1.0 - zalf
-            zalf1q       =  zalf1 * zalf1
-            !
-            zfm (ji,jj)  =  zalf  * psm(ji,jj)
-            zf0 (ji,jj)  =  zalf  * ( ps0(ji,jj) + zalf1 * ( psy(ji,jj)  + (zalf1-zalf) * psyy(ji,jj)  ) ) 
-            zfy (ji,jj)  =  zalfq *( psy(ji,jj) + 3.0*zalf1*psyy(ji,jj) )
-            zfyy(ji,jj)  =  zalf  * zalfq * psyy(ji,jj)
-            zfx (ji,jj)  =  zalf  * ( psx(ji,jj) + zalf1 * psxy(ji,jj) )
-            zfxy(ji,jj)  =  zalfq * psxy(ji,jj)
-            zfxx(ji,jj)  =  zalf  * psxx(ji,jj)
-            !
-            !  Readjust moments remaining in the box.
-            psm (ji,jj)  =  psm (ji,jj) - zfm(ji,jj)
-            ps0 (ji,jj)  =  ps0 (ji,jj) - zf0(ji,jj)
-            psy (ji,jj)  =  zalf1q * ( psy(ji,jj) -3.0 * zalf * psyy(ji,jj) )
-            psyy(ji,jj)  =  zalf1 * zalf1q * psyy(ji,jj)
-            psx (ji,jj)  =  psx (ji,jj) - zfx(ji,jj)
-            psxx(ji,jj)  =  psxx(ji,jj) - zfxx(ji,jj)
-            psxy(ji,jj)  =  zalf1q * psxy(ji,jj)
+      !-- Lateral boundary conditions
+      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm(:,:,1:jcat) , 'T',  1., ps0 , 'T',  1.   &
+         &                                , psx             , 'T', -1., psy , 'T', -1.   &   ! caution gradient ==> the sign changes
+         &                                , psxx            , 'T',  1., psyy, 'T',  1. , psxy, 'T',  1. )
+      !
+   END SUBROUTINE adv_y
+
+
+   SUBROUTINE Hsnow( pdt, pv_i, pv_s, pa_i, pa_ip, pe_s )
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE Hsnow  ***
+      !!
+      !! ** Purpose : 1- Check snow load after advection
+      !!              2- Correct pond concentration to avoid a_ip > a_i
+      !!
+      !! ** Method :  If snow load makes snow-ice interface to deplet below the ocean surface
+      !!              then put the snow excess in the ocean
+      !!
+      !! ** Notes :   This correction is crucial because of the call to routine icecor afterwards
+      !!              which imposes a mini of ice thick. (rn_himin). This imposed mini can artificially
+      !!              make the snow very thick (if concentration decreases drastically)
+      !!              This behavior has been seen in Ultimate-Macho and supposedly it can also be true for Prather
+      !!-------------------------------------------------------------------
+      REAL(wp)                    , INTENT(in   ) ::   pdt   ! tracer time-step
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip
+      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s
+      !
+      INTEGER  ::   ji, jj, jl   ! dummy loop indices
+      REAL(wp) ::   z1_dt, zvs_excess, zfra
+      !!-------------------------------------------------------------------
+      !
+      z1_dt = 1._wp / pdt
+      !
+      ! -- check snow load -- !
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( pv_i(ji,jj,jl) > 0._wp ) THEN
+                  !
+                  zvs_excess = MAX( 0._wp, pv_s(ji,jj,jl) - pv_i(ji,jj,jl) * (rau0-rhoi) * r1_rhos )
+                  !
+                  IF( zvs_excess > 0._wp ) THEN   ! snow-ice interface deplets below the ocean surface
+                     ! put snow excess in the ocean
+                     zfra = ( pv_s(ji,jj,jl) - zvs_excess ) / MAX( pv_s(ji,jj,jl), epsi20 )
+                     wfx_res(ji,jj) = wfx_res(ji,jj) + zvs_excess * rhos * z1_dt
+                     hfx_res(ji,jj) = hfx_res(ji,jj) - SUM( pe_s(ji,jj,1:nlay_s,jl) ) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0
+                     ! correct snow volume and heat content
+                     pe_s(ji,jj,1:nlay_s,jl) = pe_s(ji,jj,1:nlay_s,jl) * zfra
+                     pv_s(ji,jj,jl)          = pv_s(ji,jj,jl) - zvs_excess
+                  ENDIF
+                  !
+               ENDIF
+            END DO
          END DO
       END DO
       !
-      DO jj = 1, jpjm1                   !  Flux from j+1 to j when v LT 0.
-         DO ji = 1, jpi
-            zalf          = ( MAX(0._wp, -pvt(ji,jj) ) * zrdt * e1v(ji,jj) ) / psm(ji,jj+1) 
-            zalg  (ji,jj) = zalf
-            zalfq         = zalf * zalf
-            zalf1         = 1.0 - zalf
-            zalg1 (ji,jj) = zalf1
-            zalf1q        = zalf1 * zalf1
-            zalg1q(ji,jj) = zalf1q
-            !
-            zfm   (ji,jj) = zfm (ji,jj) + zalf  *   psm (ji,jj+1)
-            zf0   (ji,jj) = zf0 (ji,jj) + zalf  * ( ps0 (ji,jj+1) - zalf1 * (psy(ji,jj+1) - (zalf1 - zalf ) * psyy(ji,jj+1) ) )
-            zfy   (ji,jj) = zfy (ji,jj) + zalfq * ( psy (ji,jj+1) - 3.0 * zalf1 * psyy(ji,jj+1) )
-            zfyy  (ji,jj) = zfyy(ji,jj) + zalf  *   psyy(ji,jj+1) * zalfq
-            zfx   (ji,jj) = zfx (ji,jj) + zalf  * ( psx (ji,jj+1) - zalf1 * psxy(ji,jj+1) )
-            zfxy  (ji,jj) = zfxy(ji,jj) + zalfq *   psxy(ji,jj+1)
-            zfxx  (ji,jj) = zfxx(ji,jj) + zalf  *   psxx(ji,jj+1)
-         END DO
-      END DO
-
-      !  Readjust moments remaining in the box. 
-      DO jj = 2, jpj
-         DO ji = 1, jpi
-            zbt  =         zbet(ji,jj-1)
-            zbt1 = ( 1.0 - zbet(ji,jj-1) )
-            !
-            psm (ji,jj) = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) - zfm(ji,jj-1) )
-            ps0 (ji,jj) = zbt * ps0(ji,jj) + zbt1 * ( ps0(ji,jj) - zf0(ji,jj-1) )
-            psy (ji,jj) = zalg1q(ji,jj-1) * ( psy(ji,jj) + 3.0 * zalg(ji,jj-1) * psyy(ji,jj) )
-            psyy(ji,jj) = zalg1 (ji,jj-1) * zalg1q(ji,jj-1) * psyy(ji,jj)
-            psx (ji,jj) = zbt * psx (ji,jj) + zbt1 * ( psx (ji,jj) - zfx (ji,jj-1) )
-            psxx(ji,jj) = zbt * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) - zfxx(ji,jj-1) )
-            psxy(ji,jj) = zalg1q(ji,jj-1) * psxy(ji,jj)
-         END DO
-      END DO
-
-      !   Put the temporary moments into appropriate neighboring boxes.    
-      DO jj = 2, jpjm1                    !   Flux from j to j+1 IF v GT 0.
-         DO ji = 1, jpi
-            zbt  =         zbet(ji,jj-1)
-            zbt1 = ( 1.0 - zbet(ji,jj-1) )
-            psm(ji,jj)  = zbt * ( psm(ji,jj) + zfm(ji,jj-1) ) + zbt1 * psm(ji,jj) 
-            zalf        = zbt * zfm(ji,jj-1) / psm(ji,jj) 
-            zalf1       = 1.0 - zalf
-            ztemp       = zalf * ps0(ji,jj) - zalf1 * zf0(ji,jj-1)
-            !
-            ps0(ji,jj)  = zbt  * ( ps0(ji,jj) + zf0(ji,jj-1) ) + zbt1 * ps0(ji,jj)
-            psy(ji,jj)  = zbt  * ( zalf * zfy(ji,jj-1) + zalf1 * psy(ji,jj) + 3.0 * ztemp )   &
-               &                                               + zbt1 * psy(ji,jj)  
-            psyy(ji,jj) = zbt  * ( zalf * zalf * zfyy(ji,jj-1) + zalf1 * zalf1 * psyy(ji,jj)                             &
-               &                 + 5.0 * ( zalf * zalf1 * ( psy(ji,jj) - zfy(ji,jj-1) ) - ( zalf1 - zalf ) * ztemp ) )   & 
-               &                                               + zbt1 * psyy(ji,jj)
-            psxy(ji,jj) = zbt  * (  zalf * zfxy(ji,jj-1) + zalf1 * psxy(ji,jj)               &
-               &                  + 3.0 * (- zalf1 * zfx(ji,jj-1) + zalf * psx(ji,jj) )  )   &
-               &                                                + zbt1 * psxy(ji,jj)
-            psx (ji,jj) = zbt * ( psx (ji,jj) + zfx (ji,jj-1) ) + zbt1 * psx (ji,jj)
-            psxx(ji,jj) = zbt * ( psxx(ji,jj) + zfxx(ji,jj-1) ) + zbt1 * psxx(ji,jj)
-         END DO
-      END DO
-
-      DO jj = 2, jpjm1                   !  Flux from j+1 to j IF v LT 0.
-         DO ji = 1, jpi
-            zbt  =         zbet(ji,jj)
-            zbt1 = ( 1.0 - zbet(ji,jj) )
-            psm(ji,jj)  = zbt * psm(ji,jj) + zbt1 * ( psm(ji,jj) + zfm(ji,jj) )
-            zalf        = zbt1 * zfm(ji,jj) / psm(ji,jj)
-            zalf1       = 1.0 - zalf
-            ztemp       = - zalf * ps0 (ji,jj) + zalf1 * zf0(ji,jj)
-            ps0 (ji,jj) =   zbt  * ps0 (ji,jj) + zbt1  * ( ps0(ji,jj) + zf0(ji,jj) )
-            psy (ji,jj) =   zbt  * psy (ji,jj) + zbt1  * ( zalf * zfy(ji,jj) + zalf1 * psy(ji,jj) + 3.0 * ztemp )
-            psyy(ji,jj) =   zbt  * psyy(ji,jj) + zbt1 * (  zalf * zalf * zfyy(ji,jj) + zalf1 * zalf1 * psyy(ji,jj)   &
-               &                                         + 5.0 *( zalf *zalf1 *( -psy(ji,jj) + zfy(ji,jj) )          &
-               &                                         + ( zalf1 - zalf ) * ztemp )                                )
-            psxy(ji,jj) =   zbt  * psxy(ji,jj) + zbt1 * (  zalf * zfxy(ji,jj) + zalf1 * psxy(ji,jj)       &
-               &                                         + 3.0 * ( zalf1 * zfx(ji,jj) - zalf * psx(ji,jj) )  )
-            psx (ji,jj) =   zbt  * psx (ji,jj) + zbt1 * ( psx (ji,jj) + zfx (ji,jj) )
-            psxx(ji,jj) =   zbt  * psxx(ji,jj) + zbt1 * ( psxx(ji,jj) + zfxx(ji,jj) )
-         END DO
-      END DO
-
-      !-- Lateral boundary conditions
-      CALL lbc_lnk_multi( 'icedyn_adv_pra', psm , 'T',  1.,  ps0 , 'T',  1.   &
-         &              , psx , 'T', -1.,  psy , 'T', -1.   &   ! caution gradient ==> the sign changes
-         &              , psxx, 'T',  1.,  psyy, 'T',  1.   &
-         &              , psxy, 'T',  1. )
-
-      IF(ln_ctl) THEN
-         CALL prt_ctl(tab2d_1=psm  , clinfo1=' adv_y: psm  :', tab2d_2=ps0 , clinfo2=' ps0  : ')
-         CALL prt_ctl(tab2d_1=psx  , clinfo1=' adv_y: psx  :', tab2d_2=psxx, clinfo2=' psxx : ')
-         CALL prt_ctl(tab2d_1=psy  , clinfo1=' adv_y: psy  :', tab2d_2=psyy, clinfo2=' psyy : ')
-         CALL prt_ctl(tab2d_1=psxy , clinfo1=' adv_y: psxy :')
-      ENDIF
-      !
-   END SUBROUTINE adv_y
+      !-- correct pond concentration to avoid a_ip > a_i -- !
+      WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:)
+      !
+   END SUBROUTINE Hsnow
 
 
@@ -624,6 +654,5 @@
       !
       !                             !* allocate prather fields
-      ALLOCATE( sxopw(jpi,jpj)     , syopw(jpi,jpj)     , sxxopw(jpi,jpj)     , syyopw(jpi,jpj)     , sxyopw(jpi,jpj)     ,   &
-         &      sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
+      ALLOCATE( sxice(jpi,jpj,jpl) , syice(jpi,jpj,jpl) , sxxice(jpi,jpj,jpl) , syyice(jpi,jpj,jpl) , sxyice(jpi,jpj,jpl) ,   &
          &      sxsn (jpi,jpj,jpl) , sysn (jpi,jpj,jpl) , sxxsn (jpi,jpj,jpl) , syysn (jpi,jpj,jpl) , sxysn (jpi,jpj,jpl) ,   &
          &      sxa  (jpi,jpj,jpl) , sya  (jpi,jpj,jpl) , sxxa  (jpi,jpj,jpl) , syya  (jpi,jpj,jpl) , sxya  (jpi,jpj,jpl) ,   &
@@ -652,5 +681,5 @@
       !!                   ***  ROUTINE adv_pra_rst  ***
       !!                     
-      !! ** Purpose :   Read or write RHG file in restart file
+      !! ** Purpose :   Read or write file in restart file
       !!
       !! ** Method  :   use of IOM library
@@ -671,5 +700,5 @@
          !                                   !==========================!
          !
-         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxopw' , ldstop = .FALSE. )    ! file exist: id1>0
+         IF( ln_rstart ) THEN   ;   id1 = iom_varid( numrir, 'sxice' , ldstop = .FALSE. )    ! file exist: id1>0
          ELSE                   ;   id1 = 0                                                  ! no restart: id1=0
          ENDIF
@@ -689,5 +718,5 @@
             CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  )
             CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  )
-            !                                                        ! lead fraction
+            !                                                        ! ice concentration
             CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    )
             CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    )
@@ -707,10 +736,4 @@
             CALL iom_get( numrir, jpdom_autoglo, 'syyage', syyage )
             CALL iom_get( numrir, jpdom_autoglo, 'sxyage', sxyage )
-            !                                                        ! open water in sea ice
-            CALL iom_get( numrir, jpdom_autoglo, 'sxopw' , sxopw  )
-            CALL iom_get( numrir, jpdom_autoglo, 'syopw' , syopw  )
-            CALL iom_get( numrir, jpdom_autoglo, 'sxxopw', sxxopw )
-            CALL iom_get( numrir, jpdom_autoglo, 'syyopw', syyopw )
-            CALL iom_get( numrir, jpdom_autoglo, 'sxyopw', sxyopw )
             !                                                        ! snow layers heat content
             DO jk = 1, nlay_s
@@ -752,8 +775,7 @@
             sxice = 0._wp   ;   syice = 0._wp   ;   sxxice = 0._wp   ;   syyice = 0._wp   ;   sxyice = 0._wp      ! ice thickness
             sxsn  = 0._wp   ;   sysn  = 0._wp   ;   sxxsn  = 0._wp   ;   syysn  = 0._wp   ;   sxysn  = 0._wp      ! snow thickness
-            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! lead fraction
+            sxa   = 0._wp   ;   sya   = 0._wp   ;   sxxa   = 0._wp   ;   syya   = 0._wp   ;   sxya   = 0._wp      ! ice concentration
             sxsal = 0._wp   ;   sysal = 0._wp   ;   sxxsal = 0._wp   ;   syysal = 0._wp   ;   sxysal = 0._wp      ! ice salinity
             sxage = 0._wp   ;   syage = 0._wp   ;   sxxage = 0._wp   ;   syyage = 0._wp   ;   sxyage = 0._wp      ! ice age
-            sxopw = 0._wp   ;   syopw = 0._wp   ;   sxxopw = 0._wp   ;   syyopw = 0._wp   ;   sxyopw = 0._wp      ! open water in sea ice
             sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content
             sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content
@@ -786,5 +808,5 @@
          CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn  )
          CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn  )
-         !                                                           ! lead fraction
+         !                                                           ! ice concentration
          CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa    )
          CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya    )
@@ -804,10 +826,4 @@
          CALL iom_rstput( iter, nitrst, numriw, 'syyage', syyage )
          CALL iom_rstput( iter, nitrst, numriw, 'sxyage', sxyage )
-         !                                                           ! open water in sea ice
-         CALL iom_rstput( iter, nitrst, numriw, 'sxopw' , sxopw  )
-         CALL iom_rstput( iter, nitrst, numriw, 'syopw' , syopw  )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxxopw', sxxopw )
-         CALL iom_rstput( iter, nitrst, numriw, 'syyopw', syyopw )
-         CALL iom_rstput( iter, nitrst, numriw, 'sxyopw', sxyopw )
          !                                                           ! snow layers heat content
          DO jk = 1, nlay_s
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_adv_umx.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_adv_umx.F90	(revision 12165)
@@ -83,5 +83,5 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   poa_i      ! age content
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_i       ! ice concentration
-      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond concentration
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
@@ -319,13 +319,11 @@
          !
          !== Ice age ==!
-         IF( iom_use('iceage') .OR. iom_use('iceage_cat') ) THEN
-            zamsk = 1._wp
-            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
-               &                                      poa_i, poa_i )
-         ENDIF
+         zamsk = 1._wp
+         CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat, zv_cat, zcu_box, zcv_box, &
+            &                                      poa_i, poa_i )
          !
          !== melt ponds ==!
          IF ( ln_pnd_H12 ) THEN
-            ! fraction
+            ! concentration
             zamsk = 1._wp
             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zu_cat , zv_cat , zcu_box, zcv_box, &
@@ -1529,5 +1527,5 @@
       !!              3- check whether snow load deplets the snow-ice interface below sea level$
       !!                 and reduce it by sending the excess in the ocean
-      !!              4- correct pond fraction to avoid a_ip > a_i
+      !!              4- correct pond concentration to avoid a_ip > a_i
       !!
       !! ** input   : Max thickness of the surrounding 9-points
@@ -1599,5 +1597,5 @@
          END DO
       END DO 
-      !                                           !-- correct pond fraction to avoid a_ip > a_i
+      !                                           !-- correct pond concentration to avoid a_ip > a_i
       WHERE( pa_ip(:,:,:) > pa_i(:,:,:) )   pa_ip(:,:,:) = pa_i(:,:,:)
       !
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_rdgrft.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/icedyn_rdgrft.F90	(revision 12165)
@@ -86,7 +86,7 @@
       !!                ***  ROUTINE ice_dyn_rdgrft_alloc ***
       !!-------------------------------------------------------------------
-      ALLOCATE( closing_net(jpij), opning(jpij)   , closing_gross(jpij),   &
-         &      apartf(jpij,0:jpl), hrmin(jpij,jpl), hraft(jpij,jpl)    , aridge(jpij,jpl),   &
-         &      hrmax(jpij,jpl), hi_hrdg(jpij,jpl)  , araft (jpij,jpl),  &
+      ALLOCATE( closing_net(jpij)  , opning(jpij)      , closing_gross(jpij) ,               &
+         &      apartf(jpij,0:jpl) , hrmin  (jpij,jpl) , hraft(jpij,jpl) , aridge(jpij,jpl), &
+         &      hrmax (jpij,jpl)   , hi_hrdg(jpij,jpl) , araft(jpij,jpl) ,                   &
          &      ze_i_2d(jpij,nlay_i,jpl), ze_s_2d(jpij,nlay_s,jpl), STAT=ice_dyn_rdgrft_alloc )
 
@@ -137,5 +137,4 @@
       REAL(wp) ::   zfac                       ! local scalar
       INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not
-      REAL(wp), DIMENSION(jpij) ::   zdivu_adv     ! divu as implied by transport scheme  (1/s)
       REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i
       !
@@ -175,7 +174,7 @@
         
          ! just needed here
-         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i  )
          CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i )
          ! needed here and in the iteration loop
+         CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i) ! zdivu is used as a work array here (no change in divu_i)
          CALL tab_3d_2d( npti, nptidx(1:npti), a_i_2d  (1:npti,1:jpl), a_i   )
          CALL tab_3d_2d( npti, nptidx(1:npti), v_i_2d  (1:npti,1:jpl), v_i   )
@@ -187,16 +186,8 @@
             closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
             !
-            ! divergence given by the advection scheme
-            !   (which may not be equal to divu as computed from the velocity field)
-            IF    ( ln_adv_Pra ) THEN
-               zdivu_adv(ji) = ( 1._wp - ato_i_1d(ji) - SUM( a_i_2d(ji,:) ) ) * r1_rdtice
-            ELSEIF( ln_adv_UMx ) THEN
-               zdivu_adv(ji) = zdivu(ji)
-            ENDIF
-            !
-            IF( zdivu_adv(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu_adv(ji) )   ! make sure the closing rate is large enough
-            !                                                                                        ! to give asum = 1.0 after ridging
+            IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough
+            !                                                                                ! to give asum = 1.0 after ridging
             ! Opening rate (non-negative) that will give asum = 1.0 after ridging.
-            opning(ji) = closing_net(ji) + zdivu_adv(ji)
+            opning(ji) = closing_net(ji) + zdivu(ji)
          END DO
          !
@@ -215,5 +206,5 @@
                ato_i_1d   (ipti)   = ato_i_1d   (ji)
                closing_net(ipti)   = closing_net(ji)
-               zdivu_adv  (ipti)   = zdivu_adv  (ji)
+               zdivu      (ipti)   = zdivu      (ji)
                opning     (ipti)   = opning     (ji)
             ENDIF
@@ -259,7 +250,7 @@
                ELSE
                   iterate_ridging  = 1
-                  zdivu_adv  (ji) = zfac * r1_rdtice
-                  closing_net(ji) = MAX( 0._wp, -zdivu_adv(ji) )
-                  opning     (ji) = MAX( 0._wp,  zdivu_adv(ji) )
+                  zdivu      (ji) = zfac * r1_rdtice
+                  closing_net(ji) = MAX( 0._wp, -zdivu(ji) )
+                  opning     (ji) = MAX( 0._wp,  zdivu(ji) )
                ENDIF
             END DO
@@ -309,5 +300,5 @@
 
       !                       ! Ice thickness needed for rafting
-      WHERE( pa_i(1:npti,:) > epsi20 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
+      WHERE( pa_i(1:npti,:) > epsi10 )   ;   zhi(1:npti,:) = pv_i(1:npti,:) / pa_i(1:npti,:)
       ELSEWHERE                          ;   zhi(1:npti,:) = 0._wp
       END WHERE
@@ -328,5 +319,5 @@
       zasum(1:npti) = pato_i(1:npti) + SUM( pa_i(1:npti,:), dim=2 )
       !
-      WHERE( zasum(1:npti) > epsi20 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
+      WHERE( zasum(1:npti) > epsi10 )   ;   z1_asum(1:npti) = 1._wp / zasum(1:npti)
       ELSEWHERE                         ;   z1_asum(1:npti) = 0._wp
       END WHERE
@@ -454,5 +445,5 @@
       ! Based on the ITD of ridging and ridged ice, convert the net closing rate to a gross closing rate.  
       ! NOTE: 0 < aksum <= 1
-      WHERE( zaksum(1:npti) > epsi20 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
+      WHERE( zaksum(1:npti) > epsi10 )   ;   closing_gross(1:npti) = pclosing_net(1:npti) / zaksum(1:npti)
       ELSEWHERE                          ;   closing_gross(1:npti) = 0._wp
       END WHERE
@@ -537,5 +528,5 @@
             IF( apartf(ji,jl1) > 0._wp .AND. closing_gross(ji) > 0._wp ) THEN   ! only if ice is ridging
 
-               IF( a_i_2d(ji,jl1) > epsi20 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
+               IF( a_i_2d(ji,jl1) > epsi10 ) THEN   ;   z1_ai(ji) = 1._wp / a_i_2d(ji,jl1)
                ELSE                                 ;   z1_ai(ji) = 0._wp
                ENDIF
@@ -595,5 +586,5 @@
                ! virtual salt flux to keep salinity constant
                IF( nn_icesal /= 2 )  THEN
-                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )        ! ridge salinity = s_i
+                  sirdg2(ji)     = sirdg2(ji)     - vsw * ( sss_1d(ji) - s_i_1d(ji) )       ! ridge salinity = s_i
                   sfx_bri_1d(ji) = sfx_bri_1d(ji) + sss_1d(ji) * vsw * rhoi * r1_rdtice  &  ! put back sss_m into the ocean
                      &                            - s_i_1d(ji) * vsw * rhoi * r1_rdtice     ! and get  s_i  from the ocean 
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/iceitd.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/iceitd.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/ICE/iceitd.F90	(revision 12165)
@@ -211,5 +211,5 @@
                CALL itd_glinear( zhb0(1:npti)  , zhb1(1:npti)  , h_ib_1d(1:npti)  , a_i_1d(1:npti)  ,  &   ! in
                   &              g0  (1:npti,1), g1  (1:npti,1), hL     (1:npti,1), hR    (1:npti,1)   )   ! out
-                  !
+               !
                ! Area lost due to melting of thin ice
                DO ji = 1, npti
@@ -218,5 +218,5 @@
                      !
                      zdh0 =  h_i_1d(ji) - h_ib_1d(ji)                
-                     IF( zdh0 < 0.0 ) THEN      !remove area from category 1
+                     IF( zdh0 < 0.0 ) THEN      ! remove area from category 1
                         zdh0 = MIN( -zdh0, hi_max(1) )
                         !Integrate g(1) from 0 to dh0 to estimate area melted
@@ -226,8 +226,7 @@
                            zx1    = zetamax
                            zx2    = 0.5 * zetamax * zetamax 
-                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                        ! ice area removed
+                           zda0   = g1(ji,1) * zx2 + g0(ji,1) * zx1                ! ice area removed
                            zdamax = a_i_1d(ji) * (1.0 - h_i_1d(ji) / h_ib_1d(ji) ) ! Constrain new thickness <= h_i                
-                           zda0   = MIN( zda0, zdamax )                                                  ! ice area lost due to melting 
-                           !     of thin ice (zdamax > 0)
+                           zda0   = MIN( zda0, zdamax )                            ! ice area lost due to melting of thin ice (zdamax > 0)
                            ! Remove area, conserving volume
                            h_i_1d(ji) = h_i_1d(ji) * a_i_1d(ji) / ( a_i_1d(ji) - zda0 )
@@ -349,5 +348,5 @@
       DO ji = 1, npti
          !
-         IF( paice(ji) > epsi10  .AND. phice(ji) > 0._wp )  THEN
+         IF( paice(ji) > epsi10  .AND. phice(ji) > epsi10 )  THEN
             !
             ! Initialize hL and hR
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/CRS/README.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/CRS/README.rst	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/CRS/README.rst	(revision 12165)
@@ -2,4 +2,8 @@
 On line biogeochemistry coarsening
 **********************************
+
+.. todo::
+
+
 
 .. contents::
@@ -63,6 +67,6 @@
 	                           ! 1, MAX of KZ
 	                           ! 2, MIN of KZ
-	                           ! 3, 10^(MEAN(LOG(KZ)) 
-	                           ! 4, MEDIANE of KZ 
+	                           ! 3, 10^(MEAN(LOG(KZ))
+	                           ! 4, MEDIANE of KZ
 	   ln_crs_wn   = .false.   ! wn coarsened (T) or computed using horizontal divergence ( F )
 	                           !                           !
@@ -73,5 +77,5 @@
   the north-fold lateral boundary condition (ORCA025, ORCA12, ORCA36, ...).
 - ``nn_msh_crs = 1`` will activate the generation of the coarsened grid meshmask.
-- ``nn_crs_kz`` is the operator to coarsen the vertical mixing coefficient. 
+- ``nn_crs_kz`` is the operator to coarsen the vertical mixing coefficient.
 - ``ln_crs_wn``
 
@@ -80,5 +84,5 @@
   - when ``key_vvl`` is not activated,
 
-    - coarsened vertical velocities are computed using horizontal divergence (``ln_crs_wn = .false.``) 
+    - coarsened vertical velocities are computed using horizontal divergence (``ln_crs_wn = .false.``)
     - or coarsened vertical velocities are computed with an average operator (``ln_crs_wn = .true.``)
 - ``ln_crs_top = .true.``: should be activated to run BCG model in coarsened space;
@@ -97,5 +101,5 @@
 
 In the [attachment:iodef.xml iodef.xml]  file, a "nemo" context is defined and
-some variable defined in [attachment:file_def.xml file_def.xml] are writted on the ocean-dynamic grid.  
+some variable defined in [attachment:file_def.xml file_def.xml] are writted on the ocean-dynamic grid.
 To write variables on the coarsened grid, and in particular the passive tracers,
 a "nemo_crs" context should be defined in [attachment:iodef.xml iodef.xml] and
@@ -111,5 +115,5 @@
   interpolated `on-the-fly <http://forge.ipsl.jussieu.fr/nemo/wiki/Users/SetupNewConfiguration/Weight-creator>`_.
   Example of namelist for PISCES :
-  
+
 	.. code-block:: fortran
 
@@ -134,5 +138,5 @@
 	   	rn_trfac(14)  =   1.0e-06  !  -      -      -     -
 		   rn_trfac(23)  =   7.6e-06  !  -      -      -     -
-		
+
 	   	cn_dir        =  './'      !  root directory for the location of the data files
 
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/DYN/dynnxt.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/DYN/dynnxt.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/DYN/dynnxt.F90	(revision 12165)
@@ -175,4 +175,6 @@
       IF( neuler == 0 .AND. kt == nit000 ) THEN        !* Euler at first time-step: only swap
          DO jk = 1, jpkm1
+            ub(:,:,jk) = un(:,:,jk)                         ! ub <-- un
+            vb(:,:,jk) = vn(:,:,jk)
             un(:,:,jk) = ua(:,:,jk)                         ! un <-- ua
             vn(:,:,jk) = va(:,:,jk)
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/FLO/flodom.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/FLO/flodom.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/FLO/flodom.F90	(revision 12165)
@@ -433,5 +433,5 @@
       IF( ABS(dlx) > 1.0_wp ) dlx = 1.0_wp
       !
-      dld = ATAN(DSQRT( 1._wp * ( 1._wp-dlx )/( 1._wp+dlx ) )) * 222.24_wp / dls
+      dld = ATAN(SQRT( 1._wp * ( 1._wp-dlx )/( 1._wp+dlx ) )) * 222.24_wp / dls
       flo_dstnce = dld * 1000._wp
       !
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/FLO/flowri.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/FLO/flowri.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/FLO/flowri.F90	(revision 12165)
@@ -221,5 +221,5 @@
                clname=TRIM(clname)//".nc"
 
-               CALL fliocrfd( clname , (/ 'ntraj' , 't' /), (/ jpnfl , -1  /) , numflo )
+               CALL fliocrfd( clname , (/'ntraj' , '    t' /), (/ jpnfl , -1/) , numflo )
    
                CALL fliodefv( numflo, 'traj_lon'    , (/1,2/), v_t=flio_r8, long_name="Longitude"           , units="degrees_east"  )
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/LBC/mppini.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/LBC/mppini.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/LBC/mppini.F90	(revision 12165)
@@ -538,5 +538,5 @@
  9401    FORMAT('              '   ,20('   ',i3,'          ') )
  9402    FORMAT('       ',i3,' *  ',20(i3,'  x',i3,'   *   ') )
- 9404    FORMAT('           *  '   ,20('      ',i3,'   *   ') )
+ 9404    FORMAT('           *  '   ,20('     ' ,i4,'   *   ') )
       ENDIF
          
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/LDF/ldfdyn.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/LDF/ldfdyn.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/LDF/ldfdyn.F90	(revision 12165)
@@ -315,6 +315,6 @@
             DO jj = 1, jpj             ! Set local gridscale values
                DO ji = 1, jpi
-                  esqt(ji,jj) = ( e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 
-                  esqf(ji,jj) = ( e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2 
+                  esqt(ji,jj) = ( 2._wp * e1e2t(ji,jj) / ( e1t(ji,jj) + e2t(ji,jj) ) )**2 
+                  esqf(ji,jj) = ( 2._wp * e1e2f(ji,jj) / ( e1f(ji,jj) + e2f(ji,jj) ) )**2 
                END DO
             END DO
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/USR/README-setup-configs
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/USR/README-setup-configs	(revision 12164)
+++ 	(revision )
@@ -1,53 +1,0 @@
-==============================================================
-== HOW TO set up user’s configuration in NEMO version 4.0 : ==
-==============================================================
-This document describes how a configuration is set up in NEMO.
-
-In NEMO 4.0 two types of configurations can be build: global/regional and idealized.
-
-- create ./cfgs/your_config
-- case 1: regional or global 
-    => ln_read_cfg = .true. in ./cfgs/your_config/EXP00/namelist_cfg (i.e. read domain_cfg.nc)
-- case 2: idealized 
-    =>in ./cfgs/your_config/EXP00/namelist_cfg :
-      ln_read_cfg = .false. (need to define usrdef_hgr.F90 and usrdef_zgr.F90) 
-      ln_usr = .true.       (need to define usrdef_sbc.F90)    
-    => copy ./src/OCE/USR in ./cfgs/your_config/MY_SRC routines nedeed to set up domain, forcing fields and initial state
-    
-  Routines in NEMOGCM/NEMO/OCE/USR define analytically the domain, the initial state and the surface boundary conditions for the reference configuration GYRE.
-	  •	usrdef_hgr.F90 and usrdef_zgr.F90 : define horizontal and vertical grid
-	  •	usrdef_sbc.F90    : provides at each time-step the surface boundary condition, i.e. the momentum, heat and freshwater fluxes
-	  •	usrdef_istate.F90 : defines initialization of the dynamics and tracers
-	  •	usrdef_fmask.F90  : used just in ORCA CONFIGURATIONS for alteration of f-point land/ocean mask in some straits 
-	  •	usrdef_closea.F90 : used just in ORCA CONFIGURATIONS for specific treatments associated with closed seas
-	  •	usrdef_nam.F90    : set the domain characteristics of a user configuration
-
-  Other analytical examples are given in ./tests.
-
-===========================================
-== The domain_cfg.nc file must contain : ==
-===========================================
-   int ORCA  , ORCA_index            : configuration name, configuration resolution
-   int jpiglo, jpjglo , jpkglo       : global domain sizes 
-   int jperio                        : lateral global domain b.c. 
-   int ln_zco, ln_zps, ln_sco        : flags for z-coordinate, z-coordinate with partial steps and s-coordinate 
-   int ln_isfcav                     : flag  for ice shelf cavities 
-   double glamt, glamu, glamv, glamf : geographic position 
-   double gphit, gphiu, gphiv, gphif : geographic position
-   double iff, ff_f, ff_t            : Coriolis parameter (if not on the sphere)
-   double e1t, e1u, e1v, e1f         : horizontal scale factors
-   double e2t, e2u, e2v, e2f         : horizontal scale factors
-   double ie1e2u_v, e1e2u, e1e2v     : U and V surfaces (if grid size reduction in some straits)
-   double e3t_1d, e3w_1d             : reference vertical scale factors at T and W points
-   double e3t_0, e3u_0, e3v_0, e3f_0, e3w_0 : vertical scale factors 3D coordinate at T,U,V,F and W points
-   double e3uw_0,e3vw_0              : vertical scale factors 3D coordinate at UW and VW points
-   int bottom_level, top_level       : last wet T-points, 1st wet T-points (for ice shelf cavities)
-
-=======================================
-= HOW TO CREATE domain_cfg.nc file : ==
-=======================================
-Two options are possible :
-
- - users can create with his own tool the domain_cfg.nc with all mandatory fields.  
-
- - tool is available in ./tools/DOMAINcfg based on NEMO 3.6 version. See ./tools/DOMAINcfg/README.
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/USR/README.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/USR/README.rst	(revision 12165)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OCE/USR/README.rst	(revision 12165)
@@ -0,0 +1,280 @@
+******************************
+Setting up a new configuration
+******************************
+
+.. todo::
+
+
+
+.. contents::
+   :local:
+
+Starting from an existing configuration
+=======================================
+
+There are three options to build a new configuration from an existing one.
+
+Option 1: Duplicate an existing configuration
+---------------------------------------------
+
+The NEMO so-called Reference Configurations cover a number of major features for NEMO setup
+(global, regional, 1D, using embedded zoom with AGRIF...)
+
+One can create a new configuration by duplicating one of the reference configurations
+(``ORCA2_ICE_PISCES`` in the following example)
+
+.. code-block:: console
+
+   $ ./makenemo –n 'ORCA2_ICE_PISCES_MINE' -r 'ORCA2_ICE_PISCES' -m 'my_arch'
+
+Option 2: Duplicate with differences
+------------------------------------
+
+Create and compile a new configuration based on a reference configuration
+(``ORCA2_ICE_PISCES`` in the following example) but with different pre-processor options.
+For this either add ``add_key`` or ``del_key`` keys as required; e.g.
+
+.. code-block:: console
+
+   $ ./makenemo –n 'ORCA2_ICE_PISCES_MINE' -r 'ORCA2_ICE_PISCES' -m 'my_arch' del_key 'key_iomput' add_key 'key_diahth'
+
+Option 3: Use the SIREN tools to subset an existing model
+---------------------------------------------------------
+
+Define a regional configuration which is a {sub,super}-set of an existing configuration.
+
+This last option employs the SIREN software tools that are included in the standard distribution.
+The software is written in Fortran 95 and available in the :file:`./tools/SIREN` directory.
+SIREN allows you to create your own regional configuration embedded in a wider one.
+
+SIREN is a set of programs to create all the input files you need to
+run a NEMO regional configuration.
+
+:Demo:     Set of GLORYS files (GLObal ReanalYSis on the ORCA025 grid),
+           as well as examples of namelists are available `here`_.
+:Doc:      :forge:`chrome/site/doc/SIREN/html/index.html`
+:Support:  Any questions or comments regarding the use of SIREN should be posted in
+  :forge:`the corresponding forum <discussion/forum/2>`.
+
+.. _here: https://prodn.idris.fr/thredds/catalog/ipsl_public/rron463/catalog.html
+
+Creating a completely new configuration
+=======================================
+
+From NEMO version 4.0 there are two ways to build configurations from scratch.
+The appropriate method to use depends largely on the target configuration.
+Method 1 is for more complex/realistic global or regional configurations and
+method 2 is intended for simpler, idealised configurations whose
+domains and characteristics can be described in simple geometries and formulae.
+
+Option 1: Create and use a domain configuration file
+----------------------------------------------------
+
+This method is used by each of the reference configurations,
+so that downloading their input files linked to their description can help.
+Although starting from scratch,
+it is advisable to create the directory structure to house your new configuration by
+duplicating the closest reference configuration to your target application.
+For example, if your application requires both ocean ice and passive tracers,
+then use the ``ORCA2_ICE_PISCES`` as template,
+and execute following command to build your ``MY_NEW_CONFIG`` configuration:
+
+.. code-block:: sh
+
+   $ ./makenemo –n 'MY_NEW_CONFIG' -r 'ORCA2_ICE_PISCES' -m 'my_arch'
+
+where ``MY_NEW_CONFIG`` can be substituted with
+a suitably descriptive name for your new configuration.
+
+The purpose of this step is simply to create and populate the appropriate :file:`WORK`,
+:file:`MY_SRC` and :file:`EXP00` subdirectories for your new configuration.
+Other choices for the base reference configuration might be
+
+:GYRE:  If your target application is ocean-only
+:AMM12: If your target application is regional with open boundaries
+
+All the domain information for your new configuration will be contained within
+a netcdf file called :file:`domain_cfg.nc` which you will need to create and
+place in the :file:`./cfgs/MY_NEW_CONFIG/EXP00` sub-directory.
+Firstly though, ensure that your configuration is set to use such a file by checking that
+
+.. code-block:: fortran
+
+   ln_read_cfg = .true.
+
+in :file:`./cfgs/MY_NEW_CONFIG/EXP00/namelist_cfg`
+
+Create the :file:`domain_cfg.nc` file which must contain the following fields
+
+.. code-block:: c
+
+   /* configuration name, configuration resolution                 */
+   int    ORCA, ORCA_index
+   /* global domain sizes                                          */
+   int    jpiglo, jpjglo, jpkglo
+   /* lateral global domain b.c.                                   */
+   int    jperio
+   /* flags for z-coord, z-coord with partial steps and s-coord    */
+   int    ln_zco, ln_zps, ln_sco
+   /* flag  for ice shelf cavities                                 */
+   int    ln_isfcav
+   /* geographic position                                          */
+   double glamt, glamu, glamv, glamf
+   /* geographic position                                          */
+   double gphit, gphiu, gphiv, gphif
+   /* Coriolis parameter (if not on the sphere)                    */
+   double iff, ff_f, ff_t
+   /* horizontal scale factors                                     */
+   double e1t, e1u, e1v, e1f
+   /* horizontal scale factors                                     */
+   double e2t, e2u, e2v, e2f
+   /* U and V surfaces (if grid size reduction in some straits)    */
+   double ie1e2u_v, e1e2u, e1e2v
+   /* reference vertical scale factors at T and W points           */
+   double e3t_1d, e3w_1d
+   /* vertical scale factors 3D coordinate at T,U,V,F and W points */
+   double e3t_0, e3u_0, e3v_0, e3f_0, e3w_0
+   /* vertical scale factors 3D coordinate at UW and VW points     */
+   double e3uw_0, e3vw_0
+   /* last wet T-points, 1st wet T-points (for ice shelf cavities) */
+   int    bottom_level, top_level
+
+There are two options for creating a :file:`domain_cfg.nc` file:
+
+- Users can use tools of their own choice to build a :file:`domain_cfg.nc` with all mandatory fields.
+- Users can adapt and apply the supplied tool available in :file:`./tools/DOMAINcfg`.
+  This tool is based on code extracted from NEMO version 3.6 and will allow similar choices for
+  the horizontal and vertical grids that were available internally to that version.
+  See :ref:`tools <DOMAINcfg>` for details.
+
+Option 2: Adapt the usr_def configuration module of NEMO for you own purposes
+-----------------------------------------------------------------------------
+
+This method is intended for configuring easily simple/idealised configurations which
+are often used as demonstrators or for process evaluation and comparison.
+This method can be used whenever the domain geometry has a simple mathematical description and
+the ocean initial state and boundary forcing is described analytically.
+As a start, consider the case of starting a completely new ocean-only test case based on
+the ``LOCK_EXCHANGE`` example.
+
+.. note::
+
+   We probably need an even more basic example than this with only one namelist and
+   minimal changes to the usrdef modules
+
+Firstly, construct the directory structure, starting in the :file:`cfgs` directory:
+
+.. code-block:: console
+
+   $ ./makenemo -n 'MY_NEW_TEST' -t 'LOCK_EXCHANGE' -m 'my_arch'
+
+where the ``-t`` option has been used to locate the new configuration in
+the :file:`tests` subdirectory
+(it is recommended practice to keep full configurations and idealised cases clearly distinguishable).
+This command will create (amongst others) the following files and directories::
+
+   ./tests/MY_NEW_TEST:
+   BLD  EXP00  MY_SRC WORK  cpp_MY_NEW_TEST.fcm
+
+   ./tests/MY_NEW_TEST/EXP00:
+   context_nemo.xml  domain_def_nemo.xml  field_def_nemo-oce.xml  file_def_nemo-oce.xml  iodef.xml
+   namelist_cfg      namelist_ref
+
+   ./tests/MY_NEW_TEST/MY_SRC:
+   usrdef_hgr.F90  usrdef_nam.F90  usrdef_zgr.F90  usrdef_istate.F90  usrdef_sbc.F90  zdfini.F90
+
+The key to setting up an idealised configuration lies in
+adapting a small set of short Fortran 90 modules which
+should be dropped into the :file:`MY_SRC` directory.
+Here the ``LOCK_EXCHANGE`` example is using 5 such routines but the full set that is available in
+the :file:`src/OCE/USR` directory is::
+
+   ./src/OCE/USR:
+   usrdef_closea.F90  usrdef_fmask.F90  usrdef_hgr.F90  usrdef_istate.F90
+   usrdef_nam.F90     usrdef_sbc.F90    usrdef_zgr.F90
+
+Before discussing these in more detail it is worth noting the various namelist controls that
+engage the different user-defined aspects.
+These controls are set using two new logical switches or are implied by the settings of existing ones.
+For example, the mandatory requirement for an idealised configuration is to provide routines which
+define the horizontal and vertical domains.
+Templates for these are provided in the :file:`usrdef_hgr.F90` and :file:`usrdef_zgr.F90` modules.
+The application of these modules is activated whenever:
+
+.. code-block:: fortran
+
+   ln_read_cfg = .false.
+
+in any configuration's :file:`namelist_cfg` file.
+This setting also activates the reading of an optional ``&nam_usrdef`` namelist which can be used to
+supply configuration specific settings.
+These need to be declared and read in the :file:`usrdef_nam.F90` module.
+
+Another explicit control is available in the ``&namsbc`` namelist which
+activates the use of analytical forcing.
+With
+
+.. code-block:: fortran
+
+   ln_usr = .true.
+
+Other usrdef modules are activated by less explicit means.
+For example, code in :file:`usrdef_istate.F90` is used to
+define initial temperature and salinity fields if
+
+.. code-block:: fortran
+
+   ln_tsd_init   = .false.
+
+in the ``&namtsd`` namelist.
+The remaining modules, namely :file:`usrdef_closea.F90` :file:`usrdef_fmask.F90` are specific to
+ORCA configurations and set local variations of some specific fields for
+the various resolutions of the global models.
+They do not need to be considered here in the context of idealised cases but
+it is worth noting that all configuration specific code has now been isolated in the usrdef modules.
+In the case of these last two modules, they are activated only if an ORCA configuration is detected.
+Currently,
+this requires a specific integer variable named ``ORCA`` to be set in a :file:`domain_cfg.nc` file.
+
+.. note::
+
+   This would be less confusing if the ``cn_cfg`` string is read directly as
+   a character attribue from the :file:`domain_cfg.nc`.
+
+So, in most cases, the set up of idealised model configurations can be completed by
+copying the template routines from :file:`./src/OCE/USR` into
+your new :file:`./cfgs/MY_NEW_TEST/MY_SRC` directory and
+editing the appropriate modules as needed.
+The default set are those used for the GYRE reference configuration.
+The contents of :file:`MY_SRC` directories from other idealised configurations may provide
+more convenient templates if they share common characteristics with your target application.
+
+Whatever the starting point,
+it should not require too many changes or additional lines of code to produce routines in
+:file:`./src/OCE/USR` that define analytically the domain,
+the initial state and the surface boundary conditions for your new configuration.
+
+To summarize, the base set of modules is:
+
+:usrdef_hgr.F90:    Define horizontal grid
+:usrdef_zgr.F90:    Define vertical grid
+:usrdef_sbc.F90:    Provides at each time-step the surface boundary condition,
+   i.e. the momentum, heat and freshwater fluxes
+:usrdef_istate.F90: Defines initialization of the dynamics and tracers
+:usrdef_nam.F90:    Configuration-specific namelist processing to
+   set any associated run-time parameters
+
+with two specialised ORCA modules
+(not related to idealised configurations but used to isolate configuration specific code that
+is used in ORCA2 reference configurations and established global configurations using
+the ORCA tripolar grid):
+
+:usrdef_fmask.F90:  only used in ORCA configurations for
+   alteration of f-point land/ocean mask in some straits
+:usrdef_closea.F90: only used in ORCA configurations for
+   specific treatments associated with closed seas
+
+From version 4.0, the NEMO release includes a :file:`tests` subdirectory containing available and
+up to date :doc:`test cases <tests>` build by the community.
+These will not be fully supported as are NEMO reference configurations,
+but should provide a source of raw material.
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OFF/nemogcm.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OFF/nemogcm.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/OFF/nemogcm.F90	(revision 12165)
@@ -114,7 +114,9 @@
 #else
                                 CALL dta_dyn    ( istp )         ! Interpolation of the dynamical fields
+#endif
+                                CALL trc_stp    ( istp )         ! time-stepping
+#if ! defined key_sed_off
          IF( .NOT.ln_linssh )   CALL dta_dyn_swp( istp )         ! swap of sea  surface height and vertical scale factors
 #endif
-                                CALL trc_stp    ( istp )         ! time-stepping
                                 CALL stp_ctl    ( istp, indic )  ! Time loop: control and print
          istp = istp + 1
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/README.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/README.rst	(revision 12165)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/README.rst	(revision 12165)
@@ -0,0 +1,15 @@
+*******
+Sources
+*******
+
+.. todo::
+
+
+
+:file:`ICE`: |NEMO-ICE|
+
+:file:`NST`: AGRIF for embedded zooms
+
+:file:`OCE`: |NEMO-OCE|
+
+:file:`TOP`: |NEMO-MBG|
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/TOP/README.rst
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/TOP/README.rst	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/TOP/README.rst	(revision 12165)
@@ -3,46 +3,51 @@
 ***************
 
+.. todo::
+
+
+
 .. contents::
-	:local:
-
-TOP (Tracers in the Ocean Paradigm) is the NEMO hardwired interface toward biogeochemical models and
-provide the physical constraints/boundaries for oceanic tracers.
-It consists of a modular framework to handle multiple ocean tracers, including also a variety of built-in modules.
+   :local:
+
+TOP (Tracers in the Ocean Paradigm) is the NEMO hardwired interface toward
+biogeochemical models and provide the physical constraints/boundaries for oceanic tracers.
+It consists of a modular framework to handle multiple ocean tracers,
+including also a variety of built-in modules.
 
 This component of the NEMO framework allows one to exploit available modules (see below) and
 further develop a range of applications, spanning from the implementation of a dye passive tracer to
 evaluate dispersion processes (by means of MY_TRC), track water masses age (AGE module),
-assess the ocean interior penetration of persistent chemical compounds (e.g., gases like CFC or even PCBs),
-up to the full set of equations involving marine biogeochemical cycles.
+assess the ocean interior penetration of persistent chemical compounds
+(e.g., gases like CFC or even PCBs), up to the full set of equations involving
+marine biogeochemical cycles.
 
 Structure
 =========
 
-TOP interface has the following location in the source code ``./src/MBG/`` and
+TOP interface has the following location in the source code :file:`./src/TOP` and
 the following modules are available:
 
-``TRP``
-	Interface to NEMO physical core for computing tracers transport
-
-``CFC``
-	Inert carbon tracers (CFC11,CFC12,SF6)
-
-``C14``
-	Radiocarbon passive tracer
-
-``AGE``
-	Water age tracking
-
-``MY_TRC``
-	Template for creation of new modules and external BGC models coupling
-
-``PISCES``
-	Built in BGC model.
-	See [https://www.geosci-model-dev.net/8/2465/2015/gmd-8-2465-2015-discussion.html Aumont et al. (2015)] for
-	a throughout description. |
-
-The usage of TOP is activated i) by including in the configuration definition  the component ``MBG`` and
-ii) by adding the macro ``key_top`` in the configuration CPP file
-(see for more details [http://forge.ipsl.jussieu.fr/nemo/wiki/Users "Learn more about the model"]).
+:file:`TRP`
+   Interface to NEMO physical core for computing tracers transport
+
+:file:`CFC`
+   Inert carbon tracers (CFC11,CFC12,SF6)
+
+:file:`C14`
+   Radiocarbon passive tracer
+
+:file:`AGE`
+   Water age tracking
+
+:file:`MY_TRC`
+   Template for creation of new modules and external BGC models coupling
+
+:file:`PISCES`
+   Built in BGC model. See :cite:`gmd-8-2465-2015` for a throughout description.
+
+The usage of TOP is activated
+*i)* by including in the configuration definition the component ``TOP`` and
+*ii)* by adding the macro ``key_top`` in the configuration CPP file
+(see for more details :forge:`"Learn more about the model" <wiki/Users>`).
 
 As an example, the user can refer to already available configurations in the code,
@@ -51,5 +56,6 @@
 (see also Section 4) .
 
-Note that, since version 4.0, TOP interface core functionalities are activated by means of logical keys and
+Note that, since version 4.0,
+TOP interface core functionalities are activated by means of logical keys and
 all submodules preprocessing macros from previous versions were removed.
 
@@ -57,118 +63,131 @@
 
 ``key_iomput``
-	use XIOS I/O
+   use XIOS I/O
 
 ``key_agrif``
-	enable AGRIF coupling
+   enable AGRIF coupling
 
 ``key_trdtrc`` & ``key_trdmxl_trc``
-	trend computation for tracers
+   trend computation for tracers
 
 Synthetic Workflow
 ==================
 
-A synthetic description of the TOP interface workflow is given below to summarize the steps involved in
-the computation of biogeochemical and physical trends and their time integration and outputs,
+A synthetic description of the TOP interface workflow is given below to
+summarize the steps involved in the computation of biogeochemical and physical trends and
+their time integration and outputs,
 by reporting also the principal Fortran subroutine herein involved.
 
-**Model initialization (OPA_SRC/nemogcm.F90)**
-
-call to trc_init (trcini.F90)
-
-  ↳ call trc_nam (trcnam.F90) to initialize TOP tracers and run setting
-
-  ↳ call trc_ini_sms, to initialize each submodule
-
-  ↳ call trc_ini_trp, to initialize transport for tracers
-
-  ↳ call trc_ice_ini, to initialize tracers in seaice
-
-  ↳ call trc_ini_state, read passive tracers from a restart or input data
-
-  ↳ call trc_sub_ini, setup substepping if {{{nn_dttrc /= 1}}}
-
-**Time marching procedure (OPA_SRC/stp.F90)**
-
-call to trc_stp.F90 (trcstp.F90)
-
-  ↳ call trc_sub_stp, averaging physical variables for sub-stepping
-
-  ↳ call trc_wri, call XIOS for output of data
-
-  ↳ call trc_sms, compute BGC trends for each submodule
-
-    ↳ call trc_sms_my_trc, includes also surface and coastal BCs trends
-
-  ↳ call trc_trp (TRP/trctrp.F90), compute physical trends
-
-    ↳ call trc_sbc, get trend due to surface concentration/dilution
-
-    ↳ call trc_adv, compute tracers advection
-
-    ↳ call to trc_ldf, compute tracers lateral diffusion
-
-    ↳ call to trc_zdf, vertical mixing and after tracer fields
-
-    ↳ call to trc_nxt, tracer fields at next time step. Lateral Boundary Conditions are solved in here.
-
-    ↳ call to trc_rad, Correct artificial negative concentrations
-
-  ↳ call trc_rst_wri, output tracers restart files
+Model initialization (:file:`./src/OCE/nemogcm.F90`)
+----------------------------------------------------
+
+Call to ``trc_init`` subroutine (:file:`./src/TOP/trcini.F90`) to initialize TOP.
+
+.. literalinclude:: ../../../src/TOP/trcini.F90
+   :language:        fortran
+   :lines:           41-86
+   :emphasize-lines: 21,30-32,38-40
+   :caption:         ``trc_init`` subroutine
+
+Time marching procedure (:file:`./src/OCE/step.F90`)
+----------------------------------------------------
+
+Call to ``trc_stp`` subroutine (:file:`./src/TOP/trcstp.F90`) to compute/update passive tracers.
+
+.. literalinclude:: ../../../src/TOP/trcstp.F90
+   :language:        fortran
+   :lines:           46-125
+   :emphasize-lines: 42,55-57
+   :caption:         ``trc_stp`` subroutine
+
+BGC trends computation for each submodule (:file:`./src/TOP/trcsms.F90`)
+------------------------------------------------------------------------
+
+.. literalinclude:: ../../../src/TOP/trcsms.F90
+   :language:        fortran
+   :lines:           21
+   :caption:         :file:`trcsms` snippet
+
+Physical trends computation (:file:`./src/TOP/TRP/trctrp.F90`)
+--------------------------------------------------------------
+
+.. literalinclude:: ../../../src/TOP/TRP/trctrp.F90
+   :language:        fortran
+   :lines:           46-95
+   :emphasize-lines: 17,21,29,33-35
+   :caption:         ``trc_trp`` subroutine
 
 Namelists walkthrough
 =====================
 
-namelist_top
-------------
-
-Here below are listed the features/options of the TOP interface accessible through the namelist_top_ref and
-modifiable by means of namelist_top_cfg (as for NEMO physical ones).
-
-Note that ## is used to refer to a number in an array field.
+:file:`namelist_top`
+--------------------
+
+Here below are listed the features/options of the TOP interface accessible through
+the :file:`namelist_top_ref` and modifiable by means of :file:`namelist_top_cfg`
+(as for NEMO physical ones).
+
+Note that ``##`` is used to refer to a number in an array field.
 
 .. literalinclude:: ../../namelists/namtrc_run
+   :language: fortran
 
 .. literalinclude:: ../../namelists/namtrc
+   :language: fortran
 
 .. literalinclude:: ../../namelists/namtrc_dta
+   :language: fortran
 
 .. literalinclude:: ../../namelists/namtrc_adv
+   :language: fortran
 
 .. literalinclude:: ../../namelists/namtrc_ldf
+   :language: fortran
 
 .. literalinclude:: ../../namelists/namtrc_rad
+   :language: fortran
 
 .. literalinclude:: ../../namelists/namtrc_snk
+   :language: fortran
 
 .. literalinclude:: ../../namelists/namtrc_dmp
+   :language: fortran
 
 .. literalinclude:: ../../namelists/namtrc_ice
+   :language: fortran
 
 .. literalinclude:: ../../namelists/namtrc_trd
+   :language: fortran
 
 .. literalinclude:: ../../namelists/namtrc_bc
+   :language: fortran
 
 .. literalinclude:: ../../namelists/namtrc_bdy
+   :language: fortran
 
 .. literalinclude:: ../../namelists/namage
-
-Two main types of data structure are used within TOP interface to initialize tracer properties (1) and
+   :language: fortran
+
+Two main types of data structure are used within TOP interface
+to initialize tracer properties (1) and
 to provide related initial and boundary conditions (2).
 
-**1. TOP tracers initialization**: sn_tracer (namtrc)
+1. TOP tracers initialization: ``sn_tracer`` (``&namtrc``)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 Beside providing name and metadata for tracers,
-here are also defined the use of initial ({{{sn_tracer%llinit}}}) and
-boundary ({{{sn_tracer%llsbc, sn_tracer%llcbc, sn_tracer%llobc}}}) conditions.
-
-In the following, an example of the full structure definition is given for two idealized tracers both with
-initial conditions given, while the first has only surface boundary forcing and
+here are also defined the use of initial (``sn_tracer%llinit``) and
+boundary (``sn_tracer%llsbc, sn_tracer%llcbc, sn_tracer%llobc``) conditions.
+
+In the following, an example of the full structure definition is given for
+two idealized tracers both with initial conditions given,
+while the first has only surface boundary forcing and
 the second both surface and coastal forcings:
 
 .. code-block:: fortran
 
-	!             !    name   !           title of the field            !   units    ! initial data ! sbc   !   cbc  !   obc  !
-	sn_tracer(1)  = 'TRC1'    , 'Tracer 1 Concentration                ',   ' - '    ,  .true.      , .true., .false., .true.
-	sn_tracer(2)  = 'TRC2 '   , 'Tracer 2 Concentration                ',   ' - '    ,  .true.      , .true., .true. , .false.
+   !             !    name   !           title of the field            !   units    ! initial data ! sbc   !   cbc  !   obc  !
+   sn_tracer(1)  = 'TRC1'    , 'Tracer 1 Concentration                ',   ' - '    ,  .true.      , .true., .false., .true.
+   sn_tracer(2)  = 'TRC2 '   , 'Tracer 2 Concentration                ',   ' - '    ,  .true.      , .true., .true. , .false.
 
 As tracers in BGC models are increasingly growing,
@@ -177,94 +196,101 @@
 .. code-block:: fortran
 
-	!             !    name   !           title of the field            !   units    ! initial data !
-	sn_tracer(1)  = 'TRC1'    , 'Tracer 1 Concentration                ',   ' - '    ,   .true.
-	sn_tracer(2)  = 'TRC2 '   , 'Tracer 2 Concentration                ',   ' - '    ,   .true.
-	! sbc
-	sn_tracer(1)%llsbc = .true.
-	sn_tracer(2)%llsbc = .true.
-	! cbc
-	sn_tracer(2)%llcbc = .true.
+   !             !    name   !           title of the field            !   units    ! initial data !
+   sn_tracer(1)  = 'TRC1'    , 'Tracer 1 Concentration                ',   ' - '    ,   .true.
+   sn_tracer(2)  = 'TRC2 '   , 'Tracer 2 Concentration                ',   ' - '    ,   .true.
+   ! sbc
+   sn_tracer(1)%llsbc = .true.
+   sn_tracer(2)%llsbc = .true.
+   ! cbc
+   sn_tracer(2)%llcbc = .true.
 
 The data structure is internally initialized by code with dummy names and
-all initialization/forcing logical fields set to .false. .
-
-**2. Structures to read input initial and boundary conditions**: namtrc_dta (sn_trcdta), namtrc_bc (sn_trcsbc/sn_trccbc/sn_trcobc)
+all initialization/forcing logical fields set to ``.false.`` .
+
+2. Structures to read input initial and boundary conditions: ``&namtrc_dta`` (``sn_trcdta``), ``&namtrc_bc`` (``sn_trcsbc`` / ``sn_trccbc`` / ``sn_trcobc``)
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 The overall data structure (Fortran type) is based on the general one defined for NEMO core in the SBC component
-(see details in User Manual SBC Chapter on Input Data specification).
-
-Input fields are prescribed within namtrc_dta (with sn_trcdta structure),
-while Boundary Conditions are applied to the model by means of namtrc_bc,
-with dedicated structure fields for surface (sn_trcsbc), riverine (sn_trccbc), and
-lateral open (sn_trcobc) boundaries.
+(see details in ``SBC`` Chapter of :doc:`Reference Manual <cite>` on Input Data specification).
+
+Input fields are prescribed within ``&namtrc_dta`` (with ``sn_trcdta`` structure),
+while Boundary Conditions are applied to the model by means of ``&namtrc_bc``,
+with dedicated structure fields for surface (``sn_trcsbc``), riverine (``sn_trccbc``), and
+lateral open (``sn_trcobc``) boundaries.
 
 The following example illustrates the data structure in the case of initial condition for
-a single tracer contained in the file named tracer_1_data.nc (.nc is implicitly assumed in namelist filename),
-with a doubled initial value, and located in the usr/work/model/inputdata/ folder:
+a single tracer contained in the file named :file:`tracer_1_data.nc`
+(``.nc`` is implicitly assumed in namelist filename),
+with a doubled initial value, and located in the :file:`usr/work/model/inputdata` folder:
 
 .. code-block:: fortran
 
-	!               !  file name             ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
-	!               !                        !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
-	  sn_trcdta(1)  = 'tracer_1_data'        ,        -12        ,  'TRC1'   ,    .false.   , .true. , 'yearly'  , ''       , ''       , ''
-	  rf_trfac(1) = 2.0
-	  cn_dir = “usr/work/model/inputdata/”
-
-Note that, the Lateral Open Boundaries conditions are applied on the segments defined for the physical core of NEMO
-(see BDY description in the User Manual).
-
-namelist_trc
-------------
-
-Here below the description of namelist_trc_ref used to handle Carbon tracers modules, namely CFC and C14.
-
-|||| &'''namcfc'''     !   CFC ||
-
-|||| &'''namc14_typ'''     !  C14 - type of C14 tracer, default values of C14/C and pco2 ||
-
-|||| &'''namc14_sbc'''     !  C14 - surface BC ||
-
-|||| &'''namc14_fcg'''     !  files & dates ||
+   !               !  file name             ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+   !               !                        !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  ! filename      !
+     sn_trcdta(1)  = 'tracer_1_data'        ,        -12        ,  'TRC1'   ,    .false.   , .true. , 'yearly'  , ''       , ''       , ''
+     rf_trfac(1) = 2.0
+     cn_dir = 'usr/work/model/inputdata/'
+
+Note that, the Lateral Open Boundaries conditions are applied on
+the segments defined for the physical core of NEMO
+(see ``BDY`` description in the :doc:`Reference Manual <cite>`).
+
+:file:`namelist_trc`
+--------------------
+
+Here below the description of :file:`namelist_trc_ref` used to handle Carbon tracers modules,
+namely CFC and C14.
+
+.. literalinclude:: ../../../cfgs/SHARED/namelist_trc_ref
+   :language: fortran
+   :lines: 7,17,26,34
+   :caption: :file:`namelist_trc_ref` snippet
 
 ``MY_TRC`` interface for coupling external BGC models
 =====================================================
 
-The generalized interface is pivoted on MY_TRC module that contains template files to build the coupling between
+The generalized interface is pivoted on MY_TRC module that contains template files to
+build the coupling between
 NEMO and any external BGC model.
 
-The call to MY_TRC is activated by setting ``ln_my_trc = .true.`` (in namtrc)
+The call to MY_TRC is activated by setting ``ln_my_trc = .true.`` (in ``&namtrc``)
 
 The following 6 fortran files are available in MY_TRC with the specific purposes here described.
 
-``par_my_trc.F90``
-	This module allows to define additional arrays and public variables to be used within the MY_TRC interface
-
-``trcini_my_trc.F90``
-	Here are initialized user defined namelists and the call to the external BGC model initialization procedures to
-	populate general tracer array (trn and trb). Here are also likely to be defined suport arrays related to
-	system metrics that could be needed by the BGC model.
-
-``trcnam_my_trc.F90``
-	This routine is called at the beginning of trcini_my_trc and should contain the initialization of
-	additional namelists for the BGC model or user-defined code.
-
-``trcsms_my_trc.F90``
-	The routine performs the call to Boundary Conditions and its main purpose is to
-	contain the Source-Minus-Sinks terms due to the biogeochemical processes of the external model.
-	Be aware that lateral boundary conditions are applied in trcnxt routine.
-	IMPORTANT: the routines to compute the light penetration along the water column and
-	the tracer vertical sinking should be defined/called in here, as generalized modules are still missing in
-	the code.
-
-``trcice_my_trc.F90``
-	Here it is possible to prescribe the tracers concentrations in the seaice that will be used as
-	boundary conditions when ice melting occurs (nn_ice_tr =1 in namtrc_ice).
-	See e.g. the correspondent PISCES subroutine.
-
-``trcwri_my_trc.F90``
-	This routine performs the output of the model tracers (only those defined in namtrc) using IOM module
-	(see Manual Chapter “Output and Diagnostics”).
-	It is possible to place here the output of additional variables produced by the model,
-	if not done elsewhere in the code, using the call to iom_put.
+:file:`par_my_trc.F90`
+   This module allows to define additional arrays and public variables to
+   be used within the MY_TRC interface
+
+:file:`trcini_my_trc.F90`
+   Here are initialized user defined namelists and
+   the call to the external BGC model initialization procedures to populate general tracer array
+   (``trn`` and ``trb``).
+   Here are also likely to be defined support arrays related to system metrics that
+   could be needed by the BGC model.
+
+:file:`trcnam_my_trc.F90`
+   This routine is called at the beginning of ``trcini_my_trc`` and
+   should contain the initialization of additional namelists for the BGC model or user-defined code.
+
+:file:`trcsms_my_trc.F90`
+   The routine performs the call to Boundary Conditions and its main purpose is to
+   contain the Source-Minus-Sinks terms due to the biogeochemical processes of the external model.
+   Be aware that lateral boundary conditions are applied in trcnxt routine.
+
+   .. warning::
+      The routines to compute the light penetration along the water column and
+      the tracer vertical sinking should be defined/called in here,
+      as generalized modules are still missing in the code.
+
+:file:`trcice_my_trc.F90`
+   Here it is possible to prescribe the tracers concentrations in the sea-ice that
+   will be used as boundary conditions when ice melting occurs (``nn_ice_tr = 1`` in ``&namtrc_ice``).
+   See e.g. the correspondent PISCES subroutine.
+
+:file:`trcwri_my_trc.F90`
+   This routine performs the output of the model tracers (only those defined in ``&namtrc``) using
+   IOM module (see chapter “Output and Diagnostics” in the :doc:`Reference Manual <cite>`).
+   It is possible to place here the output of additional variables produced by the model,
+   if not done elsewhere in the code, using the call to ``iom_put``.
 
 Coupling an external BGC model using NEMO framework
@@ -273,18 +299,20 @@
 The coupling with an external BGC model through the NEMO compilation framework can be achieved in
 different ways according to the degree of coding complexity of the Biogeochemical model, like e.g.,
-the whole code is made only by one file or it has multiple modules and interfaces spread across several subfolders.
-
-Beside the 6 core files of MY_TRC module, let’s assume an external BGC model named *MYBGC* and constituted by
-a rather essential coding structure, likely few Fortran files.
+the whole code is made only by one file or
+it has multiple modules and interfaces spread across several subfolders.
+
+Beside the 6 core files of MY_TRC module, let’s assume an external BGC model named *MYBGC* and
+constituted by a rather essential coding structure, likely few Fortran files.
 The new coupled configuration name is *NEMO_MYBGC*.
 
-The best solution is to have all files (the modified ``MY_TRC`` routines and the BGC model ones) placed in
-a unique folder with root ``MYBGCPATH`` and to use the makenemo external readdressing of ``MY_SRC`` folder.
-
-The coupled configuration listed in ``work_cfgs.txt`` will look like
+The best solution is to have all files (the modified ``MY_TRC`` routines and the BGC model ones)
+placed in a unique folder with root ``MYBGCPATH`` and
+to use the makenemo external readdressing of ``MY_SRC`` folder.
+
+The coupled configuration listed in :file:`work_cfgs.txt` will look like
 
 ::
 
-	NEMO_MYBGC OPA_SRC TOP_SRC
+   NEMO_MYBGC OCE TOP
 
 and the related ``cpp_MYBGC.fcm`` content will be
@@ -292,46 +320,46 @@
 .. code-block:: perl
 
-	bld::tool::fppkeys  key_iomput key_mpp_mpi key_top
-
-the compilation with ``makenemo`` will be executed through the following syntax
+   bld::tool::fppkeys key_iomput key_mpp_mpi key_top
+
+the compilation with :file:`makenemo` will be executed through the following syntax
 
 .. code-block:: console
 
-	$ makenemo -n 'NEMO_MYBGC' -m '<arch_my_machine>' -j 8 -e '<MYBGCPATH>'
-
-The makenemo feature “-e” was introduced to readdress at compilation time the standard MY_SRC folder
-(usually found in NEMO configurations) with a user defined external one.
-
-The compilation of more articulated BGC model code & infrastructure, like in the case of BFM
-([http://www.bfm-community.eu/publications/bfmnemomanual_r1.0_201508.pdf BFM-NEMO coupling manual]),
-requires some additional features.
+   $ makenemo -n 'NEMO_MYBGC' -m '<arch_my_machine>' -j 8 -e '<MYBGCPATH>'
+
+The makenemo feature ``-e`` was introduced to
+readdress at compilation time the standard MY_SRC folder (usually found in NEMO configurations) with
+a user defined external one.
+
+The compilation of more articulated BGC model code & infrastructure,
+like in the case of BFM (|BFM man|_), requires some additional features.
 
 As before, let’s assume a coupled configuration name *NEMO_MYBGC*,
-but in this case MYBGC model root becomes ``<MYBGCPATH>`` that contains 4 different subfolders for
-biogeochemistry, named ``initialization``, ``pelagic``, and ``benthic``, and
-a separate one named ``nemo_coupling`` including the modified ``MY_SRC`` routines.
+but in this case MYBGC model root becomes :file:`MYBGC` path that
+contains 4 different subfolders for biogeochemistry,
+named :file:`initialization`, :file:`pelagic`, and :file:`benthic`,
+and a separate one named :file:`nemo_coupling` including the modified `MY_SRC` routines.
 The latter folder containing the modified NEMO coupling interface will be still linked using
-the makenemo “-e” option.
+the makenemo ``-e`` option.
 
 In order to include the BGC model subfolders in the compilation of NEMO code,
-it will be necessary to extend the configuration ``cpp_NEMO_MYBGC.fcm`` file to include the specific paths of
-``MYBGC`` folders, as in the following example
+it will be necessary to extend the configuration :file:`cpp_NEMO_MYBGC.fcm` file to include the specific paths of :file:`MYBGC` folders, as in the following example
 
 .. code-block:: perl
 
-	bld::tool::fppkeys  key_iomput key_mpp_mpi key_top
-	
-	src::MYBGC::initialization         <MYBGCPATH>/initialization
-	src::MYBGC::pelagic                <MYBGCPATH>/pelagic
-	src::MYBGC::benthic                <MYBGCPATH>/benthic
-	
-	bld::pp::MYBGC      1
-	bld::tool::fppflags::MYBGC   %FPPFLAGS
-	bld::tool::fppkeys           %bld::tool::fppkeys MYBGC_MACROS
+   bld::tool::fppkeys  key_iomput key_mpp_mpi key_top
+
+   src::MYBGC::initialization         <MYBGCPATH>/initialization
+   src::MYBGC::pelagic                <MYBGCPATH>/pelagic
+   src::MYBGC::benthic                <MYBGCPATH>/benthic
+
+   bld::pp::MYBGC      1
+   bld::tool::fppflags::MYBGC   %FPPFLAGS
+   bld::tool::fppkeys           %bld::tool::fppkeys MYBGC_MACROS
 
 where *MYBGC_MACROS* is the space delimited list of macros used in *MYBGC* model for
 selecting/excluding specific parts of the code.
-The BGC model code will be preprocessed in the configuration ``BLD`` folder as for NEMO,
-but with an independent path, like ``NEMO_MYBGC/BLD/MYBGC/<subforlders>``.
+The BGC model code will be preprocessed in the configuration :file:`BLD` folder as for NEMO,
+but with an independent path, like :file:`NEMO_MYBGC/BLD/MYBGC/<subforlders>`.
 
 The compilation will be performed similarly to in the previous case with the following
@@ -339,20 +367,26 @@
 .. code-block:: console
 
-	$ makenemo -n 'NEMO_MYBGC' -m '<arch_my_machine>' -j 8 -e '<MYBGCPATH>/nemo_coupling'
-
-Note that, the additional lines specific for the BGC model source and build paths can be written into
-a separate file, e.g. named ``MYBGC.fcm``, and then simply included in the ``cpp_NEMO_MYBGC.fcm`` as follow
-
-.. code-block:: perl
-
-	bld::tool::fppkeys  key_zdftke key_dynspg_ts key_iomput key_mpp_mpi key_top
-	inc <MYBGCPATH>/MYBGC.fcm
-
-This will enable a more portable compilation structure for all MYBGC related configurations.
-
-**Important**: the coupling interface contained in nemo_coupling cannot be added using the FCM syntax,
-as the same files already exists in NEMO and they are overridden only with the readdressing of MY_SRC contents to
-avoid compilation conflicts due to duplicate routines.
-
-All modifications illustrated above, can be easily implemented using shell or python scripting to
-edit the NEMO configuration CPP.fcm file and to create the BGC model specific FCM compilation file with code paths.
+   $ makenemo -n 'NEMO_MYBGC' -m '<arch_my_machine>' -j 8 -e '<MYBGCPATH>/nemo_coupling'
+
+.. note::
+   The additional lines specific for the BGC model source and build paths can be written into
+   a separate file, e.g. named :file:`MYBGC.fcm`,
+   and then simply included in the :file:`cpp_NEMO_MYBGC.fcm` as follow
+
+   .. code-block:: perl
+
+      bld::tool::fppkeys  key_zdftke key_dynspg_ts key_iomput key_mpp_mpi key_top
+      inc <MYBGCPATH>/MYBGC.fcm
+
+   This will enable a more portable compilation structure for all MYBGC related configurations.
+
+.. warning::
+   The coupling interface contained in :file:`nemo_coupling` cannot be added using the FCM syntax,
+   as the same files already exists in NEMO and they are overridden only with
+   the readdressing of MY_SRC contents to avoid compilation conflicts due to duplicate routines.
+
+All modifications illustrated above, can be easily implemented using shell or python scripting
+to edit the NEMO configuration :file:`CPP.fcm` file and
+to create the BGC model specific FCM compilation file with code paths.
+
+.. |BFM man| replace:: BFM-NEMO coupling manual
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/TOP/TRP/trcnxt.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/TOP/TRP/trcnxt.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/TOP/TRP/trcnxt.F90	(revision 12165)
@@ -139,5 +139,14 @@
       ENDIF
       !                                ! Leap-Frog + Asselin filter time stepping
-      IF( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) THEN    ! Euler time-stepping (only swap)
+      IF( (neuler == 0 .AND. kt == nittrc000) ) THEN
+         ! set up for leapfrog on second timestep
+         DO jn = 1, jptra
+            DO jk = 1, jpkm1
+               trb(:,:,jk,jn) = trn(:,:,jk,jn)  
+               trn(:,:,jk,jn) = tra(:,:,jk,jn)
+            END DO
+         END DO
+      ELSE IF( ln_top_euler ) THEN
+         ! always doing euler timestepping
          DO jn = 1, jptra
             DO jk = 1, jpkm1
@@ -146,4 +155,6 @@
             END DO
          END DO
+      ENDIF
+      IF( (neuler == 0 .AND. kt == nittrc000) .OR. ln_top_euler ) THEN    ! Euler time-stepping (only swap)
          IF (l_trdtrc .AND. .NOT. ln_linssh ) THEN   ! Zero Asselin filter contribution must be explicitly written out since for vvl
             !                                        ! Asselin filter is output by tra_nxt_vvl that is not called on this time step
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/src/TOP/trcbdy.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/src/TOP/trcbdy.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/src/TOP/trcbdy.F90	(revision 12165)
@@ -95,5 +95,5 @@
          END DO
          IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
-            CALL lbc_lnk( 'bdytra', tsa, 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
+            CALL lbc_lnk( 'trcbdy', tra, 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
          END IF
          !
Index: NEMO/branches/2019/dev_ASINTER-01-05_merged/tests/CANAL/MY_SRC/usrdef_nam.F90
===================================================================
--- NEMO/branches/2019/dev_ASINTER-01-05_merged/tests/CANAL/MY_SRC/usrdef_nam.F90	(revision 12164)
+++ NEMO/branches/2019/dev_ASINTER-01-05_merged/tests/CANAL/MY_SRC/usrdef_nam.F90	(revision 12165)
@@ -86,5 +86,5 @@
       REWIND( numnam_cfg )          ! Namelist namusr_def (exist in namelist_cfg only)
       READ  ( numnam_cfg, namusr_def, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist', cdtxt )
+902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namusr_def in configuration namelist' )
       !
       IF(lwm)   WRITE( numond, namusr_def )
