Advection schemes -- grid centers (Source File: adv_center.F90)

INTERFACE:

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

DESCRIPTION:

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

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

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

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

The discretized form of (232) is

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

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

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

$$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).$$ (235)

Which advection method is applied is decided by the flag mode, which gives conservative advection (234) for mode=1 and non-conservative advection (235) 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,

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

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{array}{rcll} F_{i} &=& w_{i} \left(Y_i + \dfrac{1}{2} ... ...ight) \quad & \text{for} \quad w_{i} < 0 \; , \quad \end{array}$$ (237)

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

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

where

$$\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 .$$ (239)

The upstream and downstream slope parameters are

$$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 .$$ (240)

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):

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

USES:

    use util
    IMPLICIT NONE

INPUT PARAMETERS:

 
    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

INPUT/OUTPUT PARAMETERS:

    REALTYPE                            :: Y(0:N)

DEFINED PARAMETERS:

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

REVISION HISTORY:

    Original author(s): Lars Umlauf
    $Log: adv_center.F90,v $
    Revision 1.5  2010-09-17 12:53:52  jorn
    extensive code clean-up to ensure proper initialization and clean-up of all variables
    Revision 1.4  2006-11-06 13:36:46  hb
    Option for conservative vertical advection added to adv_center
    Revision 1.3  2006-03-20 09:06:38  kbk
    removed explicit double precission dependency
    Revision 1.2  2005/11/18 10:59:34  kbk
    removed unused variables - some left in parameter lists
    Revision 1.1  2005/06/27 10:54:33  kbk
    new files needed