chap_DYN.tex

  \eta(\tau+\Delta) - \eta^{F}(\tau-\Delta) = 2\rdt \ \left[ - \nabla \cdot \textbf{U}(\tau) + \text{EMP}_w \right]
\end{equation}

The use of this "big-leap-frog" scheme for the surface height ensures compatibility between
the mass/volume budgets and the tracer budgets.
More discussion of this point is provided in Chapter 10 (see in particular Section 10.2).

In general, some form of time filter is needed to maintain integrity of the surface height field due to
the leap-frog splitting mode in equation \autoref{eq:DYN_spg_ts_ssh}.
We have tried various forms of such filtering,
with the following method discussed in \cite{griffies.pacanowski.ea_MWR01} chosen due to
its stability and reasonably good maintenance of tracer conservation properties (see ??).

\begin{equation}
  \label{eq:DYN_spg_ts_sshf}
  \eta^{F}(\tau-\Delta) =  \overline{\eta^{(b)}(\tau)}
\end{equation}
Another approach tried was

\[
  % \label{eq:DYN_spg_ts_sshf2}
  \eta^{F}(\tau-\Delta) = \eta(\tau)
  + (\alpha/2) \left[\overline{\eta^{(b)}}(\tau+\rdt)
    + \overline{\eta^{(b)}}(\tau-\rdt) -2 \;\eta(\tau) \right]
\]

which is useful since it isolates all the time filtering aspects into the term multiplied by $\alpha$.
This isolation allows for an easy check that tracer conservation is exact when
eliminating tracer and surface height time filtering (see ?? for more complete discussion).
However, in the general case with a non-zero $\alpha$,
the filter \autoref{eq:DYN_spg_ts_sshf} was found to be more conservative, and so is recommended.

}            %%end gm comment (copy of griffies book)

%% =================================================================================================
\section[Lateral diffusion term and operators (\textit{dynldf.F90})]{Lateral diffusion term and operators (\protect\mdl{dynldf})}
\label{sec:DYN_ldf}

\begin{listing}
  \nlst{namdyn_ldf}
  \caption{\forcode{&namdyn_ldf}}
  \label{lst:namdyn_ldf}
\end{listing}

Options are defined through the \nam{dyn_ldf}{dyn\_ldf} namelist variables.
The options available for lateral diffusion are to use either laplacian (rotated or not) or biharmonic operators.
The coefficients may be constant or spatially variable;
the description of the coefficients is found in the chapter on lateral physics (\autoref{chap:LDF}).
The lateral diffusion of momentum is evaluated using a forward scheme,
\ie\ the velocity appearing in its expression is the \textit{before} velocity in time,
except for the pure vertical component that appears when a tensor of rotation is used.
This latter term is solved implicitly together with the vertical diffusion term (see \autoref{chap:TD}).

At the lateral boundaries either free slip,
no slip or partial slip boundary conditions are applied according to the user's choice (see \autoref{chap:LBC}).

\cmtgm{
  Hyperviscous operators are frequently used in the simulation of turbulent flows to
  control the dissipation of unresolved small scale features.
  Their primary role is to provide strong dissipation at the smallest scale supported by
  the grid while minimizing the impact on the larger scale features.
  Hyperviscous operators are thus designed to be more scale selective than the traditional,
  physically motivated Laplace operator.
  In finite difference methods,
  the biharmonic operator is frequently the method of choice to achieve this scale selective dissipation since
  its damping time (\ie\ its spin down time) scale like $\lambda^{-4}$ for disturbances of wavelength $\lambda$
  (so that short waves damped more rapidelly than long ones),
  whereas the Laplace operator damping time scales only like $\lambda^{-2}$.
}

%% =================================================================================================
\subsection[Iso-level laplacian (\forcode{ln_dynldf_lap})]{Iso-level laplacian operator (\protect\np{ln_dynldf_lap}{ln\_dynldf\_lap})}
\label{subsec:DYN_ldf_lap}

