diff --git a/latex/NEMO/main/bibliography.bib b/latex/NEMO/main/bibliography.bib
index 62a015d322e818720ef78d3656e1d611c84d1202..a51cf1c21aad34c9534a9410a14304941002e17f 100644
--- a/latex/NEMO/main/bibliography.bib
+++ b/latex/NEMO/main/bibliography.bib
@@ -14,6 +14,18 @@
   doi           = "10.1016/j.ocemod.2003.09.003"
 }
 
+@article{         andreas_JMS15,
+  title         = "An improved bulk air-sea surface flux algorithm,
+                  including spray-mediated transfer.",
+  pages         = "642--654",
+  journal       = "Q.J.R. Meteorol. Soc.",
+  volume        = "141",
+  author        = "E.L. Andreas and L. Mahrt and D. Vickers",
+  year          = "2015",
+  doi           = "10.1002/qj.2424"
+}
+
+
 @article{         arakawa.hsu_MWR90,
   title         = "Energy Conserving and Potential-Enstrophy Dissipating
                   Schemes for the Shallow Water Equations",
@@ -248,7 +260,7 @@
 
 @article{         bignami_JGR95,
   title         = " Longwave radiation budget in the mediterranean sea",
-  pages         = "2501�2514",
+  pages         = "2501--2514",
   journal     = "Journal of Geophysical Research",
   volume        = "100",
   author        = "F. Bignami and S. Marullo and R. Santoleri and M. E. Schiano",
@@ -3247,6 +3259,20 @@ forecasting",
   doi           = "10.1029/91jc02350"
 }
 
+@article{         tsamados_JPO14,
+  title         = "Impact of Variable Atmospheric and Oceanic Form Drag on Simulations of Arctic Sea Ice",
+  pages         = "1329--1353",
+  journal       = "Journal of Physical Oceanography",
+  volume        = "44",
+  number        = "5",
+  author        = "M. Tsamados",
+  year          = "2014",
+  month         = "may",
+  publisher     = "American Meteorological Society (AMS)",
+  issn          = "0022-3670",
+  doi           = "10.1175/JPO-D-13-0215.1"
+}
+
 @article{         tu.tsuang_GRL05,
   title         = "Cool-skin simulation by a one-column ocean model",
   pages         = "L22602",
diff --git a/latex/NEMO/subfiles/chap_SBC.tex b/latex/NEMO/subfiles/chap_SBC.tex
index e73c948769830bc0932a0129d3746a6d4dbd4050..4fc1a9c682d25da93c2bed4cebf5f52a074522ef 100644
--- a/latex/NEMO/subfiles/chap_SBC.tex
+++ b/latex/NEMO/subfiles/chap_SBC.tex
@@ -2,7 +2,7 @@
 
 \begin{document}
 
-\chapter{Surface Boundary Condition (SBC, SAS, ISF, ICB, TDE)}
+\chapter{Surface Boundary Condition (SBC, SAS, TDE)}
 \label{chap:SBC}
 
 \chaptertoc
@@ -13,9 +13,10 @@
   \begin{tabularx}{\textwidth}{l||X|X}
     Release & Author(s) & Modifications \\
     \hline
-    {\em  next} & {\em A. Moulin, E. Clementi} & {\em Update of \autoref{sec:SBC_wave}}\\[2mm]
-    {\em  next} & {\em Simon M{\" u}ller} & {\em Update of \autoref{sec:SBC_TDE}; revision of \autoref{subsec:SBC_fwb}}\\[2mm]
-    {\em  next} & {\em Pierre Mathiot} & {\em update of the ice shelf section (2019 developments)}\\[2mm]  
+    {\em   5.0} & {\em A. Moulin, G. Samson and P. Mathiot} & {\em Update of the SBC chapter and moving ICB and ISF in there own chapter}\\[2mm]
+    {\em   4.2} & {\em A. Moulin, E. Clementi} & {\em Update of \autoref{sec:SBC_wave}}\\[2mm]
+    {\em   4.2} & {\em Simon M{\" u}ller} & {\em Update of \autoref{sec:SBC_TDE}; revision of \autoref{subsec:SBC_fwb}}\\[2mm]
+    {\em   4.2} & {\em Pierre Mathiot} & {\em update of the ice shelf section (2019 developments)}\\[2mm]  
     {\em   4.0} & {\em ...} & {\em ...} \\
     {\em   3.6} & {\em ...} & {\em ...} \\
     {\em   3.4} & {\em ...} & {\em ...} \\
@@ -45,11 +46,10 @@ Five different ways are available to provide these fields to the ocean. They are
 namelist \nam{sbc}{sbc} variables:
 
 \begin{itemize}
-\item a bulk formulation (\np[=.true.]{ln_blk}{ln\_blk}), featuring a selection of five bulk parameterization algorithms,
-\item an atmospheric boundary layer model (\np[=.true.]{ln_abl}{ln\_abl}) associated with the bulk formulation,
-\item a flux formulation (\np[=.true.]{ln_flx}{ln\_flx}),
-\item a coupled or mixed forced/coupled formulation (exchanges with a atmospheric model via the OASIS coupler),
-(\np{ln_cpl}{ln\_cpl} or \np[=.true.]{ln_mixcpl}{ln\_mixcpl}),
+\item a bulk formulation (\autoref{sec:SBC_blk}), featuring a selection of six bulk parameterization algorithms,
+\item an atmospheric boundary layer model (\autoref{sec:SBC_abl}) associated with the bulk formulation,
+\item a flux formulation (\autoref{sec:SBC_flx}),
+\item a coupled or mixed forced/coupled formulation (exchanges with an atmospheric model via the OASIS coupler)(\autoref{sec:SBC_cpl}),
 \item a user defined formulation (\np[=.true.]{ln_usr}{ln\_usr}).
 \end{itemize}
 
@@ -68,19 +68,19 @@ These options control:
 
 \begin{itemize}
 \item the rotation of vector components supplied relative to an east-north coordinate system onto
-  the local grid directions in the model,
-\item the use of a land/sea mask for input fields (\np[=.true.]{nn_lsm}{nn\_lsm}),
-\item the addition of a surface restoring term to observed SST and/or SSS (\np[=.true.]{ln_ssr}{ln\_ssr}),
+  the local grid directions in the model (\autoref{subsec:SBC_rotation}),
+\item the use of a land/sea mask for input fields (\autoref{subsec:SBC_iof}),
+\item the addition of a surface restoring term to observed SST and/or SSS (\autoref{subsec:SBC_ssr}),
 \item the modification of fluxes below ice-covered areas (using climatological ice-cover or a sea-ice model)
-  (\np[=0..3]{nn_ice}{nn\_ice}),
-\item the addition of river runoffs as surface freshwater fluxes or lateral inflow (\np[=.true.]{ln_rnf}{ln\_rnf}),
+  (\autoref{subsec:SBC_ice-cover}),
+\item the addition of river runoffs as surface freshwater fluxes or lateral inflow (\autoref{sec:SBC_rnf}),
 \item the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift
-  (\np[=0..2]{nn_fwb}{nn\_fwb}),
+  (\autoref{subsec:SBC_fwb}),
 \item the transformation of the solar radiation (if provided as daily mean) into an analytical diurnal cycle
-  (\np[=.true.]{ln_dm2dc}{ln\_dm2dc}),
-\item the activation of wave effects from an external wave model  (\np[=.true.]{ln_wave}{ln\_wave}),
-\item the light penetration in the ocean (\np[=.true.]{ln_traqsr}{ln\_traqsr} with \nam{tra_qsr}{tra\_qsr}),
-\item the atmospheric surface pressure gradient effect on ocean and ice dynamics (\np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} with namelist \nam{sbc_apr}{sbc\_apr}),
+  (\autoref{subsec:SBC_dcy}),
+\item the activation of wave effects from an external wave model  (\autoref{sec:SBC_wave}),
+\item the light penetration in the ocean (\autoref{subsec:TRA_qsr}),
+\item the atmospheric surface pressure gradient effect on ocean and ice dynamics (\autoref{sec:SBC_apr}),
 \item the effect of sea-ice pressure on the ocean (\np[=.true.]{ln_ice_embd}{ln\_ice\_embd}).
 \end{itemize}
 
@@ -89,8 +89,6 @@ Then we present the four ways of providing the surface boundary conditions,
 followed by the description of the atmospheric pressure and the river runoff.
 Next, the scheme for interpolation on the fly is described.
 Finally, the different options that further modify the fluxes applied to the ocean are discussed.
-One of these is modification by icebergs (see \autoref{sec:SBC_ICB_icebergs}),
-which act as drifting sources of fresh water.
 
 %% =================================================================================================
 \section{Surface boundary condition for the ocean}
@@ -221,7 +219,7 @@ where
     \centering
     \begin{tabular}{|l|c|c|c|}
       \hline
-                                  &  daily or weekLL     &  monthly           &  yearly        \\
+                                  &  daily or weekLLL     &  monthly           &  yearly        \\
       \hline
       \texttt{clim=.false.} &  fn\_yYYYYmMMdDD.nc  &  fn\_yYYYYmMM.nc   &  fn\_yYYYY.nc  \\
       \hline
@@ -260,7 +258,7 @@ where
   If you want to change this behaviour, it is possible to prepend the variable name with a '-' or a '+' sign.
   In the first case, the records will be dated at the beginning of the forcing period,
   while in the second case, the records will be dated at the end of the forcing period.
-\item [Climatological forcing]: a logical to specify if a input file contains climatological forcing which can be cycle in time,
+\item [Climatological forcing]: a logical to specify if an input file contains climatological forcing which can be cycle in time,
   or an interannual forcing which will requires additional files if
   the period covered by the simulation exceeds the one of the file.
   See the above file naming strategy which impacts the expected name of the file to be opened.
@@ -431,7 +429,7 @@ or is a copy of one from the first few columns on the opposite side of the grid
 
 % to be completed
 A set of utilities to create a weights file for a rectilinear input grid is available
-(see the directory NEMOGCM/TOOLS/WEIGHTS).
+(see the directory /tools/WEIGHTS).
 
 %% =================================================================================================
 \subsection{Standalone surface boundary condition scheme (SAS)}
@@ -457,17 +455,23 @@ For example:
 \end{itemize}
 
 The Standalone Surface scheme provides this capacity.
-Its options are defined through the \nam{sbc_sas}{sbc\_sas} namelist variables.
-A new copy of the model has to be compiled with a configuration based on ORCA2\_SAS\_LIM.
-However, no namelist parameters need be changed from the settings of the previous run (except perhaps nn\_date0).
-In this configuration, a few routines in the standard model are overriden by new versions.
+Its options are defined through the \nam{sbc_sas}{sbc\_sas} namelist variables. Here are the available options : 
+
+\begin{itemize}
+\item \np[=.true.]{l_sasread}{l\_sasread}: ocean fields are coming from netcdf files 
+\item \np[=.false.]{l_sasread}{l\_sasread}: ocean fields are coming from OASIS
+\end{itemize}
+
+To use SAS, the model has to be compiled with the ORCA2\_SAS\_ICE configuration. In standalone mode (\np[=.true.]{l_sasread}{l\_sasread}), the namelist parameters that are set in cfgs/ORCA2\_SAS\_ICE/EXPREF/namelist\_cfg can be used directly. In 'coupled mode' (\np[=.false.]{l_sasread}{l\_sasread} ) with OASIS, the user need to compile two executable (one for OCE and one for SAS). In the namelist of the first one \np{nn_components}{nn\_components} need to be set to 1 and in the sas namelist set to 2. However, the coupling of SAS and NEMO via OASIS does not seem to work properly yet. This is under investigation.
+
+In both cases, a few routines in the standard model are overwritten by SAS specific versions.
 Routines replaced are:
 
 \begin{itemize}
 \item \mdl{nemogcm}: This routine initialises the rest of the model and repeatedly calls the stp time stepping routine (\mdl{step}).
   Since the ocean state is not calculated all associated initialisations have been removed.
 \item \mdl{step}: The main time stepping routine now only needs to call the sbc routine (and a few utility functions).
-\item \mdl{sbcmod}: This has been cut down and now only calculates surface forcing and the ice model required.
+\item \mdl{sbcmod}: This has been cut down and now only calculates surface forcing, the ice model required.
   New surface modules that can function when only the surface level of the ocean state is defined can also be added
   (\eg\ icebergs).
 \item \mdl{daymod}: No ocean restarts are read or written (though the ice model restarts are retained),
@@ -478,14 +482,8 @@ Routines replaced are:
 \item \mdl{diawri}: All 3D data have been removed from the output.
   The surface temperature, salinity and velocity components (which have been read in) are written along with
   relevant forcing and ice data.
-\end{itemize}
-
-One new routine has been added:
-
-\begin{itemize}
-\item \mdl{sbcsas}: This module initialises the input files needed for reading temperature, salinity and
-  velocity arrays at the surface.
-  These filenames are supplied in namelist namsbc\_sas.
+\item \mdl{sbcssm}: When \np[=.true.]{l_sasread}{l\_sasread}, this module initialises the input files needed for reading temperature, salinity,
+  velocity arrays at the surface, surface vertical scale factor and fraction of energy absorbed by the first level (optional, \np[=.true.]{ln_read_frq}). These parameters need to supplied in the namelist namsbc\_sas parameters \np{sn_usp}{sn\_usp}, \np{sn_vsp}{sn\_vsp}, \np{sn_tem}{sn\_tem}, \np{sn_sal}{sn\_sal}, \np{sn_ssh}{sn\_ssh}, \np{sn_e3t}{sn\_e3t} and \np{sn_frq}{sn\_frq} respectively. If For velocities, either 2D velocities (\np[=.false]{ln_3d_uve}{ln\_3d\_uve}) or 3D ones (\np[=.true.]{ln_3d_uve}{ln\_3d\_uve}) can be read from the input files. 
   Unfortunately, because of limitations with the \mdl{iom} module,
   the full 3D fields from the mean files have to be read in and interpolated in time,
   before using just the top level.
@@ -624,13 +622,14 @@ parametrization} algorithm. When relevant, these algorithms also perform the
 height adjustment of humidity and temperature to the wind reference measurement
 height (from \np{rn_zqt}{rn\_zqt} to \np{rn_zu}{rn\_zu}).
 
-For the open ocean, five bulk parametrization algorithms are available in NEMO:
+For the open ocean, six bulk parametrization algorithms are available in NEMO:
 
 \begin{itemize}
 \item NCAR, formerly known as CORE, \citep{large.yeager_trpt04,large.yeager_CD09}
 \item COARE 3.0 \citep{fairall.bradley.ea_JC03}
 \item COARE 3.6 \citep{edson.jampana.ea_JPO13}
-\item ECMWF (IFS documentation, cy45)
+\item ECMWF (IFS documentation, cy45)
+\item ANDREAS \citep{andreas_JMS15}
 \item MFS \citep{castellari_JMS98}
 \end{itemize}
 
@@ -731,14 +730,14 @@ ECMWF:
 
 \begin{forlines}
   ...
-  ln_COARE3p6 = .true.
+  ln_COARE_3p6 = .true.
   ...     
   ln_skin_cs = .true. ! use the cool-skin parameterization
   ln_skin_wl = .true. ! use the warm-layer parameterization
   ...
 \end{forlines}
 
-Note: when \np[=T]{ln_COARE3p0}{ln\_COARE3p0} is selected, the selection
+Note: when \np[=T]{ln_COARE_3p0}{ln\_COARE\_3p0} is selected, the selection
 of \np{ln_skin_cs}{ln\_skin\_cs} and \np{ln_skin_wl}{ln\_skin\_wl} implicitly
 triggers the use of the COARE cool-skin and warm-layer parametrizations,
 respectively (found in \textit{sbcblk\_skin\_coare.F90}).
@@ -752,20 +751,20 @@ respectively (found in \textit{sbcblk\_skin\_coare.F90}).
 % functions of the wind speed and the near-surface stability of the atmospheric
 % surface layer (hereafter ASL), and hence, depend on $U_B$, $T_s$, $T_z$, $q_s$
 % and $q_z$.
-
-\subsubsection{MFS}
-
-The MFS bulk formulae have been developed by \cite{castellari_JMS98}. 
-They have been designed to handle ECMWF operational data for the Mediterranean \citep{oddo_OS09} and Black Sea \citep{ciliberti_JMSE22} Monitoring Forecasting Centre. 
-The bulk transfer coefficients $C_H$ and $C_E$ are computing using an empiric formula by \cite{kondo_BLM75} whereas the drag coefficient $C_D$  is computed according to \cite{hellerman_JPO83}. 
-In this bulk parametrization the net solar radiation depends on the cloud cover and 
-is computed by means of an astronomical formula \citep{reed_JPO77}. 
-Albedo monthly values are from \cite{payne_JAS72} as means of the values at 40\textdegree N and 30\textdegree N 
-for the Atlantic Ocean (hence the same latitudinal band of the Mediterranean Sea). 
-The net long-wave radiation flux are computed as a function of 
-2m air temperature, sea-surface temperature, cloud cover and 2m dew point humidity  \citep{Bignami_JGR95}. 
-
-This parametrization required as input fields the total cloud cover in \%.
+
+\subsubsection{MFS}
+
+The MFS bulk formulae have been developed by \cite{castellari_JMS98}.
+They have been designed to handle ECMWF operational data for the Mediterranean \citep{oddo_OS09} and Black Sea \citep{ciliberti_JMSE22} Monitoring Forecasting Centre. 
+The bulk transfer coefficients $C_H$ and $C_E$ are computing using an empiric formula by \cite{kondo_BLM75} whereas the drag coefficient $C_D$  is computed according to \cite{hellerman_JPO83}. 
+In this bulk parametrization the net solar radiation depends on the cloud cover and 
+is computed by means of an astronomical formula \citep{reed_JPO77}. 
+Albedo monthly values are from \cite{payne_JAS72} as means of the values at 40\textdegree N and 30\textdegree N 
+for the Atlantic Ocean (hence the same latitudinal band of the Mediterranean Sea). 
+The net long-wave radiation flux are computed as a function of 
+2m air temperature, sea-surface temperature, cloud cover and 2m dew point humidity  \citep{Bignami_JGR95}. 
+
+This parametrization required as input fields the total cloud cover in \%.
 
 The following namelist parameters must be set:
 
@@ -775,9 +774,9 @@ The following namelist parameters must be set:
   ...
   rn_zqt     = 2.     ! Air temperature & humidity reference height (m)
   rn_zu      = 10.   ! Wind vector reference height (m)
-  ...
-  ln_skin_cs = .false. ! use the cool-skin parameterization
-  ln_skin_wl = .false. ! use the warm-layer parameterization
+  ...
+  ln_skin_cs = .false. ! use the cool-skin parameterization
+  ln_skin_wl = .false. ! use the warm-layer parameterization
   ...
   ln_humi_dpt = .true. ! humidity "sn_humi" is dew point Temperature  [kg/kg]
 \end{forlines}
@@ -786,33 +785,65 @@ The following namelist parameters must be set:
 \subsection{Ice-Atmosphere Bulk formulae}
 \label{subsec:SBC_blk_ice}
 
-\texttt{\#out\_of\_place:}
- For sea-ice, three possibilities can be selected:
-a constant transfer coefficient (1.4e-3; default
-value), \citet{lupkes.gryanik.ea_JGRA12} (\np{ln_Cd_L12}{ln\_Cd\_L12}),
-and \citet{lupkes.gryanik_JGR15} (\np{ln_Cd_L15}{ln\_Cd\_L15}) parameterizations
-\texttt{\#out\_of\_place.}
-
 Surface turbulent fluxes between sea-ice and the atmosphere can be computed in three different ways:
 
 \begin{itemize}
-\item Constant value (\forcode{Cd_ice=1.4e-3}):
-  default constant value used for momentum and heat neutral transfer coefficients
-\item \citet{lupkes.gryanik.ea_JGRA12} (\np[=.true.]{ln_Cd_L12}{ln\_Cd\_L12}):
-  This scheme adds a dependency on edges at leads, melt ponds and flows
-  of the constant neutral air-ice drag. After some approximations,
-  this can be resumed to a dependency on ice concentration (A).
-  This drag coefficient has a parabolic shape (as a function of ice concentration)
-  starting at 1.5e-3 for A=0, reaching 1.97e-3 for A=0.5 and going down 1.4e-3 for A=1.
-  It is theoretically applicable to all ice conditions (not only MIZ).
-\item \citet{lupkes.gryanik_JGR15} (\np[=.true.]{ln_Cd_L15}{ln\_Cd\_L15}):
-  Alternative turbulent transfer coefficients formulation between sea-ice
-  and atmosphere with distinct momentum and heat coefficients depending
-  on sea-ice concentration and atmospheric stability (no melt-ponds effect for now).
-  The parameterization is adapted from ECHAM6 atmospheric model.
-  Compared to Lupkes2012 scheme, it considers specific skin and form drags
-  to compute neutral transfer coefficients for both heat and momentum fluxes.
-  Atmospheric stability effect on transfer coefficient is also taken into account.
+  \item Constant value (\np[=.true.]{ln_Cx_ice_cst}{ln\_Cx\_ice\_cst}):
+    Constant values are used for momentum (\np{rn_Cd_ia}{rn\_Cd\_ia}), sublimation (\np{rn_Cd_ia}{rn\_Ce\_ia})
+    and sensible heat (\np{rn_Ch_ia}{rn\_Ch\_ia}) transfer coefficients.
+    Default value for all coefficients is set to 1.4e-3.
+  \item \citet{lupkes.gryanik.ea_JGRA12} (\np[=.true.]{ln_Cx_ice_LU12}{ln\_Cx\_ice\_LU12}):
+    This scheme adds a dependency on edges at leads, melt ponds and flows
+    of the constant neutral air-ice drag. After some approximations,
+    this can be resumed to a dependency on ice concentration (A).
+    This drag coefficient has a parabolic shape (as a function of ice concentration)
+    starting at 1.5e-3 for A=0, reaching 1.97e-3 for A=0.5 and going down 1.4e-3 for A=1.
+    It is theoretically applicable to all ice conditions (not only MIZ).
+  \item \citet{lupkes.gryanik_JGR15} (\np[=.true.]{ln_Cx_ice_LG15}{ln\_Cx\_ice\_LG15}):
+    Alternative turbulent transfer coefficients formulation between sea-ice
+    and atmosphere with distinct momentum and heat coefficients depending
+    on sea-ice concentration and atmospheric stability (no melt-ponds effect for now).
+    The parameterization is adapted from ECHAM6 atmospheric model.
+    Compared to Lupkes2012 scheme, it considers specific skin and form drags
+    to compute neutral transfer coefficients for both heat and momentum fluxes.
+    Atmospheric stability effect on transfer coefficient is also taken into account.
+  \item \citet{tsamados_JPO14} (\np[=.true.]{ln_Cx_ice_frm}{ln\_Cx\_ice\_frm}):
+    Sea ice contains pressure ridges as well as floe and melt pond edges
+    that act as discrete obstructions to the flow of air or water past the ice, and are a source of form drag.
+    Here, the neutral drag coefficients are estimated from sea ice properties such as ice concentration,
+    vertical extent and area of the ridges, freeboard and floe draft, and size of floes and melt ponds.
+    The new parameterization allows the drag coefficients to be coupled to the sea ice state and
+    therefore to evolve spatially and temporally.
+    For default settings, only the transfer coefficient for the turbulent momentum fluxes are affected. The namelist allows 3 different options:
+    \begin{description}
+      \item [{\np[=1]{nn_frm}{nn\_frm}}:] affects momentum and heat transfer coefficient (ocean-ice and atm-ice)
+      \item [{\np[=2]{nn_frm}{nn\_frm}}:] affects only momentum transfer coefficient (ocean-ice and atm-ice) (default)
+      \item [{\np[=3]{nn_frm}{nn\_frm}}:] affects  momentum and heat transfer coefficient (atm-ice), and only momentum transfer coefficient (ocean-ice)
+    \end{description}
+    The total drag coefficients are derived as follows:
+    \[ zdrag_{ia} = zdrag_{ia\_skin} + zdrag_{ia\_floe} + zdrag_{ia\_rdg} + zdrag_{ia\_pond} \]
+    \[ zdrag_{io} = zdrag_{io\_skin} + zdrag_{io\_floe} + zdrag_{io\_keel} \]
+    \begin{description}
+      \item \np{zdrag_ia}{zdrag\_ia} : total drag coefficient for momentum exchange between ice and atmosphere
+      \item \np{zdrag_ia_skin}{zdrag\_ia\_skin} : skin drag coefficient (top of the ice)
+      \item \np{zdrag_ia_floe}{zdrag\_ia\_floe} : floe edge drag coefficient (top of the ice)
+      \item \np{zdrag_ia_rdg}{zdrag\_ia\_rdg}: sail drag coefficient
+      \item \np{zdrag_ia_pond}{zdrag\_ia\_rdg} : pond edge drag coefficient
+      \item \np{zdrag_io}{zdrag\_io} : total drag for momentum exhange between ice and ocean
+      \item \np{zdrag_io_skin}{zdrag\_io\_skin} : skin drag coefficient (under the ice)
+      \item \np{zdrag_io_floe}{zdrag\_io\_floe} : floe edge drag coefficient (under the ice)
+      \item \np{zdrag_io_rdg}{zdrag\_io\_rdg} : keel drag coefficient
+    \end{description}
+    The area and volume of ridged ice (required input parameters) are derived from mean ice thickness and concentration based on a polynomial fit (\citet{tournay_JPO14}).
+    In the namelist, the skin drag coefficients and factors can be modified for the strength of each contribution:
+    \begin{description}
+      \item [{\np[=0.002]{rn_Cs_io}{rn\_Cs\_io}}:] ice-ocean skin drag [0.0005,0.005]
+      \item [{\np[=0.0005]{rn_Cs_ia}{rn\_Cs\_ia}}:] ice-air skin drag [0.0001,0.001]
+      \item [{\np[=0.2]{rn_Cr_ia}{rn\_Cr\_ia}}:] factor for ridge/keel drag coefficient [0,1]
+      \item [{\np[=0.2]{rn_Cr_io}{rn\_Cr\_io}}:] factor for ridge/keel drag coefficient [0,1]
+      \item [{\np[=0.2]{rn_Cf_ia}{rn\_Cf\_ia}}:] factor for floe edge atm [0,1]
+      \item [{\np[=0.2]{rn_Cf_io}{rn\_Cf\_io}}:] factor floe edge ocean [0,1]
+    \end{description}
 \end{itemize}
 
 %% =================================================================================================
@@ -840,7 +871,8 @@ The required 9 input fields are:
     \hline
     j-component of the 10m air velocity  & wndj           & $m.s^{-1}$         & T     \\
     \hline
-    10m air temperature                  & tair           & $K$               & T     \\
+    10m absolute air temperature         & tair           & $K$               & T     \\
+    Potential air temperature            & ~              & $K$               & T     \\
     \hline
     Specific humidity                    & humi           & $-$               & T     \\
     Relative humidity                    & ~              & $\%$              & T     \\
@@ -851,6 +883,7 @@ The required 9 input fields are:
     Downwelling shortwave radiation      & qsr            & $W.m^{-2}$         & T     \\
     \hline
     Total precipitation (liquid + solid) & precip         & $Kg.m^{-2}.s^{-1}$ & T     \\
+    ~ 					 & ~              & $m$                 & T     \\
     \hline
     Solid precipitation                  & snow           & $Kg.m^{-2}.s^{-1}$ & T     \\
     \hline
@@ -883,17 +916,22 @@ Its range must be between zero and one, and it is recommended to set it to 0 at
 As for the flux parametrization, information about the input data required by the model is provided in
 the namsbc\_blk namelist (see \autoref{subsec:SBC_fldread}).
 
-\subsubsection{Air humidity}
+\subsubsection{Air humidity, temperature, precipitation parameters and units:}
+
+Air humidity can be provided as one of three parameters: \\
+specific humidity [$kg/kg$], using \np[=.true.]{ln_humi_sph}{ln\_humi\_sph}; 
+relative humidity [$\%$], using \np[=.true.]{ln_humi_rlh}{ln\_humi\_rlh}; 
+or dew-point temperature [$K$], using \np[=.true.]{ln_humi_dpt}{ln\_humi\_dpt}.
+
+Air temperature can be provided as potential temperature (\np[=.true.]{ln_tair_pot}{ln\_tair\_pot}) or as absolute temperature (\np[=.false.]{ln_tair_pot}{ln\_tair\_pot}).
 
-Air humidity can be provided as three different parameters: specific humidity
-[kg/kg], relative humidity [\%], or dew-point temperature [K] (LINK to namelist
-parameters)...
+Total precipitation can be provided in meters (\np[=.true.]{ln_prec_met}{ln\_prec\_met}) or, by default, in $kg \cdot m^{-2} \cdot s^{-1}$.
 
 %% =================================================================================================
 \section[Atmospheric Boundary Layer (ABL) model (\textit{sbcabl.F90})]{Atmospheric Boundary Layer (ABL) model (\protect\mdl{sbcabl})}
 \label{sec:SBC_abl}
 
-An atmospheric boundary layer (ABL) model is available as an alternative choice to the prescribed near-surface atmospheric forcings.
+An atmospheric boundary layer (ABL) model is available as an alternative choice to the prescribed near-surface atmospheric forcings using \np[=.true.]{ln_abl}{ln\_abl}.
 It computes the wind, air potential temperature and specific humidity evolutions in the lower atmosphere following a single-column approach
 on the same horizontal grid as the ocean component. It represents the adjustement of the air column between the large-scale
 atmospheric forcing and the surface boundary conditions over both ocean and sea-ice through vertical turbulent mixing.
@@ -938,6 +976,7 @@ Each of the three steps needed to generate the atmospheric forcings corresponds
 ABL1D model is activated by adding ABL sources directory to the sources list file (\*\_cfgs.txt) and
 by setting \np[=.true.]{ln_abl}{ln\_abl} (and \np[=.false.]{ln_blk}{ln\_bkl}) in \nam{sbc}{sbc}. \\
 It is fully compatible with Nemo Standalone Surface module and can be consequently forced by sea surface temperature and currents external data.\\
+The namelist \nam{sbc_abl}{sbc\_abl} is used to setup the ABL1D options.
 
 Atmospheric forcing files needed by ABL1D must be specified directly using the \np{sn_wndi}{sn\_wndi}, \np{sn_wndj}{sn\_wndj},
 \np{sn_tair}{sn\_tair} and \np{sn_humi}{sn\_humi} parameters from the \nam{sbc_blk}{sbc\_blk}.\\
@@ -1107,7 +1146,7 @@ included in the
   ocean tide model}: Mf, Mm, Ssa, Mtm, Msf, Msqm, Sa, K1, O1, P1, Q1, J1, S1,
 M2, S2, N2, K2, nu2, mu2, 2N2, L2, T2, eps2, lam2, R2, M3, MKS2, MN4, MS4, M4,
 N4, S4, M6, and M8; see file \textit{tide.h90} and \mdl{tide\_mod} for further
-information and references\footnote{As a legacy option \np{ln_tide_var}{ln\_tide\_var} can be
+information and references\footnote{As a legacy option \np{nn_tide_var}{nn\_tide\_var} can be
   set to \forcode{0}, in which case the 19 tidal constituents (M2, N2, 2N2, S2,
   K2, K1, O1, Q1, P1, M4, Mf, Mm, Msqm, Mtm, S1, MU2, NU2, L2, and T2; see file
   \textit{tide.h90}) and associated parameters that have been available in NEMO version
@@ -1164,12 +1203,13 @@ with a spatially uniform coefficient $\beta$, which can be configured
 via \np{rn_scal_load}{rn\_scal\_load} (default value 0.094) and is
 often tuned to minimize tidal prediction errors.\par
 
-For diagnostic purposes, the forcing potential of the individual tidal
-constituents (incl. load ptential, if activated) and the total forcing
-potential (incl. load potential, if activated) can be made available
-as diagnostic output by setting
-\np[=.true.]{ln_tide_dia}{ln\_tide\_dia} (fields
-\forcode{tide_pot_<constituent>} and \forcode{tide_pot}).\par
+%Tidal diagnostics not anymore available
+%For diagnostic purposes, the forcing potential of the individual tidal
+%constituents (incl. load ptential, if activated) and the total forcing
+%potential (incl. load potential, if activated) can be made available
+%as diagnostic output by setting
+%\np[=.true.]{ln_tide_dia}{ln\_tide\_dia} (fields
+%\forcode{tide_pot_<constituent>} and \forcode{tide_pot}).\par
 
 %% =================================================================================================
 \section[River runoffs (\textit{sbcrnf.F90})]{River runoffs (\protect\mdl{sbcrnf})}
@@ -1181,23 +1221,6 @@ as diagnostic output by setting
   \label{lst:namsbc_rnf}
 \end{listing}
 
-%River runoff generally enters the ocean at a nonzero depth rather than through the surface.
-%Many models, however, have traditionally inserted river runoff to the top model cell.
-%This was the case in \NEMO\ prior to the version 3.3. The switch toward a input of runoff
-%throughout a nonzero depth has been motivated by the numerical and physical problems
-%that arise when the top grid cells are of the order of one meter. This situation is common in
-%coastal modelling and becomes more and more often open ocean and climate modelling
-%\footnote{At least a top cells thickness of 1~meter and a 3 hours forcing frequency are
-%required to properly represent the diurnal cycle \citep{bernie.woolnough.ea_JC05}. see also \autoref{fig:SBC_dcy}.}.
-
-%To do this we need to treat evaporation/precipitation fluxes and river runoff differently in the
-%\mdl{tra\_sbc} module.  We decided to separate them throughout the code, so that the variable
-%\textit{emp} represented solely evaporation minus precipitation fluxes, and a new 2d variable
-%rnf was added which represents the volume flux of river runoff (in kg/m2s to remain consistent with
-%emp).  This meant many uses of emp and emps needed to be changed, a list of all modules which use
-%emp or emps and the changes made are below:
-
-%Rachel:
 River runoff generally enters the ocean at a nonzero depth rather than through the surface.
 Many models, however, have traditionally inserted river runoff to the top model cell.
 This was the case in \NEMO\ prior to the version 3.3,
@@ -1215,35 +1238,40 @@ and for the temperature and salinity of the river to effect the surrounding ocea
 The user is able to specify, in a NetCDF input file, the temperature and salinity of the river,
 along with the depth (in metres) which the river should be added to.
 
-Namelist variables in \nam{sbc_rnf}{sbc\_rnf}, \np{ln_rnf_depth}{ln\_rnf\_depth}, \np{ln_rnf_sal}{ln\_rnf\_sal} and
-\np{ln_rnf_temp}{ln\_rnf\_temp} control whether the river attributes (depth, salinity and temperature) are read in and used.
-If these are set as false the river is added to the surface box only, assumed to be fresh (0~psu),
-and/or taken as surface temperature respectively.
+The surface runoff is activated via the namelist parameter \np{ln_rnf}{ln\_rnf} in \nam{sbc}{sbc}. The specific options that control the runoff are described in \nam{sbc_rnf}{sbc\_rnf} (\autoref{lst:namsbc_rnf}). In case of activation, the mandatory fields is the map of surface runoff that need to be specified by the user (\np{sn_rnf}{sn\_rnf}). All the other parameters described below like the injection depth, the temperature, salinity (...) are optional.\\
 
-The runoff value and attributes are read in in sbcrnf.
-For temperature -999 is taken as missing data and the river temperature is taken to
-be the surface temperatue at the river point.
-For the depth parameter a value of -1 means the river is added to the surface box only,
-and a value of -999 means the river is added through the entire water column.
-After being read in the temperature and salinity variables are multiplied by the amount of runoff
-(converted into m/s) to give the heat and salt content of the river runoff.
-After the user specified depth is read ini,
-the number of grid boxes this corresponds to is calculated and stored in the variable \np{nz_rnf}{nz\_rnf}.
-The variable \textit{h\_dep} is then calculated to be the depth (in metres) of
+By default the surface runoff is injected in the surface level. However, by activating \np{ln_rnf_depth}{ln\_rnf\_depth}, it is possible to specify a 2D map of depth over which to inject it (\np{sn_dep_rnf}{sn\_dep\_rnf}). The depth variable is expected to be positive with two specific values. A cell value of -1 means the river is added to the surface box only,
+and a cell value of -999 means the river is added through the entire water column. There is a third way to define the runoff thickness. The depth can be automatically computed (\np{ln_rnf_depth_ini}{ln\_rnf\_depth\_ini}). The user simply need to give the depth over which the maximum runoff (\np{rn_rnf_max}{rn\_rnf\_max}) is spread (\np{rn_dep_max}{rn\_dep\_max}). THen the runoff depth is computed simply by scaling the local maximum runoff value over time by the ratio $\frac{rn\_dep\_max}{rn\_rnf\_max}$. Once computed, by setting \np{nn_rnf_depth_file}{nn\_rnf\_depth\_file} the file can be outputted for sanity check.
+
+After the depth is defined,
+the number of grid boxes this corresponds to is calculated and stored in the variable \textit{nk\_rnf}.
+The variable \textit{h\_dep} is then updated to be the depth (in metres) of
 the bottom of the lowest box the river water is being added to
-(\ie\ the total depth that river water is being added to in the model).
+(\ie\ the total depth that river water is being added to in the model).\\
 
 The mass/volume addition due to the river runoff is, at each relevant depth level, added to
-the horizontal divergence (\textit{hdivn}) in the subroutine \rou{sbc\_rnf\_div} (called from \mdl{divhor}).
-This increases the diffusion term in the vicinity of the river, thereby simulating a momentum flux.
-The sea surface height is calculated using the sum of the horizontal divergence terms,
+the horizontal divergence (\textit{hdivn}) in the subroutine \rou{sbc\_rnf\_div} (called from \mdl{divhor}) simulating a momentum flux.
+The sea surface height is then calculated using the sum of the horizontal divergence terms,
 and so the river runoff indirectly forces an increase in sea surface height.
 
 The \textit{hdivn} terms are used in the tracer advection modules to force vertical velocities.
 This causes a mass of water, equal to the amount of runoff, to be moved into the box above.
 The heat and salt content of the river runoff is not included in this step,
 and so the tracer concentrations are diluted as water of ocean temperature and salinity is moved upward out of
-the box and replaced by the same volume of river water with no corresponding heat and salt addition.
+the box and replaced by the same volume of river water with no corresponding heat and salt addition. \\
+
+For the tracers, by default (\np{ln_rnf_sal}{ln\_rnf\_sal} and
+\np{ln_rnf_tem}{ln\_rnf\_tem} set to false), the runoff is assumed to be at the river point sst and with 0 salinity (fresh).
+If set to true, the user needs to provide a runoff temperature and/or salinity field (\np{sn_t_rnf}{sn\_t\_rnf} and/or (\np{sn_s_rnf}{sn\_s\_rnf}) respectively). It is worth noting that location with a temperature value of -999 is considered as missing data and the river temperature is taken to
+be the surface temperature at the river point. Furthermore, the forcing iceberg melt flux (and associated latent heat flux) can be added as runoff by activating \np{ln_rnf_icb}{ln\_rnf\_icb} instead of using the lagrangian iceberg model (ICB, \autoref{subsec:ICB}) to simulate it.
+ In this case, the user simply need to specify a map of iceberg melt rate in the file \np{sn_i_rnf}{sn\_i\_rnf}. 
+ In this case, the iceberg fresh water flux is added to the runoff fluxes and the latent heat flux directly to the non solar heat fluxes.
+
+After being read in the temperature and salinity variables are multiplied by the amount of runoff
+(converted into m/s) to give the heat and salt content of the river runoff. These fluxes are then added to the tracer trend in \mdl{trasbc} (\autoref{subsec:TRA_sbc}) .
+This is done in the same way for both linear and non-linear free surface.
+The temperature and salinity are increased through the specified depth according to
+the heat and salt content of the river.\\
 
 For the linear free surface case, at the surface box the tracer advection causes a flux of water
 (of equal volume to the runoff) through the sea surface out of the domain,
@@ -1252,449 +1280,26 @@ As such the volume of water does not change, but the water is diluted.
 
 For the non-linear free surface case, no flux is allowed through the surface.
 Instead in the surface box (as well as water moving up from the boxes below) a volume of runoff water is added with
-no corresponding heat and salt addition and so as happens in the lower boxes there is a dilution effect.
+the corresponding heat and salt (runoff temperature at surface temperature and 0 salinity by default)  and so as happens in the lower boxes there is a dilution effect.
 (The runoff addition to the top box along with the water being moved up through
-boxes below means the surface box has a large increase in volume, whilst all other boxes remain the same size)
-
-In trasbc the addition of heat and salt due to the river runoff is added.
-This is done in the same way for both linear and non-linear free surface.
-The temperature and salinity are increased through the specified depth according to
-the heat and salt content of the river.
+boxes below means the surface box has a large increase in volume, whilst all other boxes remain the same size).
 
-In the non-linear free surface case (\np[=.false.]{ln_linssh}{ln\_linssh}),
-near the end of the time step the change in sea surface height is redistrubuted through the grid boxes,
+Furthermore, near the end of the time step the change in sea surface height is redistributed through the grid boxes,
 so that the original ratios of grid box heights are restored.
 In doing this water is moved into boxes below, throughout the water column,
-so the large volume addition to the surface box is spread between all the grid boxes.
+so the large volume addition to the surface box is spread between all the grid boxes.\\
+
+In addition, there is possibility to increase vertical mixing near river mouths (\np{ln_rnf_mouth}{ln\_rnf\_mouth}). If activated, the kz is increased by \np{rn_avt_rnf}{rn\_avt\_rnf} where
+the mask file \np{sn_cnf}{sn\_cnf} is set to 0.5 over the number of level corresponding to the depth \np{rn_hrnf}{rn\_hrnf}. This depth value can be different to the one set in \np{sn_dep_rnf}{sn\_dep\_rnf}.\\
 
-It is also possible for runnoff to be specified as a negative value for modelling flow through straits,
+Finally, it is worth nothing that it is also possible for runoff to be specified as a negative value for modelling flow through straits,
 \ie\ modelling the Baltic flow in and out of the North Sea.
 When the flow is out of the domain there is no change in temperature and salinity,
 regardless of the namelist options used,
-as the ocean water leaving the domain removes heat and salt (at the same concentration) with it.
-
-%\colorbox{yellow}{Nevertheless, Pb of vertical resolution and 3D input : increase vertical mixing near river mouths to mimic a 3D river
-
-%All river runoff and emp fluxes are assumed to be fresh water (zero salinity) and at the same temperature as the sea surface.}
-
-%\colorbox{yellow}{river mouths{\ldots}}
-
-%IF( ln_rnf ) THEN                                     ! increase diffusivity at rivers mouths
-%        DO jk = 2, nkrnf   ;   avt(:,:,jk) = avt(:,:,jk) + rn_avt_rnf * rnfmsk(:,:)   ;   END DO
-%ENDIF
-
-\cmtgm{  word doc of runoffs:
-In the current \NEMO\ setup river runoff is added to emp fluxes,
-these are then applied at just the sea surface as a volume change (in the variable volume case
-this is a literal volume change, and in the linear free surface case the free surface is moved)
-and a salt flux due to the concentration/dilution effect.
-There is also an option to increase vertical mixing near river mouths;
-this gives the effect of having a 3d river.
-All river runoff and emp fluxes are assumed to be fresh water (zero salinity) and
-at the same temperature as the sea surface.
-Our aim was to code the option to specify the temperature and salinity of river runoff,
-(as well as the amount), along with the depth that the river water will affect.
-This would make it possible to model low salinity outflow, such as the Baltic,
-and would allow the ocean temperature to be affected by river runoff.
-
-The depth option makes it possible to have the river water affecting just the surface layer,
-throughout depth, or some specified point in between.
-
-To do this we need to treat evaporation/precipitation fluxes and river runoff differently in
-the \mdl{tra_sbc} module.
-We decided to separate them throughout the code,
-so that the variable emp represented solely evaporation minus precipitation fluxes,
-and a new 2d variable rnf was added which represents the volume flux of river runoff
-(in $kg/m^2s$ to remain consistent with $emp$).
-This meant many uses of emp and emps needed to be changed,
-a list of all modules which use $emp$ or $emps$ and the changes made are below:}
-
-%% =================================================================================================
-\section[Ice Shelf (ISF)]{Interaction with ice shelves (ISF)}
-\label{sec:SBC_isf}
-
-\begin{listing}
-  \nlst{namisf}
-  \caption{\forcode{&namisf}}
-  \label{lst:namisf}
-\end{listing}
-
-The ice shelf interaction module requires the addiition of the macro \key{isf} in the configuration cpp file and 
-it is activated by means of the namelist variable \np{ln_isf}{ln\_isf} in \nam{isf}{isf}. The following interactions modes are available:
-\begin{description}
-   \item $\bullet$ representation of the ice shelf/ocean melting/freezing for opened cavity (cav, \np{ln_isfcav_mlt}{ln\_isfcav\_mlt}).
-   \item $\bullet$ parametrisation of the ice shelf/ocean melting/freezing for closed cavities (par, \np{ln_isfpar_mlt}{ln\_isfpar\_mlt}).
-   \item $\bullet$ coupling with an ice sheet model (\np{ln_isfcpl}{ln\_isfcpl}).
-\end{description}
-
-  \subsection{Ocean/Ice shelf fluxes in opened cavities}
-
-     \np{ln_isfcav_mlt}{ln\_isfcav\_mlt}\forcode{ = .true.} activates the ocean/ice shelf thermodynamics interactions at the ice shelf/ocean interface. 
-     If \np{ln_isfcav_mlt}{ln\_isfcav\_mlt}\forcode{ = .false.}, thermodynamics interactions are desctivated but the ocean dynamics inside the cavity is still active.
-     The logical flag \np{ln_isfcav}{ln\_isfcav} control whether or not the ice shelf cavities are closed. \np{ln_isfcav}{ln\_isfcav} is not defined in the namelist but in the domcfg.nc input file.\\
-
-     3 options are available to represent to ice-shelf/ocean fluxes at the interface:
-     \begin{description}
-        \item[\np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = 'spe'}]:
-        The fresh water flux is specified by a forcing fields \np{sn_isfcav_fwf}{sn\_isfcav\_fwf}. Convention of the input file is: positive toward the ocean (i.e. positive for melting and negative for freezing).
-        The latent heat fluxes is derived from the fresh water flux. 
-        The heat content flux is derived from the fwf flux assuming a temperature set to the freezing point in the top boundary layer (\np{rn_htbl}{rn\_htbl})
-
-        \item[\np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = 'oasis'}]:
-        The \forcode{'oasis'} is a prototype of what could be a method to spread precipitation on Antarctic ice sheet as ice shelf melt inside the cavity when a coupled model Atmosphere/Ocean is used. 
-        It has not been tested and therefore the model will stop if you try to use it. 
-        Actions will be undertake in 2020 to build a comprehensive interface to do so for Greenland, Antarctic and ice shelf (cav), ice shelf (par), icebergs, subglacial runoff and runoff.
-
-        \item[\np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = '2eq'}]:
-        The heat flux and the fresh water flux (negative for melting) resulting from ice shelf melting/freezing are parameterized following \citet{Grosfeld1997}. 
-        This formulation is based on a balance between the vertical diffusive heat flux across the ocean top boundary layer (\autoref{eq:ISOMIP1}) 
-        and the latent heat due to melting/freezing (\autoref{eq:ISOMIP2}):
-
-        \begin{equation}
-        \label{eq:ISOMIP1}
-        \mathcal{Q}_h = \rho c_p \gamma (T_w - T_f)
-        \end{equation}
-        \begin{equation}
-        \label{eq:ISOMIP2}
-        q = \frac{-\mathcal{Q}_h}{L_f}
-        \end{equation}
-        
-        where $\mathcal{Q}_h$($W.m^{-2}$) is the heat flux,q($kg.s^{-1}m^{-2}$) the fresh-water flux, 
-        $L_f$ the specific latent heat, $T_w$ the temperature averaged over a boundary layer below the ice shelf (explained below), 
-        $T_f$ the freezing point using  the  pressure  at  the  ice  shelf  base  and  the  salinity  of the water in the boundary layer, 
-        and $\gamma$ the thermal exchange coefficient.
-
-        \item[\np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = '3eq'}]:
-        For realistic studies, the heat and freshwater fluxes are parameterized following \citep{Jenkins2001}. This formulation is based on three equations: 
-        a balance between the vertical diffusive heat flux across the boundary layer 
-        , the latent heat due to melting/freezing of ice and the vertical diffusive heat flux into the ice shelf (\autoref{eq:3eq1}); 
-        a balance between the vertical diffusive salt flux across the boundary layer and the salt source or sink represented by the melting/freezing (\autoref{eq:3eq2}); 
-        and a linear equation for the freezing temperature of sea water (\autoref{eq:3eq3}, detailed of the linearisation coefficient in \citet{AsayDavis2016}):
-
-        \begin{equation}
-        \label{eq:3eq1}
-        c_p \rho \gamma_T (T_w-T_b) = -L_f q - \rho_i c_{p,i} \kappa \frac{T_s - T_b}{h_{isf}}
-        \end{equation}
-        \begin{equation}
-        \label{eq:3eq2}
-        \rho \gamma_S (S_w - S_b) = (S_i - S_b)q
-        \end{equation}
-        \begin{equation}
-        \label{eq:3eq3}
-        T_b = \lambda_1 S_b + \lambda_2 +\lambda_3 z_{isf}
-        \end{equation}
-
-        where $T_b$ is the temperature at the interface, $S_b$ the salinity at the interface, $\gamma_T$ and $\gamma_S$ the exchange coefficients for temperature and salt, respectively, 
-        $S_i$ the salinity of the ice (assumed to be 0), $h_{isf}$ the ice shelf thickness, $z_{isf}$ the ice shelf draft, $\rho_i$ the density of the iceshelf, 
-        $c_{p,i}$ the specific heat capacity of the ice, $\kappa$ the thermal diffusivity of the ice 
-        and $T_s$ the atmospheric surface temperature (at the ice/air interface, assumed to be -20C). 
-        The Liquidus slope ($\lambda_1$), the liquidus intercept ($\lambda_2$) and the Liquidus pressure coefficient ($\lambda_3$) 
-        for TEOS80 and TEOS10 are described in \citep{AsayDavis2016} and in \citep{Jourdain2017}.
-        The linear system formed by \autoref{eq:3eq1}, \autoref{eq:3eq2} and the linearised equation for the freezing temperature of sea water (\autoref{eq:3eq3}) can be solved for $S_b$ or $T_b$. 
-        Afterward, the freshwater flux ($q$) and the heat flux ($\mathcal{Q}_h$) can be computed.
-
-     \end{description}
-
-     \begin{table}[h]
-        \centering
-        \caption{Description of the parameters hard coded into the ISF module}
-        \label{tab:isf}
-        \begin{tabular}{|l|l|l|l|}
-        \hline
-        Symbol    & Description               & Value              & Unit               \\
-        \hline
-        $C_p$     & Ocean specific heat       & 3992               & $J.kg^{-1}.K^{-1}$ \\
-        $L_f$     & Ice latent heat of fusion & $3.34 \times 10^5$ & $J.kg^{-1}$        \\
-        $C_{p,i}$ & Ice specific heat         & 2000               & $J.kg^{-1}.K^{-1}$ \\
-        $\kappa$  & Heat diffusivity          & $1.54 \times 10^{-6}$& $m^2.s^{-1}$     \\
-        $\rho_i$  & Ice density               & 920                & $kg.m^3$           \\
-        \hline
-        \end{tabular}
-     \end{table}
-
-     Temperature and salinity used to compute the fluxes in \autoref{eq:ISOMIP1}, \autoref{eq:3eq1} and \autoref{eq:3eq2} are the average temperature in the top boundary layer \citep{losch_JGR08}. 
-     Its thickness is defined by \np{rn_htbl}{rn\_htbl}.
-     The fluxes and friction velocity are computed using the mean temperature, salinity and velocity in the first \np{rn_htbl}{rn\_htbl} m.
-     Then, the fluxes are spread over the same thickness (ie over one or several cells).
-     If \np{rn_htbl}{rn\_htbl} is larger than top $e_{3}t$, there is no more direct feedback between the freezing point at the interface and the top cell temperature.
-     This can lead to super-cool temperature in the top cell under melting condition.
-     If \np{rn_htbl}{rn\_htbl} smaller than top $e_{3}t$, the top boundary layer thickness is set to the top cell thickness.\\
-
-     Each melt formula (\np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = '3eq'} or \np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = '2eq'}) depends on an exchange coeficient ($\Gamma^{T,S}$) between the ocean and the ice.
-     Below, the exchange coeficient $\Gamma^{T}$ and $\Gamma^{S}$ are respectively defined by \np{rn_gammat0}{rn\_gammat0} and \np{rn_gammas0}{rn\_gammas0}. 
-     There are 3 different ways to compute the exchange velocity:
-
-     \begin{description}
-        \item[\np{cn_gammablk}{cn\_gammablk}\forcode{='spe'}]:
-        The salt and heat exchange coefficients are constant and defined by:
-\[
-\gamma^{T} = \Gamma^{T}
-\]
-\[
-\gamma^{S} = \Gamma^{S}
-\] 
-        This is the recommended formulation for ISOMIP.
-
-	\item[\np{cn_gammablk}{cn\_gammablk}\forcode{='vel'}]:
-        The salt and heat exchange coefficients are velocity dependent and defined as
-\[
-\gamma^{T} = \Gamma^{T} \times u_{*} 
-\]
-\[
-\gamma^{S} = \Gamma^{S} \times u_{*}
-\]
-        where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_htbl}{rn\_htbl} meters).
-        See \citet{jenkins.nicholls.ea_JPO10} for all the details on this formulation. It is the recommended formulation for realistic application and ISOMIP+/MISOMIP configuration.
-
-	\item[\np{cn_gammablk}{cn\_gammablk}\forcode{'vel\_stab'}]:
-        The salt and heat exchange coefficients are velocity and stability dependent and defined as:
-\[
-\gamma^{T,S} = \frac{u_{*}}{\Gamma_{Turb} + \Gamma^{T,S}_{Mole}} 
-\]
-        where $u_{*}$ is the friction velocity in the top boundary layer (ie first \np{rn_tbl}{rn\_htbl} meters),
-        $\Gamma_{Turb}$ the contribution of the ocean stability and
-        $\Gamma^{T,S}_{Mole}$ the contribution of the molecular diffusion.
-        See \citet{holland.jenkins_JPO99} for all the details on this formulation. 
-        This formulation has not been extensively tested in NEMO (not recommended).
-     \end{description}
-
-\subsection{Ocean/Ice shelf fluxes in parametrised cavities}
-
-  \begin{description}
-
-     \item[\np{cn_isfpar_mlt}{cn\_isfpar\_mlt}\forcode{ = 'bg03'}]:
-     The ice shelf cavities are not represented.
-     The fwf and heat flux are computed using the \citet{beckmann.goosse_OM03} parameterisation of isf melting.
-     The fluxes are distributed along the ice shelf edge between the depth of the average grounding line (GL)
-     (\np{sn_isfpar_zmax}{sn\_isfpar\_zmax}) and the base of the ice shelf along the calving front
-     (\np{sn_isfpar_zmin}{sn\_isfpar\_zmin}) as in (\np{cn_isfpar_mlt}{cn\_isfpar\_mlt}\forcode{ = 'spe'}).
-     The effective melting length (\np{sn_isfpar_Leff}{sn\_isfpar\_Leff}) is read from a file.
-     This parametrisation has not been tested since a while and based on \citet{Favier2019}, 
-     this parametrisation should probably not be used.
-
-     \item[\np{cn_isfpar_mlt}{cn\_isfpar\_mlt}\forcode{ = 'spe'}]:
-     The ice shelf cavity is not represented.
-     The fwf (\np{sn_isfpar_fwf}{sn\_isfpar\_fwf}) is prescribed and distributed along the ice shelf edge between
-     the depth of the average grounding line (GL) (\np{sn_isfpar_zmax}{sn\_isfpar\_zmax}) and
-     the base of the ice shelf along the calving front (\np{sn_isfpar_zmin}{sn\_isfpar\_min}). Convention of the input file is positive toward the ocean (i.e. positive for melting and negative for freezing).
-     The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$.
-
-     \item[\np{cn_isfpar_mlt}{cn\_isfpar\_mlt}\forcode{ = 'oasis'}]:
-     The \forcode{'oasis'} is a prototype of what could be a method to spread precipitation on Antarctic ice sheet as ice shelf melt inside the cavity when a coupled model Atmosphere/Ocean is used. 
-     It has not been tested and therefore the model will stop if you try to use it. 
-     Action will be undertake in 2020 to build a comprehensive interface to do so for Greenland, Antarctic and ice shelf (cav), ice shelf (par), icebergs, subglacial runoff and runoff.
-
-  \end{description}
-
-\np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = '2eq'}, \np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = '3eq'} and \np{cn_isfpar_mlt}{cn\_isfpar\_mlt}\forcode{ = 'bg03'} compute a melt rate based on
-the water mass properties, ocean velocities and depth.
-The resulting fluxes are thus highly dependent of the model resolution (horizontal and vertical) and 
-realism of the water masses onto the shelf.\\
-
-\np{cn_isfcav_mlt}{cn\_isfcav\_mlt}\forcode{ = 'spe'} and \np{cn_isfpar_mlt}{cn\_isfpar\_mlt}\forcode{ = 'spe'} read the melt rate from a file.
-You have total control of the fwf forcing.
-This can be useful if the water masses on the shelf are not realistic or
-the resolution (horizontal/vertical) are too coarse to have realistic melting or
-for studies where you need to control your heat and fw input. 
-However, if your forcing is not consistent with the dynamics below you can reach unrealistic low water temperature.\\
-
-The ice shelf fwf is implemented as a volume flux as for the runoff.
-The fwf addition due to the ice shelf melting is, at each relevant depth level, added to
-the horizontal divergence (\textit{hdivn}) in the subroutine \rou{isf\_hdiv}, called from \mdl{divhor}.
-See the runoff section \autoref{sec:SBC_rnf} for all the details about the divergence correction.\\
-
-Description and result of sensitivity tests to \np{ln_isfcav_mlt}{ln\_isfcav\_mlt} and \np{ln_isfpar_mlt}{ln\_isfpar\_mlt} are presented in \citet{mathiot.jenkins.ea_GMD17}. 
-The different options are illustrated in \autoref{fig:ISF}.
-
-\begin{figure}[!t]
-  \centering
-  \includegraphics[width=0.66\textwidth]{SBC_isf_v4.2}
-  \caption[Ice shelf location and fresh water flux definition]{
-    Illustration of the location where the fwf is injected and
-    whether or not the fwf is interactive or not.}
-  \label{fig:ISF}
-\end{figure}
-
-\subsection{Available outputs}
-The following outputs are availables via XIOS:
-\begin{description}
-   \item[for parametrised cavities]:
-      \begin{xmllines}
- <field id="isftfrz_par"     long_name="freezing point temperature in the parametrization boundary layer" unit="degC"     />
- <field id="fwfisf_par"      long_name="Ice shelf melt rate"                           unit="kg/m2/s"  />
- <field id="qoceisf_par"     long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     />
- <field id="qlatisf_par"     long_name="Ice shelf latent heat flux"                    unit="W/m2"     />
- <field id="qhcisf_par"      long_name="Ice shelf heat content flux of injected water" unit="W/m2"     />
- <field id="fwfisf3d_par"    long_name="Ice shelf melt rate"                           unit="kg/m2/s"  grid_ref="grid_T_3D" />
- <field id="qoceisf3d_par"   long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
- <field id="qlatisf3d_par"   long_name="Ice shelf latent heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
- <field id="qhcisf3d_par"    long_name="Ice shelf heat content flux of injected water" unit="W/m2"     grid_ref="grid_T_3D" />
- <field id="ttbl_par"        long_name="temperature in the parametrisation boundary layer" unit="degC" />
- <field id="isfthermald_par" long_name="thermal driving of ice shelf melting"          unit="degC"     />
-      \end{xmllines}
-   \item[for open cavities]:
-      \begin{xmllines}
- <field id="isftfrz_cav"     long_name="freezing point temperature at ocean/isf interface"                unit="degC"     />
- <field id="fwfisf_cav"      long_name="Ice shelf melt rate"                           unit="kg/m2/s"  />
- <field id="qoceisf_cav"     long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     />
- <field id="qlatisf_cav"     long_name="Ice shelf latent heat flux"                    unit="W/m2"     />
- <field id="qhcisf_cav"      long_name="Ice shelf heat content flux of injected water" unit="W/m2"     />
- <field id="fwfisf3d_cav"    long_name="Ice shelf melt rate"                           unit="kg/m2/s"  grid_ref="grid_T_3D" />
- <field id="qoceisf3d_cav"   long_name="Ice shelf ocean  heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
- <field id="qlatisf3d_cav"   long_name="Ice shelf latent heat flux"                    unit="W/m2"     grid_ref="grid_T_3D" />
- <field id="qhcisf3d_cav"    long_name="Ice shelf heat content flux of injected water" unit="W/m2"     grid_ref="grid_T_3D" />
- <field id="ttbl_cav"        long_name="temperature in Losch tbl"                      unit="degC"     />
- <field id="isfthermald_cav" long_name="thermal driving of ice shelf melting"          unit="degC"     />
- <field id="isfgammat"       long_name="Ice shelf heat-transfert velocity"             unit="m/s"      />
- <field id="isfgammas"       long_name="Ice shelf salt-transfert velocity"             unit="m/s"      />
- <field id="stbl"            long_name="salinity in the Losh tbl"                      unit="1e-3"     />
- <field id="utbl"            long_name="zonal current in the Losh tbl at T point"      unit="m/s"      />
- <field id="vtbl"            long_name="merid current in the Losh tbl at T point"      unit="m/s"      />
- <field id="isfustar"        long_name="ustar at T point used in ice shelf melting"    unit="m/s"      />
- <field id="qconisf"         long_name="Conductive heat flux through the ice shelf"    unit="W/m2"     />
-      \end{xmllines}
-\end{description}
-
-%% =================================================================================================
-\subsection{Ice sheet coupling}
-\label{subsec:ISF_iscpl}
-
-Ice sheet/ocean coupling is done through file exchange at the restart step.
-At each restart step, the procedure is this one:
-
-\begin{description}
-\item[Step 1]: the ice sheet model send a new bathymetry and ice shelf draft netcdf file.
-\item[Step 2]: a new domcfg.nc file is built using the DOMAINcfg tools.
-\item[Step 3]: NEMO run for a specific period and output the average melt rate over the period.
-\item[Step 4]: the ice sheet model run using the melt rate outputed in step 3.
-\item[Step 5]: go back to 1.
-\end{description}
-
-If \np{ln_iscpl}{ln\_iscpl}\forcode{ = .true.}, the isf draft is assume to be different at each restart step with
-potentially some new wet/dry cells due to the ice sheet dynamics/thermodynamics.
-The wetting and drying scheme, applied on the restart, is very simple. The 6 different possible cases for the tracer and ssh are:
-
-\begin{description}
-   \item[Thin a cell]:
-   T/S/ssh are unchanged.
-
-   \item[Enlarge  a cell]:
-   See case "Thin a cell down"
-
-   \item[Dry a cell]:
-   Mask, T/S, U/V and ssh are set to 0.
-
-   \item[Wet a cell]: 
-   Mask is set to 1, T/S is extrapolated from neighbours, $ssh_n = ssh_b$.
-   If no neighbours, T/S is extrapolated from old top cell value. 
-   If no neighbours along i,j and k (both previous tests failed), T/S/ssh and mask are set to 0.
-
-   \item[Dry a column]:
-   mask, T/S and ssh are set to 0.
-
-   \item[Wet a column]:
-   set mask to 1, T/S/ssh are extrapolated from neighbours.
-   If no neighbour, T/S/ssh and mask set to 0.
-\end{description}
-
-The method described above will strongly affect the barotropic transport under an ice shelf when the geometry change.
-In order to keep the model stable, an adjustment of the dynamics at the initialisation after the coupling step is needed. 
-The idea behind this is to keep $\pd[\eta]{t}$ as it should be without change in geometry at the initialisation. 
-This will prevent any strong velocity due to large pressure gradient. 
-To do so, we correct the horizontal divergence before $\pd[\eta]{t}$ is computed in the first time step.\\
-
-Furthermore, as the before and now fields are not compatible (modification of the geometry),
-the restart time step is prescribed to be an euler time step instead of a leap frog and $fields_b = fields_n$.\\
-
-The horizontal extrapolation to fill new cell with realistic value is called \np{nn_drown}{nn\_drown} times.
-It means that if the grounding line retreat by more than \np{nn_drown}{nn\_drown} cells between 2 coupling steps,
-the code will be unable to fill all the new wet cells properly and the model is likely to blow up at the initialisation.
-The default number is set up for the MISOMIP idealised experiments.
-This coupling procedure is able to take into account grounding line and calving front migration.
-However, it is a non-conservative proccess. 
-This could lead to a trend in heat/salt content and volume.\\
-
-In order to remove the trend and keep the conservation level as close to 0 as possible,
-a simple conservation scheme is available with \np{ln_isfcpl_cons}{ln\_isfcpl\_cons}\forcode{ = .true.}.
-The heat/salt/vol. gain/loss are diagnosed, as well as the location.
-A correction increment is computed and applied each time step during the model run.
-The corrective increment are applied into the cells itself (if it is a wet cell), the neigbouring cells or the closest wet cell (if the cell is now dry).
+as the ocean water leaving the domain removes heat and salt (at the same concentration) with it.\\
 
 %% =================================================================================================
