Changeset 6347


Ignore:
Timestamp:
2016-02-24T08:56:48+01:00 (5 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:
25 edited

Legend:

Unmodified
Added
Removed
  • branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Biblio/Biblio.bib

    r6320 r6347  
    472472} 
    473473 
     474@article{bouffard_Boegman_DAO2013, 
     475   author = {D. Bouffard and L. Boegman}, 
     476   title = {A diapycnal diffusivity model for stratified environmental flows}, 
     477   volume = {61-62}, 
     478   issn = {03770265}, 
     479   url = {http://dx.doi.org/10.1016/j.dynatmoce.2013.02.002}, 
     480   doi = {10.1016/j.dynatmoce.2013.02.002}, 
     481   journal = DAO, 
     482   year = {2013}, 
     483   pages = {14--34}, 
     484} 
     485 
    474486@ARTICLE{Bougeault1989, 
    475487  author = {P. Bougeault and P. Lacarrere}, 
     
    787799  volume = {34}, 
    788800  pages = {8--13} 
     801} 
     802 
     803@article{de_lavergne_JPO2016_mixing, 
     804   author = {C. de Lavergne and G. Madec and J. Le Sommer and A. J. G. Nurser and A. C. Naveira Garabato }, 
     805   title = {On Antarctic Bottom Water consumption in the abyssal ocean}, 
     806   issn = {0022-3670}, 
     807   url = {http://dx.doi.org/10.1175/JPO-D-14-0201.1}, 
     808   doi = {10.1175/JPO-D-14-0201.1}, 
     809   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.}, 
     810   journal = {Journal of Physical Oceanography}, 
     811   year = {2016}, 
     812   volume = {46},  pages = {635-–661} 
     813} 
     814 
     815@article{de_lavergne_JPO2016_efficiency, 
     816   author = {C. de Lavergne and G. Madec and J. Le Sommer and A. J. G. Nurser and A. C. Naveira Garabato }, 
     817   title = {The impact of a variable mixing efficiency on the abyssal overturning}, 
     818   issn = {0022-3670}, 
     819   url = {http://dx.doi.org//10.1175/JPO-D-14-0259.1}, 
     820   doi = {10.1175/JPO-D-14-0259.1}, 
     821   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.}, 
     822   journal = {Journal of Physical Oceanography}, 
     823   year = {2016}, 
     824   volume = {46},  pages = {663-–681} 
    789825} 
    790826 
     
    11601196} 
    11611197 
     1198@article{goff_JGR2010, 
     1199   author = {J. A. Goff}, 
     1200   title = {Global prediction of abyssal hill root-mean-square heights from small-scale altimetric gravity variability}, 
     1201   issn = {2156-2202}, 
     1202   url = {http://dx.doi.org/10.1029/2010JB007867}, 
     1203   doi = {10.1029/2010JB007867}, 
     1204   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.}, 
     1205   volume = {115}, 
     1206   number = {B12}, 
     1207   journal = {Journal of Geophysical Research: Solid Earth}, 
     1208   year = {2010}, 
     1209   pages = {B12104}, 
     1210} 
     1211 
    11621212@ARTICLE{Goosse_al_JGR99, 
    11631213  author = {H. Goosse and E. Deleersnijder and T. Fichefet and M. England}, 
     
    12641314 
    12651315@ARTICLE{Griffies_Hallberg_MWR00, 
    1266   author = {S.M. Griffies and R.H. Hallberg}, 
    1267   title = {Biharmonic friction with a smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models}, 
     1316  author = {S.M. Griffies and R.W. Hallberg}, 
     1317  title = {Biharmonic friction with a Smagorinsky-like viscosity for use in large-scale eddy-permitting ocean models}, 
    12681318  journal = MWR, 
    12691319  year = {2000}, 
     
    15861636  volume = {12}, 
    15871637  pages = {381--389} 
     1638} 
     1639 
     1640@article{Jackson_Rehmann_JPO2014, 
     1641   author = {P. R. Jackson and C. R. Rehmann}, 
     1642   title = {Experiments on differential scalar mixing in turbulence in a sheared, stratified flow}, 
     1643   journal = JPO, 
     1644   volume = {44}, 
     1645   issn = {0022-3670}, 
     1646   url = {http://dx.doi.org/10.1175/JPO-D-14-0027.1}, 
     1647   doi = {10.1175/JPO-D-14-0027.1}, 
     1648   number = {10}, 
     1649   year = {2014}, 
     1650   pages = {2661--2680}, 
    15881651} 
    15891652 
     
    24302493} 
    24312494 
     2495@ARTICLE{Morel_Berthon_LO89, 
     2496  author = {A. Morel and J.-F. Berthon}, 
     2497  title = {Surface pigments, algal biomass profiles, and potential production of the euphotic layer:  
     2498           Relationships reinvestigated in view of remote-sensing applications}, 
     2499  journal = {Limnol. Oceanogr.}, 
     2500  year = {1989}, 
     2501  volume = {34(8)}, 
     2502  pages = {1545--1562} 
     2503} 
     2504 
    24322505@ARTICLE{Morel_Maritorena_JGR01, 
    24332506  author = {A. Morel and S. Maritorena}, 
     
    24782551  title = {Estimates of the local rate of vertical diffusion from dissipation measurements}, 
    24792552  journal = JPO, 
     2553  year = {1980}, 
    24802554  volume = {10}, 
    24812555  pages = {83--89} 
     
    27142788} 
    27152789 
     2790@ARTICLE{Rousset_GMD2015, 
     2791  author = {C. Rousset and M. Vancoppenolle and G. Madec and T. Fichefet and S. Flavoni  
     2792            and A. Barth\'{e}lemy and R. Benshila and J. Chanut and C. L\'{e}vy and S. Masson and F. Vivier }, 
     2793  title  = {The Louvain-La-Neuve sea-ice model LIM3.6: Global and regional capabilities}, 
     2794  journal= {Geoscientific Model Development}, 
     2795  year = {2015}, 
     2796  volume = {8}, pages={2991--3005}, 
     2797  doi = {10.5194/gmd-8-2991-2015}, 
     2798  url = {http://dx.doi.org/10.5194/gmd-8-2991-2015} 
     2799} 
     2800 
    27162801@ARTICLE{Sadourny1975, 
    27172802  author = {R. Sadourny}, 
     
    27942879  year = {2004}, 
    27952880  pages = {245--263}, 
     2881} 
     2882 
     2883@INBOOK{Smagorinsky_93, 
     2884  author = {Smagorinsky, J.}, 
     2885  chapter = {Some historical remarks on the use of non-linear viscosities}, 
     2886  title = {Large Eddy Simulation of Complex Engineering and Geophysical Flows}, 
     2887  pages = {3--36}, 
     2888  year = {1993}, 
     2889  publisher = {Cambridge University Press, B. Galperin and S. A. Orszag (eds.)}, 
    27962890} 
    27972891 
  • branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_DIA.tex

    r6289 r6347  
    110110even without a parallel-enabled NetCDF4 library, simply by employing only one dedicated I/O server. 
    111111 
    112 \subsection{XIOS: the IO\_SERVER} 
     112\subsection{XIOS: the I/O server} 
    113113 
    114114\subsubsection{Attached or detached mode?} 
     
    14091409 
    14101410% ------------------------------------------------------------------------------------------------------------- 
    1411 %       25 hour mean and hourly Surface, Mid and Bed  
    1412 % ------------------------------------------------------------------------------------------------------------- 
    1413 \section{25 hour mean output for tidal models } 
    1414  
    1415 A module is available to compute a crudely detided M2 signal by obtaining a 25 hour mean. 
    1416 The 25 hour mean is available for daily runs by summing up the 25 hourly instantananeous hourly values from 
    1417 midnight at the start of the day to midight at the day end. 
    1418 This diagnostic is actived with the logical  $ln\_dia25h$ 
    1419  
    1420 %------------------------------------------nam_dia25h------------------------------------------------------ 
    1421 \namdisplay{nam_dia25h} 
    1422 %---------------------------------------------------------------------------------------------------------- 
    1423  
    1424 \section{Top Middle and Bed hourly output } 
    1425  
    1426 A module is available to output the surface (top), mid water and bed diagnostics of a set of standard variables.  
    1427 This can be a useful diagnostic when hourly or sub-hourly output is required in high resolution tidal outputs. 
    1428 The tidal signal is retained but the overall data usage is cut to just three vertical levels. Also the bottom level  
    1429 is calculated for each cell. 
    1430 This diagnostic is actived with the logical  $ln\_diatmb$ 
    1431  
    1432 %------------------------------------------nam_diatmb----------------------------------------------------- 
    1433 \namdisplay{nam_diatmb} 
    1434 %---------------------------------------------------------------------------------------------------------- 
    1435  
    1436 % ------------------------------------------------------------------------------------------------------------- 
    14371411%       Sections transports 
    14381412% ------------------------------------------------------------------------------------------------------------- 
     
    14401414\label{DIA_diag_dct} 
    14411415 
     1416%------------------------------------------namdct---------------------------------------------------- 
     1417\namdisplay{namdct} 
     1418%------------------------------------------------------------------------------------------------------------- 
     1419 
    14421420A module is available to compute the transport of volume, heat and salt through sections.  
    14431421This diagnostic is actived with \key{diadct}. 
     
    14591437and the time scales over which they are averaged, as well as the level of output for debugging: 
    14601438 
    1461 %------------------------------------------namdct---------------------------------------------------- 
    1462 \namdisplay{namdct} 
    1463 %------------------------------------------------------------------------------------------------------------- 
    1464  
    14651439\np{nn\_dct}: frequency of instantaneous transports computing 
    14661440 
     
    14691443\np{nn\_debug}: debugging of the section 
    14701444 
    1471 \subsubsection{ To create a binary file containing the pathway of each section } 
    1472  
    1473 In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run}, the file \texttt{ {list\_sections.ascii\_global}} 
     1445\subsubsection{ Creating a binary file containing the pathway of each section } 
     1446 
     1447In \texttt{NEMOGCM/TOOLS/SECTIONS\_DIADCT/run}, the file \textit{ {list\_sections.ascii\_global}} 
    14741448contains a list of all the sections that are to be computed (this list of sections is based on MERSEA project metrics). 
    14751449 
     
    15831557\texttt{=/0, =/ 1000.}   &  diagonal   & eastward  & westward  & postive: eastward  \\ \hline                 
    15841558\end{tabular} 
    1585  
    1586  
    1587  
    1588 % ------------------------------------------------------------------------------------------------------------- 
    1589 %       Other Diagnostics 
    1590 % ------------------------------------------------------------------------------------------------------------- 
    1591 \section{Other Diagnostics (\key{diahth}, \key{diaar5})} 
    1592 \label{DIA_diag_others} 
    1593  
    1594  
    1595 Aside from the standard model variables, other diagnostics can be computed  
    1596 on-line. The available ready-to-add diagnostics routines can be found in directory DIA.  
    1597 Among the available diagnostics the following ones are obtained when defining  
    1598 the \key{diahth} CPP key:  
    1599  
    1600 - the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth}) 
    1601  
    1602 - the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth}) 
    1603  
    1604 - the depth of the 20\deg C isotherm (\mdl{diahth}) 
    1605  
    1606 - the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth}) 
    1607  
    1608 The poleward heat and salt transports, their advective and diffusive component, and  
    1609 the meriodional stream function can be computed on-line in \mdl{diaptr}  
    1610 \np{ln\_diaptr} to true (see the \textit{\ngn{namptr} } namelist below).   
    1611 When \np{ln\_subbas}~=~true, transports and stream function are computed  
    1612 for the Atlantic, Indian, Pacific and Indo-Pacific Oceans (defined north of 30\deg S)  
    1613 as well as for the World Ocean. The sub-basin decomposition requires an input file  
    1614 (\ifile{subbasins}) which contains three 2D mask arrays, the Indo-Pacific mask  
    1615 been deduced from the sum of the Indian and Pacific mask (Fig~\ref{Fig_mask_subasins}).  
    1616  
    1617 %------------------------------------------namptr---------------------------------------------------- 
    1618 \namdisplay{namptr}  
    1619 %------------------------------------------------------------------------------------------------------------- 
    1620 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1621 \begin{figure}[!t]     \begin{center} 
    1622 \includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_mask_subasins.pdf} 
    1623 \caption{   \label{Fig_mask_subasins} 
    1624 Decomposition of the World Ocean (here ORCA2) into sub-basin used in to compute 
    1625 the heat and salt transports as well as the meridional stream-function: Atlantic basin (red),  
    1626 Pacific basin (green), Indian basin (bleue), Indo-Pacific basin (bleue+green).  
    1627 Note that semi-enclosed seas (Red, Med and Baltic seas) as well as Hudson Bay  
    1628 are removed from the sub-basins. Note also that the Arctic Ocean has been split  
    1629 into Atlantic and Pacific basins along the North fold line.  } 
    1630 \end{center}   \end{figure}   
    1631 %>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1632  
    1633 In addition, a series of diagnostics has been added in the \mdl{diaar5}.  
    1634 They corresponds to outputs that are required for AR5 simulations  
    1635 (see Section \ref{DIA_steric} below for one of them).  
    1636 Activating those outputs requires to define the \key{diaar5} CPP key. 
    1637 \\ 
    1638 \\ 
    1639  
    1640 \section{Courant numbers} 
    1641 Courant numbers provide a theoretical indication of the model's numerical stability. The advective Courant numbers can be calculated according to 
    1642 \begin{equation} 
    1643 \label{eq:CFL} 
    1644 C_u = |u|\frac{\rdt}{e_{1u}}, \quad C_v = |v|\frac{\rdt}{e_{2v}}, \quad C_w = |w|\frac{\rdt}{e_{3w}} 
    1645 \end{equation} 
    1646 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. 
    1647  
    1648 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.  
    16491559 
    16501560 
     
    18021712the \key{diaar5} defined to be called. 
    18031713 
     1714 
     1715 
     1716% ------------------------------------------------------------------------------------------------------------- 
     1717%       Other Diagnostics 
     1718% ------------------------------------------------------------------------------------------------------------- 
     1719\section{Other Diagnostics (\key{diahth}, \key{diaar5})} 
     1720\label{DIA_diag_others} 
     1721 
     1722 
     1723Aside from the standard model variables, other diagnostics can be computed on-line.  
     1724The available ready-to-add diagnostics modules can be found in directory DIA.  
     1725 
     1726\subsection{Depth of various quantities (\mdl{diahth})} 
     1727 
     1728Among the available diagnostics the following ones are obtained when defining  
     1729the \key{diahth} CPP key:  
     1730 
     1731- the mixed layer depth (based on a density criterion \citep{de_Boyer_Montegut_al_JGR04}) (\mdl{diahth}) 
     1732 
     1733- the turbocline depth (based on a turbulent mixing coefficient criterion) (\mdl{diahth}) 
     1734 
     1735- the depth of the 20\deg C isotherm (\mdl{diahth}) 
     1736 
     1737- the depth of the thermocline (maximum of the vertical temperature gradient) (\mdl{diahth}) 
     1738 
     1739% ----------------------------------------------------------- 
     1740%     Poleward heat and salt transports 
     1741% ----------------------------------------------------------- 
     1742 
     1743\subsection{Poleward heat and salt transports (\mdl{diaptr})} 
     1744 
     1745%------------------------------------------namptr----------------------------------------- 
     1746\namdisplay{namptr}  
     1747%----------------------------------------------------------------------------------------- 
     1748 
     1749The poleward heat and salt transports, their advective and diffusive component, and  
     1750the meriodional stream function can be computed on-line in \mdl{diaptr}  
     1751\np{ln\_diaptr} to true (see the \textit{\ngn{namptr} } namelist below).   
     1752When \np{ln\_subbas}~=~true, transports and stream function are computed  
     1753for the Atlantic, Indian, Pacific and Indo-Pacific Oceans (defined north of 30\deg S)  
     1754as well as for the World Ocean. The sub-basin decomposition requires an input file  
     1755(\ifile{subbasins}) which contains three 2D mask arrays, the Indo-Pacific mask  
     1756been deduced from the sum of the Indian and Pacific mask (Fig~\ref{Fig_mask_subasins}).  
     1757 
     1758%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     1759\begin{figure}[!t]     \begin{center} 
     1760\includegraphics[width=1.0\textwidth]{./TexFiles/Figures/Fig_mask_subasins.pdf} 
     1761\caption{   \label{Fig_mask_subasins} 
     1762Decomposition of the World Ocean (here ORCA2) into sub-basin used in to compute 
     1763the heat and salt transports as well as the meridional stream-function: Atlantic basin (red),  
     1764Pacific basin (green), Indian basin (bleue), Indo-Pacific basin (bleue+green).  
     1765Note that semi-enclosed seas (Red, Med and Baltic seas) as well as Hudson Bay  
     1766are removed from the sub-basins. Note also that the Arctic Ocean has been split  
     1767into Atlantic and Pacific basins along the North fold line.  } 
     1768\end{center}   \end{figure}   
     1769%>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
     1770 
     1771 
     1772% ----------------------------------------------------------- 
     1773%       CMIP specific diagnostics  
     1774% ----------------------------------------------------------- 
     1775\subsection{CMIP specific diagnostics (\mdl{diaar5})} 
     1776 
     1777A series of diagnostics has been added in the \mdl{diaar5}.  
     1778They corresponds to outputs that are required for AR5 simulations (CMIP5) 
     1779(see also Section \ref{DIA_steric} for one of them).  
     1780Activating those outputs requires to define the \key{diaar5} CPP key. 
     1781 
     1782 
     1783% ----------------------------------------------------------- 
     1784%       25 hour mean and hourly Surface, Mid and Bed  
     1785% ----------------------------------------------------------- 
     1786\subsection{25 hour mean output for tidal models } 
     1787 
     1788%------------------------------------------nam_dia25h------------------------------------- 
     1789\namdisplay{nam_dia25h} 
     1790%----------------------------------------------------------------------------------------- 
     1791 
     1792A module is available to compute a crudely detided M2 signal by obtaining a 25 hour mean. 
     1793The 25 hour mean is available for daily runs by summing up the 25 hourly instantananeous hourly values from 
     1794midnight at the start of the day to midight at the day end. 
     1795This diagnostic is actived with the logical  $ln\_dia25h$ 
     1796 
     1797 
     1798% ----------------------------------------------------------- 
     1799%     Top Middle and Bed hourly output 
     1800% ----------------------------------------------------------- 
     1801\subsection{Top Middle and Bed hourly output } 
     1802 
     1803%------------------------------------------nam_diatmb----------------------------------------------------- 
     1804\namdisplay{nam_diatmb} 
     1805%---------------------------------------------------------------------------------------------------------- 
     1806 
     1807A module is available to output the surface (top), mid water and bed diagnostics of a set of standard variables.  
     1808This can be a useful diagnostic when hourly or sub-hourly output is required in high resolution tidal outputs. 
     1809The tidal signal is retained but the overall data usage is cut to just three vertical levels. Also the bottom level  
     1810is calculated for each cell. 
     1811This diagnostic is actived with the logical  $ln\_diatmb$ 
     1812 
     1813 
     1814 
     1815% ----------------------------------------------------------- 
     1816%     Courant numbers 
     1817% ----------------------------------------------------------- 
     1818\subsection{Courant numbers} 
     1819Courant numbers provide a theoretical indication of the model's numerical stability. The advective Courant numbers can be calculated according to 
     1820\begin{equation} 
     1821\label{eq:CFL} 
     1822C_u = |u|\frac{\rdt}{e_{1u}}, \quad C_v = |v|\frac{\rdt}{e_{2v}}, \quad C_w = |w|\frac{\rdt}{e_{3w}} 
     1823\end{equation} 
     1824in 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. 
     1825 
     1826The 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.  
     1827 
     1828 
    18041829% ================================================================ 
    18051830 
  • branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_DOM.tex

    r6320 r6347  
    486486The last choice in terms of vertical coordinate concerns the presence (or not) in the model domain  
    487487of ocean cavities beneath ice shelves. Setting \np{ln\_isfcav} to true allows to manage ocean cavities,  
    488 otherwise they are filled in. 
     488otherwise they are filled in. This option is currently only available in $z$- or $zps$-coordinate, 
     489and partial step are also applied at the ocean/ice shelf interface.  
    489490 
    490491Contrary to the horizontal grid, the vertical grid is computed in the code and no  
     
    494495\ifile{bathy\_meter} file, so that the computation of the number of wet ocean point  
    495496in each water column is by-passed}.  
    496 If \np{ln\_isfcav}~=~true, an extra file input file describing the ice shelf draft  
    497 (in meters) (\ifile{isf\_draft\_meter}) is needed. 
    498  
    499497After reading the bathymetry, the algorithm for vertical grid definition differs  
    500498between the different options: 
     
    760758as the minimum and maximum depths at which the terrain-following vertical coordinate is calculated. 
    761759 
    762 Options for stretching the coordinate are provided as examples, but care must be taken to ensure  
    763 that the vertical stretch used is appropriate for the application. 
     760Options for stretching the coordinate are provided as examples, but care must be taken  
     761to ensure that the vertical stretch used is appropriate for the application. 
    764762 
    765763The original default NEMO s-coordinate stretching is available if neither of the other options  
     
    772770\end{equation} 
    773771 
    774 where $s_{min}$ is the depth at which the s-coordinate stretching starts and  
    775 allows a z-coordinate to placed on top of the stretched coordinate,  
    776 and z is the depth (negative down from the asea surface). 
     772where $s_{min}$ is the depth at which the $s$-coordinate stretching starts and  
     773allows a $z$-coordinate to placed on top of the stretched coordinate,  
     774and $z$ is the depth (negative down from the asea surface). 
    777775 
    778776\begin{equation} 
     
    886884that do not communicate with another ocean point at the same level are eliminated. 
    887885 
    888 In case of ice shelf cavities, as for the representation of bathymetry, a 2D integer array, misfdep, is created.  
    889 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.  
    890 By default, $misfdep(:,:)=1$ and no cells are masked. 
    891 Modifications of the model bathymetry and ice shelf draft into  
     886As for the representation of bathymetry, a 2D integer array, misfdep, is created.  
     887misfdep defines the level of the first wet $t$-point. All the cells between $k=1$ and $misfdep(i,j)-1$ are masked.  
     888By default, misfdep(:,:)=1 and no cells are masked. 
     889 
     890In case of ice shelf cavities, modifications of the model bathymetry and ice shelf draft into  
    892891the cavities are performed in the \textit{zgr\_isf} routine. The compatibility between ice shelf draft and bathymetry is checked.  
    893 All the locations where the isf cavity is thinnest than \np{rn\_isfhmin} meters are grounded ($i.e.$ masked).  
    894892If 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. 
    895893If the incompatibility is too strong (need to dig more than 1 cell), the cell is masked.\\  
    896894 
    897 From the \textit{mbathy} and \textit{misfdep} array, the mask fields are defined as follows: 
     895From the \textit{mbathy} array, the mask fields are defined as follows: 
    898896\begin{align*} 
    899897tmask(i,j,k) &= \begin{cases}   \; 0&   \text{ if $k < misfdep(i,j) $ } \\ 
     
    903901vmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j+1,k)   \\ 
    904902fmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i+1,j,k)   \\ 
    905                    & \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ 
     903             &    \ \ \, * tmask(i,j,k) \ * \ tmask(i+1,j,k) \\ 
    906904wmask(i,j,k) &=         \; tmask(i,j,k) \ * \ tmask(i,j,k-1) \text{ with } wmask(i,j,1) = tmask(i,j,1)  
    907905\end{align*} 
    908906 
    909 Note, wmask is now defined. It allows, in case of ice shelves,  
    910 to deal with the top boundary (ice shelf/ocean interface) exactly in the same way as for the bottom boundary.  
     907Note that, without ice shelves cavities, masks at $t-$ and $w-$points are identical with  
     908the numerical indexing used (\S~\ref{DOM_Num_Index}). Nevertheless, $wmask$ are required  
     909with oceean cavities to deal with the top boundary (ice shelf/ocean interface)  
     910exactly in the same way as for the bottom boundary.  
    911911 
    912912The specification of closed lateral boundaries requires that at least the first and last  
     
    916916(and so too the mask arrays) (see \S~\ref{LBC_jperio}). 
    917917 
    918 %%% 
    919 \gmcomment{   \colorbox{yellow}{Add one word on tricky trick !} mbathy in further modified in zdfbfr{\ldots}.  } 
    920 %%% 
    921918 
    922919% ================================================================ 
  • branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_DYN.tex

    r6320 r6347  
    637637($e_{3w}$). 
    638638  
    639 $\bullet$ Traditional coding with adaptation for ice shelf cavities (\np{ln\_dynhpg\_isf}=true). 
    640 This scheme need the activation of ice shelf cavities (\np{ln\_isfcav}=true). 
    641  
    642639$\bullet$ Pressure Jacobian scheme (prj) (a research paper in preparation) (\np{ln\_dynhpg\_prj}=true) 
    643640 
     
    667664 
    668665$\bullet$ The ocean load is computed using the expression \eqref{Eq_dynhpg_sco} described in \ref{DYN_hpg_sco}.  
     666A treatment of the partial cell for top and bottom similar to the one described in \ref{DYN_hpg_zps} is done  
     667to reduce the residual circulation generated by the top partial cell.  
    669668 
    670669%-------------------------------------------------------------------------------------------------------------- 
  • branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_LDF.tex

    r6289 r6347  
    397397\subsubsection{Space and Time Varying Mixing Coefficients} 
    398398 
    399 There is no default specification of space and time varying mixing coefficient.  
    400 The only case available is specific to the ORCA2 and ORCA05 global ocean  
    401 configurations. It provides only a tracer  
    402 mixing coefficient for eddy induced velocity (ORCA2) or both iso-neutral and  
    403 eddy induced velocity (ORCA05) that depends on the local growth rate of  
    404 baroclinic instability. This specification is actually used when an ORCA key  
     399There are no default specifications of space and time varying mixing coefficient.  One 
     400available case is specific to the ORCA2 and ORCA05 global ocean configurations. It 
     401provides only a tracer mixing coefficient for eddy induced velocity (ORCA2) or both 
     402iso-neutral and eddy induced velocity (ORCA05) that depends on the local growth rate of 
     403baroclinic instability. This specification is actually used when an ORCA key 
    405404and both \key{traldf\_eiv} and \key{traldf\_c2d} are defined. 
     405 
     406\subsubsection{Smagorinsky viscosity (\key{dynldf\_c3d} and \key{dynldf\_smag})} 
     407 
     408The \key{dynldf\_smag} key activates a 3D, time-varying viscosity that depends on the 
     409resolved motions. Following \citep{Smagorinsky_93} the viscosity coefficient is set 
     410proportional to a local deformation rate based on the horizontal shear and tension, 
     411namely: 
     412 
     413\begin{equation} 
     414A_{m_{Smag}} = \left(\frac{{\sf CM_{Smag}}}{\pi}\right)^2L^2\vert{D}\vert 
     415\end{equation} 
     416 
     417\noindent where the deformation rate $\vert{D}\vert$ is given by  
     418 
     419\begin{equation} 
     420\vert{D}\vert=\sqrt{\left({\frac{\partial{u}} {\partial{x}}} 
     421                         -{\frac{\partial{v}} {\partial{y}}}\right)^2 
     422                 +  \left({\frac{\partial{u}} {\partial{y}}} 
     423                         +{\frac{\partial{v}} {\partial{x}}}\right)^2}  
     424\end{equation} 
     425 
     426\noindent and $L$ is the local gridscale given by: 
     427 
     428\begin{equation} 
     429L^2 = \frac{2{e_1}^2 {e_2}^2}{\left ( {e_1}^2 + {e_2}^2 \right )} 
     430\end{equation} 
     431 
     432\citep{Griffies_Hallberg_MWR00} suggest values in the range 2.2 to 4.0 of the coefficient 
     433$\sf CM_{Smag}$ for oceanic flows. This value is set via the \np{rn\_cmsmag\_1} namelist 
     434parameter. An additional parameter: \np{rn\_cmsh} is included in NEMO for experimenting 
     435with the contribution of the shear term. A value of 1.0 (the default) calculates the 
     436deformation rate as above; a value of 0.0 will discard the shear term entirely. 
     437 
     438For numerical stability, the calculated viscosity is bounded according to the following: 
     439 
     440\begin{equation} 
     441{\rm MIN}\left ({ L^2\over {8\Delta{t}}}, rn\_ahm\_m\_lap\right ) \geq A_{m_{Smag}}  
     442                                                                  \geq rn\_ahm\_0\_lap 
     443\end{equation} 
     444 
     445\noindent with both parameters for the upper and lower bounds being provided via the 
     446indicated namelist parameters. 
     447 
     448\bigskip When $ln\_dynldf\_bilap = .true.$, a biharmonic version of the Smagorinsky 
     449viscosity is also available which sets a coefficient for the biharmonic viscosity as: 
     450 
     451\begin{equation} 
     452B_{m_{Smag}} = - \left(\frac{{\sf CM_{bSmag}}}{\pi}\right)^2 {L^4\over 8}\vert{D}\vert 
     453\end{equation} 
     454 
     455\noindent which is bounded according to: 
     456 
     457\begin{equation} 
     458{\rm MAX}\left (-{ L^4\over {64\Delta{t}}}, rn\_ahm\_m\_blp\right ) \leq B_{m_{Smag}}  
     459                                                                    \leq rn\_ahm\_0\_blp 
     460\end{equation} 
     461 
     462\noindent Note the reversal of the inequalities here because NEMO requires the biharmonic 
     463coefficients as negative numbers. $\sf CM_{bSmag}$ is set via the \np{rn\_cmsmag\_2} 
     464namelist parameter and the bounding values have corresponding entries in the namelist too. 
     465 
     466\bigskip The current implementation in NEMO also allows for 3D, time-varying diffusivities 
     467to be set using the Smagorinsky approach. Users should note that this option is not 
     468recommended for many applications since diffusivities will tend to be largest near 
     469boundaries (where shears are greatest) leading to spurious upwellings 
     470(\citep{Griffies_Bk04}, chapter 18.3.4). Nevertheless the option is there for those 
     471wishing to experiment. This choice requires both \key{traldf\_c3d} and \key{traldf\_smag} 
     472and uses the \np{rn\_chsmag} (${\sf CH_{Smag}}$), \np{rn\_smsh} and \np{rn\_aht\_m} 
     473namelist parameters in an analogous way to \np{rn\_cmsmag\_1}, \np{rn\_cmsh} and 
     474\np{rn\_ahm\_m\_lap} (see above) to set the diffusion coefficient: 
     475 
     476\begin{equation} 
     477A_{h_{Smag}} = \left(\frac{{\sf CH_{Smag}}}{\pi}\right)^2L^2\vert{D}\vert 
     478\end{equation} 
     479 
     480  
     481For numerical stability, the calculated diffusivity is bounded according to the following: 
     482 
     483\begin{equation} 
     484{\rm MIN}\left ({ L^2\over {8\Delta{t}}}, rn\_aht\_m\right ) \geq A_{h_{Smag}}  
     485                                                             \geq rn\_aht\_0 
     486\end{equation} 
     487 
    406488 
    407489$\ $\newline    % force a new ligne 
  • branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_SBC.tex

    r6320 r6347  
    5151\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) ;  
    5252\item the addition of river runoffs as surface freshwater fluxes or lateral inflow (\np{ln\_rnf}~=~true) ;  
    53 \item the addition of isf melting as lateral inflow (parameterisation) or as fluxes applied at the land-ice ocean interface (\np{ln\_isf}) ;  
     53\item the addition of isf melting as lateral inflow (parameterisation) (\np{nn\_isf}~=~2 or 3 and \np{ln\_isfcav}~=~false)  
     54or as fluxes applied at the land-ice ocean interface (\np{nn\_isf}~=~1 or 4 and \np{ln\_isfcav}~=~true) ;  
    5455\item the addition of a freshwater flux adjustment in order to avoid a mean sea-level drift (\np{nn\_fwb}~=~0,~1~or~2) ;  
    5556\item the transformation of the solar radiation (if provided as daily mean) into a diurnal cycle (\np{ln\_dm2dc}~=~true) ;  
     
    128129The ocean model provides, at each time step, to the surface module (\mdl{sbcmod})  
    129130the surface currents, temperature and salinity.   
    130 These variables are averaged over \np{nf\_sbc} time-step (\ref{Tab_ssm}),  
     131These variables are averaged over \np{nn\_fsbc} time-step (\ref{Tab_ssm}),  
    131132and it is these averaged fields which are used to computes the surface fluxes  
    132 at a frequency of \np{nf\_sbc} time-step. 
     133at a frequency of \np{nn\_fsbc} time-step. 
    133134 
    134135 
     
    144145\caption{  \label{Tab_ssm}    
    145146Ocean variables provided by the ocean to the surface module (SBC).  
    146 The variable are averaged over nf{\_}sbc time step, $i.e.$ the frequency of  
    147 computation of surface fluxes.} 
     147The variable are averaged over nn{\_}fsbc time step,  
     148$i.e.$ the frequency of computation of surface fluxes.} 
    148149\end{center}   \end{table} 
    149150%-------------------------------------------------------------------------------------------------------------- 
     
    557558reanalysis and satellite data. They use an inertial dissipative method to compute  
    558559the turbulent transfer coefficients (momentum, sensible heat and evaporation)  
    559 from the 10 metre wind speed, air temperature and specific humidity. 
     560from the 10 meters wind speed, air temperature and specific humidity. 
    560561This \citet{Large_Yeager_Rep04} dataset is available through the  
    561562\href{http://nomads.gfdl.noaa.gov/nomads/forms/mom4/CORE.html}{GFDL web site}.  
     
    592593or larger than the one of the input atmospheric fields. 
    593594 
     595The \np{sn\_wndi}, \np{sn\_wndj}, \np{sn\_qsr}, \np{sn\_qlw}, \np{sn\_tair}, \np{sn\_humi}, 
     596\np{sn\_prec}, \np{sn\_snow}, \np{sn\_tdif} parameters describe the fields  
     597and the way they have to be used (spatial and temporal interpolations).  
     598 
     599\np{cn\_dir} is the directory of location of bulk files 
     600\np{ln\_taudif} is the flag to specify if we use Hight Frequency (HF) tau information (.true.) or not (.false.) 
     601\np{rn\_zqt}: is the height of humidity and temperature measurements (m) 
     602\np{rn\_zu}: is the height of wind measurements (m) 
     603 
     604Three multiplicative factors are availables :  
     605\np{rn\_pfac} and \np{rn\_efac} allows to adjust (if necessary) the global freshwater budget  
     606by increasing/reducing the precipitations (total and snow) and or evaporation, respectively. 
     607The third one,\np{rn\_vfac}, control to which extend the ice/ocean velocities are taken into account  
     608in the calculation of surface wind stress. Its range should be between zero and one,  
     609and it is recommended to set it to 0. 
     610 
    594611% ------------------------------------------------------------------------------------------------------------- 
    595612%        CLIO Bulk formulea 
     
    926943\begin{description} 
    927944\item[\np{nn\_isf}~=~1] 
    928 The ice shelf cavity is represented (\np{ln\_isfcav}~=~true needed). The fwf and heat flux are computed. Two different bulk formula are available: 
     945The ice shelf cavities are explicitly represented. The fwf and heat flux are computed. Two different bulk formula are available: 
    929946   \begin{description} 
    930947   \item[\np{nn\_isfblk}~=~1] 
     
    934951   \item[\np{nn\_isfblk}~=~2]  
    935952   The bulk formula used to compute the melt is based the one described in \citet{Jenkins1991}. 
    936         This formulation is based on a 3 equations formulation (a heat flux budget, a salt flux budget 
    937          and a linearised freezing point temperature equation). 
     953        This formulation is based on a 3 equations formulation (a heat flux budget, a salt flux budget and a linearised freezing point temperature equation). 
    938954   \end{description} 
    939955 
     
    971987 
    972988\item[\np{nn\_isf}~=~4] 
    973 The ice shelf cavity is opened (\np{ln\_isfcav}~=~true needed). However, the fwf is not computed but specified from file \np{sn\_fwfisf}).  
     989The ice shelf cavity is opened. However, the fwf is not computed but specified from file \np{sn\_fwfisf}).  
    974990The heat flux ($Q_h$) is computed as $Q_h = fwf \times L_f$.\\ 
    975991\end{description} 
     
    9841000coarse to have realistic melting or for studies where you need to control your heat and fw input.\\  
    9851001 
    986 A namelist parameters control over how many meters the heat and fw fluxes are spread.  
    987 \np{rn\_hisf\_tbl}] is the top boundary layer thickness as defined in \citet{Losch2008}.  
     1002Two namelist parameters control how the heat and fw fluxes are passed to NEMO: \np{rn\_hisf\_tbl} and \np{ln\_divisf} 
     1003\begin{description} 
     1004\item[\np{rn\_hisf\_tbl}] is the top boundary layer thickness as defined in \citet{Losch2008}.  
    9881005This parameter is only used if \np{nn\_isf}~=~1 or \np{nn\_isf}~=~4 
     1006It allows you to control over which depth you want to spread the heat and fw fluxes.  
    9891007 
    9901008If \np{rn\_hisf\_tbl} = 0.0, the fluxes are put in the top level whatever is its tickness.  
    9911009 
    992 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).\\ 
    993  
    994 The ice shelf melt is implemented as a volume flux with in the same way as for the runoff. 
     1010If \np{rn\_hisf\_tbl} $>$ 0.0, the fluxes are spread over the first \np{rn\_hisf\_tbl} m (ie over one or several cells). 
     1011 
     1012\item[\np{ln\_divisf}] is a flag to apply the fw flux as a volume flux or as a salt flux.  
     1013 
     1014\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. 
    9951015The fw addition due to the ice shelf melting is, at each relevant depth level, added to the horizontal divergence  
    9961016(\textit{hdivn}) in the subroutine \rou{sbc\_isf\_div}, called from \mdl{divcur}.  
    9971017See the runoff section \ref{SBC_rnf} for all the details about the divergence correction.  
    9981018 
    999  
    1000 \section{ Ice sheet coupling} 
    1001 \label{SBC_iscpl} 
    1002 %------------------------------------------namsbc_iscpl---------------------------------------------------- 
    1003 \namdisplay{namsbc_iscpl} 
    1004 %-------------------------------------------------------------------------------------------------------- 
    1005 Ice sheet/ocean coupling is done through file exchange at the restart step. NEMO, at each restart step,  
    1006 read the bathymetry and ice shelf draft variable in a netcdf file.  
    1007 If \np{ln\_iscpl = ~true}, the isf draft is assume to be different at each restart step  
    1008 with potentially some new wet/dry cells due to the ice sheet dynamics/thermodynamics. 
    1009 The wetting and drying scheme applied on the restart is very simple and described below for the 6 different cases: 
    1010 \begin{description} 
    1011 \item[Thin a cell down:] 
    1012    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$). 
    1013 \item[Enlarge  a cell:] 
    1014    See case "Thin a cell down" 
    1015 \item[Dry a cell:] 
    1016    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$). 
    1017 \item[Wet a cell:]  
    1018    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. 
    1019 \item[Dry a column:] 
    1020    mask, T/S, U/V are set to 0 everywhere in the column and ssh set to 0. 
    1021 \item[Wet a column:] 
    1022    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. 
     1019\np{ln\_divisf}~=~false applies the fwf and heat flux directly on the salinity and temperature tendancy. 
     1020 
     1021\item[\np{ln\_conserve}] is a flag for \np{nn\_isf}~=~1. A conservative boundary layer scheme as described in \citet{Jenkins2001}  
     1022is 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.  
     1023It is only relevant for \np{ln\_divisf}~=~false.  
     1024If \np{ln\_divisf}~=~true, \np{ln\_conserve} has to be set to false to avoid a double counting of the contribution.  
     1025  
    10231026\end{description} 
    1024 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, 
    1025  the code will be unable to fill all the new wet cells properly. The default number is set up for the MISOMIP idealised experiments.\\ 
    1026 This coupling procedure is able to take into account grounding line and calving front migration. However, it is a non-conservative processe.  
    1027 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, 
    1028  a simple conservation scheme is available with \np{ln\_hsb = ~true}. The heat/salt/vol. gain/loss is diagnose, as well as the location.  
    1029 Based on what is done on sbcrnf to prescribed a source of heat/salt/vol., the heat/salt/vol. gain/loss is removed/added, 
    1030  over a period of \np{rn\_fiscpl} time step, into the system.  
    1031 So after \np{rn\_fiscpl} time step, all the heat/salt/vol. gain/loss due to extrapolation process is canceled.\\ 
    1032  
    1033 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$. 
    10341027% 
    10351028% ================================================================ 
  • branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_STO.tex

    r6289 r6347  
    55\label{STO} 
    66 
     7Authors: P.-A. Bouttier 
     8 
    79\minitoc 
    810 
     11\newpage 
    912 
    10 \newpage 
    11 $\ $\newline    % force a new line 
     13 
     14The stochastic parametrization module aims to explicitly simulate uncertainties in the model.  
     15More particularly, \cite{Brankart_OM2013} has shown that,  
     16because of the nonlinearity of the seawater equation of state, unresolved scales represent  
     17a major source of uncertainties in the computation of the large scale horizontal density gradient  
     18(from T/S large scale fields), and that the impact of these uncertainties can be simulated  
     19by random processes representing unresolved T/S fluctuations. 
     20 
     21The stochastic formulation of the equation of state can be written as: 
     22\begin{equation} 
     23 \label{eq:eos_sto} 
     24  \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)] \} 
     25\end{equation} 
     26where $p_o(z)$ is the reference pressure depending on the depth and,  
     27$\Delta T_i$ and $\Delta S_i$ are a set of T/S perturbations defined as the scalar product  
     28of the respective local T/S gradients with random walks $\mathbf{\xi}$: 
     29\begin{equation} 
     30 \label{eq:sto_pert} 
     31 \Delta T_i = \mathbf{\xi}_i \cdot \nabla T \qquad \hbox{and} \qquad \Delta S_i = \mathbf{\xi}_i \cdot \nabla S 
     32\end{equation} 
     33$\mathbf{\xi}_i$ are produced by a first-order autoregressive processes (AR-1) with  
     34a parametrized decorrelation time scale, and horizontal and vertical standard deviations $\sigma_s$.  
     35$\mathbf{\xi}$ are uncorrelated over the horizontal and fully correlated along the vertical. 
     36 
     37 
     38\section{Stochastic processes} 
     39\label{STO_the_details} 
     40 
     41The starting point of our implementation of stochastic parameterizations 
     42in NEMO is to observe that many existing parameterizations are based 
     43on autoregressive processes, which are used as a basic source of randomness 
     44to transform a deterministic model into a probabilistic model. 
     45A generic approach is thus to add one single new module in NEMO, 
     46generating processes with appropriate statistics 
     47to simulate each kind of uncertainty in the model 
     48(see \cite{Brankart_al_GMD2015} for more details). 
     49 
     50In practice, at every model grid point, independent Gaussian autoregressive 
     51processes~$\xi^{(i)},\,i=1,\ldots,m$ are first generated 
     52using the same basic equation: 
     53 
     54\begin{equation} 
     55\label{eq:autoreg} 
     56\xi^{(i)}_{k+1} = a^{(i)} \xi^{(i)}_k + b^{(i)} w^{(i)} + c^{(i)} 
     57\end{equation} 
     58 
     59\noindent 
     60where $k$ is the index of the model timestep; and 
     61$a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are parameters defining 
     62the mean ($\mu^{(i)}$) standard deviation ($\sigma^{(i)}$) 
     63and correlation timescale ($\tau^{(i)}$) of each process: 
     64 
     65\begin{itemize} 
     66\item for order~1 processes, $w^{(i)}$ is a Gaussian white noise, 
     67with zero mean and standard deviation equal to~1, and the parameters 
     68$a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are given by: 
     69 
     70\begin{equation} 
     71\label{eq:ord1} 
     72\left\{ 
     73\begin{array}{l} 
     74a^{(i)} = \varphi \\ 
     75b^{(i)} = \sigma^{(i)} \sqrt{ 1 - \varphi^2 }  
     76 \qquad\qquad\mbox{with}\qquad\qquad 
     77\varphi = \exp \left( - 1 / \tau^{(i)} \right) \\ 
     78c^{(i)} = \mu^{(i)} \left( 1 - \varphi \right) \\ 
     79\end{array} 
     80\right. 
     81\end{equation} 
     82 
     83\item for order~$n>1$ processes, $w^{(i)}$ is an order~$n-1$ autoregressive process, 
     84with zero mean, standard deviation equal to~$\sigma^{(i)}$; correlation timescale 
     85equal to~$\tau^{(i)}$; and the parameters $a^{(i)}$, $b^{(i)}$, $c^{(i)}$ are given by: 
     86 
     87\begin{equation} 
     88\label{eq:ord2} 
     89\left\{ 
     90\begin{array}{l} 
     91a^{(i)} = \varphi \\ 
     92b^{(i)} = \frac{n-1}{2(4n-3)} \sqrt{ 1 - \varphi^2 }  
     93 \qquad\qquad\mbox{with}\qquad\qquad 
     94\varphi = \exp \left( - 1 / \tau^{(i)} \right) \\ 
     95c^{(i)} = \mu^{(i)} \left( 1 - \varphi \right) \\ 
     96\end{array} 
     97\right. 
     98\end{equation} 
     99 
     100\end{itemize} 
     101 
     102\noindent 
     103In this way, higher order processes can be easily generated recursively using  
     104the same piece of code implementing Eq.~(\ref{eq:autoreg}),  
     105and using succesively processes from order $0$ to~$n-1$ as~$w^{(i)}$. 
     106The parameters in Eq.~(\ref{eq:ord2}) are computed so that this recursive application 
     107of Eq.~(\ref{eq:autoreg}) leads to processes with the required standard deviation 
     108and correlation timescale, with the additional condition that 
     109the $n-1$ first derivatives of the autocorrelation function 
     110are equal to zero at~$t=0$, so that the resulting processes 
     111become smoother and smoother as $n$ is increased. 
     112 
     113Overall, this method provides quite a simple and generic way of generating  
     114a wide class of stochastic processes.  
     115However, this also means that new model parameters are needed to specify each of  
     116these stochastic processes. As in any parameterization of lacking physics,  
     117a very important issues then to tune these new parameters using either first principles,  
     118model simulations, or real-world observations. 
     119 
     120\section{Implementation details} 
     121\label{STO_thech_details} 
     122 
    12123%---------------------------------------namsbc-------------------------------------------------- 
    13124\namdisplay{namsto} 
    14125%-------------------------------------------------------------------------------------------------------------- 
    15 $\ $\newline    % force a new ligne 
    16126 
     127The computer code implementing stochastic parametrisations can be found in the STO directory. 
     128It involves three modules :  
     129\begin{description} 
     130\item[\mdl{stopar}] : define the Stochastic parameters and their time evolution. 
     131\item[\mdl{storng}] : a random number generator based on (and includes) the 64-bit KISS  
     132                      (Keep It Simple Stupid) random number generator distributed by George Marsaglia  
     133                      (see \href{https://groups.google.com/forum/#!searchin/comp.lang.fortran/64-bit$20KISS$20RNGs}{here}) 
     134\item[\mdl{stopts}] : stochastic parametrisation associated with the non-linearity of the equation of seawater,  
     135 implementing Eq~\ref{eq:sto_pert} and specific piece of code in the equation of state implementing Eq~\ref{eq:eos_sto}. 
     136\end{description} 
    17137 
    18 See \cite{Brankart_OM2013} and \cite{Brankart_al_GMD2015} papers for a description of the parameterization. 
     138The \mdl{stopar} module has 3 public routines to be called by the model (in our case, NEMO): 
     139 
     140The first routine (\rou{sto\_par}) is a direct implementation of Eq.~(\ref{eq:autoreg}), 
     141applied at each model grid point (in 2D or 3D),  
     142and called at each model time step ($k$) to update 
     143every autoregressive process ($i=1,\ldots,m$). 
     144This routine also includes a filtering operator, applied to $w^{(i)}$, 
     145to introduce a spatial correlation between the stochastic processes. 
     146 
     147The second routine (\rou{sto\_par\_init}) is an initialization routine mainly dedicated 
     148to the computation of parameters $a^{(i)}, b^{(i)}, c^{(i)}$ 
     149for each autoregressive process, as a function of the statistical properties 
     150required by the model user (mean, standard deviation, time correlation, 
     151order of the process,\ldots).  
     152 
     153Parameters for the processes can be specified through the following \ngn{namsto} namelist parameters: 
     154\begin{description} 
     155   \item[\np{nn\_sto\_eos}]   : number of independent random walks  
     156   \item[\np{rn\_eos\_stdxy}] : random walk horz. standard deviation (in grid points) 
     157   \item[\np{rn\_eos\_stdz}]  : random walk vert. standard deviation (in grid points) 
     158   \item[\np{rn\_eos\_tcor}]  : random walk time correlation (in timesteps) 
     159   \item[\np{nn\_eos\_ord}]   : order of autoregressive processes 
     160   \item[\np{nn\_eos\_flt}]   : passes of Laplacian filter 
     161   \item[\np{rn\_eos\_lim}]   : limitation factor (default = 3.0) 
     162\end{description} 
     163This routine also includes the initialization (seeding) of the random number generator. 
     164 
     165The third routine (\rou{sto\_rst\_write}) writes a restart file (which suffix name is  
     166given by \np{cn\_storst\_out} namelist parameter) containing the current value of  
     167all autoregressive processes to allow restarting a simulation from where it has been interrupted. 
     168This file also contains the current state of the random number generator. 
     169When \np{ln\_rststo} is set to \textit{true}), the restart file (which suffix name is  
     170given by \np{cn\_storst\_in} namelist parameter) is read by the initialization routine  
     171(\rou{sto\_par\_init}). The simulation will continue exactly as if it was not interrupted 
     172only  when \np{ln\_rstseed} is set to \textit{true}, $i.e.$ when the state of  
     173the random number generator is read in the restart file. 
  • branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_TRA.tex

    r6320 r6347  
    734734(see \S\ref{SBC_rnf} for further detail of how it acts on temperature and salinity tendencies) 
    735735 
    736 $\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt, (see \S\ref{SBC_isf} for further details  
    737 on how the ice shelf melt is computed and applied). 
     736$\bullet$ \textit{fwfisf}, the mass flux associated with ice shelf melt,  
     737(see \S\ref{SBC_isf} for further details on how the ice shelf melt is computed and applied). 
    738738 
    739739The surface boundary condition on temperature and salinity is applied as follows: 
     
    840840($i.e.$ the inverses of the extinction length scales) are tabulated over 61 nonuniform  
    841841chlorophyll classes ranging from 0.01 to 10 g.Chl/L (see the routine \rou{trc\_oce\_rgb}  
    842 in \mdl{trc\_oce} module). Three types of chlorophyll can be chosen in the RGB formulation: 
    843 (1) a constant 0.05 g.Chl/L value everywhere (\np{nn\_chdta}=0) ; (2) an observed  
    844 time varying chlorophyll (\np{nn\_chdta}=1) ; (3) simulated time varying chlorophyll 
    845 by TOP biogeochemical model (\np{ln\_qsr\_bio}=true). In the latter case, the RGB  
    846 formulation is used to calculate both the phytoplankton light limitation in PISCES  
    847 or LOBSTER and the oceanic heating rate.  
    848  
     842in \mdl{trc\_oce} module). Four types of chlorophyll can be chosen in the RGB formulation: 
     843\begin{description}  
     844\item[\np{nn\_chdta}=0]  
     845a constant 0.05 g.Chl/L value everywhere ;  
     846\item[\np{nn\_chdta}=1]   
     847an observed time varying chlorophyll deduced from satellite surface ocean color measurement  
     848spread uniformly in the vertical direction ;  
     849\item[\np{nn\_chdta}=2]   
     850same as previous case except that a vertical profile of chlorophyl is used.  
     851Following \cite{Morel_Berthon_LO89}, the profile is computed from the local surface chlorophyll value ; 
     852\item[\np{ln\_qsr\_bio}=true]   
     853simulated time varying chlorophyll by TOP biogeochemical model.  
     854In this case, the RGB formulation is used to calculate both the phytoplankton  
     855light limitation in PISCES or LOBSTER and the oceanic heating rate.  
     856\end{description}  
    849857The trend in \eqref{Eq_tra_qsr} associated with the penetration of the solar radiation  
    850858is added to the temperature trend, and the surface heat flux is modified in routine \mdl{traqsr}.  
     
    13851393                   I've changed "derivative" to "difference" and "mean" to "average"} 
    13861394 
    1387 With partial cells (\np{ln\_zps}=true) at bottom and top (\np{ln\_isfcav}=true), in general, tracers in horizontally  
    1388 adjacent cells live at different depths. Horizontal gradients of tracers are needed  
    1389 for horizontal diffusion (\mdl{traldf} module) and for the hydrostatic pressure  
    1390 gradient (\mdl{dynhpg} module) to be active. The partial cell properties  
    1391 at the top (\np{ln\_isfcav}=true) are computed in the same way as for the bottom. So, only the bottom interpolation is shown. 
     1395With partial cells (\np{ln\_zps}=true) at bottom and top (\np{ln\_isfcav}=true), in general,  
     1396tracers in horizontally adjacent cells live at different depths.  
     1397Horizontal gradients of tracers are needed for horizontal diffusion (\mdl{traldf} module)  
     1398and for the hydrostatic pressure gradient (\mdl{dynhpg} module) to be active.  
    13921399\gmcomment{STEVEN from gm : question: not sure of  what -to be active- means} 
    1393  
    13941400Before taking horizontal gradients between the tracers next to the bottom, a linear  
    13951401interpolation in the vertical is used to approximate the deeper tracer as if it actually  
     
    14671473\gmcomment{gm :   this last remark has to be done} 
    14681474%%% 
     1475 
     1476If under ice shelf seas opened (\np{ln\_isfcav}=true), the partial cell properties  
     1477at the top are computed in the same way as for the bottom. Some extra variables are,  
     1478however, computed to reduce the flow generated at the top and bottom if $z*$ coordinates activated. 
     1479The extra variables calculated and used by \S\ref{DYN_hpg_isf} are: 
     1480 
     1481$\bullet$ $\overline{T}_k^{\,i+1/2}$ as described in \eqref{Eq_zps_hde} 
     1482 
     1483$\bullet$ $\delta _{i+1/2} Z_{T_k} = \widetilde {Z}^{\,i}_{T_k}-Z^{\,i}_{T_k}$ to compute  
     1484the pressure gradient correction term used by \eqref{Eq_dynhpg_sco} in \S\ref{DYN_hpg_isf}, 
     1485 with $\widetilde {Z}_{T_k}$ the depth of the point $\widetilde {T}_{k}$ in case of $z^*$ coordinates  
     1486(this term = 0 in z-coordinates) 
  • branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/Chapters/Chap_ZDF.tex

    r6320 r6347  
    262262\end{equation} 
    263263 
    264 At the ocean surface, a non zero length scale is set through the  \np{rn\_lmin0} namelist  
     264At the ocean surface, a non zero length scale is set through the  \np{rn\_mxl0} namelist  
    265265parameter. Usually the surface scale is given by $l_o = \kappa \,z_o$  
    266266where $\kappa = 0.4$ is von Karman's constant and $z_o$ the roughness  
    267267parameter of the surface. Assuming $z_o=0.1$~m \citep{Craig_Banner_JPO94}  
    268 leads to a 0.04~m, the default value of \np{rn\_lsurf}. In the ocean interior  
     268leads to a 0.04~m, the default value of \np{rn\_mxl0}. In the ocean interior  
    269269a minimum length scale is set to recover the molecular viscosity when $\bar{e}$  
    270270reach its minimum value ($1.10^{-6}= C_k\, l_{min} \,\sqrt{\bar{e}_{min}}$ ). 
     
    295295As the surface boundary condition on TKE is prescribed through $\bar{e}_o = e_{bb} |\tau| / \rho_o$,  
    296296with $e_{bb}$ the \np{rn\_ebb} namelist parameter, setting \np{rn\_ebb}~=~67.83 corresponds  
    297 to $\alpha_{CB} = 100$. further setting  \np{ln\_lsurf} to true applies \eqref{ZDF_Lsbc}  
    298 as surface boundary condition on length scale, with $\beta$ hard coded to the Stacet's value. 
     297to $\alpha_{CB} = 100$. Further setting  \np{ln\_mxl0} to true applies \eqref{ZDF_Lsbc}  
     298as surface boundary condition on length scale, with $\beta$ hard coded to the Stacey's value. 
    299299Note that a minimal threshold of \np{rn\_emin0}$=10^{-4}~m^2.s^{-2}$ (namelist parameters)  
    300300is applied on surface $\bar{e}$ value. 
     
    852852The bottom friction represents the friction generated by the bathymetry.  
    853853The top friction represents the friction generated by the ice shelf/ocean interface.  
    854 As the friction processes at the top and bottom are represented similarly, only the bottom friction is described in detail below.\\ 
     854As the friction processes at the top and bottom are represented similarly,  
     855only the bottom friction is described in detail below. 
    855856 
    856857 
     
    926927$H = 4000$~m, the resulting friction coefficient is $r = 4\;10^{-4}$~m\;s$^{-1}$.  
    927928This is the default value used in \NEMO. It corresponds to a decay time scale  
    928 of 115~days. It can be changed by specifying \np{rn\_bfric1} (namelist parameter). 
     929of 115~days. It can be changed by specifying \np{rn\_bfri1} (namelist parameter). 
    929930 
    930931For the linear friction case the coefficients defined in the general  
     
    936937\end{split} 
    937938\end{equation} 
    938 When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfric1}.  
     939When \np{nn\_botfr}=1, the value of $r$ used is \np{rn\_bfri1}.  
    939940Setting \np{nn\_botfr}=0 is equivalent to setting $r=0$ and leads to a free-slip  
    940941bottom boundary condition. These values are assigned in \mdl{zdfbfr}.  
     
    943944in the \ifile{bfr\_coef} input NetCDF file. The mask values should vary from 0 to 1.  
    944945Locations with a non-zero mask value will have the friction coefficient increased  
    945 by $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfric1}. 
     946by $mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri1}. 
    946947 
    947948% ------------------------------------------------------------------------------------------------------------- 
     
    963964$e_b = 2.5\;10^{-3}$m$^2$\;s$^{-2}$, while the FRAM experiment \citep{Killworth1992}  
    964965uses $C_D = 1.4\;10^{-3}$ and $e_b =2.5\;\;10^{-3}$m$^2$\;s$^{-2}$.  
    965 The CME choices have been set as default values (\np{rn\_bfric2} and \np{rn\_bfeb2}  
     966The CME choices have been set as default values (\np{rn\_bfri2} and \np{rn\_bfeb2}  
    966967namelist parameters). 
    967968 
     
    978979\end{equation} 
    979980 
    980 The coefficients that control the strength of the non-linear bottom friction are  
    981 initialised as namelist parameters: $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}.  
    982 Note for applications which treat tides explicitly a low or even zero value of  
    983 \np{rn\_bfeb2} is recommended. From v3.2 onwards a local enhancement of $C_D$  
    984 is possible via an externally defined 2D mask array (\np{ln\_bfr2d}=true).  
    985 See previous section for details. 
     981The coefficients that control the strength of the non-linear bottom friction are 
     982initialised as namelist parameters: $C_D$= \np{rn\_bfri2}, and $e_b$ =\np{rn\_bfeb2}. 
     983Note for applications which treat tides explicitly a low or even zero value of 
     984\np{rn\_bfeb2} is recommended. From v3.2 onwards a local enhancement of $C_D$ is possible 
     985via an externally defined 2D mask array (\np{ln\_bfr2d}=true).  This works in the same way 
     986as for the linear bottom friction case with non-zero masked locations increased by 
     987$mask\_value$*\np{rn\_bfrien}*\np{rn\_bfri2}. 
     988 
     989% ------------------------------------------------------------------------------------------------------------- 
     990%       Bottom Friction Log-layer 
     991% ------------------------------------------------------------------------------------------------------------- 
     992\subsection{Log-layer Bottom Friction enhancement (\np{nn\_botfr} = 2, \np{ln\_loglayer} = .true.)} 
     993\label{ZDF_bfr_loglayer} 
     994 
     995In the non-linear bottom friction case, the drag coefficient, $C_D$, can be optionally 
     996enhanced using a "law of the wall" scaling. If  \np{ln\_loglayer} = .true., $C_D$ is no 
     997longer constant but is related to the thickness of the last wet layer in each column by: 
     998 
     999\begin{equation} 
     1000C_D = \left ( {\kappa \over {\rm log}\left ( 0.5e_{3t}/rn\_bfrz0 \right ) } \right )^2 
     1001\end{equation} 
     1002 
     1003\noindent where $\kappa$ is the von-Karman constant and \np{rn\_bfrz0} is a roughness 
     1004length provided via the namelist. 
     1005 
     1006For stability, the drag coefficient is bounded such that it is kept greater or equal to 
     1007the base \np{rn\_bfri2} value and it is not allowed to exceed the value of an additional 
     1008namelist parameter: \np{rn\_bfri2\_max}, i.e.: 
     1009 
     1010\begin{equation} 
     1011rn\_bfri2 \leq C_D \leq rn\_bfri2\_max 
     1012\end{equation} 
     1013 
     1014\noindent Note also that a log-layer enhancement can also be applied to the top boundary 
     1015friction if under ice-shelf cavities are in use (\np{ln\_isfcav}=.true.).  In this case, the 
     1016relevant namelist parameters are \np{rn\_tfrz0}, \np{rn\_tfri2} 
     1017and \np{rn\_tfri2\_max}. 
    9861018 
    9871019% ------------------------------------------------------------------------------------------------------------- 
     
    12671299 
    12681300% ================================================================ 
     1301% Internal wave-driven mixing 
     1302% ================================================================ 
     1303\section{Internal wave-driven mixing (\key{zdftmx\_new})} 
     1304\label{ZDF_tmx_new} 
     1305 
     1306%--------------------------------------------namzdf_tmx_new------------------------------------------ 
     1307\namdisplay{namzdf_tmx_new} 
     1308%-------------------------------------------------------------------------------------------------------------- 
     1309 
     1310The parameterization of mixing induced by breaking internal waves is a generalization  
     1311of the approach originally proposed by \citet{St_Laurent_al_GRL02}.  
     1312A three-dimensional field of internal wave energy dissipation $\epsilon(x,y,z)$ is first constructed,  
     1313and the resulting diffusivity is obtained as  
     1314\begin{equation} \label{Eq_Kwave} 
     1315A^{vT}_{wave} =  R_f \,\frac{ \epsilon }{ \rho \, N^2 } 
     1316\end{equation} 
     1317where $R_f$ is the mixing efficiency and $\epsilon$ is a specified three dimensional distribution  
     1318of the energy available for mixing. If the \np{ln\_mevar} namelist parameter is set to false,  
     1319the mixing efficiency is taken as constant and equal to 1/6 \citep{Osborn_JPO80}.  
     1320In the opposite (recommended) case, $R_f$ is instead a function of the turbulence intensity parameter  
     1321$Re_b = \frac{ \epsilon}{\nu \, N^2}$, with $\nu$ the molecular viscosity of seawater,  
     1322following the model of \cite{Bouffard_Boegman_DAO2013}  
     1323and the implementation of \cite{de_lavergne_JPO2016_efficiency}. 
     1324Note that $A^{vT}_{wave}$ is bounded by $10^{-2}\,m^2/s$, a limit that is often reached when the mixing efficiency is constant. 
     1325 
     1326In addition to the mixing efficiency, the ratio of salt to heat diffusivities can chosen to vary  
     1327as a function of $Re_b$ by setting the \np{ln\_tsdiff} parameter to true, a recommended choice).  
     1328This parameterization of differential mixing, due to \cite{Jackson_Rehmann_JPO2014},  
     1329is implemented as in \cite{de_lavergne_JPO2016_efficiency}. 
     1330 
     1331The three-dimensional distribution of the energy available for mixing, $\epsilon(i,j,k)$, is constructed  
     1332from three static maps of column-integrated internal wave energy dissipation, $E_{cri}(i,j)$,  
     1333$E_{pyc}(i,j)$, and $E_{bot}(i,j)$, combined to three corresponding vertical structures  
     1334(de Lavergne et al., in prep): 
     1335\begin{align*} 
     1336F_{cri}(i,j,k) &\propto e^{-h_{ab} / h_{cri} }\\ 
     1337F_{pyc}(i,j,k) &\propto N^{n\_p}\\ 
     1338F_{bot}(i,j,k) &\propto N^2 \, e^{- h_{wkb} / h_{bot} } 
     1339\end{align*}  
     1340In the above formula, $h_{ab}$ denotes the height above bottom,  
     1341$h_{wkb}$ denotes the WKB-stretched height above bottom, defined by 
     1342\begin{equation*} 
     1343h_{wkb} = H \, \frac{ \int_{-H}^{z} N \, dz' } { \int_{-H}^{\eta} N \, dz'  } \; , 
     1344\end{equation*} 
     1345The $n_p$ parameter (given by \np{nn\_zpyc} in \ngn{namzdf\_tmx\_new} namelist)  controls the stratification-dependence of the pycnocline-intensified dissipation.  
     1346It can take values of 1 (recommended) or 2. 
     1347Finally, the vertical structures $F_{cri}$ and $F_{bot}$ require the specification of  
     1348the decay scales $h_{cri}(i,j)$ and $h_{bot}(i,j)$, which are defined by two additional input maps.  
     1349$h_{cri}$ is related to the large-scale topography of the ocean (etopo2)  
     1350and $h_{bot}$ is a function of the energy flux $E_{bot}$, the characteristic horizontal scale of  
     1351the abyssal hill topography \citep{Goff_JGR2010} and the latitude. 
     1352 
     1353% ================================================================ 
     1354 
     1355 
     1356 
  • branches/2016/dev_r6325_SIMPLIF_1/DOC/TexFiles/clean.sh

    r707 r6347  
    1 rm -f $( ls -1 NEMO_book.* | egrep -v "(tex|pdf)" ) *mpgraph* 
     1rm -f $( ls -1 ../NEMO_book.* | egrep -v "(tex|pdf)" ) *mpgraph* 
    22rm -f Chapters/*.aux Chapters/*.log 
  • branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/CONFIG/SHARED/field_def.xml

    r6140 r6347  
    479479        <!-- avt_tide: available with key_zdftmx --> 
    480480        <field id="av_tide"      long_name="tidal vertical diffusivity"   standard_name="ocean_vertical_tracer_diffusivity_due_to_tides"   unit="m2/s" /> 
     481 
     482        <!-- fields available with key_zdftmx_new --> 
     483        <field id="av_wave"      long_name="wave-induced vertical diffusivity"     standard_name="ocean_vertical_tracer_diffusivity_due_to_internal_waves"         unit="m2/s" /> 
     484        <field id="bn2"          long_name="squared Brunt-Vaisala frequency"       standard_name="squared_brunt_vaisala_frequency"                                 unit="s-1"  /> 
     485        <field id="bflx_tmx"     long_name="wave-induced buoyancy flux"            standard_name="buoyancy_flux_due_to_internal_waves"                             unit="W/kg" /> 
     486        <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" /> 
     487        <field id="emix_tmx"     long_name="power density available for mixing"    standard_name="power_available_for_mixing_from_breaking_internal_waves"         unit="W/kg" /> 
     488        <field id="av_ratio"     long_name="S over T diffusivity ratio"            standard_name="salinity_over_temperature_diffusivity_ratio"                     unit="1"    /> 
    481489 
    482490        <!-- variables available with key_diaar5 -->    
  • branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/CONFIG/SHARED/namelist_ref

    r6152 r6347  
    1111!!              6 - Tracer           (nameos, namtra_adv, namtra_ldf, namtra_ldfeiv, namtra_dmp) 
    1212!!              7 - dynamics         (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf) 
    13 !!              8 - Verical physics  (namzdf, namzdf_ric, namzdf_tke, namzdf_ddm, namzdf_tmx) 
     13!!              8 - Verical physics  (namzdf, namzdf_ric, namzdf_tke, namzdf_ddm, namzdf_tmx, namzdf_tmx_new) 
    1414!!              9 - diagnostics      (namnc4, namtrd, namspr, namflo, namhsb, namsto) 
    1515!!             10 - miscellaneous    (nammpp, namctl) 
     
    414414   ln_qsr_2bd  = .false.   !  2 bands              light penetration 
    415415   ln_qsr_bio  = .false.   !  bio-model light penetration 
    416    nn_chldta   =      1    !  RGB : Chl data (=1) or cst value (=0) 
     416   nn_chldta   =      1    !  RGB : 2D Chl data (=1), 3D Chl data (=2) or cst value (=0) 
    417417   rn_abs      =   0.58    !  RGB & 2 bands: fraction of light (rn_si1) 
    418418   rn_si0      =   0.35    !  RGB & 2 bands: shortess depth of extinction 
     
    516516&namsbc_alb    !   albedo parameters 
    517517!----------------------------------------------------------------------- 
    518    rn_cloud    =    0.06   !  cloud correction to snow and ice albedo 
    519    rn_albice   =    0.53   !  albedo of melting ice in the arctic and antarctic 
    520    rn_alphd    =    0.80   !  coefficients for linear interpolation used to 
    521    rn_alphc    =    0.65   !  compute albedo between two extremes values 
    522    rn_alphdi   =    0.72   !  (Pyane, 1972) 
     518   nn_ice_alb  =    0   !  parameterization of ice/snow albedo 
     519                        !     0: Shine & Henderson-Sellers (JGR 1985) 
     520                        !     1: "home made" based on Brandt et al. (J. Climate 2005) 
     521                        !                         and Grenfell & Perovich (JGR 2004) 
     522   rn_albice   =  0.53  !  albedo of bare puddled ice (values from 0.49 to 0.58) 
     523                        !     0.53 (default) => if nn_ice_alb=0 
     524                        !     0.50 (default) => if nn_ice_alb=1 
    523525/ 
    524526!----------------------------------------------------------------------- 
     
    575577!!====================================================================== 
    576578!!   namlbc        lateral momentum boundary condition 
    577 !!   namobc        open boundaries parameters                           ("key_obc") 
    578579!!   namagrif      agrif nested grid ( read by child model only )       ("key_agrif") 
    579580!!   nambdy        Unstructured open boundaries                         ("key_bdy") 
     
    584585&namlbc        !   lateral momentum boundary condition 
    585586!----------------------------------------------------------------------- 
     587   rn_shlat    =    2.     !  shlat = 0  !  0 < shlat < 2  !  shlat = 2  !  2 < shlat 
    586588   !                       !  free slip  !   partial slip  !   no slip   ! strong slip 
    587    rn_shlat    =    2.     !  shlat = 0  !  0 < shlat < 2  !  shlat = 2  !  2 < shlat 
    588589   ln_vorlat   = .false.   !  consistency of vorticity boundary condition with analytical Eqs. 
    589590/ 
     
    600601&nam_tide      !   tide parameters                                      ("key_tide") 
    601602!----------------------------------------------------------------------- 
    602    ln_tide_pot   = .true.   !  use tidal potential forcing 
    603    ln_tide_ramp  = .false.  ! 
    604    rdttideramp   =    0.    ! 
    605    clname(1)     = 'DUMMY'  !  name of constituent - all tidal components must be set in namelist_cfg 
     603   ln_tide_pot   = .true.  !  use tidal potential forcing 
     604   ln_tide_ramp  = .false. ! 
     605   rdttideramp   =    0.   ! 
     606   clname(1)     = 'DUMMY' !  name of constituent - all tidal components must be set in namelist_cfg 
    606607/ 
    607608!----------------------------------------------------------------------- 
     
    943944!!    namzdf_ddm    double diffusive mixing parameterization            ("key_zdfddm") 
    944945!!    namzdf_tmx    tidal mixing parameterization                       ("key_zdftmx") 
     946!!    namzdf_tmx_new    new internal wave mixing parameterization                       ("key_zdftmx") 
    945947!!====================================================================== 
    946948! 
     
    10361038   rn_tfe_itf  = 1.        !  ITF tidal dissipation efficiency 
    10371039/ 
     1040!----------------------------------------------------------------------- 
     1041&namzdf_tmx_new !   internal wave-driven mixing parameterization        ("key_zdftmx_new" & "key_zdfddm") 
     1042!----------------------------------------------------------------------- 
     1043   nn_zpyc     = 1         !  pycnocline-intensified dissipation scales as N (=1) or N^2 (=2) 
     1044   ln_mevar    = .true.    !  variable (T) or constant (F) mixing efficiency 
     1045   ln_tsdiff   = .true.    !  account for differential T/S mixing (T) or not (F) 
     1046/ 
     1047 
    10381048 
    10391049!!====================================================================== 
     
    10831093   nn_eos_flt   = 0        ! passes of Laplacian filter 
    10841094   rn_eos_lim   = 2.0      ! limitation factor (default = 3.0) 
     1095   ! 
    10851096   ln_rststo    = .false.  ! start from mean parameter (F) or from restart file (T) 
    1086    ln_rstseed = .true.           ! read seed of RNG from restart file 
    1087    cn_storst_in  = "restart_sto" !  suffix of stochastic parameter restart file (input) 
    1088    cn_storst_out = "restart_sto" !  suffix of stochastic parameter restart file (output) 
     1097   ln_rstseed   = .true.        ! read seed of RNG from restart file 
     1098   cn_storst_in = "restart_sto" !  suffix of stochastic parameter restart file (input) 
     1099   cn_storst_out= "restart_sto" !  suffix of stochastic parameter restart file (output) 
    10891100/ 
    10901101 
  • branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/CONFIG/cfg.txt

    r6140 r6347  
    77AMM12 OPA_SRC 
    88ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    9 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 
    109ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
    1110ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 
    1211GYRE OPA_SRC 
     12ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 
  • branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r6140 r6347  
    273273                        !--------------------- 
    274274                        ! Test 1: area conservation 
    275                         zA_cons = SUM(za_i_ini(ji,jj,:)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons ) 
    276                         IF ( zconv .LT. 1.0e-6 ) THEN 
     275                        zA_cons = SUM( za_i_ini(ji,jj,:) ) 
     276                        zconv   = ABS( zat_i_ini(ji,jj) - zA_cons ) 
     277                        IF ( zconv < 1.0e-6 ) THEN 
    277278                           ztest_1 = 1 
    278279                        ELSE  
    279                           !this write is useful 
    280                           IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(ji,jj)  
     280                          ! this write is useful 
     281                          IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,   & 
     282                           &                                                    ' zat_i_ini = ',zat_i_ini(ji,jj) 
    281283                          ztest_1 = 0 
    282284                        ENDIF 
    283285 
    284286                        ! Test 2: volume conservation 
    285                         zV_cons = SUM(zv_i_ini(ji,jj,:)) 
    286                         zconv = ABS(zvt_i_ini(ji,jj) - zV_cons) 
    287  
    288                         IF( zconv .LT. 1.0e-6 ) THEN 
     287                        zV_cons = SUM( zv_i_ini(ji,jj,:) ) 
     288                        zconv   = ABS( zvt_i_ini(ji,jj) - zV_cons ) 
     289 
     290                        IF( zconv < 1.0e-6 ) THEN 
    289291                           ztest_2 = 1 
    290292                        ELSE 
    291                            !this write is useful 
    292                            IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 
    293                                                     ' zvt_i_ini = ', zvt_i_ini(ji,jj) 
     293                           ! this write is useful 
     294                           IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons,   & 
     295                              &                                                    ' zvt_i_ini = ', zvt_i_ini(ji,jj) 
    294296                           ztest_2 = 0 
    295297                        ENDIF 
     
    300302                        ELSE 
    301303                           ! this write is useful 
    302                            IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zh_i_ini(ji,jj,i_fill) = ', & 
    303                            zh_i_ini(ji,jj,i_fill), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 
     304                           IF(lwp) WRITE(numout,*) ' * TEST3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** ',  & 
     305                              &    ' zh_i_ini(ji,jj,i_fill) = ', zh_i_ini(ji,jj,i_fill), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 
    304306                           IF(lwp) WRITE(numout,*) ' ji,jj,i_fill ',ji,jj,i_fill 
    305307                           IF(lwp) WRITE(numout,*) 'zht_i_ini ',zht_i_ini(ji,jj) 
     
    529531      !!----------------------------------------------------------------------------- 
    530532      ! 
    531       REWIND( numnam_ice_ref )              ! Namelist namiceini in reference namelist : Ice initial state 
     533      REWIND( numnam_ice_ref )         ! Namelist namiceini in reference namelist : Ice initial state 
    532534      READ  ( numnam_ice_ref, namiceini, IOSTAT = ios, ERR = 901) 
    533535901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in reference namelist', lwp ) 
    534536 
    535       REWIND( numnam_ice_cfg )              ! Namelist namiceini in configuration namelist : Ice initial state 
     537      REWIND( numnam_ice_cfg )         ! Namelist namiceini in configuration namelist : Ice initial state 
    536538      READ  ( numnam_ice_cfg, namiceini, IOSTAT = ios, ERR = 902 ) 
    537539902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp ) 
    538540      IF(lwm) WRITE ( numoni, namiceini ) 
    539541 
    540       slf_i(jp_hti) = sn_hti  ;  slf_i(jp_hts) = sn_hts 
    541       slf_i(jp_ati) = sn_ati  ;  slf_i(jp_tsu) = sn_tsu 
    542       slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_smi) = sn_smi 
    543  
    544       ! Define the initial parameters 
    545       ! ------------------------- 
    546  
    547       IF(lwp) THEN 
     542      slf_i(jp_hti) = sn_hti   ;   slf_i(jp_hts) = sn_hts 
     543      slf_i(jp_ati) = sn_ati   ;   slf_i(jp_tsu) = sn_tsu 
     544      slf_i(jp_tmi) = sn_tmi   ;   slf_i(jp_smi) = sn_smi 
     545 
     546      IF(lwp) THEN                     ! control print 
    548547         WRITE(numout,*) 
    549548         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' 
     
    571570            CALL ctl_stop( 'Ice_ini in limistate: unable to allocate si structure' )   ;   RETURN 
    572571         ENDIF 
    573  
     572         ! 
    574573         DO ifpr = 1, jpfldi 
    575574            ALLOCATE( si(ifpr)%fnow(jpi,jpj,1) ) 
    576575            ALLOCATE( si(ifpr)%fdta(jpi,jpj,1,2) ) 
    577576         END DO 
    578  
     577         ! 
    579578         ! fill si with slf_i and control print 
    580579         CALL fld_fill( si, slf_i, cn_dir, 'lim_istate', 'lim istate ini', 'numnam_ice' ) 
    581  
     580         ! 
    582581         CALL fld_read( nit000, 1, si )                ! input fields provided at the current time-step 
    583  
     582         ! 
    584583      ENDIF 
    585  
     584      ! 
    586585   END SUBROUTINE lim_istate_init 
    587586 
  • branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/DIA/dia25h.F90

    r6140 r6347  
    44   !! Harmonic analysis of tidal constituents  
    55   !!====================================================================== 
    6    !! History :  3.6  !  2014  (E O'Dea)  Original code 
     6   !! History :  3.0  !  07-2004  (A. Hines) New routine, developed from dia_wri_foam 
     7   !!            3.4  !  02-2013  (J. Siddorn) Routine taken from old dia_wri_foam 
     8   !!            3.6  !  11-2014  (E O'Dea)  Original code 
    79   !!---------------------------------------------------------------------- 
    8    USE oce             ! ocean dynamics and tracers variables 
    9    USE dom_oce         ! ocean space and time domain 
    10    USE in_out_manager  ! I/O units 
    11    USE iom             ! I/0 library 
    12    USE wrk_nemo        ! working arrays 
    13 #if defined key_zdftke  
    14    USE zdf_oce, ONLY: en 
    15 #endif 
    16    USE zdf_oce, ONLY: avt, avm 
    17 #if defined key_zdfgls 
    18    USE zdf_oce, ONLY: en 
    19    USE zdfgls, ONLY: mxln 
    20 #endif 
     10   USE oce            ! ocean dynamics and tracers variables 
     11   USE dom_oce        ! ocean space and time domain 
     12   USE zdf_oce        ! ocean vertical physics variables 
     13   USE zdfgls 
     14   ! 
     15   USE in_out_manager ! I/O units 
     16   USE iom            ! I/0 library 
     17   USE wrk_nemo       ! working arrays 
    2118 
    2219   IMPLICIT NONE 
    2320   PRIVATE 
    2421 
    25    LOGICAL , PUBLIC ::   ln_dia25h     !:  25h mean output 
    26    PUBLIC   dia_25h_init               ! routine called by nemogcm.F90 
    27    PUBLIC   dia_25h                    ! routine called by diawri.F90 
     22   PUBLIC   dia_25h_init   ! routine called by nemogcm.F90 
     23   PUBLIC   dia_25h        ! routine called by diawri.F90 
     24 
     25   LOGICAL, PUBLIC ::   ln_dia25h     !:  25h mean output 
    2826 
    2927  !! * variables for calculating 25-hourly means 
     28  REAL(wp) ::   r1_25 = 1._wp / 25.0_wp    ! factor for the mean calulation 
    3029   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:,:) ::   tn_25h  , sn_25h 
    3130   REAL(wp),SAVE, ALLOCATABLE,   DIMENSION(:,:)   ::   sshn_25h  
     
    3938#endif 
    4039   INTEGER, SAVE :: cnt_25h     ! Counter for 25 hour means 
    41  
    42  
    4340 
    4441   !!---------------------------------------------------------------------- 
     
    5653      !!         
    5754      !! ** Method : Read namelist 
    58       !!   History 
    59       !!   3.6  !  08-14  (E. O'Dea) Routine to initialize dia_25h 
     55      !! 
    6056      !!--------------------------------------------------------------------------- 
    61       !! 
    6257      INTEGER ::   ios                 ! Local integer output status for namelist read 
    6358      INTEGER ::   ierror              ! Local integer for memory allocation 
     
    6661      !!---------------------------------------------------------------------- 
    6762      ! 
    68       REWIND ( numnam_ref )              ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics 
     63      REWIND ( numnam_ref )          ! Read Namelist nam_dia25h in reference namelist : 25hour mean diagnostics 
    6964      READ   ( numnam_ref, nam_dia25h, IOSTAT=ios, ERR= 901 ) 
    7065901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in reference namelist', lwp ) 
    7166 
    72       REWIND( numnam_cfg )              ! Namelist nam_dia25h in configuration namelist  25hour diagnostics 
     67      REWIND( numnam_cfg )           ! Namelist nam_dia25h in configuration namelist  25hour diagnostics 
    7368      READ  ( numnam_cfg, nam_dia25h, IOSTAT = ios, ERR = 902 ) 
    7469902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_dia25h in configuration namelist', lwp ) 
     
    134129      ! ------------------------- ! 
    135130      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)  
    136       tn_25h(:,:,:) = tsb(:,:,:,jp_tem) 
    137       sn_25h(:,:,:) = tsb(:,:,:,jp_sal) 
    138       sshn_25h(:,:) = sshb(:,:) 
    139       un_25h(:,:,:) = ub(:,:,:) 
    140       vn_25h(:,:,:) = vb(:,:,:) 
    141       wn_25h(:,:,:) = wn(:,:,:) 
    142       avt_25h(:,:,:) = avt(:,:,:) 
    143       avm_25h(:,:,:) = avm(:,:,:) 
    144 # if defined key_zdfgls || defined key_zdftke 
    145          en_25h(:,:,:) = en(:,:,:) 
    146 #endif 
    147 # if defined key_zdfgls 
    148          rmxln_25h(:,:,:) = mxln(:,:,:) 
     131      tn_25h  (:,:,:) = tsb (:,:,:,jp_tem) 
     132      sn_25h  (:,:,:) = tsb (:,:,:,jp_sal) 
     133      sshn_25h(:,:)   = sshb(:,:) 
     134      un_25h  (:,:,:) = ub  (:,:,:) 
     135      vn_25h  (:,:,:) = vb  (:,:,:) 
     136      wn_25h  (:,:,:) = wn  (:,:,:) 
     137      avt_25h (:,:,:) = avt (:,:,:) 
     138      avm_25h (:,:,:) = avm (:,:,:) 
     139# if defined key_zdfgls || defined key_zdftke 
     140      en_25h  (:,:,:) = en  (:,:,:) 
     141#endif 
     142# if defined key_zdfgls 
     143      rmxln_25h(:,:,:) = mxln(:,:,:) 
    149144#endif 
    150145#if defined key_lim3 || defined key_lim2 
    151          CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice') 
     146      CALL ctl_stop('STOP', 'dia_25h not setup yet to do tidemean ice') 
    152147#endif  
    153  
    154       ! -------------------------- ! 
    155       ! 3 - Return to dia_wri      ! 
    156       ! -------------------------- ! 
    157  
    158  
     148      ! 
    159149   END SUBROUTINE dia_25h_init 
    160150 
     
    163153      !!---------------------------------------------------------------------- 
    164154      !!                 ***  ROUTINE dia_25h  *** 
    165       !!          
    166       !! 
    167       !!-------------------------------------------------------------------- 
    168155      !!                    
    169156      !! ** Purpose :   Write diagnostics with M2/S2 tide removed 
    170157      !! 
    171       !! ** Method  :    
    172       !!      25hr mean outputs for shelf seas 
     158      !! ** Method  :   25hr mean outputs for shelf seas 
    173159      !! 
    174       !! History : 
    175       !!   ?.0  !  07-04  (A. Hines) New routine, developed from dia_wri_foam 
    176       !!   3.4  !  02-13  (J. Siddorn) Routine taken from old dia_wri_foam 
    177       !!   3.6  !  08-14  (E. O'Dea) adapted for VN3.6 
    178       !!---------------------------------------------------------------------- 
    179       !! * Modules used 
    180  
    181       IMPLICIT NONE 
    182  
    183       !! * Arguments 
     160      !!---------------------------------------------------------------------- 
    184161      INTEGER, INTENT( in ) ::   kt      ! ocean time-step index 
    185  
    186  
    187       !! * Local declarations 
     162      ! 
    188163      INTEGER ::   ji, jj, jk 
    189  
    190164      LOGICAL ::   ll_print = .FALSE.    ! =T print and flush numout 
    191       REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi       ! temporary reals 
    192       INTEGER                          ::   i_steps                               ! no of timesteps per hour 
    193       REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                    ! temporary workspace 
    194       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                                  ! temporary workspace 
    195       REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                                 ! temporary workspace 
    196       INTEGER                          ::   iyear0, nimonth0,iday0                ! start year,imonth,day 
    197  
     165      REAL(wp)                         ::   zsto, zout, zmax, zjulian, zmdi   ! temporary reals 
     166      INTEGER                          ::   i_steps                           ! no of timesteps per hour 
     167      REAL(wp), DIMENSION(jpi,jpj    ) ::   zw2d, un_dm, vn_dm                ! temporary workspace 
     168      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zw3d                              ! temporary workspace 
     169      REAL(wp), DIMENSION(jpi,jpj,3)   ::   zwtmb                             ! temporary workspace 
     170      INTEGER                          ::   iyear0, nimonth0,iday0            ! start year,imonth,day 
    198171      !!---------------------------------------------------------------------- 
    199172 
     
    216189      ! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours 
    217190      ! every day 
    218       IF( MOD( kt, i_steps ) == 0  .and. kt .ne. nn_it000 ) THEN 
    219  
    220          IF (lwp) THEN 
    221               WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ',kt 
    222               WRITE(numout,*) '~~~~~~~~~~~~ ' 
     191      IF( MOD( kt, i_steps ) == 0  .AND. kt /= nn_it000 ) THEN 
     192         ! 
     193         IF(lwp) THEN 
     194            WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ', kt 
     195            WRITE(numout,*) '~~~~~~~~~~~~ ' 
    223196         ENDIF 
    224  
    225          tn_25h(:,:,:)        = tn_25h(:,:,:) + tsn(:,:,:,jp_tem) 
    226          sn_25h(:,:,:)        = sn_25h(:,:,:) + tsn(:,:,:,jp_sal) 
    227          sshn_25h(:,:)        = sshn_25h(:,:) + sshn (:,:) 
    228          un_25h(:,:,:)        = un_25h(:,:,:) + un(:,:,:) 
    229          vn_25h(:,:,:)        = vn_25h(:,:,:) + vn(:,:,:) 
    230          wn_25h(:,:,:)        = wn_25h(:,:,:) + wn(:,:,:) 
    231          avt_25h(:,:,:)       = avt_25h(:,:,:) + avt(:,:,:) 
    232          avm_25h(:,:,:)       = avm_25h(:,:,:) + avm(:,:,:) 
    233 # if defined key_zdfgls || defined key_zdftke 
    234          en_25h(:,:,:)        = en_25h(:,:,:) + en(:,:,:) 
    235 #endif 
    236 # if defined key_zdfgls 
    237          rmxln_25h(:,:,:)      = rmxln_25h(:,:,:) + mxln(:,:,:) 
     197         ! 
     198         tn_25h  (:,:,:) = tn_25h  (:,:,:) + tsn (:,:,:,jp_tem) 
     199         sn_25h  (:,:,:) = sn_25h  (:,:,:) + tsn (:,:,:,jp_sal) 
     200         sshn_25h(:,:)   = sshn_25h(:,:)   + sshn(:,:) 
     201         un_25h  (:,:,:) = un_25h  (:,:,:) + un  (:,:,:) 
     202         vn_25h  (:,:,:) = vn_25h  (:,:,:) + vn  (:,:,:) 
     203         wn_25h  (:,:,:) = wn_25h  (:,:,:) + wn  (:,:,:) 
     204         avt_25h (:,:,:) = avt_25h (:,:,:) + avt (:,:,:) 
     205         avm_25h (:,:,:) = avm_25h (:,:,:) + avm (:,:,:) 
     206# if defined key_zdfgls || defined key_zdftke 
     207         en_25h  (:,:,:) = en_25h  (:,:,:) + en  (:,:,:) 
     208#endif 
     209# if defined key_zdfgls 
     210         rmxln_25h(:,:,:) = rmxln_25h(:,:,:) + mxln(:,:,:) 
    238211#endif 
    239212         cnt_25h = cnt_25h + 1 
    240  
    241          IF (lwp) THEN 
    242             WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h 
     213         ! 
     214         IF(lwp) THEN 
     215            WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far', cnt_25h 
    243216         ENDIF 
    244  
     217         ! 
    245218      ENDIF ! MOD( kt, i_steps ) == 0 
    246  
    247          ! Write data for 25 hour mean output streams 
    248       IF( cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN 
    249  
    250             IF(lwp) THEN 
    251                WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt 
    252                WRITE(numout,*) '~~~~~~~~~~~~ ' 
    253             ENDIF 
    254  
    255             tn_25h(:,:,:)        = tn_25h(:,:,:) / 25.0_wp 
    256             sn_25h(:,:,:)        = sn_25h(:,:,:) / 25.0_wp 
    257             sshn_25h(:,:)        = sshn_25h(:,:) / 25.0_wp 
    258             un_25h(:,:,:)        = un_25h(:,:,:) / 25.0_wp 
    259             vn_25h(:,:,:)        = vn_25h(:,:,:) / 25.0_wp 
    260             wn_25h(:,:,:)        = wn_25h(:,:,:) / 25.0_wp 
    261             avt_25h(:,:,:)       = avt_25h(:,:,:) / 25.0_wp 
    262             avm_25h(:,:,:)       = avm_25h(:,:,:) / 25.0_wp 
    263 # if defined key_zdfgls || defined key_zdftke 
    264             en_25h(:,:,:)        = en_25h(:,:,:) / 25.0_wp 
    265 #endif 
    266 # if defined key_zdfgls 
    267             rmxln_25h(:,:,:)       = rmxln_25h(:,:,:) / 25.0_wp 
    268 #endif 
    269  
    270             IF (lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output' 
    271             zmdi=1.e+20 !missing data indicator for masking 
    272             ! write tracers (instantaneous) 
    273             zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    274             CALL iom_put("temper25h", zw3d)   ! potential temperature 
    275             zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    276             CALL iom_put( "salin25h", zw3d  )   ! salinity 
    277             zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
    278             CALL iom_put( "ssh25h", zw2d )   ! sea surface  
    279  
    280  
    281             ! Write velocities (instantaneous) 
    282             zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:)) 
    283             CALL iom_put("vozocrtx25h", zw3d)    ! i-current 
    284             zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:)) 
    285             CALL iom_put("vomecrty25h", zw3d  )   ! j-current 
    286  
    287             zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    288             CALL iom_put("vomecrtz25h", zw3d )   ! k-current 
    289             zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    290             CALL iom_put("avt25h", zw3d )   ! diffusivity 
    291             zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    292             CALL iom_put("avm25h", zw3d)   ! viscosity 
     219      ! 
     220      ! Write data for 25 hour mean output streams 
     221      IF( cnt_25h == 25 .AND. MOD( kt, i_steps*24) == 0 .AND. kt /= nn_it000 ) THEN 
     222         ! 
     223         IF(lwp) THEN 
     224            WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt 
     225            WRITE(numout,*) '~~~~~~~~~~~~ ' 
     226         ENDIF 
     227         tn_25h  (:,:,:) = tn_25h  (:,:,:) * r1_25 
     228         sn_25h  (:,:,:) = sn_25h  (:,:,:) * r1_25 
     229         sshn_25h(:,:)   = sshn_25h(:,:)   * r1_25 
     230         un_25h  (:,:,:) = un_25h  (:,:,:) * r1_25 
     231         vn_25h  (:,:,:) = vn_25h  (:,:,:) * r1_25 
     232         wn_25h  (:,:,:) = wn_25h  (:,:,:) * r1_25 
     233         avt_25h (:,:,:) = avt_25h (:,:,:) * r1_25 
     234         avm_25h (:,:,:) = avm_25h (:,:,:) * r1_25 
     235# if defined key_zdfgls || defined key_zdftke 
     236         en_25h  (:,:,:)        = en_25h(:,:,:) * r1_25 
     237#endif 
     238# if defined key_zdfgls 
     239         rmxln_25h(:,:,:)       = rmxln_25h(:,:,:) * r1_25 
     240#endif 
     241 
     242         IF(lwp)  WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output' 
     243         zmdi=1.e+20 !missing data indicator for masking 
     244         ! write tracers (instantaneous) 
     245         zw3d(:,:,:) = tn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     246         CALL iom_put("temper25h", zw3d)   ! potential temperature 
     247         zw3d(:,:,:) = sn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     248         CALL iom_put( "salin25h", zw3d  )   ! salinity 
     249         zw2d(:,:) = sshn_25h(:,:)*tmask(:,:,1) + zmdi*(1.0-tmask(:,:,1)) 
     250         CALL iom_put( "ssh25h", zw2d )   ! sea surface  
     251         ! 
     252         ! Write velocities (instantaneous) 
     253         zw3d(:,:,:) = un_25h(:,:,:)*umask(:,:,:) + zmdi*(1.0-umask(:,:,:)) 
     254         CALL iom_put("vozocrtx25h", zw3d)    ! i-current 
     255         zw3d(:,:,:) = vn_25h(:,:,:)*vmask(:,:,:) + zmdi*(1.0-vmask(:,:,:)) 
     256         CALL iom_put("vomecrty25h", zw3d  )   ! j-current 
     257         zw3d(:,:,:) = wn_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     258         CALL iom_put("vomecrtz25h", zw3d )   ! k-current 
     259         zw3d(:,:,:) = avt_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     260         CALL iom_put("avt25h", zw3d )   ! diffusivity 
     261         zw3d(:,:,:) = avm_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     262         CALL iom_put("avm25h", zw3d)   ! viscosity 
    293263#if defined key_zdftke || defined key_zdfgls  
    294             zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    295             CALL iom_put("tke25h", zw3d)   ! tke 
     264         zw3d(:,:,:) = en_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     265         CALL iom_put("tke25h", zw3d)   ! tke 
    296266#endif 
    297267#if defined key_zdfgls  
    298             zw3d(:,:,:) = rmxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
    299             CALL iom_put( "mxln25h",zw3d) 
    300 #endif 
    301  
    302             ! After the write reset the values to cnt=1 and sum values equal current value  
    303             tn_25h(:,:,:) = tsn(:,:,:,jp_tem) 
    304             sn_25h(:,:,:) = tsn(:,:,:,jp_sal) 
    305             sshn_25h(:,:) = sshn (:,:) 
    306             un_25h(:,:,:) = un(:,:,:) 
    307             vn_25h(:,:,:) = vn(:,:,:) 
    308             wn_25h(:,:,:) = wn(:,:,:) 
    309             avt_25h(:,:,:) = avt(:,:,:) 
    310             avm_25h(:,:,:) = avm(:,:,:) 
    311 # if defined key_zdfgls || defined key_zdftke 
    312             en_25h(:,:,:) = en(:,:,:) 
    313 #endif 
    314 # if defined key_zdfgls 
    315             rmxln_25h(:,:,:) = mxln(:,:,:) 
    316 #endif 
    317             cnt_25h = 1 
    318             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 
    319  
     268         zw3d(:,:,:) = rmxln_25h(:,:,:)*tmask(:,:,:) + zmdi*(1.0-tmask(:,:,:)) 
     269         CALL iom_put( "mxln25h",zw3d) 
     270#endif 
     271         ! 
     272         ! After the write reset the values to cnt=1 and sum values equal current value  
     273         tn_25h  (:,:,:) = tsn (:,:,:,jp_tem) 
     274         sn_25h  (:,:,:) = tsn (:,:,:,jp_sal) 
     275         sshn_25h(:,:)   = sshn(:,:) 
     276         un_25h  (:,:,:) = un  (:,:,:) 
     277         vn_25h  (:,:,:) = vn  (:,:,:) 
     278         wn_25h  (:,:,:) = wn  (:,:,:) 
     279         avt_25h (:,:,:) = avt (:,:,:) 
     280         avm_25h (:,:,:) = avm (:,:,:) 
     281# if defined key_zdfgls || defined key_zdftke 
     282         en_25h  (:,:,:) = en  (:,:,:) 
     283#endif 
     284# if defined key_zdfgls 
     285         rmxln_25h(:,:,:) = mxln(:,:,:) 
     286#endif 
     287         cnt_25h = 1 
     288         IF(lwp)  WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value ',   & 
     289            &                                    'and cnt_25h to one for overlapping average',cnt_25h 
     290         ! 
    320291      ENDIF !  cnt_25h .EQ. 25 .AND.  MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000 
    321  
    322  
     292      ! 
    323293   END SUBROUTINE dia_25h  
    324294 
  • branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/DIA/diacfl.F90

    r6140 r6347  
    151151            WRITE(numout,*) 'dia_cfl     : Maximum Courant number information for the run:' 
    152152            WRITE(numout,*) '~~~~~~~~~~~~' 
    153             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), ')' 
     153            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cu', cu_max,   & 
     154               &                                                           'at (i,j,k) = (', cu_loc(1), cu_loc(2), cu_loc(3), ')' 
    154155            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cu_max) 
    155             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), ')' 
     156            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cv', cv_max,   & 
     157               &                                                           'at (i,j,k) = (', cv_loc(1), cv_loc(2), cv_loc(3), ')' 
    156158            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cv_max) 
    157             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), ')' 
     159            WRITE(numout,FMT='(12x,a12,7x,f6.4,5x,a16,i4,1x,i4,1x,i4,a1)') 'Run Max Cw', cw_max,   & 
     160               &                                                           'at (i,j,k) = (', cw_loc(1), cw_loc(2), cw_loc(3), ')' 
    158161            WRITE(numout,FMT='(12x,a8,11x,f7.1)') ' => dt/C', dt*(1.0/cw_max) 
    159162 
  • branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/DOM/daymod.F90

    r6140 r6347  
    142142 
    143143      ! control print 
    144       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 = ',   & 
    145            &                   nyear, '/', nmonth, '/', nday, '  nsec_day:', nsec_day, '  nsec_week:', nsec_week, '  & 
    146            &                   nsec_month:', nsec_month , '  nsec_year:' , nsec_year 
     144      IF(lwp) WRITE(numout,'(a,i6,a,i2,a,i2,a,i8,a,i8,a,i8,a,i8)')                                  & 
     145         &  ' =======>> 1/2 time step before the start of the run DATE Y/M/D = ', nyear,            & 
     146         &             '/', nmonth, '/', nday, '  nsec_day:', nsec_day,'  nsec_week:', nsec_week,   & 
     147         &             '  nsec_month:', nsec_month ,'  nsec_year:' , nsec_year 
    147148 
    148149      ! 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

    r6140 r6347  
    215215      ! -------------------------------------- 
    216216         DO jk = 1,jpk 
    217             DO jj=1,jpj 
    218                DO ji=1,jpi 
    219                   IF (tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN 
    220                      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) ) 
     217            DO jj = 1, jpj 
     218               DO ji = 1, jpi 
     219                  IF( tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp) THEN 
     220                     e3t_n(ji,jj,jk) = e3t_0(ji,jj,jk)       & 
     221                        &            * ( 1._wp + sshn(ji,jj) / ( ht_0(ji,jj) + 1._wp - ssmask(ji,jj) ) * tmask(ji,jj,jk) ) 
    221222                  ENDIF 
    222223               END DO 
     
    386387      !PM: Is this IF needed since change to VVL by default 
    387388      IF (.NOT.ln_linssh) THEN 
    388          DO jk = 2,jpk-1 
    389             DO jj = 1,jpj 
    390                DO ji = 1,jpi 
    391                   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 
     389         DO jk = 2, jpkm1 
     390            DO jj = 1, jpj 
     391               DO ji = 1, jpi 
     392                  IF( zwmaskn(ji,jj,jk) * zwmaskb(ji,jj,jk) == 1._wp .AND.   & 
     393                     & ( tmask(ji,jj,1) == 0._wp .OR. ptmask_b(ji,jj,1) == 0._wp )  ) THEN 
    392394                     !compute weight 
    393395                     zdzp1 = MAX(0._wp,gdepw_n(ji,jj,jk+1) - pdepw_b(ji,jj,jk+1)) 
     
    436438      CALL wrk_dealloc(jpi,jpj,       zbub   , zbvb    , zbun  , zbvn        )  
    437439      CALL wrk_dealloc(jpi,jpj,       zssh0  , zssh1  , zhu1 , zhv1          )  
    438  
     440      ! 
    439441   END SUBROUTINE iscpl_rst_interpol 
    440442 
     443   !!====================================================================== 
    441444END MODULE iscplrst 
  • branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90

    r6140 r6347  
    764764         istart(idmspc+1) = itime 
    765765 
    766          IF( PRESENT(kstart) .AND. .NOT. ll_depth_spec ) THEN ; istart(1:idmspc) = kstart(1:idmspc) ; icnt(1:idmspc) = kcount(1:idmspc) 
     766         IF( PRESENT(kstart) .AND. .NOT.ll_depth_spec ) THEN  
     767            istart(1:idmspc) = kstart(1:idmspc)  
     768            icnt  (1:idmspc) = kcount(1:idmspc) 
    767769         ELSE 
    768             IF(           idom == jpdom_unknown ) THEN                                                ; icnt(1:idmspc) = idimsz(1:idmspc) 
     770            IF(          idom == jpdom_unknown ) THEN 
     771               icnt(1:idmspc) = idimsz(1:idmspc) 
    769772            ELSE  
    770773               IF( .NOT. PRESENT(pv_r1d) ) THEN   !   not a 1D array 
    771                   IF(     idom == jpdom_data    ) THEN 
     774                  IF(    idom == jpdom_data    ) THEN 
    772775                     jstartrow = 1 
    773776                     IF( luse_jattr ) THEN 
     
    776779                     ENDIF 
    777780                     istart(1:2) = (/ mig(1), mjg(1) + jstartrow - 1 /)  ! icnt(1:2) done below 
    778                   ELSEIF( idom == jpdom_global  ) THEN ; istart(1:2) = (/ nimpp , njmpp  /)  ! icnt(1:2) done below 
     781                  ELSEIF( idom == jpdom_global  ) THEN   ;  istart(1:2) = (/ nimpp , njmpp  /)  ! icnt(1:2) done below 
    779782                  ENDIF 
    780783                  ! we do not read the overlap                     -> we start to read at nldi, nldj 
    781784! JMM + SM: ugly patch before getting the new version of lib_mpp) 
    782785!                  IF( idom /= jpdom_local_noovlap )   istart(1:2) = istart(1:2) + (/ nldi - 1, nldj - 1 /) 
    783                   IF( llnoov .AND. idom /= jpdom_local_noovlap ) istart(1:2) = istart(1:2) + (/ nldi - 1, nldj - 1 /) 
     786                  IF( llnoov .AND. idom /= jpdom_local_noovlap )   istart(1:2) = istart(1:2) + (/ nldi - 1, nldj - 1 /) 
    784787                  ! we do not read the overlap and the extra-halos -> from nldi to nlei and from nldj to nlej  
    785788! JMM + SM: ugly patch before getting the new version of lib_mpp) 
     
    789792                  ENDIF 
    790793                  IF( PRESENT(pv_r3d) ) THEN 
    791                      IF( idom == jpdom_data ) THEN                                  ; icnt(3) = jpkdta 
    792                      ELSE IF( ll_depth_spec .AND. PRESENT(kstart) ) THEN            ; istart(3) = kstart(3); icnt(3) = kcount(3) 
    793                      ELSE                                                           ; icnt(3) = jpk 
     794                     IF( idom == jpdom_data ) THEN                          icnt(3) = jpkdta 
     795                     ELSE IF( ll_depth_spec .AND. PRESENT(kstart) ) THEN   ;   istart(3) = kstart(3)   ;  icnt(3) = kcount(3) 
     796                     ELSE                                                  icnt(3) = jpk 
    794797                     ENDIF 
    795798                  ENDIF 
  • branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90

    r4624 r6347  
    88   !!             -   ! 2004-11  (C. Talandier)  add albedo_init 
    99   !!             -   ! 2001-06  (M. Vancoppenolle) LIM 3.0 
    10    !!             -   ! 2006-08  (G. Madec)  cleaning for surface module 
     10   !!            3.2  ! 2006-08  (G. Madec)  cleaning for surface module 
     11   !!            3.6  ! 2016-01  (C. Rousset) new parameterization for sea ice albedo 
    1112   !!---------------------------------------------------------------------- 
    1213 
     
    1718   !!---------------------------------------------------------------------- 
    1819   USE phycst         ! physical constants 
     20   ! 
    1921   USE in_out_manager ! I/O manager 
    2022   USE lib_mpp        ! MPP library 
     
    2527   PRIVATE 
    2628 
    27    PUBLIC   albedo_ice   ! routine called sbcice_lim.F90 
    28    PUBLIC   albedo_oce   ! routine called by ??? 
    29  
    30    INTEGER  ::   albd_init = 0      !: control flag for initialization 
    31    REAL(wp) ::   zzero     = 0.e0   ! constant values 
    32    REAL(wp) ::   zone      = 1.e0   !    "       " 
    33  
    34    REAL(wp) ::   c1     = 0.05    ! constants values 
    35    REAL(wp) ::   c2     = 0.10    !    "        " 
    36    REAL(wp) ::   rmue   = 0.40    !  cosine of local solar altitude 
    37  
     29   PUBLIC   albedo_ice   ! called sbcice_lim.F90 
     30   PUBLIC   albedo_oce   ! called by sbccpl.F90 
     31 
     32   INTEGER  ::   albd_init = 0          ! control flag for initialization 
     33   
     34   REAL(wp) ::   rmue      = 0.400_wp   ! cosine of local solar altitude 
     35   REAL(wp) ::   ralb_oce  = 0.066_wp   ! ocean or lead albedo (Pegau and Paulson, Ann. Glac. 2001) 
     36   REAL(wp) ::   c1        = 0.050_wp   ! snow thickness (only for nn_ice_alb=0) 
     37   REAL(wp) ::   c2        = 0.100_wp   !  "        " 
     38   REAL(wp) ::   rcloud    = 0.060_wp   ! cloud effect on albedo (only-for nn_ice_alb=0) 
     39  
    3840   !                             !!* namelist namsbc_alb 
    39    REAL(wp) ::   rn_cloud         !  cloudiness effect on snow or ice albedo (Grenfell & Perovich, 1984) 
    40 #if defined key_lim3 
    41    REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers) 
    42 #else 
    43    REAL(wp) ::   rn_albice        !  albedo of melting ice in the arctic and antarctic (Shine & Hendersson-Sellers) 
    44 #endif 
    45    REAL(wp) ::   rn_alphd         !  coefficients for linear interpolation used to compute 
    46    REAL(wp) ::   rn_alphdi        !  albedo between two extremes values (Pyane, 1972) 
    47    REAL(wp) ::   rn_alphc         !  
     41   INTEGER  ::   nn_ice_alb       ! type of albedo parameterization  
     42   REAL(wp) ::   rn_albice        ! albedo of bare puddled ice 
    4843 
    4944   !!---------------------------------------------------------------------- 
     
    5954      !!           
    6055      !! ** Purpose :   Computation of the albedo of the snow/ice system  
    61       !!                as well as the ocean one 
    6256      !!        
    63       !! ** Method  : - Computation of the albedo of snow or ice (choose the  
    64       !!                rignt one by a large number of tests 
    65       !!              - Computation of the albedo of the ocean 
    66       !! 
    67       !! References :   Shine and Hendersson-Sellers 1985, JGR, 90(D1), 2243-2250. 
     57      !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb) 
     58      !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies 
     59      !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005) 
     60      !!                                                                           and Grenfell & Perovich (JGR 2004) 
     61      !!                Description of scheme 1: 
     62      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005) 
     63      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999 
     64      !!                     0-5cm  : linear function of ice thickness 
     65      !!                     5-150cm: log    function of ice thickness 
     66      !!                     > 150cm: constant 
     67      !!                  2) Albedo dependency on snow thickness follows the findings from Grenfell & Perovich (2004) 
     68      !!                     i.e. it increases as -EXP(-snw_thick/0.02) during freezing and -EXP(-snw_thick/0.03) during melting 
     69      !!                  3) Albedo dependency on clouds is speculated from measurements of Grenfell and Perovich (2004) 
     70      !!                     i.e. cloudy-clear albedo depend on cloudy albedo following a 2d order polynomial law 
     71      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice 
     72      !! 
     73      !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions: 
     74      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo 
     75      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger  
     76      !!                     under melting conditions than under freezing conditions 
     77      !!                  3) the evolution of ice albedo as a function of ice thickness shows   
     78      !!                     3 sharp inflexion points (at 5cm, 100cm and 150cm) that look highly unrealistic 
     79      !! 
     80      !! References :   Shine & Henderson-Sellers 1985, JGR, 90(D1), 2243-2250. 
     81      !!                Brandt et al. 2005, J. Climate, vol 18 
     82      !!                Grenfell & Perovich 2004, JGR, vol 109  
    6883      !!---------------------------------------------------------------------- 
    6984      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice      !  ice surface temperature (Kelvin) 
     
    7388      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os   !  albedo of ice under overcast sky 
    7489      !! 
    75       INTEGER  ::   ji, jj, jl    ! dummy loop indices 
    76       INTEGER  ::   ijpl          ! number of ice categories (3rd dim of ice input arrays) 
    77       REAL(wp) ::   zalbpsnm      ! albedo of ice under clear sky when snow is melting 
    78       REAL(wp) ::   zalbpsnf      ! albedo of ice under clear sky when snow is freezing 
    79       REAL(wp) ::   zalbpsn       ! albedo of snow/ice system when ice is coverd by snow 
    80       REAL(wp) ::   zalbpic       ! albedo of snow/ice system when ice is free of snow 
    81       REAL(wp) ::   zithsn        ! = 1 for hsn >= 0 ( ice is cov. by snow ) ; = 0 otherwise (ice is free of snow) 
    82       REAL(wp) ::   zitmlsn       ! = 1 freezinz snow (pt_ice >=rt0_snow) ; = 0 melting snow (pt_ice<rt0_snow) 
    83       REAL(wp) ::   zihsc1        ! = 1 hsn <= c1 ; = 0 hsn > c1 
    84       REAL(wp) ::   zihsc2        ! = 1 hsn >= c2 ; = 0 hsn < c2 
    85       !! 
    86       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalbfz    ! = rn_alphdi for freezing ice ; = rn_albice for melting ice 
    87       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zficeth   !  function of ice thickness 
     90      INTEGER  ::   ji, jj, jl   ! dummy loop indices 
     91      INTEGER  ::   ijpl         ! number of ice categories (3rd dim of ice input arrays) 
     92      REAL(wp) ::   ralb_im, ralb_sf, ralb_sm, ralb_if 
     93      REAL(wp) ::   zswitch, z1_c1, z1_c2 
     94      REAL(wp) ::   zalb_sm, zalb_sf, zalb_st   ! albedo of snow melting, freezing, total 
     95      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalb_it   ! intermediate variable & albedo of ice (snow free) 
    8896      !!--------------------------------------------------------------------- 
     97 
     98      ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories 
    8999       
    90       ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories 
    91  
    92       CALL wrk_alloc( jpi,jpj,ijpl, zalbfz, zficeth ) 
     100      CALL wrk_alloc( jpi,jpj,ijpl,   zalb, zalb_it ) 
    93101 
    94102      IF( albd_init == 0 )   CALL albedo_init      ! initialization  
    95103 
    96       !--------------------------- 
    97       !  Computation of  zficeth 
    98       !--------------------------- 
    99       ! ice free of snow and melts 
    100       WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalbfz(:,:,:) = rn_albice 
    101       ELSE WHERE                                              ;   zalbfz(:,:,:) = rn_alphdi 
    102       END  WHERE 
    103  
    104       WHERE     ( 1.5  < ph_ice                     )  ;  zficeth = zalbfz 
    105       ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zficeth = 0.472  + 2.0 * ( zalbfz - 0.472 ) * ( ph_ice - 1.0 ) 
    106       ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zficeth = 0.2467 + 0.7049 * ph_ice              & 
    107          &                                                                 - 0.8608 * ph_ice * ph_ice     & 
    108          &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice 
    109       ELSE WHERE                                       ;  zficeth = 0.1    + 3.6    * ph_ice 
    110       END WHERE 
    111  
    112 !!gm old code 
    113 !      DO jl = 1, ijpl 
    114 !         DO jj = 1, jpj 
    115 !            DO ji = 1, jpi 
    116 !               IF( ph_ice(ji,jj,jl) > 1.5 ) THEN 
    117 !                  zficeth(ji,jj,jl) = zalbfz(ji,jj,jl) 
    118 !               ELSEIF( ph_ice(ji,jj,jl) > 1.0  .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN 
    119 !                  zficeth(ji,jj,jl) = 0.472 + 2.0 * ( zalbfz(ji,jj,jl) - 0.472 ) * ( ph_ice(ji,jj,jl) - 1.0 ) 
    120 !               ELSEIF( ph_ice(ji,jj,jl) > 0.05 .AND. ph_ice(ji,jj,jl) <= 1.0 ) THEN 
    121 !                  zficeth(ji,jj,jl) = 0.2467 + 0.7049 * ph_ice(ji,jj,jl)                               & 
    122 !                     &                    - 0.8608 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl)                 & 
    123 !                     &                    + 0.3812 * ph_ice(ji,jj,jl) * ph_ice(ji,jj,jl) * ph_ice (ji,jj,jl) 
    124 !               ELSE 
    125 !                  zficeth(ji,jj,jl) = 0.1 + 3.6 * ph_ice(ji,jj,jl)  
    126 !               ENDIF 
    127 !            END DO 
    128 !         END DO 
    129 !      END DO 
    130 !!gm end old code 
    131104       
    132       !-----------------------------------------------  
    133       !    Computation of the snow/ice albedo system  
    134       !-------------------------- --------------------- 
     105      SELECT CASE ( nn_ice_alb )    ! Select the albedo parameterization 
     106      ! 
     107      !              !------------------------------------------ 
     108      CASE( 0 )      !  Shine and Henderson-Sellers (1985) 
     109         !           !------------------------------------------ 
     110         ! 
     111         ralb_sf = 0.80       ! dry snow 
     112         ralb_sm = 0.65       ! melting snow 
     113         ralb_if = 0.72       ! bare frozen ice 
     114         ralb_im = rn_albice  ! bare puddled ice  
     115         !  
     116         !  Computation of ice albedo (free of snow) 
     117         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb(:,:,:) = ralb_im 
     118         ELSE WHERE                                              ;   zalb(:,:,:) = ralb_if 
     119         END  WHERE 
    135120       
    136       !    Albedo of snow-ice for clear sky. 
    137       !-----------------------------------------------     
    138       DO jl = 1, ijpl 
    139          DO jj = 1, jpj 
    140             DO ji = 1, jpi 
    141                !  Case of ice covered by snow.              
    142                !                                        !  freezing snow         
    143                zihsc1   = 1.0 - MAX( zzero , SIGN( zone , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 
    144                zalbpsnf = ( 1.0 - zihsc1 ) * (  zficeth(ji,jj,jl)                                             & 
    145                   &                           + ph_snw(ji,jj,jl) * ( rn_alphd - zficeth(ji,jj,jl) ) / c1  )   & 
    146                   &     +         zihsc1   * rn_alphd   
    147                !                                        !  melting snow                 
    148                zihsc2   = MAX( zzero , SIGN( zone , ph_snw(ji,jj,jl) - c2 ) ) 
    149                zalbpsnm = ( 1.0 - zihsc2 ) * ( rn_albice + ph_snw(ji,jj,jl) * ( rn_alphc - rn_albice ) / c2 )   & 
    150                   &     +         zihsc2   *   rn_alphc  
    151                ! 
    152                zitmlsn  =  MAX( zzero , SIGN( zone , pt_ice(ji,jj,jl) - rt0_snow ) )    
    153                zalbpsn  =  zitmlsn * zalbpsnm + ( 1.0 - zitmlsn ) * zalbpsnf 
    154              
    155                !  Case of ice free of snow. 
    156                zalbpic  = zficeth(ji,jj,jl)  
    157              
    158                ! albedo of the system    
    159                zithsn   = 1.0 - MAX( zzero , SIGN( zone , - ph_snw(ji,jj,jl) ) ) 
    160                pa_ice_cs(ji,jj,jl) =  zithsn * zalbpsn + ( 1.0 - zithsn ) *  zalbpic 
     121         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb 
     122         ELSE WHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 ) 
     123         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )  ;  zalb_it = 0.2467 + 0.7049 * ph_ice              & 
     124            &                                                                 - 0.8608 * ph_ice * ph_ice     & 
     125            &                                                                 + 0.3812 * ph_ice * ph_ice * ph_ice 
     126         ELSE WHERE                                       ;  zalb_it = 0.1    + 3.6    * ph_ice 
     127         END WHERE 
     128         ! 
     129         DO jl = 1, ijpl 
     130            DO jj = 1, jpj 
     131               DO ji = 1, jpi 
     132                  ! freezing snow 
     133                  ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > c2 
     134                  !                                        !  freezing snow         
     135                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - c1 ) ) ) 
     136                  zalb_sf   = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  & 
     137                     &                           + ph_snw(ji,jj,jl) * ( ralb_sf - zalb_it(ji,jj,jl) ) / c1  )   & 
     138                     &        +         zswitch   * ralb_sf   
     139                  ! 
     140                  ! melting snow 
     141                  ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > c2 
     142                  zswitch   = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - c2 ) ) 
     143                  zalb_sm = ( 1._wp - zswitch ) * ( ralb_im + ph_snw(ji,jj,jl) * ( ralb_sm - ralb_im ) / c2 )   & 
     144                      &     +         zswitch   *   ralb_sm  
     145                  ! 
     146                  ! snow albedo 
     147                  zswitch  =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )    
     148                  zalb_st  =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
     149                  ! 
     150                  ! Ice/snow albedo 
     151                  zswitch   = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
     152                  pa_ice_cs(ji,jj,jl) =  zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl) 
     153                  ! 
     154               END DO 
    161155            END DO 
    162156         END DO 
    163       END DO 
    164        
    165       !    Albedo of snow-ice for overcast sky. 
    166       !----------------------------------------------   
    167       pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rn_cloud       ! Oberhuber correction 
    168       ! 
    169       CALL wrk_dealloc( jpi,jpj,ijpl, zalbfz, zficeth ) 
     157         ! 
     158         pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + rcloud       ! Oberhuber correction for overcast sky 
     159         ! 
     160         !              !------------------------------------------ 
     161      CASE( 1 )         !  New parameterization (2016) 
     162         !              !------------------------------------------ 
     163         ! 
     164         ralb_im = rn_albice  ! bare puddled ice 
     165! compilation of values from literature 
     166         ralb_sf = 0.85      ! dry snow 
     167         ralb_sm = 0.75      ! melting snow 
     168         ralb_if = 0.60      ! bare frozen ice 
     169! Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved 
     170!         ralb_sf = 0.85       ! dry snow 
     171!         ralb_sm = 0.72       ! melting snow 
     172!         ralb_if = 0.65       ! bare frozen ice 
     173! Brandt et al 2005 (East Antarctica) 
     174!         ralb_sf = 0.87      ! dry snow 
     175!         ralb_sm = 0.82      ! melting snow 
     176!         ralb_if = 0.54      ! bare frozen ice 
     177!  
     178         !  Computation of ice albedo (free of snow) 
     179         z1_c1 = 1._wp / ( LOG(1.5_wp) - LOG(0.05_wp) )  
     180         z1_c2 = 1._wp / 0.05_wp 
     181         WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0_ice )   ;   zalb = ralb_im 
     182         ELSE WHERE                                              ;   zalb = ralb_if 
     183         END  WHERE 
     184         ! 
     185         WHERE     ( 1.5  < ph_ice                     )  ;  zalb_it = zalb 
     186         ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )  ;  zalb_it = zalb     + ( 0.18_wp - zalb     ) * z1_c1 *  & 
     187            &                                                                     ( LOG(1.5_wp) - LOG(ph_ice) ) 
     188         ELSE WHERE                                       ;  zalb_it = ralb_oce + ( 0.18_wp - ralb_oce ) * z1_c2 * ph_ice 
     189         END WHERE 
     190         ! 
     191         z1_c1 = 1._wp / 0.02_wp 
     192         z1_c2 = 1._wp / 0.03_wp 
     193         !  Computation of the snow/ice albedo 
     194         DO jl = 1, ijpl 
     195            DO jj = 1, jpj 
     196               DO ji = 1, jpi 
     197                  zalb_sf = ralb_sf - ( ralb_sf - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 ); 
     198                  zalb_sm = ralb_sm - ( ralb_sm - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 ); 
     199                  ! 
     200                  ! snow albedo 
     201                  zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )    
     202                  zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
     203                  ! 
     204                  ! Ice/snow albedo    
     205                  zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
     206                  pa_ice_os(ji,jj,jl) = zswitch *  zalb_it(ji,jj,jl) + ( 1._wp - zswitch ) * zalb_st 
     207 
     208              END DO 
     209            END DO 
     210         END DO 
     211         ! Effect of the clouds (2d order polynomial) 
     212         pa_ice_cs = pa_ice_os - ( - 0.1010_wp * pa_ice_os * pa_ice_os + 0.1933_wp * pa_ice_os - 0.0148_wp )  
     213!!gm better coding ? 
     214!         pa_ice_cs = 0.1010_wp * pa_ice_os + 0.8067_wp ) * pa_ice_os + 0.0148_wp          
     215! and even better: merge this implicit do loop with the previous one 
     216!!gm end 
     217      END SELECT 
     218      ! 
     219      CALL wrk_dealloc( jpi,jpj,ijpl,   zalb, zalb_it ) 
    170220      ! 
    171221   END SUBROUTINE albedo_ice 
     
    181231      REAL(wp), DIMENSION(:,:), INTENT(out) ::   pa_oce_cs   !  albedo of ocean under clear sky 
    182232      !! 
    183       REAL(wp) ::   zcoef   ! local scalar 
    184       !!---------------------------------------------------------------------- 
    185       ! 
    186       zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )      ! Parameterization of Briegled and Ramanathan, 1982  
    187       pa_oce_cs(:,:) = zcoef                
    188       pa_oce_os(:,:)  = 0.06                         ! Parameterization of Kondratyev, 1969 and Payne, 1972 
     233      REAL(wp) :: zcoef  
     234      !!---------------------------------------------------------------------- 
     235      ! 
     236      zcoef = 0.05 / ( 1.1 * rmue**1.4 + 0.15 )   ! Parameterization of Briegled and Ramanathan, 1982 
     237      pa_oce_cs(:,:) = zcoef  
     238      pa_oce_os(:,:) = 0.06                       ! Parameterization of Kondratyev, 1969 and Payne, 1972 
    189239      ! 
    190240   END SUBROUTINE albedo_oce 
     
    200250      !!---------------------------------------------------------------------- 
    201251      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    202       NAMELIST/namsbc_alb/ rn_cloud, rn_albice, rn_alphd, rn_alphdi, rn_alphc 
     252      NAMELIST/namsbc_alb/ nn_ice_alb, rn_albice  
    203253      !!---------------------------------------------------------------------- 
    204254      ! 
     
    219269         WRITE(numout,*) '~~~~~~~' 
    220270         WRITE(numout,*) '   Namelist namsbc_alb : albedo ' 
    221          WRITE(numout,*) '      correction for snow and ice albedo                  rn_cloud  = ', rn_cloud 
    222          WRITE(numout,*) '      albedo of melting ice in the arctic and antarctic   rn_albice = ', rn_albice 
    223          WRITE(numout,*) '      coefficients for linear                             rn_alphd  = ', rn_alphd 
    224          WRITE(numout,*) '      interpolation used to compute albedo                rn_alphdi = ', rn_alphdi 
    225          WRITE(numout,*) '      between two extremes values (Pyane, 1972)           rn_alphc  = ', rn_alphc 
     271         WRITE(numout,*) '      choose the albedo parameterization                  nn_ice_alb = ', nn_ice_alb 
     272         WRITE(numout,*) '      albedo of bare puddled ice                          rn_albice  = ', rn_albice 
    226273      ENDIF 
    227274      ! 
  • branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90

    r6140 r6347  
    1111   !!            3.2  !  2009-04  (G. Madec & NEMO team)  
    1212   !!            3.6  !  2012-05  (C. Rousset) store attenuation coef for use in ice model  
    13    !!            3.7  !  2015-11  (G. Madec, A. Coward)  remove optimisation for fix volume  
     13   !!             -   !  2015-12  (O. Aumont, J. Jouanno, C. Ethe) use vertical profile of chlorophyll 
     14   !!            3.7  !  2016-01  (G. Madec, A. Coward)  remove optimisation for fix volume  
    1415   !!---------------------------------------------------------------------- 
    1516 
     
    5556   INTEGER , PUBLIC ::   nksr         !: levels below which the light cannot penetrate (depth larger than 391 m) 
    5657  
    57    INTEGER, PARAMETER ::   np_RGB  = 1   ! R-G-B     light penetration with constant Chlorophyll 
    58    INTEGER, PARAMETER ::   np_RGBc = 2   ! R-G-B     light penetration with Chlorophyll data 
    59    INTEGER, PARAMETER ::   np_2BD  = 3   ! 2 bands   light penetration 
    60    INTEGER, PARAMETER ::   np_BIO  = 4   ! bio-model light penetration 
     58   INTEGER, PARAMETER ::   np_RGB    = 1   ! R-G-B     light penetration with constant Chlorophyll 
     59   INTEGER, PARAMETER ::   np_2BD    = 2   ! 2 bands   light penetration 
     60   INTEGER, PARAMETER ::   np_BIO    = 3   ! bio-model light penetration 
    6161   ! 
    6262   INTEGER  ::   nqsr    ! user choice of the type of light penetration 
    6363   REAL(wp) ::   xsi0r   ! inverse of rn_si0 
    6464   REAL(wp) ::   xsi1r   ! inverse of rn_si1 
     65    
     66   REAL(wp) ::   rChl_0 = 0.05_wp   ! value of Chlorophyll used in case of constant Chlorophyll 
    6567   ! 
    6668   REAL(wp) , DIMENSION(3,61)           ::   rkrgb    ! tabulated attenuation coefficients for RGB absorption 
    67    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_chl   ! structure of input Chl (file informations, fields read) 
     69   TYPE(FLD), DIMENSION(:), ALLOCATABLE ::   sf_chl   ! structure of input Chl (file informations, fields read) 
    6870 
    6971   !! * Substitutions 
     
    100102      !! Reference  : Jerlov, N. G., 1968 Optical Oceanography, Elsevier, 194pp. 
    101103      !!              Lengaigne et al. 2007, Clim. Dyn., V28, 5, 503-516. 
     104      !!              Morel, A. and Berthon, J.F., 1989, Limnol Oceanogr 34(8), 1545-1562 
    102105      !!---------------------------------------------------------------------- 
    103106      INTEGER, INTENT(in) ::   kt     ! ocean time-step 
     
    107110      REAL(wp) ::   zchl, zcoef, z1_2        ! local scalars 
    108111      REAL(wp) ::   zc0 , zc1 , zc2 , zc3    !    -         - 
    109       REAL(wp) ::   zzc0, zzc1, zzc2, zzc3   !    -         - 
    110112      REAL(wp) ::   zz0 , zz1                !    -         - 
     113      REAL(wp) ::   zCb, zCmax, zze, z1_ze, zpsi, zpsimax, zdelpsi, zCtot, zCze 
     114      REAL(wp) ::   zlogc, zlogc2, zlogc3  
    111115      REAL(wp), POINTER, DIMENSION(:,:)   :: zekb, zekg, zekr 
    112       REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, ztrdt 
     116      REAL(wp), POINTER, DIMENSION(:,:,:) :: ze0, ze1, ze2, ze3, zea, zchl3d, ztrdt 
    113117      REAL(wp), POINTER, DIMENSION(:,:,:) :: zetot 
    114118      !!---------------------------------------------------------------------- 
     
    149153      !                         !--------------------------------! 
    150154      ! 
    151       CASE( np_BIO )                   !==  bio-model fluxes  ==! 
     155      CASE( np_BIO )                !==  bio-model fluxes  ==! 
    152156         ! 
    153157         DO jk = 1, nksr 
     
    155159         END DO 
    156160         ! 
    157       CASE( np_RGB , np_RGBc )         !==  R-G-B fluxes  ==! 
    158          ! 
    159          CALL wrk_alloc( jpi,jpj,       zekb, zekg, zekr        )  
    160          CALL wrk_alloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea )  
    161          ! 
    162          IF( nqsr == np_RGBc ) THEN          !*  Variable Chlorophyll 
     161      CASE( np_RGB )                !==  R-G-B fluxes  ==! 
     162         ! 
     163         CALL wrk_alloc( jpi,jpj,       zekb, zekg, zekr                )  
     164         CALL wrk_alloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d )  
     165         ! 
     166         SELECT CASE( nn_chldta )         ! set 3D chlorophyll field 
     167         ! 
     168         CASE( 0 )                           ! constant  
     169            DO jk = 1, nksr + 1 
     170               zchl3d(:,:,jk) = rChl_0 
     171            END DO 
     172            ! 
     173         CASE( 1 )                           ! surface chlorophyl data spread uniformly on the vertical 
    163174            CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step 
    164             DO jj = 2, jpjm1                       ! Separation in R-G-B depending of the surface Chl 
     175            DO jk = 1, nksr + 1                    ! uniform vertical profile 
     176               zchl3d(:,:,jk) = sf_chl(1)%fnow(:,:,1)  
     177            END DO 
     178            ! 
     179         CASE( 2 )                           ! surface chlorophyl data + Morel and Berthon (1989) profile 
     180            CALL fld_read( kt, 1, sf_chl )         ! Read Chl data and provides it at the current time step 
     181            DO jj = 2, jpjm1                       ! Chl profile = F( surface Chl value) 
    165182               DO ji = fs_2, fs_jpim1 
    166                   zchl = MIN( 10. , MAX( 0.03, sf_chl(1)%fnow(ji,jj,1) ) ) 
    167                   irgb = NINT( 41 + 20.*LOG10(zchl) + 1.e-15 ) 
    168                   zekb(ji,jj) = rkrgb(1,irgb) 
    169                   zekg(ji,jj) = rkrgb(2,irgb) 
    170                   zekr(ji,jj) = rkrgb(3,irgb) 
     183                  zchl    = sf_chl(1)%fnow(ji,jj,1) 
     184                  zCtot   =  40.6_wp *  zchl**0.459 
     185                  zze     = 568.2_wp * zCtot**(-0.746) 
     186                  IF( zze > 102. )   zze = 200.0 * zCtot**(-0.293) 
     187                  zlogc   = LOG( zchl ) 
     188!!gm : instead of this : 
     189                  zlogc2  = zlogc * zlogc 
     190                  zlogc3  = zlogc * zlogc * zlogc 
     191                  zCb     = 0.768 + 0.087 * zlogc - 0.179 * zlogc2 - 0.025 * zlogc3 
     192                  zCmax   = 0.299 - 0.289 * zlogc + 0.579 * zlogc2 
     193                  zpsimax = 0.6   - 0.640 * zlogc + 0.021 * zlogc2 + 0.115 * zlogc3 
     194                  zdelpsi = 0.710 + 0.159 * zlogc + 0.021 * zlogc2 
     195!!gm faster & more precise: 
     196!                  zCb     = 0.768 + zlogc * (  (  0.087 + zlogc * (- 0.179 - zlogc * 0.025 )  ) 
     197!                  zCmax   = 0.299 + zlogc * (   - 0.289 + zlogc *    0.579  ) 
     198!                  zpsimax = 0.6   + zlogc * (  (- 0.640 + zlogc * (  0.021 + zlogc * 0.115 )  ) 
     199!                  zdelpsi = 0.710 + zlogc * (     0.159 + zlogc *    0.021  ) 
     200!!gm end           
     201                  zCze    = 1.12_wp * (zchl)**0.803  
     202                  z1_ze   = 1._wp / zze 
     203                  DO jk = 1, nksr + 1 
     204                     zpsi = gdept_n(ji,jj,jk) * z1_ze 
     205                     zchl3d(ji,jj,jk) = zCze * ( zCb + zCmax * EXP( -( (zpsi - zpsimax) / zdelpsi )**2 ) ) 
     206                  END DO 
    171207               END DO 
    172208            END DO 
    173          ELSE                                !* constant chrlorophyll 
    174             zchl = 0.05                            ! constant chlorophyll 
    175             !                                      ! Separation in R-G-B depending of the chlorophyll 
    176             irgb = NINT( 41 + 20.*LOG10( zchl ) + 1.e-15 ) 
    177             DO jj = 2, jpjm1 
    178                DO ji = fs_2, fs_jpim1 
    179                   zekb(ji,jj) = rkrgb(1,irgb)                       
    180                   zekg(ji,jj) = rkrgb(2,irgb) 
    181                   zekr(ji,jj) = rkrgb(3,irgb) 
    182                END DO 
    183             END DO 
    184          ENDIF 
     209            ! 
     210         END SELECT 
    185211         ! 
    186212         zcoef  = ( 1. - rn_abs ) / 3._wp    !* surface equi-partition in R-G-B 
     
    195221         END DO 
    196222         ! 
    197          DO jk = 2, nksr+1                   !* interior equi-partition in R-G-B 
     223         DO jk = 2, nksr+1                   !* interior partition in R-G-B function of 3D Chl 
     224            DO jj = 2, jpjm1 
     225               DO ji = fs_2, fs_jpim1 
     226                  zchl = MIN( 10._wp , MAX( 0.03_wp , zchl3d(ji,jj,jk) ) ) 
     227                  irgb = NINT( 41._wp + 20._wp * LOG10(zchl) + 1.e-15 ) 
     228                  zekb(ji,jj) = rkrgb(1,irgb) 
     229                  zekg(ji,jj) = rkrgb(2,irgb) 
     230                  zekr(ji,jj) = rkrgb(3,irgb) 
     231               END DO 
     232            END DO 
    198233            DO jj = 2, jpjm1 
    199234               DO ji = fs_2, fs_jpim1 
     
    219254         END DO 
    220255         ! 
    221          CALL wrk_dealloc( jpi,jpj,        zekb, zekg, zekr        )  
    222          CALL wrk_dealloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea )  
     256         CALL wrk_dealloc( jpi,jpj,       zekb, zekg, zekr                )  
     257         CALL wrk_dealloc( jpi,jpj,jpk,   ze0, ze1, ze2, ze3, zea, zchl3d )  
    223258         ! 
    224259      CASE( np_2BD  )            !==  2-bands fluxes  ==! 
     
    309344      INTEGER  ::   ji, jj, jk                  ! dummy loop indices 
    310345      INTEGER  ::   ios, irgb, ierror, ioptio   ! local integer 
    311       REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars 
    312       REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      - 
     346!      REAL(wp) ::   zz0, zc0 , zc1, zcoef      ! local scalars 
     347!      REAL(wp) ::   zz1, zc2 , zc3, zchl       !   -      - 
    313348      ! 
    314349      CHARACTER(len=100) ::   cn_dir   ! Root directory for location of ssr files 
     
    321356      IF( nn_timing == 1 )   CALL timing_start('tra_qsr_init') 
    322357      ! 
    323       REWIND( numnam_ref )              ! Namelist namtra_qsr in reference     namelist 
     358      REWIND( numnam_ref )       ! Namelist namtra_qsr in reference     namelist 
    324359      READ  ( numnam_ref, namtra_qsr, IOSTAT = ios, ERR = 901) 
    325360901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in reference namelist', lwp ) 
    326361      ! 
    327       REWIND( numnam_cfg )              ! Namelist namtra_qsr in configuration namelist 
     362      REWIND( numnam_cfg )       ! Namelist namtra_qsr in configuration namelist 
    328363      READ  ( numnam_cfg, namtra_qsr, IOSTAT = ios, ERR = 902 ) 
    329364902   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtra_qsr in configuration namelist', lwp ) 
    330365      IF(lwm) WRITE ( numond, namtra_qsr ) 
    331366      ! 
    332       IF(lwp) THEN                ! control print 
     367      IF(lwp) THEN               ! control print 
    333368         WRITE(numout,*) 
    334369         WRITE(numout,*) 'tra_qsr_init : penetration of the surface solar radiation' 
     
    339374         WRITE(numout,*) '      bio-model            light penetration       ln_qsr_bio = ', ln_qsr_bio 
    340375         WRITE(numout,*) '      light penetration for ice-model (LIM3)       ln_qsr_ice = ', ln_qsr_ice 
    341          WRITE(numout,*) '      RGB : Chl data (=1) or cst value (=0)        nn_chldta  = ', nn_chldta 
     376         WRITE(numout,*) '      RGB : Chl data (=1,2) or cst value (=0)      nn_chldta  = ', nn_chldta 
    342377         WRITE(numout,*) '      RGB & 2 bands: fraction of light (rn_si1)    rn_abs     = ', rn_abs 
    343378         WRITE(numout,*) '      RGB & 2 bands: shortess depth of extinction  rn_si0     = ', rn_si0 
     
    346381      ENDIF 
    347382      ! 
    348       ioptio = 0                    ! Parameter control 
     383      ioptio = 0                 ! Parameter control 
    349384      IF( ln_qsr_rgb  )   ioptio = ioptio + 1 
    350385      IF( ln_qsr_2bd  )   ioptio = ioptio + 1 
     
    355390      ! 
    356391      IF( ln_qsr_rgb .AND. nn_chldta == 0 )   nqsr = np_RGB  
    357       IF( ln_qsr_rgb .AND. nn_chldta == 1 )   nqsr = np_RGBc 
    358392      IF( ln_qsr_2bd                      )   nqsr = np_2BD 
    359393      IF( ln_qsr_bio                      )   nqsr = np_BIO 
    360394      ! 
    361       !                             ! Initialisation 
     395      !                          ! Initialisation 
    362396      xsi0r = 1._wp / rn_si0 
    363397      xsi1r = 1._wp / rn_si1 
     
    365399      SELECT CASE( nqsr ) 
    366400      !                                
    367       CASE( np_RGB , np_RGBc )         !==  Red-Green-Blue light penetration  ==! 
     401      CASE( np_RGB )             !==  Red-Green-Blue light penetration  ==! 
    368402         !                              
    369403         IF(lwp)   WRITE(numout,*) '   R-G-B   light penetration ' 
     
    375409         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 
    376410         ! 
    377          IF( nqsr == np_RGBc ) THEN                ! Chl data : set sf_chl structure 
    378             IF(lwp) WRITE(numout,*) '        Chlorophyll read in a file' 
     411         SELECT CASE( nn_chldta )         ! set 3D chlorophyll field 
     412         CASE( 0 )                           ! constant 
     413            IF(lwp) WRITE(numout,*) '        constant Chlorophyll set to rChl_0 =', rChl_0  
     414            ! 
     415         CASE( 1 , 2 )                       ! 3D chlorophyl field : read 2D surface data 
     416            ! 
     417            IF(lwp) WRITE(numout,*) '        surface 2D Chlorophyll field read in a file' 
    379418            ALLOCATE( sf_chl(1), STAT=ierror ) 
    380419            IF( ierror > 0 ) THEN 
     
    386425            CALL fld_fill( sf_chl, (/ sn_chl /), cn_dir, 'tra_qsr_init',   & 
    387426               &           'Solar penetration function of read chlorophyll', 'namtra_qsr' ) 
    388          ENDIF 
    389          IF( nqsr == np_RGB ) THEN                 ! constant Chl 
    390             IF(lwp) WRITE(numout,*) '        Constant Chlorophyll concentration = 0.05' 
    391          ENDIF 
    392          ! 
    393       CASE( np_2BD )                   !==  2 bands light penetration  ==! 
     427            ! 
     428            IF( lwp .AND. nn_chldta == 1 )   WRITE(numout,*) '        profile of chlorophyll : Chl(z) = Chl(z=0)' 
     429            IF( lwp .AND. nn_chldta == 2 )   WRITE(numout,*) '        profile of chlorophyll : Chl(z) = Func[Chl(z=0)]' 
     430            ! 
     431         END SELECT 
     432         ! 
     433      CASE( np_2BD )             !==  2 bands light penetration  ==! 
    394434         ! 
    395435         IF(lwp)  WRITE(numout,*) '   2 bands light penetration' 
     
    398438         IF(lwp) WRITE(numout,*) '        level of light extinction = ', nksr, ' ref depth = ', gdepw_1d(nksr+1), ' m' 
    399439         ! 
    400       CASE( np_BIO )                   !==  BIO light penetration  ==! 
     440      CASE( np_BIO )                         !==  BIO light penetration  ==! 
    401441         ! 
    402442         IF(lwp) WRITE(numout,*) '   bio-model light penetration' 
  • branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfddm.F90

    r6140 r6347  
    174174                  &                             +  0.15 * zrau(ji,jj)          * zmskd2(ji,jj)  ) 
    175175               ! add to the eddy viscosity coef. previously computed 
     176# if defined key_zdftmx_new 
     177               ! key_zdftmx_new: New internal wave-driven param: use avs value computed by zdftmx 
     178               avs (ji,jj,jk) = avs(ji,jj,jk) + zavfs + zavds 
     179# else 
    176180               avs (ji,jj,jk) = avt(ji,jj,jk) + zavfs + zavds 
     181# endif 
    177182               avt (ji,jj,jk) = avt(ji,jj,jk) + zavft + zavdt 
    178183               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

    r6140 r6347  
    7878      ! 
    7979      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
    80       INTEGER  ::   iikn, iiki, ikt, imkt   ! local integer 
     80      INTEGER  ::   iikn, iiki, ikt   ! local integer 
    8181      REAL(wp) ::   zN2_c        ! local scalar 
    8282      INTEGER, POINTER, DIMENSION(:,:) ::   imld   ! 2D workspace 
     
    123123            iiki = imld(ji,jj) 
    124124            iikn = nmln(ji,jj) 
    125             imkt = mikt(ji,jj) 
    126125            hmld (ji,jj) = gdepw_n(ji,jj,iiki  ) * ssmask(ji,jj)    ! Turbocline depth  
    127126            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

    r6140 r6347  
    323323                  zwlc = zind * rn_lc * zus * SIN( rpi * gdepw_n(ji,jj,jk) / zhlc(ji,jj) ) 
    324324                  !                                           ! TKE Langmuir circulation source term 
    325                   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) 
     325                  en(ji,jj,jk) = en(ji,jj,jk) + rdt * (1._wp - fr_i(ji,jj) ) * ( zwlc * zwlc * zwlc )   & 
     326                     &                              / zhlc(ji,jj) * wmask(ji,jj,jk) * tmask(ji,jj,1) 
    326327               END DO 
    327328            END DO 
     
    732733      ! 
    733734      ri_cri   = 2._wp    / ( 2._wp + rn_ediss / rn_ediff )   ! resulting critical Richardson number 
     735# if defined key_zdftmx_new 
     736      ! key_zdftmx_new: New internal wave-driven param: specified value of rn_emin & rmxl_min are used 
     737      rn_emin  = 1.e-10_wp 
     738      rmxl_min = 1.e-03_wp 
     739      IF(lwp) THEN                  ! Control print 
     740         WRITE(numout,*) 
     741         WRITE(numout,*) 'zdf_tke_init :  New tidal mixing case: force rn_emin = 1.e-10 and rmxl_min = 1.e-3 ' 
     742         WRITE(numout,*) '~~~~~~~~~~~~' 
     743      ENDIF 
     744# else 
    734745      rmxl_min = 1.e-6_wp / ( rn_ediff * SQRT( rn_emin ) )    ! resulting minimum length to recover molecular viscosity 
     746# endif 
    735747      ! 
    736748      IF(lwp) THEN                    !* Control print 
  • branches/2016/dev_r6325_SIMPLIF_1/NEMOGCM/NEMO/OPA_SRC/ZDF/zdftmx.F90

    r6140 r6347  
    541541   END SUBROUTINE zdf_tmx_init 
    542542 
     543#elif defined key_zdftmx_new 
     544   !!---------------------------------------------------------------------- 
     545   !!   'key_zdftmx_new'               Internal wave-driven vertical mixing 
     546   !!---------------------------------------------------------------------- 
     547   !!   zdf_tmx       : global     momentum & tracer Kz with wave induced Kz 
     548   !!   zdf_tmx_init  : global     momentum & tracer Kz with wave induced Kz 
     549   !!---------------------------------------------------------------------- 
     550   USE oce            ! ocean dynamics and tracers variables 
     551   USE dom_oce        ! ocean space and time domain variables 
     552   USE zdf_oce        ! ocean vertical physics variables 
     553   USE zdfddm         ! ocean vertical physics: double diffusive mixing 
     554   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
     555   USE eosbn2         ! ocean equation of state 
     556   USE phycst         ! physical constants 
     557   USE prtctl         ! Print control 
     558   USE in_out_manager ! I/O manager 
     559   USE iom            ! I/O Manager 
     560   USE lib_mpp        ! MPP library 
     561   USE wrk_nemo       ! work arrays 
     562   USE timing         ! Timing 
     563   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     564 
     565   IMPLICIT NONE 
     566   PRIVATE 
     567 
     568   PUBLIC   zdf_tmx         ! called in step module  
     569   PUBLIC   zdf_tmx_init    ! called in nemogcm module  
     570   PUBLIC   zdf_tmx_alloc   ! called in nemogcm module 
     571 
     572   LOGICAL, PUBLIC, PARAMETER ::   lk_zdftmx = .TRUE.    !: wave-driven mixing flag 
     573 
     574   !                       !!* Namelist  namzdf_tmx : internal wave-driven mixing * 
     575   INTEGER  ::  nn_zpyc     ! pycnocline-intensified mixing energy proportional to N (=1) or N^2 (=2) 
     576   LOGICAL  ::  ln_mevar    ! variable (=T) or constant (=F) mixing efficiency 
     577   LOGICAL  ::  ln_tsdiff   ! account for differential T/S wave-driven mixing (=T) or not (=F) 
     578 
     579   REAL(wp) ::  r1_6 = 1._wp / 6._wp 
     580 
     581   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ebot_tmx     ! power available from high-mode wave breaking (W/m2) 
     582   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   epyc_tmx     ! power available from low-mode, pycnocline-intensified wave breaking (W/m2) 
     583   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ecri_tmx     ! power available from low-mode, critical slope wave breaking (W/m2) 
     584   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hbot_tmx     ! WKB decay scale for high-mode energy dissipation (m) 
     585   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hcri_tmx     ! decay scale for low-mode critical slope dissipation (m) 
     586   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   emix_tmx     ! local energy density available for mixing (W/kg) 
     587   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bflx_tmx     ! buoyancy flux Kz * N^2 (W/kg) 
     588   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   pcmap_tmx    ! vertically integrated buoyancy flux (W/m2) 
     589   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_ratio    ! S/T diffusivity ratio (only for ln_tsdiff=T) 
     590   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   zav_wave     ! Internal wave-induced diffusivity 
     591 
     592   !! * Substitutions 
     593#  include "zdfddm_substitute.h90" 
     594#  include "vectopt_loop_substitute.h90" 
     595   !!---------------------------------------------------------------------- 
     596   !! NEMO/OPA 4.0 , NEMO Consortium (2016) 
     597   !! $Id$ 
     598   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     599   !!---------------------------------------------------------------------- 
     600CONTAINS 
     601 
     602   INTEGER FUNCTION zdf_tmx_alloc() 
     603      !!---------------------------------------------------------------------- 
     604      !!                ***  FUNCTION zdf_tmx_alloc  *** 
     605      !!---------------------------------------------------------------------- 
     606      ALLOCATE(     ebot_tmx(jpi,jpj),  epyc_tmx(jpi,jpj),  ecri_tmx(jpi,jpj)    ,   & 
     607      &             hbot_tmx(jpi,jpj),  hcri_tmx(jpi,jpj),  emix_tmx(jpi,jpj,jpk),   & 
     608      &         bflx_tmx(jpi,jpj,jpk), pcmap_tmx(jpi,jpj), zav_ratio(jpi,jpj,jpk),   &  
     609      &         zav_wave(jpi,jpj,jpk), STAT=zdf_tmx_alloc     ) 
     610      ! 
     611      IF( lk_mpp             )   CALL mpp_sum ( zdf_tmx_alloc ) 
     612      IF( zdf_tmx_alloc /= 0 )   CALL ctl_warn('zdf_tmx_alloc: failed to allocate arrays') 
     613   END FUNCTION zdf_tmx_alloc 
     614 
     615 
     616   SUBROUTINE zdf_tmx( kt ) 
     617      !!---------------------------------------------------------------------- 
     618      !!                  ***  ROUTINE zdf_tmx  *** 
     619      !!                    
     620      !! ** Purpose :   add to the vertical mixing coefficients the effect of 
     621      !!              breaking internal waves. 
     622      !! 
     623      !! ** Method  : - internal wave-driven vertical mixing is given by: 
     624      !!                  Kz_wave = min(  100 cm2/s, f(  Reb = emix_tmx /( Nu * N^2 )  ) 
     625      !!              where emix_tmx is the 3D space distribution of the wave-breaking  
     626      !!              energy and Nu the molecular kinematic viscosity. 
     627      !!              The function f(Reb) is linear (constant mixing efficiency) 
     628      !!              if the namelist parameter ln_mevar = F and nonlinear if ln_mevar = T. 
     629      !! 
     630      !!              - Compute emix_tmx, the 3D power density that allows to compute 
     631      !!              Reb and therefrom the wave-induced vertical diffusivity. 
     632      !!              This is divided into three components: 
     633      !!                 1. Bottom-intensified low-mode dissipation at critical slopes 
     634      !!                     emix_tmx(z) = ( ecri_tmx / rau0 ) * EXP( -(H-z)/hcri_tmx ) 
     635      !!                                   / ( 1. - EXP( - H/hcri_tmx ) ) * hcri_tmx 
     636      !!              where hcri_tmx is the characteristic length scale of the bottom  
     637      !!              intensification, ecri_tmx a map of available power, and H the ocean depth. 
     638      !!                 2. Pycnocline-intensified low-mode dissipation 
     639      !!                     emix_tmx(z) = ( epyc_tmx / rau0 ) * ( sqrt(rn2(z))^nn_zpyc ) 
     640      !!                                   / SUM( sqrt(rn2(z))^nn_zpyc * e3w(z) ) 
     641      !!              where epyc_tmx is a map of available power, and nn_zpyc 
     642      !!              is the chosen stratification-dependence of the internal wave 
     643      !!              energy dissipation. 
     644      !!                 3. WKB-height dependent high mode dissipation 
     645      !!                     emix_tmx(z) = ( ebot_tmx / rau0 ) * rn2(z) * EXP(-z_wkb(z)/hbot_tmx) 
     646      !!                                   / SUM( rn2(z) * EXP(-z_wkb(z)/hbot_tmx) * e3w(z) ) 
     647      !!              where hbot_tmx is the characteristic length scale of the WKB bottom  
     648      !!              intensification, ebot_tmx is a map of available power, and z_wkb is the 
     649      !!              WKB-stretched height above bottom defined as 
     650      !!                    z_wkb(z) = H * SUM( sqrt(rn2(z'>=z)) * e3w(z'>=z) ) 
     651      !!                                 / SUM( sqrt(rn2(z'))    * e3w(z')    ) 
     652      !! 
     653      !!              - update the model vertical eddy viscosity and diffusivity:  
     654      !!                     avt  = avt  +    av_wave 
     655      !!                     avm  = avm  +    av_wave 
     656      !!                     avmu = avmu + mi(av_wave) 
     657      !!                     avmv = avmv + mj(av_wave) 
     658      !! 
     659      !!              - if namelist parameter ln_tsdiff = T, account for differential mixing: 
     660      !!                     avs  = avt  +    av_wave * diffusivity_ratio(Reb) 
     661      !! 
     662      !! ** Action  : - Define emix_tmx used to compute internal wave-induced mixing 
     663      !!              - avt, avs, avm, avmu, avmv increased by internal wave-driven mixing     
     664      !! 
     665      !! References :  de Lavergne et al. 2015, JPO; 2016, in prep. 
     666      !!---------------------------------------------------------------------- 
     667      INTEGER, INTENT(in) ::   kt   ! ocean time-step  
     668      ! 
     669      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     670      REAL(wp) ::   ztpc         ! scalar workspace 
     671      REAL(wp), DIMENSION(:,:)  , POINTER ::  zfact     ! Used for vertical structure 
     672      REAL(wp), DIMENSION(:,:)  , POINTER ::  zhdep     ! Ocean depth 
     673      REAL(wp), DIMENSION(:,:,:), POINTER ::  zwkb      ! WKB-stretched height above bottom 
     674      REAL(wp), DIMENSION(:,:,:), POINTER ::  zweight   ! Weight for high mode vertical distribution 
     675      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_t     ! Molecular kinematic viscosity (T grid) 
     676      REAL(wp), DIMENSION(:,:,:), POINTER ::  znu_w     ! Molecular kinematic viscosity (W grid) 
     677      REAL(wp), DIMENSION(:,:,:), POINTER ::  zReb      ! Turbulence intensity parameter 
     678      !!---------------------------------------------------------------------- 
     679      ! 
     680      IF( nn_timing == 1 )   CALL timing_start('zdf_tmx') 
     681      ! 
     682      CALL wrk_alloc( jpi,jpj,       zfact, zhdep ) 
     683      CALL wrk_alloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
     684 
     685      !                          ! ----------------------------- ! 
     686      !                          !  Internal wave-driven mixing  !  (compute zav_wave) 
     687      !                          ! ----------------------------- ! 
     688      !                              
     689      !                        !* Critical slope mixing: distribute energy over the time-varying ocean depth, 
     690      !                                                 using an exponential decay from the seafloor. 
     691      DO jj = 1, jpj                ! part independent of the level 
     692         DO ji = 1, jpi 
     693            zhdep(ji,jj) = gdepw_0(ji,jj,mbkt(ji,jj)+1)       ! depth of the ocean 
     694            zfact(ji,jj) = rau0 * (  1._wp - EXP( -zhdep(ji,jj) / hcri_tmx(ji,jj) )  ) 
     695            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ecri_tmx(ji,jj) / zfact(ji,jj) 
     696         END DO 
     697      END DO 
     698 
     699      DO jk = 2, jpkm1              ! complete with the level-dependent part 
     700         emix_tmx(:,:,jk) = zfact(:,:) * (  EXP( ( fsde3w(:,:,jk  ) - zhdep(:,:) ) / hcri_tmx(:,:) )                      & 
     701            &                             - EXP( ( fsde3w(:,:,jk-1) - zhdep(:,:) ) / hcri_tmx(:,:) )  ) * wmask(:,:,jk)   & 
     702            &                          / ( fsde3w(:,:,jk) - fsde3w(:,:,jk-1) ) 
     703      END DO 
     704 
     705      !                        !* Pycnocline-intensified mixing: distribute energy over the time-varying  
     706      !                        !* ocean depth as proportional to sqrt(rn2)^nn_zpyc 
     707 
     708      SELECT CASE ( nn_zpyc ) 
     709 
     710      CASE ( 1 )               ! Dissipation scales as N (recommended) 
     711 
     712         zfact(:,:) = 0._wp 
     713         DO jk = 2, jpkm1              ! part independent of the level 
     714            zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
     715         END DO 
     716 
     717         DO jj = 1, jpj 
     718            DO ji = 1, jpi 
     719               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     720            END DO 
     721         END DO 
     722 
     723         DO jk = 2, jpkm1              ! complete with the level-dependent part 
     724            emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
     725         END DO 
     726 
     727      CASE ( 2 )               ! Dissipation scales as N^2 
     728 
     729         zfact(:,:) = 0._wp 
     730         DO jk = 2, jpkm1              ! part independent of the level 
     731            zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
     732         END DO 
     733 
     734         DO jj= 1, jpj 
     735            DO ji = 1, jpi 
     736               IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = epyc_tmx(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     737            END DO 
     738         END DO 
     739 
     740         DO jk = 2, jpkm1              ! complete with the level-dependent part 
     741            emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zfact(:,:) * MAX( 0._wp, rn2(:,:,jk) ) * wmask(:,:,jk) 
     742         END DO 
     743 
     744      END SELECT 
     745 
     746      !                        !* WKB-height dependent mixing: distribute energy over the time-varying  
     747      !                        !* ocean depth as proportional to rn2 * exp(-z_wkb/rn_hbot) 
     748       
     749      zwkb(:,:,:) = 0._wp 
     750      zfact(:,:) = 0._wp 
     751      DO jk = 2, jpkm1 
     752         zfact(:,:) = zfact(:,:) + fse3w(:,:,jk) * SQRT(  MAX( 0._wp, rn2(:,:,jk) )  ) * wmask(:,:,jk) 
     753         zwkb(:,:,jk) = zfact(:,:) 
     754      END DO 
     755 
     756      DO jk = 2, jpkm1 
     757         DO jj = 1, jpj 
     758            DO ji = 1, jpi 
     759               IF( zfact(ji,jj) /= 0 )   zwkb(ji,jj,jk) = zhdep(ji,jj) * ( zfact(ji,jj) - zwkb(ji,jj,jk) )   & 
     760                                            &           * tmask(ji,jj,jk) / zfact(ji,jj) 
     761            END DO 
     762         END DO 
     763      END DO 
     764      zwkb(:,:,1) = zhdep(:,:) * tmask(:,:,1) 
     765 
     766      zweight(:,:,:) = 0._wp 
     767      DO jk = 2, jpkm1 
     768         zweight(:,:,jk) = MAX( 0._wp, rn2(:,:,jk) ) * hbot_tmx(:,:) * wmask(:,:,jk)                    & 
     769            &   * (  EXP( -zwkb(:,:,jk) / hbot_tmx(:,:) ) - EXP( -zwkb(:,:,jk-1) / hbot_tmx(:,:) )  ) 
     770      END DO 
     771 
     772      zfact(:,:) = 0._wp 
     773      DO jk = 2, jpkm1              ! part independent of the level 
     774         zfact(:,:) = zfact(:,:) + zweight(:,:,jk) 
     775      END DO 
     776 
     777      DO jj = 1, jpj 
     778         DO ji = 1, jpi 
     779            IF( zfact(ji,jj) /= 0 )   zfact(ji,jj) = ebot_tmx(ji,jj) / ( rau0 * zfact(ji,jj) ) 
     780         END DO 
     781      END DO 
     782 
     783      DO jk = 2, jpkm1              ! complete with the level-dependent part 
     784         emix_tmx(:,:,jk) = emix_tmx(:,:,jk) + zweight(:,:,jk) * zfact(:,:) * wmask(:,:,jk)   & 
     785            &                                / ( fsde3w(:,:,jk) - fsde3w(:,:,jk-1) ) 
     786      END DO 
     787 
     788 
     789      ! Calculate molecular kinematic viscosity 
     790      znu_t(:,:,:) = 1.e-4_wp * (  17.91_wp - 0.53810_wp * tsn(:,:,:,jp_tem) + 0.00694_wp * tsn(:,:,:,jp_tem) * tsn(:,:,:,jp_tem)  & 
     791         &                                  + 0.02305_wp * tsn(:,:,:,jp_sal)  ) * tmask(:,:,:) * r1_rau0 
     792      DO jk = 2, jpkm1 
     793         znu_w(:,:,jk) = 0.5_wp * ( znu_t(:,:,jk-1) + znu_t(:,:,jk) ) * wmask(:,:,jk) 
     794      END DO 
     795 
     796      ! Calculate turbulence intensity parameter Reb 
     797      DO jk = 2, jpkm1 
     798         zReb(:,:,jk) = emix_tmx(:,:,jk) / MAX( 1.e-20_wp, znu_w(:,:,jk) * rn2(:,:,jk) ) 
     799      END DO 
     800 
     801      ! Define internal wave-induced diffusivity 
     802      DO jk = 2, jpkm1 
     803         zav_wave(:,:,jk) = znu_w(:,:,jk) * zReb(:,:,jk) * r1_6   ! This corresponds to a constant mixing efficiency of 1/6 
     804      END DO 
     805 
     806      IF( ln_mevar ) THEN              ! Variable mixing efficiency case : modify zav_wave in the 
     807         DO jk = 2, jpkm1              ! energetic (Reb > 480) and buoyancy-controlled (Reb <10.224 ) regimes 
     808            DO jj = 1, jpj 
     809               DO ji = 1, jpi 
     810                  IF( zReb(ji,jj,jk) > 480.00_wp ) THEN 
     811                     zav_wave(ji,jj,jk) = 3.6515_wp * znu_w(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     812                  ELSEIF( zReb(ji,jj,jk) < 10.224_wp ) THEN 
     813                     zav_wave(ji,jj,jk) = 0.052125_wp * znu_w(ji,jj,jk) * zReb(ji,jj,jk) * SQRT( zReb(ji,jj,jk) ) 
     814                  ENDIF 
     815               END DO 
     816            END DO 
     817         END DO 
     818      ENDIF 
     819 
     820      DO jk = 2, jpkm1                 ! Bound diffusivity by molecular value and 100 cm2/s 
     821         zav_wave(:,:,jk) = MIN(  MAX( 1.4e-7_wp, zav_wave(:,:,jk) ), 1.e-2_wp  ) * wmask(:,:,jk) 
     822      END DO 
     823 
     824      IF( kt == nit000 ) THEN        !* Control print at first time-step: diagnose the energy consumed by zav_wave 
     825         ztpc = 0._wp 
     826         DO jk = 2, jpkm1 
     827            DO jj = 1, jpj 
     828               DO ji = 1, jpi 
     829                  ztpc = ztpc + fse3w(ji,jj,jk) * e1e2t(ji,jj)   & 
     830                     &         * MAX( 0._wp, rn2(ji,jj,jk) ) * zav_wave(ji,jj,jk) * wmask(ji,jj,jk) * tmask_i(ji,jj) 
     831               END DO 
     832            END DO 
     833         END DO 
     834         IF( lk_mpp )   CALL mpp_sum( ztpc ) 
     835         ztpc = rau0 * ztpc ! Global integral of rauo * Kz * N^2 = power contributing to mixing  
     836  
     837         IF(lwp) THEN 
     838            WRITE(numout,*) 
     839            WRITE(numout,*) 'zdf_tmx : Internal wave-driven mixing (tmx)' 
     840            WRITE(numout,*) '~~~~~~~ ' 
     841            WRITE(numout,*) 
     842            WRITE(numout,*) '      Total power consumption by av_wave: ztpc =  ', ztpc * 1.e-12_wp, 'TW' 
     843         ENDIF 
     844      ENDIF 
     845 
     846      !                          ! ----------------------- ! 
     847      !                          !   Update  mixing coefs  !                           
     848      !                          ! ----------------------- ! 
     849      !       
     850      IF( ln_tsdiff ) THEN          !* Option for differential mixing of salinity and temperature 
     851         DO jk = 2, jpkm1              ! Calculate S/T diffusivity ratio as a function of Reb 
     852            DO jj = 1, jpj 
     853               DO ji = 1, jpi 
     854                  zav_ratio(ji,jj,jk) = ( 0.505_wp + 0.495_wp *                                                                  & 
     855                      &   TANH(    0.92_wp * (   LOG10(  MAX( 1.e-20_wp, zReb(ji,jj,jk) * 5._wp * r1_6 )  ) - 0.60_wp   )    )   & 
     856                      &                 ) * wmask(ji,jj,jk) 
     857               END DO 
     858            END DO 
     859         END DO 
     860         CALL iom_put( "av_ratio", zav_ratio ) 
     861         DO jk = 2, jpkm1           !* update momentum & tracer diffusivity with wave-driven mixing 
     862            fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) * zav_ratio(:,:,jk) 
     863            avt  (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
     864            avm  (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk) 
     865         END DO 
     866         ! 
     867      ELSE                          !* update momentum & tracer diffusivity with wave-driven mixing 
     868         DO jk = 2, jpkm1 
     869            fsavs(:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
     870            avt  (:,:,jk) = avt(:,:,jk) + zav_wave(:,:,jk) 
     871            avm  (:,:,jk) = avm(:,:,jk) + zav_wave(:,:,jk) 
     872         END DO 
     873      ENDIF 
     874 
     875      DO jk = 2, jpkm1              !* update momentum diffusivity at wu and wv points 
     876         DO jj = 2, jpjm1 
     877            DO ji = fs_2, fs_jpim1  ! vector opt. 
     878               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) 
     879               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) 
     880            END DO 
     881         END DO 
     882      END DO 
     883      CALL lbc_lnk( avmu, 'U', 1. )   ;   CALL lbc_lnk( avmv, 'V', 1. )      ! lateral boundary condition 
     884 
     885      !                             !* output internal wave-driven mixing coefficient 
     886      CALL iom_put( "av_wave", zav_wave ) 
     887                                    !* output useful diagnostics: N^2, Kz * N^2 (bflx_tmx),  
     888                                    !  vertical integral of rau0 * Kz * N^2 (pcmap_tmx), energy density (emix_tmx) 
     889      IF( iom_use("bflx_tmx") .OR. iom_use("pcmap_tmx") ) THEN 
     890         bflx_tmx(:,:,:) = MAX( 0._wp, rn2(:,:,:) ) * zav_wave(:,:,:) 
     891         pcmap_tmx(:,:) = 0._wp 
     892         DO jk = 2, jpkm1 
     893            pcmap_tmx(:,:) = pcmap_tmx(:,:) + fse3w(:,:,jk) * bflx_tmx(:,:,jk) * wmask(:,:,jk) 
     894         END DO 
     895         pcmap_tmx(:,:) = rau0 * pcmap_tmx(:,:) 
     896         CALL iom_put( "bflx_tmx", bflx_tmx ) 
     897         CALL iom_put( "pcmap_tmx", pcmap_tmx ) 
     898      ENDIF 
     899      CALL iom_put( "bn2", rn2 ) 
     900      CALL iom_put( "emix_tmx", emix_tmx ) 
     901       
     902      CALL wrk_dealloc( jpi,jpj,       zfact, zhdep ) 
     903      CALL wrk_dealloc( jpi,jpj,jpk,   zwkb, zweight, znu_t, znu_w, zReb ) 
     904 
     905      IF(ln_ctl)   CALL prt_ctl(tab3d_1=zav_wave , clinfo1=' tmx - av_wave: ', tab3d_2=avt, clinfo2=' avt: ', ovlap=1, kdim=jpk) 
     906      ! 
     907      IF( nn_timing == 1 )   CALL timing_stop('zdf_tmx') 
     908      ! 
     909   END SUBROUTINE zdf_tmx 
     910 
     911 
     912   SUBROUTINE zdf_tmx_init 
     913      !!---------------------------------------------------------------------- 
     914      !!                  ***  ROUTINE zdf_tmx_init  *** 
     915      !!                      
     916      !! ** Purpose :   Initialization of the wave-driven vertical mixing, reading 
     917      !!              of input power maps and decay length scales in netcdf files. 
     918      !! 
     919      !! ** Method  : - Read the namzdf_tmx namelist and check the parameters 
     920      !! 
     921      !!              - Read the input data in NetCDF files : 
     922      !!              power available from high-mode wave breaking (mixing_power_bot.nc) 
     923      !!              power available from pycnocline-intensified wave-breaking (mixing_power_pyc.nc) 
     924      !!              power available from critical slope wave-breaking (mixing_power_cri.nc) 
     925      !!              WKB decay scale for high-mode wave-breaking (decay_scale_bot.nc) 
     926      !!              decay scale for critical slope wave-breaking (decay_scale_cri.nc) 
     927      !! 
     928      !! ** input   : - Namlist namzdf_tmx 
     929      !!              - NetCDF files : mixing_power_bot.nc, mixing_power_pyc.nc, mixing_power_cri.nc, 
     930      !!              decay_scale_bot.nc decay_scale_cri.nc 
     931      !! 
     932      !! ** Action  : - Increase by 1 the nstop flag is setting problem encounter 
     933      !!              - Define ebot_tmx, epyc_tmx, ecri_tmx, hbot_tmx, hcri_tmx 
     934      !! 
     935      !! References : de Lavergne et al. 2015, JPO; 2016, in prep. 
     936      !!          
     937      !!---------------------------------------------------------------------- 
     938      INTEGER  ::   ji, jj, jk   ! dummy loop indices 
     939      INTEGER  ::   inum         ! local integer 
     940      INTEGER  ::   ios 
     941      REAL(wp) ::   zbot, zpyc, zcri   ! local scalars 
     942      !! 
     943      NAMELIST/namzdf_tmx_new/ nn_zpyc, ln_mevar, ln_tsdiff 
     944      !!---------------------------------------------------------------------- 
     945      ! 
     946      IF( nn_timing == 1 )  CALL timing_start('zdf_tmx_init') 
     947      ! 
     948      REWIND( numnam_ref )              ! Namelist namzdf_tmx in reference namelist : Wave-driven mixing 
     949      READ  ( numnam_ref, namzdf_tmx_new, IOSTAT = ios, ERR = 901) 
     950901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in reference namelist', lwp ) 
     951      ! 
     952      REWIND( numnam_cfg )              ! Namelist namzdf_tmx in configuration namelist : Wave-driven mixing 
     953      READ  ( numnam_cfg, namzdf_tmx_new, IOSTAT = ios, ERR = 902 ) 
     954902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_tmx in configuration namelist', lwp ) 
     955      IF(lwm) WRITE ( numond, namzdf_tmx_new ) 
     956      ! 
     957      IF(lwp) THEN                  ! Control print 
     958         WRITE(numout,*) 
     959         WRITE(numout,*) 'zdf_tmx_init : internal wave-driven mixing' 
     960         WRITE(numout,*) '~~~~~~~~~~~~' 
     961         WRITE(numout,*) '   Namelist namzdf_tmx_new : set wave-driven mixing parameters' 
     962         WRITE(numout,*) '      Pycnocline-intensified diss. scales as N (=1) or N^2 (=2) = ', nn_zpyc 
     963         WRITE(numout,*) '      Variable (T) or constant (F) mixing efficiency            = ', ln_mevar 
     964         WRITE(numout,*) '      Differential internal wave-driven mixing (T) or not (F)   = ', ln_tsdiff 
     965      ENDIF 
     966       
     967      ! The new wave-driven mixing parameterization elevates avt and avm in the interior, and 
     968      ! ensures that avt remains larger than its molecular value (=1.4e-7). Therefore, avtb should  
     969      ! be set here to a very small value, and avmb to its (uniform) molecular value (=1.4e-6). 
     970      avmb(:) = 1.4e-6_wp        ! viscous molecular value 
     971      avtb(:) = 1.e-10_wp        ! very small diffusive minimum (background avt is specified in zdf_tmx)     
     972      avtb_2d(:,:) = 1.e0_wp     ! uniform  
     973      IF(lwp) THEN                  ! Control print 
     974         WRITE(numout,*) 
     975         WRITE(numout,*) '   Force the background value applied to avm & avt in TKE to be everywhere ',   & 
     976            &               'the viscous molecular value & a very small diffusive value, resp.' 
     977      ENDIF 
     978       
     979      IF( .NOT.lk_zdfddm )   CALL ctl_stop( 'STOP', 'zdf_tmx_init_new : key_zdftmx_new requires key_zdfddm' ) 
     980       
     981      !                             ! allocate tmx arrays 
     982      IF( zdf_tmx_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'zdf_tmx_init : unable to allocate tmx arrays' ) 
     983      ! 
     984      !                             ! read necessary fields 
     985      CALL iom_open('mixing_power_bot',inum)       ! energy flux for high-mode wave breaking [W/m2] 
     986      CALL iom_get  (inum, jpdom_data, 'field', ebot_tmx, 1 )  
     987      CALL iom_close(inum) 
     988      ! 
     989      CALL iom_open('mixing_power_pyc',inum)       ! energy flux for pynocline-intensified wave breaking [W/m2] 
     990      CALL iom_get  (inum, jpdom_data, 'field', epyc_tmx, 1 ) 
     991      CALL iom_close(inum) 
     992      ! 
     993      CALL iom_open('mixing_power_cri',inum)       ! energy flux for critical slope wave breaking [W/m2] 
     994      CALL iom_get  (inum, jpdom_data, 'field', ecri_tmx, 1 ) 
     995      CALL iom_close(inum) 
     996      ! 
     997      CALL iom_open('decay_scale_bot',inum)        ! spatially variable decay scale for high-mode wave breaking [m] 
     998      CALL iom_get  (inum, jpdom_data, 'field', hbot_tmx, 1 ) 
     999      CALL iom_close(inum) 
     1000      ! 
     1001      CALL iom_open('decay_scale_cri',inum)        ! spatially variable decay scale for critical slope wave breaking [m] 
     1002      CALL iom_get  (inum, jpdom_data, 'field', hcri_tmx, 1 ) 
     1003      CALL iom_close(inum) 
     1004 
     1005      ebot_tmx(:,:) = ebot_tmx(:,:) * ssmask(:,:) 
     1006      epyc_tmx(:,:) = epyc_tmx(:,:) * ssmask(:,:) 
     1007      ecri_tmx(:,:) = ecri_tmx(:,:) * ssmask(:,:) 
     1008 
     1009      ! Set once for all to zero the first and last vertical levels of appropriate variables 
     1010      emix_tmx (:,:, 1 ) = 0._wp 
     1011      emix_tmx (:,:,jpk) = 0._wp 
     1012      zav_ratio(:,:, 1 ) = 0._wp 
     1013      zav_ratio(:,:,jpk) = 0._wp 
     1014      zav_wave (:,:, 1 ) = 0._wp 
     1015      zav_wave (:,:,jpk) = 0._wp 
     1016 
     1017      zbot = glob_sum( e1e2t(:,:) * ebot_tmx(:,:) ) 
     1018      zpyc = glob_sum( e1e2t(:,:) * epyc_tmx(:,:) ) 
     1019      zcri = glob_sum( e1e2t(:,:) * ecri_tmx(:,:) ) 
     1020      IF(lwp) THEN 
     1021         WRITE(numout,*) '      High-mode wave-breaking energy:             ', zbot * 1.e-12_wp, 'TW' 
     1022         WRITE(numout,*) '      Pycnocline-intensifed wave-breaking energy: ', zpyc * 1.e-12_wp, 'TW' 
     1023         WRITE(numout,*) '      Critical slope wave-breaking energy:        ', zcri * 1.e-12_wp, 'TW' 
     1024      ENDIF 
     1025      ! 
     1026      IF( nn_timing == 1 )  CALL timing_stop('zdf_tmx_init') 
     1027      ! 
     1028   END SUBROUTINE zdf_tmx_init 
     1029 
    5431030#else 
    5441031   !!---------------------------------------------------------------------- 
Note: See TracChangeset for help on using the changeset viewer.