diff --git a/latex/NEMO/subfiles/chap_D2D.tex b/latex/NEMO/subfiles/chap_D2D.tex
index f984936ec7fe9d471ac9142fbeb20538d4be6f26..c493271c179e9454ff25a694bd207b9b63cbe945 100644
--- a/latex/NEMO/subfiles/chap_D2D.tex
+++ b/latex/NEMO/subfiles/chap_D2D.tex
@@ -2,7 +2,7 @@
 
 \begin{document}
 
-\chapter{Ocean Dynamics (D2D)}
+\chapter{Sea surface height and 2D external mode (D2D)}
 \label{chap:D2D}
 
 \chaptertoc
@@ -22,83 +22,76 @@
 
 \clearpage
 
-Using the representation described in \autoref{chap:DOM},
-several semi-discrete space forms of the dynamical equations are available depending on
-the vertical coordinate used and on the conservation properties of the vorticity term.
-In all the equations presented here, the masking has been omitted for simplicity.
-One must be aware that all the quantities are masked fields and
-that each time an average or difference operator is used, the resulting field is multiplied by a mask.
-
-The prognostic ocean dynamics equation can be summarized as follows:
-\[
-  \text{NXT} = \dbinom	{\text{VOR} + \text{KEG} + \text {ZAD} }
-  {\text{COR} + \text{ADV}                       }
-  + \text{HPG} + \text{SPG} + \text{LDF} + \text{ZDF}
-\]
-NXT stands for next, referring to the time-stepping.
-The first group of terms on the rhs of this equation corresponds to the Coriolis and advection terms that
-are decomposed into either a vorticity part (VOR), a kinetic energy part (KEG) and
-a vertical advection part (ZAD) in the vector invariant formulation,
-or a Coriolis and advection part (COR+ADV) in the flux formulation.
-The terms following these are the pressure gradient contributions
-(HPG, Hydrostatic Pressure Gradient, and SPG, Surface Pressure Gradient);
-and contributions from lateral diffusion (LDF) and vertical diffusion (ZDF),
-which are added to the rhs in the \mdl{dynldf} and \mdl{dynzdf} modules.
-The vertical diffusion term includes the surface and bottom stresses.
-The external forcings and parameterisations require complex inputs
-(surface wind stress calculation using bulk formulae, estimation of mixing coefficients)
-that are carried out in modules SBC, LDF and ZDF and are described in
-\autoref{chap:SBC}, \autoref{chap:LDF} and \autoref{chap:ZDF}, respectively.
-
-In the present chapter we describe the equations used to compute the barotropic velocities and
-the sea surface height that provides the surface pressure gradient.
-
-%% =================================================================================================
-\section[Surface pressure gradient (\textit{dynspg.F90}/\textit{stp2d.F90})]
-{Surface pressure gradient (\protect\mdl{dynspg}/\protect\mdl{stp2d})}
-\label{sec:D2D_spg}
-
-The surface pressure gradient term is related to the representation of the free surface
-(\autoref{sec:MB_hor_pg}).  The main distinction is between the fixed volume case (linear
-free surface, \ie\ with the inclusion of \key{linssh}) and the variable volume case
-(nonlinear free surface, \key{qco}).  In the linear free surface case
-(\autoref{subsec:MB_free_surface}) the vertical scale factors $e_{3}$ are fixed in time,
-while they are time-dependent in the nonlinear case (\autoref{subsec:MB_free_surface}).
-With both linear and nonlinear free surface, external gravity waves are allowed in the
-equations, which imposes a very small time step when an explicit time stepping is used
-(\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp} (MLF time-stepping only)).  To allow a
-longer time step for the three-dimensional equations, one can use a split-explicit free
-surface (\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}). In that case, only a quasi-linear
- form of 2d barotropic equations is substepped with a small time increment.
-Options are defined through the \nam{dyn_spg}{dyn\_spg} namelist variables.
-
 \begin{listing}
   \nlst{namdyn_spg}
   \caption{\forcode{&namdyn_spg}}
   \label{lst:namdyn_spg}
 \end{listing}
 