-\section{Handling of icebergs (ICB)}
-\label{sec:SBC_ICB_icebergs}
-
-\begin{listing}
-  \nlst{namberg}
-  \caption{\forcode{&namberg}}
-  \label{lst:namberg}
-\end{listing}
-
-Icebergs are modelled as lagrangian particles in \NEMO\ \citep{marsh.ivchenko.ea_GMD15}.
-Their physical behaviour is controlled by equations as described in \citet{martin.adcroft_OM10} ).
-(Note that the authors kindly provided a copy of their code to act as a basis for implementation in \NEMO).
-Icebergs are initially spawned into one of ten classes which have specific mass and thickness as
-described in the \nam{berg}{berg} namelist: \np{rn_initial_mass}{rn\_initial\_mass} and \np{rn_initial_thickness}{rn\_initial\_thickness}.
-Each class has an associated scaling (\np{rn_mass_scaling}{rn\_mass\_scaling}),
-which is an integer representing how many icebergs of this class are being described as one lagrangian point
-(this reduces the numerical problem of tracking every single iceberg).
-They are enabled by setting \np[=.true.]{ln_icebergs}{ln\_icebergs}.
-
-Two initialisation schemes are possible.
-\begin{description}
-\item [{\np{nn_test_icebergs}{nn\_test\_icebergs}~$>$~0}] In this scheme, the value of \np{nn_test_icebergs}{nn\_test\_icebergs} represents the class of iceberg to generate
-  (so between 1 and 10), and \np{nn_test_icebergs}{nn\_test\_icebergs} provides a lon/lat box in the domain at each grid point of
-  which an iceberg is generated at the beginning of the run.
-  (Note that this happens each time the timestep equals \np{nn_nit000}{nn\_nit000}.)
-  \np{nn_test_icebergs}{nn\_test\_icebergs} is defined by four numbers in \np{nn_test_box}{nn\_test\_box} representing the corners of
-  the geographical box: lonmin,lonmax,latmin,latmax
-\item [{\np[=-1]{nn_test_icebergs}{nn\_test\_icebergs}}] In this scheme, the model reads a calving file supplied in the \np{sn_icb}{sn\_icb} parameter.
-  This should be a file with a field on the configuration grid (typically ORCA)
-  representing ice accumulation rate at each model point.
-  These should be ocean points adjacent to land where icebergs are known to calve.
-  Most points in this input grid are going to have value zero.
-  When the model runs, ice is accumulated at each grid point which has a non-zero source term.
-  At each time step, a test is performed to see if there is enough ice mass to
-  calve an iceberg of each class in order (1 to 10).
-  Note that this is the initial mass multiplied by the number each particle represents (\ie\ the scaling).
-  If there is enough ice, a new iceberg is spawned and the total available ice reduced accordingly.
-\end{description}
-
-Icebergs are influenced by wind, waves and currents, bottom melt and erosion.
-The latter act to disintegrate the iceberg.
-This is either all melted freshwater,
-or (if \np{rn_bits_erosion_fraction}{rn\_bits\_erosion\_fraction}~$>$~0) into melt and additionally small ice bits
-which are assumed to propagate with their larger parent and thus delay fluxing into the ocean.
-Melt water (and other variables on the configuration grid) are written into the main \NEMO\ model output files.
-
-By default, iceberg thermodynamic and dynamic are computed using ocean surface variable (sst, ssu, ssv) and the icebergs are not sensible to the bathymetry (only to land) whatever the iceberg draft. 
-\citet{Merino_OM2016} developed an option to use vertical profiles of ocean currents and temperature instead (\np{ln_M2016}{ln\_M2016}).
-Full details on the sensitivity to this parameter in done in \citet{Merino_OM2016}. 
-If \np{ln_M2016}{ln\_M2016} activated, \np{ln_icb_grd}{ln\_icb\_grd} activate (or not) an option to prevent thick icebergs to move across shallow bank (ie shallower than the iceberg draft).
-This option need to be used with care as it could required to either change the distribution to prevent generation of icebergs with draft larger than the bathymetry 
-or to build a variable \forcode{maxclass} to prevent NEMO filling the icebergs classes too thick for the local bathymetry.
-
-Extensive diagnostics can be produced.
-Separate output files are maintained for human-readable iceberg information.
-A separate file is produced for each processor (independent of \np{ln_ctl}{ln\_ctl}).
-The amount of information is controlled by two integer parameters:
-\begin{description}
-\item [{\np{nn_verbose_level}{nn\_verbose\_level}}] takes a value between one and four and
-  represents an increasing number of points in the code at which variables are written,
-  and an increasing level of obscurity.
-\item [{\np{nn_verbose_write}{nn\_verbose\_write}}] is the number of timesteps between writes
-\end{description}
-
-Iceberg trajectories can also be written out and this is enabled by setting \np{nn_sample_rate}{nn\_sample\_rate}~$>$~0.
-A non-zero value represents how many timesteps between writes of information into the output file.
-These output files are in NETCDF format.
-When running with multiple processes, each output file contains only those icebergs in the corresponding processor.
-Trajectory points are written out in the order of their parent iceberg in the model's "linked list" of icebergs.
-So care is needed to recreate data for individual icebergs,
-since its trajectory data may be spread across multiple files.
-
-%% =================================================================================================
-\section[Interactions with waves (\textit{sbcwave.F90}, \forcode{ln_wave})]{Interactions with waves (\protect\mdl{sbcwave}, \protect\np{ln_wave}{ln\_wave})}
+\section[Interactions with waves (\textit{sbcwave.F90})]{Interactions with waves (\protect\mdl{sbcwave})}
 \label{sec:SBC_wave}
 
 \begin{listing}
