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:
- 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 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.
- 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
- 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.
- 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.
- 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
- calculates the horizontal grid spacings and related metric fields, along with the grid point locations.
- initialize_masks
- initializes the land masks.
Principal Dynamic Routines
- CorAdCalc
- calculates the Coriolis and advective accelerations.
- 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 time steps the layer thicknesses.
- 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.
- 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
- 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
- moves fluid adiabatically to horizontally diffuse interface heights.
Thermodynamic Routines
- 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.
- 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
- 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.
- 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.
- 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.
- 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
- 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
- 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
- parses a parameter specification file with the same format as init.h to (possibly) change parameters at run time.
- 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.
- 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.
- 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.
- 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.
- 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
- 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.
- 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:
- sets various constants and parameters for the simulation.
- 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.
- contains the subroutine prototypes and structure definitions for I/O related subroutines.
- 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.
