Changeset 6347

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

Files:

• DOC/TexFiles/Biblio/Biblio.bib (modified) (9 diffs)
• DOC/TexFiles/Chapters/Chap_DIA.tex (modified) (7 diffs)
• DOC/TexFiles/Chapters/Chap_DOM.tex (modified) (7 diffs)
• DOC/TexFiles/Chapters/Chap_DYN.tex (modified) (2 diffs)
• DOC/TexFiles/Chapters/Chap_LDF.tex (modified) (1 diff)
• DOC/TexFiles/Chapters/Chap_SBC.tex (modified) (9 diffs)
• DOC/TexFiles/Chapters/Chap_STO.tex (modified) (1 diff)
• DOC/TexFiles/Chapters/Chap_TRA.tex (modified) (4 diffs)
• DOC/TexFiles/Chapters/Chap_ZDF.tex (modified) (9 diffs)
• DOC/TexFiles/clean.sh (modified) (1 diff)
• NEMOGCM/CONFIG/SHARED/field_def.xml (modified) (1 diff)
• NEMOGCM/CONFIG/SHARED/namelist_ref (modified) (9 diffs)
• NEMOGCM/CONFIG/cfg.txt (modified) (1 diff)
• NEMOGCM/NEMO/LIM_SRC_3/limistate.F90 (modified) (4 diffs)
• NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90 (modified) (7 diffs)
• NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90 (modified) (1 diff)
• NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90 (modified) (1 diff)
• NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90 (modified) (3 diffs)
• NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90 (modified) (3 diffs)
• NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90 (modified) (8 diffs)
• NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90 (modified) (17 diffs)
• NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90 (modified) (1 diff)
• NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90 (modified) (2 diffs)
• NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90 (modified) (2 diffs)
• NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90 (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

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/CONFIG/SHARED/field_def.xml

@@ -479 +479 @@
         <!-- avt_tide: available with key_zdftmx -->
         <field id="av_tide"      long_name="tidal vertical diffusivity"   standard_name="ocean_vertical_tracer_diffusivity_due_to_tides"   unit="m2/s" />
+
+        <!-- fields available with key_zdftmx_new -->
+        <field id="av_wave"      long_name="wave-induced vertical diffusivity"     standard_name="ocean_vertical_tracer_diffusivity_due_to_internal_waves"         unit="m2/s" />
+        <field id="bn2"          long_name="squared Brunt-Vaisala frequency"       standard_name="squared_brunt_vaisala_frequency"                                 unit="s-1"  />
+        <field id="bflx_tmx"     long_name="wave-induced buoyancy flux"            standard_name="buoyancy_flux_due_to_internal_waves"                             unit="W/kg" />
+        <field id="pcmap_tmx"    long_name="power consumed by wave-driven mixing"  standard_name="vertically_integrated_power_consumption_by_wave_driven_mixing"   unit="W/m2"      grid_ref="grid_W_2D" />
+        <field id="emix_tmx"     long_name="power density available for mixing"    standard_name="power_available_for_mixing_from_breaking_internal_waves"         unit="W/kg" />
+        <field id="av_ratio"     long_name="S over T diffusivity ratio"            standard_name="salinity_over_temperature_diffusivity_ratio"                     unit="1"    />
 
         <!-- variables available with key_diaar5 -->   

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/CONFIG/SHARED/namelist_ref

@@ -11 +11 @@
 !!              6 - Tracer           (nameos, namtra_adv, namtra_ldf, namtra_ldfeiv, namtra_dmp)
 !!              7 - dynamics         (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf)
-!!              8 - Verical physics  (namzdf, namzdf_ric, namzdf_tke, namzdf_ddm, namzdf_tmx)
+!!              8 - Verical physics  (namzdf, namzdf_ric, namzdf_tke, namzdf_ddm, namzdf_tmx, namzdf_tmx_new)
 !!              9 - diagnostics      (namnc4, namtrd, namspr, namflo, namhsb, namsto)
 !!             10 - miscellaneous    (nammpp, namctl)
@@ -414 +414 @@
    ln_qsr_2bd  = .false.   !  2 bands              light penetration
    ln_qsr_bio  = .false.   !  bio-model light penetration
-   nn_chldta   =      1    !  RGB : Chl data (=1) or cst value (=0)
+   nn_chldta   =      1    !  RGB : 2D Chl data (=1), 3D Chl data (=2) or cst value (=0)
    rn_abs      =   0.58    !  RGB & 2 bands: fraction of light (rn_si1)
    rn_si0      =   0.35    !  RGB & 2 bands: shortess depth of extinction
@@ -516 +516 @@
 &namsbc_alb    !   albedo parameters
 !-----------------------------------------------------------------------
-   rn_cloud    =    0.06   !  cloud correction to snow and ice albedo
-   rn_albice   =    0.53   !  albedo of melting ice in the arctic and antarctic
-   rn_alphd    =    0.80   !  coefficients for linear interpolation used to
-   rn_alphc    =    0.65   !  compute albedo between two extremes values
-   rn_alphdi   =    0.72   !  (Pyane, 1972)
+   nn_ice_alb  =    0   !  parameterization of ice/snow albedo
+                        !     0: Shine & Henderson-Sellers (JGR 1985)
+                        !     1: "home made" based on Brandt et al. (J. Climate 2005)
+                        !                         and Grenfell & Perovich (JGR 2004)
+   rn_albice   =  0.53  !  albedo of bare puddled ice (values from 0.49 to 0.58)
+                        !     0.53 (default) => if nn_ice_alb=0
+                        !     0.50 (default) => if nn_ice_alb=1
 /
 !-----------------------------------------------------------------------
@@ -575 +577 @@
 !!======================================================================
 !!   namlbc        lateral momentum boundary condition
-!!   namobc        open boundaries parameters                           ("key_obc")
 !!   namagrif      agrif nested grid ( read by child model only )       ("key_agrif")
 !!   nambdy        Unstructured open boundaries                         ("key_bdy")
@@ -584 +585 @@
 &namlbc        !   lateral momentum boundary condition
 !-----------------------------------------------------------------------
+   rn_shlat    =    2.     !  shlat = 0  !  0 < shlat < 2  !  shlat = 2  !  2 < shlat
    !                       !  free slip  !   partial slip  !   no slip   ! strong slip
-   rn_shlat    =    2.     !  shlat = 0  !  0 < shlat < 2  !  shlat = 2  !  2 < shlat
    ln_vorlat   = .false.   !  consistency of vorticity boundary condition with analytical Eqs.
 /
@@ -600 +601 @@
 &nam_tide      !   tide parameters                                      ("key_tide")
 !-----------------------------------------------------------------------
-   ln_tide_pot   = .true.   !  use tidal potential forcing
-   ln_tide_ramp  = .false.  !
-   rdttideramp   =    0.    !
-   clname(1)     = 'DUMMY'  !  name of constituent - all tidal components must be set in namelist_cfg
+   ln_tide_pot   = .true.  !  use tidal potential forcing
+   ln_tide_ramp  = .false. !
+   rdttideramp   =    0.   !
+   clname(1)     = 'DUMMY' !  name of constituent - all tidal components must be set in namelist_cfg
 /
 !-----------------------------------------------------------------------
@@ -943 +944 @@
 !!    namzdf_ddm    double diffusive mixing parameterization            ("key_zdfddm")
 !!    namzdf_tmx    tidal mixing parameterization                       ("key_zdftmx")
+!!    namzdf_tmx_new    new internal wave mixing parameterization                       ("key_zdftmx")
 !!======================================================================
 !
@@ -1036 +1038 @@
    rn_tfe_itf  = 1.        !  ITF tidal dissipation efficiency
 /
+!-----------------------------------------------------------------------
+&namzdf_tmx_new !   internal wave-driven mixing parameterization        ("key_zdftmx_new" & "key_zdfddm")
+!-----------------------------------------------------------------------
+   nn_zpyc     = 1         !  pycnocline-intensified dissipation scales as N (=1) or N^2 (=2)
+   ln_mevar    = .true.    !  variable (T) or constant (F) mixing efficiency
+   ln_tsdiff   = .true.    !  account for differential T/S mixing (T) or not (F)
+/
+
 
 !!======================================================================
@@ -1083 +1093 @@
    nn_eos_flt   = 0        ! passes of Laplacian filter
    rn_eos_lim   = 2.0      ! limitation factor (default = 3.0)
+   !
    ln_rststo    = .false.  ! start from mean parameter (F) or from restart file (T)
-   ln_rstseed = .true.           ! read seed of RNG from restart file
-   cn_storst_in  = "restart_sto" !  suffix of stochastic parameter restart file (input)
-   cn_storst_out = "restart_sto" !  suffix of stochastic parameter restart file (output)
+   ln_rstseed   = .true.        ! read seed of RNG from restart file
+   cn_storst_in = "restart_sto" !  suffix of stochastic parameter restart file (input)
+   cn_storst_out= "restart_sto" !  suffix of stochastic parameter restart file (output)
 /
 

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/CONFIG/cfg.txt

@@ -7 +7 @@
 AMM12 OPA_SRC
 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC
-ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC
 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC
 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC
 GYRE OPA_SRC
+ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

@@ -273 +273 @@
                         !---------------------
                         ! Test 1: area conservation
-                        zA_cons = SUM(za_i_ini(ji,jj,:)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons )
-                        IF ( zconv .LT. 1.0e-6 ) THEN
+                        zA_cons = SUM( za_i_ini(ji,jj,:) )
+                        zconv   = ABS( zat_i_ini(ji,jj) - zA_cons )
+                        IF ( zconv < 1.0e-6 ) THEN
                            ztest_1 = 1
                         ELSE 
-                          !this write is useful
-                          IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(ji,jj) 
+                          ! this write is useful
+                          IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,   &
+                           &                                                    ' zat_i_ini = ',zat_i_ini(ji,jj)
                           ztest_1 = 0
                         ENDIF
 
                         ! Test 2: volume conservation
-                        zV_cons = SUM(zv_i_ini(ji,jj,:))
-                        zconv = ABS(zvt_i_ini(ji,jj) - zV_cons)
-
-                        IF( zconv .LT. 1.0e-6 ) THEN
+                        zV_cons = SUM( zv_i_ini(ji,jj,:) )
+                        zconv   = ABS( zvt_i_ini(ji,jj) - zV_cons )
+
+                        IF( zconv < 1.0e-6 ) THEN
                            ztest_2 = 1
                         ELSE
-                           !this write is useful
-                           IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
-                                                    ' zvt_i_ini = ', zvt_i_ini(ji,jj)
+                           ! this write is useful
+                           IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons,   &
+                              &                                                    ' zvt_i_ini = ', zvt_i_ini(ji,jj)
                            ztest_2 = 0
                         ENDIF
@@ -300 +302 @@
                         ELSE
                            ! this write is useful
-                           IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(ji,jj,i_fill) = ', &
-                           zh_i_ini(ji,jj,i_fill), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
+                           IF(lwp) WRITE(numout,*) ' * TEST3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** ',   &
+                              &    ' zh_i_ini(ji,jj,i_fill) = ', zh_i_ini(ji,jj,i_fill), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
                            IF(lwp) WRITE(numout,*) ' ji,jj,i_fill ',ji,jj,i_fill
                            IF(lwp) WRITE(numout,*) 'zht_i_ini ',zht_i_ini(ji,jj)
@@ -529 +531 @@
       !!-----------------------------------------------------------------------------
       !
-      REWIND( numnam_ice_ref )              ! Namelist namiceini in reference namelist : Ice initial state
+      REWIND( numnam_ice_ref )         ! Namelist namiceini in reference namelist : Ice initial state
       READ  ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp )
 
