Index: branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Biblio/Biblio.bib
===================================================================
--- branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Biblio/Biblio.bib	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Biblio/Biblio.bib	(revision 4644)
@@ -514,4 +514,14 @@
 }
 
+@TECHREPORT{Breivik_ECMWF13,
+      AUTHOR = "{\O} Breivik and PAEM Janssen and JR Bidlot",
+      TITLE = "{Approximate Stokes Drift Profiles in Deep Water}",
+      YEAR = "2013",
+      PAGES = "18",
+      NUMBER = "716",
+      URL = "http://www.ecmwf.int/publications/library/do/references/list/14",
+      TYPE = "ECMWF Technical Memorandum",
+      INSTITUTION = "European Centre for Medium-Range Weather Forecasts"}
+
 @ARTICLE{Brown_Campana_MWR78,
   author = {J. A. Brown and K. A. Campana},
@@ -651,4 +661,13 @@
 }
 
+@article{Charnock_QJRMS55,
+  title="{Wind stress on a water surface}",
+  author="Charnock, H",
+  journal=QJRMS,
+  volume="81",
+  number="350",
+  pages={639--640},
+  year=1955}
+
 @ARTICLE{Cox1987,
   author = {M. Cox},
@@ -742,4 +761,18 @@
   volume = {34},
   pages = {8--13}
+}
+
+@article{Dee_QJRMS11,
+  title="{The ERA-Interim reanalysis: Configuration and performance of the data assimilation system}",
+  author="DP Dee and Uppala, SM and Simmons, AJ and Berrisford, P. and 
+  Poli, P. and Kobayashi, S. and Andrae, U. and Balmaseda, MA and Balsamo, G.
+  and Bauer, P and Bechtold P and Beljaars, ACM and L van de Berg and J Bidlot
+  and N Bormann and others",
+  journal=QJRMS,
+  volume={137},
+  number={656},
+  pages={553--597, doi:10.1002/qj.828},
+  year={2011},
+  publisher={Wiley Online Library}
 }
 
@@ -895,4 +928,17 @@
   url = {http://dx.doi.org/10.1029/2005GL022463}
 }
+
+
+@article{Edson_JPO13,
+  title="{On the Exchange of Momentum over the Open Ocean}",
+  author="Edson, James and Jampana, Venkata and Weller, Robert 
+          and Bigorre, Sebastien and Plueddemann, Albert and Fairall, Christopher 
+          and Miller, Scott and Mahrt, Larry and Vickers, Dean and 
+          Hersbach, Hans",
+  journal=JPO,
+  volume="43",
+  pages = "1589--1610, doi:10.1175/JPO-D-12-0173.1",
+  doi = "10.1175/JPO-D-12-0173.1",
+  year="2013"}
 
 @ARTICLE{Egbert_Ray_JGR01,
@@ -1264,4 +1310,13 @@
 }
 
+@article{Hasselmann_GAFD70,
+  title="{Wave-driven inertial oscillations}",
+  author="Hasselmann, K",
+  journal="Geophysical and Astrophysical Fluid Dynamics",
+  volume="1",
+  number="3-4",
+  pages="463--502, doi:10.1080/03091927009365783",
+  year="1970"}
+
 @ARTICLE{Hazeleger_Drijfhout_JPO98,
   author = {W. Hazeleger and S. S. Drijfhout},
@@ -1381,4 +1436,12 @@
 }
 
+@book{Holthuijsen07,
+  title="{Waves in Oceanic and Coastal Waters}",
+  author={Holthuijsen, L.H.},
+  year={2007},
+  pages={387},
+  location="Books/holthuisjen07.pdf",
+  publisher={Cambridge University Press} }
+
 @ARTICLE{Hordoir_al_CD08,
   author = {R. Hordoir and J. Polcher and J.-C. Brun-Cottan and G. Madec},
@@ -1476,4 +1539,55 @@
   pages = {381--389}
 }