For lateral iso-level diffusion, the discrete operator is:
\begin{equation}
  \label{eq:DYN_ldf_lap}
  \left\{
    \begin{aligned}
      D_u^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{1u} }\delta_{i+1/2} \left[ {A_T^{lm}
          \;\chi } \right]-\frac{1}{e_{2u} {\kern 1pt}e_{3u} }\delta_j \left[
        {A_f^{lm} \;e_{3f} \zeta } \right] \\ \\
      D_v^{l{\mathrm {\mathbf U}}} =\frac{1}{e_{2v} }\delta_{j+1/2} \left[ {A_T^{lm}
          \;\chi } \right]+\frac{1}{e_{1v} {\kern 1pt}e_{3v} }\delta_i \left[
        {A_f^{lm} \;e_{3f} \zeta } \right]
    \end{aligned}
  \right.
\end{equation}

As explained in \autoref{subsec:MB_ldf},
this formulation (as the gradient of a divergence and curl of the vorticity) preserves symmetry and
ensures a complete separation between the vorticity and divergence parts of the momentum diffusion.

%% =================================================================================================
\subsection[Rotated laplacian (\forcode{ln_dynldf_iso})]{Rotated laplacian operator (\protect\np{ln_dynldf_iso}{ln\_dynldf\_iso})}
\label{subsec:DYN_ldf_iso}

A rotation of the lateral momentum diffusion operator is needed in several cases:
for iso-neutral diffusion in the $z$-coordinate (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) and
for either iso-neutral (\np[=.true.]{ln_dynldf_iso}{ln\_dynldf\_iso}) or
geopotential (\np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}) diffusion in the $s$-coordinate.
In the partial step case, coordinates are horizontal except at the deepest level and
no rotation is performed when \np[=.true.]{ln_dynldf_hor}{ln\_dynldf\_hor}.
The diffusion operator is defined simply as the divergence of down gradient momentum fluxes on
each momentum component.
It must be emphasized that this formulation ignores constraints on the stress tensor such as symmetry.
The resulting discrete representation is:
\begin{equation}
  \label{eq:DYN_ldf_iso}
  \begin{split}
    D_u^{l\textbf{U}} &= \frac{1}{e_{1u} \, e_{2u} \, e_{3u} }	\\
    &  \left\{\quad  {\delta_{i+1/2} \left[ {A_T^{lm}  \left(
              {\frac{e_{2t} \; e_{3t} }{e_{1t} } \,\delta_{i}[u]
                -e_{2t} \; r_{1t} \,\overline{\overline {\delta_{k+1/2}[u]}}^{\,i,\,k}}
            \right)} \right]} 	\right. \\
    & \qquad +\ \delta_j \left[ {A_f^{lm} \left( {\frac{e_{1f}\,e_{3f} }{e_{2f}
            }\,\delta_{j+1/2} [u] - e_{1f}\, r_{2f}
            \,\overline{\overline {\delta_{k+1/2} [u]}} ^{\,j+1/2,\,k}}
        \right)} \right] \\
    &\qquad +\ \delta_k \left[ {A_{uw}^{lm} \left( {-e_{2u} \, r_{1uw} \,\overline{\overline
              {\delta_{i+1/2} [u]}}^{\,i+1/2,\,k+1/2} }
        \right.} \right. \\
    &  \ \qquad \qquad \qquad \quad\
    - e_{1u} \, r_{2uw} \,\overline{\overline {\delta_{j+1/2} [u]}} ^{\,j,\,k+1/2} \\
    & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\
                +\frac{e_{1u}\, e_{2u} }{e_{3uw} }\,\left( {r_{1uw}^2+r_{2uw}^2}
                \right)\,\delta_{k+1/2} [u]} \right)} \right]\;\;\;} \right\} \\ \\
    D_v^{l\textbf{V}} &= \frac{1}{e_{1v} \, e_{2v} \, e_{3v} } \\
    &  \left\{\quad  {\delta_{i+1/2} \left[ {A_f^{lm}  \left(
              {\frac{e_{2f} \; e_{3f} }{e_{1f} } \,\delta_{i+1/2}[v]
                -e_{2f} \; r_{1f} \,\overline{\overline {\delta_{k+1/2}[v]}}^{\,i+1/2,\,k}}
            \right)} \right]} 	\right. \\
    & \qquad +\ \delta_j \left[ {A_T^{lm} \left( {\frac{e_{1t}\,e_{3t} }{e_{2t}
            }\,\delta_{j} [v] - e_{1t}\, r_{2t}
            \,\overline{\overline {\delta_{k+1/2} [v]}} ^{\,j,\,k}}
        \right)} \right] \\
    & \qquad +\ \delta_k \left[ {A_{vw}^{lm} \left( {-e_{2v} \, r_{1vw} \,\overline{\overline
              {\delta_{i+1/2} [v]}}^{\,i+1/2,\,k+1/2} }\right.} \right. \\
    &  \ \qquad \qquad \qquad \quad\
    - e_{1v} \, r_{2vw} \,\overline{\overline {\delta_{j+1/2} [v]}} ^{\,j+1/2,\,k+1/2} \\
    & \left. {\left. { \ \qquad \qquad \qquad \ \ \ \left. {\
                +\frac{e_{1v}\, e_{2v} }{e_{3vw} }\,\left( {r_{1vw}^2+r_{2vw}^2}
                \right)\,\delta_{k+1/2} [v]} \right)} \right]\;\;\;} \right\}
  \end{split}