-      REWIND( numnam_ice_cfg )              ! Namelist namiceini in configuration namelist : Ice initial state
+      REWIND( numnam_ice_cfg )         ! Namelist namiceini in configuration namelist : Ice initial state
       READ  ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp )
       IF(lwm) WRITE ( numoni, namiceini )
 
-      slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts
-      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu
-      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi
-
-      ! Define the initial parameters
-      ! -------------------------
-
-      IF(lwp) THEN
+      slf_i(jp_hti) = sn_hti   ;   slf_i(jp_hts) = sn_hts
+      slf_i(jp_ati) = sn_ati   ;   slf_i(jp_tsu) = sn_tsu
+      slf_i(jp_tmi) = sn_tmi   ;   slf_i(jp_smi) = sn_smi
+
+      IF(lwp) THEN                     ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
@@ -571 +570 @@
             CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' )   ;   RETURN
          ENDIF
-
+         !
          DO ifpr = 1, jpfldi
             ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) )
             ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) )
          END DO
-
+         !
          ! fill si with slf_i and control print
          CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' )
-
+         !
          CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step
-
+         !
       ENDIF
-
+      !
    END SUBROUTINE lim_istate_init
 

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90

@@ -4 +4 @@
    !! Harmonic analysis of tidal constituents 
    !!======================================================================
-   !! History :  3.6  !  2014  (E O'Dea)  Original code
+   !! History :  3.0  !  07-2004  (A. Hines) New routine, developed from dia_wri_foam
+   !!            3.4  !  02-2013  (J. Siddorn) Routine taken from old dia_wri_foam
+   !!            3.6  !  11-2014  (E O'Dea)  Original code
    !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers variables
-   USE dom_oce         ! ocean space and time domain
-   USE in_out_manager  ! I/O units
-   USE iom             ! I/0 library
-   USE wrk_nemo        ! working arrays
-#if defined key_zdftke 
-   USE zdf_oce, ONLY: en
-#endif
-   USE zdf_oce, ONLY: avt, avm
-#if defined key_zdfgls
-   USE zdf_oce, ONLY: en
-   USE zdfgls, ONLY: mxln
-#endif
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain
+   USE zdf_oce        ! ocean vertical physics variables
+   USE zdfgls
+   !
+   USE in_out_manager ! I/O units
+   USE iom            ! I/0 library
+   USE wrk_nemo       ! working arrays
 
    IMPLICIT NONE
    PRIVATE
 
-   LOGICAL , PUBLIC ::   ln_dia25h     !:  25h mean output
-   PUBLIC   dia_25h_init               ! routine called by nemogcm.F90
-   PUBLIC   dia_25h                    ! routine called by diawri.F90
+   PUBLIC   dia_25h_init   ! routine called by nemogcm.F90
+   PUBLIC   dia_25h        ! routine called by diawri.F90
+
+   LOGICAL, PUBLIC ::   ln_dia25h     !:  25h mean output
 
   !! * variables for calculating 25-hourly means
+  REAL(wp) ::   r1_25 = 1._wp / 25.0_wp    ! factor for the mean calulation
    REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_25h  , sn_25h
    REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_25h 
@@ -39 +38 @@
 #endif
    INTEGER, SAVE :: cnt_25h     ! Counter for 25 hour means
-
-
 
    !!----------------------------------------------------------------------
@@ -56 +53 @@
       !!        
       !! ** Method : Read namelist
-      !!   History
-      !!   3.6  !  08-14  (E. O'Dea) Routine to initialize dia_25h
+      !!
       !!---------------------------------------------------------------------------
-      !!
       INTEGER ::   ios                 ! Local integer output status for namelist read
       INTEGER ::   ierror              ! Local integer for memory allocation
@@ -66 +61 @@
       !!----------------------------------------------------------------------
       !
-      REWIND ( numnam_ref )              ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics
+      REWIND ( numnam_ref )          ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics
       READ   ( numnam_ref, nam_dia25h, IOSTAT=ios, ERR= 901 )
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in reference namelist', lwp )
 
-      REWIND( numnam_cfg )              ! Namelist nam_dia25h in configuration namelist  25hour diagnostics
+      REWIND( numnam_cfg )           ! Namelist nam_dia25h in configuration namelist  25hour diagnostics
       READ  ( numnam_cfg, nam_dia25h, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist', lwp )
@@ -134 +129 @@
       ! ------------------------- !
       cnt_25h = 1  ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible) 
-      tn_25h(:,:,:) = tsb(:,:,:,jp_tem)
-      sn_25h(:,:,:) = tsb(:,:,:,jp_sal)
-      sshn_25h(:,:) = sshb(:,:)
-      un_25h(:,:,:) = ub(:,:,:)
-      vn_25h(:,:,:) = vb(:,:,:)
-      wn_25h(:,:,:) = wn(:,:,:)
-      avt_25h(:,:,:) = avt(:,:,:)
-      avm_25h(:,:,:) = avm(:,:,:)
-# if defined key_zdfgls || defined key_zdftke
-         en_25h(:,:,:) = en(:,:,:)
-#endif
-# if defined key_zdfgls
-         rmxln_25h(:,:,:) = mxln(:,:,:)
+      tn_25h  (:,:,:) = tsb (:,:,:,jp_tem)
+      sn_25h  (:,:,:) = tsb (:,:,:,jp_sal)
+      sshn_25h(:,:)   = sshb(:,:)
+      un_25h  (:,:,:) = ub  (:,:,:)
+      vn_25h  (:,:,:) = vb  (:,:,:)
+      wn_25h  (:,:,:) = wn  (:,:,:)
+      avt_25h (:,:,:) = avt (:,:,:)
+      avm_25h (:,:,:) = avm (:,:,:)
+# if defined key_zdfgls || defined key_zdftke
+      en_25h  (:,:,:) = en  (:,:,:)
+#endif
+# if defined key_zdfgls
+      rmxln_25h(:,:,:) = mxln(:,:,:)
 #endif
 #if defined key_lim3 || defined key_lim2
-         CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice')
+      CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice')
 #endif 
-
-      ! -------------------------- !
-      ! 3 - Return to dia_wri      !
-      ! -------------------------- !
-
-
+      !
    END SUBROUTINE dia_25h_init
 
@@ -163 +153 @@
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE dia_25h  ***
-      !!         
-      !!
-      !!--------------------------------------------------------------------
       !!                   
       !! ** Purpose :   Write diagnostics with M2/S2 tide removed
       !!
-      !! ** Method  :   
-      !!      25hr mean outputs for shelf seas
+      !! ** Method  :   25hr mean outputs for shelf seas
       !!
-      !! History :
-      !!   ?.0  !  07-04  (A. Hines) New routine, developed from dia_wri_foam
-      !!   3.4  !  02-13  (J. Siddorn) Routine taken from old dia_wri_foam
-      !!   3.6  !  08-14  (E. O'Dea) adapted for VN3.6
-      !!----------------------------------------------------------------------
-      !! * Modules used
-
-      IMPLICIT NONE
-
-      !! * Arguments
+      !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt      ! ocean time-step index
-
-
-      !! * Local declarations
+      !
       INTEGER ::   ji, jj, jk
-
       LOGICAL ::   ll_print = .FALSE.    ! =T print and flush numout
-      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi       ! temporary reals
-      INTEGER                          ::   i_steps                               ! no of timesteps per hour
-      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                    ! temporary workspace
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                                  ! temporary workspace
-      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                                 ! temporary workspace
-      INTEGER                          ::   iyear0, nimonth0,iday0                ! start year,imonth,day
-
+      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi   ! temporary reals
+      INTEGER                          ::   i_steps                           ! no of timesteps per hour
+      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                ! temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                              ! temporary workspace
+      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                             ! temporary workspace
+      INTEGER                          ::   iyear0, nimonth0,iday0            ! start year,imonth,day
       !!----------------------------------------------------------------------
 
@@ -216 +189 @@
       ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours
       ! every day
-      IF( MOD( kt, i_steps ) == 0  .and. kt .ne. nn_it000 ) THEN
-
-         IF (lwp) THEN
-              WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ',kt
-              WRITE(numout,*) '~~~~~~~~~~~~ '
+      IF( MOD( kt, i_steps ) == 0  .AND. kt /= nn_it000 ) THEN
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ', kt
+            WRITE(numout,*) '~~~~~~~~~~~~ '
          ENDIF
-
-         tn_25h(:,:,:)        = tn_25h(:,:,:) + tsn(:,:,:,jp_tem)
-         sn_25h(:,:,:)        = sn_25h(:,:,:) + tsn(:,:,:,jp_sal)
-         sshn_25h(:,:)        = sshn_25h(:,:) + sshn (:,:)
-         un_25h(:,:,:)        = un_25h(:,:,:) + un(:,:,:)
-         vn_25h(:,:,:)        = vn_25h(:,:,:) + vn(:,:,:)
-         wn_25h(:,:,:)        = wn_25h(:,:,:) + wn(:,:,:)
-         avt_25h(:,:,:)       = avt_25h(:,:,:) + avt(:,:,:)
-         avm_25h(:,:,:)       = avm_25h(:,:,:) + avm(:,:,:)
-# if defined key_zdfgls || defined key_zdftke
-         en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:)
-#endif
-# if defined key_zdfgls
-         rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + mxln(:,:,:)
+         !
+         tn_25h  (:,:,:) = tn_25h  (:,:,:) + tsn (:,:,:,jp_tem)
+         sn_25h  (:,:,:) = sn_25h  (:,:,:) + tsn (:,:,:,jp_sal)
+         sshn_25h(:,:)   = sshn_25h(:,:)   + sshn(:,:)
+         un_25h  (:,:,:) = un_25h  (:,:,:) + un  (:,:,:)
+         vn_25h  (:,:,:) = vn_25h  (:,:,:) + vn  (:,:,:)
+         wn_25h  (:,:,:) = wn_25h  (:,:,:) + wn  (:,:,:)
+         avt_25h (:,:,:) = avt_25h (:,:,:) + avt (:,:,:)
+         avm_25h (:,:,:) = avm_25h (:,:,:) + avm (:,:,:)
+# if defined key_zdfgls || defined key_zdftke
+         en_25h  (:,:,:) = en_25h  (:,:,:) + en  (:,:,:)
+#endif
+# if defined key_zdfgls
+         rmxln_25h(:,:,:) = rmxln_25h(:,:,:) + mxln(:,:,:)
 #endif
          cnt_25h = cnt_25h + 1
-
-         IF (lwp) THEN
-            WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far', cnt_25h
          ENDIF
-
+         !
       ENDIF ! MOD( kt, i_steps ) == 0
-
-         ! Write data for 25 hour mean output streams
-      IF( cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN
-
-            IF(lwp) THEN
-               WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt
-               WRITE(numout,*) '~~~~~~~~~~~~ '
-            ENDIF
-
-            tn_25h(:,:,:)        = tn_25h(:,:,:) / 25.0_wp
-            sn_25h(:,:,:)        = sn_25h(:,:,:) / 25.0_wp
-            sshn_25h(:,:)        = sshn_25h(:,:) / 25.0_wp
-            un_25h(:,:,:)        = un_25h(:,:,:) / 25.0_wp
-            vn_25h(:,:,:)        = vn_25h(:,:,:) / 25.0_wp
-            wn_25h(:,:,:)        = wn_25h(:,:,:) / 25.0_wp
-            avt_25h(:,:,:)       = avt_25h(:,:,:) / 25.0_wp
-            avm_25h(:,:,:)       = avm_25h(:,:,:) / 25.0_wp
-# if defined key_zdfgls || defined key_zdftke
-            en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp
-#endif
-# if defined key_zdfgls
-            rmxln_25h(:,:,:)       = rmxln_25h(:,:,:) / 25.0_wp
-#endif
-
-            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output'
-            zmdi=1.e+20 !missing data indicator for masking
-            ! write tracers (instantaneous)
-            zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put("temper25h", zw3d)   ! potential temperature
-            zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put( "salin25h", zw3d  )   ! salinity
-            zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
-            CALL iom_put( "ssh25h", zw2d )   ! sea surface 
-
-
-            ! Write velocities (instantaneous)
-            zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
-            CALL iom_put("vozocrtx25h", zw3d)    ! i-current
-            zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
-            CALL iom_put("vomecrty25h", zw3d  )   ! j-current
-
-            zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put("vomecrtz25h", zw3d )   ! k-current
-            zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put("avt25h", zw3d )   ! diffusivity
-            zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put("avm25h", zw3d)   ! viscosity
+      !
+      ! Write data for 25 hour mean output streams
+      IF( cnt_25h == 25 .AND. MOD( kt, i_steps*24) == 0 .AND. kt /= nn_it000 ) THEN
+         !
+         IF(lwp) THEN
+            WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt
+            WRITE(numout,*) '~~~~~~~~~~~~ '
+         ENDIF
+         tn_25h  (:,:,:) = tn_25h  (:,:,:) * r1_25
+         sn_25h  (:,:,:) = sn_25h  (:,:,:) * r1_25
+         sshn_25h(:,:)   = sshn_25h(:,:)   * r1_25
+         un_25h  (:,:,:) = un_25h  (:,:,:) * r1_25
+         vn_25h  (:,:,:) = vn_25h  (:,:,:) * r1_25
+         wn_25h  (:,:,:) = wn_25h  (:,:,:) * r1_25
+         avt_25h (:,:,:) = avt_25h (:,:,:) * r1_25
+         avm_25h (:,:,:) = avm_25h (:,:,:) * r1_25
+# if defined key_zdfgls || defined key_zdftke
+         en_25h  (:,:,:)        = en_25h(:,:,:) * r1_25
+#endif
+# if defined key_zdfgls
+         rmxln_25h(:,:,:)       = rmxln_25h(:,:,:) * r1_25
+#endif
+
+         IF(lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output'
+         zmdi=1.e+20 !missing data indicator for masking
+         ! write tracers (instantaneous)
+         zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         CALL iom_put("temper25h", zw3d)   ! potential temperature
+         zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         CALL iom_put( "salin25h", zw3d  )   ! salinity
+         zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1))
+         CALL iom_put( "ssh25h", zw2d )   ! sea surface 
+         !
+         ! Write velocities (instantaneous)
+         zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:))
+         CALL iom_put("vozocrtx25h", zw3d)    ! i-current
+         zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:))
+         CALL iom_put("vomecrty25h", zw3d  )   ! j-current
+         zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         CALL iom_put("vomecrtz25h", zw3d )   ! k-current
+         zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         CALL iom_put("avt25h", zw3d )   ! diffusivity
+         zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         CALL iom_put("avm25h", zw3d)   ! viscosity
 #if defined key_zdftke || defined key_zdfgls 