@@ -1705,7 +1310,7 @@ since its trajectory data may be spread across multiple files.
 
 Ocean waves represent the interface between the ocean and the atmosphere, so \NEMO\ is extended to incorporate
 physical processes related to ocean surface waves, namely the surface stress modified by growth and
-dissipation of the oceanic wave field, the Stokes-Coriolis force, the vortex-force, the Bernoulli head J term and the Stokes drift impact on mass and tracer advection; moreover the neutral surface drag coefficient or the Charnock parameter from a wave model can be used to evaluate the wind stress. NEMO has also been extended to take into account the impact of surface waves on the vertical mixing, via the parameterization of the Langmuir turbulence, the modification of the surface boundary conditions for the Turbulent Kinetic Energy closure scheme, and the inclusion the Stokes drift contribution to the shear production term in different turbulent closure schemes (RIC, TKE, GLS).\\
+dissipation of the oceanic wave field, the Stokes-Coriolis force, the vortex-force, the Bernoulli head J term and the Stokes drift impact on mass and tracer advection; moreover the neutral surface drag coefficient or the Charnock parameter from a wave model can be used to evaluate the wind stress. NEMO has also been extended to take into account the impact of surface waves on the vertical mixing, via the parameterization of the Langmuir turbulence, the modification of the surface boundary conditions for the Turbulent Kinetic Energy closure scheme, and the inclusion of the Stokes drift contribution to the shear production term in different turbulent closure schemes (TKE, GLS).\\
 
 Physical processes related to ocean surface waves can be accounted by setting the logical variable
 \np[=.true.]{ln_wave}{ln\_wave} in \nam{sbc}{sbc} namelist. In addition, specific flags accounting for
