Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/latex/NEMO/subfiles/chap_SBC.tex
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 11831)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 12019)
@@ -1,3 +1,5 @@
 \documentclass[../main/NEMO_manual]{subfiles}
+\usepackage{fontspec}
+\usepackage{fontawesome}
 
 \begin{document}
@@ -504,4 +506,7 @@
 \label{sec:SBC_flx}
 
+% Laurent: DO NOT mix up ``bulk formulae'' (the classic equation) and the ``bulk
+% parameterization'' (i.e NCAR, COARE, ECMWF...)
+
 \begin{listing}
   \nlst{namsbc_flx}
@@ -520,5 +525,13 @@
 See \autoref{subsec:SBC_ssr} for its specification.
 
-%% =================================================================================================
+
+
+
+
+
+
+%% =================================================================================================
+\pagebreak
+\newpage
 \section[Bulk formulation (\textit{sbcblk.F90})]{Bulk formulation (\protect\mdl{sbcblk})}
 \label{sec:SBC_blk}
@@ -530,17 +543,164 @@
 \end{listing}
 
-In the bulk formulation, the surface boundary condition fields are computed with bulk formulae using atmospheric fields
-and ocean (and sea-ice) variables averaged over \np{nn_fsbc}{nn\_fsbc} time-step.
-
-The atmospheric fields used depend on the bulk formulae used.
-In forced mode, when a sea-ice model is used, a specific bulk formulation is used.
-Therefore, different bulk formulae are used for the turbulent fluxes computation
-over the ocean and over sea-ice surface.
-For the ocean, four bulk formulations are available thanks to the \href{https://brodeau.github.io/aerobulk/}{Aerobulk} package (\citet{brodeau.barnier.ea_JPO16}):
-the NCAR (formerly named CORE), COARE 3.0, COARE 3.5 and ECMWF bulk formulae.
-The choice is made by setting to true one of the following namelist variable:
- \np{ln_NCAR}{ln\_NCAR}, \np{ln_COARE_3p0}{ln\_COARE\_3p0},  \np{ln_COARE_3p5}{ln\_COARE\_3p5} and  \np{ln_ECMWF}{ln\_ECMWF}.
-For sea-ice, three possibilities can be selected:
-a constant transfer coefficient (1.4e-3; default value), \citet{lupkes.gryanik.ea_JGR12} (\np{ln_Cd_L12}{ln\_Cd\_L12}), and \citet{lupkes.gryanik_JGR15} (\np{ln_Cd_L15}{ln\_Cd\_L15}) parameterizations
+In the bulk formulation, the surface boundary condition fields are computed with
+bulk formulae using prescribed atmospheric fields and prognostic ocean (and
+sea-ice) surface variables averaged over \np{nn_fsbc}{nn\_fsbc} time-step.
+
+% Turbulent air-sea fluxes are computed using the sea surface properties and
+% atmospheric SSVs at height $z$ above the sea surface, with the traditional
+% aerodynamic bulk formulae:
+
+
+%%% Bulk formulae are this:
+\subsection{Bulk formulae}
+%
+In NEMO, when the bulk formulation is selected, surface fluxes are computed by means of the traditional bulk formulae:
+%
+\begin{subequations}\label{eq_bulk}
+  \begin{eqnarray}
+    \mathbf{\tau} &=& \rho~ C_D ~ \mathbf{U}_z  ~ U_B \label{eq_b_t} \\
+    Q_H           &=& \rho~C_H~C_P~\big[ \theta_z - T_s \big] ~ U_B \label{eq_b_qh} \\
+    E             &=& \rho~C_E    ~\big[    q_s   - q_z \big] ~ U_B \label{eq_b_e}  \\
+    Q_L           &=& -L_v \, E  \label{eq_b_qe} \\
+    %
+    Q_{sr}        &=& (1 - a) Q_{sw\downarrow} \\
+    Q_{ir}        &=& \delta (Q_{lw\downarrow} -\sigma T_s^4)
+  \end{eqnarray}
+\end{subequations}
+%lulu
+%
+From which, the the non-solar heat flux is \[ Q_{ns} = Q_L + Q_H + Q_{ir} \]
+%
+   \[ \theta_z \simeq T_z+\gamma z \]
+   \[  q_s \simeq 0.98\,q_{sat}(T_s,p_a ) \]
+
+
+
+where $\mathbf{\tau}$ is the wind stress vector, $Q_H$ the sensible heat flux,
+$E$ the evaporation, $Q_L$ the latent heat flux, and $Q_{ir}$ the net longwave
+flux.
+%
+$Q_{sw\downarrow}$ and $Q_{lw\downarrow}$ are the surface downwelling shortwave
+and longwave radiative fluxes, respectively.
+%
+Note: a positive sign of $\mathbf{\tau}$, $Q_H$, and $Q_L$ means a gain of the
+relevant quantity for the ocean, while a positive $E$ implies a freshwater loss
+for the ocean.
+%
+$\rho$ is the density of air. $C_D$, $C_H$ and $C_E$ are the BTCs for momentum,
+sensible heat, and moisture, respectively.  $C_P$ is the heat capacity of moist
+air, and $L_v$ is the latent heat of vaporization of water.  $\theta_z$, $T_z$
+and $q_z$ are the potential temperature, temperature, and specific humidity of
+air at height $z$, respectively. $\gamma z$ is a temperature correction term
+which accounts for the adiabatic lapse rate and approximates the potential
+temperature at height $z$ \citep{Josey_al_2013}.  $\mathbf{U}_z$ is the wind
+speed vector at height $z$ (possibly referenced to the surface current
+$\mathbf{u_0}$, section \ref{s_res1}.\ref{ss_current}). The bulk scalar wind
+speed, $U_B$, is the scalar wind speed, $|\mathbf{U}_z|$, with the potential
+inclusion of a gustiness contribution (section
+\ref{s_res2}.\ref{ss_calm}).
+$P_0$ is the mean sea-level pressure (SLP).
+$T_s$ is the sea surface temperature. $q_s$ is the saturation specific humidity
+of air at temperature $T_s$ and includes a 2\% reduction to account for the
+presence of salt in seawater \citep{Sverdrup_al_1942,Kraus_Businger_1996}.
+Depending on the bulk parameterization used, $T_s$ can be the temperature at the
+air-sea interface (skin temperature, hereafter SSST) or at a few tens of
+centimeters below the surface (bulk sea surface temperature, hereafter SST).
+The SSST differs from the SST due to the contributions of two effects of
+opposite sign, the \emph{cool skin} and \emph{warm layer} (hereafter CSWL). The
+\emph{cool skin} refers to the cooling of the millimeter-scale uppermost layer
+of the ocean, in which the net upward flux of heat to the atmosphere is
+ineffectively sustained by molecular diffusion. As such, a steep vertical
+gradient of temperature must exist to ensure the heat flux continuity with
+underlying layers in which the same flux is sustained by turbulence.
+The \emph{warm layer} refers to the warming of the upper few meters of the ocean
+under sunny conditions.
+The CSWL effects are most significant under weak wind conditions due to the
+absence of substancial surface vertical mixing (caused by \eg breaking waves).
+The impact of the CSWL on the computed TASFs is discussed in section
+\ref{s_res1}.\ref{ss_skin}.
+
+
+%%%% Second set of equations (rad):
+where $a$ and $\delta$ are the albedo and emissivity of the sea surface,
+respectively.
+Thus, we use the computed $Q_L$ and $Q_H$ and the 3-hourly surface downwelling
+shortwave and longwave radiative fluxes ($Q_{sw\downarrow}$ and
+$Q_{lw\downarrow}$, respectively) from ERA-Interim to correct the daily SST
+every 3 hours. Due to the implicitness of the problem implied by the dependence
+of $Q_{nsol}$ on $T_s$, this correction is done iteratively during the
+computation of the TASFs.
+
+
+\subsection{Bulk parameterizations}
+
+Accuracy of the estimate of surface turbulent fluxes by means of bulk formulae
+strongly relies on that of the bulk transfer coefficients: $C_D$, $C_H$ and
+$C_E$. They are estimated with what we refer to as a \emph{bulk
+parameterization} algorithm.
+
+... also to adjust humidity and temperature of air to the wind reference measurement
+height (generally 10\,m).
+
+Over the open ocean, four bulk parameterization algorithms are available:
+\begin{itemize}
+\item NCAR, formerly known as CORE, \citep{large.yeager_rpt04}
+\item COARE 3.0 \citep{fairall.bradley.ea_JC03}
+\item COARE 3.6 \citep{edson.jampana.ea_JPO13}
+\item ECMWF (IFS documentation, cy41)
+\end{itemize}
+
+~
+
+% In a typical bulk algorithm, the BTCs under neutral stability conditions are
+% defined using \emph{in-situ} flux measurements while their dependence on the
+% stability is accounted through the \emph{Monin-Obukhov Similarity Theory} and
+% the \emph{flux-profile} relationships \citep[\eg{}][]{Paulson_1970}. BTCs are
+% functions of the wind speed and the near-surface stability of the atmospheric
+% surface layer (hereafter ASL), and hence, depend on $U_B$, $T_s$, $T_z$, $q_s$
+% and $q_z$.
+
+
+
+
+\subsection{Cool-skin and warm-layer parameterizations}
+
+As oposed to the NCAR bulk parameterization, more advanced bulk
+parameterizations such as COARE3.x and ECMWF are meant to be used with the skin
+temperature $T_s$ rather than the bulk SST (which, in NEMO is the temperature at
+the first T-point level).
+%
+So that, technically, the cool-skin and warm-layer parameterization must be
+activated (XXX) to use COARE3.x and ECMWF in a consistant way.
+
+
+\subsection{Air humidity}
+
+Air humidity can be provided as three different parameters: specific humidity
+[kg/kg], relative humidity [\%], or dew-point temperature [K] (LINK to namelist
+parameters)...
+
+
+~\\
+
+
+
+
+The atmospheric fields used depend on the bulk formulae used.  In forced mode,
+when a sea-ice model is used, a specific bulk formulation is used.  Therefore,
+different bulk formulae are used for the turbulent fluxes computation over the
+ocean and over sea-ice surface.
+%
+
+
+thanks to the \href{https://brodeau.github.io/aerobulk/}{Aerobulk} package
+(\citet{brodeau.barnier.ea_JPO16}):
+
+The choice is made by setting to true one of the following namelist
+variable: \np{ln_NCAR}{ln\_NCAR}, \np{ln_COARE_3p0}{ln\_COARE\_3p0}, \np{ln_COARE_3p6}{ln\_COARE\_3p6}
+and \np{ln_ECMWF}{ln\_ECMWF}.  For sea-ice, three possibilities can be selected:
+a constant transfer coefficient (1.4e-3; default
+value), \citet{lupkes.gryanik.ea_JGR12} (\np{ln_Cd_L12}{ln\_Cd\_L12}),
+and \citet{lupkes.gryanik_JGR15} (\np{ln_Cd_L15}{ln\_Cd\_L15}) parameterizations
 
 Common options are defined through the \nam{sbc_blk}{sbc\_blk} namelist variables.
@@ -557,11 +717,13 @@
     j-component of the 10m air velocity  & vtau           & $m.s^{-1}$         & T     \\
     \hline
-    10m air temperature                  & tair           & \r{}$K$            & T     \\
+    10m air temperature                  & tair           & $K$               & T     \\
     \hline
-    Specific humidity                    & humi           & \%                 & T     \\
+    Specific humidity                    & humi           & $-$               & T     \\
+    Relative humidity                    & ~              & $\%$              & T     \\
+    Dew-point temperature                & ~              & $K$               & T     \\    
     \hline
-    Incoming long wave radiation         & qlw            & $W.m^{-2}$         & T     \\
+    Downwelling longwave radiation       & qlw            & $W.m^{-2}$         & T     \\
     \hline
-    Incoming short wave radiation        & qsr            & $W.m^{-2}$         & T     \\
+    Downwelling shortwave radiation      & qsr            & $W.m^{-2}$         & T     \\
     \hline
     Total precipitation (liquid + solid) & precip         & $Kg.m^{-2}.s^{-1}$ & T     \\
@@ -595,13 +757,13 @@
 Its range must be between zero and one, and it is recommended to set it to 0 at low-resolution (ORCA2 configuration).
 
-As for the flux formulation, information about the input data required by the model is provided in
+As for the flux parameterization, information about the input data required by the model is provided in
 the namsbc\_blk namelist (see \autoref{subsec:SBC_fldread}).
 
 %% =================================================================================================
-\subsection[Ocean-Atmosphere Bulk formulae (\textit{sbcblk\_algo\_coare.F90, sbcblk\_algo\_coare3p5.F90, sbcblk\_algo\_ecmwf.F90, sbcblk\_algo\_ncar.F90})]{Ocean-Atmosphere Bulk formulae (\mdl{sbcblk\_algo\_coare}, \mdl{sbcblk\_algo\_coare3p5}, \mdl{sbcblk\_algo\_ecmwf}, \mdl{sbcblk\_algo\_ncar})}
+\subsection[Ocean-Atmosphere Bulk formulae (\textit{sbcblk\_algo\_coare.F90, sbcblk\_algo\_coare3p6.F90, sbcblk\_algo\_ecmwf.F90, sbcblk\_algo\_ncar.F90})]{Ocean-Atmosphere Bulk formulae (\mdl{sbcblk\_algo\_coare}, \mdl{sbcblk\_algo\_coare3p6}, \mdl{sbcblk\_algo\_ecmwf}, \mdl{sbcblk\_algo\_ncar})}
 \label{subsec:SBC_blk_ocean}
 
 Four different bulk algorithms are available to compute surface turbulent momentum and heat fluxes over the ocean.
-COARE 3.0, COARE 3.5 and ECMWF schemes mainly differ by their roughness lenghts computation and consequently
+COARE 3.0, COARE 3.6 and ECMWF schemes mainly differ by their roughness lenghts computation and consequently
 their neutral transfer coefficients relationships with neutral wind.
 \begin{itemize}
@@ -615,5 +777,5 @@
   This is the so-called DRAKKAR Forcing Set (DFS) \citep{brodeau.barnier.ea_OM10}.
 \item COARE 3.0 (\np[=.true.]{ln_COARE_3p0}{ln\_COARE\_3p0}): See \citet{fairall.bradley.ea_JC03} for more details
-\item COARE 3.5 (\np[=.true.]{ln_COARE_3p5}{ln\_COARE\_3p5}): See \citet{edson.jampana.ea_JPO13} for more details
+\item COARE 3.6 (\np[=.true.]{ln_COARE_3p6}{ln\_COARE\_3p6}): See \citet{edson.jampana.ea_JPO13} for more details
 \item ECMWF (\np[=.true.]{ln_ECMWF}{ln\_ECMWF}): Based on \href{https://www.ecmwf.int/node/9221}{IFS (Cy31)} implementation and documentation.
   Surface roughness lengths needed for the Obukhov length are computed following \citet{beljaars_QJRMS95}.
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/latex/global/packages.tex
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/latex/global/packages.tex	(revision 11831)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/latex/global/packages.tex	(revision 12019)
@@ -18,5 +18,5 @@
 %% Issue with fontawesome pkg: path to FontAwesome.otf has to be hard-coded
 \defaultfontfeatures{
-    Path = /usr/local/texlive/2019/texmf-dist/fonts/opentype/public/fontawesome/
+      Path = /usr/share/texlive/texmf-dist/tex/latex/fontawesome/opentype/
 }
 \usepackage{academicons, fontawesome, newtxtext}
Index: NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/namelists/namsbc_blk
===================================================================
--- NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/namelists/namsbc_blk	(revision 11831)
+++ NEMO/branches/2019/dev_r11085_ASINTER-05_Brodeau_Advanced_Bulk/doc/namelists/namsbc_blk	(revision 12019)
@@ -1,21 +1,29 @@
 !-----------------------------------------------------------------------
-&namsbc_blk    !   namsbc_blk  generic Bulk formula                     (ln_blk =T)
+&namsbc_blk    !   namsbc_blk  generic Bulk formula          (ln_blk =T)
 !-----------------------------------------------------------------------
    !                    !  bulk algorithm :
-   ln_NCAR     = .false.   ! "NCAR"      algorithm   (Large and Yeager 2008)
+   ln_NCAR      = .false.   ! "NCAR"      algorithm   (Large and Yeager 2008)
    ln_COARE_3p0 = .false.   ! "COARE 3.0" algorithm   (Fairall et al. 2003)
-   ln_COARE_3p5 = .false.   ! "COARE 3.5" algorithm   (Edson et al. 2013)
-   ln_ECMWF    = .false.   ! "ECMWF"     algorithm   (IFS cycle 31)
+   ln_COARE_3p6 = .false.   ! "COARE 3.6" algorithm   (Edson et al. 2013)
+   ln_ECMWF     = .false.   ! "ECMWF"     algorithm   (IFS cycle 31)
       !
-      rn_zqt      = 10.       !  Air temperature & humidity reference height (m)
-      rn_zu       = 10.       !  Wind vector reference height (m)
-      ln_Cd_L12   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2012)
-      ln_Cd_L15   = .false.   !  air-ice drags = F(ice concentration) (Lupkes et al. 2015)
-      ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data
-      rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow)
-      rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.)
-      rn_vfac     = 0.        !  multiplicative factor for ocean & ice velocity used to
-      !                       !  calculate the wind stress (0.=absolute or 1.=relative winds)
-
+      rn_zqt     = 10.      !  Air temperature & humidity reference height (m)
+      rn_zu      = 10.      !  Wind vector reference height (m)
+      ln_Cd_L12  = .false.  !  air-ice drags = F(ice conc.) (Lupkes et al. 2012)
+      ln_Cd_L15  = .false.  !  air-ice drags = F(ice conc.) (Lupkes et al. 2015)
+      ln_taudif  = .false.  !  HF tau contribution: use "mean of stress module 
+      !                     !  - module of the mean stress" data
+      rn_pfac    = 1.       !  multipl. factor for precipitation (total & snow)
+      rn_efac    = 1.       !  multipl. factor for evaporation (0. or 1.)
+      rn_vfac    = 0.       !  multipl. factor for ocean & ice velocity 
+      !                     !  used to calculate the wind stress
+      !                     ! (0. => absolute or 1. => relative winds)
+      ln_skin_cs = .false.  !  use the cool-skin parameterization
+      ln_skin_wl = .false.  !  use the warm-layer parameterization
+      !                     !   ==> only available in ECMWF and COARE algorithms
+      ln_humi_sph = .true.  !  humidity "sn_humi" is specific humidity  [kg/kg]
+      ln_humi_dpt = .false. !  humidity "sn_humi" is dew-point temperature [K]
+      ln_humi_rlh = .false. !  humidity "sn_humi" is relative humidity     [%]
+   !
    cn_dir      = './'      !  root directory for the bulk data location
    !___________!_________________________!___________________!___________!_____________!________!___________!______________________________________!__________!_______________!
