Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_top_cfg
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_top_cfg	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_top_cfg	(revision 10377)
@@ -90,4 +90,8 @@
 /
 !-----------------------------------------------------------------------
+&namtrc_snk      !  sedimentation of particles
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
 &namtrc_dmp      !   passive tracer newtonian damping   
 !-----------------------------------------------------------------------
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_top_cfg
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_top_cfg	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_top_cfg	(revision 10377)
@@ -89,4 +89,8 @@
 /
 !-----------------------------------------------------------------------
+&namtrc_snk      !  sedimentation of particles
+!-----------------------------------------------------------------------
+/
+!-----------------------------------------------------------------------
 &namtrc_dmp      !   passive tracer newtonian damping   
 !-----------------------------------------------------------------------
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/SHARED/namelist_pisces_ref
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/SHARED/namelist_pisces_ref	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/SHARED/namelist_pisces_ref	(revision 10377)
@@ -51,11 +51,9 @@
    wsbio2max  =  50.      ! Big particles maximum sinking speed
    wsbio2scale =  5000.    ! Big particles length scale of sinking
-   niter1max  =  1        ! Maximum number of iterations for POC
-   niter2max  =  2        ! Maximum number of iterations for GOC
 !                         !  ln_ligand enabled
    wfep       =  0.2      ! FeP sinking speed 
    ldocp      =  1.E-4    ! Phyto ligand production per unit doc 
    ldocz      =  1.E-4    ! Zoo ligand production per unit doc 
-   lthet      =  0.5      ! Proportional loss of ligands due to Fe uptake 
+   lthet      =  1.0      ! Proportional loss of ligands due to Fe uptake 
 !                         !  ln_p5z enabled
    no3rat3    =  0.182    ! N/C ratio in zooplankton
@@ -310,5 +308,5 @@
    xlam1        =  0.005  ! scavenging rate of Iron
    xlamdust     =  150.0  ! Scavenging rate of dust
-   ligand       =  0.6E-9 ! Ligands concentration 
+   ligand       =  0.7E-9 ! Ligands concentration 
    kfep         =  0.01   ! Nanoparticle formation rate constant
 /
@@ -322,5 +320,5 @@
    xsilab    =  0.5       ! Fraction of labile biogenic silica
    feratb    =  10.E-6    ! Fe/C quota in bacteria
-   xkferb    =  2.5E-10   ! Half-saturation constant for bacteria Fe/C
+   xkferb    =  3E-10     ! Half-saturation constant for bacteria Fe/C
 !                         ! ln_p5z
    xremikc   =  0.25      ! remineralization rate of DOC
@@ -372,5 +370,5 @@
    ln_ironsed  =  .true.   ! boolean for Fe input from sediments
    ln_ironice  =  .true.   ! boolean for Fe input from sea ice
-   ln_hydrofe  =  .false.  ! boolean for from hydrothermal vents
+   ln_hydrofe  =  .true.   ! boolean for from hydrothermal vents
    sedfeinput  =  2.e-9    ! Coastal release of Iron
    distcoast   =  5.e3     ! Distance off the coast for Iron from sediments
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/SHARED/namelist_top_ref
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/SHARED/namelist_top_ref	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/cfgs/SHARED/namelist_top_ref	(revision 10377)
@@ -95,4 +95,9 @@
 /
 !-----------------------------------------------------------------------
+&namtrc_snk      !  Sedimentation of particles
+!-----------------------------------------------------------------------
+   nitermax      =  2   !  number of iterations for sedimentation
+/
+!-----------------------------------------------------------------------
 &namtrc_dmp      !   passive tracer newtonian damping                   (ln_trcdmp=T)
 !-----------------------------------------------------------------------
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/main/NEMO_manual.bib
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/main/NEMO_manual.bib	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/main/NEMO_manual.bib	(revision 10377)
@@ -547,4 +547,15 @@
 PAGES = {1285--1297},
 DOI = {10.5194/gmd-8-1285-2015}
+}
+
+@ARTICLE{Breivik_al_JPO2014,
+   AUTHOR = {{\O}yvind Breivik and Peter A.E.M. Janssen and Jean-Raymond Bidlot},
+   YEAR = {2014},
+   TITLE = "{Approximate Stokes Drift Profiles in Deep Water}",
+   JOURNAL = {JPO},
+   VOLUME = {44},
+   NUMBER = {9},
+   DOI = {10.1175/JPO-D-14-0020.1.},
+   PAGES = {2433--2445, arXiv:1406.5039}
 }
 
@@ -1655,4 +1666,14 @@
 }
 
+@TECHREPORT{Janssen_al_TM13,
+ author = {P.A.E.M. Janssen and {\O}. Breivik and K. Mogensen 
+           and F. Vitart and  M. Balmaseda and J.B. Bidlot and 
+           S. Keeley and M. Leut-becher and L. Magnusson and F. Molteni},
+ title = {Air-Sea Interaction and Surface Waves},
+ year = {2013},
+ volume = {712},
+ institution = {ECMWF},
+}
+
 @ARTICLE{Jayne_St_Laurent_GRL01,
   author = {S.R. Jayne and L.C. {St. Laurent}},
@@ -2992,4 +3013,13 @@
   volume = {359},
   pages = {123--129}
+}
+
+@INCOLLECTION{Stokes_1847,
+  author    = {G.G. Stokes},
+  title     = {On the theory of oscillatory waves},
+  booktitle = {Transactions of the Cambridge Philosophy Society},
+  year      = {1847},
+  volume    = {8},
+  Pages     = {441--455}
 }
 
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/chap_SBC.tex
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/chap_SBC.tex	(revision 10377)
@@ -32,12 +32,10 @@
 \end{itemize}
 
-Five different ways to provide the first six fields to the ocean are available which are controlled by
+Four different ways to provide the first six fields to the ocean are available which are controlled by
 namelist \ngn{namsbc} variables:
 an analytical formulation (\np{ln\_ana}\forcode{ = .true.}),
 a flux formulation (\np{ln\_flx}\forcode{ = .true.}),
 a bulk formulae formulation (CORE (\np{ln\_blk\_core}\forcode{ = .true.}),
-CLIO (\np{ln\_blk\_clio}\forcode{ = .true.}) or
-MFS \footnote { Note that MFS bulk formulae compute fluxes only for the ocean component}
-(\np{ln\_blk\_mfs}\forcode{ = .true.}) bulk formulae) and
+CLIO (\np{ln\_blk\_clio}\forcode{ = .true.}) bulk formulae) and
 a coupled or mixed forced/coupled formulation (exchanges with a atmospheric model via the OASIS coupler)
 (\np{ln\_cpl} or \np{ln\_mixcpl}\forcode{ = .true.}). 
@@ -76,7 +74,13 @@
   the transformation of the solar radiation (if provided as daily mean) into a diurnal cycle
   (\np{ln\_dm2dc}\forcode{ = .true.});
-  and a neutral drag coefficient can be read from an external wave model (\np{ln\_cdgw}\forcode{ = .true.}). 
+\item
+  a neutral drag coefficient can be read from an external wave model (\np{ln\_cdgw}\forcode{ = .true.});
+\item
+  the Stokes drift rom an external wave model can be accounted (\np{ln\_sdw}\forcode{ = .true.}); 
+\item
+  the Stokes-Coriolis term can be included (\np{ln\_stcor}\forcode{ = .true.});
+\item
+  the surface stress felt by the ocean can be modified by surface waves (\np{ln\_tauwoc}\forcode{ = .true.}).
 \end{itemize}
-The latter option is possible only in case core or mfs bulk formulas are selected.
 
 In this chapter, we first discuss where the surface boundary condition appears in the model equations.
@@ -593,6 +597,6 @@
 % Bulk formulation
 % ================================================================
-\section[Bulk formulation {(\textit{sbcblk\{\_core,\_clio,\_mfs\}.F90})}]
-			{Bulk formulation {(\protect\mdl{sbcblk\_core}, \protect\mdl{sbcblk\_clio}, \protect\mdl{sbcblk\_mfs})}}
+\section[Bulk formulation {(\textit{sbcblk\{\_core,\_clio\}.F90})}]
+                        {Bulk formulation {(\protect\mdl{sbcblk\_core}, \protect\mdl{sbcblk\_clio})}}
 \label{sec:SBC_blk}
 
@@ -600,8 +604,8 @@
 
 The atmospheric fields used depend on the bulk formulae used.
-Three bulk formulations are available:
-the CORE, the CLIO and the MFS bulk formulea.
+Two bulk formulations are available:
+the CORE and the CLIO bulk formulea.
 The choice is made by setting to true one of the following namelist variable:
-\np{ln\_core} ; \np{ln\_clio} or  \np{ln\_mfs}.
+\np{ln\_core} or \np{ln\_clio}.
 
 Note:
@@ -712,41 +716,4 @@
 the namsbc\_blk\_core or namsbc\_blk\_clio namelist (see \autoref{subsec:SBC_fldread}). 
 
-% -------------------------------------------------------------------------------------------------------------
-%        MFS Bulk formulae
-% -------------------------------------------------------------------------------------------------------------
-\subsection{MFS formulea (\protect\mdl{sbcblk\_mfs}, \protect\np{ln\_mfs}\forcode{ = .true.})}
-\label{subsec:SBC_blk_mfs}
-%------------------------------------------namsbc_mfs----------------------------------------------------
-%
-%\nlst{namsbc_mfs}
-%----------------------------------------------------------------------------------------------------------
-
-The MFS (Mediterranean Forecasting System) bulk formulae have been developed by \citet{Castellari_al_JMS1998}. 
-They have been designed to handle the ECMWF operational data and are currently in use in
-the MFS operational system \citep{Tonani_al_OS08}, \citep{Oddo_al_OS09}.
-The wind stress computation uses a drag coefficient computed according to \citet{Hellerman_Rosenstein_JPO83}.
-The surface boundary condition for temperature involves the balance between
-surface solar radiation, net long-wave radiation, the latent and sensible heat fluxes.
-Solar radiation is dependent on cloud cover and is computed by means of an astronomical formula \citep{Reed_JPO77}.
-Albedo monthly values are from \citet{Payne_JAS72} as means of the values at $40^{o}N$ and $30^{o}N$ for
-the Atlantic Ocean (hence the same latitudinal band of the Mediterranean Sea).
-The net long-wave radiation flux \citep{Bignami_al_JGR95} is a function of
-air temperature, sea-surface temperature, cloud cover and relative humidity.
-Sensible heat and latent heat fluxes are computed by classical bulk formulae parameterised according to
-\citet{Kondo1975}.
-Details on the bulk formulae used can be found in \citet{Maggiore_al_PCE98} and \citet{Castellari_al_JMS1998}.
-
-Options are defined through the \ngn{namsbc\_mfs} namelist variables.
-The required 7 input fields must be provided on the model Grid-T and are:
-\begin{itemize}
-\item          Zonal Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windi})
-\item          Meridional Component of the 10m wind ($ms^{-1}$)  (\np{sn\_windj})
-\item          Total Claud Cover (\%)  (\np{sn\_clc})
-\item          2m Air Temperature ($K$) (\np{sn\_tair})
-\item          2m Dew Point Temperature ($K$)  (\np{sn\_rhm})
-\item          Total Precipitation ${Kg} m^{-2} s^{-1}$ (\np{sn\_prec})
-\item          Mean Sea Level Pressure (${Pa}$) (\np{sn\_msl})
-\end{itemize}
-% -------------------------------------------------------------------------------------------------------------
 % ================================================================
 % Coupled formulation
@@ -1203,4 +1170,191 @@
 since its trajectory data may be spread across multiple files.
 
