Index: /NEMO/trunk/doc/latex/NEMO/main/bibliography.bib
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/main/bibliography.bib	(revision 15521)
+++ /NEMO/trunk/doc/latex/NEMO/main/bibliography.bib	(revision 15522)
@@ -563,4 +563,15 @@
 }
 
+@article{        couvelard_2020,
+  author        = "X. Couvelard and F. Lemari{\'e} and G. Samson and J.-L. Redelsperger and F. Ardhuin and R. Benshila and G. Madec",
+  doi           = "10.5194/gmd-13-3067-2020",
+  journal       = "Geosci. Model Dev",
+  month         = "Jul",
+  pages         = "3067--3090",
+  title         = "Development of a two-way-coupled ocean--wave model: assessment on a global NEMO(v3.6)--WW3(v6.02) coupled configuration",
+  volume        = "13",
+  year          = "2020",
+}
+
 @article{         cox_OM87,
   title         = "Isopycnal diffusion in a z-coordinate ocean model",
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 15521)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 15522)
@@ -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.
 
 %% =================================================================================================
Index: /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex
===================================================================
--- /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex	(revision 15521)
+++ /NEMO/trunk/doc/latex/NEMO/subfiles/chap_ZDF.tex	(revision 15522)
@@ -278,5 +278,6 @@
 
 %% =================================================================================================
-\subsubsection{Surface wave breaking parameterization}
+\subsubsection{Surface wave breaking parameterization (No information from an external wave model)}
+\label{subsubsec:ZDF_tke_wave} 
 
 Following \citet{mellor.blumberg_JPO04}, the TKE turbulence closure model has been modified to
@@ -306,11 +307,34 @@
 with $e_{bb}$ the \np{rn_ebb}{rn\_ebb} namelist parameter, setting \np[=67.83]{rn_ebb}{rn\_ebb} corresponds
 to $\alpha_{CB} = 100$.
-Further setting  \np[=.true.]{ln_mxl0}{ln\_mxl0},  applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale,
-with $\beta$ hard coded to the Stacey's value.
-Note that a minimal threshold of \np{rn_emin0}{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the
-surface $\bar{e}$ value.
+
+Further setting  \np[=.true.]{ln_mxl0}{ln\_mxl0},  applies \autoref{eq:ZDF_Lsbc} as the surface boundary condition on the length scale, with $\beta$ hard coded to the Stacey's value. Note that a minimal threshold of \np{rn_emin0}{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on the surface $\bar{e}$ value.\\
+
+\subsubsection{Surface wave breaking parameterization (using information from an external wave model)}
+\label{subsubsec:ZDF_tke_waveco} 
+
+Surface boundary conditions for the turbulent kinetic energy, the mixing length scale and the dissipative length scale can be defined using wave fields provided from an external wave model (see \autoref{chap:SBC}, \autoref{sec:SBC_wave}). 
+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 ($\Phi_o$) from an external wave model can be provided and then converted into an ocean turbulence source by setting ln\_phioc=.true.
+
+The surface TKE can be defined by a Dirichlet boundary condition setting $nn\_bc\_surf=0$ in \nam{zdf}{tke} namelist:
+\begin{equation}
+  \bar{e}_o  = \frac{1}{2}\,\left( 15.8 \, \frac{\Phi_o}{\rho_o}\right) ^{2/3}
+\end{equation}
+
+Nevertheless, due to the definition of the computational grid, the TKE flux is not applied at the free surface but at the centre of the topmost grid cell ($z = z1$). To be more accurate, a Neumann boundary condition amounting to interpreter the half-grid cell at the top as a constant flux layer (consistent with the surface layer Monin–Obukhov theory) can be applied setting $nn\_bc\_surf=1$ in  \nam{zdf}{tke} namelist \citep{couvelard_2020}:
+
+\begin{equation}
+  \left(\frac{Km}{e_3}\,\partial_k e \right)_{z=z1} = \frac{\Phi_o}{\rho_o}
+\end{equation}
+
+
+The mixing length scale surface value $l_0$ can be estimated from the surface roughness length z0:
+\begin{equation}
+  l_o = \kappa \, \frac{ \left( C_k\,C_\epsilon \right) ^{1/4}}{C_k}\, z0
+\end{equation}
+where $z0$ is directly estimated from the significant wave height ($Hs$) provided by the external wave model as $z0=1.6Hs$. To use this option ln\_mxhsw as well as ln\_wave and ln\_sdw have to be set to .true.
 
 %% =================================================================================================
 \subsubsection{Langmuir cells}
+\label{subsubsec:ZDF_tke_langmuir}
 
 Langmuir circulations (LC) can be described as ordered large-scale vertical motions in
@@ -335,11 +359,5 @@
 \]
 where $w_{LC}(z)$ is the vertical velocity profile of LC, and $H_{LC}$ is the LC depth.
