Model Physics

Initial Conditions
Lateral Boundary Conditions
Surface Boundary Conditions
De-Checkerboarding Operator
Convection Schemes
Implicit Vertical Mixing
Pacanowski & Philander Vertical Mixing
Kraus-Turner Mixed Layer Model
QUICK Tracer Advection
Isopycnal Tracer Mixing

NB. In general many of the optional physics are selected by setting the PHYSICS variable in the Makefile. Although only one option is referred to at a time, there is no reason why a combination of options can not be selected, in which case multiple options should be specified, eg. PHYSICS = -Dopt1 -Dopt2 ...

Initial Conditions

A model run may be started from rest or from a restart, with the init variable (in the ocean.in file) indicating which it is to be (.true. implies that the model is to be started from rest, .false. implies that a restart dataset is to be used).

When the model is started from rest a density distribution is specified, by constructing an idealised potential temperature field, as a function of latitude and depth, based on the zonally averaged annual means from Levitus (1982). Salinity is set to a constant 34.9 ppt, whilst all velocities (both baroclinic and barotropic) and the free surface height field are set to zero.

To change this initialisation requires modification of the init_ocean subroutine, contained within the src/init.F module. For example, it might be preferred to read in the temperature and salinity fields from prepared data files.

To run from a restart, the file must be in a suitable format to be read in by the model and users are referred to the Data Archiving section for more details.

Lateral Boundary Conditions

Land cells are defined by kmt(i,j) = 0. For heat and salt a no-flux condition is assumed across lateral boundaries, whilst a no-slip condition is assumed for momentum.

By default the domain is defined to be finite and closed, such that solid walls are implemented on the edge of all four horizontal boundaries.

kmt(1,j)   = 0                (closed western boundary)
kmt(imt,j) = 0                (closed eastern boundary)
kmt(i,1)   = 0                (closed southern boundary)
kmt(i,jmt) = 0                (closed northern boundary)

An option to remove the solid walls and implement cyclic boundary conditions (in the longitude direction only) is provided by the DOMAIN = -Dcyclic compiler option. This implies that whatever flows out of the east/west boundary will flow into the respective west/east boundary.

kmt(1,j)   = kmt(imt-1,j)     (cyclic western boundary)
kmt(imt,j) = kmt(2,j)         (cyclic eastern boundary)
kmt(i,1)   = 0                (closed southern boundary)
kmt(i,jmt) = 0                (closed northern boundary)

Due to problems with mass conservation when using a free surface formulation, open boundary conditions are not catered for by the SEA model.

Surface Boundary Conditions

In the simplest (default) case the surface stress is based on the zonal averages of Hellerman and Rosenstein (1983) annual mean wind stresses. For surface tracers, the default case is to use zonally averaged Levitus (1982) annual means. The variable gamma (initialised in subroutine init_surface_forcing in the module src/setvbc.F) then governs the rate at which the tracer fields are relaxed toward these values. This relaxation term is converted to a tracer flux using the vertical thickness of the top layer dz(1).

To use alternate sea surface wind stress data or surface tracers, the following Makefile directive is specified.

FORCING = -Dreal_forcing

The model will now expect the surface forcing data to be read in from data files, pre processed into a suitable format and stored in the data directory. This data may either be constant or periodic data. A detailed description of the surface forcing module and how the various data files should be set up is given in the Surface Forcing section.

De-Checkerboarding Operator

Grid splitting of the barotropic mode can, under certain conditions (such as strong horizontal flow or abrupt topography), affect the baroclinic mode, resulting in a +/- checkerboard behaviour (Killworth et al., 1991). To remove this effect, a simple del-plus minus del-cross operator can be applied to the free surface height field.

PHYSICS = -Dde_checkbd

For computational efficiency this operator is applied just once per baroclinic timestep, on the second mini barotropic timestep.

Convection Schemes

Two schemes exist for the explicit treatment of convection. These will only be performed if the implicit vertical mixing option is not also selected.

By default the explicit scheme used is the complete convection scheme, as described by Rahmstorf (1993). This has the effect of mixing all unstable vertical levels in one pass. This is much faster than the old iterative scheme used by Cox (1984), which repetitively mixes pairs of unstable levels. However, the old Cox scheme may be selected instead by setting:

PHYSICS = -Dcox_convect

When using the Cox scheme, the number of iterations is limited by the ncon variable, initialised in src/init.F.

Implicit Vertical Mixing

To solve the vertical diffusion of momentum and tracers using an implicit scheme, the implicit vertical mixing scheme can be selected by setting:

PHYSICS = -Dimplicitvmix

The degree of implicitness is then governed by the aidif variable in ocean.in. When "aidif=1.0" treatment will be fully implicit, whereas "aidif=0.0" will result in fully explicit treatment. Setting "aidif=0.5" results in half the vertical diffusion being solved implicitly and half explicitly (this is the Crank-Nicholson scheme and is a good recommendation).

Note, that when implicit vertical mixing is enabled, no explicit convection is performed, as this is achieved by setting the diffusivity to a large value in regions of static instability.

Pacanowski & Philander Vertical Mixing

The basic vertical mixing parameterisation is to use constant values for the coefficients of viscosity and diffusivity (fkpm and fkph in ocean.in). These are independent of both space and time, such that for all k=1,km,

vvc(k) = fkpm
vdc(k) = fkph

Ocean observations indicate intense mixing processes in surface layers above the thermocline, but weak mixing below, except for the region below the core of the Equatorial Undercurrent. Pacanowski and Philander (1981) show that models with constant mixing coefficients behave poorly in this and go on to describe a parameterisation of this mixing process, by means of Richardson-number dependent coefficients of eddy mixing, which are therefore functions of space and time.

For a more realistic representation of heat and momentum in equatorial regions, the Pacanowski and Philander vertical mixing option should be selected.

PHYSICS = -Dppvmix

The coefficients of viscosity and diffusivity, that are Richardson-number dependent, are then calculated as follows.

if (riu >= 0)
  vvc = MAX ( [(fricmx / (1 + 5 * riu)**2) + bvvc], [wndmix] )
else
  vvc = vvclim
  
if (rit >= 0)  
  vdc = MAX ( [(fricmx / (1 + 5 * rit)**3) + bvdc], [wndmix] )
else
  vdc = vdclim

For each timestep, Richardson numbers are first calulated at bottom of u grid points to give riu, then averaged to give rit at bottom of t grid points. Six parameters, used by this option, are declared in include/cvmix.h and initialised in the subroutine init_sea in the src/init.F module. They have the following interpretation:

fricmx = max vertical mixing coefficient (cm**2/sec)
         (typically in range: 50.0 .. 100.0)
wndmix = min vertical mixing in level 1 to simulate wind mixing
         (arbitrary value meant to simulate wind bursts)
bvvc   = background vertical viscosity constant (cm**2/sec)
         (typically in range: 0.0134 .. 1.0)
bvdc   = background vertical diffusion constant (cm**2/sec)
         (typically in range: 0.00134 .. 0.1)
vvclim = vertical viscosity coefficient limit
vdclim = vertical diffusion coefficient limit

Kraus-Turner Mixed Layer Model

In the highly turbulent surface region of the ocean, there is a considerable generation of turbulent kinetic energy (tke), which may lead to vertical mixing in the water column. The Kraus-Turner (1967) module produces a mixed layer that includes the mechanical effects of wind stirring and of entrainment across the lower boundary of the mixed layer. This results in a model that responds more realistically to changing winds, which has been shown to be particularly important when studying seasonal variability.

PHYSICS = -Dkraust

Suggested values for the Kraus-Turner parameters (and an overview of the model equations) are given in the U.K. Met. Office Technical Report No.135 (Notes on how to choose parameter values for the Cox numerical ocean circulation model, ed. M.J.Bell, November 1994):

kt_lamda  = efficiency of the wind mixing
            (typically: 0.5 .. 1.0 / recommend: 0.7)
kt_epsln  = efficiency of buoyancy driven mixing
            (typically: 0.14 .. 0.17 / recommend: 0.15)
kt_wndh   = wind mixing scale depth
            (typically: 2000 .. 20,000cm)

QUICK tracer advection

In many regions of the ocean, advection plays a dominant role in the evolution of tracer fields, but in regions were large gradients and/or velocities are present, the central difference advection scheme used in many Bryan-Semtner-Cox type models (such as the SEA model) can lead to significant under/over shooting of tracer values. The problem tends to become more apparent when finer spatial grid sizes are used (since for coarse resolution models large coefficients of eddy diffusivity are used, resulting in highly diffusive solutions with limited advection effects). To reduce these effects in fine resolution models Farrow and Stevens (1995) adapted the QUICK tracer advection scheme (of Leonard, 1979) to use within the tracer component of Bryan-Semtner-Cox type models. It is this scheme that was originally implemented into the SEA model, but this has since been superceeded by a modified version, as described by Webb et al. (1998).

PHYSICS = -Dquick

The third order predict-correct QUICK scheme is split into a fourth order central difference representation of advection and a velocity dependent biharmonic mixing term (lagged by one timestep), leading to a stable one-step scheme.

Isopycnal Tracer Mixing

It is believed that mixing due to mesoscale motions takes place primarily along isopycnal surfaces, but these motions are not resolved in the standard model. Particular critisism is directed towards the lateral diffusion which is aligned with the model grid. Redi (1982) introduces a more complex mixing tensor, which senses the local isopycnal slope and conducts diffusion in that direction. The implementation of this in SEA is based on that of Cox (1987). For computational efficiency, densities are referenced to a single depth level (see isop_ref in src/init.F to specify which level), rather than using as many reference levels as there are vertical levels. When using the isopycnal mixing scheme, implicit vertical mixing must also be enabled, using the following directives.

PHYSICS = -Dimplicitvmix -Disopycmix

To calculate the effect of mesoscale eddies on isopycnals the Gent & McWilliams scheme (1990) will do so in terms of eddy-induced transport (advective) velocities. Select the following directives to include this parameterisation.

PHYSICS = -Dimplicitvmix -Disopycmix -Dgent_mcw

Various parameters, declared in include/isop_gm90.h may be specified for these modules, which are set in the subroutine init_sea in the src/init.F:

isop_ah    = isopycnal tracer diffusivity (cm**2/sec)
isop_slmxr = maximum isopycnal slope allowed in mixing scheme
             (a value of 100.0 translates to a slope of 1:100)
isop_ref   = level at which to reference densities
gm90_kappa = isopycnal thickness diffusivity (cm**2/sec)