+% -------------------------------------------------------------------------------------------------------------
+%        Interactions with waves (sbcwave.F90, ln_wave)
+% -------------------------------------------------------------------------------------------------------------
+\section{Interactions with waves (\protect\mdl{sbcwave}, \protect\np{ln\_wave})}
+\label{sec:SBC_wave}
+%------------------------------------------namsbc_wave--------------------------------------------------------
+
+\nlst{namsbc_wave}
+%-------------------------------------------------------------------------------------------------------------
+
+Ocean waves represent the interface between the ocean and the atmosphere, so NEMO is extended to incorporate 
+physical processes related to ocean surface waves, namely the surface stress modified by growth and 
+dissipation of the oceanic wave field, the Stokes-Coriolis force and the Stokes drift impact on mass and 
+tracer advection; moreover the neutral surface drag coefficient from a wave model can be used to evaluate 
+the wind stress.
+
+Physical processes related to ocean surface waves can be accounted by setting the logical variable 
+\np{ln\_wave}\forcode{= .true.} in \ngn{namsbc} namelist. In addition, specific flags accounting for 
+different porcesses should be activated as explained in the following sections.
+
+Wave fields can be provided either in forced or coupled mode:
+\begin{description}
+\item[forced mode]: wave fields should be defined through the \ngn{namsbc\_wave} namelist 
+for external data names, locations, frequency, interpolation and all the miscellanous options allowed by 
+Input Data generic Interface (see \autoref{sec:SBC_input}). 
+\item[coupled mode]: NEMO and an external wave model can be coupled by setting \np{ln\_cpl} \forcode{= .true.} 
+in \ngn{namsbc} namelist and filling the \ngn{namsbc\_cpl} namelist.
+\end{description}
+
+
+% ================================================================
+% Neutral drag coefficient from wave model (ln_cdgw)
+
+% ================================================================
+\subsection{Neutral drag coefficient from wave model (\protect\np{ln\_cdgw})}
+\label{subsec:SBC_wave_cdgw}
+
+The neutral surface drag coefficient provided from an external data source ($i.e.$ a wave
+model), 
+can be used by setting the logical variable \np{ln\_cdgw} \forcode{= .true.} in \ngn{namsbc} namelist. 
+Then using the routine \rou{turb\_ncar} and starting from the neutral drag coefficent provided, 
+the drag coefficient is computed according to the stable/unstable conditions of the 
+air-sea interface following \citet{Large_Yeager_Rep04}. 
+
+
+% ================================================================
+% 3D Stokes Drift (ln_sdw, nn_sdrift)
+% ================================================================
+\subsection{3D Stokes Drift (\protect\np{ln\_sdw, nn\_sdrift})}
+\label{subsec:SBC_wave_sdw}
+
+The Stokes drift is a wave driven mechanism of mass and momentum transport \citep{Stokes_1847}. 
+It is defined as the difference between the average velocity of a fluid parcel (Lagrangian velocity) 
+and the current measured at a fixed point (Eulerian velocity). 
+As waves travel, the water particles that make up the waves travel in orbital motions but 
+without a closed path. Their movement is enhanced at the top of the orbit and slowed slightly 
+at the bottom so the result is a net forward motion of water particles, referred to as the Stokes drift. 
+An accurate evaluation of the Stokes drift and the inclusion of related processes may lead to improved 
+representation of surface physics in ocean general circulation models.
+The Stokes drift velocity $\mathbf{U}_{st}$ in deep water can be computed from the wave spectrum and may be written as: 
+
+\begin{equation} \label{eq:sbc_wave_sdw}
+\mathbf{U}_{st} = \frac{16{\pi^3}} {g} 
+                \int_0^\infty \int_{-\pi}^{\pi} (cos{\theta},sin{\theta}) {f^3}
+                \mathrm{S}(f,\theta) \mathrm{e}^{2kz}\,\mathrm{d}\theta {d}f
+\end{equation}
+
+where: ${\theta}$ is the wave direction, $f$ is the wave intrinsic frequency, 
+$\mathrm{S}($f$,\theta)$ is the 2D frequency-direction spectrum, 
+$k$ is the mean wavenumber defined as: 
+$k=\frac{2\pi}{\lambda}$ (being $\lambda$ the wavelength). \\
+
+In order to evaluate the Stokes drift in a realistic ocean wave field the wave spectral shape is required 
+and its computation quickly becomes expensive as the 2D spectrum must be integrated for each vertical level. 
+To simplify, it is customary to use approximations to the full Stokes profile.
+Three possible parameterizations for the calculation for the approximate Stokes drift velocity profile 
+are included in the code through the \np{nn\_sdrift} parameter once provided the surface Stokes drift 
+$\mathbf{U}_{st |_{z=0}}$ which is evaluated by an external wave model that accurately reproduces the wave spectra 
+and makes possible the estimation of the surface Stokes drift for random directional waves in 
+realistic wave conditions:
+
+\begin{description}
+\item[\np{nn\_sdrift} = 0]: exponential integral profile parameterization proposed by 
+\citet{Breivik_al_JPO2014}:
+
+\begin{equation} \label{eq:sbc_wave_sdw_0a}
+\mathbf{U}_{st} \cong \mathbf{U}_{st |_{z=0}} \frac{\mathrm{e}^{-2k_ez}} {1-8k_ez} 
+\end{equation}
+
+where $k_e$ is the effective wave number which depends on the Stokes transport $T_{st}$ defined as follows:
+
+\begin{equation} \label{eq:sbc_wave_sdw_0b}
+k_e = \frac{|\mathbf{U}_{\left.st\right|_{z=0}}|} {|T_{st}|} 
+\quad \text{and }\
+T_{st} = \frac{1}{16} \bar{\omega} H_s^2 
+\end{equation}
+
+where $H_s$ is the significant wave height and $\omega$ is the wave frequency.
+
+\item[\np{nn\_sdrift} = 1]: velocity profile based on the Phillips spectrum which is considered to be a 
+reasonable estimate of the part of the spectrum most contributing to the Stokes drift velocity near the surface
+\citep{Breivik_al_OM2016}:
+
+\begin{equation} \label{eq:sbc_wave_sdw_1}
+\mathbf{U}_{st} \cong \mathbf{U}_{st |_{z=0}} \Big[exp(2k_pz)-\beta \sqrt{-2 \pi k_pz} 
+\textit{ erf } \Big(\sqrt{-2 k_pz}\Big)\Big]
+\end{equation}
+
+where $erf$ is the complementary error function and $k_p$ is the peak wavenumber.
+
+\item[\np{nn\_sdrift} = 2]: velocity profile based on the Phillips spectrum as for \np{nn\_sdrift} = 1 
+but using the wave frequency from a wave model.
+
+\end{description}
+
+The Stokes drift enters the wave-averaged momentum equation, as well as the tracer advection equations 
+and its effect on the evolution of the sea-surface height ${\eta}$ is considered as follows: 
+
+\begin{equation} \label{eq:sbc_wave_eta_sdw}
+\frac{\partial{\eta}}{\partial{t}} = 
+-\nabla_h \int_{-H}^{\eta} (\mathbf{U} + \mathbf{U}_{st}) dz 
+\end{equation}
+
+The tracer advection equation is also modified in order for Eulerian ocean models to properly account 
+for unresolved wave effect. The divergence of the wave tracer flux equals the mean tracer advection 
+that is induced by the three-dimensional Stokes velocity. 
+The advective equation for a tracer $c$ combining the effects of the mean current and sea surface waves 
+can be formulated as follows: 
+
+\begin{equation} \label{eq:sbc_wave_tra_sdw}
+\frac{\partial{c}}{\partial{t}} = 
+- (\mathbf{U} + \mathbf{U}_{st}) \cdot \nabla{c}
+\end{equation}
+
+
+% ================================================================
+% Stokes-Coriolis term (ln_stcor)
+% ================================================================
+\subsection{Stokes-Coriolis term (\protect\np{ln\_stcor})}
+\label{subsec:SBC_wave_stcor}
+
+In a rotating ocean, waves exert a wave-induced stress on the mean ocean circulation which results 
+in a force equal to $\mathbf{U}_{st}$×$f$, where $f$ is the Coriolis parameter. 
+This additional force may have impact on the Ekman turning of the surface current. 
+In order to include this term, once evaluated the Stokes drift (using one of the 3 possible 
+approximations described in \autoref{subsec:SBC_wave_sdw}), 
+\np{ln\_stcor}\forcode{ = .true.} has to be set.
+
+
+% ================================================================
+% Waves modified stress (ln_tauwoc, ln_tauw)
+% ================================================================
+\subsection{Wave modified sress (\protect\np{ln\_tauwoc, ln\_tauw})} 
+\label{subsec:SBC_wave_tauw}
+
+The surface stress felt by the ocean is the atmospheric stress minus the net stress going 
+into the waves \citep{Janssen_al_TM13}. Therefore, when waves are growing, momentum and energy is spent and is not 
+available for forcing the mean circulation, while in the opposite case of a decaying sea 
+state more momentum is available for forcing the ocean. 
+Only when the sea state is in equilibrium the ocean is forced by the atmospheric stress, 
+but in practice an equilibrium sea state is a fairly rare event. 
+So the atmospheric stress felt by the ocean circulation $\tau_{oc,a}$ can be expressed as: 
+
+\begin{equation} \label{eq:sbc_wave_tauoc}
+\tau_{oc,a} = \tau_a - \tau_w
+\end{equation}
+
+where $\tau_a$ is the atmospheric surface stress;
+$\tau_w$ is the atmospheric stress going into the waves defined as:
+
+\begin{equation} \label{eq:sbc_wave_tauw}
+\tau_w = \rho g \int {\frac{dk}{c_p} (S_{in}+S_{nl}+S_{diss})}
+\end{equation}
+
+where: $c_p$ is the phase speed of the gravity waves,
+$S_{in}$, $S_{nl}$ and $S_{diss}$ are three source terms that represent 
+the physics of ocean waves. The first one, $S_{in}$, describes the generation 
+of ocean waves by wind and therefore represents the momentum and energy transfer 
+from air to ocean waves; the second term $S_{nl}$ denotes 
+the nonlinear transfer by resonant four-wave interactions; while the third term $S_{diss}$ 
+describes the dissipation of waves by processes such as white-capping, large scale breaking 
+eddy-induced damping.
+
+The wave stress derived from an external wave model can be provided either through the normalized 
+wave stress into the ocean by setting \np{ln\_tauwoc}\forcode{ = .true.}, or through the zonal and 
+meridional stress components by setting \np{ln\_tauw}\forcode{ = .true.}.
+
 
 % ================================================================
@@ -1421,24 +1575,4 @@
 \end{description}
 