@@ -1722,6 +1327,7 @@ Input Data generic Interface (see \autoref{sec:SBC_input}).
 in \nam{sbc}{sbc} namelist and filling the \nam{sbc_cpl}{sbc\_cpl} namelist. NEMO can receive the significant wave height, the mean wave period ($T0M1$), the mean wavenumber, the Charnock parameter, the neutral drag coefficient, the two components of the surface Stokes drift and the associated transport, the wave to ocean momentum flux, the net wave-supported stress, the Bernoulli head $J$ term and the wave to ocean energy flux term.
 \end{description}
 
+The option \np[=.true.]{ln_wave_test}{ln\_wave\_test} enables testing of ocean-wave interactions using an idealized constant wave field. The fields define are the surface Stokes drift (in the x-direction), as well as the significant wave height and mean wave period.
 
 %% =================================================================================================
 \subsection[Neutral drag coefficient from wave model (\forcode{ln_cdgw})]{Neutral drag coefficient from wave model (\protect\np{ln_cdgw}{ln\_cdgw})}
@@ -1729,9 +1335,9 @@ in \nam{sbc}{sbc} namelist and filling the \nam{sbc_cpl}{sbc\_cpl} namelist. NEM
 
 The neutral surface drag coefficient provided from an external data source (\ie\ forced or coupled wave model),
 can be used by setting the logical variable \np[=.true.]{ln_cdgw}{ln\_cdgw} in \nam{sbc_wave}{sbc\_wave} namelist.