\end{equation}
where $r_1$ and $r_2$ are the slopes between the surface along which the diffusion operator acts and
the surface of computation ($z$- or $s$-surfaces).
The way these slopes are evaluated is given in the lateral physics chapter (\autoref{chap:LDF}).

%% =================================================================================================
\subsection[Iso-level bilaplacian (\forcode{ln_dynldf_bilap})]{Iso-level bilaplacian operator (\protect\np{ln_dynldf_bilap}{ln\_dynldf\_bilap})}
\label{subsec:DYN_ldf_bilap}

The lateral fourth order operator formulation on momentum is obtained by applying \autoref{eq:DYN_ldf_lap} twice.
It requires an additional assumption on boundary conditions:
the first derivative term normal to the coast depends on the free or no-slip lateral boundary conditions chosen,
while the third derivative terms normal to the coast are set to zero (see \autoref{chap:LBC}).
\cmtgm{add a remark on the the change in the position of the coefficient}

%% =================================================================================================
\section[Vertical diffusion term (\textit{dynzdf.F90})]{Vertical diffusion term (\protect\mdl{dynzdf})}
\label{sec:DYN_zdf}

Options are defined through the \nam{zdf}{zdf} namelist variables.
The large vertical diffusion coefficient found in the surface mixed layer together with high vertical resolution implies that in the case of explicit time stepping there would be too restrictive a constraint on the time step.
Two time stepping schemes can be used for the vertical diffusion term:
$(a)$ a forward time differencing scheme
(\np[=.true.]{ln_zdfexp}{ln\_zdfexp}) using a time splitting technique (\np{nn_zdfexp}{nn\_zdfexp} $>$ 1) or
$(b)$ a backward (or implicit) time differencing scheme (\np[=.false.]{ln_zdfexp}{ln\_zdfexp})
(see \autoref{chap:TD}).
Note that namelist variables \np{ln_zdfexp}{ln\_zdfexp} and \np{nn_zdfexp}{nn\_zdfexp} apply to both tracers and dynamics.

