chap_SBC.tex

When an external wave model (see \autoref{sec:SBC_wave}) is used in the coupled system, wave parameters, surface currents and sea surface height can be exchanged between both models (and need to be activated in \nam{sbc_cpl}{sbc\_cpl} ).


The namelist above allows control of various aspects of the coupling fields (particularly for vectors) and
now allows for any coupling fields to have multiple sea ice categories (as required by SI3).
When indicating a multi-category coupling field in \nam{sbc_cpl}{sbc\_cpl}, the number of categories will be determined by
the number used in the sea ice model.
In some limited cases, it may be possible to specify single category coupling fields even when
the sea ice model is running with multiple categories -
in this case, the user should examine the code to be sure the assumptions made are satisfactory.
In cases where this is definitely not possible, the model should abort with an error message.

%% =================================================================================================
\section[Atmospheric pressure (\textit{sbcapr.F90})]{Atmospheric pressure (\protect\mdl{sbcapr})}
\label{sec:SBC_apr}

\begin{listing}
  \nlst{namsbc_apr}
  \caption{\forcode{&namsbc_apr}}
  \label{lst:namsbc_apr}
\end{listing}

The optional atmospheric pressure can be used to force ocean and ice dynamics
(\np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn}, \nam{sbc}{sbc} namelist).
The input atmospheric forcing defined via \np{sn_apr}{sn\_apr} structure (\nam{sbc_apr}{sbc\_apr} namelist)
can be interpolated in time to the model time step, and even in space when the interpolation on-the-fly is used.
When used to force the dynamics, the atmospheric pressure is further transformed into
an equivalent inverse barometer sea surface height, $\eta_{ib}$, using:
\[
  % \label{eq:SBC_ssh_ib}
  \eta_{ib} = -  \frac{1}{g\,\rho_o}  \left( P_{atm} - P_o \right)
\]
where $P_{atm}$ is the atmospheric pressure and $P_o$ a reference atmospheric pressure.
A value of $101,000~N/m^2$ is used unless \np{ln_ref_apr}{ln\_ref\_apr} is set to true.
In this case, $P_o$ is set to the value of $P_{atm}$ averaged over the ocean domain,
\ie\ the mean value of $\eta_{ib}$ is kept to zero at all time steps.

The gradient of $\eta_{ib}$ is added to the RHS of the ocean momentum equation (see \mdl{dynspg} for the ocean).
For sea-ice, the sea surface height, $\eta_m$, which is provided to the sea ice model is set to $\eta - \eta_{ib}$
(see \mdl{sbcssr} module).
$\eta_{ib}$ can be written in the output.
This can simplify altimetry data and model comparison as
inverse barometer sea surface height is usually removed from these date prior to their distribution.

When using time-splitting and BDY package for open boundaries conditions,
the equivalent inverse barometer sea surface height $\eta_{ib}$ can be added to BDY ssh data:
\np{ln_apr_obc}{ln\_apr\_obc}  might be set to true.

%% =================================================================================================
\section{Surface tides (TDE)}
\label{sec:SBC_TDE}

\begin{listing}
  \nlst{nam_tide}
  \caption{\forcode{&nam_tide}}
  \label{lst:nam_tide}
\end{listing}

\subsection{Tidal constituents}
Ocean model component TDE provides the common functionality for tidal forcing
and tidal analysis in the model framework. This includes the computation of the gravitational
surface forcing, as well as support for lateral forcing at open boundaries (see
\autoref{subsec:LBC_bdy_tides}) and tidal harmonic analysis \iffalse (see
\autoref{subsec:DIA_diamlr?} and \autoref{subsec:DIA_diadetide?}) \fi . The module is
activated with \np[=.true.]{ln_tide}{ln\_tide} in namelist
\nam{_tide}{\_tide}. It provides the same 34 tidal constituents that are
included in the
\href{https://www.aviso.altimetry.fr/en/data/products/auxiliary-products/global-tide-fes.html}{FES2014
  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
  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
  4.0 and earlier are available}. Constituents to be included in the tidal forcing
(surface and lateral boundaries) are selected by enumerating their respective
names in namelist array \np{sn_tide_cnames}{sn\_tide\_cnames}.\par

\subsection{Surface tidal forcing}
Surface tidal forcing can be represented in the model through an additional
barotropic force in the momentum equation (\autoref{eq:MB_PE_dyn}) such that:
\[
  \frac{\partial {\mathrm {\mathbf U}}_h }{\partial t} = \ldots +g\nabla (\gamma
  \Pi_{eq} + \Pi_{sal})
\]
where $\gamma \Pi_{eq}$ stands for the equilibrium tidal forcing scaled by a spatially
uniform tilt factor $\gamma$, and $\Pi_{sal}$ is an optional
self-attraction and loading term (SAL). These additional terms are enabled when,
in addition to \np[=.true.]{ln_tide}{ln\_tide}, the parameter
\np[=.true.]{ln_tide_pot}{ln\_tide\_pot}.\par

The equilibrium tidal forcing is expressed as a sum over the subset of
constituents listed in \np{sn_tide_cnames}{sn\_tide\_cnames} of
\nam{_tide}{\_tide} (e.g.,
\begin{forlines}
      sn_tide_cnames(1) = 'M2'
      sn_tide_cnames(2) = 'K1'
      sn_tide_cnames(3) = 'S2'
      sn_tide_cnames(4) = 'O1'
\end{forlines}
to select the four tidal constituents of strongest equilibrium tidal
potential). The tidal tilt factor $\gamma = 1 + k - h$ includes the
Love numbers $k$ and $h$ \citep{love_PRSL09}; this factor is
configurable using \np{rn_tide_gamma}{rn\_tide\_gamma} (default value 0.7). Optionally,
when \np[=.true.]{ln_tide_ramp}{ln\_tide\_ramp}, the equilibrium tidal
forcing can be ramped up linearly from zero during the initial
\np{rn_tide_ramp_dt}{rn\_tide\_ramp\_dt} days of the model run.\par

The SAL term should in principle be computed online as it depends on
the model tidal prediction itself (see \citet{arbic.garner.ea_DSR04} for a
discussion about the practical implementation of this term). The complex
calculations involved in such computations, however, are computationally very
expensive. Here, two mutually exclusive simpler variants are available:
amplitudes generated by an external model for oscillatory $\Pi_{sal}$
contributions from each of the selected tidal constituents can be read in
(\np[=.true.]{ln_read_load}{ln\_read\_load}) from the file specified in
\np{cn_tide_load}{cn\_tide\_load} (the variable names are comprised of the
tidal-constituent name and suffixes \forcode{_z1} and \forcode{_z2} for the two
orthogonal components, respectively); alternatively, a ``scalar approximation''
can be used (\np[=.true.]{ln_scal_load}{ln\_scal\_load}), where
\[
  \Pi_{sal} = \beta \eta,
\]
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

%% =================================================================================================
\section[River runoffs (\textit{sbcrnf.F90})]{River runoffs (\protect\mdl{sbcrnf})}
\label{sec:SBC_rnf}

\begin{listing}
  \nlst{namsbc_rnf}
  \caption{\forcode{&namsbc_rnf}}
  \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,
and was combined with an option to increase vertical mixing near the river mouth.

However, with this method numerical and physical problems arise when the top grid cells are of the order of one meter.
This situation is common in coastal modelling and is becoming more common in 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}.}.

As such from V~3.3 onwards it is possible to add river runoff through a non-zero depth,
and for the temperature and salinity of the river to effect the surrounding ocean.
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 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
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).

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

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,
which causes a salt and heat flux out of the model.
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 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.

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

It is also possible for runnoff 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).

%% =================================================================================================
\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})}
\label{sec:SBC_wave}

\begin{listing}
  \nlst{namsbc_wave}
  \caption{\forcode{&namsbc_wave}}
  \label{lst:namsbc_wave}
\end{listing}

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).\\

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
different processes should be activated as explained in the following sections.\\

Wave fields can be provided either in forced or coupled mode:
\begin{description}
\item [forced mode]: the neutral drag coefficient, the two components of the surface Stokes drift, the significant wave height, the mean wave period, the mean wave number, and the normalized
wave stress going into the ocean can be read from external files. Wave fields should be defined through the \nam{sbc_wave}{sbc\_wave} namelist
for external data names, locations, frequency, interpolation and all the miscellanous options allowed by
Input Data generic Interface (see \autoref{sec:SBC_input}).

