Changeset 6347 for branches/2016/dev_r6325_SIMPLIF_1/DOC

Timestamp:

2016-02-24T08:56:48+01:00 (8 years ago)

Author:

gm

Message:

#1683: SIMPLIF-1 : Phase with the v3.6_Stable (DOC+ZDF+traqsr+lbedo)

Location:

branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles

Files:

• Biblio/Biblio.bib (modified) (9 diffs)
• Chapters/Chap_DIA.tex (modified) (7 diffs)
• Chapters/Chap_DOM.tex (modified) (7 diffs)
• Chapters/Chap_DYN.tex (modified) (2 diffs)
• Chapters/Chap_LDF.tex (modified) (1 diff)
• Chapters/Chap_SBC.tex (modified) (9 diffs)
• Chapters/Chap_STO.tex (modified) (1 diff)
• Chapters/Chap_TRA.tex (modified) (4 diffs)
• Chapters/Chap_ZDF.tex (modified) (9 diffs)
• clean.sh (modified) (1 diff)

branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Biblio/Biblio.bib

@@ -472 +472 @@
 }
 
+@article{bouffard_Boegman_DAO2013,
+   author = {D. Bouffard and L. Boegman},
+   title = {A diapycnal diffusivity model for stratified environmental flows},
+   volume = {61-62},
+   issn = {03770265},
+   url = {http://dx.doi.org/10.1016/j.dynatmoce.2013.02.002},
+   doi = {10.1016/j.dynatmoce.2013.02.002},
+   journal = DAO,
+   year = {2013},
+   pages = {14--34},
+}
+
 @ARTICLE{Bougeault1989,
   author = {P. Bougeault and P. Lacarrere},
@@ -787 +799 @@
   volume = {34},
   pages = {8--13}
+}
+
+@article{de_lavergne_JPO2016_mixing,
+   author = {C. de Lavergne and G. Madec and J. Le Sommer and A. J. G. Nurser and A. C. Naveira Garabato },
+   title = {On Antarctic Bottom Water consumption in the abyssal ocean},
+   issn = {0022-3670},
+   url = {http://dx.doi.org/10.1175/JPO-D-14-0201.1},
+   doi = {10.1175/JPO-D-14-0201.1},
+   abstract = {In studies of ocean mixing, it is generally assumed that small-scale turbulent overturns lose 15-20 \% of their energy in eroding the background stratification. Accumulating evidence that this energy fraction, or mixing efficiency Rf, significantly varies depending on flow properties challenges this assumption, however. Here, we examine the implications of a varying mixing efficiency for ocean energetics and deep water mass transformation. Combining current parameterizations of internal wave-driven mixing with a recent model expressing Rf as a function of a turbulence intensity parameter Reb = εν/νN2, we show that accounting for reduced mixing efficiencies in regions of weak stratification or energetic turbulence (high Reb) strongly limits the ability of breaking internal waves to supply oceanic potential energy and drive abyssal upwelling. Moving from a fixed Rf = 1/6 to a variable efficiency Rf(Reb) causes Antarctic Bottom Water upwelling induced by locally-dissipating internal tides and lee waves to fall from 9 to 4 Sv, and the corresponding potential energy source to plunge from 97 to 44 GW. When adding the contribution of remotely-dissipating internal tides under idealized distributions of energy dissipation, the total rate of Antarctic Bottom Water upwelling is reduced by about a factor of 2, reaching 5-15 Sv compared to 10-33 Sv for a fixed efficiency. Our results suggest that distributed mixing, overflow-related boundary processes and geothermal heating are more effective in consuming abyssal waters than topographically-enhanced mixing by breaking internal waves. Our calculations also point to the importance of accurately constraining Rf(Reb) and including the effect in ocean models.},
+   journal = {Journal of Physical Oceanography},
+   year = {2016},
+   volume = {46},  pages = {635-–661}
+}
+
+@article{de_lavergne_JPO2016_efficiency,
+   author = {C. de Lavergne and G. Madec and J. Le Sommer and A. J. G. Nurser and A. C. Naveira Garabato },
+   title = {The impact of a variable mixing efficiency on the abyssal overturning},
+   issn = {0022-3670},
+   url = {http://dx.doi.org//10.1175/JPO-D-14-0259.1},
+   doi = {10.1175/JPO-D-14-0259.1},
+   abstract = {In studies of ocean mixing, it is generally assumed that small-scale turbulent overturns lose 15-20 \% of their energy in eroding the background stratification. Accumulating evidence that this energy fraction, or mixing efficiency Rf, significantly varies depending on flow properties challenges this assumption, however. Here, we examine the implications of a varying mixing efficiency for ocean energetics and deep water mass transformation. Combining current parameterizations of internal wave-driven mixing with a recent model expressing Rf as a function of a turbulence intensity parameter Reb = εν/νN2, we show that accounting for reduced mixing efficiencies in regions of weak stratification or energetic turbulence (high Reb) strongly limits the ability of breaking internal waves to supply oceanic potential energy and drive abyssal upwelling. Moving from a fixed Rf = 1/6 to a variable efficiency Rf(Reb) causes Antarctic Bottom Water upwelling induced by locally-dissipating internal tides and lee waves to fall from 9 to 4 Sv, and the corresponding potential energy source to plunge from 97 to 44 GW. When adding the contribution of remotely-dissipating internal tides under idealized distributions of energy dissipation, the total rate of Antarctic Bottom Water upwelling is reduced by about a factor of 2, reaching 5-15 Sv compared to 10-33 Sv for a fixed efficiency. Our results suggest that distributed mixing, overflow-related boundary processes and geothermal heating are more effective in consuming abyssal waters than topographically-enhanced mixing by breaking internal waves. Our calculations also point to the importance of accurately constraining Rf(Reb) and including the effect in ocean models.},
+   journal = {Journal of Physical Oceanography},
+   year = {2016},
+   volume = {46},  pages = {663-–681}
 }
 
@@ -1160 +1196 @@
 }
 
+@article{goff_JGR2010,
+   author = {J. A. Goff},
+   title = {Global prediction of abyssal hill root-mean-square heights from small-scale altimetric gravity variability},
+   issn = {2156-2202},
+   url = {http://dx.doi.org/10.1029/2010JB007867},
+   doi = {10.1029/2010JB007867},
+   abstract = {Abyssal hills, which are pervasive landforms on the seafloor of the Earth's oceans, represent a potential tectonic record of the history of mid-ocean ridge spreading. However, the most detailed global maps of the seafloor, derived from the satellite altimetry-based gravity field, cannot be used to deterministically characterize such small-scale ({\textless}10 km) morphology. Nevertheless, the small-scale variability of the gravity field can be related to the statistical properties of abyssal hill morphology using the upward continuation formulation. In this paper, I construct a global prediction of abyssal hill root-mean-square (rms) heights from the small-scale variability of the altimetric gravity field. The abyssal hill-related component of the gravity field is derived by first masking distinct features, such as seamounts, mid-ocean ridges, and continental margins, and then applying a newly designed adaptive directional filter algorithm to remove fracture zone/discontinuity fabric. A noise field is derived empirically by correlating the rms variability of the small-scale gravity field to the altimetric noise field in regions of very low relief, and the noise variance is subtracted from the small-scale gravity variance. Suites of synthetically derived, abyssal hill formed gravity fields are generated as a function of water depth, basement rms heights, and sediment thickness and used to predict abyssal hill seafloor rms heights from corrected small-scale gravity rms height. The resulting global prediction of abyssal hill rms heights is validated qualitatively by comparing against expected variations in abyssal hill morphology and quantitatively by comparing against actual measurements of rms heights. Although there is scatter, the prediction appears unbiased.},
+   volume = {115},
+   number = {B12},
+   journal = {Journal of Geophysical Research: Solid Earth},
+   year = {2010},
+   pages = {B12104},
+}
+
 @ARTICLE{Goosse_al_JGR99,
   author = {H. Goosse and E. Deleersnijder and T. Fichefet and M. England},
@@ -1264 +1314 @@
 
 @ARTICLE{Griffies_Hallberg_MWR00,
-  author = {S.M. Griffies and R.H. Hallberg},
-  title = {Biharmonic friction with a smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models},
+  author = {S.M. Griffies and R.W. Hallberg},
+  title = {Biharmonic friction with a Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models},
   journal = MWR,
   year = {2000},
@@ -1586 +1636 @@
   volume = {12},
   pages = {381--389}
+}
+
+@article{Jackson_Rehmann_JPO2014,
+   author = {P. R. Jackson and C. R. Rehmann},
+   title = {Experiments on differential scalar mixing in turbulence in a sheared, stratified flow},
+   journal = JPO,
+   volume = {44},
+   issn = {0022-3670},
+   url = {http://dx.doi.org/10.1175/JPO-D-14-0027.1},
+   doi = {10.1175/JPO-D-14-0027.1},
+   number = {10},
+   year = {2014},
+   pages = {2661--2680},
 }
 
@@ -2430 +2493 @@
 }
 
+@ARTICLE{Morel_Berthon_LO89,
+  author = {A. Morel and J.-F. Berthon},
+  title = {Surface pigments, algal biomass profiles, and potential production of the euphotic layer: 
+           Relationships reinvestigated in view of remote-sensing applications},
+  journal = {Limnol. Oceanogr.},
+  year = {1989},
+  volume = {34(8)},
+  pages = {1545--1562}
+}
+
 @ARTICLE{Morel_Maritorena_JGR01,
   author = {A. Morel and S. Maritorena},
@@ -2478 +2551 @@
   title = {Estimates of the local rate of vertical diffusion from dissipation measurements},
   journal = JPO,
+  year = {1980},
   volume = {10},
   pages = {83--89}
@@ -2714 +2788 @@
 }
 
+@ARTICLE{Rousset_GMD2015,
+  author = {C. Rousset and M. Vancoppenolle and G. Madec and T. Fichefet and S. Flavoni 
+            and A. Barth\'{e}lemy and R. Benshila and J. Chanut and C. L\'{e}vy and S. Masson and F. Vivier },
+  title  = {The Louvain-La-Neuve sea-ice model LIM3.6: Global and regional capabilities},
+  journal= {Geoscientific Model Development},
+  year = {2015},
+  volume = {8}, pages={2991--3005},
+  doi = {10.5194/gmd-8-2991-2015},
+  url = {http://dx.doi.org/10.5194/gmd-8-2991-2015}
+}
+
 @ARTICLE{Sadourny1975,
   author = {R. Sadourny},
@@ -2794 +2879 @@
   year = {2004},
   pages = {245--263},
+}
+
+@INBOOK{Smagorinsky_93,
+  author = {Smagorinsky, J.},
+  chapter = {Some historical remarks on the use of non-linear viscosities},
+  title = {Large Eddy Simulation of Complex Engineering and Geophysical Flows},
+  pages = {3--36},
+  year = {1993},
+  publisher = {Cambridge University Press, B. Galperin and S. A. Orszag (eds.)},
 }
 

branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_DIA.tex

@@ -110 +110 @@
 even without a parallel-enabled NetCDF4 library, simply by employing only one dedicated I/O server.
 
-\subsection{XIOS: the IO\_SERVER}
+\subsection{XIOS: the I/O server}
 
 \subsubsection{Attached or detached mode?}
@@ -1409 +1409 @@
 
 % -------------------------------------------------------------------------------------------------------------
-%       25 hour mean and hourly Surface, Mid and Bed 
-% -------------------------------------------------------------------------------------------------------------
-\section{25 hour mean output for tidal models }
-
-A module is available to compute a crudely detided M2 signal by obtaining a 25 hour mean.
-The 25 hour mean is available for daily runs by summing up the 25 hourly instantananeous hourly values from
-midnight at the start of the day to midight at the day end.
-This diagnostic is actived with the logical  $ln\_dia25h$
-
-%------------------------------------------nam_dia25h------------------------------------------------------
-\namdisplay{nam_dia25h}
-%----------------------------------------------------------------------------------------------------------
-
-\section{Top Middle and Bed hourly output }
-
-A module is available to output the surface (top), mid water and bed diagnostics of a set of standard variables. 
-This can be a useful diagnostic when hourly or sub-hourly output is required in high resolution tidal outputs.
-The tidal signal is retained but the overall data usage is cut to just three vertical levels. Also the bottom level 
-is calculated for each cell.
-This diagnostic is actived with the logical  $ln\_diatmb$
-
-%------------------------------------------nam_diatmb-----------------------------------------------------
-\namdisplay{nam_diatmb}
-%----------------------------------------------------------------------------------------------------------
-
-% -------------------------------------------------------------------------------------------------------------
 %       Sections transports
 % -------------------------------------------------------------------------------------------------------------
@@ -1440 +1414 @@
 \label{DIA_diag_dct}
 
+%------------------------------------------namdct----------------------------------------------------
+\namdisplay{namdct}
+%-------------------------------------------------------------------------------------------------------------
+
 A module is available to compute the transport of volume, heat and salt through sections. 
 This diagnostic is actived with \key{diadct}.
@@ -1459 +1437 @@
 and the time scales over which they are averaged, as well as the level of output for debugging:
 
-%------------------------------------------namdct----------------------------------------------------
-\namdisplay{namdct}
-%-------------------------------------------------------------------------------------------------------------
-
 \np{nn\_dct}: frequency of instantaneous transports computing
 
@@ -1469 +1443 @@
 \np{nn\_debug}: debugging of the section
 
-\subsubsection{ To create a binary file containing the pathway of each section }
-
-In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run}, the file \texttt{ {list\_sections.ascii\_global}}
+\subsubsection{ Creating a binary file containing the pathway of each section }
+
+In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run}, the file \textit{ {list\_sections.ascii\_global}}
 contains a list of all the sections that are to be computed (this list of sections is based on MERSEA project metrics).
 
@@ -1583 +1557 @@
 \texttt{=/0, =/ 1000.}   &  diagonal   & eastward  & westward  & postive: eastward  \\ \hline                
 \end{tabular}
-
-
-
-% -------------------------------------------------------------------------------------------------------------
-%       Other Diagnostics
-% -------------------------------------------------------------------------------------------------------------
-\section{Other Diagnostics (\key{diahth}, \key{diaar5})}
-\label{DIA_diag_others}
-
-
-Aside from the standard model variables, other diagnostics can be computed 
-on-line. The available ready-to-add diagnostics routines can be found in directory DIA. 
-Among the available diagnostics the following ones are obtained when defining 
-the \key{diahth} CPP key: 
-
-- the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth})
-
-- the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth})
-
-- the depth of the 20\deg C isotherm (\mdl{diahth})
-
-- the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth})
-
-The poleward heat and salt transports, their advective and diffusive component, and 
-the meriodional stream function can be computed on-line in \mdl{diaptr} 
-\np{ln\_diaptr} to true (see the \textit{\ngn{namptr} } namelist below).  
-When \np{ln\_subbas}~=~true, transports and stream function are computed 
-for the Atlantic, Indian, Pacific and Indo-Pacific Oceans (defined north of 30\deg S) 
-as well as for the World Ocean. The sub-basin decomposition requires an input file 
-(\ifile{subbasins}) which contains three 2D mask arrays, the Indo-Pacific mask 
-been deduced from the sum of the Indian and Pacific mask (Fig~\ref{Fig_mask_subasins}). 
-
-%------------------------------------------namptr----------------------------------------------------
-\namdisplay{namptr} 
-%-------------------------------------------------------------------------------------------------------------
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-\begin{figure}[!t]     \begin{center}
-\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_mask_subasins.pdf}
-\caption{   \label{Fig_mask_subasins}
-Decomposition of the World Ocean (here ORCA2) into sub-basin used in to compute
-the heat and salt transports as well as the meridional stream-function: Atlantic basin (red), 
-Pacific basin (green), Indian basin (bleue), Indo-Pacific basin (bleue+green). 
-Note that semi-enclosed seas (Red, Med and Baltic seas) as well as Hudson Bay 
-are removed from the sub-basins. Note also that the Arctic Ocean has been split 
-into Atlantic and Pacific basins along the North fold line.  }
-\end{center}   \end{figure}  
-%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
-
-In addition, a series of diagnostics has been added in the \mdl{diaar5}. 
-They corresponds to outputs that are required for AR5 simulations 
-(see Section \ref{DIA_steric} below for one of them). 
-Activating those outputs requires to define the \key{diaar5} CPP key.
-\\
-\\
-
-\section{Courant numbers}
-Courant numbers provide a theoretical indication of the model's numerical stability. The advective Courant numbers can be calculated according to
-\begin{equation}
-\label{eq:CFL}
-C_u = |u|\frac{\rdt}{e_{1u}}, \quad C_v = |v|\frac{\rdt}{e_{2v}}, \quad C_w = |w|\frac{\rdt}{e_{3w}}
-\end{equation}
-in the zonal, meridional and vertical directions respectively. The vertical component is included although it is not strictly valid as the vertical velocity is calculated from the continuity equation rather than as a prognostic variable. Physically this represents the rate at which information is propogated across a grid cell. Values greater than 1 indicate that information is propagated across more than one grid cell in a single time step.
-
-The variables can be activated by setting the \np{nn\_diacfl} namelist parameter to 1 in the \ngn{namctl} namelist. The diagnostics will be written out to an ascii file named cfl\_diagnostics.ascii. In this file the maximum value of $C_u$, $C_v$, and $C_w$ are printed at each timestep along with the coordinates of where the maximum value occurs. At the end of the model run the maximum value of $C_u$, $C_v$, and $C_w$ for the whole model run is printed along with the coordinates of each. The maximum values from the run are also copied to the ocean.output file. 
 
 
@@ -1802 +1712 @@
 the \key{diaar5} defined to be called.
 