-Then using the routine \rou{sbcblk\_algo\_ncar} and starting from the neutral drag coefficient provided,
-the drag coefficient is computed according to the stable/unstable conditions of the
-air-sea interface following \citet{large.yeager_trpt04}.
+The drag coefficient is computed according to the stable/unstable conditions of the
+air-sea interface following \citet{large.yeager_trpt04}, starting from the neutral drag coefficient provided.
+This option is only available for the \np{ln_NCAR}{ln\_NCAR} and \np{ln_MFS}{ln\_MFS} bulk formulae.
 %% =================================================================================================
 
 
@@ -1739,7 +1345,7 @@ air-sea interface following \citet{large.yeager_trpt04}.
 \label{subsec:SBC_wave_charn}
 
 The dimensionless Charnock parameter characterising the sea surface roughness provided from an external wave model, can be used by setting the logical variable \np[=.true.]{ln_charn}{ln\_charn} in \nam{sbc_wave}{sbc\_wave} namelist. Then using the routine \rou{sbcblk\_algo\_ecmwf}, the roughness length that enters the definition of the drag coefficient is computed according to the Charnock parameter depending on the sea state.
-Note that this option is only available in coupled mode.\\
+Note that this option is only available in coupled mode and for \np[=.true.]{ln_ECMWF}{ln\_ECMWF}.\\
 
 %% =================================================================================================
 