+
+@article{Janssen_JPO89,
+  title="{Wave-induced stress and the drag of air flow over sea waves}",
+  author="Janssen, PAEM",
+  year="1989",
+  journal=JPO,
+  volume={19},
+  number={6},
+  pages={745--754, doi:10/fsz7vd}
+}
+
+@inproceedings{Janssen_AH04,
+  title="{Impact of the sea state on the atmosphere and ocean}",
+  author="Janssen, P.A.E.M. and Saetra, O. and Wettre, C. and 
+           Hersbach, H. and Bidlot, J.",
+  booktitle="Annales hydrographiques",
+  volume={3-772},
+  pages={3.1--3.23},
+  year={2004},
+  organization={Service hydrographique et oc{\'e}anographique de la marine}
+}
+
+@article{Janssen_Rep08,
+  title="{Progress in ocean wave forecasting}",
+  author="Janssen, PAEM",
+  journal="Journal of Computational Physics",
+  volume="227",
+  number="7",
+  pages="3572--3594, doi:10.1016/j.jcp.2007.04.029",
+  year="2008"}
+
+
+@article{Janssen_JGR12,
+  author="Janssen, PAEM",
+  title="{Ocean Wave Effects on the Daily Cycle in SST}",
+  year="2012",
+  journal=JGR,
+  volume="117",
+  pages="C00J32, 24 pp, doi:10/mth"}
+
+@TECHREPORT{Janssen_ECMWF13,
+      AUTHOR = "PAEM Janssen and {\O} Breivik and K Mogensen and F Vitart and 
+                M Balmaseda and JR Bidlot and S Keeley and M Leutbecher and 
+                L Magnusson and F Molteni",
+      TITLE = "{Air-Sea Interaction and Surface Waves}",
+      YEAR = "2013",
+      PAGES = "36",
+      NUMBER = "712",
+      URL = "http://www.ecmwf.int/publications/library/do/references/list/14",
+      TYPE = "ECMWF Technical Memorandum",
+      INSTITUTION = "European Centre for Medium-Range Weather Forecasts"}
 
 @ARTICLE{Jayne_St_Laurent_GRL01,
@@ -2681,4 +2795,12 @@
 }
 
+@ARTICLE{Stokes_TCPS47,
+      AUTHOR = "G~G Stokes",
+      YEAR = "1847",
+      TITLE = "{On the theory of oscillatory waves}",
+      JOURNAL = "Trans Cambridge Philos Soc",
+      VOLUME = "8",
+      PAGES = "441--455"}
+
 @ARTICLE{Talagrand_JAS72,
   author = {O. Talagrand},
@@ -2859,4 +2981,12 @@
 }
 
+@TECHREPORT{wam38r1,
+      AUTHOR = "ECMWF",
+      TITLE = "{IFS Documentation CY38r1, Part VII: ECMWF Wave Model}",
+      YEAR = "2012",
+      PAGES = "77 pp, available at http://ecmwf.int/research/ifsdocs/CY38r1/",
+      TYPE = "{ECMWF Model Documentation}",
+      INSTITUTION = "European Centre for Medium-Range Weather Forecasts"}
+
 @ARTICLE{Warner_al_OM05,
   author = {J. C. Warner and C. R. Sherwood and H. G. Arango and R. P. Signell},
@@ -2956,4 +3086,13 @@
   pages = {593--611}
 }
+
+@MANUAL{wmo98,
+      AUTHOR = "{World Meteorological Organization}",
+      TITLE = "{Guide to wave analysis and forecasting}",
+      YEAR = "1998",
+      ADDRESS = "Geneva, Switzerland",
+      NUMBER = "702",
+      EDITION = "2",
+      ORGANIZATION = "World Meteorological Organization"}
 
 @ARTICLE{Zalesak_JCP79,
Index: branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Chapters/Chap_DYN.tex
===================================================================
--- branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Chapters/Chap_DYN.tex	(revision 4644)
@@ -1,4 +1,4 @@
 % ================================================================
-% Chapter � Ocean Dynamics (DYN)
+% Chapter Ocean Dynamics (DYN)
 % ================================================================
 \chapter{Ocean Dynamics (DYN)}
@@ -795,5 +795,5 @@
 %>   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >   >
 \begin{figure}[!t]    \begin{center}
-\includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_time_split.pdf}
+\includegraphics[width=0.7\textwidth]{./TexFiles/Figures/Fig_DYN_dynspg_ts.pdf}
 \caption{  \label{Fig_DYN_dynspg_ts}
 Schematic of the split-explicit time stepping scheme for the external 
@@ -1293,2 +1293,91 @@
 
 % ================================================================
+% Coriolis-Stokes force
+% ================================================================
+\section  [Coriolis-Stokes Force (\textit{dynstcor})]
+                {Coriolis-Stokes Force (\mdl{dynstcor})}
+\label{DYN_stcor}
+Waves set up a Lagrangian drift in the down-wave direction known
+as the Stokes drift \citep{Stokes_TCPS47}. Although its drift speed
+$\mathbf{v}_\mathrm{s}$ decays rapidly with depth, it can be substantial
+near the surface ($v_\mathrm{s} {\sim}0.7\, \mathrm{m/s}$). In combination
+with the earth's rotation it adds an additional veering to the upper-ocean
+currents known as the Coriolis-Stokes force \citep{Hasselmann_GAFD70},
+\begin{equation}
+   \frac{D\mathbf{u}}{Dt} = -\frac{1}{\rho} \nabla p 
+  + (\mathbf{u} + \mathbf{v}_\mathrm{s}) \times f\hat{\mathbf{z}}
+  + \frac{1}{\rho} \frac{\partial \tau}{\partial z}.
+  \label{Eq_dynstcor_stcor}
+\end{equation}
+It requires integration of the full two-dimensional spectrum to get the
+Stokes profile \citep{Janssen_AH04,Janssen_JGR12},
+\begin{equation}
+   \mathbf{v}_\mathrm{s}(z) = 4\pi \int_0^{2\pi} \int_0^{\infty} 
+                              f \mathbf{k} e^{2kz} F(f,\theta) \, df\, d\theta,
+   \label{Eq_dynstcor_uvfth}
+\end{equation}
+This is computationally demanding and requires access to the full
+two-dimensional wave spectra from a numerical wave model (see e.g. the
+ECMWF WAM implementation, ECWAM, \citealt{wam38r1}), so
+we introduce a parameterized Stokes drift velocity profile
+\citep{Janssen_ECMWF13,Breivik_ECMWF13},
+\begin{equation}
+   \mathbf{v}_\mathrm{e} = \mathbf{v}_0 
+   \frac{e^{2k_\mathrm{e}z}}{1-8k_\mathrm{e}z}.
+   \label{Eq_dynstcor_uve1}
+\end{equation}
+The surface velocity vector $\mathbf{v}_0$ is computed by ECWAM and is
+available both in ERA-Interim \citep{Dee_QJRMS11} and operationally.
+
+The transport under such a profile involves the exponential integral $E_1$ and 
+can be solved analytically \citep{Breivik_ECMWF13} to yield
+\begin{equation}
+   {T}_\mathrm{s} = \frac{{v}_0 e^{1/4} E_1(1/4)}{8 k_\mathrm{e}}.
+   \label{Eq_dynstcor_UVe} 
+\end{equation}
+This imposes the following constraint on the wavenumber,
+\begin{equation}
+   k_\mathrm{e} = \frac{{v}_0 e^{1/4} E_1(1/4)}{8
+   {T}_\mathrm{s}}.
+   \label{Eq_dynstcor_ke} 
+\end{equation}
+Here $E_1(1/4) \approx 1.34$, thus
+\begin{equation}
+   k_\mathrm{e} \approx \frac{{v}_0}{5.97{T}_\mathrm{s}}.
+   \label{Eq_dynstcor_keapprox} 
+\end{equation}
+The $n$-th order spectral moment is defined as
+\begin{equation}
+   m_{n} = \int_0^{2\pi} \int_0^{\infty} 
+           f^{n} F(f,\theta) \, df\, d\theta.
+   \label{Eq_dynstcor_moment}
+\end{equation}
+The mean frequency is defined as $\overline{f} = m_1/m_0$
+\citep{wmo98,Holthuijsen07} and the significant wave height $H_{m_0} =
+4\sqrt{m_0}$.  We can derive the first moment from the integrated parameters
+of a wave model or from wave observations and find an estimate for the
+Stokes transport,
+\begin{equation}
+  \mathbf{T}_\mathrm{s} \approx \frac{2\pi}{16} \overline{f} H_{m_0}^2
+  \hat{\mathbf{k}}_\mathrm{s}.
+  \label{Eq_dynstcor_UVHsf}
+\end{equation}
+Here $\hat{\mathbf{k}}_\mathrm{s} = (\sin \theta_\mathrm{s},
+\cos \theta_\mathrm{s})$ is the unit vector in the
+direction $\theta_\mathrm{s}$ of the Stokes transport.  From
+Eqs~(\ref{Eq_dynstcor_keapprox})-(\ref{Eq_dynstcor_UVHsf}) it is clear that
+in order to compute the Stokes drift velocity profile at the desired vertical
+levels we need $H_\mathrm{s}$, $\overline{f}$ and $\mathbf{v}_0$.
+
+The Coriolis-Stokes effect is enabled when \np{ln\_stcor} = true (default =
+false). All wave-related switches are found in \ngn{namsbc}.
+The surface Stokes drift velocity vectors (east and north components) are
+archived in ERA-Interim as GRIB parameters 215 and 216 respectively (table
+140).
+%\smallskip
+%%----------------------------------------------namsbc----------------------------------------------------
+%\namdisplay{namsbc}
+%%--------------------------------------------------------------------------------------------------------
+%\smallskip
+%
+% ================================================================
Index: branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Chapters/Chap_SBC.tex
===================================================================
--- branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Chapters/Chap_SBC.tex	(revision 4644)
@@ -15,11 +15,11 @@
 The ocean needs six fields as surface boundary condition:
 \begin{itemize}
-	\item the two components of the surface ocean stress $\left( {\tau _u \;,\;\tau _v} \right)$
-	\item the incoming solar and non solar heat fluxes $\left( {Q_{ns} \;,\;Q_{sr} } \right)$
-	\item the surface freshwater budget $\left( {\textit{emp},\;\textit{emp}_S } \right)$
+    \item the two components of the surface ocean stress $\left( {\tau _u \;,\;\tau _v} \right)$
+    \item the incoming solar and non solar heat fluxes $\left( {Q_{ns} \;,\;Q_{sr} } \right)$
+    \item the surface freshwater budget $\left( {\textit{emp},\;\textit{emp}_S } \right)$
 \end{itemize}
 plus an optional field:
 \begin{itemize}
-	\item the atmospheric pressure at the ocean surface $\left( p_a \right)$
+    \item the atmospheric pressure at the ocean surface $\left( p_a \right)$
 \end{itemize}
 
@@ -75,5 +75,5 @@
 \begin{equation} \label{Eq_sbc_dynzdf}
 \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1}
-	 = \frac{1}{\rho _o} \binom{\tau _u}{\tau _v }
+     = \frac{1}{\rho _o} \binom{\tau _u}{\tau _v }
 \end{equation}
 where $(\tau _u ,\;\tau _v )=(utau,vtau)$ are the two components of the wind 
@@ -348,6 +348,6 @@
 horizontal and vertical dimensions of the associated variable and should 
 be equal to 1 over land and 0 elsewhere.
-The procedure can be recursively applied setting nn_lsm > 1 in namsbc namelist. 
-Note that nn_lsm=0 forces the code to not apply the procedure even if a file for land/sea mask is supplied.
+The procedure can be recursively applied setting nn\_lsm > 1 in namsbc namelist. 
+Note that nn\_lsm=0 forces the code to not apply the procedure even if a file for land/sea mask is supplied.
 
 \subsubsection{Bilinear Interpolation}
@@ -565,16 +565,16 @@
 
 The atmospheric fields used depend on the bulk formulae used. Three bulk formulations 
-are available : the CORE, the CLIO and the MFS bulk formulea. The choice is made by setting to true
+are available : the CORE, the CLIO and the MFS bulk formulae. The choice is made by setting to true
 one of the following namelist variable : \np{ln\_core} ; \np{ln\_clio} or  \np{ln\_mfs}.
 
 Note : in forced mode, when a sea-ice model is used, a bulk formulation (CLIO or CORE) have to be used. 
-Therefore the two bulk (CLIO and CORE) formulea include the computation of the fluxes over both 
+Therefore the two bulk (CLIO and CORE) formulae include the computation of the fluxes over both 
 an ocean and an ice surface. 
 
 % -------------------------------------------------------------------------------------------------------------
-%        CORE Bulk formulea
-% -------------------------------------------------------------------------------------------------------------
-\subsection    [CORE Bulk formulea (\np{ln\_core}=true)]
-		      {CORE Bulk formulea (\np{ln\_core}=true, \mdl{sbcblk\_core})}
+%        CORE Bulk formulae
+% -------------------------------------------------------------------------------------------------------------
+\subsection    [CORE Bulk formulae (\np{ln\_core}=true)]
+		      {CORE Bulk formulae (\np{ln\_core}=true, \mdl{sbcblk\_core})}
 \label{SBC_blk_core}
 %------------------------------------------namsbc_core----------------------------------------------------
@@ -591,5 +591,5 @@
 
 Note that substituting ERA40 to NCEP reanalysis fields 
-does not require changes in the bulk formulea themself. 
+does not require changes in the bulk formulae themself. 
 This is the so-called DRAKKAR Forcing Set (DFS) \citep{Brodeau_al_OM09}. 
 
@@ -621,9 +621,87 @@
 or larger than the one of the input atmospheric fields.
 
-% -------------------------------------------------------------------------------------------------------------
-%        CLIO Bulk formulea
-% -------------------------------------------------------------------------------------------------------------
-\subsection    [CLIO Bulk formulea (\np{ln\_clio}=true)]
-		      {CLIO Bulk formulea (\np{ln\_clio}=true, \mdl{sbcblk\_clio})}
+\subsubsection    [The ECMWF parametric drag law (\np{ln\_cdec}=true)]
+                  {The ECMWF parametric drag law (\np{ln\_cdec}=true)}
+As an alternative to the \citet{Large_Yeager_Rep04} drag law the
+parameterization used operationally by ECMWF \citep{Janssen_Rep08,Edson_JPO13} is
+included,
+\begin{equation}
+   C_\mathrm{D}(z=10 \, \mathrm{m}) = \left(a + bU_{10}^{p_1}\right)/U_{10}^{p_2}.
+   \label{Eq_blk_core_cdec}
+\end{equation}
+The coefficients are $a = 1.03 \times 10^{-3}$, $b = 0.04\times 10^{-3}$,
+$p_1 = 1.48$ and $p_2 = 0.21$. 
+
+\subsubsection    [Wave-modified air-side stress (\np{ln\_cdgw}=true)]
+                  {Wave-modified air-side stress (\np{ln\_cdgw}=true)}
+The atmospheric momentum flux to the ocean is denoted $\tau_\mathrm{a}$. It is
+customary to define an air-side friction velocity as $u_*^2 = 
+\tau_\mathrm{a}/\rho_\mathrm{a}$.
+\citet{Charnock_QJRMS55} was the first to relate the roughness of the sea
+surface to the friction velocity,
+\begin{equation}
+   z_0 = \alpha_\mathrm{CH} \frac{u_{*}^2}{g},
+\end{equation}
+where $\alpha_\mathrm{CH}$ is known as the Charnock constant.
+\citet{Janssen_JPO89} showed that $\alpha$ is not constant but varies with
+the sea state,
+\begin{equation}
+   \alpha_\mathrm{CH} =
+   \frac{\hat{\alpha}_\mathrm{CH}}
+        {\sqrt{1-\tau_\mathrm{w}/\tau_\mathrm{a}}},
+\end{equation}
+where $\hat{\alpha}_\mathrm{CH} = 0.01$ and the wave-induced stress,
+$\tau_\mathrm{in}$, is related to the wind input as
+\begin{equation}
+   \boldsymbol{\tau}_\mathrm{in} = \rho_\mathrm{w}g \int_0^{2\pi} 
+      \int_0^{\infty} \frac{\mathbf{k}}{\omega} S_\mathrm{in} \, 
+      d\omega \, d\theta.
+   \label{Eq_blk_core_tauin}
+\end{equation}
+The wave-modified drag coefficient is then
+\begin{equation}
+   C_\mathrm{D} = \frac{\kappa^2}{\log^2(10/z_0)}.
+\end{equation}
+This parameter is stored as CDWW (GRIB parameter 233, table 140) in ERA-Interim and operationally by ECMWF.
+Note that it is used in conjunction with the 10-m \emph{neutral} wind speed,
+$U_\mathrm{10N}$, also archived. The wind direction is taken from the 10-m
+wind vector as before, and only the wind \emph{speed} is changed.  Note also
+that where there is a discrepancy between the ice cover of the wave model
+and NEMO, a drag parametric drag law should used.  Where the wave model
+has ice (as $C_\mathrm{D} = 0$ under ice), a drag law such as the one put
+forward by \citet{Large_Yeager_Rep04} or the one used operationally by ECMWF,
+see Eq~(\ref{Eq_blk_core_cdec}), must be used to pad the fields.
+
+\subsubsection    [Wave-modified water-side stress (\np{ln\_tauoc}=true)]
+                  {Wave-modified water-side stress (\np{ln\_tauoc}=true)}
+As waves break they feed momentum
+into the currents. If wind input and dissipation in the wave field were in
+equilibrium, the air-side stress would be equal to the total water-side
+stress. However, most of the time waves are not in equilibrium
+\citep{Janssen_JGR12,Janssen_ECMWF13}, giving
+differences in air-side and water-side stress of the order of 5-10\%.
+The water-side stress is the total
+atmospheric stress minus the momentum absorbed by the wave field (positive)
+minus the momentum injected from breaking waves to the ocean (negative), 
+$\boldsymbol{\tau}_\mathrm{oc} = \boldsymbol{\tau}_\mathrm{a} -
+\boldsymbol{\tau}_\mathrm{in} - \boldsymbol{\tau}_\mathrm{ds}$. This can be
+written \citep{wam38r1}
+\begin{equation}
+   \boldsymbol{\tau}_\mathrm{oc} = \boldsymbol{\tau}_\mathrm{a} -
+   \rho_\mathrm{w}g \int_0^{2\pi} \int_0^{\infty} 
+   \frac{\mathbf{k}}{\omega}(S_\mathrm{in} + S_\mathrm{ds})\, d\omega d\theta.
+   \label{eq:tauoc}
+\end{equation}
+This parameter is known as TAUOC (GRIB parameter 214, table 140) is stored in
+normalized form, $\tilde{\tau} = \tau_\mathrm{oc}/\tau_\mathrm{a}$, in
+ERA-Interim and operationally by ECMWF.  It is controlled by the namelist
+parameter \np{ln\_tauoc} in namelist \ngn{namsbc}.
+
+
+% -------------------------------------------------------------------------------------------------------------
+%        CLIO Bulk formulae
+% -------------------------------------------------------------------------------------------------------------
+\subsection    [CLIO Bulk formulae (\np{ln\_clio}=true)]
+		      {CLIO Bulk formulae (\np{ln\_clio}=true, \mdl{sbcblk\_clio})}
 \label{SBC_blk_clio}
 %------------------------------------------namsbc_clio----------------------------------------------------
@@ -665,6 +743,6 @@
 %        MFS Bulk formulae
 % -------------------------------------------------------------------------------------------------------------
-\subsection    [MFS Bulk formulea (\np{ln\_mfs}=true)]
-		      {MFS Bulk formulea (\np{ln\_mfs}=true, \mdl{sbcblk\_mfs})}
+\subsection    [MFS Bulk formulae (\np{ln\_mfs}=true)]
+		      {MFS Bulk formulae (\np{ln\_mfs}=true, \mdl{sbcblk\_mfs})}
 \label{SBC_blk_mfs}
 %------------------------------------------namsbc_mfs----------------------------------------------------
@@ -1048,5 +1126,5 @@
 of incident SWF. The \cite{Bernie_al_CD07} reconstruction algorithm is available
 in \NEMO by setting \np{ln\_dm2dc}~=~true (a \textit{\ngn{namsbc}} namelist variable) when using 
-CORE bulk formulea (\np{ln\_blk\_core}~=~true) or the flux formulation (\np{ln\_flx}~=~true). 
+CORE bulk formulae (\np{ln\_blk\_core}~=~true) or the flux formulation (\np{ln\_flx}~=~true). 
 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_al_CD07}. The algorithm preserve the daily 
Index: branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Chapters/Chap_ZDF.tex
===================================================================
--- branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Chapters/Chap_ZDF.tex	(revision 4644)
@@ -197,5 +197,5 @@
 instabilities associated with too weak vertical diffusion. They must be 
 specified at least larger than the molecular values, and are set through 
-\np{rn\_avm0} and \np{rn\_avt0} (namzdf namelist, see \S\ref{ZDF_cst}).
+\np{rn\_avm0} and \np{rn\_avt0} (\ngn{namzdf} namelist, see \S\ref{ZDF_cst}).
 
 \subsubsection{Turbulent length scale}
@@ -262,11 +262,12 @@
 \end{equation}
 
-At the ocean surface, a non zero length scale is set through the  \np{rn\_lmin0} namelist 
-parameter. Usually the surface scale is given by $l_o = \kappa \,z_o$ 
-where $\kappa = 0.4$ is von Karman's constant and $z_o$ the roughness 
-parameter of the surface. Assuming $z_o=0.1$~m \citep{Craig_Banner_JPO94} 
-leads to a 0.04~m, the default value of \np{rn\_lsurf}. In the ocean interior 
-a minimum length scale is set to recover the molecular viscosity when $\bar{e}$ 
-reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ).
+At the ocean surface, a non zero length scale is set through the
+\np{rn\_lmin0} namelist parameter. Usually the surface scale is given
+by $l_o = \kappa \,z_o$ where $\kappa = 0.4$ is von Karman's constant
+and $z_o$ the roughness parameter of the surface. Assuming $z_o=0.1$~m
+\citep{Craig_Banner_JPO94} leads to a 0.04~m, the default value of
+\np{rn\_lsurf}. In the ocean interior a minimum length scale is set to
+recover the molecular viscosity when $\bar{e}$ reach its minimum value
+($1\times 10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$).
 
 
@@ -283,20 +284,93 @@
 \bar{e}_o = \frac{1}{2}\,\left(  15.8\,\alpha_{CB} \right)^{2/3} \,\frac{|\tau|}{\rho_o}
 \end{equation}
