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
 /
Index: /branches/2014/dev_r4642_WavesWG/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg
===================================================================
--- /branches/2014/dev_r4642_WavesWG/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg	(revision 4643)
+++ /branches/2014/dev_r4642_WavesWG/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist_cfg	(revision 4644)
@@ -188,2 +188,18 @@
 !-----------------------------------------------------------------------
 /
+!-----------------------------------------------------------------------
+&namsbc_wave  ! External fields from wave model
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_wave_ecmwf   ! ECWMF external wave model drag coefficient
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_wavepar ! namsbc_wavepar  parameters from wave model
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
+&namsbc_waveparlim
+!-----------------------------------------------------------------------
+/
Index: /branches/2014/dev_r4642_WavesWG/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2014/dev_r4642_WavesWG/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4643)
+++ /branches/2014/dev_r4642_WavesWG/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4644)
@@ -240,4 +240,7 @@
    ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
    ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave)
+   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
    cn_iceflx = 'linear'    !  redistribution of solar input into ice categories during coupling ice/atm.
 /
@@ -1164,4 +1167,12 @@
 /
 !-----------------------------------------------------------------------
+&namsbc_wave_ecmwf   ! External fields from wave model
+!-----------------------------------------------------------------------
+!              !  file name  ! frequency (hours) ! variable     ! time interp. !  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask !
+!              !             !  (if <0  months)  !   name       !   (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      !
+   sn_cdg      =  'cdww' ,        6          , 'cdww'       ,     .true.   , .false. , 'monthly'  , ''      , ''       , '' !
+   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+/
+!-----------------------------------------------------------------------
 &namdyn_nept  !   Neptune effect (simplified: lateral and vertical diffusions removed)
 !-----------------------------------------------------------------------
@@ -1178,2 +1189,30 @@
    rn_htrmax         =  200.0   ! max. depth of transition range
 /
+!-----------------------------------------------------------------------
+&namsbc_wavepar ! namsbc_wavepar  parameters from wave model (CCC not identical to namsbc_wave)
+!-----------------------------------------------------------------------
+!              !  file name      ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!                                !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_ust      = 'ust'            ,         6         , 'ust'     ,   .true.    , .false. , 'monthly'  , ''       , ''
+   sn_vst      = 'vst'            ,         6         , 'vst'     ,   .true.    , .false. , 'monthly'  , ''       , ''
+   sn_swh      = 'swh'            ,         6         , 'swh'     ,   .true.    , .false. , 'monthly'  , ''       , ''
+   sn_mwp      = 'mwp'            ,         6         , 'mwp'     ,   .true.    , .false. , 'monthly'  , ''       , ''
+   sn_wspd     = 'wspd'           ,         6         , 'wspd'    ,   .true.    , .false. , 'monthly'  , ''       , ''
+   sn_phioc    = 'physphioc'      ,         6         , 'physphioc',  .true.    , .false. , 'monthly'  , ''       , ''
+   sn_tauoc    = 'tauoc'          ,         6         , 'tauoc'   ,   .true.    , .false. , 'monthly'  , ''       , ''
+   cn_dir_wavepar  = './'      !  root directory for the location of the ECWAM files
+/
+!-----------------------------------------------------------------------
+&namsbc_waveparlim
+!-----------------------------------------------------------------------
+!              !  file name      ! frequency (hours) ! variable  ! time interp. !  clim  ! 'yearly'/ ! weights  ! rotation !
+!                                !  (if <0  months)  !   name    !   (logical)  !  (T/F) ! 'monthly' ! filename ! pairing  !
+   sn_strn   = 'istr'            ,         6         , 'icestrain'     ,   .true.    , .false. , 'monthly'  , ''       , ''
+   cn_dir_waveparlim  = './'      !  root directory for the location of the ECWAM files
+/
+!-----------------------------------------------------------------------
+&namsbc_wave   ! External fields from wave model
+!-----------------------------------------------------------------------
+sn_cdg     = 'cdww' ,        6          , 'cdww' , .true.   , .false. , 'monthly'  ,''         , ''  !
+   cn_dir_cdg  = './'  !  root directory for the location of drag coefficient files
+/
Index: /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90
===================================================================
--- /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 4644)
+++ /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/DYN/dynstcor.F90	(revision 4644)
@@ -0,0 +1,116 @@
+MODULE dynstcor
+   !!======================================================================
+   !!                       ***  MODULE  dynstcor  ***
+   !! Ocean dynamics: Stokes-Coriolis effect.
+   !!
+   !!======================================================================
+   !! History :  0.1  !  2012-10  (Oyvind Breivik)
+   !!----------------------------------------------------------------------
+   !! dyn_stcor        : Add the Stokes-Coriolis forcing to the momentum equation
+   !!----------------------------------------------------------------------
+   USE oce              ! ocean dynamics and tracers
+   USE dom_oce          ! ocean space and time domain
+   USE in_out_manager   ! I/O manager
+   USE lib_mpp          ! distributed memory computing
+   USE prtctl           ! Print control
+   USE phycst
+   USE lbclnk
+   USE wrk_nemo         ! Memory Allocation
+   USE sbcmod           ! Access to ln_stcor (sbc_oce) and wave parameters (sbc_wave)
+   USE sbcwave_ecmwf    ! Wave module
+
+   IMPLICIT NONE
+
+   REAL(wp) :: rn_deptmaxstcor = 150.0_wp ! maximum depth [m] to be affected by Stokes-Coriolis effect
+   !PRIVATE
+
+   !! * Routine accessibility
+   PUBLIC dyn_stcor           ! routine called by step.F90
+
+   !! * Shared module variables
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ustc, vstc  ! Stokes-Coriolis u and v
+
+   !! * Module variables
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+#  include "domzgr_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   OPA 9.0 , implemented by Bedford Institute of Oceanography
+   !!----------------------------------------------------------------------
+
+ CONTAINS
+
+   INTEGER FUNCTION dynstcor_alloc()
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE dynstcor_alloc  ***
+      !!
+      !! History :  1.0  !  2012-10  (Oyvind Breivik)
+      !!----------------------------------------------------------------------
+      ALLOCATE( ustc(jpi,jpj,jpk) , vstc(jpi,jpj,jpk) , STAT=dynstcor_alloc )
+         !
+      IF( dynstcor_alloc /= 0 )   CALL ctl_warn('dynstcor_alloc: array allocate failed.')
+   END FUNCTION dynstcor_alloc
+
+   SUBROUTINE dyn_stcor( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE dyn_stcor ***
+      !!
+      !! ** Purpose:  Add Stokes-Coriolis forcing to horizontal momentum equation.
+      !!
+      !! ** History:  0.1  !   2012-10  oyvind.breivik@ecmwf.int
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt  ! ocean time-step index
+      !!
+      INTEGER  ::  ji, jj, jk        ! dummy loop indices
+      REAL(wp) ::  ztransp, zsp0, zk, zfac
+      REAL(wp) ::  zus, zvs
+      !
+
+      ! Allocation at first time step.
+
+      IF ( kt == nit000 ) THEN
+         IF( dynstcor_alloc() /= 0 ) CALL ctl_stop('dyn_stcor: array allocate failed.')
+      ENDIF
+
+      !!----------------------------------------------------------------------
+      !
+      ! Update velocity tendencies ua, va by adding the Stokes-Coriolis velocities ustc, vstc
+      !
+      DO jk = 1, jpk
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ! Skip deep levels where the Stokes-Coriolis effect is negligible
+               IF (fsdept(ji,jj,jk) <= rn_deptmaxstcor) THEN
+                  ! Stokes transport speed estimated from Hs and Tmean
+                  ztransp = 2.0_wp*rpi*swh_wavepar(ji,jj)**2.0_wp/(16.0_wp*MAX(mwp_wavepar(ji,jj),0.0000001_wp))
+
+                  ! Stokes surface speed
+                  zsp0 = SQRT(ust_wavepar(ji,jj)**2 + vst_wavepar(ji,jj)**2)
+
+                  ! Wavenumber scale
+                  zk = ABS(zsp0)/MAX(ABS(5.97_wp*ztransp),0.0000001_wp)
+
+                  ! Depth attenuation
+                  zfac = EXP(-2.0_wp*zk*fsdept(ji,jj,jk))/(1.0_wp+8.0_wp*zk*fsdept(ji,jj,jk))
+
+                  ! The Stokes-Coriolis forcing
+                  zus =  ff(ji,jj)*vst_wavepar(ji,jj)*zfac
+                  zvs = -ff(ji,jj)*ust_wavepar(ji,jj)*zfac
+
+                  ! Store arrays of tendencies for diagnostics output
+                  ! This may be removed later for efficiency
+                  ustc(ji,jj,jk) = zus
+                  vstc(ji,jj,jk) = zvs
+
+                  ua(ji,jj,jk) = ua(ji,jj,jk) + zus * umask(ji,jj,jk)
+                  va(ji,jj,jk) = va(ji,jj,jk) + zvs * vmask(ji,jj,jk)
+               ENDIF
+            ENDDO
+         ENDDO
+      ENDDO
+      !
+   END SUBROUTINE dyn_stcor
+
+
+END MODULE dynstcor
Index: /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 4643)
+++ /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 4644)
@@ -53,4 +53,8 @@
    LOGICAL , PUBLIC ::   ln_sdw         !: true if 3d stokes drift from wave model
    !
