next up previous
Next: saltflux.F Up: The Proudman Oceanographic Laboratory Previous: metset.F

Subsections


pgrad.F


subroutine pgrad

Description

This subroutine calculates the buoyancy contribution to the horizontal pressure gradient term, $-\nabla P$ (where $\nabla$ is the horizontal gradient operator). To avoid errors with steep topography, we calculate this term by interpolating the buoyancy field to z-coordinates then integrating to calculate the pressure. Fig. 49.1 shows the arrangement at a grid point. This is then differenced to give the pressure gradient at the central u-point.

The pressure is written

\begin{displaymath}
P= P_a+\rho_0(\psi + g \zeta - gz)+F(Z)
\end{displaymath} (42)

Where $F(Z)=0.002282gZ^2$ is an arbitrary function of $Z$ included in the equation of state calcucalculation for accuracy (value here is determined from the initial condition of the CS3 model - model results should not be sensitive to this). Here we are estimating $- \nabla \psi$ where $\psi= H \int_0 ^{\sigma} b d \sigma $. The depth of the plane on which the pressure is to be calculated (Figure 49.1) is $z= \sigma HU_q + \zeta U_q $. The four corners are labeled $q=1..4$, so $b$-column, $q$, this is at $\sigma_q =(z-\zeta_q)/H_q $. The nearest $\sigma $-levels, $k_q$, above the plane (shown as double circles in Fig. 49.1) lie at a fractional distance, $r=(\sigma_{k_q}
- \sigma_q)/(\sigma_{k_q} - \sigma_{k_q-1})$, from the plane. The pressure at $\sigma_q$ can be estimated to be
$\displaystyle \psi_q=$ $\textstyle -0.5 \delta \sigma
[\sum\limits_{k=k_q}^{N-1} (b_{q,k} +b_{q,k+1})$    
  $\textstyle +(1-2r+r^2)b_{q,k_q-1}+(1-r^2)b_{q,k_q } ],$   (43)

with no summation for $k_q=N-1$. In this calculation level $N$ is at the sea surface, where the buoyancy is defined by $b_{q,N-1}=b_{0 q,N-2}$.

The pressure gradients along the edges of the plane are

\begin{displaymath}
\Delta \psi_1 = \psi_1-\psi_2, \quad, \Delta \psi_2 =
\psi_2...
...a \psi_3 = \psi_4-\psi_3, \quad \Delta \psi_4
= \psi_1-\psi_4
\end{displaymath} (44)

and the corresponding changes in velocity are
$\displaystyle \Delta u_{i,j,k} =-0.5 \frac{\Delta
t}{Rcos \phi \Delta \chi} (\Delta \psi_1 +\Delta \psi_3)$      
$\displaystyle \Delta v_{i,j,k} =-0.5 \frac{\Delta t}{R \Delta \phi} (\Delta \psi_2
+\Delta \psi_4)$     (45)

Figure 4: Arrangement for the calculation of pressure gradients on a horizontal plane
\includegraphics[scale=1]{pgradgridpoint.ps}

This technique allows a straightforward treatment of cases where the sea bed lies above the plane at one or more of the surrounding $b$-points; if either of the corners on an edge are below the sea bed then the corresponding pressure gradient takes the last defined value above it. Alternatively these edges can be removed from the calculation.

Subroutine Arguments

none

Local variables

i,j,k, ip, jp, kp, ks local indices
icg,jcg global indices
i4 loop count
icliff, it which points are below the sea bed /above the surface
inextp,jnextp(4), ipnt, jpnt, kpnt direction of b-point from u-point
tagh,tagzet,tagb,tagdi tags for exchanges
press pressure contribution between $\sigma $-level and z-level.
di pressure contribution at a $\sigma $-level
dp pressure gradient along an edge
dp0last non-zero pressure gradient along an edge
r relative vertical coordinate within a grid cell
sigp zz vertical position of the horizontal plane
first .true. the first time the subroutine is called

Global variables changed

u v
kpgrad
b, pressure

Logical units

none

Order of Things

  1. add pressure dependence to equation of state and iterate pressure (if used)
  2. set b at sea surface and sea bed
  3. integrate density on sigma levels
  4. calculate extra pressure between $\sigma $-level and z-level.
  5. calculate pressure difference along edges
  6. remember last non-zero difference
  7. update velocity