The formulation of the vertical subgrid scale physics is the same whatever the vertical coordinate is.
The vertical diffusion operators given by \autoref{eq:MB_zdf} take the following semi-discrete space form:
\[
  % \label{eq:DYN_zdf}
  \left\{
    \begin{aligned}
      D_u^{vm} &\equiv \frac{1}{e_{3u}} \ \delta_k \left[ \frac{A_{uw}^{vm} }{e_{3uw} }
        \ \delta_{k+1/2} [\,u\,]         \right]     \\
      \\
      D_v^{vm} &\equiv \frac{1}{e_{3v}} \ \delta_k \left[ \frac{A_{vw}^{vm} }{e_{3vw} }
        \ \delta_{k+1/2} [\,v\,]         \right]
    \end{aligned}
  \right.
\]
where $A_{uw}^{vm} $ and $A_{vw}^{vm} $ are the vertical eddy viscosity and diffusivity coefficients.
The way these coefficients are evaluated depends on the vertical physics used (see \autoref{chap:ZDF}).

The surface boundary condition on momentum is the stress exerted by the wind.
At the surface, the momentum fluxes are prescribed as the boundary condition on
the vertical turbulent momentum fluxes,
\begin{equation}
  \label{eq:DYN_zdf_sbc}
  \left.{\left( {\frac{A^{vm} }{e_3 }\ \frac{\partial \textbf{U}_h}{\partial k}} \right)} \right|_{z=1}
  = \frac{1}{\rho_o} \binom{\tau_u}{\tau_v }
\end{equation}
where $\left( \tau_u ,\tau_v \right)$ are the two components of the wind stress vector in
the (\textbf{i},\textbf{j}) coordinate system.
The high mixing coefficients in the surface mixed layer ensure that the surface wind stress is distributed in
the vertical over the mixed layer depth.
If the vertical mixing coefficient is small (when no mixed layer scheme is used)
the surface stress enters only the top model level, as a body force.
The surface wind stress is calculated in the surface module routines (SBC, see \autoref{chap:SBC}).

The turbulent flux of momentum at the bottom of the ocean is specified through a bottom friction parameterisation
(see \autoref{sec:ZDF_drg})

%% =================================================================================================
\section{External forcings}
\label{sec:DYN_forcing}

Besides the surface and bottom stresses (see the above section)
which are introduced as boundary conditions on the vertical mixing,
three other forcings may enter the dynamical equations by affecting the surface pressure gradient.

(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.

\cmtgm{ missing : the lateral boundary condition !!!   another external forcing
 }

%% =================================================================================================
\section{Wetting and drying }
\label{sec:DYN_wetdry}

There are two main options for wetting and drying code (wd):
(a) an iterative limiter (il) and (b) a directional limiter (dl).
The directional limiter is based on the scheme developed by \cite{warner.defne.ea_CG13} for RO
MS
which was in turn based on ideas developed for POM by \cite{oey_OM06}. The iterative
limiter is a new scheme.  The iterative limiter is activated by setting $\mathrm{ln\_wd\_il} = \mathrm{.true.}$
and $\mathrm{ln\_wd\_dl} = \mathrm{.false.}$. The directional limiter is activated
by setting $\mathrm{ln\_wd\_dl} = \mathrm{.true.}$ and $\mathrm{ln\_wd\_il} = \mathrm{.false.}$.

\begin{listing}
  \nlst{namwad}
  \caption{\forcode{&namwad}}
  \label{lst:namwad}
\end{listing}

The following terminology is used. The depth of the topography (positive downwards)
at each $(i,j)$ point is the quantity stored in array $\mathrm{ht\_wd}$ in the \NEMO\ code.
The height of the free surface (positive upwards) is denoted by $ \mathrm{ssh}$. Given the sign
conventions used, the water depth, $h$, is the height of the free surface plus the depth of the
topography (i.e. $\mathrm{ssh} + \mathrm{ht\_wd}$).

Both wd schemes take all points in the domain below a land elevation of $\mathrm{rn\_wdld}$ to be
covered by water. They require the topography specified with a model
configuration to have negative depths at points where the land is higher than the
topography's reference sea-level. The vertical grid in \NEMO\ is normally computed relative to an
initial state with zero sea surface height elevation.
The user can choose to compute the vertical grid and heights in the model relative to
a non-zero reference height for the free surface. This choice affects the calculation of the metrics and depths
(i.e. the $\mathrm{e3t\_0, ht\_0}$ etc. arrays).

Points where the water depth is less than $\mathrm{rn\_wdmin1}$ are interpreted as ``dry''.
$\mathrm{rn\_wdmin1}$ is usually chosen to be of order $0.05$m but extreme topographies
with very steep slopes require larger values for normal choices of time-step. Surface fluxes
are also switched off for dry cells to prevent freezing, boiling etc. of very thin water layers.
The fluxes are tappered down using a $\mathrm{tanh}$ weighting function
to no flux as the dry limit $\mathrm{rn\_wdmin1}$ is approached. Even wet cells can be very shallow.
The depth at which to start tapering is controlled by the user by setting $\mathrm{rn\_wd\_sbcdep}$.
The fraction $(<1)$ of sufrace fluxes to use at this depth is set by $\mathrm{rn\_wd\_sbcfra}$.

Both versions of the code have been tested in six test cases provided in the WAD\_TEST\_CASES configuration
and in ``realistic'' configurations covering parts of the north-west European shelf.
All these configurations have used pure sigma coordinates. It is expected that
the wetting and drying code will work in domains with more general s-coordinates provided
the coordinates are pure sigma in the region where wetting and drying actually occurs.

The next sub-section descrbies the directional limiter and the following sub-section the iterative limiter.
The final sub-section covers some additional considerations that are relevant to both schemes.

%   Iterative limiters
%% =================================================================================================
\subsection[Directional limiter (\textit{wet\_dry.F90})]{Directional limiter (\mdl{wet\_dry})}
\label{subsec:DYN_wd_directional_limiter}

The principal idea of the directional limiter is that
water should not be allowed to flow out of a dry tracer cell (i.e. one whose water depth is less than \np{rn_wdmin1}{rn\_wdmin1}).

All the changes associated with this option are made to the barotropic solver for the non-linear
free surface code within dynspg\_ts.
On each barotropic sub-step the scheme determines the direction of the flow across each face of all the tracer cells
and sets the flux across the face to zero when the flux is from a dry tracer cell. This prevents cells
whose depth is rn\_wdmin1 or less from drying out further. The scheme does not force $h$ (the water depth) at tracer cells
to be at least the minimum depth and hence is able to conserve mass / volume.

The flux across each $u$-face of a tracer cell is multiplied by a factor zuwdmask (an array which depends on ji and jj).
If the user sets \np[=.false.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} then zuwdmask is 1 when the
flux is from a cell with water depth greater than \np{rn_wdmin1}{rn\_wdmin1} and 0 otherwise. If the user sets
\np[=.true.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} the flux across the face is ramped down as the water depth decreases
from 2 * \np{rn_wdmin1}{rn\_wdmin1} to \np{rn_wdmin1}{rn\_wdmin1}. The use of this ramp reduced grid-scale noise in idealised test cases.

At the point where the flux across a $u$-face is multiplied by zuwdmask , we have chosen
also to multiply the corresponding velocity on the ``now'' step at that face by zuwdmask. We could have
chosen not to do that and to allow fairly large velocities to occur in these ``dry'' cells.
The rationale for setting the velocity to zero is that it is the momentum equations that are being solved
and the total momentum of the upstream cell (treating it as a finite volume) should be considered
to be its depth times its velocity. This depth is considered to be zero at ``dry'' $u$-points consistent with its
treatment in the calculation of the flux of mass across the cell face.

\cite{warner.defne.ea_CG13} state that in their scheme the velocity masks at the cell faces for the baroclinic
timesteps are set to 0 or 1 depending on whether the average of the masks over the barotropic sub-steps is respectively less than
or greater than 0.5. That scheme does not conserve tracers in integrations started from constant tracer
fields (tracers independent of $x$, $y$ and $z$). Our scheme conserves constant tracers because
the velocities used at the tracer cell faces on the baroclinic timesteps are carefully calculated by dynspg\_ts
to equal their mean value during the barotropic steps. If the user sets \np[=.true.]{ln_wd_dl_bc}{ln\_wd\_dl\_bc}, the
baroclinic velocities are also multiplied by a suitably weighted average of zuwdmask.

%   Iterative limiters

%% =================================================================================================
\subsection[Iterative limiter (\textit{wet\_dry.F90})]{Iterative limiter (\mdl{wet\_dry})}
\label{subsec:DYN_wd_iterative_limiter}

%% =================================================================================================
\subsubsection[Iterative flux limiter (\textit{wet\_dry.F90})]{Iterative flux limiter (\mdl{wet\_dry})}
\label{subsec:DYN_wd_il_spg_limiter}

The iterative limiter modifies the fluxes across the faces of cells that are either already ``dry''
or may become dry within the next time-step using an iterative method.

The flux limiter for the barotropic flow (devised by Hedong Liu) can be understood as follows:

The continuity equation for the total water depth in a column
\begin{equation}
  \label{eq:DYN_wd_continuity}
  \frac{\partial h}{\partial t} + \mathbf{\nabla.}(h\mathbf{u}) = 0 .
\end{equation}
can be written in discrete form  as

\begin{align}
  \label{eq:DYN_wd_continuity_2}
  \frac{e_1 e_2}{\Delta t} ( h_{i,j}(t_{n+1}) - h_{i,j}(t_e) )
  &= - ( \mathrm{flxu}_{i+1,j} - \mathrm{flxu}_{i,j}  + \mathrm{flxv}_{i,j+1} - \mathrm{flxv}_{i,j} ) \\
  &= \mathrm{zzflx}_{i,j} .
\end{align}

In the above $h$ is the depth of the water in the column at point $(i,j)$,
$\mathrm{flxu}_{i+1,j}$ is the flux out of the ``eastern'' face of the cell and
$\mathrm{flxv}_{i,j+1}$ the flux out of the ``northern'' face of the cell; $t_{n+1}$ is
the new timestep, $t_e$ is the old timestep (either $t_b$ or $t_n$) and $ \Delta t =
t_{n+1} - t_e$; $e_1 e_2$ is the area of the tracer cells centred at $(i,j)$ and
$\mathrm{zzflx}$ is the sum of the fluxes through all the faces.

The flux limiter splits the flux $\mathrm{zzflx}$ into fluxes that are out of the cell
(zzflxp) and fluxes that are into the cell (zzflxn).  Clearly

\begin{equation}
  \label{eq:DYN_wd_zzflx_p_n_1}
  \mathrm{zzflx}_{i,j} = \mathrm{zzflxp}_{i,j} + \mathrm{zzflxn}_{i,j} .
\end{equation}

The flux limiter iteratively adjusts the fluxes $\mathrm{flxu}$ and $\mathrm{flxv}$ until
none of the cells will ``dry out''. To be precise the fluxes are limited until none of the
cells has water depth less than $\mathrm{rn\_wdmin1}$ on step $n+1$.

Let the fluxes on the $m$th iteration step be denoted by $\mathrm{flxu}^{(m)}$ and
$\mathrm{flxv}^{(m)}$.  Then the adjustment is achieved by seeking a set of coefficients,
$\mathrm{zcoef}_{i,j}^{(m)}$ such that:

\begin{equation}
  \label{eq:DYN_wd_continuity_coef}
  \begin{split}
    \mathrm{zzflxp}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxp}^{(0)}_{i,j} \\
    \mathrm{zzflxn}^{(m)}_{i,j} =& \mathrm{zcoef}_{i,j}^{(m)} \mathrm{zzflxn}^{(0)}_{i,j}
  \end{split}
\end{equation}

where the coefficients are $1.0$ generally but can vary between $0.0$ and $1.0$ around
cells that would otherwise dry.

The iteration is initialised by setting

\begin{equation}
  \label{eq:DYN_wd_zzflx_initial}
  \mathrm{zzflxp^{(0)}}_{i,j} = \mathrm{zzflxp}_{i,j} , \quad  \mathrm{zzflxn^{(0)}}_{i,j} = \mathrm{zzflxn}_{i,j} .
\end{equation}

The fluxes out of cell $(i,j)$ are updated at the $m+1$th iteration if the depth of the
cell on timestep $t_e$, namely $h_{i,j}(t_e)$, is less than the total flux out of the cell
times the timestep divided by the cell area. Using (\autoref{eq:DYN_wd_continuity_2}) this
condition is

\begin{equation}
  \label{eq:DYN_wd_continuity_if}
  h_{i,j}(t_e)  - \mathrm{rn\_wdmin1} <  \frac{\Delta t}{e_1 e_2} ( \mathrm{zzflxp}^{(m)}_{i,j} + \mathrm{zzflxn}^{(m)}_{i,j} ) .
\end{equation}

Rearranging (\autoref{eq:DYN_wd_continuity_if}) we can obtain an expression for the maximum
outward flux that can be allowed and still maintain the minimum wet depth:

\begin{equation}
  \label{eq:DYN_wd_max_flux}
  \begin{split}
    \mathrm{zzflxp}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\
    \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big]
  \end{split}