-            zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put("tke25h", zw3d)   ! tke
+         zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         CALL iom_put("tke25h", zw3d)   ! tke
 #endif
 #if defined key_zdfgls 
-            zw3d(:,:,:) = rmxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
-            CALL iom_put( "mxln25h",zw3d)
-#endif
-
-            ! After the write reset the values to cnt=1 and sum values equal current value 
-            tn_25h(:,:,:) = tsn(:,:,:,jp_tem)
-            sn_25h(:,:,:) = tsn(:,:,:,jp_sal)
-            sshn_25h(:,:) = sshn (:,:)
-            un_25h(:,:,:) = un(:,:,:)
-            vn_25h(:,:,:) = vn(:,:,:)
-            wn_25h(:,:,:) = wn(:,:,:)
-            avt_25h(:,:,:) = avt(:,:,:)
-            avm_25h(:,:,:) = avm(:,:,:)
-# if defined key_zdfgls || defined key_zdftke
-            en_25h(:,:,:) = en(:,:,:)
-#endif
-# if defined key_zdfgls
-            rmxln_25h(:,:,:) = mxln(:,:,:)
-#endif
-            cnt_25h = 1
-            IF (lwp)  WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average',cnt_25h
-
+         zw3d(:,:,:) = rmxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:))
+         CALL iom_put( "mxln25h",zw3d)
+#endif
+         !
+         ! After the write reset the values to cnt=1 and sum values equal current value 
+         tn_25h  (:,:,:) = tsn (:,:,:,jp_tem)
+         sn_25h  (:,:,:) = tsn (:,:,:,jp_sal)
+         sshn_25h(:,:)   = sshn(:,:)
+         un_25h  (:,:,:) = un  (:,:,:)
+         vn_25h  (:,:,:) = vn  (:,:,:)
+         wn_25h  (:,:,:) = wn  (:,:,:)
+         avt_25h (:,:,:) = avt (:,:,:)
+         avm_25h (:,:,:) = avm (:,:,:)
+# if defined key_zdfgls || defined key_zdftke
+         en_25h  (:,:,:) = en  (:,:,:)
+#endif
+# if defined key_zdfgls
+         rmxln_25h(:,:,:) = mxln(:,:,:)
+#endif
+         cnt_25h = 1
+         IF(lwp)  WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value ',   &
+            &                                    'and cnt_25h to one for overlapping average',cnt_25h
+         !
       ENDIF !  cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000
-
-
+      !
    END SUBROUTINE dia_25h 
 

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90

@@ -151 +151 @@
             WRITE(numout,*) 'dia_cfl     : Maximum Courant number information for the run:'
             WRITE(numout,*) '~~~~~~~~~~~~'
-            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cu', cu_max, 'at (i, j, k) = (', cu_loc(1), cu_loc(2), cu_loc(3), ')'
+            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cu', cu_max,   &
+               &                                                           'at (i,j,k) = (', cu_loc(1), cu_loc(2), cu_loc(3), ')'
             WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max)
-            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cv', cv_max, 'at (i, j, k) = (', cv_loc(1), cv_loc(2), cv_loc(3), ')'
+            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cv', cv_max,   &
+               &                                                           'at (i,j,k) = (', cv_loc(1), cv_loc(2), cv_loc(3), ')'
             WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max)
-            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cw', cw_max, 'at (i, j, k) = (', cw_loc(1), cw_loc(2), cw_loc(3), ')'
+            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cw', cw_max,   &
+               &                                                           'at (i,j,k) = (', cw_loc(1), cw_loc(2), cw_loc(3), ')'
             WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max)
 

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90

@@ -142 +142 @@
 
       ! control print
-      IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i8,a,i8,a,i8,a,i8)')' =======>> 1/2 time step before the start of the run DATE Y/M/D = ',   &
-           &                   nyear, '/', nmonth, '/', nday, '  nsec_day:', nsec_day, '  nsec_week:', nsec_week, '  &
-           &                   nsec_month:', nsec_month , '  nsec_year:' , nsec_year
+      IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i8,a,i8,a,i8,a,i8)')                                  &
+         &  ' =======>> 1/2 time step before the start of the run DATE Y/M/D = ', nyear,            &
+         &             '/', nmonth, '/', nday, '  nsec_day:', nsec_day,'  nsec_week:', nsec_week,   &
+         &             '  nsec_month:', nsec_month ,'  nsec_year:' , nsec_year
 
       ! Up to now, calendar parameters are related to the end of previous run (nit000-1)

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/DOM/iscplrst.F90

@@ -215 +215 @@
       ! --------------------------------------
          DO jk = 1,jpk
-            DO jj=1,jpj
-               DO ji=1,jpi
-                  IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
-                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk) * ( 1._wp + sshn(ji,jj) / ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN
+                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk)       &
+                        &            * ( 1._wp + sshn(ji,jj) / ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) )
                   ENDIF
                END DO
@@ -386 +387 @@
       !PM: Is this IF needed since change to VVL by default
       IF (.NOT.ln_linssh) THEN
-         DO jk = 2,jpk-1
-            DO jj = 1,jpj
-               DO ji = 1,jpi
-                  IF (zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND. (tmask(ji,jj,1)==0._wp .OR. ptmask_b(ji,jj,1)==0._wp) ) THEN
+         DO jk = 2, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND.   &
+                     & ( tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp )  ) THEN
                      !compute weight
                      zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1))
@@ -436 +438 @@
       CALL wrk_dealloc(jpi,jpj,       zbub   , zbvb    , zbun  , zbvn        ) 
       CALL wrk_dealloc(jpi,jpj,       zssh0  , zssh1  , zhu1 , zhv1          ) 
-
+      !
    END SUBROUTINE iscpl_rst_interpol
 
+   !!======================================================================
 END MODULE iscplrst

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