Calls

bcalcP
bcalc
Exch3DS
Exch3DR

Called By

baroc

Options - Logical

compress use compressibility in equation of state

Options - Compiler

SCOORD use horizontal varying vertical coordinates
NOTEVENSIG $\sigma $-level are not evenly spaced

Known Issues

Bottom boundary condition does not reflect the correct no normal gradient condition. This is a finite difference calculation in a different coordinate system from the rest of the model, so there might be conservation issues (p.e. /p.v. ?), but none have been identified as yet.


subroutine pgrad_sigma

Description

No documentation yet. If you require information about this subroutine please contact POL.

Subroutine Arguments

none

Local variables

Global variables changed

Logical units

Order of Things

Calls

Called By

Options - Logical

none

Options - Compiler

none

Known Issues


subroutine pgrad_spline

Description

This subroutine calculates the buoyancy contribution to the horizontal pressure gradient term, $-\nabla P$ (where $\nabla$ is the horizontal gradient operator). To avoid errors with steep topography, we calculate this term by interpolating the buoyancy field to z-coordinates then integrating to calculate the pressure. Here cubic spline interpolation is used.

In this case the pressure at $\sigma_q$ can be estimated to be

$\displaystyle \psi_q=$ $\textstyle -0.5 \delta \sigma
[\sum\limits_{k=k_q}^{N-1} (b_{q,k} +b_{q,k+1}-\delta \sigma ^2 /12
(b''_{q,k} +b''_{q,k+1}))$    
  $\textstyle +(1-2r+r^2)b_{q,k_q-1}+(1-r^2)b_{q,k_q }$    
  $\textstyle + \delta \sigma ^2
/6[(0.5(1-r)^4-1+r^2)b''_{q,k-1}+(0.5-0.5r^4-1+r^2)b''_{q,k}] ],$   (46)

with no summation for $k_q=N-1$. In this calculation level $N$ is at the sea surface, where the buoyancy is defined by $b_{q,N-1}=b_{0 q,N-2}$.

Apart from this, the subroutine is identical to pgrad, where a full description can be found.

Subroutine Arguments

none

Local variables

i,j,k, ip, jp, kp, ks local indices
icg,jcg global indices
i4 loop count
icliff, it which points are below the sea bed /above the surface
inextp,jnextp(4), ipnt, jpnt, kpnt direction of b-point from u-point
tagh,tagzet,tagb,tagdi,tagbpp tags for exchanges
press pressure contribution between $\sigma $-level and z-level.
di pressure contribution at a $\sigma $-level
Bpp second deivative of b from spline calculation
dp pressure gradient along an edge
dp0last non-zero pressure gradient along an edge
r relative vertical coordinate within a grid cell
sigp zz vertical position of the horizontal plane
first .true. the first time the subroutine is called

Global variables changed

u v
kpgrad
b, pressure

Logical units

none

Order of Things

  1. add pressure dependence to equation of state and iterate pressure (if used)
  2. set b at sea surface and sea bed
  3. fitspline to b
  4. integrate density on sigma levels
  5. calculate extra pressure between $\sigma $-level and z-level.
  6. calculate pressure difference along edges
  7. remember last non-zero difference
  8. update velocity

Calls

bcalcP
bcalc
Exch3DS
Exch3DR

Called By

baroc

Options - Logical

compress use compressibility in equation of state

Options - Compiler

SCOORD use horizontal varying vertical coordinates
NOTEVENSIG $\sigma $-level are not evenly spaced

Known Issues

Bottom boundary condition does not reflect the correct no normal gradient condition. This is a finite difference calculation in a different coordinate system from the rest of the model, so there might be conservation issues (p.e. /p.v. ?), but none have been identified as yet.


subroutine spline_fit(c,cpp)

Description

No documentation yet. If you require information about this subroutine please contact POL.

Subroutine Arguments

none

Local variables

Global variables changed

Logical units

Order of Things

Calls

Called By

Options - Logical

none

Options - Compiler

none

Known Issues


next up previous
Next: saltflux.F Up: The Proudman Oceanographic Laboratory Previous: metset.F
The AMMP Project 2005-04-20