\item [coupled mode]: \NEMO\ and an external wave model can be coupled by setting \np[=.true.]{ln_cpl}{ln\_cpl}
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}


%% =================================================================================================
\subsection[Neutral drag coefficient from wave model (\forcode{ln_cdgw})]{Neutral drag coefficient from wave model (\protect\np{ln_cdgw}{ln\_cdgw})}
\label{subsec:SBC_wave_cdgw}

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


\subsection[Charnok coefficient from wave model (\forcode{ln_charn})]{ Charnok coefficient from wave model (\protect\np{ln_charn}{ln\_charn})}
\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.\\

%% =================================================================================================


\subsection[3D Stokes Drift (\forcode{ln_sdw})]{3D Stokes Drift (\protect\np{ln_sdw}{ln\_sdw}) }
\label{subsec:SBC_wave_sdw}

The Stokes drift is a wave driven mechanism of mass and momentum transport \citep{stokes_ibk09}.
It is defined as the difference between the average velocity of a fluid parcel (Lagrangian velocity)
and the current measured at a fixed point (Eulerian velocity).
As waves travel, the water particles that make up the waves travel in orbital motions but
without a closed path. Their movement is enhanced at the top of the orbit and slowed slightly
at the bottom, so the result is a net forward motion of water particles, referred to as the Stokes drift.
An accurate evaluation of the Stokes drift and the inclusion of related processes may lead to improved
representation of surface physics in ocean general circulation models. %GS: reference needed
The Stokes drift velocity $\mathbf{U}_{st}$ in deep water can be computed from the wave spectrum and may be written as:

\[
  % \label{eq:SBC_wave_sdw}
  \mathbf{U}_{st} = \frac{16{\pi^3}} {g}
  \int_0^\infty \int_{-\pi}^{\pi} (cos{\theta},sin{\theta}) {f^3}
  \mathrm{S}(f,\theta) \mathrm{e}^{2kz}\,\mathrm{d}\theta {d}f
\]

where: ${\theta}$ is the wave direction, $f$ is the wave intrinsic frequency,
$\mathrm{S}($f$,\theta)$ is the 2D frequency-direction spectrum,
$k$ is the mean wavenumber defined as:
$k=\frac{2\pi}{\lambda}$ (being $\lambda$ the wavelength). \\

In order to evaluate the Stokes drift in a realistic ocean wave field, the wave spectral shape is required
and its computation quickly becomes expensive as the 2D spectrum must be integrated for each vertical level.
To simplify, it is customary to use approximations to the full Stokes profile.
Two possible parameterizations for the calculation for the approximate Stokes drift velocity profile
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.})]:\\ 
An exponential integral profile parameterization proposed by \citet{breivik.janssen.ea_JPO14} is used:

\[
  % \label{eq:SBC_wave_sdw_0a}
  \mathbf{U}_{st} \cong \mathbf{U}_{st |_{z=0}} \frac{\mathrm{e}^{-2k_ez}} {1-8k_ez}
\]

where $k_e$ is the effective wave number which depends on the Stokes transport $T_{st}$ defined as follows:

\[
  % \label{eq:SBC_wave_sdw_0b}
  k_e = \frac{|\mathbf{U}_{\left.st\right|_{z=0}}|} {5.97|T_{st}|}
  \quad \text{and }\
  T_{st} = \frac{1}{16} \bar{\omega} H_s^2
\]

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.} ]: \\

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:

\[
  % \label{eq:SBC_wave_sdw_1}
  \mathbf{U}_{st} \cong \mathbf{U}_{st |_{z=0}} \Big[exp(2k_pz)-\beta \sqrt{-2 \pi k_pz}
  \textit{ erf } \Big(\sqrt{-2 k_pz}\Big)\Big]
\]

where $erf$ is the complementary error function , $ \beta =1$ and $k_p$ is the peak wavenumber defined as:
\[
  % \label{eq:SBC_wave_kp}
  k_p = \frac{|\mathbf{U}_{\left.st\right|_{z=0}}|}{2 |T_{st}| } (1-2 \beta /3) 
\]

$|T_{st}|$ is estimated from integral wave parameters (Hs and Tm) in forced mode and is provided directly from an external wave model in coupled mode. 

