GFDL - Geophysical Fluid Dynamics Laboratory

The Hallberg Isopycnal Model (HIM)

by Robert Hallberg

Introduction


This program (HIM) simulates the ocean by numerically solving the Boussinesq
primitive equations in isopycnal vertical coordinates and general orthogonal horizontal
coordinates. These equations are horizontally discretized on an Arakawa C-grid.
There are a range of options for the physical parameterizations, from those most
appropriate to highly idealized models for studies of geophysical fluid dynamics
to a rich suite of processes appropriate for realistic ocean simulations. The thermodynamic
options range from an adiabatic model with fixed density layers to a model with
temperature and salinity as state variables and using a full nonlinear equation
of state. The uppermost few layers may be used to describe a bulk mixed layer, including
the effects of penetrating shortwave radiation. Either a split- explicit time stepping
scheme or a non-split scheme may be used for the dynamics, while the time stepping
may be split (and use different numbers of steps to cover the same interval) for
the forcing, the thermodynamics, and for the dynamics. Most of the numerics are
second order accurate in space. HIM can run with an absurdly thin minimum layer
thickness.

Details of the numerics and physical parameterizations are provided in the appropriate
source files. Most of the available options are selected by the settings in init.h,
although some (such as the equation of state) are selected by specifying which file
to use in the make file (Makefile).

The file HIM.c contains the main time stepping loops. One time integration
option for the dynamics uses a split explicit time stepping scheme to rapidly step
the barotropic pressure and velocity fields. The barotropic velocities are averaged
over the baroclinic time step before they are used to advect thickness and determine
the baroclinic accelerations. At the end of every time step, the free surface height
perturbation is determining by adding up the layer thicknesses; this perturbation
is used to drive the free surface heights from the barotropic calculation and from
the sum of the layer thicknesses toward each other over subsequent time steps. The
barotropic and baroclinic velocities are synchronized as part of the vertical viscosity
algorithm and be recalculating. the barotropic velocities from the baroclinic velocities
each time step. This scheme is described in Hallberg, 1997, J. Comp. Phys. 135,
54-65.

The other time integration option uses a non-split time stepping scheme based
on the 3-step third order Runge-Kutta scheme described in Matsuno, 1966, J. Met.
Soc. Japan, 44, 85-88. For problems with more than a very few layers, this non-split
scheme is much less efficient than the split scheme, but because it is much simpler
and formally more accurate, it is useful to have to verify that temporal truncation
errors are not significant.

There are a range of closure options available in HIM. Horizontal velocities
are subject to a combination of horizontal biharmonic and Laplacian friction (based
on a stress tensor formalism) and a vertical Fickian viscosity (perhaps using the
kinematic viscosity of water). The horizontal viscosities may be constant, spatially
varying or may be dynamically calculated using Smagorinsky’s approach. A diapycnal
diffusion of density and thermodynamic quantities is also allowed, but not required,
as is horizontal diffusion of interface heights (akin to the Gent-McWilliams closure
of geopotential coordinate models). The diapycnal mixing may use a fixed diffusivity
or it may use the shear Richardson number dependent closure described in Hallberg
(MWR, 2000). When there is diapycnal diffusion, it applies to momentum as well.
As this is in addition to the vertical viscosity, the vertical Prandtl always exceeds
1.

HIM has a number of noteworthy debugging capabilities. Excessively large velocities
are truncated and HIM will stop itself after a number of such instances to keep
the model from crashing altogether, and the model state is output with a reported
time of 9.9e9. This is useful in diagnosing failures, or (by accepting some truncations)
it may be useful for getting the model past the adjustment from an ill-balanced
initial condition. In addition, all of the accelerations in the columns with excessively
large velocities may be directed to a text file. Parallelization errors may be diagnosed
with the CHECK_PARALLEL option, whereby ostensibly identical model incarnations
are run simultaneously on one and multiple processors and any differences are reported.

About 35 other files of source code and 6 header files must be used to make HIM
work. A makefile is included to make compiling easy. Type "make HIM" to
compile. Some run time input is required, but this is prompted for. It may be convenient
to direct a small file containing the needed information to standard input, and
direct the output to another file.