\end{equation}

Note a small tolerance ($\mathrm{rn\_wdmin2}$) has been introduced here {\itshape [Q: Why is
this necessary/desirable?]}. Substituting from (\autoref{eq:DYN_wd_continuity_coef}) gives an
expression for the coefficient needed to multiply the outward flux at this cell in order
to avoid drying.

\begin{equation}
  \label{eq:DYN_wd_continuity_nxtcoef}
  \begin{split}
    \mathrm{zcoef}^{(m+1)}_{i,j} = \Big[ (h_{i,j}(t_e) & - \mathrm{rn\_wdmin1} - \mathrm{rn\_wdmin2})  \frac{e_1 e_2}{\Delta t} \phantom{]} \\
    \phantom{[} & -  \mathrm{zzflxn}^{(m)}_{i,j} \Big] \frac{1}{ \mathrm{zzflxp}^{(0)}_{i,j} }
  \end{split}
\end{equation}

Only the outward flux components are altered but, of course, outward fluxes from one cell
are inward fluxes to adjacent cells and the balance in these cells may need subsequent
adjustment; hence the iterative nature of this scheme.  Note, for example, that the flux
across the ``eastern'' face of the $(i,j)$th cell is only updated at the $m+1$th iteration
if that flux at the $m$th iteration is out of the $(i,j)$th cell. If that is the case then
the flux across that face is into the $(i+1,j)$ cell and that flux will not be updated by
the calculation for the $(i+1,j)$th cell. In this sense the updates to the fluxes across
the faces of the cells do not ``compete'' (they do not over-write each other) and one
would expect the scheme to converge relatively quickly. The scheme is flux based so
conserves mass. It also conserves constant tracers for the same reason that the
directional limiter does.