\end{description}

The Stokes drift enters the wave-averaged momentum equation, as well as the tracer advection equations
and its effect on the evolution of the sea-surface height ${\eta}$ by including the barotropic Stokes transport horizontal divergence in the term $D$ of Eq.\ref{eq:MB_ssh}

The tracer advection equation is also modified in order for Eulerian ocean models to properly account
for unresolved wave effect. The divergence of the wave tracer flux equals the mean tracer advection
that is induced by the three-dimensional Stokes velocity.
The advective equation for a tracer $c$ combining the effects of the mean current and sea surface waves
can be formulated as follows:

\[
  % \label{eq:SBC_wave_tra_sdw}
  \frac{\partial{c}}{\partial{t}} =
  - (\mathbf{U} + \mathbf{U}_{st}) \cdot \nabla{c}
\]

%% =================================================================================================
\subsection[Stokes-Coriolis term (\forcode{ln_stcor})]{Stokes-Coriolis term (\protect\np{ln_stcor}{ln\_stcor})}
\label{subsec:SBC_wave_stcor}

In a rotating ocean, waves exert a wave-induced stress on the mean ocean circulation which results
in a force equal to $\mathbf{U}_{st}$×$f$, where $f$ is the Coriolis parameter.
This additional force may have impact on the Ekman turning of the surface current.
In order to include this term, once evaluated the Stokes drift (using one of the 2 possible
approximations described in \autoref{subsec:SBC_wave_sdw}),
\np[=.true.]{ln_stcor}{ln\_stcor} has to be set.

%% =================================================================================================

%% =================================================================================================
\subsection[Vortex-force term (\forcode{ln_vortex_force})]{Vortex-force term (\protect\np{ln_vortex_force}{ln\_vortex\_force})}
\label{subsec:SBC_wave_vf}

The vortex-force term arises from the interaction of the mean flow vorticity with the Stokes drift. 
It results in a force equal to $\mathbf{U}_{st}$×$\zeta$, where $\zeta$ is the mean flow vorticity.
In order to include this term, once evaluated the Stokes drift (using one of the 2 possible
approximations described in \autoref{subsec:SBC_wave_sdw}), \np[=.true.]{ln_vortex_force}{ln\_vortex\_force} has to be set.

%% =================================================================================================

%% =================================================================================================
\subsection[Wave-induced pressure term (\forcode{ln_bern_srfc})]{ Wave-induced pressure term (\protect\np{ln_bern_srfc}{ln\_bern\_srfc})}
\label{subsec:SBC_wave_bhd}
An adjustment in the mean pressure arises to accommodate for the presence of waves.
The mean pressure is corrected adding a depth-uniform wave-induced kinematic pressure term named Bernoulli head $J$ term. The Bernoulli head $J$ term is provided to NEMO from an external wave model where it is defined as:
\[
  % \label{eq:SBC_wave_tauw}
  J = g \iint {\frac{k}{sinh(2kd)} S(k,\theta) d\theta dk}
\]
with $d$ the water depth. \\
In order to include this term, \np[=.true.]{ln_bern_srfc}{ln\_bern\_srfc} has to be defined as well as the Stokes drift option (\autoref{subsec:SBC_wave_sdw}) and the coupling with an external wave model (\autoref{sec:SBC_wave}).

%% =================================================================================================


\subsection[Wave modified stress (\forcode{ln_tauoc} \& \forcode{ln_taw})]{Wave modified stress (\protect\np{ln_tauoc}{ln\_tauoc} \& \np{ln_taw}{ln\_taw})}
\label{subsec:SBC_wave_taw}

The surface stress felt by the ocean is the atmospheric stress minus the net stress going
into the waves \citep{janssen.breivik.ea_trpt13}. Therefore, when waves are growing, momentum and energy is spent and is not
available for forcing the mean circulation, while in the opposite case of a decaying sea
state, more momentum is available for forcing the ocean.
Only when the sea state is in equilibrium, the ocean is forced by the atmospheric stress,
but in practice, an equilibrium sea state is a fairly rare event.
So the atmospheric stress felt by the ocean circulation $\tau_{oc,a}$ can be expressed as:

\[
  % \label{eq:SBC_wave_tauoc}
  \tau_{oc,a} = \tau_a - \tau_w
\]

where $\tau_a$ is the atmospheric surface stress; $\tau_w$ is the atmospheric stress going into the waves defined as:

\[
  % \label{eq:SBC_wave_tauw}
  \tau_w = \rho g \int_{0}^{2\pi} \int {\frac{1}{c_p} (S_{in}+S_{nl}+S_{diss})}dkd\theta
\]

%% ∫2π0∫∞0kω(Sin+Sds) dωdθ

where: $c_p$ is the phase speed of the gravity waves,
$S_{in}$, $S_{nl}$ and $S_{diss}$ are three source terms that represent
the physics of ocean waves. The first one, $S_{in}$, describes the generation of ocean waves by wind and therefore represents the momentum and energy transfer from air to ocean waves; the second term $S_{nl}$ denotes
the nonlinear transfer by resonant four-wave interactions; while the third term $S_{diss}$ describes the dissipation of waves by processes such as white-capping, large scale breaking eddy-induced damping. Note that the $S_{nl}$ is not always taken into account for the calculation of the atmospheric stress going into the waves, depending on the external wave model.
The wave stress derived from an external wave model can be provided either through the normalized wave stress into the ocean by setting \np[=.true.]{ln_tauoc}{ln\_tauoc}, or through the zonal and meridional stress components by setting
 \np[=.true.]{ln_taw}{ln\_taw} .  In coupled mode both options can be used while in forced mode only the first option is included.
 
If the normalized wave stress into the ocean ($\widetilde{\tau}$) is provided (\np[=.true.]{ln_tauoc}{ln\_tauoc}) the atmospheric stress felt by the ocean circulation is expressed as:
\[
  % \label{eq:SBC_wave_tauoc}
  \tau_{oc,a} = \tau_a \times \widetilde{\tau}
\]

If  \np[=.true.]{ln_taw}{ln\_taw} , the zonal and meridional stress fields components from the coupled wave model have to be sent directly to u-grid and v-grid through OASIS.


%% =================================================================================================

\subsection[Waves impact vertical mixing  (\forcode{ln_phioc} \& \forcode{ln_stshear})]{Waves impact vertical mixing (\protect\np{ln_phioc}{ln\_phioc} \& \protect\np{ln_stshear}{ln\_stshear})}
\label{subsec:SBC_wave_TKE}


The vortex-force vertical term gives rise to extra terms in the turbulent kinetic energy (TKE) prognostic \citep{couvelard_2020}. The first term corresponds to a modification of the shear production term. 
The Stokes Drift shear contribution can be included, in coupled mode, by setting \np[=.true.]{ln_stshear}{ln\_stshear}.


In addition, waves affect the surface boundary condition for the turbulent kinetic energy, the mixing length scale and the dissipative length scale of the TKE closure scheme.
The injection of turbulent kinetic energy at the surface can be given by the dissipation of the wave field usually dominated by wave breaking.

In coupled mode, the wave to ocean energy flux term from an external wave model ($ \Phi_o $) can be provided to NEMO and then converted into an ocean turbulence source by setting \np[=.true.]{ln_phioc}{ln\_phioc}.
The boundary condition for the turbulent kinetic energy is implemented in the \rou{zdftke} as a Dirichlet or as a Neumann boundary condition (see \autoref{subsubsec:ZDF_tke_waveco}). The boundary condition for the mixing length scale and the dissipative length scale can also account for surface waves (see \autoref{subsubsec:ZDF_tke_waveco})

Some improvements are introduced in the Langmuir turbulence parameterization (see \autoref{chap:ZDF} \autoref{subsubsec:ZDF_tke_langmuir}) if wave coupled mode is activated.

%% =================================================================================================
\section{Miscellaneous options}
\label{sec:SBC_misc}

%% =================================================================================================
\subsection[Diurnal cycle (\textit{sbcdcy.F90})]{Diurnal cycle (\protect\mdl{sbcdcy})}
\label{subsec:SBC_dcy}

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.66\textwidth]{SBC_diurnal}
  \caption[Reconstruction of the diurnal cycle variation of short wave flux]{
    Example of reconstruction of the diurnal cycle variation of short wave flux from
    daily mean values.
    The reconstructed diurnal cycle (black line) is chosen as
    the mean value of the analytical cycle (blue line) over a time step,
    not as the mid time step value of the analytically cycle (red square).
    From \citet{bernie.guilyardi.ea_CD07}.}
  \label{fig:SBC_diurnal}