+   LOGICAL , PUBLIC ::   ln_wavetke  = .FALSE.   !: true if wave parameters are read
+   LOGICAL , PUBLIC ::   ln_stcor    = .FALSE.   !: true if Stokes-Coriolis forcing is included 
+   LOGICAL , PUBLIC ::   ln_tauoc    = .FALSE.   !: true if wave-modified water-side stress is used
+
    LOGICAL , PUBLIC ::   ln_icebergs    !: Icebergs
    !
@@ -69,4 +73,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vtau   , vtau_b   !: sea surface j-stress (ocean referential)     [N/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   taum              !: module of sea surface stress (at T-point)    [N/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   cdn_wave          !: wave model drag coefficient [N/m2] 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wspd_wavepar      !: wave model 10-m neutral wind [m/s] CCC
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tauoc_wavepar     !: normalized stress into the ocean from wave model  CCC
    !! wndm is used onmpute surface gases exchanges in ice-free ocean or leads
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wndm              !: wind speed module at T-point (=|U10m-Uoce|)  [m/s]
@@ -131,4 +138,11 @@
       ALLOCATE( rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     &
          &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) )
+         ! Initialize these since it may not done elsewhere in the code.
+                rnf        (:,:) = 0.0_wp
+                sbc_tsc  (:,:,:) = 0.0_wp
+                qsr_hc   (:,:,:) = 0.0_wp
+                rnf_b      (:,:) = 0.0_wp
+                sbc_tsc_b(:,:,:) = 0.0_wp
+                qsr_hc_b (:,:,:) = 0.0_wp
          !
       ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     &
@@ -138,4 +152,6 @@
          &      ssu_m  (jpi,jpj) , sst_m(jpi,jpj) ,                       &
          &      ssv_m  (jpi,jpj) , sss_m  (jpi,jpj), ssh_m(jpi,jpj) , STAT=ierr(4) )
+      ALLOCATE( cdn_wave(jpi,jpj) )
+      ALLOCATE( wspd_wavepar(jpi,jpj) )
          !
 #if defined key_vvl
Index: /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 4643)
+++ /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 4644)
@@ -38,5 +38,5 @@
    USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
    USE prtctl          ! Print control
-   USE sbcwave,ONLY :  cdn_wave !wave module 
+   !USE sbcwave,ONLY :  cdn_wave !wave module  ! moved to sbc_oce 
 #if defined key_lim3 || defined key_cice
    USE sbc_ice         ! Surface boundary condition: ice fields
@@ -64,4 +64,5 @@
    
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf   ! structure of input fields (file informations, fields read)
+   REAL(wp), PUBLIC, DIMENSION(:,:), ALLOCATABLE :: Cd_n10 ! 10m neutral drag coefficient for output
          
    !                                             !!! CORE bulk parameters