-With no information about the wave field, $w_{LC}$ is assumed to be proportional to
-the Stokes drift $u_s = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module
-\footnote{Following \citet{li.garrett_JMR93}, the surface Stoke drift velocity may be expressed as
-  $u_s =  0.016 \,|U_{10m}|$.
-  Assuming an air density of $\rho_a=1.22 \,Kg/m^3$ and a drag coefficient of
-  $1.5~10^{-3}$ give the expression used of $u_s$ as a function of the module of surface stress
-}.
+
 For the vertical variation, $w_{LC}$ is assumed to be zero at the surface as well as at
 a finite depth $H_{LC}$ (which is often close to the mixed layer depth),
@@ -349,10 +367,23 @@
   w_{LC}  =
   \begin{cases}
-    c_{LC} \,u_s \,\sin(- \pi\,z / H_{LC} )    &      \text{if $-z \leq H_{LC}$} 	\\
+    c_{LC} \,\|u_s^{LC}\| \,\sin(- \pi\,z / H_{LC} )    &      \text{if $-z \leq H_{LC}$} 	\\
     0                 				 &      \text{otherwise}
   \end{cases}
 \]
-where $c_{LC} = 0.15$ has been chosen by \citep{axell_JGR02} as a good compromise to fit LES data.
-The chosen value yields maximum vertical velocities $w_{LC}$ of the order of a few centimeters per second.
+
+
+In the absence of information about the wave field, $w_{LC}$ is assumed to be proportional to
+the surface Stokes drift ($u_s^{LC}=u_{s0} $) empirically estimated by $ u_{s0} = 0.377\,\,|\tau|^{1/2}$, where $|\tau|$ is the surface wind stress module
+\footnote{Following \citet{li.garrett_JMR93}, the surface Stoke drift velocity may be expressed as
+  $u_{s0} =  0.016 \,|U_{10m}|$.
+  Assuming an air density of $\rho_a=1.22 \,Kg/m^3$ and a drag coefficient of
+  $1.5~10^{-3}$ give the expression used of $u_{s0}$ as a function of the module of surface stress
+}.
+
+In case of online coupling with an external wave model (see \autoref{chap:SBC} \autoref{sec:SBC_wave}), $w_{LC}$ is proportional to the component of the Stokes drift aligned with the wind \citep{couvelard_2020} and $ u_s^{LC}  = \max(u_{s0}.e_\tau,0)$ where $e_\tau$ is the unit vector in the wind stress direction and $u_{s0}$ is the surface Stokes drift provided by the external wave model.
+
+
+$c_{LC} = 0.15$ has been chosen by \citep{axell_JGR02} as a good compromise to fit LES data.
+The chosen value yields maximum vertical velocities $w_{LC}$ of the order of a few centimetres per second.
 The value of $c_{LC}$ is set through the \np{rn_lc}{rn\_lc} namelist parameter,
 having in mind that it should stay between 0.15 and 0.54 \citep{axell_JGR02}.
@@ -362,5 +393,5 @@
 converting its kinetic energy to potential energy, according to
 \[
-- \int_{-H_{LC}}^0 { N^2\;z  \;dz} = \frac{1}{2} u_s^2
+- \int_{-H_{LC}}^0 { N^2\;z  \;dz} = \frac{1}{2} \|u_s^{LC}\|^2
 \]
 
@@ -1427,5 +1458,5 @@
 the Stokes Drift can be evaluated by setting \forcode{ln_sdw=.true.}
 (see \autoref{subsec:SBC_wave_sdw})
-and the needed wave fields can be provided either in forcing or coupled mode
+and the needed wave fields (significant wave height and mean wave number) can be provided either in forcing or coupled mode
 (for more information on wave parameters and settings see \autoref{sec:SBC_wave})
 
Index: /NEMO/trunk/doc/namelists/namsbc
===================================================================
--- /NEMO/trunk/doc/namelists/namsbc	(revision 15521)
+++ /NEMO/trunk/doc/namelists/namsbc	(revision 15522)
@@ -33,13 +33,4 @@
    ln_isf      = .false.   !  ice shelf                                 (T   => fill namsbc_isf & namsbc_iscpl)
    ln_wave     = .false.   !  Activate coupling with wave  (T => fill namsbc_wave)
