Module ocean_neutral_physics_mod
OVERVIEW
Thickness weighted time tendency for tracer
from Laplacian neutral diffusion + Laplacian GM skew-diffusion.
This module computes the cell thickness weighted tracer tendency from
small angle Laplacian neutral diffusion + Laplacian GM skew-diffusion.
It uses the "triad" algorithm also coded in MOM2.2, MOM3, and MOM3.1.
The MOM4 algorithm accounts for partial bottom cells and generalized
orthogonal horizontal coordinates.
OTHER MODULES USED
constants_mod
diag_manager_mod
fms_mod
mpp_domains_mod
mpp_mod
time_manager_mod
ocean_density_mod
ocean_domains_mod
ocean_operators_mod
ocean_sigma_diffuse_mod
ocean_types_mod
ocean_workspace_mod
PUBLIC INTERFACE
PUBLIC DATA
None.
PUBLIC ROUTINES
-
ocean_neutral_physics_init
-
DESCRIPTION
- Initialize the neutral physics module by registering fields for
diagnostic output and performing some numerical checks to see
that namelist settings are appropriate.
-
neutral_physics
-
DESCRIPTION
- This function computes the thickness weighted time tendency for tracer
from neutral physics. Full discussion and details are provided by
Griffies (2004). Here is a brief summary.
---How the neutral diffusive flux components are computed:
The vertical flux component is split into diagonal (3,3) and
off-diagonal (3,1) and (3,2) terms. The off-diagonal (3,1) and (3,2)
terms are included explicitly in time. The main contribution from the
(3,3) term to the time tendency is included implicitly in time
along with the usual contribution from diapycnal processes
(vertical mixing schemes). This is the K33_implicit term.
This approach is necessary with high vertical resolution, as
noted by Cox (1987). However, splitting the vertical flux into
an implicit and explicit piece compromises the
integrity of the vertical flux component (see Griffies et al. 1998).
So to minimize the disparity engendered by this split, the portion of
K33 that can be stably included explicitly in time is computed along
with the (3,1) and (3,2) terms.
All other terms in the mixing tensor are included explicitly in time.
The off-diagonal terms in the horizontal flux components, and all terms
in the vertical flux component, are tapered in regions of steep neutral
slope according to the requirements of linear stability. MOM4 allows for
choice of two tapering schemes:
(a) the tanh taper of Danabasoglu and McWilliams (1995)
(b) the quadratic scheme of Gerdes, Koberle, and Willebrand (1991)
Linear stability is far less stringent on the diagonal (1,1) and (2,2)
part of the horizontal flux. Indeed, these terms in practice need
not be tapered in steep sloped regions. The namelist
neutral_taper_diagonal=.false. keeps the diagnonal terms maintained
for all neutral slopes. This approach assists in reducing numerical
noise in regions where the physical system experiences a lot of
diapycnal mixing anyhow.
---How the skew diffusive flux components are computed:
The GM skew flux components are purely off-diagonal.
They are generally tapered when neutral slope
is large (neutral_physics_simple=.false).
Doing so maintains a nontrivial GM slumping effect even when the
neutral slopes are vertical. The alternative neutral_physics_simple=.true.
is the approach used in MOM3, whereby GM effects are removed
in steep sloped regions. neutral_physics_simple=.false. is
less efficient, but has been seen to yield superior simulations.
Comments: Much research as of 2003-2004 aims to clarify the
schemes used to meld neutral physics with physics occuring
in steep sloped regions (i.e., boundary layers). Methods used
in mom4 are not expected to be the final word on these matters.
-
neutral_slopes
-
DESCRIPTION
- Subroutine computes the neutral slopes for the triads associated
with the vertical flux component. No tapering applied here.
Array tensor_31 initially holds the x-slope used for flux component fz.
Array tensor_32 initially holds the y-slope used for flux component fz.
-
neutral_blayer
-
DESCRIPTION
- Subroutine computes the boundary layer as determined by
1. steep neutral slopes
2. depth within which typical mesoscale eddies are partially outcropped
3. depth within which vertical mixing scheme (e.g., kpp) computes a boundary layer
Note: Only consider surface boundary layers here.
-
fz_terms
-
DESCRIPTION
- Subroutine computes the tracer independent pieces of the vertical
flux component. As a result of this routine,
Array tensor_31 = x-diffusivity*slope (m^2/sec) for fz
Array tensor_32 = y-diffusivity*slope (m^2/sec) for fz
K33 is the (3,3) term in small angle Redi diffusion tensor.
It is broken into an explicit in time piece and implicit
in time piece.
K33 has units m^2/sec.
Also will compute the squared Eady growth rate, with the maximum slope
contributing to this growth rate set by smax.
-
fx_flux
-
DESCRIPTION
- Subroutine computes the zonal neutral physics tracer flux component.
Compute this component for all tracers at level k.
fx has physical dimensions (area*diffusivity*tracer gradient)
-
fy_flux
-
DESCRIPTION
- Subroutine computes the meridional neutral physics tracer flux component.
Compute this component for all tracers at level k.
fy has physical dimensions (area*diffusivity*tracer gradient)
-
fz_flux
-
DESCRIPTION
- Subroutine computes the vertical neutral physics tracer flux component.
Compute this component for all tracers at level k.
Surface and bottom boundary condition fz(k=0)=fz(k=kmt(i,j))=0
fz has physical dimensions (diffusivity*tracer gradient) (usual flux units)
-
compute_eady_rate
-
DESCRIPTION
- Finish computing eady growth rate.
-
compute_baroclinicity
-
DESCRIPTION
- Finish computing baroclinicity, which is defined to be the vertically
averaged magnitude of the horizontal density gradient.
-
compute_rossby_radius
-
DESCRIPTION
- Subroutine computes the first baroclinic Rossby radius of deformation.
Employ WKB approach described by Chelton et al. In particular,
use formulae (2.2), (2.3a) and (2.3b) from their paper.
Place a max and min value on the Rossby radius, as such is useful
for computing the diffusivity as well as for the sine_taper scheme.
Compute buoyancy frequency in terms of vertical gradient of
locally referenced potential density. Place the reference point
at the interface between the tracer cells, which is also where
the vertical derivative of neutral density is located. This amounts
to a centered difference computation similar to that used by
Chelton et al. equation (B.4).
-
compute_bczone_radius
-
DESCRIPTION
- Subroutine computes the radius of the baroclinic zone in a manner
suggested by the Hadley Centre approach (Malcolm Roberts, personal
communication). Algorithm was used in MOM3.
-
compute_diffusivity
-
DESCRIPTION
- Subroutine computes flow dependent diffusivity.
Allow for an added dimensionless tuning factor as well as a
minimum and maximum diffusivity.
-
gm_velocity
-
DESCRIPTION
- Subroutine computes GM eddy-induced velocity field for diagnostics.
Compute ustar and vstar at U-cell point, and wstar at T-cell bottom.
Do a two-point average rather than more democratic four-point avg
in order to avoid having to call mpp_update domains on tensor_31 and
tensor_32. The 0.5 factor is due to the two-point average.
Note that this algorithm is ad hoc. Researchers interested in this
field may wish to test alternatives.
-
slope_function_gm
-
DESCRIPTION
- Function for defining effective slope in diagnostic GM velocity
calculation. Used only for diagnostic purposes.
-
ocean_neutral_physics_end
-
DESCRIPTION
- Write to restart for case when running with agm_closure
NAMELIST
&ocean_neutral_physics_nml
-
neutral_physics_on
Must be true to use this module.
[logical]
-
neutral_physics_debug
For printing starting and ending checksums for restarts
[logical]
-
diffusion_all_explicit
To compute all contributions from neutral diffusion explicitly in time, including
the K33 diagonal piece. This approach is available only when have small time
steps and/or running with just a single tracer. It is for testing purposes.
[logical, units: dimensionless]
-
neutral_physics_limit
When tracer falls outside a specified range, revert to horizontal
diffusive fluxes at this cell. This is an ad hoc and incomplete attempt
to maintain monotonicity with the neutral physics scheme.
[logical]
-
tmask_neutral_on
If .true. then this logical reduces the neutral fluxes to
horizontal/vertical diffusion next to boundaries.
This approach has been found to reduce spurious
extrema resulting from truncation of triads used to compute
a neutral flux component.
[logical]
-
dm_taper
Set to true to use the tanh tapering scheme of Danabasoglu and McWilliams.
Default is true.
[logical]
-
gkw_taper
Set to true to use the quadradic tapering scheme of Gerdes, Koberle, and Willebrand.
Default is false.
[logical]
-
smax
Value of the maximum neutral direction slope above which the neutral fluxes are
either tapered to zero or saturated. Typical value is smax=0.01 or smaller.
[real]
-
swidth
Width in slope over which use tanh to taper fluxes in steep sloped regions.
Typical value swidth=0.1*smax
[real]
-
aredi
Neutral diffusivity used for experiments using a constant diffusivity.
[real]
-
agm
GM-skew diffusivity used for experiments using a constant diffusivity.
[real]
-
agm_lat_zones
If true, will set agm_array as constant within two latitudinal zones.
The idea is that one may wish to use a larger agm in the ACC than
elsewhere.
[logical]
-
agm_lat_zones_boundary
Boundary between agm in the south and north zones.
[real]
-
agm_lat_zones_ratio
Ratio between the large agm used in the southern latitudinal zone
to that used in the north.
agm_array(north) = agm
agm_array(south) = agm*agm_lat_zones_ratio
[real]
-
tracer_mix_micom
If .true., then the GM-skew diffusivity is set according to a velocity scale
times the grid spacing.
[logical]
-
vel_micom
Velocity scale that is used for computing the MICOM diffusivity.
[real, units: m/sec]
-
neutral_horz_mix_bdy
If .true., then use a horizontal diffusivity in the neutral boundary layer.
[logical]
-
vel_micom_bdy
Velocity scale that is used for computing the MICOM horizontal diffusivity
within the neutral boundary layer.
[real, units: m/sec]
-
ah_bdy
Constant horizontal diffusivity for the boundary layer.
[real, units: m^2/sec]
-
bryan_lewis_aredi
Set bryan_lewis_aredi=.true. when want to have aredi a function of depth
according to the Bryan and Lewis (1979) profile.
[logical]
-
ahs
ahs = adjustable parameter at the surface for bryan_lewis_aredi
[real]
-
ahb
ahb = adjustable parameter at the bottom for bryan_lewis_aredi
[real]
-
agm_closure
If .true. then will compute the GM-skew diffusivity as a function of the flow.
The length scale is determined by the Rossby radius and the time scale is
determined by the Eady growth rate. Diffusivities are depth independent.
[logical]
-
agm_closure_max
Maximum GM diffusivity allowed when using agm_closure=.true.
[real, units: m^2/sec]
-
agm_closure_min
Minimum GM diffusivity allowed when using agm_closure=.true.
[real, units: m^2/sec]
-
agm_closure_scaling
Dimensionless tuning parameter for computing flow dependent diffusivities.
[logical, units: dimensionless]
-
agm_closure_upper_depth
Upper depth where start the depth integration to compute the Eady
growth rate and/or baroclinicity.
[real, units: m]
-
agm_closure_lower_depth
Deeper depth where finish the depth integration to compute the Eady
growth rate and/or baroclinicity.
[real, units: m]
-
agm_closure_baroclinic
For computing the agm coefficient using only the vertically
averaged magnitude of the horizontal density gradient
(i.e., the "baroclinicity").
[logical]
-
agm_closure_buoy_freq
For computing the agm coefficient using only the vertically
averaged horizontal density gradient, use a buoyancy frequency,
which is fixed for all space-time.
[real, units: sec^-1]
-
agm_closure_length_fixed
Used fixed length scale for computing agm_closure diffusivity
[logical]
-
agm_closure_length
Fixed length scale for use with agm_closure_fixed_length
[real, units: meter]
-
agm_closure_length_rossby
For computing the agm_closure length scale according to Rossby radius.
[logical]
-
agm_closure_length_bczone
For computing the agm_closure length scale according to radius of baroclinic zone.
[logical]
-
bczone_max_pts
Max number of horizontal grid points for use in computing the baroclinic zone radius.
[integer]
-
agm_closure_bczone_crit_rate
Critical growth rate for determining width of the baroclinic zone.
[real, units: sec^-1]
-
agm_closure_growth_scale
Dimensionless scaling used to set a maximum for agm_growth.
[real, units: dimensionless]
-
rossby_radius_max
Maximum Rossby Radius used for agm_closure_length_rossby and
the sine_taper scheme. Default = 100e3 m.
[real, units: meter]
-
rossby_radius_min
Minimum Rossby Radius used for agm_closure_length_rossby and
the sine_taper scheme. Default = 15e3 m.
[real, units: meter]
-
neutral_physics_simple
If .true. then must have aredi=agm. The horizontal flux components are then
computed as horizontal downgradient diffusive fluxes regardless the neutral slope.
This approach improves algorithm efficiency. However, a .true. setting precludes
being able to have the GM-skew fluxes remain active in the steep sloped regions,
thus shutting off their effects to reduce the slopes of isopycnals in convective and
mixed layer regimes.
[logical]
-
neutral_linear_gm_taper
If .true. then with neutral_physics_simple=.false., will linearly taper GM
skew fluxes towards the surface within regions of steep neutral slopes.
This approach leads to a constant horizontal eddy-induced velocity in
the steeply sloping regions.
[logical]
-
neutral_taper_diagonal
For cases with neutral_physics_simple=.false., then neutral_taper_diagonal=.true.
will taper the diagonal pieces of the horizontal flux components when neutral slopes
are steep. With neutral_taper_diagonal=.false., then the horizontal flux components will
remain enabled for all slopes, thus producing horizontal downgradient diffusion in
regions of vertical neutral directions.
[logical, units: dimensionless]
-
neutral_sine_taper
If .true. then with neutral_physics_simple=.false., will apply a sine-taper
to GM and neutral diffusive fluxes in regions where the penetration depth
of eddies is deeper than the grid point. This method is essential for
fine vertical resolution grids.
[logical]
-
gm_velocity_save
If .true. then will compute the GM eddy-induce advection velocity components
for purposes of diagnostics. This computation is nontrivial as it requires
some extra calls to mpp_update_domains. The algorithm used to compute the
velocities may require modifications for efficiency, and suggestions are welcome.
[logical]
-
neutral_blayer_diagnose
Diagnose properties of the neutral physics boundary layer, whether have
neutral_linear_gm_taper or neutral_sine_taper true or not.
[logical]
DATA SETS
None.
ERROR MESSAGES
None.
REFERENCES
- D.B. Chelton, R.A. deSzoeke, M.G. Schlax, K.E. Naggar, N. Siwertz
Geographical Variability of the First Baroclinic Rossby Radius of Deformation
Journal of Physical Oceanography (1998) vol 28 pages 433-460
- G. Danabasoglu and J. C. McWilliams
Sensitivity of the global ocean circulation to
parameterizations of mesoscale tracer transports
Journal of Climate (1995) vol 8 pages 2967--2987
- S.M. Griffies, A. Gnanadesikan, R.C. Pacanowski, V. Larichev,
J.K. Dukowicz, and R.D. Smith
Isoneutral diffusion in a z-coordinate ocean model
Journal of Physical Oceanography (1998) vol 28 pages 805-830
- S.M. Griffies
The Gent-McWilliams Skew-flux
Journal of Physical Oceanography (1998) vol 28 pages 831-841
- S.M. Griffies
Fundamentals of Ocean Climate Models (2003)
NOAA/Geophysical Fluid Dynamics Laboratory
- I.M. Held and V.D. Larichev
A scaling theory for horizontally homogeneous baroclinically
unstable flow on a beta plane
Journal of Atmospheric Sciences (1996) vol 53 pages 946-952
- M. Visbeck, J.C. Marshall, T. Haine and M. Spall
Specification of eddy transfer coefficients in coarse resolution ocean
circulation models
Journal of Physical Oceanography (1997) vol 27 pages 381--402
COMPILER SPECIFICS
None.
PRECOMPILER OPTIONS
None.
LOADER OPTIONS
None.
TEST PROGRAM
None.
KNOWN BUGS
None.
NOTES
This module represents perhaps the most complicated piece of code in
mom4. It is strongly suggested that some time be taken prior to
making major changes. If you find problems, either in documentation,
algorithm, or implementation, please feel free to contact Griffies.
Numerical implementation follows the triad approach documented in the
references and implemented in MOM2 and MOM3.
Diffusivities can be flow dependent, where the length scale is related
to the first baroclinic Rossby radius and the time scale to the Eady
growth rate. Within 5deg of equator, use is made of the equatorial
Rossby radius. Length scale can also be determined by the width of
the baroclinic zone, as done in the Hadley Centre model.
If the namelist settings agm and aredi are not equal,
then flow dependent diffusivities will be computed ONLY for the GM
skew-fluxes.
neutral_physics_simple=.true. requires the neutral diffusive
diffusivity to be set equal to the GM skew-diffusive diffusivity
(aredi=agm). neutral_physics_simple=.true. results in down-gradient
horizontal flux components. This setting reduces the overall cost
of the neutral physics scheme significantly.
In steep slope regions, neutral fluxes are generally tapered to zero
with the tanh taper of Danabasoglu and McWilliams (1995) or the
quadratic scheme of Gerdes, Koberle, and Willebrand. However, if
neutral_physics_simple=.false., the GM skew-diffusive fluxes
can remain nonzero if have neutral_linear_gm_taper=.true.
FUTURE PLANS
None.