@@ -83,4 +84,6 @@
    REAL(wp) ::   rn_zqt      ! z(q,t) : height of humidity and temperature measurements
    REAL(wp) ::   rn_zu       ! z(u)   : height of wind measurements
+   LOGICAL  ::   ln_cdec   = .FALSE.   ! logical flag for using ec neutral wind drag coef 
+   ! References : P. Janssen, 2008, ECMWF Workshop on Ocean-Atmosphere Interactions, 10-12 November p47-60
 
    !! * Substitutions
@@ -137,6 +140,6 @@
       TYPE(FLD_N) ::   sn_tdif                                 !   "                                 "
       NAMELIST/namsbc_core/ cn_dir , ln_2m  , ln_taudif, rn_pfac, rn_efac, rn_vfac,  &
-         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr ,           &
-         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,           &
+         &                  sn_wndi, sn_wndj, sn_humi  , sn_qsr , ln_cdec,           &
+         &                  sn_qlw , sn_tair, sn_prec  , sn_snow,                    &
          &                  sn_tdif, rn_zqt , ln_bulk2z, rn_zu
       !!---------------------------------------------------------------------
@@ -163,4 +166,5 @@
             sn_qsr%ln_tint = .false.
          ENDIF
+         IF( ln_cdec .AND. lwp ) WRITE(numout,*)'Using Hans Hersbach formula for drag.'
          !                                         ! store namelist information in an array
          slf_i(jp_wndi) = sn_wndi   ;   slf_i(jp_wndj) = sn_wndj
@@ -182,4 +186,7 @@
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' )
          !
+         ! Drag coefficent so we can write to disk
+         ALLOCATE(Cd_n10(jpi,jpj))
+
          sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
          !
@@ -239,4 +246,5 @@
       INTEGER  ::   ji, jj               ! dummy loop indices
       REAL(wp) ::   zcoef_qsatw, zztmp   ! local variable
+      REAL(wp) ::   ztheta               ! local variable, wind direction
       REAL(wp), DIMENSION(:,:), POINTER ::   zwnd_i, zwnd_j    ! wind speed components at T-point
       REAL(wp), DIMENSION(:,:), POINTER ::   zqsatw            ! specific humidity at pst
@@ -283,10 +291,29 @@
 !CDIR COLLAPSE
 #endif
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vect. opt.
-            zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pu(ji-1,jj  ) + pu(ji,jj) )  )
-            zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - rn_vfac * 0.5 * ( pv(ji  ,jj-1) + pv(ji,jj) )  )
+      ! Wind vector at T points
+      IF ( ln_cdgw ) THEN
+         ! Use neutral 10-m wind speed from wave model if available
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vect. opt.
+               ! Local wind direction from non-neutral 10-m wind vector (relative to grid and no current)
+               ! since we do not have the directional information from the wave model
+               ztheta = ATAN2(sf(jp_wndi)%fnow(ji,jj,1), sf(jp_wndj)%fnow(ji,jj,1))
+               ! Wind vector magnitude from 10-m neutral wind speed from wave model
+               zwnd_i(ji,jj) = wspd_wavepar(ji,jj) * SIN(ztheta)
+               zwnd_j(ji,jj) = wspd_wavepar(ji,jj) * COS(ztheta)
+               ! Correct for surface current, 0.0 <= rn_vfac <= 1.0
+               zwnd_i(ji,jj) = (zwnd_i(ji,jj) - 0.5*rn_vfac*(pu(ji-1,jj  ) + pu(ji,jj)))
+               zwnd_j(ji,jj) = (zwnd_j(ji,jj) - 0.5*rn_vfac*(pv(ji  ,jj-1) + pv(ji,jj)))
+            END DO
          END DO
-      END DO
+      ELSE
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vect. opt.
+               ! Correct for surface current, 0.0 <= rn_vfac <= 1.0
+               zwnd_i(ji,jj) = (  sf(jp_wndi)%fnow(ji,jj,1) - 0.5*rn_vfac*( pu(ji-1,jj  ) + pu(ji,jj) )  )
+               zwnd_j(ji,jj) = (  sf(jp_wndj)%fnow(ji,jj,1) - 0.5*rn_vfac*( pv(ji  ,jj-1) + pv(ji,jj) )  )
+            END DO
+         END DO
+      ENDIF
       CALL lbc_lnk( zwnd_i(:,:) , 'T', -1. )
       CALL lbc_lnk( zwnd_j(:,:) , 'T', -1. )
@@ -353,12 +380,27 @@
 
       ! ... tau module, i and j component
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj)
-            taum  (ji,jj) = zztmp * wndm  (ji,jj)
-            zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
-            zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
-         END DO
-      END DO
+      ! Use cdn_wave instead of Cd if we have wave parameters
+      IF ( ln_cdgw ) THEN
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !zztmp = rhoa * wspd_wavepar(ji,jj) * Cd_n10(ji,jj)
+               zztmp = rhoa * wndm(ji,jj) * cdn_wave(ji,jj)
+               !taum(ji,jj) = tauoc_wavepar(ji,jj)
+               !taum(ji,jj) = zztmp * wspd_wavepar(ji,jj)
+               taum  (ji,jj) = zztmp * wndm  (ji,jj)
+               zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
+               zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+            ENDDO
+         ENDDO
+      ELSE
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zztmp = rhoa * wndm(ji,jj) * Cd(ji,jj)
+               taum  (ji,jj) = zztmp * wndm  (ji,jj)
+               zwnd_i(ji,jj) = zztmp * zwnd_i(ji,jj)
+               zwnd_j(ji,jj) = zztmp * zwnd_j(ji,jj)
+            ENDDO
+         ENDDO
+      ENDIF
 
       ! ... add the HF tau contribution to the wind stress module?
@@ -379,4 +421,10 @@
       CALL lbc_lnk( utau(:,:), 'U', -1. )
       CALL lbc_lnk( vtau(:,:), 'V', -1. )
+
+      IF (ln_tauoc) THEN 
+         utau(:,:) = utau(:,:)*tauoc_wavepar(:,:)
+         vtau(:,:) = vtau(:,:)*tauoc_wavepar(:,:)
+         taum(:,:) = taum(:,:)*tauoc_wavepar(:,:)
+      ENDIF
 
       !  Turbulent fluxes over ocean