@@ -764 +764 @@
          istart(idmspc+1) = itime
 
-         IF( PRESENT(kstart) .AND. .NOT. ll_depth_spec ) THEN ; istart(1:idmspc) = kstart(1:idmspc) ; icnt(1:idmspc) = kcount(1:idmspc)
+         IF( PRESENT(kstart) .AND. .NOT.ll_depth_spec ) THEN 
+            istart(1:idmspc) = kstart(1:idmspc) 
+            icnt  (1:idmspc) = kcount(1:idmspc)
          ELSE
-            IF(           idom == jpdom_unknown ) THEN                                                ; icnt(1:idmspc) = idimsz(1:idmspc)
+            IF(          idom == jpdom_unknown ) THEN
+               icnt(1:idmspc) = idimsz(1:idmspc)
             ELSE 
                IF( .NOT. PRESENT(pv_r1d) ) THEN   !   not a 1D array
-                  IF(     idom == jpdom_data    ) THEN
+                  IF(    idom == jpdom_data    ) THEN
                      jstartrow = 1
                      IF( luse_jattr ) THEN
@@ -776 +779 @@
                      ENDIF
                      istart(1:2) = (/ mig(1), mjg(1) + jstartrow - 1 /)  ! icnt(1:2) done below
-                  ELSEIF( idom == jpdom_global  ) THEN ; istart(1:2) = (/ nimpp , njmpp  /)  ! icnt(1:2) done below
+                  ELSEIF( idom == jpdom_global  ) THEN   ;   istart(1:2) = (/ nimpp , njmpp  /)  ! icnt(1:2) done below
                   ENDIF
                   ! we do not read the overlap                     -> we start to read at nldi, nldj
 ! JMM + SM: ugly patch before getting the new version of lib_mpp)
 !                  IF( idom /= jpdom_local_noovlap )   istart(1:2) = istart(1:2) + (/ nldi - 1, nldj - 1 /)
-                  IF( llnoov .AND. idom /= jpdom_local_noovlap ) istart(1:2) = istart(1:2) + (/ nldi - 1, nldj - 1 /)
+                  IF( llnoov .AND. idom /= jpdom_local_noovlap )   istart(1:2) = istart(1:2) + (/ nldi - 1, nldj - 1 /)
                   ! we do not read the overlap and the extra-halos -> from nldi to nlei and from nldj to nlej 
 ! JMM + SM: ugly patch before getting the new version of lib_mpp)
@@ -789 +792 @@
                   ENDIF
                   IF( PRESENT(pv_r3d) ) THEN
-                     IF( idom == jpdom_data ) THEN                                  ; icnt(3) = jpkdta
-                     ELSE IF( ll_depth_spec .AND. PRESENT(kstart) ) THEN            ; istart(3) = kstart(3); icnt(3) = kcount(3)
-                     ELSE                                                           ; icnt(3) = jpk
+                     IF( idom == jpdom_data ) THEN                         ;   icnt(3) = jpkdta
+                     ELSE IF( ll_depth_spec .AND. PRESENT(kstart) ) THEN   ;   istart(3) = kstart(3)   ;   icnt(3) = kcount(3)
+                     ELSE                                                  ;   icnt(3) = jpk
                      ENDIF
                   ENDIF

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

@@ -8 +8 @@
    !!             -   ! 2004-11  (C. Talandier)  add albedo_init
    !!             -   ! 2001-06  (M. Vancoppenolle) LIM 3.0
-   !!             -   ! 2006-08  (G. Madec)  cleaning for surface module
+   !!            3.2  ! 2006-08  (G. Madec)  cleaning for surface module
+   !!            3.6  ! 2016-01  (C. Rousset) new parameterization for sea ice albedo
    !!----------------------------------------------------------------------
 
@@ -17 +18 @@
    !!----------------------------------------------------------------------
    USE phycst         ! physical constants
+   !
    USE in_out_manager ! I/O manager
    USE lib_mpp        ! MPP library
@@ -25 +27 @@
    PRIVATE
 
-   PUBLIC   albedo_ice   ! routine called sbcice_lim.F90
-   PUBLIC   albedo_oce   ! routine called by ???
-
-   INTEGER  ::   albd_init = 0      !: control flag for initialization
-   REAL(wp) ::   zzero     = 0.e0   ! constant values
-   REAL(wp) ::   zone      = 1.e0   !    "       "
-
-   REAL(wp) ::   c1     = 0.05    ! constants values
-   REAL(wp) ::   c2     = 0.10    !    "        "
-   REAL(wp) ::   rmue   = 0.40    !  cosine of local solar altitude
-
+   PUBLIC   albedo_ice   ! called sbcice_lim.F90
+   PUBLIC   albedo_oce   ! called by sbccpl.F90
+
+   INTEGER  ::   albd_init = 0          ! control flag for initialization
+  
+   REAL(wp) ::   rmue      = 0.400_wp   ! cosine of local solar altitude
+   REAL(wp) ::   ralb_oce  = 0.066_wp   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001)
+   REAL(wp) ::   c1        = 0.050_wp   ! snow thickness (only for nn_ice_alb=0)
+   REAL(wp) ::   c2        = 0.100_wp   !  "        "
+   REAL(wp) ::   rcloud    = 0.060_wp   ! cloud effect on albedo (only-for nn_ice_alb=0)
+ 
    !                             !!* namelist namsbc_alb
-   REAL(wp) ::   rn_cloud         !  cloudiness effect on snow or ice albedo (Grenfell & Perovich, 1984)
-#if defined key_lim3
-   REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
-#else
-   REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers)
-#endif
-   REAL(wp) ::   rn_alphd         !  coefficients for linear interpolation used to compute
-   REAL(wp) ::   rn_alphdi        !  albedo between two extremes values (Pyane, 1972)
-   REAL(wp) ::   rn_alphc         ! 
+   INTEGER  ::   nn_ice_alb       ! type of albedo parameterization 
+   REAL(wp) ::   rn_albice        ! albedo of bare puddled ice
 
    !!----------------------------------------------------------------------
@@ -59 +54 @@
       !!          
       !! ** Purpose :   Computation of the albedo of the snow/ice system 
-      !!                as well as the ocean one
       !!       
-      !! ** Method  : - Computation of the albedo of snow or ice (choose the 
-      !!                rignt one by a large number of tests
-      !!              - Computation of the albedo of the ocean
-      !!
-      !! References :   Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250.
+      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb)
+      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies
+      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005)
+      !!                                                                           and Grenfell & Perovich (JGR 2004)
+      !!                Description of scheme 1:
+      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005)
+      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999
+      !!                     0-5cm  : linear function of ice thickness
+      !!                     5-150cm: log    function of ice thickness
+      !!                     > 150cm: constant
+      !!                  2) Albedo dependency on snow thickness follows the findings from Grenfell & Perovich (2004)
+      !!                     i.e. it increases as -EXP(-snw_thick/0.02) during freezing and -EXP(-snw_thick/0.03) during melting
+      !!                  3) Albedo dependency on clouds is speculated from measurements of Grenfell and Perovich (2004)
+      !!                     i.e. cloudy-clear albedo depend on cloudy albedo following a 2d order polynomial law
+      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice
+      !!
+      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions:
+      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo
+      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger 
+      !!                     under melting conditions than under freezing conditions
+      !!                  3) the evolution of ice albedo as a function of ice thickness shows  
+      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic
+      !!
+      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250.
+      !!                Brandt et al. 2005, J. Climate, vol 18
+      !!                Grenfell & Perovich 2004, JGR, vol 109 
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin)
@@ -73 +88 @@
       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os   !  albedo of ice under overcast sky
       !!
-      INTEGER  ::   ji, jj, jl    ! dummy loop indices
-      INTEGER  ::   ijpl          ! number of ice categories (3rd dim of ice input arrays)
-      REAL(wp) ::   zalbpsnm      ! albedo of ice under clear sky when snow is melting
-      REAL(wp) ::   zalbpsnf      ! albedo of ice under clear sky when snow is freezing
-      REAL(wp) ::   zalbpsn       ! albedo of snow/ice system when ice is coverd by snow
-      REAL(wp) ::   zalbpic       ! albedo of snow/ice system when ice is free of snow
-      REAL(wp) ::   zithsn        ! = 1 for hsn >= 0 ( ice is cov. by snow ) ; = 0 otherwise (ice is free of snow)
-      REAL(wp) ::   zitmlsn       ! = 1 freezinz snow (pt_ice >=rt0_snow) ; = 0 melting snow (pt_ice<rt0_snow)
-      REAL(wp) ::   zihsc1        ! = 1 hsn <= c1 ; = 0 hsn > c1
-      REAL(wp) ::   zihsc2        ! = 1 hsn >= c2 ; = 0 hsn < c2
-      !!
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalbfz    ! = rn_alphdi for freezing ice ; = rn_albice for melting ice
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zficeth   !  function of ice thickness
+      INTEGER  ::   ji, jj, jl   ! dummy loop indices
+      INTEGER  ::   ijpl         ! number of ice categories (3rd dim of ice input arrays)
+      REAL(wp) ::   ralb_im, ralb_sf, ralb_sm, ralb_if
+      REAL(wp) ::   zswitch, z1_c1, z1_c2
+      REAL(wp) ::   zalb_sm, zalb_sf, zalb_st   ! albedo of snow melting, freezing, total
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalb_it   ! intermediate variable & albedo of ice (snow free)
       !!---------------------------------------------------------------------
+
+      ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
       
-      ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
-
-      CALL wrk_alloc( jpi,jpj,ijpl, zalbfz, zficeth )
+      CALL wrk_alloc( jpi,jpj,ijpl,   zalb, zalb_it )
 
       IF( albd_init == 0 )   CALL albedo_init      ! initialization 
 
-      !---------------------------
-      !  Computation of  zficeth
-      !---------------------------
-      ! ice free of snow and melts
-      WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalbfz(:,:,:) = rn_albice
-      ELSE WHERE                                              ;   zalbfz(:,:,:) = rn_alphdi
-      END  WHERE
-
-      WHERE     ( 1.5  < ph_ice                     )  ;  zficeth = zalbfz
-      ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zficeth = 0.472  + 2.0 * ( zalbfz - 0.472 ) * ( ph_ice - 1.0 )
-      ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zficeth = 0.2467 + 0.7049 * ph_ice              &
-         &                                                                 - 0.8608 * ph_ice * ph_ice     &
-         &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
-      ELSE WHERE                                       ;  zficeth = 0.1    + 3.6    * ph_ice
-      END WHERE
-
-!!gm old code
-!      DO jl = 1, ijpl
-!         DO jj = 1, jpj
-!            DO ji = 1, jpi
-!               IF( ph_ice(ji,jj,jl) > 1.5 ) THEN
-!                  zficeth(ji,jj,jl) = zalbfz(ji,jj,jl)
-!               ELSEIF( ph_ice(ji,jj,jl) > 1.0  .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN
-!                  zficeth(ji,jj,jl) = 0.472 + 2.0 * ( zalbfz(ji,jj,jl) - 0.472 ) * ( ph_ice(ji,jj,jl) - 1.0 )
-!               ELSEIF( ph_ice(ji,jj,jl) > 0.05 .AND. ph_ice(ji,jj,jl) <= 1.0 ) THEN
-!                  zficeth(ji,jj,jl) = 0.2467 + 0.7049 * ph_ice(ji,jj,jl)                               &
-!                     &                    - 0.8608 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl)                 &
-!                     &                    + 0.3812 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl) * ph_ice (ji,jj,jl)
-!               ELSE
-!                  zficeth(ji,jj,jl) = 0.1 + 3.6 * ph_ice(ji,jj,jl) 
-!               ENDIF
-!            END DO
-!         END DO
-!      END DO
-!!gm end old code
       