+
+
+% -------------------------------------------------------------------------------------------------------------
+%       Other Diagnostics
+% -------------------------------------------------------------------------------------------------------------
+\section{Other Diagnostics (\key{diahth}, \key{diaar5})}
+\label{DIA_diag_others}
+
+
+Aside from the standard model variables, other diagnostics can be computed on-line. 
+The available ready-to-add diagnostics modules can be found in directory DIA. 
+
+\subsection{Depth of various quantities (\mdl{diahth})}
+
+Among the available diagnostics the following ones are obtained when defining 
+the \key{diahth} CPP key: 
+
+- the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth})
+
+- the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth})
+
+- the depth of the 20\deg C isotherm (\mdl{diahth})
+
+- the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth})
+
+% -----------------------------------------------------------
+%     Poleward heat and salt transports
+% -----------------------------------------------------------
+
+\subsection{Poleward heat and salt transports (\mdl{diaptr})}
+
+%------------------------------------------namptr-----------------------------------------
+\namdisplay{namptr} 
+%-----------------------------------------------------------------------------------------
+
+The poleward heat and salt transports, their advective and diffusive component, and 
+the meriodional stream function can be computed on-line in \mdl{diaptr} 
+\np{ln\_diaptr} to true (see the \textit{\ngn{namptr} } namelist below).  
+When \np{ln\_subbas}~=~true, transports and stream function are computed 
+for the Atlantic, Indian, Pacific and Indo-Pacific Oceans (defined north of 30\deg S) 
+as well as for the World Ocean. The sub-basin decomposition requires an input file 
+(\ifile{subbasins}) which contains three 2D mask arrays, the Indo-Pacific mask 
+been deduced from the sum of the Indian and Pacific mask (Fig~\ref{Fig_mask_subasins}). 
+
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+\begin{figure}[!t]     \begin{center}
+\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_mask_subasins.pdf}
+\caption{   \label{Fig_mask_subasins}
+Decomposition of the World Ocean (here ORCA2) into sub-basin used in to compute
+the heat and salt transports as well as the meridional stream-function: Atlantic basin (red), 
+Pacific basin (green), Indian basin (bleue), Indo-Pacific basin (bleue+green). 
+Note that semi-enclosed seas (Red, Med and Baltic seas) as well as Hudson Bay 
+are removed from the sub-basins. Note also that the Arctic Ocean has been split 
+into Atlantic and Pacific basins along the North fold line.  }
+\end{center}   \end{figure}  
+%>>>>>>>>>>>>>>>>>>>>>>>>>>>>
+
+
+% -----------------------------------------------------------
+%       CMIP specific diagnostics 
+% -----------------------------------------------------------
+\subsection{CMIP specific diagnostics (\mdl{diaar5})}
+
+A series of diagnostics has been added in the \mdl{diaar5}. 
+They corresponds to outputs that are required for AR5 simulations (CMIP5)
+(see also Section \ref{DIA_steric} for one of them). 
+Activating those outputs requires to define the \key{diaar5} CPP key.
+
+
+% -----------------------------------------------------------
+%       25 hour mean and hourly Surface, Mid and Bed 
+% -----------------------------------------------------------
+\subsection{25 hour mean output for tidal models }
+
+%------------------------------------------nam_dia25h-------------------------------------
+\namdisplay{nam_dia25h}
+%-----------------------------------------------------------------------------------------
+
+A module is available to compute a crudely detided M2 signal by obtaining a 25 hour mean.
+The 25 hour mean is available for daily runs by summing up the 25 hourly instantananeous hourly values from
+midnight at the start of the day to midight at the day end.
+This diagnostic is actived with the logical  $ln\_dia25h$
+
+
+% -----------------------------------------------------------
+%     Top Middle and Bed hourly output
+% -----------------------------------------------------------
+\subsection{Top Middle and Bed hourly output }
+
+%------------------------------------------nam_diatmb-----------------------------------------------------
+\namdisplay{nam_diatmb}
+%----------------------------------------------------------------------------------------------------------
+
+A module is available to output the surface (top), mid water and bed diagnostics of a set of standard variables. 
+This can be a useful diagnostic when hourly or sub-hourly output is required in high resolution tidal outputs.
+The tidal signal is retained but the overall data usage is cut to just three vertical levels. Also the bottom level 
+is calculated for each cell.
+This diagnostic is actived with the logical  $ln\_diatmb$
+
+
+
+% -----------------------------------------------------------
+%     Courant numbers
+% -----------------------------------------------------------
+\subsection{Courant numbers}
+Courant numbers provide a theoretical indication of the model's numerical stability. The advective Courant numbers can be calculated according to
+\begin{equation}
+\label{eq:CFL}
+C_u = |u|\frac{\rdt}{e_{1u}}, \quad C_v = |v|\frac{\rdt}{e_{2v}}, \quad C_w = |w|\frac{\rdt}{e_{3w}}
+\end{equation}
+in the zonal, meridional and vertical directions respectively. The vertical component is included although it is not strictly valid as the vertical velocity is calculated from the continuity equation rather than as a prognostic variable. Physically this represents the rate at which information is propogated across a grid cell. Values greater than 1 indicate that information is propagated across more than one grid cell in a single time step.
+
+The variables can be activated by setting the \np{nn\_diacfl} namelist parameter to 1 in the \ngn{namctl} namelist. The diagnostics will be written out to an ascii file named cfl\_diagnostics.ascii. In this file the maximum value of $C_u$, $C_v$, and $C_w$ are printed at each timestep along with the coordinates of where the maximum value occurs. At the end of the model run the maximum value of $C_u$, $C_v$, and $C_w$ for the whole model run is printed along with the coordinates of each. The maximum values from the run are also copied to the ocean.output file. 
+
+
 % ================================================================
 

branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_DOM.tex

@@ -486 +486 @@
 The last choice in terms of vertical coordinate concerns the presence (or not) in the model domain 
 of ocean cavities beneath ice shelves. Setting \np{ln\_isfcav} to true allows to manage ocean cavities, 
-otherwise they are filled in.
+otherwise they are filled in. This option is currently only available in $z$- or $zps$-coordinate,
+and partial step are also applied at the ocean/ice shelf interface. 
 
 Contrary to the horizontal grid, the vertical grid is computed in the code and no 
@@ -494 +495 @@
 \ifile{bathy\_meter} file, so that the computation of the number of wet ocean point 
 in each water column is by-passed}. 
-If \np{ln\_isfcav}~=~true, an extra file input file describing the ice shelf draft 
-(in meters) (\ifile{isf\_draft\_meter}) is needed.
-
 After reading the bathymetry, the algorithm for vertical grid definition differs 
 between the different options:
@@ -760 +758 @@
 as the minimum and maximum depths at which the terrain-following vertical coordinate is calculated.
 
-Options for stretching the coordinate are provided as examples, but care must be taken to ensure 
-that the vertical stretch used is appropriate for the application.
+Options for stretching the coordinate are provided as examples, but care must be taken 
+to ensure that the vertical stretch used is appropriate for the application.
 
 The original default NEMO s-coordinate stretching is available if neither of the other options 
@@ -772 +770 @@
 \end{equation}
 
-where $s_{min}$ is the depth at which the s-coordinate stretching starts and 
-allows a z-coordinate to placed on top of the stretched coordinate, 
-and z is the depth (negative down from the asea surface).
+where $s_{min}$ is the depth at which the $s$-coordinate stretching starts and 
+allows a $z$-coordinate to placed on top of the stretched coordinate, 
+and $z$ is the depth (negative down from the asea surface).
 
 \begin{equation}
@@ -886 +884 @@
 that do not communicate with another ocean point at the same level are eliminated.
 
-In case of ice shelf cavities, as for the representation of bathymetry, a 2D integer array, misfdep, is created. 
-misfdep defines the level of the first wet $t$-point (ie below the ice-shelf/ocean interface). All the cells between $k=1$ and $misfdep(i,j)-1$ are masked. 
-By default, $misfdep(:,:)=1$ and no cells are masked.
-Modifications of the model bathymetry and ice shelf draft into 
+As for the representation of bathymetry, a 2D integer array, misfdep, is created. 
+misfdep defines the level of the first wet $t$-point. All the cells between $k=1$ and $misfdep(i,j)-1$ are masked. 
+By default, misfdep(:,:)=1 and no cells are masked.
+
+In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into 
 the cavities are performed in the \textit{zgr\_isf} routine. The compatibility between ice shelf draft and bathymetry is checked. 
-All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded ($i.e.$ masked). 
 If only one cell on the water column is opened at $t$-, $u$- or $v$-points, the bathymetry or the ice shelf draft is dug to fit this constrain.
 If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked.\\ 
 
-From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows:
+From the \textit{mbathy} array, the mask fields are defined as follows:
 \begin{align*}
 tmask(i,j,k) &= \begin{cases}   \; 0&   \text{ if $k < misfdep(i,j) $ } \\
@@ -903 +901 @@
 vmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j+1,k)   \\
 fmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)   \\
-                   & \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) \\
+             &    \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) \\
 wmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j,k-1) \text{ with } wmask(i,j,1) = tmask(i,j,1) 
 \end{align*}
 
-Note, wmask is now defined. It allows, in case of ice shelves, 
-to deal with the top boundary (ice shelf/ocean interface) exactly in the same way as for the bottom boundary. 
+Note that, without ice shelves cavities, masks at $t-$ and $w-$points are identical with 
+the numerical indexing used (\S~\ref{DOM_Num_Index}). Nevertheless, $wmask$ are required 
+with oceean cavities to deal with the top boundary (ice shelf/ocean interface) 
+exactly in the same way as for the bottom boundary. 
 
 The specification of closed lateral boundaries requires that at least the first and last 
@@ -916 +916 @@
 (and so too the mask arrays) (see \S~\ref{LBC_jperio}).
 
-%%%
-\gmcomment{   \colorbox{yellow}{Add one word on tricky trick !} mbathy in further modified in zdfbfr{\ldots}.  }
-%%%
 
 % ================================================================

branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_DYN.tex

@@ -637 +637 @@
 ($e_{3w}$).
  
-$\bullet$ Traditional coding with adaptation for ice shelf cavities (\np{ln\_dynhpg\_isf}=true).
-This scheme need the activation of ice shelf cavities (\np{ln\_isfcav}=true).
-
 $\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}=true)
 
@@ -667 +664 @@
 
 $\bullet$ The ocean load is computed using the expression \eqref{Eq_dynhpg_sco} described in \ref{DYN_hpg_sco}. 
+A treatment of the partial cell for top and bottom similar to the one described in \ref{DYN_hpg_zps} is done 
+to reduce the residual circulation generated by the top partial cell. 
 
 %--------------------------------------------------------------------------------------------------------------

branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_LDF.tex

@@ -397 +397 @@
 \subsubsection{Space and Time Varying Mixing Coefficients}
 
-There is no default specification of space and time varying mixing coefficient. 
-The only case available is specific to the ORCA2 and ORCA05 global ocean 
-configurations. It provides only a tracer 
-mixing coefficient for eddy induced velocity (ORCA2) or both iso-neutral and 
-eddy induced velocity (ORCA05) that depends on the local growth rate of 
-baroclinic instability. This specification is actually used when an ORCA key 
+There are no default specifications of space and time varying mixing coefficient.  One
+available case is specific to the ORCA2 and ORCA05 global ocean configurations. It
+provides only a tracer mixing coefficient for eddy induced velocity (ORCA2) or both
+iso-neutral and eddy induced velocity (ORCA05) that depends on the local growth rate of
+baroclinic instability. This specification is actually used when an ORCA key
 and both \key{traldf\_eiv} and \key{traldf\_c2d} are defined.
+
+\subsubsection{Smagorinsky viscosity (\key{dynldf\_c3d} and \key{dynldf\_smag})}
+
+The \key{dynldf\_smag} key activates a 3D, time-varying viscosity that depends on the
+resolved motions. Following \citep{Smagorinsky_93} the viscosity coefficient is set
+proportional to a local deformation rate based on the horizontal shear and tension,
+namely:
+
+\begin{equation}
+A_{m_{Smag}} = \left(\frac{{\sf CM_{Smag}}}{\pi}\right)^2L^2\vert{D}\vert
+\end{equation}
+
+\noindent where the deformation rate $\vert{D}\vert$ is given by 
+
+\begin{equation}
+\vert{D}\vert=\sqrt{\left({\frac{\partial{u}} {\partial{x}}}
+                         -{\frac{\partial{v}} {\partial{y}}}\right)^2
+                 +  \left({\frac{\partial{u}} {\partial{y}}}
+                         +{\frac{\partial{v}} {\partial{x}}}\right)^2} 
+\end{equation}
+
+\noindent and $L$ is the local gridscale given by:
+
+\begin{equation}
+L^2 = \frac{2{e_1}^2 {e_2}^2}{\left ( {e_1}^2 + {e_2}^2 \right )}
+\end{equation}
+
+\citep{Griffies_Hallberg_MWR00} suggest values in the range 2.2 to 4.0 of the coefficient
+$\sf CM_{Smag}$ for oceanic flows. This value is set via the \np{rn\_cmsmag\_1} namelist
+parameter. An additional parameter: \np{rn\_cmsh} is included in NEMO for experimenting
+with the contribution of the shear term. A value of 1.0 (the default) calculates the
+deformation rate as above; a value of 0.0 will discard the shear term entirely.
+
+For numerical stability, the calculated viscosity is bounded according to the following:
+
+\begin{equation}
+{\rm MIN}\left ({ L^2\over {8\Delta{t}}}, rn\_ahm\_m\_lap\right ) \geq A_{m_{Smag}} 
+                                                                  \geq rn\_ahm\_0\_lap
+\end{equation}
+
+\noindent with both parameters for the upper and lower bounds being provided via the
+indicated namelist parameters.
+
+\bigskip When $ln\_dynldf\_bilap = .true.$, a biharmonic version of the Smagorinsky
+viscosity is also available which sets a coefficient for the biharmonic viscosity as:
+
+\begin{equation}
+B_{m_{Smag}} = - \left(\frac{{\sf CM_{bSmag}}}{\pi}\right)^2 {L^4\over 8}\vert{D}\vert
+\end{equation}
+
+\noindent which is bounded according to:
+
+\begin{equation}
+{\rm MAX}\left (-{ L^4\over {64\Delta{t}}}, rn\_ahm\_m\_blp\right ) \leq B_{m_{Smag}} 
+                                                                    \leq rn\_ahm\_0\_blp
+\end{equation}
+
+\noindent Note the reversal of the inequalities here because NEMO requires the biharmonic
+coefficients as negative numbers. $\sf CM_{bSmag}$ is set via the \np{rn\_cmsmag\_2}
+namelist parameter and the bounding values have corresponding entries in the namelist too.
+
+\bigskip The current implementation in NEMO also allows for 3D, time-varying diffusivities
+to be set using the Smagorinsky approach. Users should note that this option is not
+recommended for many applications since diffusivities will tend to be largest near
+boundaries (where shears are greatest) leading to spurious upwellings
+(\citep{Griffies_Bk04}, chapter 18.3.4). Nevertheless the option is there for those
+wishing to experiment. This choice requires both \key{traldf\_c3d} and \key{traldf\_smag}
+and uses the \np{rn\_chsmag} (${\sf CH_{Smag}}$), \np{rn\_smsh} and \np{rn\_aht\_m}
+namelist parameters in an analogous way to \np{rn\_cmsmag\_1}, \np{rn\_cmsh} and
+\np{rn\_ahm\_m\_lap} (see above) to set the diffusion coefficient:
+
+\begin{equation}
+A_{h_{Smag}} = \left(\frac{{\sf CH_{Smag}}}{\pi}\right)^2L^2\vert{D}\vert
+\end{equation}
+
+ 
+For numerical stability, the calculated diffusivity is bounded according to the following:
+
+\begin{equation}
+{\rm MIN}\left ({ L^2\over {8\Delta{t}}}, rn\_aht\_m\right ) \geq A_{h_{Smag}} 
+                                                             \geq rn\_aht\_0
+\end{equation}
+
 
 $\ $\newline    % force a new ligne

branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_SBC.tex

@@ -51 +51 @@
 \item the modification of fluxes below ice-covered areas (using observed ice-cover or a sea-ice model) (\np{nn\_ice}~=~0,1, 2 or 3) ; 
 \item the addition of river runoffs as surface freshwater fluxes or lateral inflow (\np{ln\_rnf}~=~true) ; 
-\item the addition of isf melting as lateral inflow (parameterisation) or as fluxes applied at the land-ice ocean interface (\np{ln\_isf}) ; 
+\item the addition of isf melting as lateral inflow (parameterisation) (\np{nn\_isf}~=~2 or 3 and \np{ln\_isfcav}~=~false) 
+or as fluxes applied at the land-ice ocean interface (\np{nn\_isf}~=~1 or 4 and \np{ln\_isfcav}~=~true) ; 
 \item the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2) ; 
 \item the transformation of the solar radiation (if provided as daily mean) into a diurnal cycle (\np{ln\_dm2dc}~=~true) ; 
@@ -128 +129 @@
 The ocean model provides, at each time step, to the surface module (\mdl{sbcmod}) 
 the surface currents, temperature and salinity.  
-These variables are averaged over \np{nf\_sbc} time-step (\ref{Tab_ssm}), 
+These variables are averaged over \np{nn\_fsbc} time-step (\ref{Tab_ssm}), 
 and it is these averaged fields which are used to computes the surface fluxes 
-at a frequency of \np{nf\_sbc} time-step.
+at a frequency of \np{nn\_fsbc} time-step.
 
 
@@ -144 +145 @@
 \caption{  \label{Tab_ssm}   
 Ocean variables provided by the ocean to the surface module (SBC). 
-The variable are averaged over nf{\_}sbc time step, $i.e.$ the frequency of 
-computation of surface fluxes.}
+The variable are averaged over nn{\_}fsbc time step, 
+$i.e.$ the frequency of computation of surface fluxes.}
 \end{center}   \end{table}
 %--------------------------------------------------------------------------------------------------------------
@@ -557 +558 @@
 reanalysis and satellite data. They use an inertial dissipative method to compute 
 the turbulent transfer coefficients (momentum, sensible heat and evaporation) 
-from the 10 metre wind speed, air temperature and specific humidity.
+from the 10 meters wind speed, air temperature and specific humidity.
 This \citet{Large_Yeager_Rep04} dataset is available through the 
 \href{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}{GFDL web site}. 
@@ -592 +593 @@
 or larger than the one of the input atmospheric fields.
 
+The \np{sn\_wndi}, \np{sn\_wndj}, \np{sn\_qsr}, \np{sn\_qlw}, \np{sn\_tair}, \np{sn\_humi},
+\np{sn\_prec}, \np{sn\_snow}, \np{sn\_tdif} parameters describe the fields 
+and the way they have to be used (spatial and temporal interpolations). 
+
+\np{cn\_dir} is the directory of location of bulk files
+\np{ln\_taudif} is the flag to specify if we use Hight Frequency (HF) tau information (.true.) or not (.false.)
+\np{rn\_zqt}: is the height of humidity and temperature measurements (m)
+\np{rn\_zu}: is the height of wind measurements (m)
+
+Three multiplicative factors are availables : 
+\np{rn\_pfac} and \np{rn\_efac} allows to adjust (if necessary) the global freshwater budget 
+by increasing/reducing the precipitations (total and snow) and or evaporation, respectively.
+The third one,\np{rn\_vfac}, control to which extend the ice/ocean velocities are taken into account 
+in the calculation of surface wind stress. Its range should be between zero and one, 
+and it is recommended to set it to 0.
+
 % -------------------------------------------------------------------------------------------------------------
 %        CLIO Bulk formulea
@@ -926 +943 @@
 \begin{description}
 \item[\np{nn\_isf}~=~1]
-The ice shelf cavity is represented (\np{ln\_isfcav}~=~true needed). The fwf and heat flux are computed. Two different bulk formula are available:
+The ice shelf cavities are explicitly represented. The fwf and heat flux are computed. Two different bulk formula are available:
    \begin{description}
    \item[\np{nn\_isfblk}~=~1]
@@ -934 +951 @@
    \item[\np{nn\_isfblk}~=~2] 
    The bulk formula used to compute the melt is based the one described in \citet{Jenkins1991}.
-        This formulation is based on a 3 equations formulation (a heat flux budget, a salt flux budget
-         and a linearised freezing point temperature equation).
+        This formulation is based on a 3 equations formulation (a heat flux budget, a salt flux budget and a linearised freezing point temperature equation).
    \end{description}
 
@@ -971 +987 @@
 
 \item[\np{nn\_isf}~=~4]
-The ice shelf cavity is opened (\np{ln\_isfcav}~=~true needed). However, the fwf is not computed but specified from file \np{sn\_fwfisf}). 
+The ice shelf cavity is opened. However, the fwf is not computed but specified from file \np{sn\_fwfisf}). 
 The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$.\\
 \end{description}
@@ -984 +1000 @@
 coarse to have realistic melting or for studies where you need to control your heat and fw input.\\ 
 
-A namelist parameters control over how many meters the heat and fw fluxes are spread. 
-\np{rn\_hisf\_tbl}] is the top boundary layer thickness as defined in \citet{Losch2008}. 
+Two namelist parameters control how the heat and fw fluxes are passed to NEMO: \np{rn\_hisf\_tbl} and \np{ln\_divisf}
+\begin{description}
+\item[\np{rn\_hisf\_tbl}] is the top boundary layer thickness as defined in \citet{Losch2008}. 
 This parameter is only used if \np{nn\_isf}~=~1 or \np{nn\_isf}~=~4
+It allows you to control over which depth you want to spread the heat and fw fluxes. 
 
 If \np{rn\_hisf\_tbl} = 0.0, the fluxes are put in the top level whatever is its tickness. 
 
-If \np{rn\_hisf\_tbl} $>$ 0.0, the fluxes are spread over the first \np{rn\_hisf\_tbl} m (ie over one or several cells).\\
-
-The ice shelf melt is implemented as a volume flux with in the same way as for the runoff.
+If \np{rn\_hisf\_tbl} $>$ 0.0, the fluxes are spread over the first \np{rn\_hisf\_tbl} m (ie over one or several cells).
+
+\item[\np{ln\_divisf}] is a flag to apply the fw flux as a volume flux or as a salt flux. 
+
+\np{ln\_divisf}~=~true applies the fwf as a volume flux. This volume flux is implemented with in the same way as for the runoff.
 The fw addition due to the ice shelf melting is, at each relevant depth level, added to the horizontal divergence 
 (\textit{hdivn}) in the subroutine \rou{sbc\_isf\_div}, called from \mdl{divcur}. 
 See the runoff section \ref{SBC_rnf} for all the details about the divergence correction. 
 
