Advection schemes -- grid centers


   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

$\displaystyle \dfrac{\partial {Y}}{\partial {t}} = - \dfrac{\partial {F}}{\partial {z}} \; , \quad$ (234)

where $ F=wY$ is the flux caused by the advective velocity, $ w$.

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

$\displaystyle \dfrac{\partial {Y}}{\partial {t}} = - w\dfrac{\partial {Y}}{\par...
...frac{\partial {F}}{\partial {z}} - Y\dfrac{\partial {w}}{\partial {z}} \right).$ (235)

The discretized form of (234) is

$\displaystyle Y_i^{n+1} = Y_i^n - \dfrac{\Delta t}{h_i} \left( F^n_{i} - F^n_{i-1} \right) \; , \quad$ (236)

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

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

$\displaystyle Y_i^{n+1} = Y_i^n - \dfrac{\Delta t}{h_i} \left( F^n_{i} - F^n_{i-1} -Y_i^n \left(w_k-w_{k-1} \right)\right).$ (237)

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 $ Y_i$ is defined at the grid centers. The fluxes are computed in an upstream-biased way,

$\displaystyle F^n_{i} = \dfrac{1}{\Delta t} \int_{z^\text{Face}_{i} - w \Delta t}^{z^\text{Face}_{i}} Y(z') dz' \quad .$ (238)

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

\begin{displaymath}\begin{array}{rcll} F_{i} &=& w_{i} \left(Y_i + \dfrac{1}{2} ...
...ight) \quad & \text{for} \quad w_{i} < 0 \; , \quad \end{array}\end{displaymath} (239)

where $ c_{i} = 2 w_{i} \Delta t / (h_i+h_{i+1})$ is the Courant number. The factors appearing in (239) are defined as

$\displaystyle \Phi^+_{i} = \alpha_{i} + \beta_{i} r^+_{i} \; , \quad \Phi^-_{i} = \alpha_{i} + \beta_{i} r^-_{i} \; , \quad$ (240)


$\displaystyle \alpha_{i} = \dfrac{1}{2} + \dfrac{1}{6} \left( 1- 2 \vert c_{i} ...
...{i} = \dfrac{1}{2} - \dfrac{1}{6} \left( 1- 2 \vert c_{i} \vert \right) \quad .$ (241)

The upstream and downstream slope parameters are

$\displaystyle r^+_{i} = \dfrac{Y_i - Y_{i-1}}{Y_{i+1}-Y_{i}} \; , \quad r^-_{i} = \dfrac{Y_{i+2} - Y_{i+1}}{Y_{i+1}-Y_{i}} \quad .$ (242)

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

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 $ Y$ 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 $ Y$ 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.


   use util
   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