-% -------------------------------------------------------------------------------------------------------------
-%        Neutral Drag Coefficient from external wave model
-% -------------------------------------------------------------------------------------------------------------
-\subsection[Neutral drag coeff. from external wave model (\protect\mdl{sbcwave})]
-				{Neutral drag coefficient from external wave model (\protect\mdl{sbcwave})}
-\label{subsec:SBC_wave}
-%------------------------------------------namwave----------------------------------------------------
-
-\nlst{namsbc_wave}
-%-------------------------------------------------------------------------------------------------------------
-
-In order to read a neutral drag coefficient, from an external data source ($i.e.$ a wave model),
-the logical variable \np{ln\_cdgw} in \ngn{namsbc} namelist must be set to \forcode{.true.}.
-The \mdl{sbcwave} module containing the routine \np{sbc\_wave} reads the namelist \ngn{namsbc\_wave}
-(for external data names, locations, frequency, interpolation and all the miscellanous options allowed by
-Input Data generic Interface see \autoref{sec:SBC_input}) and
-a 2D field of neutral drag coefficient.
-Then using the routine TURB\_CORE\_1Z or TURB\_CORE\_2Z, and starting from the neutral drag coefficent provided, 
-the drag coefficient is computed according to stable/unstable conditions of the air-sea interface following
-\citet{Large_Yeager_Rep04}.
 
 
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/introduction.tex
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/introduction.tex	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/doc/latex/NEMO/subfiles/introduction.tex	(revision 10377)
@@ -311,4 +311,5 @@
 \item new definition of configurations ;
 \item bulk formulation ;
+\item NEMO-wave large scale interactions ;
 \item ... ; 
 \end{enumerate}
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/IOM/iom_nf90.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/IOM/iom_nf90.F90	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/IOM/iom_nf90.F90	(revision 10377)
@@ -691,5 +691,5 @@
       INTEGER, DIMENSION(4) :: idimid               ! dimensions id
       CHARACTER(LEN=256)    :: clinfo               ! info character
-      CHARACTER(LEN= 12), DIMENSION(4) :: cltmp     ! temporary character
+      CHARACTER(LEN= 12), DIMENSION(5) :: cltmp     ! temporary character
       INTEGER               :: if90id               ! nf90 file identifier
       INTEGER               :: idmy                 ! dummy variable
@@ -716,9 +716,5 @@
          ENDIF
          ! define the dimension variables if it is not already done
-         IF(iom_file(kiomid)%nlev == jpk ) THEN
-          cltmp = (/ 'nav_lon     ', 'nav_lat     ', 'nav_lev     ', 'time_counter' /)
-         ELSE
-          cltmp = (/ 'nav_lon     ', 'nav_lat     ', 'numcat      ', 'time_counter' /)
-         ENDIF
+         cltmp = (/ 'nav_lon', 'nav_lat', 'nav_lev', 'time_counter', 'numcat' /)
          CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(1)), NF90_FLOAT , (/ 1, 2 /), iom_file(kiomid)%nvid(1) ), clinfo)
          CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(2)), NF90_FLOAT , (/ 1, 2 /), iom_file(kiomid)%nvid(2) ), clinfo)
@@ -728,6 +724,13 @@
          iom_file(kiomid)%nvars       = 4
          iom_file(kiomid)%luld(1:4)   = (/ .FALSE., .FALSE., .FALSE., .TRUE. /)
-         iom_file(kiomid)%cn_var(1:4) = cltmp
-         iom_file(kiomid)%ndims(1:4)  = (/ 2, 2, 1, 1 /)  
+         iom_file(kiomid)%cn_var(1:4) = cltmp(1:4)
+         iom_file(kiomid)%ndims(1:4)  = (/ 2, 2, 1, 1 /)
+         IF( NF90_INQ_DIMID( if90id, 'numcat', idmy ) == nf90_noerr ) THEN   ! add a 5th variable corresponding to the 5th dimension
+            CALL iom_nf90_check(NF90_DEF_VAR( if90id, TRIM(cltmp(5)), NF90_FLOAT , (/ 5 /), iom_file(kiomid)%nvid(5) ), clinfo)
+            iom_file(kiomid)%nvars     = 5
+            iom_file(kiomid)%luld(5)   = .FALSE.
+            iom_file(kiomid)%cn_var(5) = cltmp(5)
+            iom_file(kiomid)%ndims(5)  = 1
+         ENDIF
          ! trick: defined to 0 to say that dimension variables are defined but not yet written
          iom_file(kiomid)%dimsz(1, 1)  = 0   
@@ -841,10 +844,7 @@
                CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lat'     , idmy )         , clinfo )
                CALL iom_nf90_check( NF90_PUT_VAR  ( if90id, idmy, gphit(ix1:ix2, iy1:iy2) ), clinfo )
-               IF(iom_file(kiomid)%nlev == jpk ) THEN 
-                  !NEMO
-                  CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lev'     , idmy ), clinfo )
-                  CALL iom_nf90_check( NF90_PUT_VAR  ( if90id, idmy, gdept_1d       ), clinfo )
-               ELSE
-                  CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'numcat'     , idmy ), clinfo)
+               CALL iom_nf90_check( NF90_INQ_VARID( if90id, 'nav_lev'     , idmy ), clinfo )
+               CALL iom_nf90_check( NF90_PUT_VAR  ( if90id, idmy, gdept_1d       ), clinfo )
+               IF( NF90_INQ_VARID( if90id, 'numcat', idmy ) == nf90_noerr ) THEN
                   CALL iom_nf90_check( NF90_PUT_VAR  ( if90id, idmy, (/ (idlv, idlv = 1,iom_file(kiomid)%nlev) /)), clinfo )
                ENDIF
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/geo2ocean.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/geo2ocean.F90	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/OCE/SBC/geo2ocean.F90	(revision 10377)
@@ -10,12 +10,4 @@
    !!            3.7  !  11-2015  (G. Madec)  remove the unused repere and repcmo routines
    !!----------------------------------------------------------------------
-#if defined key_agrif
-!clem: these lines do not seem necessary anymore
-!!DIR$ OPTIMIZE (-O 1)  ! cray formulation
-# if defined __INTEL_COMPILER
-!acc: still breaks on at least one Ivybridge cluster with ifort 17.0.4 without this directive
-!DIR$ OPTIMIZE:1        ! intel formulation
-# endif
-#endif
    !!----------------------------------------------------------------------
    !!   rot_rep       : Rotate the Repere: geographic grid <==> stretched coordinates grid
@@ -81,5 +73,5 @@
          IF(lwp) WRITE(numout,*) ' ~~~~~~~~    '
          !
-         CALL angle       ! initialization of the transformation
+         CALL angle( glamt, gphit, glamu, gphiu, glamv, gphiv, glamf, gphif )       ! initialization of the transformation
          lmust_init = .FALSE.
       ENDIF
@@ -126,5 +118,5 @@
 
 
-   SUBROUTINE angle
+   SUBROUTINE angle( plamt, pphit, plamu, pphiu, plamv, pphiv, plamf, pphif )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE angle  ***
@@ -138,4 +130,8 @@
       !! ** Action  : - gsint, gcost, gsinu, gcosu, gsinv, gcosv, gsinf, gcosf
       !!----------------------------------------------------------------------
+      ! WARNING: for an unexplained reason, we need to pass all glam, gphi arrays as input parameters in
+      !          order to get AGRIF working with -03 compilation option
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) :: plamt, pphit, plamu, pphiu, plamv, pphiv, plamf, pphif  
+      !
       INTEGER  ::   ji, jj   ! dummy loop indices
       INTEGER  ::   ierr     ! local integer
@@ -167,32 +163,32 @@
          DO ji = fs_2, jpi   ! vector opt.
             !                  
-            zlam = glamt(ji,jj)     ! north pole direction & modulous (at t-point)
-            zphi = gphit(ji,jj)
+            zlam = plamt(ji,jj)     ! north pole direction & modulous (at t-point)
+            zphi = pphit(ji,jj)
             zxnpt = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             zynpt = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             znnpt = zxnpt*zxnpt + zynpt*zynpt
             !