@@ -755,5 +803,4 @@
       REAL(wp), DIMENSION(:,:), POINTER  ::   dT            ! air/sea temperature differeence      [K]
       REAL(wp), DIMENSION(:,:), POINTER  ::   dq            ! air/sea humidity difference          [K]
-      REAL(wp), DIMENSION(:,:), POINTER  ::   Cd_n10        ! 10m neutral drag coefficient
       REAL(wp), DIMENSION(:,:), POINTER  ::   Ce_n10        ! 10m neutral latent coefficient
       REAL(wp), DIMENSION(:,:), POINTER  ::   Ch_n10        ! 10m neutral sensible coefficient
@@ -775,10 +822,10 @@
       !
       CALL wrk_alloc( jpi,jpj, stab )   ! integer
-      CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
+      CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
       CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m )
 
       !! * Start
       !! Air/sea differences
-      dU10 = max(0.5, dU)     ! we don't want to fall under 0.5 m/s
+      dU10 = MAX(0.5, dU)     ! we don't want to fall under 0.5 m/s
       dT = T_a - sst          ! assuming that T_a is allready the potential temp. at zzu
       dq = q_a - q_sat
@@ -788,10 +835,17 @@
       !!
       !! Neutral Drag Coefficient
-      stab    = 0.5 + sign(0.5,dT)    ! stable : stab = 1 ; unstable : stab = 0 
-      IF  ( ln_cdgw ) THEN
-        cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
-        Cd_n10(:,:) =   cdn_wave
+      stab   = 0.5 + SIGN(0.5,dT)    ! stable : stab = 1 ; unstable : stab = 0 
+
+      ! Use drag coefficient from a wave model ...
+      IF (ln_cdgw) THEN
+         cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
+         Cd_n10(:,:) =  cdn_wave(:,:)
       ELSE
-        Cd_n10  = 1.e-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )    !   L & Y eq. (6a)
+      ! ... or choose a drag law
+         IF (ln_cdec) THEN
+            Cd_n10  = 1E-3 * ( 1.03 + 0.04 * dU10**1.48 ) / (dU10**0.21)    ! PJ (6)
+         ELSE
+            Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )              ! L & Y eq. (6a)
+         ENDIF
       ENDIF
       sqrt_Cd_n10 = sqrt(Cd_n10)
@@ -824,5 +878,9 @@
 
            !! Updating the neutral 10m transfer coefficients :
-           Cd_n10  = 1.e-3 * (2.7/U_n10 + 0.142 + U_n10/13.09)              !  L & Y eq. (6a)
+           IF (ln_cdec) THEN
+             Cd_n10  = 1E-3 * ( 1.03 + 0.04 * U_n10**1.48 ) / (U_n10**0.21)    ! PJ (6)
+           ELSE
+             Cd_n10  = 1E-3 * ( 2.7/U_n10 + 0.142 + U_n10/13.09 )              ! L & Y eq. (6a)
+           ENDIF
            sqrt_Cd_n10 = sqrt(Cd_n10)
            Ce_n10  = 1.e-3 * (34.6 * sqrt_Cd_n10)                           !  L & Y eq. (6b)
@@ -847,5 +905,5 @@
       !!
       CALL wrk_dealloc( jpi,jpj, stab )   ! integer
-      CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
+      CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
       CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta, U_n10, xlogt, xct, zpsi_h, zpsi_m )
       !
@@ -891,5 +949,4 @@
       REAL(wp), DIMENSION(:,:), POINTER ::   dT            ! air/sea temperature differeence   [K]
       REAL(wp), DIMENSION(:,:), POINTER ::   dq            ! air/sea humidity difference       [K]
-      REAL(wp), DIMENSION(:,:), POINTER ::   Cd_n10        ! 10m neutral drag coefficient
       REAL(wp), DIMENSION(:,:), POINTER ::   Ce_n10        ! 10m neutral latent coefficient
       REAL(wp), DIMENSION(:,:), POINTER ::   Ch_n10        ! 10m neutral sensible coefficient
@@ -911,5 +968,5 @@
       IF( nn_timing == 1 )  CALL timing_start('TURB_CORE_2Z')
       !
-      CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
+      CALL wrk_alloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
       CALL wrk_alloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 )
       CALL wrk_alloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m )
@@ -923,9 +980,16 @@
       !! Neutral Drag Coefficient :
       stab = 0.5 + sign(0.5,dT)                 ! stab = 1  if dT > 0  -> STABLE
-      IF( ln_cdgw ) THEN
-        cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
-        Cd_n10(:,:) =   cdn_wave
+
+      ! Use drag coefficient from a wave model ...
+      IF (ln_cdgw) THEN
+         cdn_wave = cdn_wave - rsmall*(tmask(:,:,1)-1)
+         Cd_n10(:,:) =  cdn_wave
       ELSE
-        Cd_n10  = 1.e-3*( 2.7/dU10 + 0.142 + dU10/13.09 ) 
+      ! ... or choose a drag law
+         IF (ln_cdec) THEN
+            Cd_n10  = 1E-3 * ( 1.03 + 0.04 * dU10**1.48 ) / (dU10**0.21)    ! PJ (6)
+         ELSE
+            Cd_n10  = 1E-3 * ( 2.7/dU10 + 0.142 + dU10/13.09 )              ! L & Y eq. (6a)
+         ENDIF
       ENDIF
       sqrt_Cd_n10 = sqrt(Cd_n10)
@@ -996,5 +1060,5 @@
       END DO
       !!
-      CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Cd_n10, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
+      CALL wrk_dealloc( jpi,jpj, dU10, dT, dq, Ce_n10, Ch_n10, sqrt_Cd_n10, sqrt_Cd, L )
       CALL wrk_dealloc( jpi,jpj, T_vpot, T_star, q_star, U_star, zeta_u, zeta_t, U_n10 )
       CALL wrk_dealloc( jpi,jpj, xlogt, xct, zpsi_hu, zpsi_ht, zpsi_m )
Index: /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4643)
+++ /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4644)
@@ -50,4 +50,5 @@
    USE timing           ! Timing
    USE sbcwave          ! Wave module