-where $\alpha_{CB}$ is the \citet{Craig_Banner_JPO94} constant of proportionality 
-which depends on the ''wave age'', ranging from 57 for mature waves to 146 for 
-younger waves \citep{Mellor_Blumberg_JPO04}. 
-The boundary condition on the turbulent length scale follows the Charnock's relation:
+where $\alpha_{CB}$ is the \citet{Craig_Banner_JPO94} constant of
+proportionality which depends on the ''wave age'', ranging from 57 for mature
+waves to 146 for younger waves \citep{Mellor_Blumberg_JPO04}.  The boundary
+condition on the turbulent length scale follows Charnock's relation:
 \begin{equation} \label{ZDF_Lsbc}
 l_o = \kappa \beta \,\frac{|\tau|}{g\,\rho_o}
 \end{equation}
-where $\kappa=0.40$ is the von Karman constant, and $\beta$ is the Charnock's constant.
-\citet{Mellor_Blumberg_JPO04} suggest $\beta = 2.10^{5}$ the value chosen by \citet{Stacey_JPO99}
-citing observation evidence, and $\alpha_{CB} = 100$ the Craig and Banner's value.
-As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$, 
-with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}~=~67.83 corresponds 
-to $\alpha_{CB} = 100$. further setting  \np{ln\_lsurf} to true applies \eqref{ZDF_Lsbc} 
-as surface boundary condition on length scale, with $\beta$ hard coded to the Stacet's value.
-Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) 
-is applied on surface $\bar{e}$ value.
+where $\kappa=0.40$ is the von Karman constant, and $\beta$ is Charnock's
+constant.  \citet{Mellor_Blumberg_JPO04} suggest $\beta = 2\times10^{5}$ the value
+chosen by \citet{Stacey_JPO99} citing observation evidence, and $\alpha_{CB}
+= 100$ the Craig and Banner's value.  As the surface boundary condition
+on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$, with
+$e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}~=~67.83
+corresponds to $\alpha_{CB} = 100$. further setting  \np{ln\_lsurf} to true
+applies \eqref{ZDF_Lsbc} as surface boundary condition on length scale, with
+$\beta$ hard coded to Stacey's value.  Note that a minimal threshold
+of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters) is applied on
+surface $\bar{e}$ value.
+
+\subsubsection{Surface wave breaking flux from a wave model \np{ln\_wavetke}}
+%-----------------------------------------------------------------------%
+The constant of proportionality $\alpha_{CB}$ in Eq~\eqref{ZDF_Esbc} relates
+the water-side friction velocity $w_*$ to the turbulent energy flux as follows,
+\begin{equation}
+   \Phi_\mathrm{oc} = \rho_\mathrm{w} \alpha_\mathrm{CB} w_*^3.
+   \label{Eq_ZDF_alpha}
+\end{equation}
+The default option in NEMO \eqref{ZDF_Esbc} is to assume $\alpha_{CB}
+= 100$ as explained in the previous section.
+However, the energy flux can be computed from the dissipation source term
+of a wave model \citep{Janssen_AH04,Janssen_JGR12,Janssen_ECMWF13},
+\begin{equation}
+   \Phi_\mathrm{oc} = \Phi_\mathrm{in} - \rho_\mathrm{w}g \int_0^{2\pi} 
+                      \int_0^{\infty} (S_\mathrm{in} + S_\mathrm{ds})\, 
+                        d\omega d\theta.
+   \label{Eq_ZDF_phioc}
+\end{equation}
+Assuming high-frequency equilibrium and ignoring the direct turbulent energy
+flux from the atmosphere to the ocean we get
+\begin{equation}
+   \Phi_\mathrm{oc} = -\rho_\mathrm{w}g \int_0^{2\pi} 
+                      \int_0^{\omega_\mathrm{c}} S_\mathrm{ds}\, 
+                        d\omega d\theta = -\rho_\mathrm{a} m u_*^3.
+   \label{Eq_ZDF_m}
+\end{equation}
+Here, $m \approx -\sqrt{\rho_\mathrm{a}/\rho_\mathrm{w}} \alpha_\mathrm{CB}$
+is the energy flux \emph{from} the waves (thus always negative) normalized by
+the air friction velocity $u_*$. It
+is archived as PHIOC (GRIB parameter 212, table 140) in ERA-Interim and also
+operationally by ECMWF.  The namelist parameter \np{ln\_wavetke} controls
+the wave TKE flux.  We assume that the flux has been converted to physical
+units following \eqref{Eq_ZDF_m} before ingested by NEMO.
+
+NEMO computes the upper boundary condition following
+\citet{Mellor_Blumberg_JPO04}, see \eqref{ZDF_Esbc}.  Since $e$ varies
+rapidly with depth, we want to weight the surface value $\overline{e}_o$
+by the thickness of the uppermost level to attain a value representative
+for the turbulence level of the uppermost level,
+\begin{equation}
+   \overline{e}_1 = \frac{\overline{e}_o}{L} \int_{-L}^{0} e(z) \,dz.
+   \label{Eq_ZDF_eavg}
+\end{equation}
+Here $L = \Delta z_1/2$ is the depth of the $T$-point of the uppermost level.
+This adjustment is crucial with model configurations with a thick uppermost
+level, e.g. ORCA1L42.  If we assume, as \citet{Mellor_Blumberg_JPO04} do, that in the
+wave-affected layer the roughness length can be set to a constant which we
+choose to be $z_\mathrm{w} = 0.5H_\mathrm{s}$ and that in this near-surface
+region diffusion balances dissipation, the TKE equation attains the simple
+exponential solution \citep{Mellor_Blumberg_JPO04}
+\begin{equation}
+  e(z) = \overline{e}_o \exp(2\lambda z/3).
+   \label{Eq_ZDF_phimb}
+\end{equation}
+Here, the length scale $\lambda^{-1}$ is sea-state dependent, see Eq (10) by
+\citet{Mellor_Blumberg_JPO04},
+\begin{equation}
+   \lambda = [3/(S_q B_1 \kappa^2)]^{1/2}z_\mathrm{w}^{-1} \approx
+   \frac{2.38}{z_\mathrm{w}}.
+\end{equation}
+We have assumed $S_q=0.2$ and $B=16.6$ \citep{Mellor_Yamada_1982}, as
+used in NEMO.  For a wave height of 2.5 m, which is close to the global
+mean, $\lambda^{-1} \approx 0.5\, \mathrm{m}$.  Integrating \eqref{Eq_ZDF_phimb}
+is straightforward, and the average TKE in \eqref{Eq_ZDF_eavg} becomes
+\begin{equation}
+  \overline{e}_1 = \overline{e}_o \frac{3}{2\lambda L} \left[1 -
+  \exp(-2\lambda L/3)\right].
+\end{equation}
+The wave model energy flux is controlled by \np{ln\_wavetke} in namelist
+\ngn{namsbc}.
 
 
@@ -318,5 +392,5 @@
  
 By making an analogy with the characteristic convective velocity scale 
-($e.g.$, \citet{D'Alessio_al_JPO98}), $P_{LC}$ is assumed to be : 
+($e.g.$, \citet{D'Alessio_al_JPO98}), $P_{LC}$ is assumed to be: 
 \begin{equation}
 P_{LC}(z) = \frac{w_{LC}^3(z)}{H_{LC}}
Index: branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Namelist/namsbc
===================================================================
--- branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Namelist/namsbc	(revision 4643)
+++ branches/2014/dev_r4642_WavesWG/DOC/TexFiles/Namelist/namsbc	(revision 4644)
@@ -30,3 +30,6 @@
                            !       is left empty in namelist)  ,
                            !  =1:n number of iterations of land/sea mask application for input fields
+   ln_stcor    = .false.   !   Stokes drift read from wave model
+   ln_wavetke  = .false.   !   Wave parameters from wave model for the TKE BC
+   ln_tauoc    = .false.   !   Wave-modified stress from wave model
 /