Source File Subroutines


The ~35 source files contain the following subroutines:


    HIM.c

    step_HIM
    steps HIM over a specified interval of time.
    HIM_initialize
    calls initialize and does other initialization that does not warrant user
    modification.
    set_restart_fields
    is used to specify those fields that are written to and read from the restart
    file.
    calculate_surface_state
    determines the surface (mixed layer) properties of the current model state
    and packages pointers to these fields into an exported structure.

    HIM_checkparallel.c

    HIM_checkparallel.c is identical to HIM.c, except that it has a number
    of diagnostic calls to assess whether a parallel run is giving identical
    results to a single processor run.
    HIM_driver.c

    main
    is where HIM starts. Inside of main are the calls that set up the
    run, step the model, and orchestrate output and normal termination of the run.

Initial Condition, Forcing, and Domain Specification Routines



    initialize.c

    initialize
    does just that to all of the fields that are needed to specify the initial
    conditions of the model. initialize calls a number of other subroutines
    in initialize.c, each of which initializes a single field (or a few
    closely related fields) that are indicated by the subroutine name.

    initialize_output.c

    USER_set_output
    specifies what fields are output into which files and when, and whether
    they are averaged in time.
    USER_set_diagnostics
    specifies what diagnostic quantities are calculated.
    GetInputLines
    gets 4 controlling inputs from stdin, including the directories
    for I/O and the parameter specification file.
    write_grid_file
    writes out a file describing the model grid.

    surface_forcing.c

    set_forcing
    sets the current values of surface forcing fields.
    wind_forcing
    sets the current surface wind stresses.
    buoyancy_forcing
    sets the current surface heat, fresh water, buoyancy or other appropriate
    tracer fluxes.
    set_forcing_output
    sets up the output of any forcing fields.
    average_forcing
    accumulates time averages of indicated forcing fields.
    register_forcing_restarts
    is used to specify the forcing-related fields that are written to and read
    from the restart file.

    set_metrics.c

    set_metrics
    calculates the horizontal grid spacings and related metric fields, along
    with the grid point locations.
    initialize_masks
    initializes the land masks.

Principal Dynamic Routines



    CoriolisAdv.c

    CorAdCalc
    calculates the Coriolis and advective accelerations.

    PressureForce.c

    PressureForce
    calculates the pressure acceleration.
    register_compress
    is used to specify the reference profile of potential temperature and salinity
    that is used to compensate for compressibility.
    uncompress_e_rho
    makes internally consistent changes to profiles of interface height and
    density to offset compressibility and to minimize the non-solenoidal pressure
    gradient term.

    continuity.c, continuity_mpdata.c, or continuity_FCT.c

    continuity time steps the layer thicknesses.

    barotropic.c

    btstep
    time steps the linearized barotropic equations for use with the split explicit
    time stepping scheme.
    btcalc
    calculates the barotropic velocities from the layer velocities.
    barotropic_init
    initializes several split-related variables and calculates several static
    quantities for use by btstep.
    register_barotropic_restarts
    indicates those time splitting-related fields that are to be in the restart
    file.

    hor_visc.c

    horizontal_viscosity
    calculates the convergence of momentum due to Laplacian or biharmonic horizontal
    viscosity.
    set_up_hor_visc
    calculates combinations of metric coefficients and other static quantities
    used in horizontal_viscosity.

    vertvisc.c

    vertvisc
    changes the velocity due to vertical viscosity, including application of
    a surface stress and bottom drag.
    set_viscous_BBL
    determines the bottom boundary layer thickness and viscosity according to
    a linear or quadratic drag law.

    thickness_diffuse.c

    thickness_diffuse
    moves fluid adiabatically to horizontally diffuse interface heights.

