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

Subsections


tubulence.F


subroutine turbulence_closure

Description

Calculate turbulent kinetic energy, mixing length and hence vertical mixing coeficients eddy viscosity and diffusivity. Uses a Mellor-Yamada level 2.5 closure scheme. It has two compiling options one to calculate the mixing length l (al) algebraically, but limited by stratification using the Bakhmetev scale; or to solve the $q^2l$ differential equation. It first calculate turbulent kinetic energy, eddy viscosity, etc. at u points and the at b points.

The subroutine uses the Galperin et al (1988)stability relations. The surface boundary conditions includes the wind waves effect as in Craig and Banner (1994), while the bottom boundary condition assumes a steady state balance between shear production and dissipation.

Numerical Equations


\begin{displaymath}
\frac{\partial q^2}{\partial t} = 2K_ZM^2 - 2\left( {K_H - \alpha K_Z}
\right)N^2 - \frac{2q^3}{B_1 l} + D(q^2)
\end{displaymath}

where

\begin{displaymath}
\begin{array}{l}
M^2 = H^{ - 2}\left[ {\left( {\frac{\parti...
...H^{ - 1}\frac{\partial b_0 }{\partial \sigma } \\
\end{array}\end{displaymath}

and $D(q^2)$ is the TKE diffusive term. With B$_{1}$ = 16.6. The term $\alpha K_Z$ (with $\alpha$ = 0.7) is a simple representation of vertical mixing by long wavelength internal waves following Mellor (1989).

The algebraic mixing l length used here is:

\begin{displaymath}
l = l_\sigma H
\end{displaymath}

where l$_{\sigma }$ is:

\begin{displaymath}
l_\sigma = \kappa \left( {1 - \sigma } \right)\left( { - \si...
...right)^{\frac{1}{{\begin{array}{l}
2 \\
\\
\end{array}}}}
\end{displaymath}

or in case of deep water:


\begin{displaymath}
\begin{array}{l}
l_\sigma = \kappa \left( {\frac{h_c }{h}} ...
...d \quad \quad
\sigma \le \frac{2h_c }{3h} - 1 \\
\end{array}\end{displaymath}

where h$_{c}$ is a critical depth set to 150 m.

\begin{displaymath}
\begin{array}{l}
l = l\left( {\left( {Ri^2 + 1} \right)^{\f...
...al z}}}
\right)\quad \quad \quad \quad other. \\
\end{array}\end{displaymath}

The bed boundary condition following steady state balance between shear production and dissipation:

\begin{displaymath}
q^2 = B_1 ^{\frac{2}{3}}\left( {F_b ^2 + G_b ^2} \right)^{\frac{1}{2}}
\end{displaymath}

While the surface boundary condition is the uses the Craig and Banner formulation, which includes the effect of surface breaking waves and assumes a balance between diffusion, shear production and dissipation.

\begin{displaymath}
q^2 = \alpha _{wave}^{\frac{2}{3}} \left( {3\frac{B_1 }{S_q ...
...t)^{\frac{1}{3}}\left( {F_s ^2 + G_s ^2} \right)^{\frac{1}{2}}
\end{displaymath}

with $\alpha_{wave}$ = 100.

The system is then closed by

\begin{displaymath}
K_Z = S_M lq\quad \quad K_H = S_H lq\quad \quad K_q = S_q lq
\end{displaymath}

which uses the stability relations from Galperin et al.(1988)

\begin{displaymath}
\begin{array}{l}
S_M = \frac{0.3933 - 3.086G}{\left( {1 - 3...
... S_H = \frac{0.4939}{1 - 34.68G}, \\
S_q = 2 \\
\end{array}\end{displaymath}

With G = -l$^{2}$N$^{2}$ / q$^{2}$.

Subroutine Arguments

None

Local variables

i,j,k,iig,jjg Grid indeces
tagb,tagh,tagqsq,tagaln Communication tags for MPI
first Logical

Global variables changed

qsq,aln twice the turbulent kinetic energy and micing length

Logical units

none

Order of Things

  1. Exchange values of b and h.
  2. First time called allocate variables and set initial value of al.
  3. Looping through u-points:
  4. Calculates bulk gradients.
  5. Calculates the mixing length if using the algebraic option (default).
  6. Calculates Stability funtions Sm and Sh
  7. Calculates the Changes of TKE due to shear strees and buoyancy.
  8. Calculates the Changes of TKE due to diffusion.
  9. Updates TKE
  10. Updates Mixing length scale either algebraically or using the second differential equation.
  11. Calculates stability functions.
  12. Calculates $K_Z$
  13. Exchanges qsq and al
  14. Looping through b-points:
  15. Calculates qsq and al
  16. Calculates stability functions
  17. Calculates $K_H$

Calls

Exch3DS,Exch3DR
al_init
bulk_gradients
stability_functions
source_terms
KP_Q2
update_TKE
update_AL
second_eqn
Q_L_on_B
Q_L_on_B_Q2L

Called By

b3drun

Options - Logical

first

Options - Compiler

none

Known Issues

none


subroutine al_init

Description

Sets the initial value of the Mixing lenght in terms of $\sigma $ $(l_\sigma)$ al using algebraic formulation of Bakhmetev. But limited in case of deep water.

Numerical Equations


\begin{displaymath}
l_\sigma = \kappa \left( {1 - \sigma } \right)\left( { - \si...
...right)^{\frac{1}{{\begin{array}{l}
2 \\
\\
\end{array}}}}
\end{displaymath}

or in case of deep water:

\begin{displaymath}
\begin{array}{l}
l_\sigma = \kappa \left( {\frac{h_c }{h}} ...
...d \quad \quad
\sigma \le \frac{2h_c }{3h} - 1 \\
\end{array}\end{displaymath}

Subroutine Arguments

none

Local variables

i,j,k grid indices
tagaln Message passing tag
al1,al2,almax,hc max l, critical h

Global variables changed

aln

Logical units

none

Order of Things

  1. Sets the depth critical value.
  2. Sets maximum l
  3. Calculates length scale

Calls

Exch3DS,Exch3DR

Called By

turbulence_closure

Options - Logical

none

Options - Compiler

SCOORD

Known Issues

none


subroutine bulk_gradients

Description

Calculates the vertical shear and buoyancy gradients.

Numerical Equations


\begin{displaymath}
\begin{array}{l}
M^2 = H^{ - 2}\left[ {\left( {\frac{\parti...
...H^{ - 1}\frac{\partial b_0 }{\partial \sigma } \\
\end{array}\end{displaymath}

Subroutine Arguments

i,j Horizontal grid indeces

Local variables

none

Global variables changed

none

Logical units

none

Order of Things

  1. Calculates vertical buoyancy gradient $N^2$
  2. Calculates Shear

Called By

turbulence_closure

Options - Compiler

SCOORD
NO_CONVADJ

Known Issues

none


subroutine update_AL

Description

Updates the Mixing length using initial algebraic $\sigma $-formulation from al_init and converts it l by multiplying it by the local depth.

Subroutine Arguments

i,j Horizontal grid indeces.

Local variables

k Vertical grid index

Global variables changed

al

Order of Things

  1. Calculates al by multiplying aln by h.
  2. Limits al.

Calls

limitAL

Called By

turbulence_closure

Options - Compiler

SCOORD

Known Issues

none


subroutine limit_AL

Description

Setup a limiter for the mixing length whe stratification is present.

Numerical Equations


\begin{displaymath}
\begin{array}{l}
l = l\left( {\left( {Ri^2 + 1} \right)^{\f...
...al z}}}
\right)\quad \quad \quad \quad other. \\
\end{array}\end{displaymath}

Subroutine Arguments

none

Local variables

k vertical grid index
alim, Ri length scale scaling factor, Richardson Number.

Global variables changed

al

Order of Things

  1. Check stratification and shear
  2. Calculates limiting factor for mixing length.

Calls

none

Called By

update_AL
update_Q2L
Q_L_on_B_Q2L
update_Q2L

Options - Compiler

none

Known Issues

none


subroutine stability_functions

Description

Calculates the stability functions S_M and S_H

Subroutine Arguments

none

Local variables

k Vertical grid index
gg Gaperin Stability function

Global variables changed

none

Logical units

none

Order of Things

  1. Sets the boundary values of S_H and S_M.
  2. Calculates gg for the internal points.
  3. Calculates S_M and S_H for the internal points.

Called By

turbulence_closure

Options - Logical

none

Options - Compiler

none

Known Issues

none


subroutine Q_L_on_B

Description

Interpolates $q^2$ and $l$ to b-points.

Subroutine Arguments

i,j,iig,jjg Horizontal grid indeces.

Local variables

k vertical grid indeces.

Global variables changed

qsq,al Twice TKE and mixing scale.

Order of Things

  1. Calculates vertical buoyancy gradients and length scale in case of S-coordinates.
  2. If no convective adjsutment the buoyancy gradients are small they are set to zero.
  3. Interpolates qsq.
  4. interpolates al.
  5. Calculates sqs.
  6. Limits al.

Calls

limit_AL

Called By

turbulence_closure

Options - Compiler

SCOORD
NO_CONVADJ

Known Issues

none


subroutine KP_Q2

Description

Calculates the Eddy diffusivity for TKE

Subroutine Arguments

None

Local variables

k Vertical grid index
qq Eddy diffusivity for TKE

Global variables changed

none

Order of Things

  1. Sets value of $S_q$.
  2. Calculates Eddy diffusivity.

Calls

none

Called By

turbulence_closure

Options - Logical

none

Options - Compiler

none

Known Issues

none


subroutine update_TKE

Description

Solves the TKE differential equation. With the sources of order 0 s0 explicitly and the sources of order 1 s1 and diffusion imlicitly.

Subroutine Arguments

i,j Horizontal grdi indeces.

Local variables

k vertical grid.
hh depth squate
qk diffusion

Global variables changed

qsq

Order of Things

  1. Evaluates Bottom Boundary Condition
  2. Evaluates Surface Boundary Condition
  3. Set coefficients for (n-1)*(n-1) tridiagonal matrix: acr*qsq(k-1,i,j)+bcr*qsq(k,i,j)+ccr*qsq(k+1,i,j)=rcr
  4. Solves tridiagonal matrix for qsq
  5. Checks that qsq is greater than the minimum qsq

Calls

tridag

Called By

turbulence_closure

Options - Logical

none

Options - Compiler

SCOORD

Known Issues

none


subroutine second_eqn

Description

Solves the differential equation for $q^2l$

Numerical Equations


\begin{displaymath}
\frac{\partial q^2l}{\partial t} = 2E_1 l\left( {K_ZM^2 + K_...
...right) +
D\left( {q^2l} \right) - \frac{2q^3}{B_1 }\tilde {W}
\end{displaymath}

where $\tilde{W}$ is the parabolic wall proximity function:

\begin{displaymath}
\tilde {W} = 1 + E_2 \left( {\frac{l}{\kappa }\left( {\frac{1}{d_s } +
\frac{1}{d_b }} \right)} \right)^2
\end{displaymath}

Subroutine Arguments

i,j Horizontal grid indeces.

Local variables

k Vertical grid index

Global variables changed

al mixing length scale

Order of Things

  1. Calculates balance between production and dissipation
  2. Updates $q^2l$ by solving diferential equation including diffusion terms

Calls

source_terms_q2l
update_q2l

Called By

turbulence_closure

Options - Compiler

none

Known Issues

Not tested in anger


subroutine source_terms_q2l

Description

Calculates the changes of $q^2l$ due to production buoyancy and dissipation.

Numerical Equations


\begin{displaymath}
\frac{\partial q^2l}{\partial t} = 2E_1 l\left( {K_ZM^2 + K_...
...right) +
D\left( {q^2l} \right) - \frac{2q^3}{B_1 }\tilde {W}
\end{displaymath}

Subroutine Arguments

i,j Vertical grid indeces.

Local variables

k Vertical grid index
gamma Use to store the shear and buoyancy components.
sigop Use in the wall proximity function.

Global variables changed

none

Order of Things

  1. Calculates 0-order changes in $q^2l$.
  2. Calculates 1-order changes in $q^2l$.

Calls

none

Called By

second_equation

Options - Compiler

SCOORD

Known Issues

Not tested in anger


subroutine update_Q2L

Description

Solves the $q^2l$ differential equation. sl0 explicit, sl1 implicit.

Subroutine Arguments

i,j Horizontal indeces.

Local variables

k Vertical index.
hh Depth square.
qk Eddy diffusivity for turbulence.

Global variables changed

al

Order of Things

  1. Set coefficients for (n-1)*(n-1) tridiagonal matrix acr*qsq(k-1,i,j)+bcr*qsq(k,i,j)+ccr*qsq(k+1,i,j)=rcr
  2. Solves matrix.
  3. Calculates al
  4. limits al

Calls

tridiag
limit_AL

Called By

second_equation

Options - Compiler

SCOORD

Known Issues

Not used in Anger


subroutine Q_L_on_B_Q2L

Description

Interpolates $q^2$ and $l$ to b points.

Subroutine Arguments

i,j,iig,jjg Horizontal grid indeces.

Local variables

k vertical grid indeces.

Global variables changed

qsq,al Twice TKE and mixing scale.

Order of Things

  1. Calculates vertical buoyancy gradients and length scale in case of S-coordinates.
  2. If no convective adjsutment the buoyancy gradients are small they are set to zero.
  3. Interpolates qsq.
  4. interpolates al.
  5. Calculates sqs.
  6. Limits al.

Calls

limit_AL

Called By

turbulence_closure

Options - Compiler

SCOORD
NO_CONVADJ

Known Issues

none


subroutine aakcon

Description

Set the eddy viscosity and eddy diffusivity as constant

Subroutine Arguments

none

Local variables

i,j,k Grid indeces

Global variables changed

aa,ak Eddy viscosity and eddy diffusivity

Logical units

none

Order of Things

  1. Sets eddy diffusitity on wet points to 0.1 m$^2$s$^{-1}$ as default.
  2. Set eddy viscosity on wet points to 0.01 m$^2$s$^{-1}$ as default.

Calls

none

Called By

b3drun

Options - Logical

none

Options - Compiler

none

Known Issues

none


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