Index: NEMO/branches/2021/dev_r14886_VLD-03_Aimie_Moulin_Wave_Coupling_TestCase/doc/latex/NEMO/subfiles/chap_SBC.tex
===================================================================
--- NEMO/branches/2021/dev_r14886_VLD-03_Aimie_Moulin_Wave_Coupling_TestCase/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 14896)
+++ NEMO/branches/2021/dev_r14886_VLD-03_Aimie_Moulin_Wave_Coupling_TestCase/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 14968)
@@ -78,11 +78,5 @@
   (\np[=.true.]{ln_dm2dc}{ln\_dm2dc}),
 \item the activation of wave effects from an external wave model  (\np[=.true.]{ln_wave}{ln\_wave}),
-\item a neutral drag coefficient is read from an external wave model (\np[=.true.]{ln_cdgw}{ln\_cdgw}),
-\item the Stokes drift from an external wave model is accounted for (\np[=.true.]{ln_sdw}{ln\_sdw}),
-\item the choice of the Stokes drift profile parameterization (\np[=0..2]{nn_sdrift}{nn\_sdrift}),
-\item the surface stress given to the ocean is modified by surface waves (\np[=.true.]{ln_tauwoc}{ln\_tauwoc}),
-\item the surface stress given to the ocean is read from an external wave model (\np[=.true.]{ln_tauw}{ln\_tauw}),
-\item the Stokes-Coriolis term is included (\np[=.true.]{ln_stcor}{ln\_stcor}),
-\item the light penetration in the ocean (\np[=.true.]{ln_traqsr}{ln\_traqsr} with namelist \nam{tra_qsr}{tra\_qsr}),
+\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}),
 \item the effect of sea-ice pressure on the ocean (\np[=.true.]{ln_ice_embd}{ln\_ice\_embd}).
@@ -171,4 +165,5 @@
 
 %\colorbox{yellow}{Penser a} mettre dans le restant l'info nn\_fsbc ET nn\_fsbc*rdt de sorte de reinitialiser la moyenne si on change la frequence ou le pdt
+
 
 %% =================================================================================================
@@ -720,4 +715,6 @@
 respectively (found in \textit{sbcblk\_skin\_ecmwf.F90}).
 
+
+
 \subsubsection{COARE 3.x}
 
@@ -906,4 +903,8 @@
 In this case, CO$_2$ fluxes will be exchanged between the atmosphere and the ice-ocean system
 (and need to be activated in \nam{sbc_cpl}{sbc\_cpl} ).
+
+
+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
@@ -1572,20 +1573,21 @@
 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 and the Stokes drift impact on mass and
-tracer advection; moreover the neutral surface drag coefficient from a wave model can be used to evaluate
-the wind stress.
+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.
+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]: wave fields should be defined through the \nam{sbc_wave}{sbc\_wave} namelist
+\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.
+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}
+
 
 %% =================================================================================================
@@ -1593,12 +1595,22 @@
 \label{subsec:SBC_wave_cdgw}
 
-The neutral surface drag coefficient provided from an external data source (\ie\ a wave model),
-can be used by setting the logical variable \np[=.true.]{ln_cdgw}{ln\_cdgw} in \nam{sbc}{sbc} namelist.
-Then using the routine \rou{sbcblk\_algo\_ncar} and starting from the neutral drag coefficent provided,
+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[3D Stokes Drift (\forcode{ln_sdw} \& \forcode{nn_sdrift})]{3D Stokes Drift (\protect\np{ln_sdw}{ln\_sdw} \& \np{nn_sdrift}{nn\_sdrift})}
+%% =================================================================================================
+
+
+\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}
 
@@ -1628,13 +1640,10 @@
 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.
-Three possible parameterizations for the calculation for the approximate Stokes drift velocity profile
-are included in the code through the \np{nn_sdrift}{nn\_sdrift} parameter 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:
+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 [{\np{nn_sdrift}{nn\_sdrift} = 0}]: exponential integral profile parameterization proposed by
-\citet{breivik.janssen.ea_JPO14}:
+\item [By default (\forcode{ln_breivikFV_2016=.true.})]:\\ 
+An exponential integral profile parameterization proposed by \citet{breivik.janssen.ea_JPO14} is used:
 
 \[
@@ -1647,14 +1656,14 @@
 \[
   % \label{eq:SBC_wave_sdw_0b}
-  k_e = \frac{|\mathbf{U}_{\left.st\right|_{z=0}}|} {|T_{st}|}
+  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 $\omega$ is the wave frequency.
-
-\item [{\np{nn_sdrift}{nn\_sdrift} = 1}]: 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}:
+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:
 
 \[
@@ -1664,19 +1673,16 @@
 \]
 
-where $erf$ is the complementary error function and $k_p$ is the peak wavenumber.
-
-\item [{\np{nn_sdrift}{nn\_sdrift} = 2}]: velocity profile based on the Phillips spectrum as for \np{nn_sdrift}{nn\_sdrift} = 1
-but using the wave frequency from a wave model.
+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}$ is considered as follows:
-
-\[
-  % \label{eq:SBC_wave_eta_sdw}
-  \frac{\partial{\eta}}{\partial{t}} =
-  -\nabla_h \int_{-H}^{\eta} (\mathbf{U} + \mathbf{U}_{st}) dz
-\]
+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
@@ -1699,11 +1705,38 @@
 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 3 possible
+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[Wave modified stress (\forcode{ln_tauwoc} \& \forcode{ln_tauw})]{Wave modified sress (\protect\np{ln_tauwoc}{ln\_tauwoc} \& \np{ln_tauw}{ln\_tauw})}
-\label{subsec:SBC_wave_tauw}
+
+%% =================================================================================================
+\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
@@ -1720,24 +1753,46 @@
 \]
 
-where $\tau_a$ is the atmospheric surface stress;
-$\tau_w$ is the atmospheric stress going into the waves defined as:
+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 {\frac{dk}{c_p} (S_{in}+S_{nl}+S_{diss})}
+  \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.
-
-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_tauwoc}{ln\_tauwoc}, or through the zonal and
-meridional stress components by setting \np[=.true.]{ln_tauw}{ln\_tauw}.
+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.
 
 %% =================================================================================================
