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}}
