PUBLIC INTERFACE ~ PUBLIC DATA ~ PUBLIC ROUTINES ~ NAMELIST ~ DIAGNOSTIC FIELDS ~ ERROR MESSAGES ~ REFERENCES ~ NOTES

Module ocean_neutral_physics_mod

Contact:  Stephen M. Griffies Russell Fiedler
Reviewers: 
Change History: WebCVS Log


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

ocean_neutral_physics_init:
neutral_physics:
neutral_slopes:
neutral_blayer:
fz_terms:
fx_flux:
fy_flux:
fz_flux:
compute_eady_rate:
compute_baroclinicity:
compute_rossby_radius:
compute_bczone_radius:
compute_diffusivity:
gm_velocity:
slope_function_gm:
ocean_neutral_physics_end:


PUBLIC DATA

None.


PUBLIC ROUTINES

  1. 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.


  2. 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.



  3. 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.


  4. 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.



  5. 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.


  6. 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)


  7. 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)


  8. 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)


  9. compute_eady_rate

    DESCRIPTION
    Finish computing eady growth rate.


  10. compute_baroclinicity

    DESCRIPTION
    Finish computing baroclinicity, which is defined to be the vertically averaged magnitude of the horizontal density gradient.


  11. 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).


  12. 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.


  13. compute_diffusivity

    DESCRIPTION
    Subroutine computes flow dependent diffusivity. Allow for an added dimensionless tuning factor as well as a minimum and maximum diffusivity.


  14. 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.


  15. slope_function_gm

    DESCRIPTION
    Function for defining effective slope in diagnostic GM velocity calculation. Used only for diagnostic purposes.


  16. 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

  1. 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
  2. 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
  3. 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
  4. S.M. Griffies The Gent-McWilliams Skew-flux Journal of Physical Oceanography (1998) vol 28 pages 831-841
  5. S.M. Griffies Fundamentals of Ocean Climate Models (2003) NOAA/Geophysical Fluid Dynamics Laboratory
  6. 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
  7. 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.


top