Commit 3fbfd134 authored by hb's avatar hb
Browse files

Source code documentation extended

parent 7e106788
!$Id: bottom_friction.F90,v 1.4 2003-04-23 12:09:43 kbk Exp $
!$Id: bottom_friction.F90,v 1.5 2006-02-04 11:21:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -10,6 +10,28 @@
!
! !DESCRIPTION:
!
! In this routine the bottom friction for the external (vertically integrated)
! mode is calculated. This is done separately for the $U$-equation in the
! U-points and for the $V$-equation in the V-points.
! The drag coefficient $R$ for the external mode is given in eq.\
! (\ref{bottom_vert}) on page \pageref{bottom_vert}.
! For {\tt runtype=1} (only vertically integrated calculations), the
! bottom roughness length is depending on the bed friction
! velocity $u_*^b$ and the molecular viscosity $\nu$:
!
! \begin{equation}\label{Defz0b}
! z_0^b = 0.1 \frac{\nu}{u_*^b} + \left(z^b_0\right)_{\min},
! \end{equation}
!
! see e.g.\ \cite{KAGAN95}, i.e.\ the given roughness may be increased
! by viscous effects.
! After this, the drag coefficient is multiplied by the absolute value of the
! local velocity, which is alculated by dividing the local transports by the
! local water depths and by properly interpolating these velocities
! to the U- and V-points. The resulting fields are {\tt ru}, representing
! $R\sqrt{u^2+v^2}$ on the U-points and {\tt rv}, representing
! this quantity on the V-points.
!
! !USES:
use parameters, only: kappa,avmmol
use domain, only: imin,imax,jmin,jmax,au,av,min_depth
......@@ -27,6 +49,9 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: bottom_friction.F90,v $
! Revision 1.5 2006-02-04 11:21:52 hb
! Source code documentation extended
!
! Revision 1.4 2003-04-23 12:09:43 kbk
! cleaned code + TABS to spaces
!
......
!$Id: cfl_check.F90,v 1.4 2003-04-23 12:09:43 kbk Exp $
!$Id: cfl_check.F90,v 1.5 2006-02-04 11:21:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: cfl_check() - barotropic residual currents.
! !ROUTINE: cfl_check() - check for explicit barotropic time step.
!
! !INTERFACE:
subroutine cfl_check()
!
! !DESCRIPTION:
!
! This routine loops over all horizontal grid points and calculated the
! maximum time step according to the shallow water criterium by
! \cite{BECKERSea93}:
!
! \begin{equation}
! \Delta t_{\max} = \min_{i,j} \left\{\frac{\Delta x_{i,j} \Delta y_{i,j}}
! {\sqrt{2} c_{i,j} \sqrt{\Delta x_{i,j}^2+ \Delta y_{i,j}^2}}\right\}
! \end{equation}
!
! with the local Courant number
!
! \begin{equation}
! c_{i,j}=\sqrt{g H_{i,j}},
! \end{equation}
!
! where $g$ is the gravitational acceleration and $H_{i,j}$ is the local
! bathymetry value. In case that the chosen micro time step $\Delta t_m$
! is larger than $\Delta t_{\max}$, the program will be aborted. In any
! the CFL diagnostics will be written to standard output.
!
!
! !USES:
use parameters, only: g
......@@ -31,6 +53,9 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: cfl_check.F90,v $
! Revision 1.5 2006-02-04 11:21:52 hb
! Source code documentation extended
!
! Revision 1.4 2003-04-23 12:09:43 kbk
! cleaned code + TABS to spaces
!
......
!$Id: depth_update.F90,v 1.6 2005-10-06 09:54:00 hb Exp $
!$Id: depth_update.F90,v 1.7 2006-02-04 11:21:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -10,6 +10,17 @@
!
! !DESCRIPTION:
!
! This routine which is called at every micro time step updates all
! necessary depth related information. These are the water depths in the
! T-, U- and V-points, {\tt D}, {\tt DU} and {\tt DV}, respectively,
! the sea surface elevation in the U- and V-points, {\tt zu} and {\tt zv},
! respectively, and the drying value $\alpha$ defined in equation (\ref{alpha})
! on page \pageref{alpha} in the T-, the U- and the V-points
! ({\tt dry\_z}, {\tt dry\_u} and {\tt dry\_v}).
!
! When working with the option {\tt SLICE\_MODEL}, the water depths in the
! V-points are mirrored from $j=2$ to $j=1$ and $j=3$.
!
! !USES:
use domain, only: imin,imax,jmin,jmax,H,HU,HV,min_depth,crit_depth
use domain, only: az,au,av,dry_z,dry_u,dry_v
......@@ -26,6 +37,9 @@
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: depth_update.F90,v $
! Revision 1.7 2006-02-04 11:21:52 hb
! Source code documentation extended
!
! Revision 1.6 2005-10-06 09:54:00 hb
! added support for vertical slice model - via -DSLICE_MODEL
!
......
!$Id: have_bdy.F90,v 1.5 2004-02-23 15:14:29 kbk Exp $
!$Id: have_bdy.F90,v 1.6 2006-02-04 11:21:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -10,6 +10,10 @@
!
! !DESCRIPTION:
!
! This routine which is called in {\tt domain.F90} checks whether the present
! node has open lateral boundaries. The integer field {\tt bdy\_index}
! is then set accordingly.
!
! !USES:
use domain
use m2d, only: have_boundaries
......@@ -25,6 +29,9 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: have_bdy.F90,v $
! Revision 1.6 2006-02-04 11:21:52 hb
! Source code documentation extended
!
! Revision 1.5 2004-02-23 15:14:29 kbk
! correct mapping of eastern boundary - Staneva
!
......
!$Id: m2d.F90,v 1.12 2006-01-27 16:10:20 kbk Exp $
!$Id: m2d.F90,v 1.13 2006-02-04 11:21:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -9,13 +9,14 @@
module m2d
!
! !DESCRIPTION:
! This modules contains declarations for all variables related to 2D
! This module contains declarations for all variables related to 2D
! hydrodynamical calculations. Information about the calculation domain
! is included from the \emph{domain.F90} module.
! is included from the {\tt domain} module.
! The module contains public subroutines for initialisation, integration
! and clean up of the 2D model component.
! The actual calculation routines are called in integrate\_2d and is linked
! in from the library lib2d.a.
! The actual calculation routines are called in {\tt integrate\_2d}
! and are linked
! in from the library {\tt lib2d.a}.
!
! !USES:
use time, only: julianday,secondsofday
......@@ -41,6 +42,9 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: m2d.F90,v $
! Revision 1.13 2006-02-04 11:21:52 hb
! Source code documentation extended
!
! Revision 1.12 2006-01-27 16:10:20 kbk
! only call do_residual() when specified
!
......@@ -124,7 +128,7 @@
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: init_2d - initialise 2D relatedstuff.
! !IROUTINE: init_2d - initialise 2D related stuff.
!
! !INTERFACE:
subroutine init_2d(runtype,timestep,hotstart)
......@@ -140,7 +144,12 @@
! !OUTPUT PARAMETERS:
!
! !DESCRIPTION:
! Allocates memiory for 2D related fields.
! Here, the {\tt m2d} namelist is read from {\tt getm.inp}, and the check
! for the fulfilment of the CFL criterium for shallow water theory
! {\tt cfl\_check} is called. A major part of this subroutine deals
! then with the setting of local bathymetry values and initial surface
! elevations in $u$- and $v$-points, also by calls to the subroutines
! {\tt uv\_depths} and {\tt depth\_update}.
!
! !REVISION HISTORY:
!
......@@ -192,13 +201,6 @@
LEVEL2 'Format=',bdyfmt_2d
end if
! Boundary related information
if (bdy2d) then
! call have_bdy()
! call print_bdy('Local Boundary Information')
!kbk if (have_boundaries) call init_2d_bdy(bdyfmt_2d,bdyfile_2d)
end if
call uv_depths(vel_depth_method)
where ( -H+min_depth .gt. _ZERO_ )
......@@ -246,7 +248,26 @@
! !OUTPUT PARAMETERS:
!
! !DESCRIPTION:
! A wrapper to call all 2D related subroutines in one subroutine.
! Here, all 2D related subroutines are called. The major calls and their
! meaning are:
!
! \vspace{0.5cm}
!
! \begin{tabular}{ll}
! {\tt call update\_2d\_bdy} & read in new lateral boundary conditions \\
! {\tt call bottom\_friction} & update bottom friction\\
! {\tt call uv\_advect} & calculate 2D advection terms\\
! {\tt call uv\_diffusion} & calculate 2D diffusion terms\\
! {\tt call momentum} & iterate 2D momemtum equations\\
! {\tt call sealevel} & update sea surface elevation\\
! {\tt call depth\_update} & update water depths\\
! {\tt call do\_residual} & calculate intermdediate values for residual currents
! \end{tabular}
!
! \vspace{0.5cm}
!
! It should be noted that some of these calls may be excluded for certain
! compiler options set in the {\tt Makefile} of the application.
!
! !REVISION HISTORY:
! 22Nov Author name Initial code
......@@ -315,7 +336,8 @@
! !OUTPUT PARAMETERS:
!
! !DESCRIPTION:
! This routine cleans up after a 2D integration. Close open files etc.
! This routine executes a final call to {\tt do\_residual} where the residual
! current calculations are finished.
!
! !REVISION HISTORY:
! 22Nov Author name Initial code
......
!$Id: momentum.F90,v 1.9 2006-01-28 20:07:52 hb Exp $
!$Id: momentum.F90,v 1.10 2006-02-04 11:21:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -10,6 +10,11 @@
!
! !DESCRIPTION:
!
! This small routine calls the $U$-equation and the $V$-equation in an
! alternating sequence (UVVUUVVUUVVU), in order to provide higher
! accuracy and energy conservation for the explicit numerical treatment
! of the Coriolis term.
!
! !USES:
use domain, only: imin,imax,jmin,jmax
IMPLICIT NONE
......@@ -28,6 +33,9 @@
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: momentum.F90,v $
! Revision 1.10 2006-02-04 11:21:52 hb
! Source code documentation extended
!
! Revision 1.9 2006-01-28 20:07:52 hb
! Extensions to compiler option SLICE_MODEL for better representation of zero gradients in y-direction
!
......@@ -109,6 +117,50 @@
!
! !DESCRIPTION:
!
! Here, the vertically integrated $U$-momentum equation (\ref{UMOM}) given
! on page \pageref{UMOM} including a
! number of slow terms is calculated. One slight modification is that
! for better stability of drying and flooding processes the slow friction
! term $S^x_F$ is now also multiplied with the parameter $\alpha$ defined
! in eq.\ (\ref{alpha}).
!
! Furthermore, the horizontal pressure gradient $\partial^*_x\zeta$ is
! modified in order to
! support drying and flooding, see figure \ref{figpressgrad} on page
! \pageref{figpressgrad} and the explanations in section \ref{Section_dry}.
! $\partial^*_x\zeta$ is now also considering the atmospheric pressure
! gradient at sea surface height.
!
! For numerical stability reasons, the $U$-momentum equation is here
! discretised in time such that the
! bed friction is treated explicitely:
!
! \begin{equation}\label{Umom_discrete}
! \displaystyle
! U^{n+1}=\frac{U^n-\Delta t_m(gD\partial^*_x\zeta
! +\alpha(-\frac{\tau_x^s}{\rho_0}-fV^n+U_{Ex}+S_A^x-S_D^x+S_B^x+S_F^x))}
! {1+\Delta t_m\frac{R}{D^2}\sqrt{\left(U^n\right)^2+\left(V^n\right)^2}},
! \end{equation}
!
! with $U_{Ex}$ combining advection and diffusion of $U$, see routines
! {\tt uv\_advect} (section \ref{sec-uv-advect} on page
! \pageref{sec-uv-advect}) and {\tt uv\_diffusion}
! (section \ref{sec-uv-diffusion} on page
! \pageref{sec-uv-diffusion}). The slow terms
! are calculated in the routine {\tt slow\_terms} documented in section
! \ref{sec-slow-terms} on page \pageref{sec-slow-terms}.
! In (\ref{Umom_discrete}), $U^{n+1}$ denotes the transport on the
! new and $U^n$ and $V^n$ the transports on the old time level.
!
! The Coriolis term $fU$ for the subsequent $V$-momentum is also calculated
! here, by directly interpolating the $U$-transports to the V-points
! or by a method suggested by \cite{ESPELIDea00} which takes the
! varying water depths into account.
!
! Some provisions for proper behaviour of the $U$-transports when
! GETM runs as slice model are made as well, see section
! \ref{Section_GETM_Slice} on page \pageref{Section_GETM_Slice}.
!
! !USES:
use parameters, only: g,rho_0
use domain, only: imin,imax,jmin,jmax
......@@ -226,6 +278,50 @@
!
! !DESCRIPTION:
!
! Here, the vertically integrated $V$-momentum equation (\ref{VMOM}) given
! on page \pageref{VMOM} including a
! number of slow terms is calculated. One slight modification is that
! for better stability of drying and flooding processes the slow friction
! term $S^y_F$ is now also multiplied with the parameter $\alpha$ defined
! in eq.\ (\ref{alpha}).
!
! Furthermore, the horizontal pressure gradient $\partial^*_y\zeta$ is
! modified in order to
! support drying and flooding, see figure \ref{figpressgrad} on page
! \pageref{figpressgrad} and the explanations in section \ref{Section_dry}.
! $\partial^*_y\zeta$ is now also considering the atmospheric pressure
! gradient at sea surface height.
!
! For numerical stability reasons, the $V$-momentum equation is here
! discretised in time such that the
! bed friction is treated explicitely:
!
! \begin{equation}\label{Vmom_discrete}
! \displaystyle
! V^{n+1}=\frac{V^n-\Delta t_m(gD\partial^*_y\zeta
! +\alpha(-\frac{\tau_y^s}{\rho_0}+fU^n+V_{Ex}+S_A^y-S_D^y+S_B^y+S_F^y))}
! {1+\Delta t_m\frac{R}{D^2}\sqrt{\left(U^n\right)^2+\left(V^n\right)^2}},
! \end{equation}
!
! with $V_{Ex}$ combining advection and diffusion of $V$, see routines
! {\tt uv\_advect} (section \ref{sec-uv-advect} on page
! \pageref{sec-uv-advect}) and {\tt uv\_diffusion}
! (section \ref{sec-uv-diffusion} on page
! \pageref{sec-uv-diffusion}). The slow terms
! are calculated in the routine {\tt slow\_terms} documented in section
! \ref{sec-slow-terms} on page \pageref{sec-slow-terms}.
! In (\ref{Vmom_discrete}), $V^{n+1}$ denotes the transport on the
! new and $U^n$ and $V^n$ the transports on the old time level.
!
! The Coriolis term $fV$ for the subsequent $U$-momentum is also calculated
! here, by directly interpolating the $U$-transports to the U-points
! or by a method suggested by \cite{ESPELIDea00} which takes the
! varying water depths into account.
!
! Some provisions for proper behaviour of the $V$-transports when
! GETM runs as slice model are made as well, see section
! \ref{Section_GETM_Slice} on page \pageref{Section_GETM_Slice}.
!
! !USES:
use parameters, only: g,rho_0
use domain, only: imin,imax,jmin,jmax
......
!$Id: residual.F90,v 1.4 2004-01-03 16:40:28 kbk Exp $
!$Id: residual.F90,v 1.5 2006-02-04 11:21:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -10,6 +10,23 @@
!
! !DESCRIPTION:
!
! Here, the residual transports and depths are integrated up every time step.
! At the end of the simulation, the Eulerian residual currents are
! calculated from:
!
! \begin{equation}
! u_{res} =
! \frac{\displaystyle\int_{t_0}^{t_1}U\,\mbox{d}\tau}
! {\displaystyle\int_{t_0}^{t_1}D^u\,\mbox{d}\tau}, \quad
! v_{res} =
! \frac{\displaystyle\int_{t_0}^{t_1}V\,\mbox{d}\tau}
! {\displaystyle\int_{t_0}^{t_1}D^v\,\mbox{d}\tau},
! \end{equation}
!
! where $t_0$ is the time when the residual calculation begins (to be
! chosen from namelist) and $t_1$ is the finishing time of the model simulation.
!
!
! !USES:
use variables_2d, only: u,v,res_du,res_u,res_dv,res_v,du,dv
IMPLICIT NONE
......@@ -25,6 +42,9 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: residual.F90,v $
! Revision 1.5 2006-02-04 11:21:52 hb
! Source code documentation extended
!
! Revision 1.4 2004-01-03 16:40:28 kbk
! bug fix
!
......
!$Id: sealevel.F90,v 1.8 2006-01-29 12:25:20 kbk Exp $
!$Id: sealevel.F90,v 1.9 2006-02-04 11:21:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -10,6 +10,13 @@
!
! !DESCRIPTION:
!
! Here, the sea surface elevation is iterated according to the vertically
! integrated continuity equation given in (\ref{Elevation}) on page
! \pageref{Elevation}.
!
! When working with the option {\tt SLICE\_MODEL}, the elevations
! at $j=2$ are copied to $j=3$.
!
! !USES:
use domain, only: imin,imax,jmin,jmax,az,H
#if defined(SPHERICAL) || defined(CURVILINEAR)
......@@ -32,6 +39,9 @@
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: sealevel.F90,v $
! Revision 1.9 2006-02-04 11:21:52 hb
! Source code documentation extended
!
! Revision 1.8 2006-01-29 12:25:20 kbk
! NOMADS -> FRESHWATER_LENSE
!
......
!$Id: update_2d_bdy.F90,v 1.5 2006-01-29 12:25:20 kbk Exp $
!$Id: update_2d_bdy.F90,v 1.6 2006-02-04 11:21:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -10,6 +10,15 @@
!
! !DESCRIPTION:
!
! In this routine sea surface elevation boundary conditions are read
! in from a file, interpolated to the actual time step, and distributed
! to the open boundary grid boxes.
! Only for a special test case ({\tt SYLT\_TEST}), ascii data reading is
! supported. For a few special simple cases, analytical calculation
! of boundary elevations is supported. The generic way is reading in
! boundary data from a netcdf file, which is managed in
! {\tt get\_2d\_bdy} via {\tt get\_2d\_bdy\_ncdf}.
!
! !USES:
use domain, only: NWB,NNB,NEB,NSB,H,min_depth,imin,imax,jmin,jmax,az
use domain, only: wi,wfj,wlj,nj,nfi,nli,ei,efj,elj,sj,sfi,sli
......@@ -29,6 +38,9 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: update_2d_bdy.F90,v $
! Revision 1.6 2006-02-04 11:21:52 hb
! Source code documentation extended
!
! Revision 1.5 2006-01-29 12:25:20 kbk
! NOMADS -> FRESHWATER_LENSE
!
......@@ -143,7 +155,6 @@
first = .false.
if (bdyfmt_2d .eq. 1) then
! open(BDYDATA,file='databoun.dat')
open(BDYDATA,file='bdy_data.dat')
t1 = 0.
do i=1,nsbv
......
!$Id: uv_advect.F90,v 1.7 2005-10-06 09:54:00 hb Exp $
!$Id: uv_advect.F90,v 1.8 2006-02-04 11:21:52 hb Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !IROUTINE: uv_advect() - advection of momentum.
! !IROUTINE: uv_advect() - 2D advection of momentum \label{sec-uv-advect}
!
! !INTERFACE:
subroutine uv_advect
!
! !DESCRIPTION:
!
! The advective terms in the vertically integrated
! momentum equation are discretised in
! a momentum-conservative form. This is carried out here for the
! advective terms in the $U$-equation (\ref{UMOM}) and the
! $V$-equation (\ref{VMOM}) (after applying the curvilinear
! coordinate transformationand multiplying these
! equations with $mn$).
!
! First advection term in (\ref{UMOM}):
! \begin{equation}
! \begin{array}{l}
! \displaystyle
! \left(mn\,\partial_{\cal X}\left(\frac{U^2}{Dn}\right)\right)_{i,j}\approx \\ \\
! \quad
! \displaystyle
! \frac{
! \frac12(U_{i+1,j}+U_{i,j})\tilde u_{i+1,j}\Delta y^c_{i+1,j}-
! \frac12(U_{i,j}+U_{i-1,j})\tilde u_{i,j}\Delta y^c_{i,j}
! }{\Delta x^u_{i,j}\Delta y^u_{i,j}}
! \end{array}
! \end{equation}
!
! For the upwind scheme used here, the inter-facial velocities which are defined
! on T-points are here
! calculated as:
!
! \begin{equation}
! \tilde u_{i,j}=
! \left\{
! \begin{array}{ll}
! \displaystyle
! \frac{U_{i-1,j}}{D^u_{i-1,j}} & \mbox{ for } \frac12(U_{i,j}+U_{i-1,j})>0\\ \\
! \displaystyle
! \frac{U_{i,j}}{D^u_{i,j}} & \mbox{ else. }
! \end{array}
! \right.
! \end{equation}
!
! Second advection term in (\ref{UMOM}):
! \begin{equation}
! \begin{array}{l}
! \displaystyle
! \left(mn\,\partial_{\cal Y}y\left(\frac{UV}{Dm}\right)\right)_{i,j,k}\approx \\ \\
! \displaystyle
! \quad
! \frac{
! \frac12(V_{i+1,j}+V_{i,j})\tilde u_{i,j}\Delta x^+_{i,j}-
! \frac12(V_{i+1,j-1}+V_{i,j-1})\tilde u_{i,j-1}\Delta x^+_{i,j-1}
! }{\Delta x^u_{i,j}\Delta y^u_{i,j}}
! \end{array}
! \end{equation}
!
! For the upwind scheme used here, the inter-facial
! velocities which are defined on
! X-points are here
! calculated as:
!
! \begin{equation}
! \tilde u_{i,j}=
! \left\{
! \begin{array}{ll}
! \displaystyle
! \frac{U_{i,j}}{D^u_{i,j}} & \mbox{ for } \frac12(V_{i+1,j,k}+V_{i,j,k})>0\\ \\
! \displaystyle
! \frac{U_{i,j+1}}{D^u_{i,j+1}} & \mbox{ else. }
! \end{array}
! \right.
! \end{equation}
!
! First advection term in (\ref{VMOM}):
! \begin{equation}
! \begin{array}{l}
! \displaystyle
! \left(mn\,\partial_{\cal X}\left(\frac{UV}{Dn}\right)\right)_{i,j,k}\approx \\ \\
! \displaystyle
! \quad
! \frac{
! \frac12(U_{i,j+1}+U_{i,j})\tilde v_{i,j}\Delta y^+_{i,j}-
! \frac12(U_{i-1,j+1}+U_{i-1,j})\tilde v_{i-1,j}\Delta y^+_{i-1,j}
! }{\Delta x^v_{i,j}\Delta y^v_{i,j}}
! \end{array}
! \end{equation}
!
! For the upwind scheme used here, the interfacial
! velocities which are defined on
! X-points are here
! calculated as:
!
! \begin{equation}
! \tilde v_{i,j}=
! \left\{
! \begin{array}{ll}
! \displaystyle
! \frac{V_{i,j}}{D^v_{i,j}} & \mbox{ for } \frac12(U_{i+1,j}+U_{i,j})>0\\ \\
! \displaystyle
! \frac{V_{i+1,j}}{D^v_{i+1,j}} & \mbox{ else. }
! \end{array}
! \right.
! \end{equation}
!
! Second advection term in (\ref{VMOM}):
! \begin{equation}
! \begin{array}{l}
! \displaystyle
! \left(mn\,\partial_{\cal Y}\left(\frac{V^2}{Dm}\right)\right)_{i,j,k}\approx \\ \\
! \quad
! \displaystyle
! \frac{
! \frac12(V_{i,j+1}+V_{i,j})\tilde v_{i,j+1}\Delta x^c_{i,j+1}-
! \frac12(V_{i,j}+V_{i,j-1})\tilde v_{i,j}\Delta x^c_{i,j}
! }{\Delta x^v_{i,j}\Delta y^v_{i,j}}
! \end{array}
! \end{equation}
!
! For the upwind scheme used here, the interfacial velocities which are defined
! on T-points are here
! calculated as:
!
! \begin{equation}
! \tilde v_{i,j}=
! \left\{
! \begin{array}{ll}
! \displaystyle
! \frac{V_{i,j-1}}{D^v_{i,j-1}} & \mbox{ for } \frac12(V_{i,j}+V_{i,j-1})>0\\ \\
! \displaystyle
! \frac{V_{i,j}}{D^v_{i,j}} & \mbox{ else. }
! \end{array}
! \right.
! \end{equation}
!
! When working with the option {\tt SLICE\_MODEL}, the calculation of
! all gradients in $y$-direction is suppressed.
!
! !USES:
use domain,