-
-\section{ Ice sheet coupling}
-\label{SBC_iscpl}
-%------------------------------------------namsbc_iscpl----------------------------------------------------
-\namdisplay{namsbc_iscpl}
-%--------------------------------------------------------------------------------------------------------
-Ice sheet/ocean coupling is done through file exchange at the restart step. NEMO, at each restart step, 
-read the bathymetry and ice shelf draft variable in a netcdf file. 
-If \np{ln\_iscpl = ~true}, the isf draft is assume to be different at each restart step 
-with potentially some new wet/dry cells due to the ice sheet dynamics/thermodynamics.
-The wetting and drying scheme applied on the restart is very simple and described below for the 6 different cases:
-\begin{description}
-\item[Thin a cell down:]
-   T/S/ssh are unchanged and U/V in the top cell are corrected to keep the barotropic transport (bt) constant ($bt_b=bt_n$).
-\item[Enlarge  a cell:]
-   See case "Thin a cell down"
-\item[Dry a cell:]
-   mask, T/S, U/V and ssh are set to 0. Furthermore, U/V into the water column are modified to satisfy ($bt_b=bt_n$).
-\item[Wet a cell:] 
-   mask is set to 1, T/S is extrapolated from neighbours, $ssh_n = ssh_b$ and U/V set to 0. If no neighbours along i,j and k, T/S/U/V and mask are set to 0.
-\item[Dry a column:]
-   mask, T/S, U/V are set to 0 everywhere in the column and ssh set to 0.
-\item[Wet a column:]
-   set mask to 1, T/S is extrapolated from neighbours, ssh is extrapolated from neighbours and U/V set to 0. If no neighbour, T/S/U/V and mask set to 0.
+\np{ln\_divisf}~=~false applies the fwf and heat flux directly on the salinity and temperature tendancy.
+
+\item[\np{ln\_conserve}] is a flag for \np{nn\_isf}~=~1. A conservative boundary layer scheme as described in \citet{Jenkins2001} 
+is used if \np{ln\_conserve}=true. It takes into account the fact that the melt water is at freezing T and needs to be warm up to ocean temperature. 
+It is only relevant for \np{ln\_divisf}~=~false. 
+If \np{ln\_divisf}~=~true, \np{ln\_conserve} has to be set to false to avoid a double counting of the contribution. 
+ 
 \end{description}
-The extrapolation is call \np{nn\_drown} times. It means that if the grounding line retreat by more than \np{nn\_drown} cells between 2 coupling steps,
- the code will be unable to fill all the new wet cells properly. The default number is set up for the MISOMIP idealised experiments.\\
-This coupling procedure is able to take into account grounding line and calving front migration. However, it is a non-conservative processe. 
-This could lead to a trend in heat/salt content and volume. In order to remove the trend and keep the conservation level as close to 0 as possible,
- a simple conservation scheme is available with \np{ln\_hsb = ~true}. The heat/salt/vol. gain/loss is diagnose, as well as the location. 
-Based on what is done on sbcrnf to prescribed a source of heat/salt/vol., the heat/salt/vol. gain/loss is removed/added,
- over a period of \np{rn\_fiscpl} time step, into the system. 
-So after \np{rn\_fiscpl} time step, all the heat/salt/vol. gain/loss due to extrapolation process is canceled.\\
-
-As the before and now fields are not compatible (modification of the geometry), the restart time step is prescribed to be an euler time step instead of a leap frog and $fields_b = fields_n$.
 %
 % ================================================================

branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_STO.tex

@@ -5 +5 @@
 \label{STO}
 
+Authors: P.-A. Bouttier
+
 \minitoc
 
+\newpage
 