-            zlam = glamu(ji,jj)     ! north pole direction & modulous (at u-point)
-            zphi = gphiu(ji,jj)
+            zlam = plamu(ji,jj)     ! north pole direction & modulous (at u-point)
+            zphi = pphiu(ji,jj)
             zxnpu = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             zynpu = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             znnpu = zxnpu*zxnpu + zynpu*zynpu
             !
-            zlam = glamv(ji,jj)     ! north pole direction & modulous (at v-point)
-            zphi = gphiv(ji,jj)
+            zlam = plamv(ji,jj)     ! north pole direction & modulous (at v-point)
+            zphi = pphiv(ji,jj)
             zxnpv = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             zynpv = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             znnpv = zxnpv*zxnpv + zynpv*zynpv
             !
-            zlam = glamf(ji,jj)     ! north pole direction & modulous (at f-point)
-            zphi = gphif(ji,jj)
+            zlam = plamf(ji,jj)     ! north pole direction & modulous (at f-point)
+            zphi = pphif(ji,jj)
             zxnpf = 0. - 2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             zynpf = 0. - 2. * SIN( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )
             znnpf = zxnpf*zxnpf + zynpf*zynpf
             !
-            zlam = glamv(ji,jj  )   ! j-direction: v-point segment direction (around t-point)
-            zphi = gphiv(ji,jj  )
-            zlan = glamv(ji,jj-1)
-            zphh = gphiv(ji,jj-1)
+            zlam = plamv(ji,jj  )   ! j-direction: v-point segment direction (around t-point)
+            zphi = pphiv(ji,jj  )
+            zlan = plamv(ji,jj-1)
+            zphh = pphiv(ji,jj-1)
             zxvvt =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
                &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
@@ -202,8 +198,8 @@
             znvvt = MAX( znvvt, 1.e-14 )
             !
-            zlam = glamf(ji,jj  )   ! j-direction: f-point segment direction (around u-point)
-            zphi = gphif(ji,jj  )
-            zlan = glamf(ji,jj-1)
-            zphh = gphif(ji,jj-1)
+            zlam = plamf(ji,jj  )   ! j-direction: f-point segment direction (around u-point)
+            zphi = pphif(ji,jj  )
+            zlan = plamf(ji,jj-1)
+            zphh = pphif(ji,jj-1)
             zxffu =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
                &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
@@ -213,8 +209,8 @@
             znffu = MAX( znffu, 1.e-14 )
             !
-            zlam = glamf(ji  ,jj)   ! i-direction: f-point segment direction (around v-point)
-            zphi = gphif(ji  ,jj)
-            zlan = glamf(ji-1,jj)
-            zphh = gphif(ji-1,jj)
+            zlam = plamf(ji  ,jj)   ! i-direction: f-point segment direction (around v-point)
+            zphi = pphif(ji  ,jj)
+            zlan = plamf(ji-1,jj)
+            zphh = pphif(ji-1,jj)
             zxffv =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
                &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
@@ -224,8 +220,8 @@
             znffv = MAX( znffv, 1.e-14 )
             !
-            zlam = glamu(ji,jj+1)   ! j-direction: u-point segment direction (around f-point)
-            zphi = gphiu(ji,jj+1)
-            zlan = glamu(ji,jj  )
-            zphh = gphiu(ji,jj  )
+            zlam = plamu(ji,jj+1)   ! j-direction: u-point segment direction (around f-point)
+            zphi = pphiu(ji,jj+1)
+            zlan = plamu(ji,jj  )
+            zphh = pphiu(ji,jj  )
             zxuuf =  2. * COS( rad*zlam ) * TAN( rpi/4. - rad*zphi/2. )   &
                &  -  2. * COS( rad*zlan ) * TAN( rpi/4. - rad*zphh/2. )
@@ -257,17 +253,17 @@
       DO jj = 2, jpjm1
          DO ji = fs_2, jpi   ! vector opt.