\end{figure}

\cite{bernie.woolnough.ea_JC05} have shown that to capture 90$\%$ of the diurnal variability of SST requires a vertical resolution in upper ocean of 1~m or better and a temporal resolution of the surface fluxes of 3~h or less.
%Unfortunately high frequency forcing fields are rare, not to say inexistent. GS: not true anymore !
Nevertheless, it is possible to obtain a reasonable diurnal cycle of the SST knowning only short wave flux (SWF) at high frequency \citep{bernie.guilyardi.ea_CD07}.
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
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.
The detail of the algoritm used can be found in the appendix~A of \cite{bernie.guilyardi.ea_CD07}.
The algorithm preserves the daily mean incoming SWF as the reconstructed SWF at
a given time step is the mean value of the analytical cycle over this time step (\autoref{fig:SBC_diurnal}).
The use of diurnal cycle reconstruction requires the input SWF to be daily
(\ie\ a frequency of 24 hours and a time interpolation set to true in \np{sn_qsr}{sn\_qsr} namelist parameter).
Furthermore, it is recommended to have a least 8 surface module time steps per day,
that is  $\rdt \ nn\_fsbc < 10,800~s = 3~h$.
An example of recontructed SWF is given in \autoref{fig:SBC_dcy} for a 12 reconstructed diurnal cycle,
one every 2~hours (from 1am to 11pm).

\begin{figure}[!t]
  \centering
  \includegraphics[width=0.66\textwidth]{SBC_dcy}
  \caption[Reconstruction of the diurnal cycle variation of short wave flux on an ORCA2 grid]{
    Example of reconstruction of the diurnal cycle variation of short wave flux from
    daily mean values on an ORCA2 grid with a time sampling of 2~hours (from 1am to 11pm).
    The display is on (i,j) plane.}
  \label{fig:SBC_dcy}
\end{figure}

Note also that the setting a diurnal cycle in SWF is highly recommended when
the top layer thickness approach 1~m or less, otherwise large error in SST can appear due to
an inconsistency between the scale of the vertical resolution and the forcing acting on that scale.

%% =================================================================================================
\subsection{Rotation of vector pairs onto the model grid directions}
\label{subsec:SBC_rotation}

When using a flux (\np[=.true.]{ln_flx}{ln\_flx}) or bulk (\np[=.true.]{ln_blk}{ln\_blk}) formulation,
pairs of vector components can be rotated from east-north directions onto the local grid directions.
This is particularly useful when interpolation on the fly is used since here any vectors are likely to
be defined relative to a rectilinear grid.
To activate this option, a non-empty string is supplied in the rotation pair column of the relevant namelist.
The eastward component must start with "U" and the northward component with "V".
The remaining characters in the strings are used to identify which pair of components go together.
So for example, strings "U1" and "V1" next to "utau" and "vtau" would pair the wind stress components together and
rotate them on to the model grid directions;
"U2" and "V2" could be used against a second pair of components, and so on.
The extra characters used in the strings are arbitrary.
The rot\_rep routine from the \mdl{geo2ocean} module is used to perform the rotation.

%% =================================================================================================
\subsection[Surface restoring to observed SST and/or SSS (\textit{sbcssr.F90})]{Surface restoring to observed SST and/or SSS (\protect\mdl{sbcssr})}
\label{subsec:SBC_ssr}

\begin{listing}
  \nlst{namsbc_ssr}
  \caption{\forcode{&namsbc_ssr}}
  \label{lst:namsbc_ssr}
\end{listing}

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}
  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),
$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$.
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.
Converted into an equivalent freshwater flux, it takes the following expression :

\begin{equation}
  \label{eq:SBC_dmp_emp}
  \textit{emp} = \textit{emp}_o + \gamma_s^{-1} e_{3t}  \frac{  \left(\left.S\right|_{k=1}-SSS_{Obs}\right)}
  {\left.S\right|_{k=1}}
\end{equation}

where $\textit{emp}_{o }$ is a net surface fresh water flux