chap_DYN.tex

      w_i \, \delta_{k} [u] &= \frac{1}{e_{1u}\,e_{2u}} \overline{\overline{ e_{1t}\,e_{2t}\,w_i }}^{\,i,k+1/2} \, \delta_{k} \left[\frac{u}{e_{3uw}}\right]^{up1}, \\[10pt]
      w_i \, \delta_{k} [v] &= \frac{1}{e_{1v}\,e_{2v}} \overline{\overline{ e_{1t}\,e_{2t}\,w_i }}^{\,j,k+1/2} \, \delta_{k} \left[\frac{v}{e_{3vw}}\right]^{up1}.
    \end{aligned}
  \right.
\end{equation}

In flux form case :
\begin{equation}
  \left\{
    \begin{aligned}
      \delta_{k} [w_i u] &= \frac{1}{e_{1u}\,e_{2u}\,e_{3u}} \left( \left( \overline{ e_{1t}\,e_{2t}\,w_i }^{\,i} \, u \right)^{k,up1} - \left( \overline{ e_{1t}\,e_{2t}\,w_i }^{\,i} \, u \right)^{k+1,up1}\right), \\[10pt]
      \delta_{k} [w_i v] &= \frac{1}{e_{1v}\,e_{2v}\,e_{3v}} \left(\left( \overline{ e_{1t}\,e_{2t}\,w_i }^{\,j} \, v \right)^{k,up1} - \left( \overline{ e_{1t}\,e_{2t}\,w_i }^{\,j} \, v \right)^{k+1,up1}\right).
    \end{aligned}
  \right.
\end{equation}

where $w_i$ is the part of the vertical velocity to be treated implicitly
and all other variables are implicit in time (\textit{after}). Note vertical derivatives are 
done with a 1$^\mathrm{st}$ order upstream scheme which provides stability but is diffusive. It is 
therefore important, when using this option, to confirm that the active partitioning of the 
vertical velocity is not frequently required in areas critical for the study being undertaken.

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

There is currently only one choice of limiter for the wetting and drying code
(wd): a directional limiter (dl). Previous versions also provided an iterative
limiter (il) but this has been removed due to performance and robustness issues.
The framework for providing alternatives has been retained in case of future
interest so the directional limiter has to be exlicitly selected despite being
the only choice.

The directional limiter is based on the scheme developed by \cite{warner.defne.ea_CG13}
for ROMS which was in turn based on ideas developed for POM by \cite{oey_OM06}.  The
directional limiter is activated by setting \np[=.true.]{ln_wd_dl}{ln\_wd\_dl}. 

\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 \forcode{ht_wd} in the \NEMO\ code.
The height of the free surface (positive upwards) is denoted by \forcode{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. \forcode{ssh} $+$ \forcode{ht_wd}).

Wetting and Drying schemes take all points in the domain below a land elevation of
\forcode{rn_wdld} to be covered by water. The topography specified with a model
configuration is required 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 "\forcode{ssh_ref}" value may be supplied as the 
\forcode{rn_wd_ref_depth} variable in the domain configuration file. Otherwise it is 
assumed to be zero. This choice affects the calculation of the metrics and depths
(i.e. the \forcode{e3t_0, ht_0} etc. arrays).

Points where the water depth is less than \forcode{rn_wdmin1} are interpreted as ``dry''.
\forcode{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 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 \forcode{rn_wdmin1} is approached. Even wet cells
can be very shallow and may need their surface fluxes reduced.  The depth at which to
start tapering is controlled by the user by setting \forcode{rn_wd_sbcdep}.  The
fraction $(<1)$ of suface fluxes to use at this depth is set by
\forcode{rn_wd_sbcfra}.

The code has been tested in seven 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 describes the directional limiter.  The final sub-section covers some
additional considerations that are relevant to all possible limiting schemes.

\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 \forcode{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 \forcode{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
\forcode{zuwdmask} (an array which depends on ji and jj).  If the user sets
\np[=.false.]{ln_wd_dl_ramp}{ln\_wd\_dl\_ramp} then \forcode{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 \forcode{zuwdmask} , we
have chosen also to multiply the corresponding velocity on the ``now'' step at that face
by \forcode{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 \forcode{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 \forcode{zuwdmask}.

%% =================================================================================================
\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 using
the \forcode{rn_wd_sbcdep} and \forcode{rn_wd_sbcfra} options discussed above.

%      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 - leapfrog (\textit{dynatf.F90})]{Time evolution term - leapfrog (\protect\mdl{dynatf})}
\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 (\key{qco}),
where it has to be applied to the thickness weighted velocity (see \autoref{sec:SCOORD_momentum})
The time integration is done within \mdl{dynzdf}

$\bullet$ vector invariant form or linear free surface
(\np[=.true.]{ln_dynadv_vec}{ln\_dynadv\_vec} or \key{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_dynadv_vec}{ln\_dynadv\_vec} and \key{qco}):
\[
  % \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}