%      Surface pressure gradients
%% =================================================================================================
\subsubsection[Modification of surface pressure gradients (\textit{dynhpg.F90})]{Modification of surface pressure gradients (\mdl{dynhpg})}
\label{subsec:DYN_wd_il_spg}

At ``dry'' points the water depth is usually close to $\mathrm{rn\_wdmin1}$. If the
topography is sloping at these points the sea-surface will have a similar slope and there
will hence be very large horizontal pressure gradients at these points. The WAD modifies
the magnitude but not the sign of the surface pressure gradients (zhpi and zhpj) at such
points by mulitplying them by positive factors (zcpx and zcpy respectively) that lie
between $0$ and $1$.

We describe how the scheme works for the ``eastward'' pressure gradient, zhpi, calculated
at the $(i,j)$th $u$-point. The scheme uses the ht\_wd depths and surface heights at the
neighbouring $(i+1,j)$ and $(i,j)$ tracer points.  zcpx is calculated using two logicals
variables, $\mathrm{ll\_tmp1}$ and $\mathrm{ll\_tmp2}$ which are evaluated for each grid
column.  The three possible combinations are illustrated in \autoref{fig:DYN_WAD_dynhpg}.

\begin{figure}[!ht]
  \centering
  \includegraphics[width=0.66\textwidth]{DYN_WAD_dynhpg}
  \caption[Combinations controlling the limiting of the horizontal pressure gradient in
  wetting and drying regimes]{
    Three possible combinations of the logical variables controlling the
    limiting of the horizontal pressure gradient in wetting and drying regimes}
  \label{fig:DYN_WAD_dynhpg}