-            IF( MOD( ABS( glamv(ji,jj) - glamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
+            IF( MOD( ABS( plamv(ji,jj) - plamv(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
                gsint(ji,jj) = 0.
                gcost(ji,jj) = 1.
             ENDIF
-            IF( MOD( ABS( glamf(ji,jj) - glamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
+            IF( MOD( ABS( plamf(ji,jj) - plamf(ji,jj-1) ), 360. ) < 1.e-8 ) THEN
                gsinu(ji,jj) = 0.
                gcosu(ji,jj) = 1.
             ENDIF
-            IF(      ABS( gphif(ji,jj) - gphif(ji-1,jj) )         < 1.e-8 ) THEN
+            IF(      ABS( pphif(ji,jj) - pphif(ji-1,jj) )         < 1.e-8 ) THEN
                gsinv(ji,jj) = 0.
                gcosv(ji,jj) = 1.
             ENDIF
-            IF( MOD( ABS( glamu(ji,jj) - glamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN
+            IF( MOD( ABS( plamu(ji,jj) - plamu(ji,jj+1) ), 360. ) < 1.e-8 ) THEN
                gsinf(ji,jj) = 0.
                gcosf(ji,jj) = 1.
@@ -457,5 +453,5 @@
          IF(lwp) WRITE(numout,*) ' obs_rot : geographic <--> stretched'
          IF(lwp) WRITE(numout,*) ' ~~~~~~~   coordinate transformation'
-         CALL angle       ! initialization of the transformation
+         CALL angle( glamt, gphit, glamu, gphiu, glamv, gphiv, glamf, gphif )       ! initialization of the transformation
          lmust_init = .FALSE.
       ENDIF
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zmeso.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zmeso.F90	(revision 10377)
@@ -252,4 +252,6 @@
          DEALLOCATE( zw3d )
       ENDIF
+      !
+      IF (ln_ligand)  DEALLOCATE( zz2ligprod )
       !
       IF(ln_ctl)   THEN  ! print mean trends (used for debugging)
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zmicro.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zmicro.F90	(revision 10377)
@@ -205,4 +205,6 @@
       ENDIF
       !
+      IF (ln_ligand)  DEALLOCATE( zzligprod )
+      !
       IF(ln_ctl) THEN      ! print mean trends (used for debugging)
          WRITE(charout, FMT="('micro')")
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsink.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsink.F90	(revision 10377)
@@ -16,4 +16,5 @@
    USE trc             !  passive tracers common variables 
    USE sms_pisces      !  PISCES Source Minus Sink variables
+   USE trcsink         !  General routine to compute sedimentation
    USE prtctl_trc      !  print control for debugging
    USE iom             !  I/O manager
@@ -59,10 +60,7 @@
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) :: kt, knt
-      INTEGER  ::   ji, jj, jk, jit
-      INTEGER  ::   iiter1, iiter2
-      REAL(wp) ::   zagg1, zagg2, zagg3, zagg4
-      REAL(wp) ::   zagg , zaggfe, zaggdoc, zaggdoc2, zaggdoc3
-      REAL(wp) ::   zfact, zwsmax, zmax
+      INTEGER  ::   ji, jj, jk
       CHARACTER (len=25) :: charout
+      REAL(wp) :: zmax, zfact
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: zw3d
       REAL(wp), ALLOCATABLE, DIMENSION(:,:  ) :: zw2d
@@ -70,5 +68,4 @@
       !
       IF( ln_timing )   CALL timing_start('p4z_sink')
-
 
       ! Initialization of some global variables
@@ -97,50 +94,4 @@
 
       !
-      ! OA This is (I hope) a temporary solution for the problem that may 
-      ! OA arise in specific situation where the CFL criterion is broken 
-      ! OA for vertical sedimentation of particles. To avoid this, a time
-      ! OA splitting algorithm has been coded. A specific maximum
-      ! OA iteration number is provided and may be specified in the namelist 
-      ! OA This is to avoid very large iteration number when explicit free
-      ! OA surface is used (for instance). When niter?max is set to 1, 
-      ! OA this computation is skipped. The crude old threshold method is 
-      ! OA then applied. This also happens when niter exceeds nitermax.
-      IF( MAX( niter1max, niter2max ) == 1 ) THEN
-        iiter1 = 1
-        iiter2 = 1
-      ELSE
-        iiter1 = 1
-        iiter2 = 1
-        DO jk = 1, jpkm1
-          DO jj = 1, jpj
-             DO ji = 1, jpi
-                IF( tmask(ji,jj,jk) == 1) THEN
-                   zwsmax =  0.5 * e3t_n(ji,jj,jk) / xstep
-                   iiter1 =  MAX( iiter1, INT( wsbio3(ji,jj,jk) / zwsmax ) )
-                   iiter2 =  MAX( iiter2, INT( wsbio4(ji,jj,jk) / zwsmax ) )
-                ENDIF
-             END DO
-          END DO
-        END DO
-        IF( lk_mpp ) THEN
-           CALL mpp_max( 'p4zsink', iiter1 )
-           CALL mpp_max( 'p4zsink', iiter2 )
-        ENDIF
-        iiter1 = MIN( iiter1, niter1max )
-        iiter2 = MIN( iiter2, niter2max )
-      ENDIF
-
-      DO jk = 1,jpkm1
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( tmask(ji,jj,jk) == 1 ) THEN
-                 zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep
-                 wsbio3(ji,jj,jk) = MIN( wsbio3(ji,jj,jk), zwsmax * REAL( iiter1, wp ) )
-                 wsbio4(ji,jj,jk) = MIN( wsbio4(ji,jj,jk), zwsmax * REAL( iiter2, wp ) )
-               ENDIF
-            END DO
-         END DO
-      END DO
-
       !  Initializa to zero all the sinking arrays 
       !   -----------------------------------------
@@ -154,15 +105,10 @@
       !   Compute the sedimentation term using p4zsink2 for all the sinking particles
       !   -----------------------------------------------------
-      DO jit = 1, iiter1
-        CALL p4z_sink2( wsbio3, sinking , jppoc, iiter1 )
-        CALL p4z_sink2( wsbio3, sinkfer , jpsfe, iiter1 )
-      END DO
-
-      DO jit = 1, iiter2
-        CALL p4z_sink2( wsbio4, sinking2, jpgoc, iiter2 )
-        CALL p4z_sink2( wsbio4, sinkfer2, jpbfe, iiter2 )
-        CALL p4z_sink2( wsbio4, sinksil , jpgsi, iiter2 )
-        CALL p4z_sink2( wsbio4, sinkcal , jpcal, iiter2 )
-      END DO
+      CALL trc_sink( kt, wsbio3, sinking , jppoc, rfact2 )
+      CALL trc_sink( kt, wsbio3, sinkfer , jpsfe, rfact2 )
+      CALL trc_sink( kt, wsbio4, sinking2, jpgoc, rfact2 )
+      CALL trc_sink( kt, wsbio4, sinkfer2, jpbfe, rfact2 )
+      CALL trc_sink( kt, wsbio4, sinksil , jpgsi, rfact2 )
+      CALL trc_sink( kt, wsbio4, sinkcal , jpcal, rfact2 )
 
       IF( ln_p5z ) THEN
@@ -174,32 +120,15 @@
          !   Compute the sedimentation term using p4zsink2 for all the sinking particles
          !   -----------------------------------------------------
-         DO jit = 1, iiter1
-           CALL p4z_sink2( wsbio3, sinkingn , jppon, iiter1 )
-           CALL p4z_sink2( wsbio3, sinkingp , jppop, iiter1 )
-         END DO
-
-         DO jit = 1, iiter2
-           CALL p4z_sink2( wsbio4, sinking2n, jpgon, iiter2 )
-           CALL p4z_sink2( wsbio4, sinking2p, jpgop, iiter2 )
-         END DO
+         CALL trc_sink( kt, wsbio3, sinkingn , jppon, rfact2 )
+         CALL trc_sink( kt, wsbio3, sinkingp , jppop, rfact2 )
+         CALL trc_sink( kt, wsbio4, sinking2n, jpgon, rfact2 )
+         CALL trc_sink( kt, wsbio4, sinking2p, jpgop, rfact2 )
       ENDIF
 
       IF( ln_ligand ) THEN
          wsfep (:,:,:) = wfep
-         DO jk = 1,jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  IF( tmask(ji,jj,jk) == 1 ) THEN
-                    zwsmax = 0.5 * e3t_n(ji,jj,jk) / xstep
-                    wsfep(ji,jj,jk) = MIN( wsfep(ji,jj,jk), zwsmax * REAL( iiter1, wp ) )
-                  ENDIF
-               END DO
-            END DO
-         END DO
          !
          sinkfep(:,:,:) = 0.e0
-         DO jit = 1, iiter1
-           CALL p4z_sink2( wsfep, sinkfep , jpfep, iiter1 )
-         END DO
+         CALL trc_sink( kt, wsfep, sinkfep , jpfep, rfact2 )
       ENDIF
 
@@ -281,112 +210,4 @@
    END SUBROUTINE p4z_sink_init
 
-
-   SUBROUTINE p4z_sink2( pwsink, psinkflx, jp_tra, kiter )
-      !!---------------------------------------------------------------------
-      !!                     ***  ROUTINE p4z_sink2  ***
-      !!
-      !! ** Purpose :   Compute the sedimentation terms for the various sinking
-      !!     particles. The scheme used to compute the trends is based
-      !!     on MUSCL.
-      !!
-      !! ** Method  : - this ROUTINE compute not exactly the advection but the
-      !!      transport term, i.e.  div(u*tra).
-      !!---------------------------------------------------------------------
-      INTEGER , INTENT(in   )                         ::   jp_tra    ! tracer index index      
-      INTEGER , INTENT(in   )                         ::   kiter     ! number of iterations for time-splitting 
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed
-      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
-      !
-      INTEGER  ::   ji, jj, jk, jn
-      REAL(wp) ::   zigma,zew,zign, zflx, zstep
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz, zwsink2, ztrb 
-      !!---------------------------------------------------------------------
-      !
-      IF( ln_timing )   CALL timing_start('p4z_sink2')
-      !
-      zstep = rfact2 / REAL( kiter, wp ) / 2.
-
-      ztraz(:,:,:) = 0.e0
-      zakz (:,:,:) = 0.e0
-      ztrb (:,:,:) = trb(:,:,:,jp_tra)
-
-      DO jk = 1, jpkm1
-         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) 
-      END DO
-      zwsink2(:,:,1) = 0.e0
-
-
-      ! Vertical advective flux
-      DO jn = 1, 2
-         !  first guess of the slopes interior values
-         DO jk = 2, jpkm1
-            ztraz(:,:,jk) = ( trb(:,:,jk-1,jp_tra) - trb(:,:,jk,jp_tra) ) * tmask(:,:,jk)
-         END DO
-         ztraz(:,:,1  ) = 0.0
-         ztraz(:,:,jpk) = 0.0
-
-         ! slopes
-         DO jk = 2, jpkm1
-            DO jj = 1,jpj
-               DO ji = 1, jpi
-                  zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
-                  zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
-               END DO
-            END DO
-         END DO
-         
-         ! Slopes limitation
-         DO jk = 2, jpkm1
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) *        &
-                     &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
-               END DO
-            END DO
-         END DO
-         
-         ! vertical advective flux
-         DO jk = 1, jpkm1
-            DO jj = 1, jpj      
-               DO ji = 1, jpi    
-                  zigma = zwsink2(ji,jj,jk+1) * zstep / e3w_n(ji,jj,jk+1)
-                  zew   = zwsink2(ji,jj,jk+1)
-                  psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
-               END DO
-            END DO
-         END DO
-         !
-         ! Boundary conditions
-         psinkflx(:,:,1  ) = 0.e0
-         psinkflx(:,:,jpk) = 0.e0
-         
-         DO jk=1,jpkm1
-            DO jj = 1,jpj
-               DO ji = 1, jpi
-                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
-                  trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx
-               END DO
-            END DO
-         END DO
-
-      ENDDO
-
-      DO jk = 1,jpkm1
-         DO jj = 1,jpj
-            DO ji = 1, jpi
-               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
-               ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
-            END DO
-         END DO
-      END DO
-
-      trb(:,:,:,jp_tra) = ztrb(:,:,:)
-      psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
-      !
-      IF( ln_timing )  CALL timing_stop('p4z_sink2')
-      !
-   END SUBROUTINE p4z_sink2
-
-
    INTEGER FUNCTION p4z_sink_alloc()
       !!----------------------------------------------------------------------
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsms.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/P4Z/p4zsms.F90	(revision 10377)
@@ -74,5 +74,4 @@
             CALL p4z_che                              ! initialize the chemical constants
             CALL ahini_for_at(hi)   !  set PH at kt=nit000
-            t_oce_co2_flx_cum = 0._wp
         ELSE
             CALL p4z_rst( nittrc000, 'READ' )  !* read or initialize all required fields
@@ -189,5 +188,5 @@
       !!
       NAMELIST/nampisbio/ nrdttrc, wsbio, xkmort, ferat3, wsbio2, wsbio2max, wsbio2scale,    &
-         &                   niter1max, niter2max, wfep, ldocp, ldocz, lthet,  &
+         &                   wfep, ldocp, ldocz, lthet,  &
          &                   no3rat3, po4rat3
          !
@@ -223,6 +222,4 @@
          WRITE(numout,*) '      Big particles maximum sinking speed       wsbio2max   =', wsbio2max
          WRITE(numout,*) '      Big particles sinking speed length scale  wsbio2scale =', wsbio2scale
-         WRITE(numout,*) '      Maximum number of iterations for POC      niter1max   =', niter1max
-         WRITE(numout,*) '      Maximum number of iterations for GOC      niter2max   =', niter2max
          IF( ln_ligand ) THEN
             WRITE(numout,*) '      FeP sinking speed                              wfep   =', wfep
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/sms_pisces.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/sms_pisces.F90	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/PISCES/sms_pisces.F90	(revision 10377)
@@ -36,5 +36,4 @@
 
    !!*  Biological parameters 
-   INTEGER  ::   niter1max, niter2max   !: Maximum number of iterations for sinking
    REAL(wp) ::   rno3              !: ???
    REAL(wp) ::   o2ut              !: ???
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/TRP/trcsink.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/TRP/trcsink.F90	(revision 10377)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/TRP/trcsink.F90	(revision 10377)
@@ -0,0 +1,235 @@
+MODULE trcsink
+   !!======================================================================
+   !!                         ***  MODULE trcsink  ***
+   !! TOP :  vertical flux of particulate matter due to gravitational sinking
+   !!======================================================================
+   !! History :   1.0  !  2004     (O. Aumont) Original code
+   !!             2.0  !  2007-12  (C. Ethe, G. Madec)  F90
+   !!             3.4  !  2011-06  (O. Aumont, C. Ethe) Change aggregation formula
+   !!             3.5  !  2012-07  (O. Aumont) Introduce potential time-splitting
+   !!             4.0  !  2018-12  (O. Aumont) Generalize the PISCES code to make it usable by any model
+   !!----------------------------------------------------------------------
+   !!   trc_sink       :  Compute vertical flux of particulate matter due to gravitational sinking
+   !!----------------------------------------------------------------------
+   USE oce_trc         !  shared variables between ocean and passive tracers
+   USE trc             !  passive tracers common variables 
+   USE lib_mpp
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC trc_sink
+   PUBLIC trc_sink_ini
+
+   INTEGER, PUBLIC :: nitermax      !: Maximum number of iterations for sinking
+
+   !!----------------------------------------------------------------------
+   !! NEMO/TOP 4.0 , NEMO Consortium (2018)
+   !! $Id: trcsink.F90 10069 2018-08-28 14:12:24Z nicolasmartin $ 
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   !!----------------------------------------------------------------------
+   !!   'standard sinking parameterisation'                  ???
+   !!----------------------------------------------------------------------
+
+   SUBROUTINE trc_sink ( kt, pwsink, psinkflx, jp_tra, rsfact )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_sink  ***
+      !!
+      !! ** Purpose :   Compute vertical flux of particulate matter due to 
+      !!                gravitational sinking
+      !!
+      !! ** Method  : - ???
+      !!---------------------------------------------------------------------
+      INTEGER , INTENT(in)  :: kt
+      INTEGER , INTENT(in)  :: jp_tra    ! tracer index index      
+      REAL(wp), INTENT(in)  :: rsfact    ! time step duration
+      REAL(wp), INTENT(in)   , DIMENSION(jpi,jpj,jpk) :: pwsink
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) :: psinkflx
+      INTEGER  ::   ji, jj, jk
+      INTEGER, DIMENSION(jpi, jpj) ::   iiter
+      REAL(wp) ::   zfact, zwsmax, zmax
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zwsink
+      !!---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('trc_sink')
+      !
+      !
+      ! OA This is (I hope) a temporary solution for the problem that may 
+      ! OA arise in specific situation where the CFL criterion is broken 
+      ! OA for vertical sedimentation of particles. To avoid this, a time
+      ! OA splitting algorithm has been coded. A specific maximum
+      ! OA iteration number is provided and may be specified in the namelist 
+      ! OA This is to avoid very large iteration number when explicit free
+      ! OA surface is used (for instance). When niter?max is set to 1, 
+      ! OA this computation is skipped. The crude old threshold method is 
+      ! OA then applied. This also happens when niter exceeds nitermax.
+      IF( nitermax == 1 ) THEN
+         iiter(:,:) = 1
+      ELSE
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               iiter(ji,jj) = 1
+               DO jk = 1, jpkm1
+                  IF( tmask(ji,jj,jk) == 1.0 ) THEN
+                      zwsmax =  0.5 * e3t_n(ji,jj,jk) * rday / rsfact
+                      iiter(ji,jj) =  MAX( iiter(ji,jj), INT( pwsink(ji,jj,jk) / zwsmax ) )
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         iiter(ji,jj) = MIN( iiter(ji,jj), nitermax )
+      ENDIF
+
+      DO jk = 1,jpkm1
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( tmask(ji,jj,jk) == 1 ) THEN
+                 zwsmax = 0.5 * e3t_n(ji,jj,jk) * rday / rsfact
+                 zwsink(ji,jj,jk) = MIN( pwsink(ji,jj,jk), zwsmax * REAL( iiter(ji,jj), wp ) )
+               ENDIF
+            END DO
+         END DO
+      END DO
+
+      !  Initializa to zero all the sinking arrays 
+      !  -----------------------------------------
+      psinkflx(:,:,:) = 0.e0
+
+      !   Compute the sedimentation term using trc_sink2 for the considered sinking particle
+      !   -----------------------------------------------------
+      CALL trc_sink2( zwsink, psinkflx, jp_tra, iiter, rsfact )
+      !
+      IF( ln_timing )   CALL timing_stop('trc_sink')
+      !
+   END SUBROUTINE trc_sink
+
+   SUBROUTINE trc_sink2( pwsink, psinkflx, jp_tra, kiter, rsfact )
+      !!---------------------------------------------------------------------
+      !!                     ***  ROUTINE trc_sink2  ***
+      !!
+      !! ** Purpose :   Compute the sedimentation terms for the various sinking
+      !!     particles. The scheme used to compute the trends is based
+      !!     on MUSCL.
+      !!
+      !! ** Method  : - this ROUTINE compute not exactly the advection but the
+      !!      transport term, i.e.  div(u*tra).
+      !!---------------------------------------------------------------------
+      INTEGER,  INTENT(in   )                         ::   jp_tra    ! tracer index index      
+      REAL(wp), INTENT(in   )                         ::   rsfact    ! duration of time step
+      INTEGER,  INTENT(in   ), DIMENSION(jpi,jpj)     ::   kiter     ! number of iterations for time-splitting 
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj,jpk) ::   pwsink    ! sinking speed
+      REAL(wp), INTENT(inout), DIMENSION(jpi,jpj,jpk) ::   psinkflx  ! sinking fluxe
+      !
+      INTEGER  ::   ji, jj, jk, jn
+      REAL(wp) ::   zigma,zew,zign, zflx, zstep
+      REAL(wp), DIMENSION(jpi,jpj,jpk) :: ztraz, zakz, zwsink2, ztrb 
+      !!---------------------------------------------------------------------
+      !
+      IF( ln_timing )   CALL timing_start('trc_sink2')
+      !
+      ztraz(:,:,:) = 0.e0
+      zakz (:,:,:) = 0.e0
+      ztrb (:,:,:) = trb(:,:,:,jp_tra)
+
+      DO jk = 1, jpkm1
+         zwsink2(:,:,jk+1) = -pwsink(:,:,jk) / rday * tmask(:,:,jk+1) 
+      END DO
+      zwsink2(:,:,1) = 0.e0
+
+
+      ! Vertical advective flux
+      DO jn = 1, 2
+         !  first guess of the slopes interior values
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               !
+               zstep = rsfact / REAL( kiter(ji,jj), wp ) / 2.
+               !              
+               DO jk = 2, jpkm1
+                  ztraz(ji,jj,jk) = ( trb(ji,jj,jk-1,jp_tra) - trb(ji,jj,jk,jp_tra) ) * tmask(ji,jj,jk)
+               END DO
+               ztraz(ji,jj,1  ) = 0.0
+               ztraz(ji,jj,jpk) = 0.0
+
+               ! slopes
+               DO jk = 2, jpkm1
+                  zign = 0.25 + SIGN( 0.25, ztraz(ji,jj,jk) * ztraz(ji,jj,jk+1) )
+                  zakz(ji,jj,jk) = ( ztraz(ji,jj,jk) + ztraz(ji,jj,jk+1) ) * zign
+               END DO
+         
+               ! Slopes limitation
+               DO jk = 2, jpkm1
+                  zakz(ji,jj,jk) = SIGN( 1., zakz(ji,jj,jk) ) *        &
+                     &             MIN( ABS( zakz(ji,jj,jk) ), 2. * ABS(ztraz(ji,jj,jk+1)), 2. * ABS(ztraz(ji,jj,jk) ) )
+               END DO
+         
+               ! vertical advective flux
+               DO jk = 1, jpkm1
+                  zigma = zwsink2(ji,jj,jk+1) * zstep / e3w_n(ji,jj,jk+1)
+                  zew   = zwsink2(ji,jj,jk+1)
+                  psinkflx(ji,jj,jk+1) = -zew * ( trb(ji,jj,jk,jp_tra) - 0.5 * ( 1 + zigma ) * zakz(ji,jj,jk) ) * zstep
+               END DO
+               !
+               ! Boundary conditions
+               psinkflx(ji,jj,1  ) = 0.e0
+               psinkflx(ji,jj,jpk) = 0.e0
+         
+               DO jk=1,jpkm1
+                  zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
+                  trb(ji,jj,jk,jp_tra) = trb(ji,jj,jk,jp_tra) + zflx
+               END DO
+            END DO
+         END DO
+      END DO
+
+      DO jk = 1,jpkm1
+         DO jj = 1,jpj
+            DO ji = 1, jpi
+               zflx = ( psinkflx(ji,jj,jk) - psinkflx(ji,jj,jk+1) ) / e3t_n(ji,jj,jk)
+               ztrb(ji,jj,jk) = ztrb(ji,jj,jk) + 2. * zflx
+            END DO
+         END DO
+      END DO
+
+      trb(:,:,:,jp_tra) = ztrb(:,:,:)
+      psinkflx(:,:,:)   = 2. * psinkflx(:,:,:)
+      !
+      IF( ln_timing )  CALL timing_stop('trc_sink2')
+      !
+   END SUBROUTINE trc_sink2
+
+  SUBROUTINE trc_sink_ini
+      !!---------------------------------------------------------------------
+      !!                  ***  ROUTINE trc_sink_ini ***
+      !!
+      !! ** Purpose :   read  namelist options 
+      !!----------------------------------------------------------------------
+      INTEGER ::   ios   ! Local integer output status for namelist read
+      !!
+      NAMELIST/namtrc_snk/ nitermax
+      !!----------------------------------------------------------------------
+      !
+      REWIND( numnat_ref )              ! namtrc_rad in reference namelist 
+      READ  ( numnat_ref, namtrc_snk, IOSTAT = ios, ERR = 907)
+907   IF( ios /= 0 )   CALL ctl_nam ( ios , 'namtrc_snk in reference namelist', lwp )
+      REWIND( numnat_cfg )              ! namtrc_rad in configuration namelist 
+      READ  ( numnat_cfg, namtrc_snk, IOSTAT = ios, ERR = 908 )
+908   IF( ios > 0 )   CALL ctl_nam ( ios , 'namtrc_snk in configuration namelist', lwp )
+      IF(lwm) WRITE( numont, namtrc_snk )
+
+      IF(lwp) THEN                     !   ! Control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'trc_sink : Sedimentation of particles '
+         WRITE(numout,*) '~~~~~~~ '
+         WRITE(numout,*) '   Namelist namtrc_snk : sedimentation of particles'
+         WRITE(numout,*) '   Maximum number of iterations   nitermax = ', nitermax
+         WRITE(numout,*)
+      ENDIF
+
+   END SUBROUTINE trc_sink_ini
+
+   !!======================================================================
+END MODULE trcsink
Index: NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/trcini.F90
===================================================================
--- NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/trcini.F90	(revision 10372)
+++ NEMO/branches/2018/dev_r10164_HPC09_ESIWACE_PREP_MERGE/src/TOP/trcini.F90	(revision 10377)
@@ -205,4 +205,5 @@
       USE trcldf , ONLY:  trc_ldf_ini
       USE trcrad , ONLY:  trc_rad_ini
+      USE trcsink, ONLY:  trc_sink_ini
       !
       INTEGER :: ierr
@@ -214,4 +215,5 @@
                        !                          ! vertical diffusion: always implicit time stepping scheme
                        CALL  trc_rad_ini          ! positivity of passive tracers 
+                       CALL  trc_sink_ini         ! Vertical sedimentation of particles
       !
    END SUBROUTINE trc_ini_trp