-   ln_cdgw     = .false.   !  Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave)
-   ln_sdw      = .false.   !  Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave) 
-   nn_sdrift   =  0        !  Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift
-      !                    !   = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]
-      !                    !   = 1 Phillips:                      v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))]
-      !                    !   = 2 Phillips as (1) but using the wave frequency from a wave model
-   ln_tauwoc   = .false.   !  Activate ocean stress modified by external wave induced stress (T => ln_wave=.true. & fill namsbc_wave)
-   ln_tauw     = .false.   !  Activate ocean stress components from wave model
-   ln_stcor    = .false.   !  Activate Stokes Coriolis term (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave)
    nn_lsm      = 0         !  =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) ,
                            !  =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field)
Index: /NEMO/trunk/doc/namelists/namsbc_cpl
===================================================================
--- /NEMO/trunk/doc/namelists/namsbc_cpl	(revision 15521)
+++ /NEMO/trunk/doc/namelists/namsbc_cpl	(revision 15522)
@@ -48,6 +48,9 @@
    sn_rcv_isf    =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
    sn_rcv_icb    =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
-   sn_rcv_tauwoc =   'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_tauw   =   'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
-   sn_rcv_wdrag  =   'none'                 ,    'no'    ,     ''      ,         ''          ,   ''
+   sn_rcv_wdrag  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_charn  =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_taw    =   'none'                 ,    'no'    ,     ''      ,         ''           ,   'U,V'
+   sn_rcv_bhd    =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_tusd   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
+   sn_rcv_tvsd   =   'none'                 ,    'no'    ,     ''      ,         ''           ,   ''
 /
Index: /NEMO/trunk/doc/namelists/namsbc_wave
===================================================================
--- /NEMO/trunk/doc/namelists/namsbc_wave	(revision 15521)
+++ /NEMO/trunk/doc/namelists/namsbc_wave	(revision 15522)
@@ -2,4 +2,18 @@
 &namsbc_wave   ! External fields from wave model                        (ln_wave=T)
 !-----------------------------------------------------------------------
+   ln_sdw      = .false.       !  get the 2D Surf Stokes Drift & Compute the 3D stokes drift
+   ln_stcor    = .false.       !  add Stokes Coriolis and tracer advection terms
+   ln_cdgw     = .false.       !  Neutral drag coefficient read from wave model
+   ln_tauoc    = .false.       !  ocean stress is modified by wave induced stress
+   ln_wave_test= .false.       !  Test case with constant wave fields
+!
+   ln_charn    = .false.       !  Charnock coefficient read from wave model (IFS only)
+   ln_taw      = .false.       !  ocean stress is modified by wave induced stress (coupled mode)
+   ln_phioc    = .false.       !  TKE flux from wave model
+   ln_bern_srfc= .false.       !  wave induced pressure. Bernoulli head J term
+   ln_breivikFV_2016 = .false. !  breivik 2016 vertical stokes profile
+   ln_vortex_force = .false.   !  Vortex Force term
+   ln_stshear  = .false.       !  include stokes shear in EKE computation
+!
    cn_dir      = './'      !  root directory for the waves data location
    !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________!
@@ -11,8 +25,5 @@
    sn_hsw      =  'sdw_ecwaves_orca2'    ,        6.         , 'hs'         ,  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
    sn_wmp      =  'sdw_ecwaves_orca2'    ,        6.         , 'wmp'        ,  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
-   sn_wfr      =  'sdw_ecwaves_orca2'    ,        6.         , 'wfr'        ,  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
    sn_wnum     =  'sdw_ecwaves_orca2'    ,        6.         , 'wave_num'   ,  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
-   sn_tauwoc   =  'sdw_ecwaves_orca2'    ,        6.         , 'wave_stress',  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
-   sn_tauwx    =  'sdw_ecwaves_orca2'    ,        6.         , 'wave_stress',  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
-   sn_tauwy    =  'sdw_ecwaves_orca2'    ,        6.         , 'wave_stress',  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
+   sn_tauoc    =  'sdw_ecwaves_orca2'    ,        6.         , 'wave_stress',  .true.  , .true. , 'yearly'  ,  ''              , ''       , ''
 /
Index: /NEMO/trunk/doc/namelists/namzdf_tke
===================================================================
--- /NEMO/trunk/doc/namelists/namzdf_tke	(revision 15521)
+++ /NEMO/trunk/doc/namelists/namzdf_tke	(revision 15522)
@@ -30,3 +30,5 @@
    !                       !           = 2 weighted by 1-fr_i
    !                       !           = 3 weighted by 1-MIN(1,4*fr_i)
+   nn_bc_surf   =     1    !  surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling (ln_cplwave=1)
+   nn_bc_bot    =     1    !  bottom condition (0/1=Dir/Neum) ! Only applicable for wave coupling (ln_cplwave=1)
 /