+In the present chapter we describe the equations used to compute the the sea surface
+height. Recall from \autoref{subsec:MB_free_surface} that allowing a free surface permits
+the existence of External Gravity Waves (EGWs). Resolving these fast moving waves explicitly
+imposes a severe limit on the model timestep to avoid breaching the CFL condition.
+Alternatively, EGWs can be filtered by discretisation of the temporal derivatives using a
+split-expilcit scheme.  This chapter describes these options controlled via the
+\nam{dyn_spg}{dyn_spg} namelist.
+
+%% =================================================================================================
+\section[Sea surface height  (\textit{sshwzv.F90}/\textit{stp2d.F90})]
+{Sea surface height (\protect\mdl{sshwzv}/\protect\mdl{stp2d})}
+\label{sec:D2D_ssh}
+
+The sea surface height is given by the vertical average of the kinematic surface condition 
+(\autoref{subsec:MB_free_surface}):
+\begin{equation}
+  \label{eq:D2D_spg_ssh}
+  \begin{aligned}
+    \frac{\partial \eta }{\partial t}
+    &\equiv    \frac{1}{e_{1t} e_{2t} }\sum\limits_k { \left\{  \delta_i \left[ {e_{2u}\,e_{3u}\;u} \right]
+        +\delta_j \left[ {e_{1v}\,e_{3v}\;v} \right]  \right\} }
+    -    \frac{\textit{emp}}{\rho_w }   \\
+    &\equiv    \sum\limits_k {\chi \ e_{3t}}  -  \frac{\textit{emp}}{\rho_w }
+  \end{aligned}
+\end{equation}
+
+where \textit{emp} is the surface freshwater budget (i.e. evaporation minus precipitation
++ optional contributions (see \autoref{subsec:D2D_forcing})), expressed in Kg/m$^2$/s (which is
+equal to mm/s), and $\rho_w$=1,026~Kg/m$^3$ is the reference density of sea water
+(Boussinesq approximation).  The sea-surface height is evaluated using exactly the same
+time stepping scheme as the tracer \autoref{eq:TRA_nxt}.  E.g., in the case of
+MLF this is a leapfrog scheme in combination with an Asselin time filter, \ie\ the
+velocity appearing in \autoref{eq:D2D_spg_ssh} is centred in time (\textit{now} velocity).
+This is of paramount importance.  Replacing $T$ by the number $1$ in the tracer equation
+and summing over the water column must lead to the sea surface height equation otherwise
+tracer content will not be conserved \citep{griffies.pacanowski.ea_MWR01,
+leclair.madec_OM09}.
+
 %% =================================================================================================