@@ -1776,7 +1382,7 @@ Two possible parameterizations for the calculation for the approximate Stokes dr
 are included in the code once provided the surface Stokes drift $\mathbf{U}_{st |_{z=0}}$ which is evaluated by an external wave model that accurately reproduces the wave spectra and makes possible the estimation of the surface Stokes drift for random directional waves in realistic wave conditions. To evaluate the 3D Stokes drift, the surface Stokes drift (zonal and meridional components), the Stokes transport or alternatively the significant wave height and the mean wave period should be provided either in forced or coupled mode.
 
 \begin{description}
-\item [By default (\forcode{ln_breivikFV_2016=.true.})]:\\ 
+\item By default \np[=.false.]{ln_breivikFV_2016}{ln\_breivikFV\_2016}:\\
 An exponential integral profile parameterization proposed by \citet{breivik.janssen.ea_JPO14} is used:
 
 \[
@@ -1795,7 +1401,7 @@ where $k_e$ is the effective wave number which depends on the Stokes transport $
 
 where $H_s$ is the significant wave height and $\bar{\omega}$ is the wave frequency defined as: $\bar{\omega}=\frac{2\pi}{T_m}$ (being $T_m$ the mean wave period).
 
-\item [If \forcode{ln_breivikFV_2016= .true.} ]: \\
+\item If \np[=.true.]{ln_breivikFV_2016}{ln\_breivikFV\_2016}: \\
 
 A velocity profile based on the Phillips spectrum which is considered to be a reasonable estimate of the part of the spectrum mostly contributing to the Stokes drift velocity near the surface \citep{breivik.bidlot.ea_OM16} is used:
 
@@ -1956,7 +1562,7 @@ Furthermore, only the knowledge of daily mean value of SWF is needed,
 as higher frequency variations can be reconstructed from them,
 assuming that the diurnal cycle of SWF is a scaling of the top of the atmosphere diurnal cycle of incident SWF.
 The \cite{bernie.guilyardi.ea_CD07} reconstruction algorithm is available in \NEMO\ by
-setting \np[=.true.]{ln_dm2dc}{ln\_dm2dc} (a \textit{\nam{sbc}{sbc}} namelist variable) when
+setting \np[=.true.]{ln_dm2dc}{ln\_dm2dc} (a \nam{sbc}{sbc} namelist variable) when
 using a bulk formulation (\np[=.true.]{ln_blk}{ln\_blk}) or
 the flux formulation (\np[=.true.]{ln_flx}{ln\_flx}).
 The reconstruction is performed in the \mdl{sbcdcy} module.
@@ -2011,20 +1617,23 @@ The rot\_rep routine from the \mdl{geo2ocean} module is used to perform the rota
   \label{lst:namsbc_ssr}
 \end{listing}
 
+The addition of a surface restoring term to observed SST and/or SSS can be activated defining \np[=.true.]{ln_ssr}{ln\_ssr} in \nam{sbc}{sbc}.
 Options are defined through the \nam{sbc_ssr}{sbc\_ssr} namelist variables.
 On forced mode using a flux formulation (\np[=.true.]{ln_flx}{ln\_flx}),
-a feedback term \emph{must} be added to the surface heat flux $Q_{ns}^o$:
-\[
-  % \label{eq:SBC_dmp_q}
+a feedback term \emph{must} be added to the surface heat flux $Q_{ns}^o$ (\np[= 1]{nn_sstr}{nn\_sstr}):
+
+\begin{equation}
+  \label{eq:SBC_dmp_q}
   Q_{ns} = Q_{ns}^o + \frac{dQ}{dT} \left( \left. T \right|_{k=1} - SST_{Obs} \right)
-\]
-where SST is a sea surface temperature field (observed or climatological),
+\end{equation}
+
+where SST is a sea surface temperature field (observed or climatological, \np{sn_sst}{sn\_sst}),
 $T$ is the model surface layer temperature and
-$\frac{dQ}{dT}$ is a negative feedback coefficient usually taken equal to $-40~W/m^2/K$.
+$\frac{dQ}{dT}$  (\np{rn_dqdt}{rn\_dqdt}) is a negative feedback coefficient usually taken equal to $-40~W/m^2/K$.
 For a $50~m$ mixed-layer depth, this value corresponds to a relaxation time scale of two months.
 This term ensures that if $T$ perfectly matches the supplied SST, then $Q$ is equal to $Q_o$.
 
-In the fresh water budget, a feedback term can also be added.
+In the fresh water budget, a feedback term can also be added (\np[= 1]{nn_sssr}{nn\_sssr}).
 Converted into an equivalent freshwater flux, it takes the following expression :
 
 \begin{equation}
@@ -2033,12 +1642,20 @@ Converted into an equivalent freshwater flux, it takes the following expression
   {\left.S\right|_{k=1}}
 \end{equation}
 
-where $\textit{emp}_{o }$ is a net surface fresh water flux
+where $\textit{emp}_{o}$ is a net surface fresh water flux
 (observed, climatological or an atmospheric model product),
-\textit{SSS}$_{Obs}$ is a sea surface salinity
-(usually a time interpolation of the monthly mean Polar Hydrographic Climatology \citep{steele.morley.ea_JC01}),
+$\textit{SSS}_{Obs}$ is a sea surface salinity estimate ( \np{sn_sss}{sn\_sss}),
 $\left.S\right|_{k=1}$ is the model surface layer salinity and
-$\gamma_s$ is a negative feedback coefficient which is provided as a namelist parameter.
+$\gamma_s$ is a negative feedback coefficient which is provided as \np{rn_deds}{rn\_deds}. At high resolution is is worth bounding the restoring term (\np{ln_sssr_bnd}{ln\_sssr\_bnd}) 
+because it can be extremely high in structured not resolved in the sss observation (filament or eddies). If activated the maximum value of the correction term is set to \np{rn_sssr_bnd}{rn\_sssr\_bnd}.
+There is also an option to deactivate the salinity restoring under sea ice (\np{nn_sssr_ice}{nn\_sssr\_ice}). There is 3 options:
+\begin{description}
+\item [nn\_sssr\_ice = 0]: no restoring under sea ice. This can be justify by the fact that the observation under the sea ice less reliable as in the open ocean.
+\item [nn\_sssr\_ice = 1]: same restoring under sea ice than within the open ocean
+\item [nn\_sssr\_ice > 1]: strong restoring under the ice. In this case, the user trust more the sss observation that its sea ice model as the restoring will tend to strongly constrain the surface salinity that in reality is strongly driven by the sea ice melt / formation in these area.
+\end{description}
+If the sea surface salinity estimates do not resolve the river fresh water plumes, the restoring will tend to counter balance the effect of the surface runoff in the model. To avoid such effect, the salinity restoring is not applied at the river mouth mask defined when the runoff parameter \np{ln_rnf_mouth}{ln\_rnf\_mouth} is activated.
+
 Unlike heat flux, there is no physical justification for the feedback term in \autoref{eq:SBC_dmp_emp} as
 the atmosphere does not care about ocean surface salinity \citep{madec.delecluse_IWN97}.
 The SSS restoring term should be viewed as a flux correction on freshwater fluxes to
@@ -2088,12 +1705,11 @@ the value of the \np{nn_ice}{nn\_ice} namelist parameter found in \nam{sbc}{sbc}
 
 For global ocean simulations, it can be useful to introduce a control of the
 mean sea level in order to prevent unrealistic drifting of the sea surface
-height due to unbalanced freshwater fluxes. In \NEMO, two options for
-controlling the freshwater budget are proposed.
+height due to unbalanced freshwater fluxes. In \NEMO, four options for
+controlling the freshwater budget are proposed:
 
 \begin{description}
-\item [{\np[=0]{nn_fwb}{nn\_fwb}}:] No control at all; the mean sea level is
-  free to drift, and will certainly do so.
+\item [{\np[=0]{nn_fwb}{nn\_fwb}}:] No control at all; the mean sea level is free to drift, and will certainly do so.
 \item [{\np[=1]{nn_fwb}{nn\_fwb}}:] The global mean \textit{emp} is set to zero at each model time step.
   %GS: comment below still relevant ?
   %Note that with a sea-ice model, this technique only controls the mean sea level with linear free surface and no mass flux between ocean and ice (as it is implemented in the current ice-ocean coupling).
@@ -2103,7 +1719,22 @@ controlling the freshwater budget are proposed.
   approximation, the freshwater budget can be evaluated from the change in the
   mean sea level and in the ice and snow mass after the end of each simulation
   year; at the start of the model run, an initial adjustment flux can be set
-  using parameter \np{rn_rwb0}{rn\_fwb0} in namelist \nam{sbc_fwb}{sbc\_fwb}.
+  using parameter \np{rn_fwb0}{rn\_fwb0} in namelist \nam{sbc_fwb}{sbc\_fwb}.
+\item [{\np[=3]{nn_fwb}{nn\_fwb}}:] volume adjusted at each time step and spread out over >0 erp (sea surface salinity restoring term) area to increase evaporation or 
+  spread out over <0 erp area to increase precipitation.
+\item [{\np[=4]{nn_fwb}{nn\_fwb}}:] volume adjusted at each time step with a correction on non solar heat flux (\textit{qns}) and salt flux (\textit{sflx})
+  to avoid any surface buoyancy flux associated with the emp correction.
+\end{description}
+
+If \np{nn_fwb}{nn\_fwb} is set > 0, \nam{sbc_fwb}{sbc\_fwb} block must be filled accordingly:
+\begin{description}
+  \item [rn\_fwb0 :] if {\np[=2]{nn_fwb}{nn\_fwb}}, it defines the initial freshwater adjustment flux.
+  \item [nn\_fwb\_voltype :] it refers to the variable considered as the global value of volume (in equivalent liquid height in m) to be conserved by the adjustment process.
+  \item [nn\_fwb\_voltype :] if set to 1, the total ocean+equivalent liquid sea ice volume water budget is controled, or if set to 0, only the ocean volume is controled.
+    The former is now the default and recommended, being obviously more accurate on a physical point of view.
+  \item [ln\_hvolg\_var :] if set to .true., an analytical variation of the global liquid height can be specified by the user.
+    It is defined as the sum of an annual harmonic signal (with a peak to peak amplitude \np{rn_hvolg_amp}{rn\_hvolg\_amp} in m,
+    and a zero crossing at the beginning of month \np{nn_hvolg_mth}{nn\_hvolg\_mth}) and a linear trend given by \np{rn_hvolg_trd}{rn\_hvolg\_trd} (in m/s).
 \end{description}
 
 % Griffies doc:
@@ -2116,7 +1747,7 @@ controlling the freshwater budget are proposed.
 % How often the normalization is done is a matter of choice. In mom4p1, we choose to do so at each model time step,
 % so that there is always a zero net input of water to the ocean-ice system.
 % Others choose to normalize over an annual cycle, in which case the net imbalance over an annual cycle is used
-% to alter the subsequent year�s water budget in an attempt to damp the annual water imbalance.
+% to alter the subsequent years water budget in an attempt to damp the annual water imbalance.
 % Note that the annual budget approach may be inappropriate with interannually varying precipitation forcing.
 % When running ocean-ice coupled models, it is incorrect to include the water transport between the ocean
 % and ice models when aiming to balance the hydrological cycle.