-      !----------------------------------------------- 
-      !    Computation of the snow/ice albedo system 
-      !-------------------------- ---------------------
+      SELECT CASE ( nn_ice_alb )    ! Select the albedo parameterization
+      !
+      !              !------------------------------------------
+      CASE( 0 )      !  Shine and Henderson-Sellers (1985)
+         !           !------------------------------------------
+         !
+         ralb_sf = 0.80       ! dry snow
+         ralb_sm = 0.65       ! melting snow
+         ralb_if = 0.72       ! bare frozen ice
+         ralb_im = rn_albice  ! bare puddled ice 
+         ! 
+         !  Computation of ice albedo (free of snow)
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im
+         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if
+         END  WHERE
       
-      !    Albedo of snow-ice for clear sky.
-      !-----------------------------------------------    
-      DO jl = 1, ijpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               !  Case of ice covered by snow.             
-               !                                        !  freezing snow        
-               zihsc1   = 1.0 - MAX( zzero , SIGN( zone , - ( ph_snw(ji,jj,jl) - c1 ) ) )
-               zalbpsnf = ( 1.0 - zihsc1 ) * (  zficeth(ji,jj,jl)                                             &
-                  &                           + ph_snw(ji,jj,jl) * ( rn_alphd - zficeth(ji,jj,jl) ) / c1  )   &
-                  &     +         zihsc1   * rn_alphd  
-               !                                        !  melting snow                
-               zihsc2   = MAX( zzero , SIGN( zone , ph_snw(ji,jj,jl) - c2 ) )
-               zalbpsnm = ( 1.0 - zihsc2 ) * ( rn_albice + ph_snw(ji,jj,jl) * ( rn_alphc - rn_albice ) / c2 )   &
-                  &     +         zihsc2   *   rn_alphc 
-               !
-               zitmlsn  =  MAX( zzero , SIGN( zone , pt_ice(ji,jj,jl) - rt0_snow ) )   
-               zalbpsn  =  zitmlsn * zalbpsnm + ( 1.0 - zitmlsn ) * zalbpsnf
-            
-               !  Case of ice free of snow.
-               zalbpic  = zficeth(ji,jj,jl) 
-            
-               ! albedo of the system   
-               zithsn   = 1.0 - MAX( zzero , SIGN( zone , - ph_snw(ji,jj,jl) ) )
-               pa_ice_cs(ji,jj,jl) =  zithsn * zalbpsn + ( 1.0 - zithsn ) *  zalbpic
+         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
+         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 )
+         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              &
+            &                                                                 - 0.8608 * ph_ice * ph_ice     &
+            &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice
+         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice
+         END WHERE
+         !
+         DO jl = 1, ijpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! freezing snow
+                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2
+                  !                                        !  freezing snow        
+                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) )
+                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  &
+                     &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   &
+                     &        +         zswitch   * ralb_sf  
+                  !
+                  ! melting snow
+                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2
+                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) )
+                  zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   &
+                      &     +         zswitch   *   ralb_sm 
+                  !
+                  ! snow albedo
+                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
+                  !
+                  ! Ice/snow albedo
+                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
+                  pa_ice_cs(ji,jj,jl) =  zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl)
+                  !
+               END DO
             END DO
          END DO
-      END DO
-      
-      !    Albedo of snow-ice for overcast sky.
-      !----------------------------------------------  
-      pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rn_cloud       ! Oberhuber correction
-      !
-      CALL wrk_dealloc( jpi,jpj,ijpl, zalbfz, zficeth )
+         !
+         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky
+         !
+         !              !------------------------------------------
+      CASE( 1 )         !  New parameterization (2016)
+         !              !------------------------------------------
+         !
+         ralb_im = rn_albice  ! bare puddled ice
+! compilation of values from literature
+         ralb_sf = 0.85      ! dry snow
+         ralb_sm = 0.75      ! melting snow
+         ralb_if = 0.60      ! bare frozen ice
+! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved
+!         ralb_sf = 0.85       ! dry snow
+!         ralb_sm = 0.72       ! melting snow
+!         ralb_if = 0.65       ! bare frozen ice
+! Brandt et al 2005 (East Antarctica)
+!         ralb_sf = 0.87      ! dry snow
+!         ralb_sm = 0.82      ! melting snow
+!         ralb_if = 0.54      ! bare frozen ice
+! 
+         !  Computation of ice albedo (free of snow)
+         z1_c1 = 1._wp / ( LOG(1.5_wp) - LOG(0.05_wp) ) 
+         z1_c2 = 1._wp / 0.05_wp
+         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im
+         ELSE WHERE                                              ;   zalb = ralb_if
+         END  WHERE
+         !
+         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb
+         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18_wp - zalb     ) * z1_c1 *  &
+            &                                                                     ( LOG(1.5_wp) - LOG(ph_ice) )
+         ELSE WHERE                                       ;  zalb_it = ralb_oce + ( 0.18_wp - ralb_oce ) * z1_c2 * ph_ice
+         END WHERE
+         !
+         z1_c1 = 1._wp / 0.02_wp
+         z1_c2 = 1._wp / 0.03_wp
+         !  Computation of the snow/ice albedo
+         DO jl = 1, ijpl
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zalb_sf = ralb_sf - ( ralb_sf - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 );
+                  zalb_sm = ralb_sm - ( ralb_sm - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 );
+                  !
+                  ! snow albedo
+                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )   
+                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf
+                  !
+                  ! Ice/snow albedo   
+                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) )
+                  pa_ice_os(ji,jj,jl) = zswitch *  zalb_it(ji,jj,jl) + ( 1._wp - zswitch ) * zalb_st
+
+              END DO
+            END DO
+         END DO
+         ! Effect of the clouds (2d order polynomial)
+         pa_ice_cs = pa_ice_os - ( - 0.1010_wp * pa_ice_os * pa_ice_os + 0.1933_wp * pa_ice_os - 0.0148_wp ) 
+!!gm better coding ?
+!         pa_ice_cs = 0.1010_wp * pa_ice_os + 0.8067_wp ) * pa_ice_os + 0.0148_wp         
+! and even better: merge this implicit do loop with the previous one
+!!gm end
+      END SELECT
+      !
+      CALL wrk_dealloc( jpi,jpj,ijpl,   zalb, zalb_it )
       !
    END SUBROUTINE albedo_ice
@@ -181 +231 @@
       REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky
       !!
-      REAL(wp) ::   zcoef   ! local scalar
-      !!----------------------------------------------------------------------
-      !
-      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )      ! Parameterization of Briegled and Ramanathan, 1982 
-      pa_oce_cs(:,:) = zcoef               
-      pa_oce_os(:,:)  = 0.06                         ! Parameterization of Kondratyev, 1969 and Payne, 1972
+      REAL(wp) :: zcoef 
+      !!----------------------------------------------------------------------
+      !
+      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982
+      pa_oce_cs(:,:) = zcoef 
+      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972
       !
    END SUBROUTINE albedo_oce
@@ -200 +250 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namsbc_alb/ rn_cloud, rn_albice, rn_alphd, rn_alphdi, rn_alphc
+      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice 
       !!----------------------------------------------------------------------
       !
@@ -219 +269 @@
          WRITE(numout,*) '~~~~~~~'
          WRITE(numout,*) '   Namelist namsbc_alb : albedo '
-         WRITE(numout,*) '      correction for snow and ice albedo                  rn_cloud  = ', rn_cloud
-         WRITE(numout,*) '      albedo of melting ice in the arctic and antarctic   rn_albice = ', rn_albice
-         WRITE(numout,*) '      coefficients for linear                             rn_alphd  = ', rn_alphd
-         WRITE(numout,*) '      interpolation used to compute albedo                rn_alphdi = ', rn_alphdi
-         WRITE(numout,*) '      between two extremes values (Pyane, 1972)           rn_alphc  = ', rn_alphc
+         WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb
+         WRITE(numout,*) '      albedo of bare puddled ice                          rn_albice  = ', rn_albice
       ENDIF
       !

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

@@ -11 +11 @@
    !!            3.2  !  2009-04  (G. Madec & NEMO team) 
    !!            3.6  !  2012-05  (C. Rousset) store attenuation coef for use in ice model 
-   !!            3.7  !  2015-11  (G. Madec, A. Coward)  remove optimisation for fix volume 
+   !!             -   !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll
+   !!            3.7  !  2016-01  (G. Madec, A. Coward)  remove optimisation for fix volume 
    !!----------------------------------------------------------------------
 
@@ -55 +56 @@
    INTEGER , PUBLIC ::   nksr         !: levels below which the light cannot penetrate (depth larger than 391 m)
  
-   INTEGER, PARAMETER ::   np_RGB  = 1   ! R-G-B     light penetration with constant Chlorophyll
-   INTEGER, PARAMETER ::   np_RGBc = 2   ! R-G-B     light penetration with Chlorophyll data
-   INTEGER, PARAMETER ::   np_2BD  = 3   ! 2 bands   light penetration
-   INTEGER, PARAMETER ::   np_BIO  = 4   ! bio-model light penetration
+   INTEGER, PARAMETER ::   np_RGB    = 1   ! R-G-B     light penetration with constant Chlorophyll
+   INTEGER, PARAMETER ::   np_2BD    = 2   ! 2 bands   light penetration
+   INTEGER, PARAMETER ::   np_BIO    = 3   ! bio-model light penetration
    !
    INTEGER  ::   nqsr    ! user choice of the type of light penetration
    REAL(wp) ::   xsi0r   ! inverse of rn_si0
    REAL(wp) ::   xsi1r   ! inverse of rn_si1
+   
+   REAL(wp) ::   rChl_0 = 0.05_wp   ! value of Chlorophyll used in case of constant Chlorophyll
    !
    REAL(wp) , DIMENSION(3,61)           ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read)
