Subroutine VVELOC

Given an initial distribution of layer thicknesses (H) and some function (F) calculate the amount of mass that must be moved across the interfaces (W) in order to bring the interfacial values of F back toward the desired values ( F_desired )

COORDINATE DEFINITION FUNCTION:

(1) The coordinate definition function is presumed to be a monotonic decreasing function of the state variables. The value passed as F(:,:,:) is the average value of this function over each layer.

(2) The F_Desired value is the value of the coordinate definition function that we would like to have as an interpolated value at the bottom of the interface.

(3) The damping to F_Desired is done over a time F_DecayTime(k). This means that F_DecayTime is a function of the layer, so that if we want layers near the top to be more strongly damped than layers near the bottom, we can do so.

The algorithm does not care what function you use for F so long as F and F_Desired are consistent and F is a monotonic decreasing function over k.

The exact nature of the function would be determined by the calling program.

Examples:

     F = buoyancy
     F = depth
     F = EXP( -k*Pressure ) + Specific Volume anomaly

MIXED LAYER AND BUFFER REGION CONSIDERATION:

if MLK > 1 then the layers 1..mlk represent the "mixed" layer

MATERL is the first layer which should be considered for isopycnal treatment

if MATERL > MLK+1, then the region between mlk and k=materl is subdivided into equally spaced layers.

The mixed layer gets precedence over all deeper layers. Pmax(:,:,1) will prevent the mixed layer from becoming so deep that there is no room for other layers

ISO-SURFACE REGION:

In the deep layers, the model attempts to keep the layers along iso- surfaces. Depending on the choice of the coordinate definition function, this may be isopycnal iso-neutral, isothermal, etc.

The algorithm attempts to move mass so that a linearly interpolated value of the coordinate definition function at the layer interfaces is kept close to the desired interfacial value through simple relaxation.

w = %DELTA F {partial z}over{partial F} / T

where then gradient is computed as

{partial F}over{partial z} = %delta_k F / overline{H}

and the difference is

%DELTA F = { F_k H_{k+1} + F_{k+1} H_k }over{h_k + h_{k+1}} - F_d

The timescale T is layer dependent. The argument F_DecayTime gives the decay time for each layer. This permits layers near the surface to be treated differently from layers further down.

LIMITS:

Although the layers try to behave acording to our algorithms, it is possible that no solution exists for the desired coordinate definiton function, or that layers would lie beneath topography or out of the surface.

The model imposes a few limits on the layer behavior:

First, impose max,min limits on W.

Second, go down and enforce new hmin < H < hmax. This may push layers down from the mixed layer, and cause them to bump into the bottom.

Third, go up and enforce P < Pmax(:,:,k) where Pmax(:,:,k) is the bottom pressure minus the sum of Hmin_b betwen our layer and the bottom. This enforces the condition that at least hmin_b(k) will remain in the layer

Syntax:

call VVELOC(g,extmode,H,F,F_Desired,W,materl,mlk,Entrainment,M_O_Depth,TimeInterval,&

F_DecayTime)

INTENT(IN) :: g,extmode,H,F,F_Desired,materl,mlk,Entrainment,M_O_Depth,TimeInterval,F_DecayTime
INTENT(OUT) :: W

Subroutines Called: