The Hallberg Isopycnal Model (HIM)
by Robert Hallberg
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
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,
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
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:
- steps HIM over a specified interval of time.
- calls initialize and does other initialization that does not warrant user
- is used to specify those fields that are written to and read from the restart
- determines the surface (mixed layer) properties of the current model state
and packages pointers to these fields into an exported structure.
- 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.
- 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
- 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.
- specifies what fields are output into which files and when, and whether
they are averaged in time.
- specifies what diagnostic quantities are calculated.
- gets 4 controlling inputs from stdin, including the directories
for I/O and the parameter specification file.
- writes out a file describing the model grid.
- sets the current values of surface forcing fields.
- sets the current surface wind stresses.
- sets the current surface heat, fresh water, buoyancy or other appropriate
- sets up the output of any forcing fields.
- accumulates time averages of indicated forcing fields.
- is used to specify the forcing-related fields that are written to and read
from the restart file.
- calculates the horizontal grid spacings and related metric fields, along
with the grid point locations.
- initializes the land masks.
Principal Dynamic Routines
- calculates the Coriolis and advective accelerations.
- calculates the pressure acceleration.
- is used to specify the reference profile of potential temperature and salinity
that is used to compensate for compressibility.
- makes internally consistent changes to profiles of interface height and
density to offset compressibility and to minimize the non-solenoidal pressure
- continuity time steps the layer thicknesses.
continuity.c, continuity_mpdata.c, or continuity_FCT.c
- time steps the linearized barotropic equations for use with the split explicit
time stepping scheme.
- calculates the barotropic velocities from the layer velocities.
- initializes several split-related variables and calculates several static
quantities for use by btstep.
- indicates those time splitting-related fields that are to be in the restart
- calculates the convergence of momentum due to Laplacian or biharmonic horizontal
- calculates combinations of metric coefficients and other static quantities
used in horizontal_viscosity.
- changes the velocity due to vertical viscosity, including application of
a surface stress and bottom drag.
- determines the bottom boundary layer thickness and viscosity according to
a linear or quadratic drag law.
- moves fluid adiabatically to horizontally diffuse interface heights.
- 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
- calculates the diapycnal mass fluxes due to interior diapycnal mixing processes,
which may include a Richardson number dependent entrainment.
- estimates the Richardson number dependent entrainment in the absence of
interactions between layers, from which the full interacting entrainment can
- estimates what the velocities at thickness points will be after entrainment.
- 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.
- is called to indicate a field that is to be advected by advect_tracer
and diffused by tracer_hordiff.
- does along-isopycnal advection of tracer fields.
- diffuses tracers along isopycnals.
- registers a user-specified tracer initialization subroutine.
- calls any user-specified tracer initialization subroutines that have been
- registers a user-specified tracer column processes subroutine.
- calls any user-specified tracer column processes subroutines that have been
- damps fields back to reference profiles.
- stores the damping rates and allocates the memory for the reference profiles.
- registers reference profiles and associates them with the fields to be damped.
- calculates a list of densities at given potential temperatures, salinities
- calculates a list of the partial derivatives with temperature and salinity
at the given potential temperatures, salinities and pressures.
- calculates a list of the compressibilities (partial derivatives of density
with pressure) at the given potential temperatures, salinities and pressures.
- calculates a list of the densities at two specified reference pressures
at the given potential temperatures and salinities.
eqn_of_state.c, eqn_of_state_linear.c, or eqn_of_state_UNESCO.c
- 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.
- saves a restart file (or multiple files if they would otherwise be too large).
- is called to specify a field that is to written to and read from the restart
- reads the model state from output or restart files.
- indicates whether a specific field or all restart fields have been read
from the restart files.
- parses a parameter specification file with the same format as init.h
to (possibly) change parameters at run time.
- writes the layer energies and masses and other spatially integrated quantities
and monitors CPU time use.
- generates a list of the volumes of fluid below various depths.
- sets up a file for output, including when and how frequently the file is
to be written.
- sets up a field to be saved and also perhaps to be averaged in time.
- 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.
- enables averaging for a time interval.
- disables the accumulation of averages.
- indicates whether averaging is currently enabled.
- returns a reference integer for a field that is to be averaged, or 0 if
it is not to be averaged.
- accumulates the average of a field.
- t_ge, t_gt, t_le, t_lt, t_eq, and t_ne
- make logical comparisons between times.
- adds two times together.
- finds the absolute value of the difference between times. (i.e. |t1-t2|)
- multiplies a time by an integer.
- returns the real ratio of two times.
- divides a time by an integer.
- finds the signed, real ratio of the interval between two times to a third
time. (i.e. (t1-t2)/t3)
- converts a real number to a time.
- converts a time to a real number.
timetype.c (Partial C implementation of the FMS time
- creates a new file, set up structures that are needed for subsequent output,
and write the coordinates.
- reopens an existing file for writing and set up structures that are needed
for subsequent output.
- opens the indicated file for reading only.
- closes an open file.
- flushes the buffers for a file, completing all pending output.
- writes a field to an open file.
- writes a value of the time axis to an open file.
- provides a simpler interface to read_field.
- reads a field from an open file.
- reads a time from an open file.
- provides a name for an output file based on a name root and the time of
- finds a file that has been previously written by HIM and named by name_output_file
and opens it for reading.
- writes an error code and quits.
HIM_io.c (Input/Output utility subroutines)
- starts parallel runs and calculates subdomain sizes and neighboring domains.
- passes a 3-D variable to neighboring processors or applies corresponding
- passes a 2-D variable to neighboring processors or applies corresponding
- sums fields across all processors.
- sums integer fields across all processors.
- collects a full layer of a field to processor 0.
- distributes a character string from processor 0.
- distributes copies of a 1-D array from processor 0.
- 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.
- causes all processors to stop execution.
HIM_parallel.c (Parallelization utility subroutines)
Purely Diagnostic Routines
- is used to calculate several diagnostic fields that are not naturally calculated
- is used to register the information needed for diagnostically calculating
a time derivative.
- calculates any registered time derivatives.
- 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.
- 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.
In addition there are 5 header files:
- sets various constants and parameters for the simulation.
- contains the subroutine prototypes and the definitions of the following
- 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.
- contains the subroutine prototypes and structure definitions for I/O
- contains the descriptions for a number of metric terms.
- 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.