+   TYPE(FLD), DIMENSION(:), ALLOCATABLE ::   sf_chl   ! structure of input Chl (file informations, fields read)
 
    !! * Substitutions
@@ -100 +102 @@
       !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp.
       !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516.
+      !!              Morel, A. and Berthon, J.F., 1989, Limnol Oceanogr 34(8), 1545-1562
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! ocean time-step
@@ -107 +110 @@
       REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars
       REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         -
-      REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         -
       REAL(wp) ::   zz0 , zz1                !    -         -
+      REAL(wp) ::   zCb, zCmax, zze, z1_ze, zpsi, zpsimax, zdelpsi, zCtot, zCze
+      REAL(wp) ::   zlogc, zlogc2, zlogc3 
       REAL(wp), POINTER, DIMENSION(:,:)   :: zekb, zekg, zekr
-      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt
+      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, zchl3d, ztrdt
       REAL(wp), POINTER, DIMENSION(:,:,:) :: zetot
       !!----------------------------------------------------------------------
@@ -149 +153 @@
       !                         !--------------------------------!
       !
-      CASE( np_BIO )                   !==  bio-model fluxes  ==!
+      CASE( np_BIO )                !==  bio-model fluxes  ==!
          !
          DO jk = 1, nksr
@@ -155 +159 @@
          END DO
          !
-      CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==!
-         !
-         CALL wrk_alloc( jpi,jpj,       zekb, zekg, zekr        ) 
-         CALL wrk_alloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea ) 
-         !
-         IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll
+      CASE( np_RGB )                !==  R-G-B fluxes  ==!
+         !
+         CALL wrk_alloc( jpi,jpj,       zekb, zekg, zekr                ) 
+         CALL wrk_alloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d ) 
+         !
+         SELECT CASE( nn_chldta )         ! set 3D chlorophyll field
+         !
+         CASE( 0 )                           ! constant 
+            DO jk = 1, nksr + 1
+               zchl3d(:,:,jk) = rChl_0
+            END DO
+            !
+         CASE( 1 )                           ! surface chlorophyl data spread uniformly on the vertical
             CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
-            DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl
+            DO jk = 1, nksr + 1                    ! uniform vertical profile
+               zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1) 
+            END DO
+            !
+         CASE( 2 )                           ! surface chlorophyl data + Morel and Berthon (1989) profile
+            CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step
+            DO jj = 2, jpjm1                       ! Chl profile = F( surface Chl value)
                DO ji = fs_2, fs_jpim1
-                  zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) )
-                  irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 )
-                  zekb(ji,jj) = rkrgb(1,irgb)
-                  zekg(ji,jj) = rkrgb(2,irgb)
-                  zekr(ji,jj) = rkrgb(3,irgb)
+                  zchl    = sf_chl(1)%fnow(ji,jj,1)
+                  zCtot   =  40.6_wp *  zchl**0.459
+                  zze     = 568.2_wp * zCtot**(-0.746)
+                  IF( zze > 102. )   zze = 200.0 * zCtot**(-0.293)
+                  zlogc   = LOG( zchl )
+!!gm : instead of this :
+                  zlogc2  = zlogc * zlogc
+                  zlogc3  = zlogc * zlogc * zlogc
+                  zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3
+                  zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2
+                  zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3
+                  zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2
+!!gm faster & more precise:
+!                  zCb     = 0.768 + zlogc * (  (  0.087 + zlogc * (- 0.179 - zlogc * 0.025 )  )
+!                  zCmax   = 0.299 + zlogc * (   - 0.289 + zlogc *    0.579  )
+!                  zpsimax = 0.6   + zlogc * (  (- 0.640 + zlogc * (  0.021 + zlogc * 0.115 )  )
+!                  zdelpsi = 0.710 + zlogc * (     0.159 + zlogc *    0.021  )
+!!gm end          
+                  zCze    = 1.12_wp * (zchl)**0.803 
+                  z1_ze   = 1._wp / zze
+                  DO jk = 1, nksr + 1
+                     zpsi = gdept_n(ji,jj,jk) * z1_ze
+                     zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) )
+                  END DO
                END DO
             END DO
-         ELSE                                !* constant chrlorophyll
-            zchl = 0.05                            ! constant chlorophyll
-            !                                      ! Separation in R-G-B depending of the chlorophyll
-            irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 )
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zekb(ji,jj) = rkrgb(1,irgb)                      
-                  zekg(ji,jj) = rkrgb(2,irgb)
-                  zekr(ji,jj) = rkrgb(3,irgb)
-               END DO
-            END DO
-         ENDIF
+            !
+         END SELECT
          !
          zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B
@@ -195 +221 @@
          END DO
          !
-         DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B
+         DO jk = 2, nksr+1                   !* interior partition in R-G-B function of 3D Chl
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                  zchl = MIN( 10._wp , MAX( 0.03_wp , zchl3d(ji,jj,jk) ) )
+                  irgb = NINT( 41._wp + 20._wp * LOG10(zchl) + 1.e-15 )
+                  zekb(ji,jj) = rkrgb(1,irgb)
+                  zekg(ji,jj) = rkrgb(2,irgb)
+                  zekr(ji,jj) = rkrgb(3,irgb)
+               END DO
+            END DO
             DO jj = 2, jpjm1
                DO ji = fs_2, fs_jpim1
@@ -219 +254 @@
          END DO
          !
-         CALL wrk_dealloc( jpi,jpj,        zekb, zekg, zekr        ) 
-         CALL wrk_dealloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea ) 
+         CALL wrk_dealloc( jpi,jpj,       zekb, zekg, zekr                ) 
+         CALL wrk_dealloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d ) 
          !
       CASE( np_2BD  )            !==  2-bands fluxes  ==!
@@ -309 +344 @@
       INTEGER  ::   ji, jj, jk                  ! dummy loop indices
       INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer
-      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars
-      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      -
+!      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars
+!      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      -
       !
       CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files
@@ -321 +356 @@
       IF( nn_timing == 1 )   CALL timing_start('tra_qsr_init')
       !
-      REWIND( numnam_ref )              ! Namelist namtra_qsr in reference     namelist
+      REWIND( numnam_ref )       ! Namelist namtra_qsr in reference     namelist
       READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp )
       !
-      REWIND( numnam_cfg )              ! Namelist namtra_qsr in configuration namelist
+      REWIND( numnam_cfg )       ! Namelist namtra_qsr in configuration namelist
       READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp )
       IF(lwm) WRITE ( numond, namtra_qsr )
       !
-      IF(lwp) THEN                ! control print
+      IF(lwp) THEN               ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation'
@@ -339 +374 @@
          WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio
          WRITE(numout,*) '      light penetration for ice-model (LIM3)       ln_qsr_ice = ', ln_qsr_ice
-         WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta
+         WRITE(numout,*) '      RGB : Chl data (=1,2) or cst value (=0)      nn_chldta  = ', nn_chldta
          WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs
          WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0
@@ -346 +381 @@
       ENDIF
       !
-      ioptio = 0                    ! Parameter control
+      ioptio = 0                 ! Parameter control
       IF( ln_qsr_rgb  )   ioptio = ioptio + 1
       IF( ln_qsr_2bd  )   ioptio = ioptio + 1
@@ -355 +390 @@
       !
       IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB 
-      IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc
       IF( ln_qsr_2bd                      )   nqsr = np_2BD
       IF( ln_qsr_bio                      )   nqsr = np_BIO
       !
-      !                             ! Initialisation
+      !                          ! Initialisation
       xsi0r = 1._wp / rn_si0
       xsi1r = 1._wp / rn_si1
@@ -365 +399 @@
       SELECT CASE( nqsr )
       !                               
-      CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==!
+      CASE( np_RGB )             !==  Red-Green-Blue light penetration  ==!
          !                             
          IF(lwp)   WRITE(numout,*) '   R-G-B   light penetration '
@@ -375 +409 @@
          IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
          !
-         IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure
-            IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file'
+         SELECT CASE( nn_chldta )         ! set 3D chlorophyll field
+         CASE( 0 )                           ! constant
+            IF(lwp) WRITE(numout,*) '        constant Chlorophyll set to rChl_0 =', rChl_0 
+            !
+         CASE( 1 , 2 )                       ! 3D chlorophyl field : read 2D surface data
+            !
+            IF(lwp) WRITE(numout,*) '        surface 2D Chlorophyll field read in a file'
             ALLOCATE( sf_chl(1), STAT=ierror )
             IF( ierror > 0 ) THEN
@@ -386 +425 @@
             CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   &
                &           'Solar penetration function of read chlorophyll', 'namtra_qsr' )
-         ENDIF
-         IF( nqsr == np_RGB ) THEN                 ! constant Chl
-            IF(lwp) WRITE(numout,*) '        Constant Chlorophyll concentration = 0.05'
-         ENDIF
-         !
-      CASE( np_2BD )                   !==  2 bands light penetration  ==!
+            !
+            IF( lwp .AND. nn_chldta == 1 )   WRITE(numout,*) '        profile of chlorophyll : Chl(z) = Chl(z=0)'
+            IF( lwp .AND. nn_chldta == 2 )   WRITE(numout,*) '        profile of chlorophyll : Chl(z) = Func[Chl(z=0)]'
+            !
+         END SELECT
+         !
+      CASE( np_2BD )             !==  2 bands light penetration  ==!
          !
          IF(lwp)  WRITE(numout,*) '   2 bands light penetration'
@@ -398 +438 @@
          IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m'
          !
-      CASE( np_BIO )                   !==  BIO light penetration  ==!
+      CASE( np_BIO )                         !==  BIO light penetration  ==!
          !
          IF(lwp) WRITE(numout,*) '   bio-model light penetration'

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

@@ -174 +174 @@
                   &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  )
                ! add to the eddy viscosity coef. previously computed
+# if defined key_zdftmx_new
+               ! key_zdftmx_new: New internal wave-driven param: use avs value computed by zdftmx
+               avs (ji,jj,jk) = avs(ji,jj,jk) + zavfs + zavds
+# else
                avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds
+# endif
                avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt
                avm (ji,jj,jk) = avm(ji,jj,jk) + MAX( zavft + zavdt, zavfs + zavds )

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfmxl.F90

@@ -78 +78 @@
       !
       INTEGER  ::   ji, jj, jk   ! dummy loop indices
-      INTEGER  ::   iikn, iiki, ikt, imkt   ! local integer
+      INTEGER  ::   iikn, iiki, ikt   ! local integer
       REAL(wp) ::   zN2_c        ! local scalar
       INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace
@@ -123 +123 @@
             iiki = imld(ji,jj)
             iikn = nmln(ji,jj)