Thermodynamic Routines



    diabatic_driver.c

    diabatic
    orchestrates the calculation of vertical advection and diffusion of momentum
    and tracers due to diapycnal mixing and mixed layer (or other diabatic) processes.
    mixedlayer, Calculate_Entrainment, apply_sponge,
    and any user-specified tracer column physics routines are all called by
    diabatic
    .

    diabatic_entrain.c

    Calculate_Entrainment
    calculates the diapycnal mass fluxes due to interior diapycnal mixing processes,
    which may include a Richardson number dependent entrainment.
    Calculate_Rino_flux
    estimates the Richardson number dependent entrainment in the absence of
    interactions between layers, from which the full interacting entrainment can
    be found.
    Estimate_u_h
    estimates what the velocities at thickness points will be after entrainment.

    mixed_layer.c

    mixed_layer
    implements a bulk mixed layer, including entrainment and detrainment, related
    advection of dynamically active tracers, and buffer layer splitting. The bulk
    mixed layer may consist of several layers.

    tracer.c

    register_tracer
    is called to indicate a field that is to be advected by advect_tracer
    and diffused by tracer_hordiff.
    advect_tracer
    does along-isopycnal advection of tracer fields.
    tracer_hordiff
    diffuses tracers along isopycnals.
    register_tracer_init_fn
    registers a user-specified tracer initialization subroutine.
    call_tracer_init_fns
    calls any user-specified tracer initialization subroutines that have been
    registered.
    register_tracer_column_fn
    registers a user-specified tracer column processes subroutine.
    call_tracer_column_fns
    calls any user-specified tracer column processes subroutines that have been
    registered.

    sponge.c

    apply_sponge
    damps fields back to reference profiles.
    initialize_sponge
    stores the damping rates and allocates the memory for the reference profiles.
    set_up_sponge_field
    registers reference profiles and associates them with the fields to be damped.

    eqn_of_state.c, eqn_of_state_linear.c, or eqn_of_state_UNESCO.c

    calculate_density
    calculates a list of densities at given potential temperatures, salinities
    and pressures.
    calculate_density_derivs
    calculates a list of the partial derivatives with temperature and salinity
    at the given potential temperatures, salinities and pressures.
    calculate_compress
    calculates a list of the compressibilities (partial derivatives of density
    with pressure) at the given potential temperatures, salinities and pressures.
    calculate_2_densities
    calculates a list of the densities at two specified reference pressures
    at the given potential temperatures and salinities.

    fit_compressibility.c

    fit_compressibility
    determines the best fit of compressibility with pressure using a fixed 5-coefficient
    functional form, based on a provided reference profile of potential temperature
    and salinity with depth. This fit is used in PressureForce.