+   USE sbcwave_ecmwf    ! ECMWF wave module
 
    IMPLICIT NONE
@@ -82,7 +83,8 @@
       INTEGER ::   icpt   ! local integer
       !!
-      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx,  ln_blk_clio, ln_blk_core, ln_cpl,   &
-         &             ln_blk_mfs, ln_apr_dyn, nn_ice,  nn_ice_embd, ln_dm2dc   , ln_rnf,   &
-         &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave , ln_sdw, nn_lsm, cn_iceflx
+      NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx ,  ln_blk_clio, ln_blk_core, ln_cpl,   &
+         &             ln_blk_mfs, ln_apr_dyn, nn_ice ,  nn_ice_embd, ln_dm2dc   , ln_rnf,   &
+         &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave    , ln_sdw     , nn_lsm,   &
+         &             cn_iceflx , ln_stcor  , ln_tauoc, ln_wavetke 
       INTEGER  ::   ios
       !!----------------------------------------------------------------------
@@ -135,4 +137,8 @@
          WRITE(numout,*) '              closed sea (=0/1) (set in namdom)          nn_closea   = ', nn_closea
          WRITE(numout,*) '              n. of iterations if land-sea-mask applied  nn_lsm      = ', nn_lsm
+         WRITE(numout,*) '              Drag coefficient from wave model           ln_cdgw     = ', ln_cdgw  
+         WRITE(numout,*) '              Stokes-Coriolis forcing from wave model    ln_stcor    = ', ln_stcor  
+         WRITE(numout,*) '              TKE wave breaking BC from wave model       ln_wavetke  = ', ln_wavetke
+         WRITE(numout,*) '              Wave-modified stress from wave model       ln_tauoc    = ', ln_tauoc
       ENDIF
 
@@ -230,5 +236,5 @@
       ELSE
       IF ( ln_cdgw .OR. ln_sdw  )                                         & 
-         &   CALL ctl_stop('Not Activated Wave Module (ln_wave=F) but     &
+         &   CALL ctl_warn('Not Activated Wave Module (ln_wave=F) but     &
          & asked coupling with drag coefficient (ln_cdgw =T) or Stokes drift (ln_sdw=T) ')
       ENDIF 
@@ -314,4 +320,11 @@
       CALL sbc_ssm( kt )                                 ! ocean sea surface variables (sst_m, sss_m, ssu_m, ssv_m)
       !                                                  ! averaged over nf_sbc time-step
+
+      IF (ln_cdgw) CALL sbc_wave_ecmwf( kt )
+      ! Read wave parameters if Stokes-Coriolis forcing OR wave parameters for TKE are needed
+      ! Also read it if drag coefficient to ensure that we use the neutral 10m from WAM.
+      IF (ln_stcor .OR. ln_wavetke .OR. ln_tauoc .OR. ln_cdgw) THEN
+         CALL sbc_wavepar( kt )
+      ENDIF
 
       IF (ln_wave) CALL sbc_wave( kt )
Index: /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90
===================================================================
--- /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 4643)
+++ /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave.F90	(revision 4644)
@@ -31,5 +31,4 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_cd	  ! structure of input fields (file informations, fields read) Drag Coefficient
    TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_sd	  ! structure of input fields (file informations, fields read) Stokes Drift
-   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: cdn_wave 
    REAL(wp),ALLOCATABLE,DIMENSION (:,:)              :: usd2d,vsd2d,uwavenum,vwavenum 
    REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:,:)     :: usd3d,vsd3d,wsd3d 
@@ -61,4 +60,5 @@
       USE divcur
       USE wrk_nemo
+      USE sbc_oce, ONLY : cdn_wave
 #if defined key_bdy
       USE bdy_oce, ONLY : bdytmask
@@ -100,6 +100,4 @@
             IF( sn_cdg%ln_tint )   ALLOCATE( sf_cd(1)%fdta(jpi,jpj,1,2) )
             CALL fld_fill( sf_cd, (/ sn_cdg /), cn_dir, 'sbc_wave', 'Wave module ', 'namsbc_wave' )
-            ALLOCATE( cdn_wave(jpi,jpj) )
-            cdn_wave(:,:) = 0.0
         ENDIF
          IF ( ln_sdw ) THEN