-            imkt = mikt(ji,jj)
             hmld (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth 
             hmlp (ji,jj) = gdepw_n(ji,jj,iikn  ) * ssmask(ji,jj)    ! Mixed layer depth

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftke.F90

@@ -323 +323 @@
                   zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) )
                   !                                           ! TKE Langmuir circulation source term
-                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * (1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc ) / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
+                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * (1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   &
+                     &                              / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1)
                END DO
             END DO
@@ -732 +733 @@
       !
       ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number
+# if defined key_zdftmx_new
+      ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used
+      rn_emin  = 1.e-10_wp
+      rmxl_min = 1.e-03_wp
+      IF(lwp) THEN                  ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_tke_init :  New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 '
+         WRITE(numout,*) '~~~~~~~~~~~~'
+      ENDIF
+# else
       rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity
+# endif
       !
       IF(lwp) THEN                    !* Control print

branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

@@ -541 +541 @@
    END SUBROUTINE zdf_tmx_init
 
+#elif defined key_zdftmx_new
+   !!----------------------------------------------------------------------
+   !!   'key_zdftmx_new'               Internal wave-driven vertical mixing
+   !!----------------------------------------------------------------------
+   !!   zdf_tmx       : global     momentum & tracer Kz with wave induced Kz
+   !!   zdf_tmx_init  : global     momentum & tracer Kz with wave induced Kz
+   !!----------------------------------------------------------------------
+   USE oce            ! ocean dynamics and tracers variables
+   USE dom_oce        ! ocean space and time domain variables
+   USE zdf_oce        ! ocean vertical physics variables
+   USE zdfddm         ! ocean vertical physics: double diffusive mixing
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
+   USE eosbn2         ! ocean equation of state
+   USE phycst         ! physical constants
+   USE prtctl         ! Print control
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O Manager
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! work arrays
+   USE timing         ! Timing
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   zdf_tmx         ! called in step module 
+   PUBLIC   zdf_tmx_init    ! called in nemogcm module 
+   PUBLIC   zdf_tmx_alloc   ! called in nemogcm module
+
+   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: wave-driven mixing flag
+
+   !                       !!* Namelist  namzdf_tmx : internal wave-driven mixing *
+   INTEGER  ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2)
+   LOGICAL  ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency
+   LOGICAL  ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F)
+
+   REAL(wp) ::  r1_6 = 1._wp / 6._wp
+
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ebot_tmx     ! power available from high-mode wave breaking (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   epyc_tmx     ! power available from low-mode, pycnocline-intensified wave breaking (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ecri_tmx     ! power available from low-mode, critical slope wave breaking (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbot_tmx     ! WKB decay scale for high-mode energy dissipation (m)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcri_tmx     ! decay scale for low-mode critical slope dissipation (m)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   emix_tmx     ! local energy density available for mixing (W/kg)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bflx_tmx     ! buoyancy flux Kz * N^2 (W/kg)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   pcmap_tmx    ! vertically integrated buoyancy flux (W/m2)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_ratio    ! S/T diffusivity ratio (only for ln_tsdiff=T)
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_wave     ! Internal wave-induced diffusivity
+
+   !! * Substitutions
+#  include "zdfddm_substitute.h90"
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 4.0 , NEMO Consortium (2016)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   INTEGER FUNCTION zdf_tmx_alloc()
+      !!----------------------------------------------------------------------
+      !!                ***  FUNCTION zdf_tmx_alloc  ***
+      !!----------------------------------------------------------------------
+      ALLOCATE(     ebot_tmx(jpi,jpj),  epyc_tmx(jpi,jpj),  ecri_tmx(jpi,jpj)    ,   &
+      &             hbot_tmx(jpi,jpj),  hcri_tmx(jpi,jpj),  emix_tmx(jpi,jpj,jpk),   &
+      &         bflx_tmx(jpi,jpj,jpk), pcmap_tmx(jpi,jpj), zav_ratio(jpi,jpj,jpk),   & 
+      &         zav_wave(jpi,jpj,jpk), STAT=zdf_tmx_alloc     )
+      !
+      IF( lk_mpp             )   CALL mpp_sum ( zdf_tmx_alloc )
+      IF( zdf_tmx_alloc /= 0 )   CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays')
+   END FUNCTION zdf_tmx_alloc
+
+
+   SUBROUTINE zdf_tmx( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_tmx  ***
+      !!                   
+      !! ** Purpose :   add to the vertical mixing coefficients the effect of
+      !!              breaking internal waves.
+      !!
+      !! ** Method  : - internal wave-driven vertical mixing is given by:
+      !!                  Kz_wave = min(  100 cm2/s, f(  Reb = emix_tmx /( Nu * N^2 )  )
+      !!              where emix_tmx is the 3D space distribution of the wave-breaking 
+      !!              energy and Nu the molecular kinematic viscosity.
+      !!              The function f(Reb) is linear (constant mixing efficiency)
+      !!              if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T.
+      !!
+      !!              - Compute emix_tmx, the 3D power density that allows to compute
+      !!              Reb and therefrom the wave-induced vertical diffusivity.
+      !!              This is divided into three components:
+      !!                 1. Bottom-intensified low-mode dissipation at critical slopes
+      !!                     emix_tmx(z) = ( ecri_tmx / rau0 ) * EXP( -(H-z)/hcri_tmx )
+      !!                                   / ( 1. - EXP( - H/hcri_tmx ) ) * hcri_tmx
+      !!              where hcri_tmx is the characteristic length scale of the bottom 
+      !!              intensification, ecri_tmx a map of available power, and H the ocean depth.
+      !!                 2. Pycnocline-intensified low-mode dissipation
+      !!                     emix_tmx(z) = ( epyc_tmx / rau0 ) * ( sqrt(rn2(z))^nn_zpyc )
+      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) )
+      !!              where epyc_tmx is a map of available power, and nn_zpyc
+      !!              is the chosen stratification-dependence of the internal wave
+      !!              energy dissipation.
+      !!                 3. WKB-height dependent high mode dissipation
+      !!                     emix_tmx(z) = ( ebot_tmx / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_tmx)
+      !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_tmx) * e3w(z) )
+      !!              where hbot_tmx is the characteristic length scale of the WKB bottom 
+      !!              intensification, ebot_tmx is a map of available power, and z_wkb is the
+      !!              WKB-stretched height above bottom defined as
+      !!                    z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) )
+      !!                                 / SUM( sqrt(rn2(z'))    * e3w(z')    )
+      !!
+      !!              - update the model vertical eddy viscosity and diffusivity: 
+      !!                     avt  = avt  +    av_wave
+      !!                     avm  = avm  +    av_wave
+      !!                     avmu = avmu + mi(av_wave)
+      !!                     avmv = avmv + mj(av_wave)
+      !!
+      !!              - if namelist parameter ln_tsdiff = T, account for differential mixing:
+      !!                     avs  = avt  +    av_wave * diffusivity_ratio(Reb)
+      !!
+      !! ** Action  : - Define emix_tmx used to compute internal wave-induced mixing
+      !!              - avt, avs, avm, avmu, avmv increased by internal wave-driven mixing    
+      !!
+      !! References :  de Lavergne et al. 2015, JPO; 2016, in prep.
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step 
+      !
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      REAL(wp) ::   ztpc         ! scalar workspace
+      REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact     ! Used for vertical structure
+      REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep     ! Ocean depth
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  zwkb      ! WKB-stretched height above bottom
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  zweight   ! Weight for high mode vertical distribution
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_t     ! Molecular kinematic viscosity (T grid)
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_w     ! Molecular kinematic viscosity (W grid)
+      REAL(wp), DIMENSION(:,:,:), POINTER ::  zReb      ! Turbulence intensity parameter
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('zdf_tmx')
+      !
+      CALL wrk_alloc( jpi,jpj,       zfact, zhdep )
+      CALL wrk_alloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb )
+
+      !                          ! ----------------------------- !
+      !                          !  Internal wave-driven mixing  !  (compute zav_wave)
+      !                          ! ----------------------------- !
+      !                             
+      !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth,
+      !                                                 using an exponential decay from the seafloor.
+      DO jj = 1, jpj                ! part independent of the level
+         DO ji = 1, jpi
+            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean
+            zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_tmx(ji,jj) )  )
+            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ecri_tmx(ji,jj) / zfact(ji,jj)
+         END DO
+      END DO
+
+      DO jk = 2, jpkm1              ! complete with the level-dependent part
+         emix_tmx(:,:,jk) = zfact(:,:) * (  EXP( ( fsde3w(:,:,jk  ) - zhdep(:,:) ) / hcri_tmx(:,:) )                      &
+            &                             - EXP( ( fsde3w(:,:,jk-1) - zhdep(:,:) ) / hcri_tmx(:,:) )  ) * wmask(:,:,jk)   &
+            &                          / ( fsde3w(:,:,jk) - fsde3w(:,:,jk-1) )
+      END DO
+
+      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying 
+      !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc
+
+      SELECT CASE ( nn_zpyc )
+
+      CASE ( 1 )               ! Dissipation scales as N (recommended)
+
+         zfact(:,:) = 0._wp
+         DO jk = 2, jpkm1              ! part independent of the level
+            zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
+         END DO
+
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) )
+            END DO
+         END DO
+
+         DO jk = 2, jpkm1              ! complete with the level-dependent part
+            emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
+         END DO
+
+      CASE ( 2 )               ! Dissipation scales as N^2
+
+         zfact(:,:) = 0._wp
+         DO jk = 2, jpkm1              ! part independent of the level
+            zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
+         END DO
+
+         DO jj= 1, jpj
+            DO ji = 1, jpi
+               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) )
+            END DO
+         END DO
+
+         DO jk = 2, jpkm1              ! complete with the level-dependent part
+            emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk)
+         END DO
+
+      END SELECT
+
+      !                        !* WKB-height dependent mixing: distribute energy over the time-varying 
+      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot)
+      
+      zwkb(:,:,:) = 0._wp
+      zfact(:,:) = 0._wp
+      DO jk = 2, jpkm1
+         zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk)
+         zwkb(:,:,jk) = zfact(:,:)
+      END DO
+
+      DO jk = 2, jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   &
+                                            &           * tmask(ji,jj,jk) / zfact(ji,jj)
+            END DO
+         END DO
+      END DO
+      zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1)
+
+      zweight(:,:,:) = 0._wp
+      DO jk = 2, jpkm1
+         zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_tmx(:,:) * wmask(:,:,jk)                    &
+            &   * (  EXP( -zwkb(:,:,jk) / hbot_tmx(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_tmx(:,:) )  )
+      END DO
+
+      zfact(:,:) = 0._wp
+      DO jk = 2, jpkm1              ! part independent of the level
+         zfact(:,:) = zfact(:,:) + zweight(:,:,jk)
+      END DO
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_tmx(ji,jj) / ( rau0 * zfact(ji,jj) )
+         END DO
+      END DO
+
+      DO jk = 2, jpkm1              ! complete with the level-dependent part
+         emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   &
+            &                                / ( fsde3w(:,:,jk) - fsde3w(:,:,jk-1) )
+      END DO
+
+
+      ! Calculate molecular kinematic viscosity
+      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  &
+         &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0
+      DO jk = 2, jpkm1
+         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk)
+      END DO
+
+      ! Calculate turbulence intensity parameter Reb
+      DO jk = 2, jpkm1
+         zReb(:,:,jk) = emix_tmx(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) )
+      END DO
+
+      ! Define internal wave-induced diffusivity
+      DO jk = 2, jpkm1
+         zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6
+      END DO
+
+      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the
+         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  IF( zReb(ji,jj,jk) > 480.00_wp ) THEN
+                     zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+                  ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN
+                     zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) )
+                  ENDIF
+               END DO
+            END DO
+         END DO
+      ENDIF
+
+      DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s
+         zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk)
+      END DO
+
+      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave
+         ztpc = 0._wp
+         DO jk = 2, jpkm1
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ztpc = ztpc + fse3w(ji,jj,jk) * e1e2t(ji,jj)   &
+                     &         * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj)
+               END DO
+            END DO
+         END DO
+         IF( lk_mpp )   CALL mpp_sum( ztpc )
+         ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing 
+ 
+         IF(lwp) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'zdf_tmx : Internal wave-driven mixing (tmx)'
+            WRITE(numout,*) '~~~~~~~ '
+            WRITE(numout,*)
+            WRITE(numout,*) '      Total power consumption by av_wave: ztpc =  ', ztpc * 1.e-12_wp, 'TW'
+         ENDIF
+      ENDIF
+
+      !                          ! ----------------------- !
+      !                          !   Update  mixing coefs  !                          
+      !                          ! ----------------------- !
+      !      
+      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature
+         DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  &
+                      &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   &
+                      &                 ) * wmask(ji,jj,jk)
+               END DO
+            END DO
+         END DO
+         CALL iom_put( "av_ratio", zav_ratio )
+         DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing
+            fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk)
+            avt  (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
+            avm  (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk)
+         END DO
+         !
+      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing
+         DO jk = 2, jpkm1
+            fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
+            avt  (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk)
+            avm  (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk)
+         END DO
+      ENDIF
+
+      DO jk = 2, jpkm1              !* update momentum diffusivity at wu and wv points
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1  ! vector opt.
+               avmu(ji,jj,jk) = avmu(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji+1,jj  ,jk) ) * wumask(ji,jj,jk)
+               avmv(ji,jj,jk) = avmv(ji,jj,jk) + 0.5_wp * ( zav_wave(ji,jj,jk) + zav_wave(ji  ,jj+1,jk) ) * wvmask(ji,jj,jk)
+            END DO
+         END DO
+      END DO
+      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! lateral boundary condition
+
+      !                             !* output internal wave-driven mixing coefficient
+      CALL iom_put( "av_wave", zav_wave )
+                                    !* output useful diagnostics: N^2, Kz * N^2 (bflx_tmx), 
+                                    !  vertical integral of rau0 * Kz * N^2 (pcmap_tmx), energy density (emix_tmx)
+      IF( iom_use("bflx_tmx") .OR. iom_use("pcmap_tmx") ) THEN
+         bflx_tmx(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:)
+         pcmap_tmx(:,:) = 0._wp
+         DO jk = 2, jpkm1
+            pcmap_tmx(:,:) = pcmap_tmx(:,:) + fse3w(:,:,jk) * bflx_tmx(:,:,jk) * wmask(:,:,jk)
+         END DO
+         pcmap_tmx(:,:) = rau0 * pcmap_tmx(:,:)
+         CALL iom_put( "bflx_tmx", bflx_tmx )
+         CALL iom_put( "pcmap_tmx", pcmap_tmx )
+      ENDIF
+      CALL iom_put( "bn2", rn2 )
+      CALL iom_put( "emix_tmx", emix_tmx )
+      
+      CALL wrk_dealloc( jpi,jpj,       zfact, zhdep )
+      CALL wrk_dealloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb )
+
+      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' tmx - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk)
+      !
+      IF( nn_timing == 1 )   CALL timing_stop('zdf_tmx')
+      !
+   END SUBROUTINE zdf_tmx
+
+
+   SUBROUTINE zdf_tmx_init
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE zdf_tmx_init  ***
+      !!                     
+      !! ** Purpose :   Initialization of the wave-driven vertical mixing, reading
+      !!              of input power maps and decay length scales in netcdf files.
+      !!
+      !! ** Method  : - Read the namzdf_tmx namelist and check the parameters
+      !!
+      !!              - Read the input data in NetCDF files :
+      !!              power available from high-mode wave breaking (mixing_power_bot.nc)
+      !!              power available from pycnocline-intensified wave-breaking (mixing_power_pyc.nc)
+      !!              power available from critical slope wave-breaking (mixing_power_cri.nc)
+      !!              WKB decay scale for high-mode wave-breaking (decay_scale_bot.nc)
+      !!              decay scale for critical slope wave-breaking (decay_scale_cri.nc)
+      !!
+      !! ** input   : - Namlist namzdf_tmx
+      !!              - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc,
+      !!              decay_scale_bot.nc decay_scale_cri.nc
+      !!
+      !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter
+      !!              - Define ebot_tmx, epyc_tmx, ecri_tmx, hbot_tmx, hcri_tmx
+      !!
+      !! References : de Lavergne et al. 2015, JPO; 2016, in prep.
+      !!         
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   inum         ! local integer
+      INTEGER  ::   ios
+      REAL(wp) ::   zbot, zpyc, zcri   ! local scalars
+      !!
+      NAMELIST/namzdf_tmx_new/ nn_zpyc, ln_mevar, ln_tsdiff
+      !!----------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )  CALL timing_start('zdf_tmx_init')
+      !
+      REWIND( numnam_ref )              ! Namelist namzdf_tmx in reference namelist : Wave-driven mixing
+      READ  ( numnam_ref, namzdf_tmx_new, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist', lwp )
+      !
+      REWIND( numnam_cfg )              ! Namelist namzdf_tmx in configuration namelist : Wave-driven mixing
+      READ  ( numnam_cfg, namzdf_tmx_new, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp )
+      IF(lwm) WRITE ( numond, namzdf_tmx_new )
+      !
+      IF(lwp) THEN                  ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'zdf_tmx_init : internal wave-driven mixing'
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namzdf_tmx_new : set wave-driven mixing parameters'
+         WRITE(numout,*) '      Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc
+         WRITE(numout,*) '      Variable (T) or constant (F) mixing efficiency            = ', ln_mevar
+         WRITE(numout,*) '      Differential internal wave-driven mixing (T) or not (F)   = ', ln_tsdiff
+      ENDIF
+      
+      ! The new wave-driven mixing parameterization elevates avt and avm in the interior, and
+      ! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should 
+      ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6).
+      avmb(:) = 1.4e-6_wp        ! viscous molecular value
+      avtb(:) = 1.e-10_wp        ! very small diffusive minimum (background avt is specified in zdf_tmx)    
+      avtb_2d(:,:) = 1.e0_wp     ! uniform 
+      IF(lwp) THEN                  ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) '   Force the background value applied to avm & avt in TKE to be everywhere ',   &
+            &               'the viscous molecular value & a very small diffusive value, resp.'
+      ENDIF
+      
+      IF( .NOT.lk_zdfddm )   CALL ctl_stop( 'STOP', 'zdf_tmx_init_new : key_zdftmx_new requires key_zdfddm' )
+      
+      !                             ! allocate tmx arrays
+      IF( zdf_tmx_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' )
+      !
+      !                             ! read necessary fields
+      CALL iom_open('mixing_power_bot',inum)       ! energy flux for high-mode wave breaking [W/m2]
+      CALL iom_get  (inum, jpdom_data, 'field', ebot_tmx, 1 ) 
+      CALL iom_close(inum)
+      !
+      CALL iom_open('mixing_power_pyc',inum)       ! energy flux for pynocline-intensified wave breaking [W/m2]
+      CALL iom_get  (inum, jpdom_data, 'field', epyc_tmx, 1 )
+      CALL iom_close(inum)
+      !
+      CALL iom_open('mixing_power_cri',inum)       ! energy flux for critical slope wave breaking [W/m2]
+      CALL iom_get  (inum, jpdom_data, 'field', ecri_tmx, 1 )
+      CALL iom_close(inum)
+      !
+      CALL iom_open('decay_scale_bot',inum)        ! spatially variable decay scale for high-mode wave breaking [m]
+      CALL iom_get  (inum, jpdom_data, 'field', hbot_tmx, 1 )
+      CALL iom_close(inum)
+      !
+      CALL iom_open('decay_scale_cri',inum)        ! spatially variable decay scale for critical slope wave breaking [m]
+      CALL iom_get  (inum, jpdom_data, 'field', hcri_tmx, 1 )
+      CALL iom_close(inum)
+
+      ebot_tmx(:,:) = ebot_tmx(:,:) * ssmask(:,:)
+      epyc_tmx(:,:) = epyc_tmx(:,:) * ssmask(:,:)
+      ecri_tmx(:,:) = ecri_tmx(:,:) * ssmask(:,:)
+
+      ! Set once for all to zero the first and last vertical levels of appropriate variables
+      emix_tmx (:,:, 1 ) = 0._wp
+      emix_tmx (:,:,jpk) = 0._wp
+      zav_ratio(:,:, 1 ) = 0._wp
+      zav_ratio(:,:,jpk) = 0._wp
+      zav_wave (:,:, 1 ) = 0._wp
+      zav_wave (:,:,jpk) = 0._wp
+
+      zbot = glob_sum( e1e2t(:,:) * ebot_tmx(:,:) )
+      zpyc = glob_sum( e1e2t(:,:) * epyc_tmx(:,:) )
+      zcri = glob_sum( e1e2t(:,:) * ecri_tmx(:,:) )
+      IF(lwp) THEN
+         WRITE(numout,*) '      High-mode wave-breaking energy:             ', zbot * 1.e-12_wp, 'TW'
+         WRITE(numout,*) '      Pycnocline-intensifed wave-breaking energy: ', zpyc * 1.e-12_wp, 'TW'
+         WRITE(numout,*) '      Critical slope wave-breaking energy:        ', zcri * 1.e-12_wp, 'TW'
+      ENDIF
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx_init')
+      !
+   END SUBROUTINE zdf_tmx_init
+
 #else
    !!----------------------------------------------------------------------