Next: advpbu_n.F Up: The Proudman Oceanographic Laboratory Previous: advect_vel.F

Subsections

subroutine advpbu (ba, hin, hout )
subroutine advpbu (ba, hin, hout )

advpbu.F

subroutine advpbu (ba, hin, hout )

Description

The subroutine performs advection of any scalar variable in the

(or

) direction by PPM (the piecewise parabolic method). This method ensures positivity and scalar conservation with a minimum of numerical diffusion. The routine is directionally split, first advecting in the

direction, then in the vertical. The routine returns new values of the scalar at the original $\sigma$ or

levels and a new total depth. The vertical advection is performed by considering control volumes, which are remapped on to the original $\sigma$ or

levels. An alternative (`new') method is chosen if the CFL condition is violated in the vertical; this method allows any value of vertical Courant number.

When followed by advpbv, which advects in the (or ) direction, followed by a further vertical advection, a complete directionally split advection cycle has been performed, and a final total depth has been returned. Since advection is performed on the longer baroclinic time step DT, the final depth, which is given by the volume fluxes, should be equal to the final depth over the corresponding number of barotropic time steps. This is ensured by making the volume fluxes used equal to the sum of barotropic volume fluxes. The routine advect_sca controls the order of advection in the two horizontal directions: this may alternate on successive time steps.

A laplacian horizontal diffusive flux can be included at the intermediate stage, with either a fixed or Smagorinsky diffusity (see htmlrefhorizdiffusivityhorizdiffusivity)

Numerical Equations

PPM involves fitting a parabola to the variable within each grid interval. The parabola may reduce to a straight line or constant near boundaries where there are not enough points to fit a parabola. This is appropriate for one dimension. Directional splitting allows the sequential use of one-dimensional routines. Fluxes of the scalar variable

across a coordinate face are derived by integrating the parabola over a distance calculated from the given volume flux across the face.

Consider an interval to $x_{i+1}$ , with $x_{i+1} - x_i = \Delta x$ . Put $\xi=(x-x_i) / \Delta x$ . Then the parabola is

$\begin{displaymath} b(\xi) = b_L + \xi(\Delta b + b_6 (1-\xi)) \end{displaymath}$

where $\Delta b = b_R - b_L$ and $b_6 = 6\bar{b} - 3(b_L +b_R)$ , where $\bar{b}$ is the mean value of

over the interval.

The edge values $b_{i+1/2} = b_{R,i} = b_{L,i+1}$ of the parabolas are given for equal intervals by $b_{i+1/2} = (b_i + b_{i+1})/2 + (\delta b_i - \delta b_{i+1})/6$ where $\delta b_i=(b_{i+1} - b_{i-1})/2$ . For the formula for unequal intervals see Collela P. and Woodward P.R. 1984 Journal of Computational Physics 54:174-201. This formula is used in the vertical advection.

$\delta b_i$ is replaced by

$\begin{displaymath} \delta_m b_i = \mathrm{min} (\vert\delta b_i\vert,2\vert b_i... ...i-1}\vert,2\vert b_{i+1} - b_i\vert) \mathrm{sgn} (\delta b_i) \end{displaymath}$

if $(b_i - b_{i-1}).(b_{i+1} - b_i) > 0$ , zero otherwise.

Steepening (optional) tests for the presence of a discontinuity. This is deemed to occur if the third derivative is large, the second derivative changes sign and the first and third derivatives have opposite signs.

If $\delta^2 b_i = (b_{i+1}-2b_i+b_{i-1})/6$ ,
let $\tilde{\eta_i} = -(\delta^2 b_{i+1} - \delta^2 b_{i-1})/(b_{i+1} - b_i)$ if $\delta^2 b_{i+1}.\delta^2 b_{i-1} < 0$ and $\vert\delta b_i\vert> \epsilon \mathrm{min} (\vert b_{i+1}\vert, \vert b_{i-1}\vert)$ , otherwise $\tilde{\eta_i} = 0$ .
Put $\eta_i = \mathrm{max} (0, \mathrm{min} [\eta_1(\tilde{\eta_i}-\eta_2),1])$ , where $\eta_1$ , $\eta_2$ and $\epsilon$ are variable parameters. Then, set the discontinuity values $b^d_{L,i}= b_{i-1} + \delta_m b_{i-1}/2$ and $b^d_{R,i}= b_{i+1} - \delta_m b_{i+1}/2$ , and for a smooth change replace $b_{L,i}$ by $b_{L,i}(1-\eta_i) + b^d_{L,i}\eta_i$ and $b_{R,i}$ by $b_{R,i}(1-\eta_i) + b^d_{R,i}\eta_i$ .

Monotonisation is enforced by the following conditions:
put $b_{L,i}=b_{R,i}=b_i$ if $(b_{R,i}-b_i)(b_i-b_{L,i}) \le 0$
put $b_{R,i}=3b_i - 2b_{L,i}$ if $\Delta b^2 < -\Delta b.b_6$ .

Subroutine Arguments

: ba value of scalar variable, returns new value
: hin initial depth
: hout final depth

Local variables

: b2 intermediate value of advected scalar
: bl,br value of parabola at left and right of interval
: d2b value of $\delta^2 b$
: dba difference of
: dfl2 flux difference
: dmb value of $\delta_m b$
: ds2 difference between intermediate sigma levels
: dsig difference between intermediate and original sigma levels
: fba scalar flux
: fuu volume flux
: sigo3 intermediate sigma levels
: i,j,k,icg,jcg,ilb,iub grid indices
: tagba,tagdba,tagdmb,tagd2b,tagu,tagh,tagbr,tagbl ,tagfu,tagfba integers used in Exch3DS and Exch3DR routines
: a1,a2,a3,b6,bdl,bdr,delb,delb2,delb6 variables used in parabola calculations (see above description)
: eti,ett variables used in steepening option
: e1,e2,eps parameters used in steepening option
: aip mask variable used in flux calculations
: rcdal length of grid box side in lat-long coordinates
: xdb used in setting parabolas in vertical with uneven grid spacing
: xi length over which parabola is integrated to give flux
: advmask logical variable based on mask
: ivi 0 or 1 whether or not vertical Courant number is
: k1,k0,ka,kb,kc vertical indices used in `new' vertical advection scheme
: dsi sigma coordinate differences used in `new' vertical advection scheme
: baa,bab,bat variables used in `new' vertical advection scheme
: inv3,inv6 1/3,1/6
: HD imposed horizontal diffusive flux
: pri turbulent Prandlt number (inverse)