Advection schemes -- grid centers

*INTERFACE:*

subroutine adv_center(N,dt,h,ho,ww,Bcup,Bcdw,Yup,Ydw,method,mode,Y)

This subroutine solves a one-dimensional advection equation. There are two options, depending whether the advection should be conservative or not. Conservative advection has to be applied when settling of sediment or rising of phytoplankton is considered. In this case the advection is of the form

where is the flux caused by the advective velocity, .

Non-conservative advective transport has to be applied, when the water has a non-zero vertical velocity. In three-dimensional applications, this transport would be conservative, since vertical divergence would be compensated by horizontal convergence and vice versa. However, the key assumption of one-dimensional modelling is horizontal homogeneity, such that we indeed have to apply a vertically non-conservative method, which is of the form

The discretized form of (234) is

where the integers and correspond to the present time and space level, respectively.

For the non-conservative form (235), an extra term needs to be included:

Which advection method is applied is decided by the flag `mode`,
which gives conservative advection (236)
for `mode=1` and
non-conservative advection (237) for `mode=0`.

Fluxes are defined at the grid faces, the variable is defined at the grid centers. The fluxes are computed in an upstream-biased way,

For a third-order polynomial approximation of (see Pietrzak (1998)), these fluxes can be written the in so-called Lax-Wendroff form as

where is the Courant number. The factors appearing in (239) are defined as

where

The upstream and downstream slope parameters are

To obtain monotonic and positive schemes also in the presence of strong
gradients, so-called slope limiters are aplied for the factors
and
. The two most obvious cases are
the first-order upstream discretisation with
and the Lax-Wendroff scheme with
.
The subroutine `adv_center.F90` provides six different slope-limiters,
all discussed in detail by Pietrzak (1998):

- first-order upstream (
`method=UPSTREAM`) - second-order upstream-biased polynomial scheme (
`method=P1`, not yet implemented) - third-order upstream-biased polynomial scheme (
`method=P2`) - third-order scheme (TVD) with Superbee limiter (
`method=Superbee`) - third-order scheme (TVD) with MUSCL limiter (
`method=MUSCL`) - third-order scheme (TVD) with ULTIMATE QUICKEST limiter
(
`method=P2_PDM`)

If during a certain time step the maximum Courant number is larger than one, a split iteration will be carried out which guarantees that the split step Courant numbers are just smaller than 1.

Several kinds of boundary conditions are implemented for the upper
and lower boundaries. They are set by the integer values `Bcup`
and `Bcdw`, that have to correspond to the parameters defined
in the module `util`, see section 8.2. The
following choices exist at the moment:

For the value `flux`, the boundary values `Yup` and `Ydw` are
interpreted as specified fluxes at the uppermost and lowest interface.
Fluxes into the boundary cells are counted positive by convention.
For the value `value`, `Yup` and `Ydw` specify the value
of at the interfaces, and the flux is computed by multiplying with
the (known) speed at the interface. For the value `oneSided`,
`Yup` and `Ydw` are ignored and the flux is computed
from a one-sided first-order upstream discretisation using the speed
at the interface and the value of at the center of the boundary cell.
For the value `zeroDivergence`, the fluxes into and out of the
respective boundary cell are set equal.
This corresponds to a zero-gradient formulation, or to zero
flux divergence in the boundary cells.

Be careful that your boundary conditions are mathematically well defined. For example, specifying an inflow into the boundary cell with the speed at the boundary being directed outward does not make sense.

*USES:*

use util IMPLICIT NONE

number of vertical layers integer, intent(in) :: N time step (s) REALTYPE, intent(in) :: dt layer thickness (m) REALTYPE, intent(in) :: h(0:N) old layer thickness (m) REALTYPE, intent(in) :: ho(0:N) vertical advection speed REALTYPE, intent(in) :: ww(0:N) type of upper BC integer, intent(in) :: Bcup type of lower BC integer, intent(in) :: Bcdw value of upper BC REALTYPE, intent(in) :: Yup value of lower BC REALTYPE, intent(in) :: Ydw type of advection scheme integer, intent(in) :: method advection mode (0: non-conservative, 1: conservative) integer, intent(in) :: mode

REALTYPE :: Y(0:N)

REALTYPE, parameter :: one6th=1.0d0/6.0d0 integer, parameter :: itmax=100

Original author(s): Lars Umlauf

Karsten Bolding 2012-12-28