Index: /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave_ecmwf.F90
===================================================================
--- /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave_ecmwf.F90	(revision 4644)
+++ /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/SBC/sbcwave_ecmwf.F90	(revision 4644)
@@ -0,0 +1,217 @@
+MODULE sbcwave_ecmwf
+   !!======================================================================
+   !!                       ***  MODULE  sbcwave  ***
+   !! Wave module 
+   !!======================================================================
+   !! History :  3.3.1  !   2011-09  (Adani M)  Original code
+   !!----------------------------------------------------------------------
+   USE iom             ! I/O manager library
+   USE in_out_manager  ! I/O manager
+   USE lib_mpp         ! distribued memory computing library
+   USE fldread         ! read input fields
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE phycst          ! physical constants
+   USE sbcblk_core, ONLY: Cd_n10
+   
+   !!----------------------------------------------------------------------
+   !!   sbc_wave_ecmwf : read drag coefficient from wave model in netcdf files 
+   !!----------------------------------------------------------------------
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   sbc_wavepar       ! 
+   PUBLIC   sbc_wave_ecmwf    ! routine called in sbc_blk_core or sbc_blk_mfs
+   
+   INTEGER , PARAMETER ::   jpfld_wavepar = 7    ! maximum number of fields to be read 
+   INTEGER , PARAMETER ::   jp_ust = 1           ! index of Stokes velocity (east component) (m/s)  at T-point
+   INTEGER , PARAMETER ::   jp_vst = 2           ! index of Stokes velocity (north component) (m/s) at T-point
+   INTEGER , PARAMETER ::   jp_swh = 3           ! index of significant wave height (m)
+   INTEGER , PARAMETER ::   jp_mwp = 4           ! index of mean wave period (s)
+   INTEGER , PARAMETER ::   jp_phioc = 5         ! index of normalized energy flux into the ocean (non-dim)
+   INTEGER , PARAMETER ::   jp_tauoc = 6         ! index of normalized wave stress into the ocean (non-dim)
+   INTEGER , PARAMETER ::   jp_wspd = 7          ! index of 10 m neutral wind speed (m/s)
+
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wave    ! structure of input fields (file informations, fields read)
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:)  :: sf_wavepar ! structure of input fields (file informations, fields read)
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: ust_wavepar, vst_wavepar, swh_wavepar, mwp_wavepar
+   REAL(wp),PUBLIC,ALLOCATABLE,DIMENSION (:,:)       :: phioc_wavepar
+
+   !! * Substitutions
+#  include "domzgr_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2011) 
+   !! $Id: $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE sbc_wave_ecmwf( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_apr  ***
+      !!
+      !! ** Purpose :   read drag coefficient from wave model  in netcdf files.
+      !!
+      !! ** Method  : - Read namelist namsbc_wave
+      !!              - Read Cd_n10 fields in netcdf files 
+      !! ** action  :   
+      !!               
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in  ) ::  kt       ! ocean time step
+      INTEGER                ::  ierror   ! return error code
+      INTEGER                ::  ios      ! io return code
+      CHARACTER(len=100)     ::  cn_dir_cdg                      ! Root directory for location of drag coefficient files
+      TYPE(FLD_N)            ::  sn_cdg                          ! informations about the fields to be read
+      !!---------------------------------------------------------------------
+      NAMELIST/namsbc_wave_ecmwf/  sn_cdg, cn_dir_cdg
+      !!---------------------------------------------------------------------
+
+      !!----------------------------------------------------------------------
+      !
+      !
+      !                                         ! -------------------- !
+      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
+         !                                      ! -------------------- !
+         !                                            !* set file information (default values)
+         ! ... default values (NB: frequency positive => hours, negative => months)
+         !              !   file   ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation !
+         !              !   name   !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
+         sn_cdg = FLD_N( 'cdg_wave'  ,    1.    ,'drag_coeff',  .true.    , .false. ,   'daily'   , ''       , ''      , '' )
+         cn_dir_cdg = './'          ! directory in which the Patm data are 
+         
+
+         REWIND( numnam_ref )              !  Namelist namsbc_wavepar in reference namelist 
+         READ  ( numnam_ref, namsbc_wave_ecmwf, IOSTAT = ios, ERR = 901)
+901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave_ecmwf in reference namelist', lwp )
+
+         REWIND( numnam_cfg )              !  Namelist namsbc_wave_ecmwf in configuration namelist
+         READ  ( numnam_cfg, namsbc_wave_ecmwf, IOSTAT = ios, ERR = 902 )
+902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wave_ecmwf in configuration namelist', lwp )
+         IF(lwm) WRITE ( numond, namsbc_wave_ecmwf )
+         !
+
+         ALLOCATE( sf_wave(1), STAT=ierror )           !* allocate and fill sf_wave with sn_cdg
+         IF( ierror > 0 )   CALL ctl_stop( 'STOP', 'sbc_wave_ecmwf: unable to allocate sf_wave structure' )
+         !
+         CALL fld_fill( sf_wave, (/ sn_cdg /), cn_dir_cdg, 'sbc_wave_ecmwf', 'Wave module ', 'namsbc_wave_ecmwf' )
+                                ALLOCATE( sf_wave(1)%fnow(jpi,jpj,1)   )
+         IF( sn_cdg%ln_tint )   ALLOCATE( sf_wave(1)%fdta(jpi,jpj,1,2) )
+         ! Allocation done by sbc_oce 
+         cdn_wave(:,:) = 0.0
+      ENDIF
+         !
+         !
+      CALL fld_read( kt, nn_fsbc, sf_wave )      !* read drag coefficient from external forcing
+      cdn_wave(:,:) = sf_wave(1)%fnow(:,:,1)
+
+   END SUBROUTINE sbc_wave_ecmwf
+
+
+   SUBROUTINE sbc_wavepar( kt )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE sbc_wavepar  ***
+      !!
+      !! ** Purpose :   Provide at each time step wave model parameters, including the
+      !!                Stokes drift (east and north components), significant wave height and the
+      !!                mean wave period as well as the normalized stress and energy flux into the 
+      !!                ocean for TKE.
+      !!
+      !!
+      !! ** Method  : - Read namelist namsbc_wavepar
+      !!              - Read fields in NetCDF files 
+      !! ** action  :   
+      !!               
+      !!---------------------------------------------------------------------
+      INTEGER, INTENT( in  ) ::  kt  ! ocean time step
+      !!
+      INTEGER  ::  ierror ! return error code
+      INTEGER  ::   ji      ! dummy loop index
+      INTEGER  ::   jfld    ! dummy loop arguments
+      INTEGER  ::   ios     ! io return code
+      !!
+      CHARACTER(len=100)     ::  cn_dir_wavepar ! Root directory for location of ECWAM wave parameter fields
+      TYPE(FLD_N), DIMENSION(jpfld_wavepar) :: slf_i ! array of namelist informations on the fields to read
+      TYPE(FLD_N)            ::  sn_ust, sn_vst, sn_swh, sn_mwp ! information about the fields to be read
+      TYPE(FLD_N)            ::  sn_phioc, sn_tauoc
+      TYPE(FLD_N)            ::  sn_wspd
+      !TYPE(FLD_N)            ::  sn_cdww
+      !!---------------------------------------------------------------------
+      !NAMELIST/namsbc_wavepar/  sn_ust, sn_vst, sn_swh, sn_mwp, sn_phioc, sn_tauoc, cn_dir_wavepar
+      NAMELIST/namsbc_wavepar/  sn_ust, sn_vst, sn_swh, sn_mwp, sn_wspd, sn_phioc, sn_tauoc, cn_dir_wavepar
+      !!---------------------------------------------------------------------
+
+      !!----------------------------------------------------------------------
+      !
+      !
+      !                                         ! -------------------- !
+      IF( kt == nit000 ) THEN                   ! First call kt=nit000 !
+         !                                      ! -------------------- !
+         ! set file information (default values)
+         ! ... default values (NB: frequency positive => hours, negative => months)
+         !            !   file   ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation !
+         !            !   name   !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
+         sn_ust = FLD_N( 'ust'   ,     6     ,  'ust'     ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         sn_vst = FLD_N( 'vst'   ,     6     ,  'vst'     ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         sn_swh = FLD_N( 'swh'   ,     6     ,  'swh'     ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         sn_mwp = FLD_N( 'mwp'   ,     6     ,  'mwp'     ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         !sn_cdww = FLD_N( 'cdww'  ,    6     ,  'cdww'    ,  .true.    , .false. ,   'monthly' , ''       , ''      ,  '' )
+         sn_wspd = FLD_N( 'wspd'   ,   6     ,  'wspd'    ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         sn_phioc = FLD_N( 'phioc'   , 6     ,  'phioc'   ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         sn_tauoc = FLD_N( 'tauoc'   , 6     ,  'tauoc'   ,  .true.    , .false. ,   'monthly' , ''       , ''       ,  '' )
+         cn_dir_wavepar = './'          ! directory in which the wave data are found
+
+         REWIND( numnam_ref )              !  Namelist namsbc_wavepar in reference namelist 
+         READ  ( numnam_ref, namsbc_wavepar, IOSTAT = ios, ERR = 901)
+901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wavepar in reference namelist', lwp )
+
+         REWIND( numnam_cfg )              !  Namelist namsbc_wavepar in configuration namelist
+         READ  ( numnam_cfg, namsbc_wavepar, IOSTAT = ios, ERR = 902 )
+902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namsbc_wavepar in configuration namelist', lwp )
+         IF(lwm) WRITE ( numond, namsbc_wavepar )
+         !
+         slf_i(jp_ust) = sn_ust 
+         slf_i(jp_vst) = sn_vst
+         slf_i(jp_swh) = sn_swh
+         slf_i(jp_mwp) = sn_mwp
+         !slf_i(jp_cdww) = sn_cdww
+         slf_i(jp_wspd) = sn_wspd
+         slf_i(jp_phioc) = sn_phioc
+         slf_i(jp_tauoc) = sn_tauoc
+
+         ALLOCATE( sf_wavepar(jpfld_wavepar), STAT=ierror )          !* set sf_wavepar structure
+         IF ( ierror > 0 ) THEN
+            CALL ctl_stop( 'STOP', 'sbc_wavepar: unable to allocate sf_wavepar structure' ) ; RETURN
+         ENDIF
+         !
+         jfld = jpfld_wavepar
+         DO ji = 1, jpfld_wavepar
+            ALLOCATE( sf_wavepar(ji)%fnow(jpi,jpj,1)   )
+            IF ( slf_i(ji)%ln_tint ) ALLOCATE( sf_wavepar(ji)%fdta(jpi,jpj,1,2) )
+         ENDDO
+         ALLOCATE( ust_wavepar(jpi,jpj) )
+         ALLOCATE( vst_wavepar(jpi,jpj) )
+         ALLOCATE( swh_wavepar(jpi,jpj) )
+         ALLOCATE( mwp_wavepar(jpi,jpj) )
+         !ALLOCATE( wspd_wavepar(jpi,jpj) ) ! allocation done by sbc_oce
+         ALLOCATE( phioc_wavepar(jpi,jpj) )
+         ALLOCATE( tauoc_wavepar(jpi,jpj) )
+
+         !
+         CALL fld_fill( sf_wavepar, slf_i, cn_dir_wavepar, 'sbc_wavepar', 'Wave module ', 'namsbc_wavepar' )
+         !
+      ENDIF
+         !
+      CALL fld_read( kt, nn_fsbc, sf_wavepar )      !* read wave parameters from ECWAM NetCDF file
+      ust_wavepar(:,:) = sf_wavepar(jp_ust)%fnow(:,:,1)
+      vst_wavepar(:,:) = sf_wavepar(jp_vst)%fnow(:,:,1)
+      swh_wavepar(:,:) = sf_wavepar(jp_swh)%fnow(:,:,1)
+      mwp_wavepar(:,:) = sf_wavepar(jp_mwp)%fnow(:,:,1)
+      !cdww_wavepar(:,:) = sf_wavepar(jp_cdww)%fnow(:,:,1)
+      wspd_wavepar(:,:) = sf_wavepar(jp_wspd)%fnow(:,:,1)
+      phioc_wavepar(:,:) = sf_wavepar(jp_phioc)%fnow(:,:,1)
+      tauoc_wavepar(:,:) = sf_wavepar(jp_tauoc)%fnow(:,:,1)
+
+   END SUBROUTINE sbc_wavepar
+
+END MODULE sbcwave_ecmwf
Index: /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90
===================================================================
--- /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 4643)
+++ /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90	(revision 4644)
@@ -26,4 +26,5 @@
    !!                 !                                + cleaning of the parameters + bugs correction
    !!            3.3  !  2010-10  (C. Ethe, G. Madec) reorganisation of initialisation phase
+   !!                 !  2012-12  (oyvind.breivik@ecwmf.int) Changed to wave breaking source term
    !!----------------------------------------------------------------------
 #if defined key_zdftke   ||   defined key_esopa
@@ -42,4 +43,6 @@
    USE domvvl         ! domain: variable volume layer
    USE sbc_oce        ! surface boundary condition: ocean
+   USE sbcwave        ! wave parameters
+   USE sbcwave_ecmwf  ! ECMWF wave parameters
    USE zdf_oce        ! vertical physics: ocean variables
    USE zdfmxl         ! vertical physics: mixed layer
@@ -78,4 +81,9 @@
    LOGICAL  ::   ln_lc     ! Langmuir cells (LC) as a source term of TKE or not
    REAL(wp) ::   rn_lc     ! coef to compute vertical velocity of Langmuir cells
+   REAL(wp) ::   rn_alpavg = 100.0_wp      ! Average alpha (normalized, water-side energy flux, only used in
+                                           ! case we use the source term formulation with no wave model information
+                                           ! NB! The option of running with source term and no wave model has
+                                           ! not been properly tested yet but is included for consistency (Oyvind Breivik, 2013-07-16)
+   REAL(wp) ::   rn_deptmaxtkew = 50.0_wp   ! maximum depth [m] to be affected by TKE from breaking waves 
 
    REAL(wp) ::   ri_cri    ! critic Richardson number (deduced from rn_ediff and rn_ediss values)
@@ -85,4 +93,6 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   en             !: now turbulent kinetic energy   [m2/s2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   e0             !: top level turbulent kinetic energy   [m2/s2] from waves
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   dfac           !: wave depth scaling factor   
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   htau           ! depth of tke penetration (nn_htau)
    REAL(wp)        , ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dissl          ! now mixing lenght of dissipation
@@ -109,4 +119,5 @@
       !!                ***  FUNCTION zdf_tke_alloc  ***
       !!----------------------------------------------------------------------
+      INTEGER :: iexalloc
       ALLOCATE(                                                                    &
 #if defined key_c1d
@@ -118,4 +129,12 @@
          &      avmu_k(jpi,jpj,jpk) , avmv_k(jpi,jpj,jpk), STAT= zdf_tke_alloc      )
          !
+      IF ( ln_wavetke ) THEN
+         ALLOCATE( e0(jpi,jpj) , dfac(jpi,jpj), STAT= iexalloc )
+         ! Setting arrays to zero initially.
+         e0(:,:) = 0.0_wp
+         dfac(:,:) = 0.0_wp
+         !
+         zdf_tke_alloc = zdf_tke_alloc + iexalloc
+      ENDIF
       IF( lk_mpp             )   CALL mpp_sum ( zdf_tke_alloc )
       IF( zdf_tke_alloc /= 0 )   CALL ctl_warn('zdf_tke_alloc: failed to allocate arrays')
@@ -219,4 +238,8 @@
       REAL(wp) ::   zus   , zwlc  , zind            !    -         -
       REAL(wp) ::   zzd_up, zzd_lw                  !    -         -
+      REAL(wp) ::   lambda !    -         -
+      REAL(wp) ::   zphio_flux              !    -         -
+      REAL(wp) ::   sbr, z0, fb                     !    -         -
+
 !!bfr      REAL(wp) ::   zebot                           !    -         -
       INTEGER , POINTER, DIMENSION(:,:  ) :: imlc
@@ -239,10 +262,32 @@
       !                     !  Surface boundary condition on tke
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
-      DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
-         END DO
-      END DO
-      
+      ! Wave parameters read?
+      IF ( ln_wavetke ) THEN
+         ! Breaking-wave TKE injection as flux
+         DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! Roughness length proportional to wave height, capped at Hs = 0.1 m
+               z0 = 0.5_wp*MAX(swh_wavepar(ji,jj), 0.1_wp)
+               lambda = 2.38_wp/z0
+               ! The energy flux from the wave model in physical units
+               zphio_flux = MAX(phioc_wavepar(ji,jj), 0.0_wp)
+               ! TKE at surface (Mellor and Blumberg, 2004)
+               e0(ji,jj) = 0.5_wp*( 15.8_wp*zphio_flux/rau0 )**0.67_wp
+               ! Reducing the surface TKE by scaling with T-depth of first vertical level
+               dfac(ji,jj) = MAX( (2.0_wp*lambda*fsdept(ji,jj,1) )/3.0_wp, 0.00001_wp )
+               ! TKE surface boundary condition weighted by the thickness of the first level
+               en(ji,jj,1) = MAX( rn_emin0, e0(ji,jj) * ( 1.0_wp - EXP(-dfac(ji,jj)) ) / dfac(ji,jj) ) * tmask(ji,jj,1)
+            ENDDO
+         ENDDO
+         ! If wave parameters are not read, use rn_ebb, corresponding to alpha = 100.0
+      ELSE
+         DO jj = 2, jpjm1            ! en(1)   = rn_ebb taum / rau0  (min value rn_emin0)
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ! Estimate average energy flux from std value found in rn_ebb (corresponds to alpha=100.0)
+               en(ji,jj,1) = MAX( rn_emin0, zbbrau * taum(ji,jj) ) * tmask(ji,jj,1)
+            ENDDO
+         ENDDO
+      ENDIF
+ 
 !!bfr   - start commented area
       !                     !<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
@@ -443,4 +488,8 @@
       ENDIF
       CALL lbc_lnk( en, 'W', 1. )      ! Lateral boundary conditions (sign unchanged)
+      IF ( ln_wavetke ) THEN
+         CALL lbc_lnk( e0,     'W', 1. ) ! Lateral boundary conditions (sign unchanged)
+         CALL lbc_lnk( dfac,   'W', 1. ) ! Lateral boundary conditions (sign unchanged)
+      ENDIF
       !
       CALL wrk_dealloc( jpi,jpj, imlc )    ! integer
Index: /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4643)
+++ /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 4644)
@@ -286,4 +286,5 @@
                                CALL dyn_ldf( kstp )         ! lateral mixing
         IF( ln_neptsimp )      CALL dyn_nept_cor( kstp )    ! add Neptune velocities (simplified)
+        IF( ln_stcor )         CALL dyn_stcor   ( kstp )    ! Stokes-Coriolis forcing (ln_stcor set in SBC/sbc_oce.F90)
 #if defined key_agrif
         IF(.NOT. Agrif_Root()) CALL Agrif_Sponge_dyn        ! momemtum sponge
Index: /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step_oce.F90
===================================================================
--- /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 4643)
+++ /branches/2014/dev_r4642_WavesWG/NEMOGCM/NEMO/OPA_SRC/step_oce.F90	(revision 4644)
@@ -23,4 +23,5 @@
 
    USE sbcmod           ! surface boundary condition       (sbc     routine)
+   USE sbc_oce , ONLY: ln_stcor      ! seems to be necessary but should be included in sbcmod????
    USE sbcrnf           ! surface boundary condition: runoff variables
    USE sbccpl           ! surface boundary condition: coupled formulation (call send at end of step)
@@ -51,4 +52,5 @@
    USE dynspg           ! surface pressure gradient        (dyn_spg routine)
    USE dynnept          ! simp. form of Neptune effect(dyn_nept_cor routine)
+   USE dynstcor         ! Stokes-Coriolis forcing          (dyn_stcor routine)
 
    USE dynnxt           ! time-stepping                    (dyn_nxt routine)