-\subsection[Explicit free surface (\forcode{ln_dynspg_exp})]
+\subsection[Explicit free surface (\forcode{ln_dynspg_exp}, \textit{sshwzv.F90})]
 {Explicit free surface (\protect\np{ln_dynspg_exp}{ln\_dynspg\_exp})}
 \label{subsec:D2D_spg_exp}
 
 With MLF time-stepping there is an explicit free surface formulation
 (\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp}), for which the model time step must be
 chosen to be small enough to resolve the external gravity waves (typically a few tens of
-seconds).  The surface pressure gradient, evaluated using a leap-frog scheme (\ie\
-centered in time), is thus simply given by :
-\begin{equation}
-  \label{eq:DYN_spg_exp}
-  \left\{
-    \begin{aligned}
-      - \frac{1}{e_{1u}\,\rho_o} \;	\delta_{i+1/2} \left[  \,\rho \,\eta\,  \right] 	\\
-      - \frac{1}{e_{2v}\,\rho_o} \;	\delta_{j+1/2} \left[  \,\rho \,\eta\,  \right]
-    \end{aligned}
-  \right.
-\end{equation}
-
-Note that this option is not yet available with RK3 time-stepping.
+seconds).  The sea surface height is evaluated from \autoref{eq:D2D_spg_ssh}
+using a leap-frog scheme (\ie\ centered in time) with ssh at the now timestep coming from
+the previous time step computation.
+
+%, is thus simply given by :
+%\begin{equation}
+%  \label{eq:DYN_spg_exp}
+%  \left\{
+%    \begin{aligned}
+%      - \frac{1}{e_{1u}\,\rho_o} \;	\delta_{i+1/2} \left[  \,\rho \,\eta\,  \right] 	\\
+%      - \frac{1}{e_{2v}\,\rho_o} \;	\delta_{j+1/2} \left[  \,\rho \,\eta\,  \right]
+%    \end{aligned}
+%  \right.
+%\end{equation}
+
+Note that this option is not yet available with RK3 time-stepping but, if implemented, ssh 
+would be computed at each stage from \autoref{eq:D2D_spg_ssh} using values at $n$,
+$n+{1\over3}$ and $n+{1\over2}$ respectively.
 
 %% =================================================================================================
 \subsection[Split-explicit free surface (\forcode{ln_dynspg_ts})]{Split-explicit free surface (\protect\np{ln_dynspg_ts}{ln\_dynspg\_ts})}
@@ -136,7 +129,7 @@ The barotropic mode solves the following equations:
 where $\mathrm {\overline{\mathbf G}}$ is a forcing term held constant, containing
 coupling term between modes, surface atmospheric forcing as well as slowly varying
 barotropic terms not explicitly computed to gain efficiency.  The third term on the right
-hand side of \autoref{eq:D2D_BT_dyn} represents the bottom stress (see section
+hand side of \autoref{eq:D2D_BT_dyn} represents the bottom stress (see 
 \autoref{sec:ZDF_drg}), explicitly accounted for at each barotropic iteration.  
 
 For 'single 1$^{\mathrm st}$', RK3 time-stepping, barotropic velocities and sea surface
@@ -151,7 +144,7 @@ For clarity, in RK3 the RHS forcing terms on momentum are vertically averaged 3D
 \[
   \text{HPG} + \text{LDF} + \text {(COR + RVO)} + \text{KEG} + \text{ZAD} 
 \]
-for the Vector invariant formulation (for which the 3D fields are also the RHS terms for
+for the Vector invariant form (for which the 3D fields are also the RHS terms for
 the 1$^{\mathrm st}$ stage RK3 time stepping) and:
 \[
   \text{HPG} + \text{LDF} + \text {(COR + MET)} + \text{ADV}
@@ -179,7 +172,7 @@ defined in \citet{shchepetkin.mcwilliams_ibk09} (see their figure 12, lower left
   \centering
   \includegraphics[width=0.66\textwidth]{DYN_dynspg_ts}
   \caption[Split-explicit time stepping scheme for the external and internal modes]{
-    Schematic of the split-explicit time stepping scheme for the external and internal modes.
+    Schematic of the split-explicit time stepping scheme for the external and internal modes with MLF.
     Time increases to the right.
     In this particular exemple,
     a boxcar averaging window over \np{nn_e}{nn\_e} barotropic time steps is used
@@ -244,26 +237,6 @@ diffusion.
 \subsection{External forcings}
 \label{subsec:D2D_forcing}
 
-Besides the surface and bottom baroclinic stresses four other forcings may enter the barotropic
-momentum equations.
-
-(1) When \np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} (see \autoref{sec:SBC_apr}),
-the atmospheric pressure is taken into account when computing the surface pressure gradient.
-
-(2) When \np[=.true.]{ln_tide_pot}{ln\_tide\_pot} and \np[=.true.]{ln_tide}{ln\_tide} (see
-\autoref{sec:SBC_TDE}), the tidal potential is taken into account when computing the
-surface pressure gradient.
-
-(3) When \np[=2]{nn_ice_embd}{nn\_ice\_embd} and SI3 is used (\ie\ when the sea-ice is
-embedded in the ocean), the snow-ice mass is taken into account when computing the surface
-pressure gradient.
-
-(4) When \np[=.true.]{ln_wave .AND. ln_bern_srfc}{ln\_wave .AND. ln\_bern\_srfc} (see
-\autoref{subsec:SBC_wave_bhd}) , a depth-uniform wave-induced kinematic pressure term (the
-Bernoulli head) is added to the mean pressure.
-
-\vskip 0.5cm
-
 \noindent The net column-average freshwater flux applied to the ssh equation can also be modified by
 the following options:
 
@@ -286,7 +259,25 @@ computing the net freshwater flux.
 ln\_sshinc .AND. ln\_asmiau} (see \autoref{sec:ASM_IAU}), the contribution from IAU
 weighted ssh increments is taken into account when computing the net freshwater flux.
 
+\vskip 0.5cm
+
+\noindent Besides the surface and bottom baroclinic stresses four other forcings may enter the barotropic
+momentum equations.
+
+(1) When \np[=.true.]{ln_apr_dyn}{ln\_apr\_dyn} (see \autoref{sec:SBC_apr}),
+the atmospheric pressure is taken into account when computing the surface pressure gradient.
 
+(2) When \np[=.true.]{ln_tide_pot}{ln\_tide\_pot} and \np[=.true.]{ln_tide}{ln\_tide} (see
+\autoref{sec:SBC_TDE}), the tidal potential is taken into account when computing the
+surface pressure gradient.
+
+(3) When \np[=2]{nn_ice_embd}{nn\_ice\_embd} and SI3 is used (\ie\ when the sea-ice is
+embedded in the ocean), the snow-ice mass is taken into account when computing the surface
+pressure gradient.
+
+(4) When \np[=.true.]{ln_wave .AND. ln_bern_srfc}{ln\_wave .AND. ln\_bern\_srfc} (see
+\autoref{subsec:SBC_wave_bhd}) , a depth-uniform wave-induced kinematic pressure term (the
+Bernoulli head) is added to the mean pressure.
 
 \cmtgm{ missing : the lateral boundary condition !!!   another external forcing
  }
diff --git a/latex/NEMO/subfiles/chap_DYN.tex b/latex/NEMO/subfiles/chap_DYN.tex
index 591dd3c37f0850590d72946114a7943f6c19689e..2f3be49dd9a6311b3a8916e966c71871c38a8e45 100644
--- a/latex/NEMO/subfiles/chap_DYN.tex
+++ b/latex/NEMO/subfiles/chap_DYN.tex
@@ -158,14 +158,15 @@ This option can be useful when the value of the timestep is limited by vertical
   \label{lst:namdyn_adv}
 \end{listing}
 
-The vector invariant form of the momentum equation is most commonly used in coarse-resolution (1°) applications of the \NEMO\ ocean model.
-For higher resolutions it requires the activation of Hollingsworth correction (\np[=1]{nn_dynkeg}{nn\_dynkeg}) following Arakawa (2001) 
-The flux form option (see next section) has been present since version $2$.
-By structuring the equations in vector invariant form, the dynamics are expressed in terms of
-intrinsic geometric properties like gradients, curls, and divergences. 
-This ensures that the physics remain consistent and interpretable regardless of 
-the underlying curvilinear grid or coordinate system.
-It highlights key physical terms like the kinetic energy advection and the relative vorticity. 
+The vector invariant form of the momentum equation is most commonly used in
+coarse-resolution (1°) applications of the \NEMO\ ocean model.  For higher resolutions it
+requires the activation of Hollingsworth correction (\np[=1]{nn_dynkeg}{nn\_dynkeg})
+following Arakawa (2001) The flux form option (see next section) has been present since
+version $2$.  By structuring the equations in vector invariant form, the dynamics are
+expressed in terms of intrinsic geometric properties like gradients, curls, and
+divergences.  This ensures that the physics remain consistent and interpretable regardless
+of the underlying curvilinear grid or coordinate system.  It highlights key physical terms
+like the kinetic energy advection and the relative vorticity.
 
 \vskip 0.5cm
 
@@ -766,6 +767,40 @@ The pressure gradient due to ocean load is computed using the expression \autore
 %% an extra three-dimensional field, in the restart file to restart the model with exact reproducibility.
 %% This option is controlled by  \np{nn_dynhpg_rst}{nn\_dynhpg\_rst}, a namelist parameter.
 %% 
+\section[Surface pressure gradient (\textit{dynspg.F90}/\textit{stp2d.F90})]
+{Surface pressure gradient (\protect\mdl{dynspg}/\protect\mdl{stp2d})}
+\label{sec:DYN_spg}
+
+The surface pressure gradient term is related to the representation of the free surface
+(\autoref{sec:MB_hor_pg}).  The main distinction is between the fixed volume case (linear
+free surface, \ie\ with the inclusion of \key{linssh}) and the variable volume case
+(nonlinear free surface, \key{qco}).  In the linear free surface case
+(\autoref{subsec:MB_free_surface}) the vertical scale factors $e_{3}$ are fixed in time,
+while they are time-dependent in the nonlinear case (\autoref{subsec:MB_free_surface}).
+With both linear and nonlinear free surface, external gravity waves are allowed in the
+equations, which imposes a very small time step when an explicit time stepping is used
+(\np[=.true.]{ln_dynspg_exp}{ln\_dynspg\_exp} (MLF time-stepping only)). With explicit 
+time-stepping, the surface pressure gradient is evaluated using the leap-frog scheme 
+(\ie centred in time) and is thus simply given by:  
+
+\begin{equation}
+  \label{eq:DYN_spg_exp}
+  \left\{
+    \begin{aligned}
+      - \frac{1}{e_{1u}\,\rho_o} \;    \delta_{i+1/2} \left[  \,\rho \,\eta\,  \right]         \\
+      - \frac{1}{e_{2v}\,\rho_o} \;    \delta_{j+1/2} \left[  \,\rho \,\eta\,  \right]
+    \end{aligned}
+  \right.
+\end{equation}
+
+where values are evaluated at the now timestep (\textit{dynspg\_exp.F90}). 
+
+To allow a longer time step for the three-dimensional equations, one can use a
+split-explicit free surface (\np[=.true.]{ln_dynspg_ts}{ln\_dynspg\_ts}). In that case,
+a quasi-linear form of 2d barotropic equations is substepped with a small time
+increment. Details of this were provided in \autoref{chap:D2D}.  Options are defined
+through the \nam{dyn_spg}{dyn\_spg} namelist variables.
+
 %% =================================================================================================
 \section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})}
 \label{sec:DYN_ldf}