-\newpage
-$\ $\newline    % force a new line
+
+The stochastic parametrization module aims to explicitly simulate uncertainties in the model. 
+More particularly, \cite{Brankart_OM2013} has shown that, 
+because of the nonlinearity of the seawater equation of state, unresolved scales represent 
+a major source of uncertainties in the computation of the large scale horizontal density gradient 
+(from T/S large scale fields), and that the impact of these uncertainties can be simulated 
+by random processes representing unresolved T/S fluctuations.
+
+The stochastic formulation of the equation of state can be written as:
+\begin{equation}
+ \label{eq:eos_sto}
+  \rho = \frac{1}{2} \sum_{i=1}^m\{ \rho[T+\Delta T_i,S+\Delta S_i,p_o(z)] + \rho[T-\Delta T_i,S-\Delta S_i,p_o(z)] \}
+\end{equation}
+where $p_o(z)$ is the reference pressure depending on the depth and, 
+$\Delta T_i$ and $\Delta S_i$ are a set of T/S perturbations defined as the scalar product 
+of the respective local T/S gradients with random walks $\mathbf{\xi}$:
+\begin{equation}
+ \label{eq:sto_pert}
+ \Delta T_i = \mathbf{\xi}_i \cdot \nabla T \qquad \hbox{and} \qquad \Delta S_i = \mathbf{\xi}_i \cdot \nabla S
+\end{equation}
+$\mathbf{\xi}_i$ are produced by a first-order autoregressive processes (AR-1) with 
+a parametrized decorrelation time scale, and horizontal and vertical standard deviations $\sigma_s$. 
+$\mathbf{\xi}$ are uncorrelated over the horizontal and fully correlated along the vertical.
+
+
+\section{Stochastic processes}
+\label{STO_the_details}
+
+The starting point of our implementation of stochastic parameterizations
+in NEMO is to observe that many existing parameterizations are based
+on autoregressive processes, which are used as a basic source of randomness
+to transform a deterministic model into a probabilistic model.
+A generic approach is thus to add one single new module in NEMO,
+generating processes with appropriate statistics
+to simulate each kind of uncertainty in the model
+(see \cite{Brankart_al_GMD2015} for more details).
+
+In practice, at every model grid point, independent Gaussian autoregressive
+processes~$\xi^{(i)},\,i=1,\ldots,m$ are first generated
+using the same basic equation:
+
+\begin{equation}
+\label{eq:autoreg}
+\xi^{(i)}_{k+1} = a^{(i)} \xi^{(i)}_k + b^{(i)} w^{(i)} + c^{(i)}
+\end{equation}
+
+\noindent
+where $k$ is the index of the model timestep; and
+$a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are parameters defining
+the mean ($\mu^{(i)}$) standard deviation ($\sigma^{(i)}$)
+and correlation timescale ($\tau^{(i)}$) of each process:
+
+\begin{itemize}
+\item for order~1 processes, $w^{(i)}$ is a Gaussian white noise,
+with zero mean and standard deviation equal to~1, and the parameters
+$a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are given by:
+
+\begin{equation}
+\label{eq:ord1}
+\left\{
+\begin{array}{l}
+a^{(i)} = \varphi \\
+b^{(i)} = \sigma^{(i)} \sqrt{ 1 - \varphi^2 } 
+ \qquad\qquad\mbox{with}\qquad\qquad
+\varphi = \exp \left( - 1 / \tau^{(i)} \right) \\
+c^{(i)} = \mu^{(i)} \left( 1 - \varphi \right) \\
+\end{array}
+\right.
+\end{equation}
+
+\item for order~$n>1$ processes, $w^{(i)}$ is an order~$n-1$ autoregressive process,
+with zero mean, standard deviation equal to~$\sigma^{(i)}$; correlation timescale
+equal to~$\tau^{(i)}$; and the parameters $a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are given by:
+
+\begin{equation}
+\label{eq:ord2}
+\left\{
+\begin{array}{l}
+a^{(i)} = \varphi \\
+b^{(i)} = \frac{n-1}{2(4n-3)} \sqrt{ 1 - \varphi^2 } 
+ \qquad\qquad\mbox{with}\qquad\qquad
+\varphi = \exp \left( - 1 / \tau^{(i)} \right) \\
+c^{(i)} = \mu^{(i)} \left( 1 - \varphi \right) \\
+\end{array}
+\right.
+\end{equation}
+
+\end{itemize}
+
+\noindent
+In this way, higher order processes can be easily generated recursively using 
+the same piece of code implementing Eq.~(\ref{eq:autoreg}), 
+and using succesively processes from order $0$ to~$n-1$ as~$w^{(i)}$.
+The parameters in Eq.~(\ref{eq:ord2}) are computed so that this recursive application
+of Eq.~(\ref{eq:autoreg}) leads to processes with the required standard deviation
+and correlation timescale, with the additional condition that
+the $n-1$ first derivatives of the autocorrelation function
+are equal to zero at~$t=0$, so that the resulting processes
+become smoother and smoother as $n$ is increased.
+
+Overall, this method provides quite a simple and generic way of generating 
+a wide class of stochastic processes. 
+However, this also means that new model parameters are needed to specify each of 
+these stochastic processes. As in any parameterization of lacking physics, 
+a very important issues then to tune these new parameters using either first principles, 
+model simulations, or real-world observations.
+
+\section{Implementation details}
+\label{STO_thech_details}
+
 %---------------------------------------namsbc--------------------------------------------------
 \namdisplay{namsto}
 %--------------------------------------------------------------------------------------------------------------
-$\ $\newline    % force a new ligne
 
+The computer code implementing stochastic parametrisations can be found in the STO directory.
+It involves three modules : 
+\begin{description}
+\item[\mdl{stopar}] : define the Stochastic parameters and their time evolution.
+\item[\mdl{storng}] : a random number generator based on (and includes) the 64-bit KISS 
+                      (Keep It Simple Stupid) random number generator distributed by George Marsaglia 
+                      (see \href{https://groups.google.com/forum/#!searchin/comp.lang.fortran/64-bit$20KISS$20RNGs}{here})
+\item[\mdl{stopts}] : stochastic parametrisation associated with the non-linearity of the equation of seawater, 
+ implementing Eq~\ref{eq:sto_pert} and specific piece of code in the equation of state implementing Eq~\ref{eq:eos_sto}.
+\end{description}
 
-See \cite{Brankart_OM2013} and \cite{Brankart_al_GMD2015} papers for a description of the parameterization.
+The \mdl{stopar} module has 3 public routines to be called by the model (in our case, NEMO):
+
+The first routine (\rou{sto\_par}) is a direct implementation of Eq.~(\ref{eq:autoreg}),
+applied at each model grid point (in 2D or 3D), 
+and called at each model time step ($k$) to update
+every autoregressive process ($i=1,\ldots,m$).
+This routine also includes a filtering operator, applied to $w^{(i)}$,
+to introduce a spatial correlation between the stochastic processes.
+
+The second routine (\rou{sto\_par\_init}) is an initialization routine mainly dedicated
+to the computation of parameters $a^{(i)}, b^{(i)}, c^{(i)}$
+for each autoregressive process, as a function of the statistical properties
+required by the model user (mean, standard deviation, time correlation,
+order of the process,\ldots). 
+
+Parameters for the processes can be specified through the following \ngn{namsto} namelist parameters:
+\begin{description}
+   \item[\np{nn\_sto\_eos}]   : number of independent random walks 
+   \item[\np{rn\_eos\_stdxy}] : random walk horz. standard deviation (in grid points)
+   \item[\np{rn\_eos\_stdz}]  : random walk vert. standard deviation (in grid points)
+   \item[\np{rn\_eos\_tcor}]  : random walk time correlation (in timesteps)
+   \item[\np{nn\_eos\_ord}]   : order of autoregressive processes
+   \item[\np{nn\_eos\_flt}]   : passes of Laplacian filter
+   \item[\np{rn\_eos\_lim}]   : limitation factor (default = 3.0)
+\end{description}
+This routine also includes the initialization (seeding) of the random number generator.
+
+The third routine (\rou{sto\_rst\_write}) writes a restart file (which suffix name is 
+given by \np{cn\_storst\_out} namelist parameter) containing the current value of 
+all autoregressive processes to allow restarting a simulation from where it has been interrupted.
+This file also contains the current state of the random number generator.
+When \np{ln\_rststo} is set to \textit{true}), the restart file (which suffix name is 
+given by \np{cn\_storst\_in} namelist parameter) is read by the initialization routine 
+(\rou{sto\_par\_init}). The simulation will continue exactly as if it was not interrupted
+only  when \np{ln\_rstseed} is set to \textit{true}, $i.e.$ when the state of 
+the random number generator is read in the restart file.

branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_TRA.tex

@@ -734 +734 @@
 (see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies)
 
-$\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt, (see \S\ref{SBC_isf} for further details 
-on how the ice shelf melt is computed and applied).
+$\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt, 
+(see \S\ref{SBC_isf} for further details on how the ice shelf melt is computed and applied).
 
 The surface boundary condition on temperature and salinity is applied as follows:
@@ -840 +840 @@
 ($i.e.$ the inverses of the extinction length scales) are tabulated over 61 nonuniform 
 chlorophyll classes ranging from 0.01 to 10 g.Chl/L (see the routine \rou{trc\_oce\_rgb} 
-in \mdl{trc\_oce} module). Three types of chlorophyll can be chosen in the RGB formulation:
-(1) a constant 0.05 g.Chl/L value everywhere (\np{nn\_chdta}=0) ; (2) an observed 
-time varying chlorophyll (\np{nn\_chdta}=1) ; (3) simulated time varying chlorophyll
-by TOP biogeochemical model (\np{ln\_qsr\_bio}=true). In the latter case, the RGB 
-formulation is used to calculate both the phytoplankton light limitation in PISCES 
-or LOBSTER and the oceanic heating rate. 
-
+in \mdl{trc\_oce} module). Four types of chlorophyll can be chosen in the RGB formulation:
+\begin{description} 
+\item[\np{nn\_chdta}=0] 
+a constant 0.05 g.Chl/L value everywhere ; 
+\item[\np{nn\_chdta}=1]  
+an observed time varying chlorophyll deduced from satellite surface ocean color measurement 
+spread uniformly in the vertical direction ; 
+\item[\np{nn\_chdta}=2]  
+same as previous case except that a vertical profile of chlorophyl is used. 
+Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value ;
+\item[\np{ln\_qsr\_bio}=true]  
+simulated time varying chlorophyll by TOP biogeochemical model. 
+In this case, the RGB formulation is used to calculate both the phytoplankton 
+light limitation in PISCES or LOBSTER and the oceanic heating rate. 
+\end{description} 
 The trend in \eqref{Eq_tra_qsr} associated with the penetration of the solar radiation 
 is added to the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}. 
@@ -1385 +1393 @@
                    I've changed "derivative" to "difference" and "mean" to "average"}
 
-With partial cells (\np{ln\_zps}=true) at bottom and top (\np{ln\_isfcav}=true), in general, tracers in horizontally 
-adjacent cells live at different depths. Horizontal gradients of tracers are needed 
-for horizontal diffusion (\mdl{traldf} module) and for the hydrostatic pressure 
-gradient (\mdl{dynhpg} module) to be active. The partial cell properties 
-at the top (\np{ln\_isfcav}=true) are computed in the same way as for the bottom. So, only the bottom interpolation is shown.
+With partial cells (\np{ln\_zps}=true) at bottom and top (\np{ln\_isfcav}=true), in general, 
+tracers in horizontally adjacent cells live at different depths. 
+Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module) 
+and for the hydrostatic pressure gradient (\mdl{dynhpg} module) to be active. 
 \gmcomment{STEVEN from gm : question: not sure of  what -to be active- means}
-
 Before taking horizontal gradients between the tracers next to the bottom, a linear 
 interpolation in the vertical is used to approximate the deeper tracer as if it actually 
@@ -1467 +1473 @@
 \gmcomment{gm :   this last remark has to be done}
 %%%
+
+If under ice shelf seas opened (\np{ln\_isfcav}=true), the partial cell properties 
+at the top are computed in the same way as for the bottom. Some extra variables are, 
+however, computed to reduce the flow generated at the top and bottom if $z*$ coordinates activated.
+The extra variables calculated and used by \S\ref{DYN_hpg_isf} are:
+
+$\bullet$ $\overline{T}_k^{\,i+1/2}$ as described in \eqref{Eq_zps_hde}
+
+$\bullet$ $\delta _{i+1/2} Z_{T_k} = \widetilde {Z}^{\,i}_{T_k}-Z^{\,i}_{T_k}$ to compute 
+the pressure gradient correction term used by \eqref{Eq_dynhpg_sco} in \S\ref{DYN_hpg_isf},
+ with $\widetilde {Z}_{T_k}$ the depth of the point $\widetilde {T}_{k}$ in case of $z^*$ coordinates 
+(this term = 0 in z-coordinates)

branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_ZDF.tex

@@ -262 +262 @@
 \end{equation}
 
-At the ocean surface, a non zero length scale is set through the  \np{rn\_lmin0} namelist 
+At the ocean surface, a non zero length scale is set through the  \np{rn\_mxl0} 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 
+leads to a 0.04~m, the default value of \np{rn\_mxl0}. 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}}$ ).
@@ -295 +295 @@
 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.
+to $\alpha_{CB} = 100$. Further setting  \np{ln\_mxl0} to true applies \eqref{ZDF_Lsbc} 
+as surface boundary condition on length scale, with $\beta$ hard coded to the 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.
@@ -852 +852 @@
 The bottom friction represents the friction generated by the bathymetry. 
 The top friction represents the friction generated by the ice shelf/ocean interface. 
-As the friction processes at the top and bottom are represented similarly, only the bottom friction is described in detail below.\\
+As the friction processes at the top and bottom are represented similarly, 
+only the bottom friction is described in detail below.
 
 
@@ -926 +927 @@
 $H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$. 
 This is the default value used in \NEMO. It corresponds to a decay time scale 
-of 115~days. It can be changed by specifying \np{rn\_bfric1} (namelist parameter).
+of 115~days. It can be changed by specifying \np{rn\_bfri1} (namelist parameter).
 
 For the linear friction case the coefficients defined in the general 
@@ -936 +937 @@
 \end{split}
 \end{equation}
-When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfric1}. 
+When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfri1}. 
 Setting \np{nn\_botfr}=0 is equivalent to setting $r=0$ and leads to a free-slip 
 bottom boundary condition. These values are assigned in \mdl{zdfbfr}. 
@@ -943 +944 @@
 in the \ifile{bfr\_coef} input NetCDF file. The mask values should vary from 0 to 1. 
 Locations with a non-zero mask value will have the friction coefficient increased 
-by $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfric1}.
+by $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri1}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -963 +964 @@
 $e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{Killworth1992} 
 uses $C_D = 1.4\;10^{-3}$ and $e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$. 
-The CME choices have been set as default values (\np{rn\_bfric2} and \np{rn\_bfeb2} 
+The CME choices have been set as default values (\np{rn\_bfri2} and \np{rn\_bfeb2} 
 namelist parameters).
 
@@ -978 +979 @@
 \end{equation}
 
-The coefficients that control the strength of the non-linear bottom friction are 
-initialised as namelist parameters: $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}. 
-Note for applications which treat tides explicitly a low or even zero value of 
-\np{rn\_bfeb2} is recommended. From v3.2 onwards a local enhancement of $C_D$ 
-is possible via an externally defined 2D mask array (\np{ln\_bfr2d}=true). 
-See previous section for details.
+The coefficients that control the strength of the non-linear bottom friction are
+initialised as namelist parameters: $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}.
+Note for applications which treat tides explicitly a low or even zero value of
+\np{rn\_bfeb2} is recommended. From v3.2 onwards a local enhancement of $C_D$ is possible
+via an externally defined 2D mask array (\np{ln\_bfr2d}=true).  This works in the same way
+as for the linear bottom friction case with non-zero masked locations increased by
+$mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri2}.
+
+% -------------------------------------------------------------------------------------------------------------
+%       Bottom Friction Log-layer
+% -------------------------------------------------------------------------------------------------------------
+\subsection{Log-layer Bottom Friction enhancement (\np{nn\_botfr} = 2, \np{ln\_loglayer} = .true.)}
+\label{ZDF_bfr_loglayer}
+
+In the non-linear bottom friction case, the drag coefficient, $C_D$, can be optionally
+enhanced using a "law of the wall" scaling. If  \np{ln\_loglayer} = .true., $C_D$ is no
+longer constant but is related to the thickness of the last wet layer in each column by:
+
+\begin{equation}
+C_D = \left ( {\kappa \over {\rm log}\left ( 0.5e_{3t}/rn\_bfrz0 \right ) } \right )^2
+\end{equation}
+
+\noindent where $\kappa$ is the von-Karman constant and \np{rn\_bfrz0} is a roughness
+length provided via the namelist.
+
+For stability, the drag coefficient is bounded such that it is kept greater or equal to
+the base \np{rn\_bfri2} value and it is not allowed to exceed the value of an additional
+namelist parameter: \np{rn\_bfri2\_max}, i.e.:
+
+\begin{equation}
+rn\_bfri2 \leq C_D \leq rn\_bfri2\_max
+\end{equation}
+
+\noindent Note also that a log-layer enhancement can also be applied to the top boundary
+friction if under ice-shelf cavities are in use (\np{ln\_isfcav}=.true.).  In this case, the
+relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2}
+and \np{rn\_tfri2\_max}.
 
 % -------------------------------------------------------------------------------------------------------------
@@ -1267 +1299 @@
 
 % ================================================================
+% Internal wave-driven mixing
+% ================================================================
+\section{Internal wave-driven mixing (\key{zdftmx\_new})}
+\label{ZDF_tmx_new}
+
+%--------------------------------------------namzdf_tmx_new------------------------------------------
+\namdisplay{namzdf_tmx_new}
+%--------------------------------------------------------------------------------------------------------------
+
+The parameterization of mixing induced by breaking internal waves is a generalization 
+of the approach originally proposed by \citet{St_Laurent_al_GRL02}. 
+A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed, 
+and the resulting diffusivity is obtained as 
+\begin{equation} \label{Eq_Kwave}
+A^{vT}_{wave} =  R_f \,\frac{ \epsilon }{ \rho \, N^2 }
+\end{equation}
+where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution 
+of the energy available for mixing. If the \np{ln\_mevar} namelist parameter is set to false, 
+the mixing efficiency is taken as constant and equal to 1/6 \citep{Osborn_JPO80}. 
+In the opposite (recommended) case, $R_f$ is instead a function of the turbulence intensity parameter 
+$Re_b = \frac{ \epsilon}{\nu \, N^2}$, with $\nu$ the molecular viscosity of seawater, 
+following the model of \cite{Bouffard_Boegman_DAO2013} 
+and the implementation of \cite{de_lavergne_JPO2016_efficiency}.
+Note that $A^{vT}_{wave}$ is bounded by $10^{-2}\,m^2/s$, a limit that is often reached when the mixing efficiency is constant.
+
+In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary 
+as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to true, a recommended choice). 
+This parameterization of differential mixing, due to \cite{Jackson_Rehmann_JPO2014}, 
+is implemented as in \cite{de_lavergne_JPO2016_efficiency}.
+
+The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$, is constructed 
+from three static maps of column-integrated internal wave energy dissipation, $E_{cri}(i,j)$, 
+$E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures 
+(de Lavergne et al., in prep):
+\begin{align*}
+F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\
+F_{pyc}(i,j,k) &\propto N^{n\_p}\\
+F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} }
+\end{align*} 
+In the above formula, $h_{ab}$ denotes the height above bottom, 
+$h_{wkb}$ denotes the WKB-stretched height above bottom, defined by
+\begin{equation*}
+h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz'  } \; ,
+\end{equation*}
+The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_tmx\_new} namelist)  controls the stratification-dependence of the pycnocline-intensified dissipation. 
+It can take values of 1 (recommended) or 2.
+Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of 
+the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps. 
+$h_{cri}$ is related to the large-scale topography of the ocean (etopo2) 
+and $h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of 
+the abyssal hill topography \citep{Goff_JGR2010} and the latitude.
+
+% ================================================================
+
+
+

branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/clean.sh

@@ -1 @@
-rm -f $( ls -1 NEMO_book.* | egrep -v "(tex|pdf)" ) *mpgraph*
+rm -f $( ls -1 ../NEMO_book.* | egrep -v "(tex|pdf)" ) *mpgraph*
 rm -f Chapters/*.aux Chapters/*.log