Infrastructural Routines



    HIM_restart.c

    save_restart
    saves a restart file (or multiple files if they would otherwise be too large).
    register_restart_field
    is called to specify a field that is to written to and read from the restart
    file.
    restore_state
    reads the model state from output or restart files.
    query_initialized
    indicates whether a specific field or all restart fields have been read
    from the restart files.

    HIM_parser.c

    HIM_parser
    parses a parameter specification file with the same format as init.h
    to (possibly) change parameters at run time.

    HIM_sum_output.c

    write_energy
    writes the layer energies and masses and other spatially integrated quantities
    and monitors CPU time use.
    depth_list_setup
    generates a list of the volumes of fluid below various depths.

    HIM_field_output.c

    specify_saves_file
    sets up a file for output, including when and how frequently the file is
    to be written.
    set_up_save_field
    sets up a field to be saved and also perhaps to be averaged in time.
    save_fields
    saves all fields that are to be saved within the given time interval and
    indicates the next scheduled time that a field is to be saved.
    enable_averaging
    enables averaging for a time interval.
    disable_averaging
    disables the accumulation of averages.
    query_averaging_enabled
    indicates whether averaging is currently enabled.
    get_field_ref
    returns a reference integer for a field that is to be averaged, or 0 if
    it is not to be averaged.
    average_field
    accumulates the average of a field.

    timetype.c (Partial C implementation of the FMS time
    manager.)

    t_ge, t_gt, t_le, t_lt, t_eq, and t_ne
    make logical comparisons between times.
    t_plus
    adds two times together.
    t_minus
    finds the absolute value of the difference between times. (i.e. |t1-t2|)
    t_mult
    multiplies a time by an integer.
    t_div
    returns the real ratio of two times.
    t_div_int
    divides a time by an integer.
    t_interval
    finds the signed, real ratio of the interval between two times to a third
    time. (i.e. (t1-t2)/t3)
    double_to_time
    converts a real number to a time.
    time_to_double
    converts a time to a real number.

    HIM_io.c (Input/Output utility subroutines)

    create_file
    creates a new file, set up structures that are needed for subsequent output,
    and write the coordinates.
    reopen_file
    reopens an existing file for writing and set up structures that are needed
    for subsequent output.
    open_input_file
    opens the indicated file for reading only.
    close_file
    closes an open file.
    synch_file
    flushes the buffers for a file, completing all pending output.
    write_field
    writes a field to an open file.
    write_time
    writes a value of the time axis to an open file.
    Read_Field
    provides a simpler interface to read_field.
    read_field
    reads a field from an open file.
    read_time
    reads a time from an open file.
    name_output_file
    provides a name for an output file based on a name root and the time of
    the output.
    find_input_file
    finds a file that has been previously written by HIM and named by name_output_file
    and opens it for reading.
    handle_error
    writes an error code and quits.

    HIM_parallel.c (Parallelization utility subroutines)

    set_up_parallel
    starts parallel runs and calculates subdomain sizes and neighboring domains.
    pass_var
    passes a 3-D variable to neighboring processors or applies corresponding
    boundary conditions.
    pass_var_2d
    passes a 2-D variable to neighboring processors or applies corresponding
    boundary conditions.
    sum_fields
    sums fields across all processors.
    sum_int_fields
    sums integer fields across all processors.
    collect_double_fields
    collects a full layer of a field to processor 0.
    spread_string
    distributes a character string from processor 0.
    spread_double_vector
    distributes copies of a 1-D array from processor 0.
    check_field
    checks whether a field is identical on different processors, and reports
    any differences. This is useful for checking parallelization or for debugging
    coding changes that are not expected to change answers.
    quit
    causes all processors to stop execution.

Purely Diagnostic Routines



    diagnostics.c

    calculate_diagnostic_fields
    is used to calculate several diagnostic fields that are not naturally calculated
    elsewhere.
    register_time_deriv
    is used to register the information needed for diagnostically calculating
    a time derivative.
    calculate_derivs
    calculates any registered time derivatives.

    PointAccel.c

    write_u_accel
    writes a long list of zonal accelerations and related quantities for one
    column out to a file. This is typically called for diagnostic purposes from
    vertvisc when a zonal velocity exceeds the specified threshold.
    write_v_accel
    writes a long list of meridional accelerations and related quantities for
    one column out to a file. This is typically called for diagnostic purposes from
    vertvisc when a meridional velocity exceeds the specified threshold.

Header Files


In addition there are 5 header files:


    init.h

    sets various constants and parameters for the simulation.

    HIM.h

    contains the subroutine prototypes and the definitions of the following
    6 structures:
    params (HIM_params) contains changeable simulation parameters.
    forcing (fluxes) contains pointers to surface forcing fields.
    surface (state) contains pointers to ocean surface fields.
    diag_fld (diag) contains pointers to diagnostic quantities.
    thermo_var_ptrs (tv) contains pointers to thermodynamic fields,
    such as potential temperature and salinity.
    bt_vars_ptrs (bt) contains pointers to variables related to the
    baroclinic/barotropic split time stepping.

    HIM_io.h

    contains the subroutine prototypes and structure definitions for I/O
    related subroutines.

    metrics.h

    contains the descriptions for a number of metric terms.

    hor_visc.h

    contains the descriptions for a number of metric-related fields that
    are only used by the horizontal viscosity.

Most simulations can be set up by modifying only the files init.h,
initialize.c, and surface_forcing.c. In addition, the file
initialize_output.c
will commonly be modified to tailor the output to the needs
of the question at hand. These altered files and the makefile might reside in the
same directory as the executable file. All of the other (unaltered) source code
should probably remain in some central directory. The makefile will, of course,
need to reflect the locations of the desired source files.