\end{figure}

The first logical, $\mathrm{ll\_tmp1}$, is set to true if and only if the water depth at
both neighbouring points is greater than $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$ and
the minimum height of the sea surface at the two points is greater than the maximum height
of the topography at the two points:

\begin{equation}
  \label{eq:DYN_ll_tmp1}
  \begin{split}
    \mathrm{ll\_tmp1}  = & \mathrm{MIN(sshn(ji,jj), sshn(ji+1,jj))} > \\
                     & \quad \mathrm{MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj))\  .and.} \\
                     & \mathrm{MAX(sshn(ji,jj) + ht\_wd(ji,jj),} \\
                     & \mathrm{\phantom{MAX(}sshn(ji+1,jj) + ht\_wd(ji+1,jj))} >\\
                     & \quad\quad\mathrm{rn\_wdmin1 + rn\_wdmin2 }
  \end{split}
\end{equation}

The second logical, $\mathrm{ll\_tmp2}$, is set to true if and only if the maximum height
of the sea surface at the two points is greater than the maximum height of the topography
at the two points plus $\mathrm{rn\_wdmin1} + \mathrm{rn\_wdmin2}$

\begin{equation}
  \label{eq:DYN_ll_tmp2}
  \begin{split}
    \mathrm{ ll\_tmp2 } = & \mathrm{( ABS( sshn(ji,jj) - sshn(ji+1,jj) ) > 1.E-12 )\ .AND.}\\
    & \mathrm{( MAX(sshn(ji,jj), sshn(ji+1,jj)) > } \\
    & \mathrm{\phantom{(} MAX(-ht\_wd(ji,jj), -ht\_wd(ji+1,jj)) + rn\_wdmin1 + rn\_wdmin2}) .
  \end{split}
\end{equation}

If $\mathrm{ll\_tmp1}$ is true then the surface pressure gradient, zhpi at the $(i,j)$
point is unmodified. If both logicals are false zhpi is set to zero.

If $\mathrm{ll\_tmp1}$ is true and $\mathrm{ll\_tmp2}$ is false then the surface pressure
gradient is multiplied through by zcpx which is the absolute value of the difference in
the water depths at the two points divided by the difference in the surface heights at the
two points. Thus the sign of the sea surface height gradient is retained but the magnitude
of the pressure force is determined by the difference in water depths rather than the
difference in surface height between the two points. Note that dividing by the difference
between the sea surface heights can be problematic if the heights approach parity. An
additional condition is applied to $\mathrm{ ll\_tmp2 }$ to ensure it is .false. in such
conditions.

%% =================================================================================================
\subsubsection[Additional considerations (\textit{usrdef\_zgr.F90})]{Additional considerations (\mdl{usrdef\_zgr})}
\label{subsec:DYN_WAD_additional}

In the very shallow water where wetting and drying occurs the parametrisation of
bottom drag is clearly very important. In order to promote stability
it is sometimes useful to calculate the bottom drag using an implicit time-stepping approach.

Suitable specifcation of the surface heat flux in wetting and drying domains in forced and
coupled simulations needs further consideration. In order to prevent freezing or boiling
in uncoupled integrations the net surface heat fluxes need to be appropriately limited.

%      The WAD test cases
%% =================================================================================================
\subsection[The WAD test cases (\textit{usrdef\_zgr.F90})]{The WAD test cases (\mdl{usrdef\_zgr})}
\label{subsec:DYN_WAD_test_cases}

See the WAD tests MY\_DOC documention for details of the WAD test cases.

%% =================================================================================================
\section[Time evolution term (\textit{dynnxt.F90})]{Time evolution term (\protect\mdl{dynnxt})}
\label{sec:DYN_nxt}

Options are defined through the \nam{dom}{dom} namelist variables.
The general framework for dynamics time stepping is a leap-frog scheme,
\ie\ a three level centred time scheme associated with an Asselin time filter (cf. \autoref{chap:TD}).
The scheme is applied to the velocity, except when
using the flux form of momentum advection (cf. \autoref{sec:DYN_adv_cor_flux})
in the variable volume case (\np[=.false.]{ln_linssh}{ln\_linssh}),
where it has to be applied to the thickness weighted velocity (see \autoref{sec:SCOORD_momentum})

$\bullet$ vector invariant form or linear free surface
(\np[=.true.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} or \np[=.true.]{ln_linssh}{ln\_linssh}):
\[
  % \label{eq:DYN_nxt_vec}
  \left\{
    \begin{aligned}
      &u^{t+\rdt} = u_f^{t-\rdt} + 2\rdt  \ \text{RHS}_u^t  	\\
      &u_f^t \;\quad = u^t+\gamma \,\left[ {u_f^{t-\rdt} -2u^t+u^{t+\rdt}} \right]
    \end{aligned}
  \right.
\]

$\bullet$ flux form and nonlinear free surface
(\np[=.false.]{ln_dynhpg_vec}{ln\_dynhpg\_vec} and \np[=.false.]{ln_linssh}{ln\_linssh}):
\[
  % \label{eq:DYN_nxt_flux}
  \left\{
    \begin{aligned}
      &\left(e_{3u}\,u\right)^{t+\rdt} = \left(e_{3u}\,u\right)_f^{t-\rdt} + 2\rdt \; e_{3u} \;\text{RHS}_u^t  	\\
      &\left(e_{3u}\,u\right)_f^t \;\quad = \left(e_{3u}\,u\right)^t
      +\gamma \,\left[ {\left(e_{3u}\,u\right)_f^{t-\rdt} -2\left(e_{3u}\,u\right)^t+\left(e_{3u}\,u\right)^{t+\rdt}} \right]
    \end{aligned}
  \right.
\]
where RHS is the right hand side of the momentum equation,
the subscript $f$ denotes filtered values and $\gamma$ is the Asselin coefficient.
$\gamma$ is initialized as \np{nn_atfp}{nn\_atfp} (namelist parameter).
Its default value is \np[=10.e-3]{nn_atfp}{nn\_atfp}.
In both cases, the modified Asselin filter is not applied since perfect conservation is not an issue for
the momentum equations.

\subinc{\input{../../global/epilogue}}

\end{document}