MOM.F90

Go to the documentation of this file.

!> This is the main routine for MOM
module mom

! This file is part of MOM6. See LICENSE.md for the license.

use mom_variables, only : vertvisc_type
use mom_open_boundary, only : ocean_obc_type

! A Structure with pointers to forcing fields to drive MOM;
! all fluxes are positive downward.
use mom_forcing_type, only : forcing

use mom_variables, only: accel_diag_ptrs, cont_diag_ptrs, ocean_internal_state

! A structure containing pointers to various fields
! to describe surface state, and will be returned
! to the calling program.
use mom_variables, only : surface

! A structure containing pointers to an assortment of
! thermodynamic fields, including potential/Conservative
! temperature, salinity and mixed layer density.
use mom_variables, only: thermo_var_ptrs

! Infrastructure modules
use mom_debugging,            only : mom_debugging_init, hchksum, uvchksum
use mom_checksum_packages,    only : mom_thermo_chksum, mom_state_chksum, mom_accel_chksum
use mom_cpu_clock,            only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock,            only : clock_component, clock_subcomponent
use mom_cpu_clock,            only : clock_module_driver, clock_module, clock_routine
use mom_coms,                 only : reproducing_sum
use mom_coord_initialization, only : mom_initialize_coord
use mom_diag_mediator,        only : diag_mediator_init, enable_averaging
use mom_diag_mediator,        only : diag_mediator_infrastructure_init
use mom_diag_mediator,        only : diag_register_area_ids
use mom_diag_mediator,        only : diag_set_state_ptrs, diag_update_remap_grids
use mom_diag_mediator,        only : disable_averaging, post_data, safe_alloc_ptr
use mom_diag_mediator,        only : register_diag_field, register_static_field
use mom_diag_mediator,        only : register_scalar_field
use mom_diag_mediator,        only : set_axes_info, diag_ctrl, diag_masks_set
use mom_domains,              only : mom_domains_init, clone_mom_domain
use mom_domains,              only : sum_across_pes, pass_var, pass_vector
use mom_domains,              only : to_north, to_east, to_south, to_west
use mom_domains,              only : to_all, omit_corners, cgrid_ne, scalar_pair
use mom_domains,              only : create_group_pass, do_group_pass, group_pass_type
use mom_domains,              only : start_group_pass, complete_group_pass, omit_corners
use mom_error_handler,        only : mom_error, fatal, warning, is_root_pe
use mom_error_handler,        only : mom_set_verbosity, calltree_showquery
use mom_error_handler,        only : calltree_enter, calltree_leave, calltree_waypoint
use mom_file_parser,          only : read_param, get_param, log_version, param_file_type
use mom_fixed_initialization, only : mom_initialize_fixed
use mom_forcing_type,         only : mom_forcing_chksum
use mom_get_input,            only : get_mom_input, directories
use mom_io,                   only : mom_io_init, vardesc, var_desc
use mom_io,                   only : slasher, file_exists, read_data
use mom_obsolete_params,      only : find_obsolete_params
use mom_restart,              only : register_restart_field, query_initialized, save_restart
use mom_restart,              only : restart_init, mom_restart_cs
use mom_spatial_means,        only : global_area_mean, global_area_integral
use mom_state_initialization, only : mom_initialize_state
use mom_time_manager,         only : time_type, set_time, time_type_to_real, operator(+)
use mom_time_manager,         only : operator(-), operator(>), operator(*), operator(/)
use mom_time_manager,         only : increment_date
use mom_unit_tests,           only : unit_tests

! MOM core modules
use mom_ale,                   only : ale_init, ale_end, ale_main, ale_cs, adjustgridforintegrity
use mom_ale,                   only : ale_getcoordinate, ale_getcoordinateunits, ale_writecoordinatefile
use mom_ale,                   only : ale_updateverticalgridtype, ale_remap_init_conds, ale_register_diags
use mom_boundary_update,       only : call_obc_register, obc_register_end, update_obc_cs
use mom_continuity,            only : continuity, continuity_init, continuity_cs, continuity_stencil
use mom_coriolisadv,           only : coradcalc, coriolisadv_init, coriolisadv_cs
use mom_diabatic_driver,       only : diabatic, diabatic_driver_init, diabatic_cs
use mom_diabatic_driver,       only : adiabatic, adiabatic_driver_init, diabatic_driver_end
use mom_diagnostics,           only : calculate_diagnostic_fields, mom_diagnostics_init
use mom_diagnostics,           only : diagnostics_cs
use mom_diag_to_z,             only : calculate_z_diag_fields, calculate_z_transport
use mom_diag_to_z,             only : mom_diag_to_z_init, register_z_tracer, diag_to_z_cs
use mom_diag_to_z,             only : mom_diag_to_z_end
use mom_dynamics_unsplit,      only : step_mom_dyn_unsplit, register_restarts_dyn_unsplit
use mom_dynamics_unsplit,      only : initialize_dyn_unsplit, end_dyn_unsplit
use mom_dynamics_unsplit,      only : mom_dyn_unsplit_cs
use mom_dynamics_split_rk2,    only : step_mom_dyn_split_rk2, register_restarts_dyn_split_rk2
use mom_dynamics_split_rk2,    only : initialize_dyn_split_rk2, end_dyn_split_rk2
use mom_dynamics_split_rk2,    only : mom_dyn_split_rk2_cs
use mom_dynamics_unsplit_rk2,  only : step_mom_dyn_unsplit_rk2, register_restarts_dyn_unsplit_rk2
use mom_dynamics_unsplit_rk2,  only : initialize_dyn_unsplit_rk2, end_dyn_unsplit_rk2
use mom_dynamics_unsplit_rk2,  only : mom_dyn_unsplit_rk2_cs
use mom_dynamics_legacy_split, only : step_mom_dyn_legacy_split, register_restarts_dyn_legacy_split
use mom_dynamics_legacy_split, only : initialize_dyn_legacy_split, end_dyn_legacy_split
use mom_dynamics_legacy_split, only : adjustments_dyn_legacy_split, mom_dyn_legacy_split_cs
use mom_dyn_horgrid,           only : dyn_horgrid_type, create_dyn_horgrid, destroy_dyn_horgrid
use mom_eos,                   only : eos_init
use mom_eos,                   only : gsw_sp_from_sr, gsw_pt_from_ct
use mom_eos,                   only : calculate_density
use mom_debugging,             only : check_redundant
use mom_grid,                  only : ocean_grid_type, set_first_direction
use mom_grid,                  only : mom_grid_init, mom_grid_end
use mom_hor_index,             only : hor_index_type, hor_index_init
use mom_hor_visc,              only : horizontal_viscosity, hor_visc_init
use mom_interface_heights,     only : find_eta
use mom_lateral_mixing_coeffs, only : calc_slope_functions, varmix_init
use mom_lateral_mixing_coeffs, only : calc_resoln_function, varmix_cs
use mom_meke,                  only : meke_init, meke_alloc_register_restart, step_forward_meke, meke_cs
use mom_meke_types,            only : meke_type
use mom_mixed_layer_restrat,   only : mixedlayer_restrat, mixedlayer_restrat_init, mixedlayer_restrat_cs
use mom_mixed_layer_restrat,   only : mixedlayer_restrat_register_restarts
use mom_neutral_diffusion,     only : neutral_diffusion_cs, neutral_diffusion_diag_init
use mom_obsolete_diagnostics,  only : register_obsolete_diagnostics
use mom_open_boundary,         only : obc_registry_type
use mom_pressureforce,         only : pressureforce, pressureforce_init, pressureforce_cs
use mom_set_visc,              only : set_viscous_bbl, set_viscous_ml, set_visc_init
use mom_set_visc,              only : set_visc_register_restarts, set_visc_cs
use mom_sponge,                only : init_sponge_diags, sponge_cs
use mom_ale_sponge,            only : init_ale_sponge_diags, ale_sponge_cs
use mom_thickness_diffuse,     only : thickness_diffuse, thickness_diffuse_init, thickness_diffuse_cs
use mom_tidal_forcing,         only : tidal_forcing_init, tidal_forcing_cs
use mom_tracer_advect,         only : advect_tracer, tracer_advect_init
use mom_tracer_advect,         only : tracer_advect_end, tracer_advect_cs
use mom_tracer_hor_diff,       only : tracer_hordiff, tracer_hor_diff_init
use mom_tracer_hor_diff,       only : tracer_hor_diff_end, tracer_hor_diff_cs
use mom_tracer_registry,       only : register_tracer, tracer_registry_init
use mom_tracer_registry,       only : add_tracer_diagnostics, tracer_registry_type
use mom_tracer_registry,       only : lock_tracer_registry, tracer_registry_end
use mom_tracer_flow_control,   only : call_tracer_register, tracer_flow_control_cs
use mom_tracer_flow_control,   only : tracer_flow_control_init, call_tracer_surface_state
use mom_tracer_flow_control,   only : tracer_flow_control_end
use mom_transcribe_grid,       only : copy_dyngrid_to_mom_grid, copy_mom_grid_to_dyngrid
use mom_vert_friction,         only : vertvisc, vertvisc_remnant
use mom_vert_friction,         only : vertvisc_limit_vel, vertvisc_init
use mom_verticalgrid,          only : verticalgrid_type, verticalgridinit, verticalgridend
use mom_verticalgrid,          only : get_thickness_units, get_flux_units, get_tr_flux_units
! Offline modules
use mom_offline_main,          only : offline_transport_cs, offline_transport_init, update_offline_fields
use mom_offline_main,          only : insert_offline_main, extract_offline_main, post_offline_convergence_diags
use mom_offline_main,          only : register_diags_offline_transport, offline_advection_ale
use mom_offline_main,          only : offline_redistribute_residual, offline_diabatic_ale
use mom_offline_main,          only : offline_fw_fluxes_into_ocean, offline_fw_fluxes_out_ocean
use mom_offline_main,          only : offline_advection_layer, offline_transport_end
use mom_ale,                   only : ale_offline_tracer_final, ale_main_offline


implicit none ; private

#include <MOM_memory.h>

!> Control structure for this module
type, public :: mom_control_struct
  real allocable_, dimension(NIMEM_,NJMEM_,NKMEM_) :: &
    h, &      !< layer thickness (m or kg/m2 (H))
    t, &      !< potential temperature (degrees C)
    s         !< salinity (ppt)
  real allocable_, dimension(NIMEMB_PTR_,NJMEM_,NKMEM_) :: &
    u,  &     !< zonal velocity component (m/s)
    uh, &     !< uh = u * h * dy at u grid points (m3/s or kg/s)
    uhtr      !< accumulated zonal thickness fluxes to advect tracers (m3 or kg)
  real allocable_, dimension(NIMEM_,NJMEMB_PTR_,NKMEM_) :: &
    v,  &     !< meridional velocity (m/s)
    vh, &     !< vh = v * h * dx at v grid points (m3/s or kg/s)
    vhtr      !< accumulated meridional thickness fluxes to advect tracers (m3 or kg)
  real allocable_, dimension(NIMEM_,NJMEM_) :: &
    ave_ssh   !< time-averaged (ave over baroclinic time steps) sea surface height (meter)
  real, pointer, dimension(:,:)   :: hml => null() !< active mixed layer depth, in m
  real, pointer, dimension(:,:,:) :: &
    u_prev => null(), &  !< previous value of u stored for diagnostics
    v_prev => null()     !< previous value of v stored for diagnostics

  type(ocean_grid_type) :: g       !< structure containing metrics and grid info
  type(verticalgrid_type), pointer :: gv => null() !< structure containing vertical grid info
  type(thermo_var_ptrs) :: tv      !< structure containing pointers to available
                                   !! thermodynamic fields
  type(diag_ctrl)       :: diag    !< structure to regulate diagnostic output timing
  type(vertvisc_type)   :: visc    !< structure containing vertical viscosities,
                                   !! bottom drag viscosities, and related fields
  type(meke_type), pointer :: meke => null()  !<  structure containing fields
                                   !! related to the Mesoscale Eddy Kinetic Energy
  type(accel_diag_ptrs) :: adp     !< structure containing pointers to accelerations,
                                   !! for derived diagnostics (e.g., energy budgets)
  type(cont_diag_ptrs)  :: cdp     !< structure containing pointers continuity equation
                                   !! terms, for derived diagnostics (e.g., energy budgets)

  logical :: split                   !< If true, use the split time stepping scheme.
  logical :: legacy_split            !< If true, use the legacy split time stepping
                                     !! code with all the options that were available in
                                     !! the predecessor isopycnal model "GOLD".
  logical :: use_rk2                 !< If true, use RK2 instead of RK3 in unsplit mode
                                     !! (i.e., no split between barotropic and baroclinic).
  logical :: adiabatic               !< If true, then no diapycnal mass fluxes, with no calls
                                     !! to routines to calculate or apply diapycnal fluxes.
  logical :: use_temperature         !< If true, temp and saln used as state variables.
  logical :: calc_rho_for_sea_lev    !< If true, calculate rho to convert pressure to sea level
  logical :: use_frazil              !< If true, liquid seawater freezes if temp below freezing,
                                     !! with accumulated heat deficit returned to surface ocean.
  logical :: bound_salinity          !< If true, salt is added to keep salinity above
                                     !! a minimum value, and the deficit is reported.
  logical :: bulkmixedlayer          !< If true, a refined bulk mixed layer scheme is used
                                     !! with nkml sublayers and nkbl buffer layer.
  logical :: diabatic_first          !< If true, apply diabatic and thermodynamic
                                     !! processes before time stepping the dynamics.
  logical :: use_cont_abss           !< If true, , the prognostics T&S are the conservative temperature
                                     !! and absolute salinity. Care should be taken to convert them
                                     !! to potential temperature and practical salinity before
                                     !! exchanging them with the coupler and/or reporting T&S diagnostics.
  logical :: thickness_diffuse       !< If true, diffuse interface height w/ a diffusivity KHTH.
  logical :: thickness_diffuse_first !< If true, diffuse thickness before dynamics.
  logical :: mixedlayer_restrat      !< If true, use submesoscale mixed layer restratifying scheme.
  logical :: usemeke                 !< If true, call the MEKE parameterization.
  logical :: debug                   !< If true, write verbose checksums for debugging purposes.
  logical :: debug_truncations       !< If true, turn on diagnostics useful for debugging truncations.
  logical :: use_ale_algorithm       !< If true, use the ALE algorithm rather than layered
                                     !! isopycnal/stacked shallow water mode. This logical is
                                     !! set by calling the function useRegridding() from the
                                     !! MOM_regridding module.
  logical :: do_dynamics             !< If false, does not call step_MOM_dyn_*. This is an
                                     !! undocumented run-time flag that is fragile.
  logical :: offline_tracer_mode = .false.
                                     !< If true, step_offline() is called instead of step_MOM().
                                     !! This is intended for running MOM6 in offline tracer mode
  logical :: advect_ts               !< If false, then no horizontal advection of temperature
                                     !! and salnity is performed
  real    :: dt                      !< (baroclinic) dynamics time step (seconds)
  real    :: dt_therm                !< thermodynamics time step (seconds)
  logical :: thermo_spans_coupling   !< If true, thermodynamic and tracer time
                                     !! steps can span multiple coupled time steps.
  real    :: t_dyn_rel_adv           !< The time of the dynamics relative to tracer
                                     !! advection and lateral mixing (in seconds), or
                                     !! equivalently the elapsed time since advectively
                                     !! updating the tracers.  t_dyn_rel_adv is invariably
                                     !! positive and may span multiple coupling timesteps.
  real    :: t_dyn_rel_thermo        !< The time of the dynamics relative to diabatic
                                     !! processes and remapping (in seconds).  t_dyn_rel_thermo
                                     !! can be negative or positive depending on whether
                                     !! the diabatic processes are applied before or after
                                     !! the dynamics and may span multiple coupling timesteps.
  real    :: t_dyn_rel_diag          !< The time of the diagnostics relative to diabatic
                                     !! processes and remapping (in seconds).  t_dyn_rel_diag
                                     !! is always positive, since the diagnostics must lag.
  type(time_type) :: z_diag_interval !< amount of time between calculating Z-space diagnostics
  type(time_type) :: z_diag_time     !< next time to compute Z-space diagnostics
  type(time_type), pointer :: time   !< pointer to ocean clock
  real :: rel_time = 0.0             !< relative time (sec) since start of current execution
  real :: dtbt_reset_period          !< The time interval in seconds between dynamic
                                     !! recalculation of the barotropic time step.  If
                                     !! this is negative, it is never calculated, and
                                     !! if it is 0, it is calculated every step.

  logical :: interp_p_surf           !< If true, linearly interpolate surface pressure
                                     !! over the coupling time step, using specified value
                                     !! at the end of the coupling step. False by default.
  logical :: p_surf_prev_set         !< If true, p_surf_prev has been properly set from
                                     !! a previous time-step or the ocean restart file.
                                     !! This is only valid when interp_p_surf is true.

  real :: hmix                       !< Diagnostic mixed layer thickness (meter) when
                                     !! bulk mixed layer is not used.
  real :: hmix_uv                    !< Depth scale over which to average surface flow to
                                     !! feedback to the coupler/driver (m).
                                     !! bulk mixed layer is not used.
  real :: missing=-1.0e34            !< missing data value for masked fields

  ! Flags needed to reach between start and finish phases of initialization
  logical :: write_ic                !< If true, then the initial conditions will be written to file
  character(len=120) :: ic_file      !< A file into which the initial conditions are
                                     !! written in a new run if SAVE_INITIAL_CONDS is true.

  integer :: ntrunc                  !< number u,v truncations since last call to write_energy
  logical :: check_bad_surface_vals  !< If true, scan surface state for ridiculous values.
  real    :: bad_val_ssh_max         !< Maximum SSH before triggering bad value message
  real    :: bad_val_sst_max         !< Maximum SST before triggering bad value message
  real    :: bad_val_sst_min         !< Minimum SST before triggering bad value message
  real    :: bad_val_sss_max         !< Maximum SSS before triggering bad value message
  real    :: bad_val_column_thickness!< Minimum column thickness before triggering bad value message

  real, pointer, dimension(:,:) :: &
    p_surf_prev  => null(), & !< surface pressure (Pa) at end  previous call to step_MOM
    p_surf_begin => null(), & !< surface pressure (Pa) at start of step_MOM_dyn_...
    p_surf_end   => null()    !< surface pressure (Pa) at end   of step_MOM_dyn_...

  type(vardesc) :: &
    vd_t, &   !< vardesc array describing potential temperature
    vd_s      !< vardesc array describing salinity

  real, pointer, dimension(:,:,:) :: &  !< diagnostic arrays of advective/diffusive tracer fluxes
    t_adx => null(), t_ady => null(), t_diffx => null(), t_diffy => null(), &
    s_adx => null(), s_ady => null(), s_diffx => null(), s_diffy => null()

  real, pointer, dimension(:,:) :: & !< diagnostic arrays of vertically integrated advective/diffusive fluxes
    t_adx_2d => null(), t_ady_2d => null(), t_diffx_2d => null(), t_diffy_2d => null(), &
    s_adx_2d => null(), s_ady_2d => null(), s_diffx_2d => null(), s_diffy_2d => null()

  real, pointer, dimension(:,:,:) :: &  !< diagnostic arrays for advection tendencies and total tendencies
    t_advection_xy => null(), s_advection_xy => null(),  &
    t_prev         => null(), s_prev         => null(),  &
    th_prev        => null(), sh_prev        => null()

  real, pointer, dimension(:,:,:) :: & !< diagnostic arrays for variance decay through ALE
    t_squared => null(), s_squared => null()

  logical :: tendency_diagnostics = .false.

  ! diagnostic ids

  ! 3-d state fields
  integer :: id_u  = -1
  integer :: id_v  = -1
  integer :: id_h  = -1
  integer :: id_t  = -1
  integer :: id_s  = -1
  integer :: id_tcon  = -1
  integer :: id_sabs  = -1

  ! 2-d surface and bottom fields
  integer :: id_zos      = -1
  integer :: id_zossq    = -1
  integer :: id_volo     = -1
  integer :: id_ssh      = -1
  integer :: id_ssh_ga   = -1
  integer :: id_sst      = -1
  integer :: id_sst_sq   = -1
  integer :: id_sss      = -1
  integer :: id_sss_sq   = -1
  integer :: id_ssu      = -1
  integer :: id_ssv      = -1
  integer :: id_speed    = -1
  integer :: id_ssh_inst = -1
  integer :: id_tob      = -1
  integer :: id_sob      = -1
  integer :: id_sstcon   = -1
  integer :: id_sssabs   = -1

  ! heat and salt flux fields
  integer :: id_fraz         = -1
  integer :: id_salt_deficit = -1
  integer :: id_heat_pme     = -1
  integer :: id_intern_heat  = -1

  ! transport of temperature and salinity
  integer :: id_tadx      = -1
  integer :: id_tady      = -1
  integer :: id_tdiffx    = -1
  integer :: id_tdiffy    = -1
  integer :: id_sadx      = -1
  integer :: id_sady      = -1
  integer :: id_sdiffx    = -1
  integer :: id_sdiffy    = -1
  integer :: id_tadx_2d   = -1
  integer :: id_tady_2d   = -1
  integer :: id_tdiffx_2d = -1
  integer :: id_tdiffy_2d = -1
  integer :: id_sadx_2d   = -1
  integer :: id_sady_2d   = -1
  integer :: id_sdiffx_2d = -1
  integer :: id_sdiffy_2d = -1

  ! tendencies for temp/heat and saln/salt
  integer :: id_t_advection_xy    = -1
  integer :: id_t_advection_xy_2d = -1
  integer :: id_t_tendency        = -1
  integer :: id_th_tendency       = -1
  integer :: id_th_tendency_2d    = -1
  integer :: id_s_advection_xy    = -1
  integer :: id_s_advection_xy_2d = -1
  integer :: id_s_tendency        = -1
  integer :: id_sh_tendency       = -1
  integer :: id_sh_tendency_2d    = -1

  ! variance decay for temp and heat
  integer :: id_t_vardec = -1
  integer :: id_s_vardec = -1

  ! fields prior to doing dynamics
  integer :: id_h_pre_dyn = -1

  ! diagnostic for fields prior to applying diapycnal physics
  integer :: id_u_predia = -1
  integer :: id_v_predia = -1
  integer :: id_h_predia = -1
  integer :: id_t_predia = -1
  integer :: id_s_predia = -1
  integer :: id_e_predia = -1
  integer :: id_u_preale = -1
  integer :: id_v_preale = -1
  integer :: id_h_preale = -1
  integer :: id_t_preale = -1
  integer :: id_s_preale = -1
  integer :: id_e_preale = -1

  ! Diagnostics for tracer horizontal transport
  integer :: id_uhtr = -1, id_umo = -1, id_umo_2d = 1
  integer :: id_vhtr = -1, id_vmo = -1, id_vmo_2d = 1

  ! The remainder provides pointers to child module control structures.
  type(mom_dyn_unsplit_cs),      pointer :: dyn_unsplit_csp      => null()
  type(mom_dyn_unsplit_rk2_cs),  pointer :: dyn_unsplit_rk2_csp  => null()
  type(mom_dyn_split_rk2_cs),    pointer :: dyn_split_rk2_csp    => null()
  type(mom_dyn_legacy_split_cs), pointer :: dyn_legacy_split_csp => null()

  type(set_visc_cs),             pointer :: set_visc_csp           => null()
  type(diabatic_cs),             pointer :: diabatic_csp           => null()
  type(thickness_diffuse_cs),    pointer :: thickness_diffuse_csp  => null()
  type(mixedlayer_restrat_cs),   pointer :: mixedlayer_restrat_csp => null()
  type(meke_cs),                 pointer :: meke_csp               => null()
  type(varmix_cs),               pointer :: varmix                 => null()
  type(tracer_registry_type),    pointer :: tracer_reg             => null()
  type(tracer_advect_cs),        pointer :: tracer_adv_csp         => null()
  type(tracer_hor_diff_cs),      pointer :: tracer_diff_csp        => null()
  type(neutral_diffusion_cs),    pointer :: neutral_diffusion_csp  => null()
  type(tracer_flow_control_cs),  pointer :: tracer_flow_csp        => null()
  type(diagnostics_cs),          pointer :: diagnostics_csp        => null()
  type(diag_to_z_cs),            pointer :: diag_to_z_csp          => null()
  type(mom_restart_cs),          pointer :: restart_csp            => null()
  type(update_obc_cs),           pointer :: update_obc_csp         => null()
  type(ocean_obc_type),          pointer :: obc                    => null()
  type(sponge_cs),               pointer :: sponge_csp             => null()
  type(ale_sponge_cs),           pointer :: ale_sponge_csp         => null()
  type(ale_cs),                  pointer :: ale_csp                => null()
  type(offline_transport_cs),    pointer :: offline_csp            => null()

  ! These are used for group halo updates.
  type(group_pass_type) :: pass_tau_ustar_psurf
  type(group_pass_type) :: pass_ray
  type(group_pass_type) :: pass_bbl_thick_kv_bbl
  type(group_pass_type) :: pass_t_s_h
  type(group_pass_type) :: pass_t_s
  type(group_pass_type) :: pass_kv_turb
  type(group_pass_type) :: pass_uv_t_s_h
  type(group_pass_type) :: pass_ssh

end type mom_control_struct

public initialize_mom
public finish_mom_initialization
public step_mom
public step_offline
public mom_end
public calculate_surface_state

integer :: id_clock_ocean
integer :: id_clock_dynamics
integer :: id_clock_thermo
integer :: id_clock_tracer
integer :: id_clock_diabatic
integer :: id_clock_continuity  ! also in dynamics s/r
integer :: id_clock_thick_diff
integer :: id_clock_bbl_visc
integer :: id_clock_ml_restrat
integer :: id_clock_diagnostics
integer :: id_clock_z_diag
integer :: id_clock_init
integer :: id_clock_mom_init
integer :: id_clock_pass       ! also in dynamics d/r
integer :: id_clock_pass_init  ! also in dynamics d/r
integer :: id_clock_ale
integer :: id_clock_other
integer :: id_clock_offline_tracer

contains


!> This subroutine orchestrates the time stepping of MOM.  The adiabatic
!! dynamics are stepped by calls to one of the step_MOM_dyn_...routines.
!! The action of lateral processes on tracers occur in calls to
!! advect_tracer and tracer_hordiff.  Vertical mixing and possibly remapping
!! occur inside of diabatic.
subroutine step_mom(fluxes, state, Time_start, time_interval, CS)
  type(forcing),    intent(inout)    :: fluxes        !< pointers to forcing fields
  type(surface),    intent(inout)    :: state         !< surface ocean state
  type(time_type),  intent(in)       :: Time_start    !< starting time of a segment, as a time type
  real,             intent(in)       :: time_interval !< time interval covered by this run segment, in s.
  type(mom_control_struct), pointer  :: CS            !< control structure from initialize_MOM

  ! local
  type(ocean_grid_type), pointer :: G ! pointer to a structure containing
                                      ! metrics and related information
  type(verticalgrid_type),  pointer :: GV => null()
  integer, save :: nt_debug = 1 ! running number of iterations, for debugging only.
  integer       :: ntstep ! time steps between tracer updates or diabatic forcing
  integer       :: n_max  ! number of steps to take in this call

  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, n
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB

  real :: dt              ! baroclinic time step (sec)
  real :: dtth            ! time step for thickness diffusion (sec)
  real :: dtdia           ! time step for diabatic processes (sec)
  real :: dt_therm        ! a limited and quantized version of CS%dt_therm (sec)
  real :: dtbt_reset_time ! value of CS%rel_time when DTBT was last calculated (sec)

  real :: mass_src_time   ! The amount of time for the surface mass source from
                          ! precipitation-evaporation, rivers, etc., that should
                          ! be applied to the start of the barotropic solver to
                          ! avoid generating tsunamis, in s.  This is negative
                          ! if the precipation has already been applied to the
                          ! layers, and positive if it will be applied later.

  real :: wt_end, wt_beg
  real :: bbl_time_int    ! The amount of time over which the calculated BBL
                          ! properties will apply, for use in diagnostics.

  logical :: calc_dtbt                 ! Indicates whether the dynamically adjusted
                                       ! barotropic time step needs to be updated.
  logical :: do_advection              ! If true, it is time to advect tracers.
  logical :: do_calc_bbl               ! If true, calculate the boundary layer properties.
  logical :: thermo_does_span_coupling ! If true, thermodynamic forcing spans
                                       ! multiple dynamic timesteps.
  real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: &
    eta_av, &   ! average sea surface height or column mass over a timestep (meter or kg/m2)
    ssh         ! sea surface height based on eta_av (meter or kg/m2)

  real, pointer, dimension(:,:,:) :: &
    u, & ! u : zonal velocity component (m/s)
    v, & ! v : meridional velocity component (m/s)
    h    ! h : layer thickness (meter (Bouss) or kg/m2 (non-Bouss))

  ! Store the layer thicknesses before any changes by the dynamics. This is necessary for remapped
  ! mass transport diagnostics
  real, dimension(SZI_(CS%G),SZJ_(CS%G),SZK_(CS%G)) :: h_pre_dyn
  real :: tot_wt_ssh, Itot_wt_ssh

  type(time_type) :: Time_local
  logical :: showCallTree
  ! These are used for group halo passes.
  logical :: do_pass_kv_turb, do_pass_Ray, do_pass_kv_bbl_thick

  g => cs%G ; gv => cs%GV
  is   = g%isc  ; ie   = g%iec  ; js   = g%jsc  ; je   = g%jec ; nz = g%ke
  isq  = g%IscB ; ieq  = g%IecB ; jsq  = g%JscB ; jeq  = g%JecB
  isd  = g%isd  ; ied  = g%ied  ; jsd  = g%jsd  ; jed  = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB
  u => cs%u ; v => cs%v ; h => cs%h

  call cpu_clock_begin(id_clock_ocean)
  call cpu_clock_begin(id_clock_other)

  if (cs%debug) then
    call mom_state_chksum("Beginning of step_MOM ", u, v, h, cs%uh, cs%vh, g, gv)
    call hchksum(cs%h,"CS%h beginning of step_MOM",g%HI, scale=gv%H_to_m)
  endif

  showcalltree = calltree_showquery()
  if (showcalltree) call calltree_enter("step_MOM(), MOM.F90")

  ! First determine the time step that is consistent with this call.
  ! It is anticipated that the time step will almost always coincide
  ! with dt. In addition, ntstep is determined, subject to the constraint
  ! that ntstep cannot exceed n_max.
  if (time_interval <= cs%dt) then
    n_max = 1
  else
    n_max = ceiling(time_interval/cs%dt - 0.001)
  endif

  dt    = time_interval / real(n_max)
  dtdia = 0.0
  thermo_does_span_coupling = (cs%thermo_spans_coupling .and. &
                              (cs%dt_therm > 1.5*time_interval))
  if (thermo_does_span_coupling) then
    ! Set dt_therm to be an integer multiple of the coupling time step.
    dt_therm = time_interval * floor(cs%dt_therm / time_interval + 0.001)
    ntstep = floor(dt_therm/dt + 0.001)
  else
    ntstep = max(1,min(n_max,floor(cs%dt_therm/dt + 0.001)))
    dt_therm = dt*ntstep
  endif

  if (.not.ASSOCIATED(fluxes%p_surf)) cs%interp_p_surf = .false.

  !---------- Begin setup for group halo pass

  call cpu_clock_begin(id_clock_pass)
  call create_group_pass(cs%pass_tau_ustar_psurf, fluxes%taux, fluxes%tauy, g%Domain)
  if (ASSOCIATED(fluxes%ustar)) &
    call create_group_pass(cs%pass_tau_ustar_psurf, fluxes%ustar(:,:), g%Domain)
  if (ASSOCIATED(fluxes%p_surf)) &
    call create_group_pass(cs%pass_tau_ustar_psurf, fluxes%p_surf(:,:), g%Domain)

  do_pass_ray = .false.
  if ((.not.g%Domain%symmetric) .and. &
      associated(cs%visc%Ray_u) .and. associated(cs%visc%Ray_v)) then
    call create_group_pass(cs%pass_ray, cs%visc%Ray_u, cs%visc%Ray_v, g%Domain, &
                           to_north+to_east+scalar_pair+omit_corners, cgrid_ne, halo=1)
    do_pass_ray = .true.
  endif
  do_pass_kv_bbl_thick = .false.
  if (associated(cs%visc%bbl_thick_u) .and. associated(cs%visc%bbl_thick_v)) then
    call create_group_pass(cs%pass_bbl_thick_kv_bbl, cs%visc%bbl_thick_u, &
                           cs%visc%bbl_thick_v, g%Domain, &
                           to_north+to_east+scalar_pair+omit_corners, cgrid_ne, halo=1)
    do_pass_kv_bbl_thick = .true.
  endif
  if (associated(cs%visc%kv_bbl_u) .and. associated(cs%visc%kv_bbl_v)) then
    call create_group_pass(cs%pass_bbl_thick_kv_bbl, cs%visc%kv_bbl_u, &
                           cs%visc%kv_bbl_v, g%Domain, &
                           to_north+to_east+scalar_pair+omit_corners, cgrid_ne, halo=1)
    do_pass_kv_bbl_thick = .true.
  endif
  do_pass_kv_turb = associated(cs%visc%Kv_turb)
  if (associated(cs%visc%Kv_turb)) &
    call create_group_pass(cs%pass_kv_turb, cs%visc%Kv_turb, g%Domain, to_all+omit_corners, halo=1)

  if (.not.cs%adiabatic .AND. cs%use_ALE_algorithm ) then
    if (cs%use_temperature) then
      call create_group_pass(cs%pass_T_S_h, cs%tv%T, g%Domain, to_all+omit_corners, halo=1)
      call create_group_pass(cs%pass_T_S_h, cs%tv%S, g%Domain, to_all+omit_corners, halo=1)
    endif
    call create_group_pass(cs%pass_T_S_h, h, g%Domain, to_all+omit_corners, halo=1)
  endif

  if ((cs%adiabatic .OR. cs%diabatic_first) .AND. cs%use_temperature) then
    call create_group_pass(cs%pass_T_S, cs%tv%T, g%Domain, halo=1) ! Could also omit corners?
    call create_group_pass(cs%pass_T_S, cs%tv%S, g%Domain, halo=1)
  endif

  !---------- End setup for group halo pass


  if (g%nonblocking_updates) then
    call start_group_pass(cs%pass_tau_ustar_psurf, g%Domain)
  else
    call do_group_pass(cs%pass_tau_ustar_psurf, g%Domain)
  endif
  call cpu_clock_end(id_clock_pass)

  if (ASSOCIATED(cs%tv%frazil))        cs%tv%frazil(:,:)        = 0.0
  if (ASSOCIATED(cs%tv%salt_deficit))  cs%tv%salt_deficit(:,:)  = 0.0
  if (ASSOCIATED(cs%tv%TempxPmE))      cs%tv%TempxPmE(:,:)      = 0.0
  if (ASSOCIATED(cs%tv%internal_heat)) cs%tv%internal_heat(:,:) = 0.0

  cs%rel_time = 0.0

  tot_wt_ssh = 0.0
  do j=js,je ; do i=is,ie ; cs%ave_ssh(i,j) = 0.0 ; ssh(i,j) = cs%missing; enddo ; enddo

  if (associated(cs%VarMix)) then
    call enable_averaging(time_interval, time_start+set_time(int(time_interval)), &
                          cs%diag)
    call calc_resoln_function(h, cs%tv, g, gv, cs%VarMix)
    call disable_averaging(cs%diag)
  endif

  if (g%nonblocking_updates) then
    call cpu_clock_begin(id_clock_pass)
    call complete_group_pass(cs%pass_tau_ustar_psurf, g%Domain)
    call cpu_clock_end(id_clock_pass)
  endif

  if (cs%interp_p_surf) then
    if (.not.ASSOCIATED(cs%p_surf_end))   allocate(cs%p_surf_end(isd:ied,jsd:jed))
    if (.not.ASSOCIATED(cs%p_surf_begin)) allocate(cs%p_surf_begin(isd:ied,jsd:jed))
    if (.not.cs%p_surf_prev_set) then
      do j=jsd,jed ; do i=isd,ied
        cs%p_surf_prev(i,j) = fluxes%p_surf(i,j)
      enddo ; enddo
      cs%p_surf_prev_set = .true.
    endif
  else
    cs%p_surf_end  => fluxes%p_surf
  endif

  if (cs%debug) then
    call mom_state_chksum("Before steps ", u, v, h, cs%uh, cs%vh, g, gv)
    call mom_forcing_chksum("Before steps", fluxes, g, haloshift=0)
    call check_redundant("Before steps ", u, v, g)
    call check_redundant("Before steps ", fluxes%taux, fluxes%tauy, g)
  endif
  call cpu_clock_end(id_clock_other)

  do n=1,n_max

    nt_debug = nt_debug + 1

    ! Set the universally visible time to the middle of the time step
    cs%Time = time_start + set_time(int(floor(cs%rel_time+0.5*dt+0.5)))
    cs%rel_time = cs%rel_time + dt

    ! Set the local time to the end of the time step.
    time_local = time_start + set_time(int(floor(cs%rel_time+0.5)))
    if (showcalltree) call calltree_enter("DT cycles (step_MOM) n=",n)

    !===========================================================================
    ! This is the first place where the diabatic processes and remapping could occur.
    if (cs%diabatic_first .and. (cs%t_dyn_rel_adv==0.0)) then ! do thermodynamics.

      if (thermo_does_span_coupling) then
        dtdia = dt_therm
        if ((fluxes%dt_buoy_accum > 0.0) .and. (dtdia > time_interval) .and. &
            (abs(fluxes%dt_buoy_accum - dtdia) > 1e-6*dtdia)) then
          call mom_error(fatal, "step_MOM: Mismatch between long thermodynamic "//&
            "timestep and time over which buoyancy fluxes have been accumulated.")
        endif
        call mom_error(fatal, "MOM is not yet set up to have restarts that work "//&
          "with THERMO_SPANS_COUPLING and DIABATIC_FIRST.")
      else
        dtdia = dt*min(ntstep,n_max-(n-1))
      endif

      ! The end-time of the diagnostic interval needs to be set ahead if there
      ! are multiple dynamic time steps worth of thermodynamics applied here.
      call enable_averaging(dtdia, time_local + &
                                   set_time(int(floor(dtdia-dt+0.5))), cs%diag)

      if (cs%debug) then
        call uvchksum("Pre set_viscous_BBL [uv]", u, v, g%HI, haloshift=1)
        call hchksum(h,"Pre set_viscous_BBL h", g%HI, haloshift=1, scale=gv%H_to_m)
        if (associated(cs%tv%T)) call hchksum(cs%tv%T, "Pre set_viscous_BBL T", g%HI, haloshift=1)
        if (associated(cs%tv%S)) call hchksum(cs%tv%S, "Pre set_viscous_BBL S", g%HI, haloshift=1)
      endif

      !   Calculate the BBL properties and store them inside visc (u,h).
      ! This is here so that CS%visc is updated before diabatic() when
      ! DIABATIC_FIRST=True. Otherwise diabatic() is called after the dynamics
      ! and set_viscous_BBL is called as a part of the dynamic stepping.
      call cpu_clock_begin(id_clock_bbl_visc)
      call set_viscous_bbl(u, v, h, cs%tv, cs%visc, g, gv, cs%set_visc_CSp)
      call cpu_clock_end(id_clock_bbl_visc)

      call cpu_clock_begin(id_clock_pass)
      if (do_pass_ray) call do_group_pass(cs%pass_ray, g%Domain )
      if (do_pass_kv_bbl_thick) call do_group_pass(cs%pass_bbl_thick_kv_bbl, g%Domain)
      call cpu_clock_end(id_clock_pass)
      if (showcalltree) call calltree_waypoint("done with set_viscous_BBL (diabatic_first)")

      call cpu_clock_begin(id_clock_thermo)

      ! These adjustments are a part of an archaic time stepping algorithm.
      if (cs%split .and. cs%legacy_split .and. .not.cs%adiabatic) then
        call adjustments_dyn_legacy_split(u, v, h, dt, g, gv, cs%dyn_legacy_split_CSp)
      endif

      ! Apply diabatic forcing, do mixing, and regrid.
      call step_mom_thermo(cs, g, gv, u, v, h, cs%tv, fluxes, dtdia)

      ! The diabatic processes are now ahead of the dynamics by dtdia.
      cs%t_dyn_rel_thermo = -dtdia
      if (showcalltree) call calltree_waypoint("finished diabatic_first (step_MOM)")

      call disable_averaging(cs%diag)
      call cpu_clock_end(id_clock_thermo)

    endif ! end of block "(CS%diabatic_first .and. (CS%t_dyn_rel_adv==0.0))"

    !===========================================================================
    ! This is the start of the dynamics stepping part of the algorithm.

    call cpu_clock_begin(id_clock_dynamics)
    call disable_averaging(cs%diag)

    if (cs%thickness_diffuse .and. cs%thickness_diffuse_first) then
      if (cs%t_dyn_rel_adv == 0.0) then
        if (thermo_does_span_coupling) then
          dtth = dt_therm
        else
          dtth = dt*min(ntstep,n_max-n+1)
        endif

        call enable_averaging(dtth,time_local+set_time(int(floor(dtth-dt+0.5))), cs%diag)
        call cpu_clock_begin(id_clock_thick_diff)
        if (associated(cs%VarMix)) &
          call calc_slope_functions(h, cs%tv, dt, g, gv, cs%VarMix)
        call thickness_diffuse(h, cs%uhtr, cs%vhtr, cs%tv, dtth, g, gv, &
                               cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp)
        call cpu_clock_end(id_clock_thick_diff)
        call cpu_clock_begin(id_clock_pass)
        call pass_var(h, g%Domain) !###, halo=max(2,cont_stensil))
        call cpu_clock_end(id_clock_pass)
        call disable_averaging(cs%diag)
        if (showcalltree) call calltree_waypoint("finished thickness_diffuse_first (step_MOM)")

        ! Whenever thickness changes let the diag manager know, target grids
        ! for vertical remapping may need to be regenerated.
        call diag_update_remap_grids(cs%diag)

      endif
    endif

    ! The bottom boundary layer properties are out-of-date and need to be
    ! recalculated.  This always occurs at the start of a coupling time
    ! step because the externally prescribed stresses may have changed.
    do_calc_bbl = ((cs%t_dyn_rel_adv == 0.0) .or. (n==1))
    if (do_calc_bbl) then
      ! Calculate the BBL properties and store them inside visc (u,h).
      call cpu_clock_begin(id_clock_bbl_visc)
      bbl_time_int = max(dt, min(dt_therm - cs%t_dyn_rel_adv, dt*(1+n_max-n)) )
      call enable_averaging(bbl_time_int, &
                time_local+set_time(int(bbl_time_int-dt+0.5)), cs%diag)
      call set_viscous_bbl(u, v, h, cs%tv, cs%visc, g, gv, cs%set_visc_CSp)
      call disable_averaging(cs%diag)
      call cpu_clock_end(id_clock_bbl_visc)
      if (showcalltree) call calltree_waypoint("done with set_viscous_BBL (step_MOM)")
    endif

    call cpu_clock_begin(id_clock_pass)
    if (do_pass_kv_turb) call do_group_pass(cs%pass_kv_turb, g%Domain)
    call cpu_clock_end(id_clock_pass)

    if (do_calc_bbl) then
      call cpu_clock_begin(id_clock_pass)
      if (g%nonblocking_updates) then
        if (do_pass_ray) call start_group_pass(cs%pass_Ray, g%Domain)
        if (do_pass_kv_bbl_thick) call start_group_pass(cs%pass_bbl_thick_kv_bbl, g%Domain)
        ! do_calc_bbl will be set to .false. when the message passing is complete.
      else
        if (do_pass_ray) call do_group_pass(cs%pass_Ray, g%Domain)
        if (do_pass_kv_bbl_thick) call do_group_pass(cs%pass_bbl_thick_kv_bbl, g%Domain)
      endif
      call cpu_clock_end(id_clock_pass)
    endif

    if (cs%interp_p_surf) then
      wt_end = real(n) / real(n_max)
      wt_beg = real(n-1) / real(n_max)
      do j=jsd,jed ; do i=isd,ied
        cs%p_surf_end(i,j) = wt_end * fluxes%p_surf(i,j) + &
                        (1.0-wt_end) * cs%p_surf_prev(i,j)
        cs%p_surf_begin(i,j) = wt_beg * fluxes%p_surf(i,j) + &
                        (1.0-wt_beg) * cs%p_surf_prev(i,j)
      enddo ; enddo
    endif

    if (associated(cs%u_prev) .and. associated(cs%v_prev)) then
      do k=1,nz ; do j=jsd,jed ; do i=isdb,iedb
        cs%u_prev(i,j,k) = u(i,j,k)
      enddo ; enddo ; enddo
      do k=1,nz ; do j=jsdb,jedb ; do i=isd,ied
        cs%v_prev(i,j,k) = u(i,j,k)
      enddo ; enddo ; enddo
    endif

    ! Store pre-dynamics layer thicknesses so that mass fluxes are remapped correctly
    do k=1,nz ; do j=jsd,jed ; do i=isd,ied
      h_pre_dyn(i,j,k) = h(i,j,k)
    enddo ; enddo ; enddo

    if (g%nonblocking_updates) then ; if (do_calc_bbl) then
      call cpu_clock_begin(id_clock_pass)
        if (do_pass_ray) call complete_group_pass(cs%pass_Ray, g%Domain)
        if (do_pass_kv_bbl_thick) call complete_group_pass(cs%pass_bbl_thick_kv_bbl, g%Domain)
      call cpu_clock_end(id_clock_pass)
    endif ; endif

    if (cs%do_dynamics .and. cs%split) then !--------------------------- start SPLIT
      ! This section uses a split time stepping scheme for the dynamic equations,
      ! basically the stacked shallow water equations with viscosity.

      calc_dtbt = .false.
      if ((cs%dtbt_reset_period >= 0.0) .and. &
          ((n==1) .or. (cs%dtbt_reset_period == 0.0) .or. &
           (cs%rel_time >= dtbt_reset_time + 0.999*cs%dtbt_reset_period))) then
        calc_dtbt = .true.
        dtbt_reset_time = cs%rel_time
      endif

      mass_src_time = cs%t_dyn_rel_thermo
      if (cs%legacy_split) then
        call step_mom_dyn_legacy_split(u, v, h, cs%tv, cs%visc, &
                    time_local, dt, fluxes, cs%p_surf_begin, cs%p_surf_end, &
                    mass_src_time, dt_therm, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
                    eta_av, g, gv, cs%dyn_legacy_split_CSp, calc_dtbt, cs%VarMix, cs%MEKE)
      else
        call step_mom_dyn_split_rk2(u, v, h, cs%tv, cs%visc, &
                    time_local, dt, fluxes, cs%p_surf_begin, cs%p_surf_end, &
                    mass_src_time, dt_therm, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
                    eta_av, g, gv, cs%dyn_split_RK2_CSp, calc_dtbt, cs%VarMix, cs%MEKE)
      endif
      if (showcalltree) call calltree_waypoint("finished step_MOM_dyn_split (step_MOM)")


    elseif (cs%do_dynamics) then ! --------------------------------------------------- not SPLIT
      !   This section uses an unsplit stepping scheme for the dynamic
      ! equations; basically the stacked shallow water equations with viscosity.
      ! Because the time step is limited by CFL restrictions on the external
      ! gravity waves, the unsplit is usually much less efficient that the split
      ! approaches. But because of its simplicity, the unsplit method is very
      ! useful for debugging purposes.

      if (cs%use_RK2) then
        call step_mom_dyn_unsplit_rk2(u, v, h, cs%tv, cs%visc, time_local, dt, fluxes, &
                 cs%p_surf_begin, cs%p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
                 eta_av, g, gv, cs%dyn_unsplit_RK2_CSp, cs%VarMix, cs%MEKE)
      else
        call step_mom_dyn_unsplit(u, v, h, cs%tv, cs%visc, time_local, dt, fluxes, &
                 cs%p_surf_begin, cs%p_surf_end, cs%uh, cs%vh, cs%uhtr, cs%vhtr, &
                 eta_av, g, gv, cs%dyn_unsplit_CSp, cs%VarMix, cs%MEKE)
      endif
      if (showcalltree) call calltree_waypoint("finished step_MOM_dyn_unsplit (step_MOM)")

    endif ! -------------------------------------------------- end SPLIT



    if (cs%thickness_diffuse .and. .not.cs%thickness_diffuse_first) then
      call cpu_clock_begin(id_clock_thick_diff)

      if (cs%debug) call hchksum(h,"Pre-thickness_diffuse h", g%HI, haloshift=0, scale=gv%H_to_m)

      if (associated(cs%VarMix)) &
        call calc_slope_functions(h, cs%tv, dt, g, gv, cs%VarMix)
      call thickness_diffuse(h, cs%uhtr, cs%vhtr, cs%tv, dt, g, gv, &
                             cs%MEKE, cs%VarMix, cs%CDp, cs%thickness_diffuse_CSp)

      if (cs%debug) call hchksum(h,"Post-thickness_diffuse h", g%HI, haloshift=1, scale=gv%H_to_m)
      call cpu_clock_end(id_clock_thick_diff)
      call cpu_clock_begin(id_clock_pass)
      call pass_var(h, g%Domain) !###, halo=max(2,cont_stensil))
      call cpu_clock_end(id_clock_pass)
      if (showcalltree) call calltree_waypoint("finished thickness_diffuse (step_MOM)")
    endif

    ! apply the submesoscale mixed layer restratification parameterization
    if (cs%mixedlayer_restrat) then
      if (cs%debug) then
        call hchksum(h,"Pre-mixedlayer_restrat h", g%HI, haloshift=1, scale=gv%H_to_m)
        call uvchksum("Pre-mixedlayer_restrat uhtr", &
                      cs%uhtr, cs%vhtr, g%HI, haloshift=0)
      endif
      call cpu_clock_begin(id_clock_ml_restrat)
      call mixedlayer_restrat(h, cs%uhtr ,cs%vhtr, cs%tv, fluxes, dt, cs%visc%MLD, &
                              cs%VarMix, g, gv, cs%mixedlayer_restrat_CSp)
      call cpu_clock_end(id_clock_ml_restrat)
      call cpu_clock_begin(id_clock_pass)
      call pass_var(h, g%Domain) !###, halo=max(2,cont_stensil))
      call cpu_clock_end(id_clock_pass)
      if (cs%debug) then
        call hchksum(h,"Post-mixedlayer_restrat h", g%HI, haloshift=1, scale=gv%H_to_m)
        call uvchksum("Post-mixedlayer_restrat [uv]htr", &
                      cs%uhtr, cs%vhtr, g%HI, haloshift=0)
      endif
    endif

    ! Whenever thickness changes let the diag manager know, target grids
    ! for vertical remapping may need to be regenerated.
    call diag_update_remap_grids(cs%diag)

    if (cs%useMEKE) call step_forward_meke(cs%MEKE, h, cs%VarMix%SN_u, cs%VarMix%SN_v, &
                                           cs%visc, dt, g, gv, cs%MEKE_CSp, cs%uhtr, cs%vhtr)
    call disable_averaging(cs%diag)

    ! Advance the dynamics time by dt.
    cs%t_dyn_rel_adv = cs%t_dyn_rel_adv + dt
    cs%t_dyn_rel_thermo = cs%t_dyn_rel_thermo + dt
    cs%t_dyn_rel_diag = cs%t_dyn_rel_diag + dt

    call cpu_clock_end(id_clock_dynamics)

    !===========================================================================
    ! This is the start of the tracer advection part of the algorithm.

    if (thermo_does_span_coupling) then
      do_advection = (cs%t_dyn_rel_adv + 0.5*dt > dt_therm)
    else
      do_advection = ((mod(n,ntstep) == 0) .or. (n==n_max))
    endif

    if (do_advection) then ! Do advective transport and lateral tracer mixing.

      if (cs%debug) then
        call cpu_clock_begin(id_clock_other)
        call uvchksum("Pre-advection [uv]", u, v, g%HI, haloshift=2)
        call hchksum(h,"Pre-advection h", g%HI, haloshift=1, scale=gv%H_to_m)
        call uvchksum("Pre-advection uhtr", cs%uhtr, cs%vhtr, g%HI, &
                      haloshift=0, scale=gv%H_to_m)
      ! call MOM_state_chksum("Pre-advection ", u, v, &
      !                       h, CS%uhtr, CS%vhtr, G, GV, haloshift=1)
          if (associated(cs%tv%T)) call hchksum(cs%tv%T, "Pre-advection T", g%HI, haloshift=1)
          if (associated(cs%tv%S)) call hchksum(cs%tv%S, "Pre-advection S", g%HI, haloshift=1)
          if (associated(cs%tv%frazil)) call hchksum(cs%tv%frazil, &
                         "Pre-advection frazil", g%HI, haloshift=0)
          if (associated(cs%tv%salt_deficit)) call hchksum(cs%tv%salt_deficit, &
                         "Pre-advection salt deficit", g%HI, haloshift=0)
      ! call MOM_thermo_chksum("Pre-advection ", CS%tv, G)
        call check_redundant("Pre-advection ", u, v, g)
        call cpu_clock_end(id_clock_other)
      endif

      call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer)
      call enable_averaging(cs%t_dyn_rel_adv, time_local, cs%diag)

      call advect_tracer(h, cs%uhtr, cs%vhtr, cs%OBC, cs%t_dyn_rel_adv, g, gv, &
                         cs%tracer_adv_CSp, cs%tracer_Reg)
      call tracer_hordiff(h, cs%t_dyn_rel_adv, cs%MEKE, cs%VarMix, g, gv, &
                          cs%tracer_diff_CSp, cs%tracer_Reg, cs%tv)
      if (showcalltree) call calltree_waypoint("finished tracer advection/diffusion (step_MOM)")
      call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo)

      call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)
      call post_transport_diagnostics(g, gv, cs, cs%diag, cs%t_dyn_rel_adv, h, h_pre_dyn)
      ! Rebuild the remap grids now that we've posted the fields which rely on thicknesses
      ! from before the dynamics calls
      call diag_update_remap_grids(cs%diag)

      call disable_averaging(cs%diag)
      call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)

      ! Reset the accumulated transports to 0 and record that the dynamics
      ! and advective times now agree.
      call cpu_clock_begin(id_clock_thermo) ; call cpu_clock_begin(id_clock_tracer)
      cs%uhtr(:,:,:) = 0.0
      cs%vhtr(:,:,:) = 0.0
      cs%t_dyn_rel_adv = 0.0
      call cpu_clock_end(id_clock_tracer) ; call cpu_clock_end(id_clock_thermo)

    endif

    !===========================================================================
    ! This is the second place where the diabatic processes and remapping could occur.
    if (cs%t_dyn_rel_adv == 0.0) then
      call cpu_clock_begin(id_clock_thermo)

      if (.not.cs%diabatic_first) then
        dtdia = cs%t_dyn_rel_thermo
        if (thermo_does_span_coupling .and. (abs(dt_therm - dtdia) > 1e-6*dt_therm)) then
          call mom_error(fatal, "step_MOM: Mismatch between dt_therm and dtdia "//&
                         "before call to diabatic.")
        endif

        call enable_averaging(cs%t_dyn_rel_thermo, time_local, cs%diag)

        ! These adjustments are a part of an archaic time stepping algorithm.
        if (cs%split .and. cs%legacy_split .and. .not.cs%adiabatic) then
          call adjustments_dyn_legacy_split(u, v, h, dt, g, gv, cs%dyn_legacy_split_CSp)
        endif

        ! Apply diabatic forcing, do mixing, and regrid.
        call step_mom_thermo(cs, g, gv, u, v, h, cs%tv, fluxes, dtdia)

        call disable_averaging(cs%diag)

      else  !  "else branch for if (.not.CS%diabatic_first) then"
        if (abs(cs%t_dyn_rel_thermo) > 1e-6*dt) call mom_error(fatal, &
              "step_MOM: Mismatch between the dynamics and diabatic times "//&
              "with DIABATIC_FIRST.")

        ! Tracers have been advected and diffused, and need halo updates.
        if (cs%use_temperature) then
          call cpu_clock_begin(id_clock_pass)
          call do_group_pass(cs%pass_T_S, g%Domain)
          call cpu_clock_end(id_clock_pass)
        endif
      endif ! close of "if (.not.CS%diabatic_first) then ; if (.not.CS%adiabatic)"

      ! Record that the dynamics and diabatic processes are synchronized.
      cs%t_dyn_rel_thermo = 0.0
      call cpu_clock_end(id_clock_thermo)
    endif

    call cpu_clock_begin(id_clock_dynamics)

    ! Determining the time-average sea surface height is part of the algorithm.
    ! This may be eta_av if Boussinesq, or need to be diagnosed if not.
    tot_wt_ssh = tot_wt_ssh + dt
    call find_eta(h, cs%tv, gv%g_Earth, g, gv, ssh, eta_av)
    do j=js,je ; do i=is,ie
      cs%ave_ssh(i,j) = cs%ave_ssh(i,j) + dt*ssh(i,j)
    enddo ; enddo
    call cpu_clock_end(id_clock_dynamics)

    !===========================================================================
    ! Calculate diagnostics at the end of the time step.
    call cpu_clock_begin(id_clock_other) ; call cpu_clock_begin(id_clock_diagnostics)

    call enable_averaging(dt,time_local, cs%diag)
    ! These diagnostics are available every time step.
    if (cs%id_u > 0) call post_data(cs%id_u, u, cs%diag)
    if (cs%id_v > 0) call post_data(cs%id_v, v, cs%diag)
    if (cs%id_h > 0) call post_data(cs%id_h, h, cs%diag)
    if (cs%id_ssh_inst > 0) call post_data(cs%id_ssh_inst, ssh, cs%diag)
    call disable_averaging(cs%diag)

    if (cs%t_dyn_rel_adv == 0.0) then
      ! Diagnostics that require the complete state to be up-to-date can be calculated.

      call enable_averaging(cs%t_dyn_rel_diag, time_local, cs%diag)
      call calculate_diagnostic_fields(u, v, h, cs%uh, cs%vh, cs%tv, cs%ADp, &
                          cs%CDp, fluxes, cs%t_dyn_rel_diag, g, gv, cs%diagnostics_CSp)
      call post_ts_diagnostics(cs, g, gv, cs%tv, cs%diag, cs%t_dyn_rel_diag)
      if (showcalltree) call calltree_waypoint("finished calculate_diagnostic_fields (step_MOM)")
      call disable_averaging(cs%diag)
      cs%t_dyn_rel_diag = 0.0

      call cpu_clock_begin(id_clock_z_diag)
      if (time_local + set_time(int(0.5*dt_therm)) > cs%Z_diag_time) then
        call enable_averaging(real(time_type_to_real(CS%Z_diag_interval)), &
                              CS%Z_diag_time, CS%diag)
        call calculate_z_diag_fields(u, v, h, ssh, fluxes%frac_shelf_h, &
                                     g, gv, cs%diag_to_Z_CSp)
        cs%Z_diag_time = cs%Z_diag_time + cs%Z_diag_interval
        call disable_averaging(cs%diag)
        if (showcalltree) call calltree_waypoint("finished calculate_Z_diag_fields (step_MOM)")
      endif
      call cpu_clock_end(id_clock_z_diag)
    endif
    call cpu_clock_end(id_clock_diagnostics) ; call cpu_clock_end(id_clock_other)

    if (showcalltree) call calltree_leave("DT cycles (step_MOM)")

  enddo ! complete the n loop

  call cpu_clock_begin(id_clock_other)

  itot_wt_ssh = 1.0/tot_wt_ssh
  do j=js,je ; do i=is,ie
    cs%ave_ssh(i,j) = cs%ave_ssh(i,j)*itot_wt_ssh
  enddo ; enddo

  call enable_averaging(dt*n_max, time_local, cs%diag)
  call post_integrated_diagnostics(cs, g, gv, cs%diag, dt*n_max, cs%tv, fluxes)
  call disable_averaging(cs%diag)

  if (showcalltree) call calltree_waypoint("calling calculate_surface_state (step_MOM)")
  call adjust_ssh_for_p_atm(cs, g, gv, cs%ave_ssh, fluxes%p_surf_SSH)
  call calculate_surface_state(state, u, v, h, cs%ave_ssh, g, gv, cs)

  call enable_averaging(dt*n_max, time_local, cs%diag)
  call post_surface_diagnostics(cs, g, cs%diag, state)
  call disable_averaging(cs%diag)

  if (cs%interp_p_surf) then ; do j=jsd,jed ; do i=isd,ied
    cs%p_surf_prev(i,j) = fluxes%p_surf(i,j)
  enddo ; enddo ; endif

  call cpu_clock_end(id_clock_other)

  if (showcalltree) call calltree_leave("step_MOM()")
  call cpu_clock_end(id_clock_ocean)

end subroutine step_mom

!> MOM_step_thermo orchestrates the thermodynamic time stepping and vertical
!! remapping, via calls to diabatic (or adiabatic) and ALE_main.
subroutine step_mom_thermo(CS, G, GV, u, v, h, tv, fluxes, dtdia)
  type(mom_control_struct), intent(inout) :: CS     !< control structure
  type(ocean_grid_type),    intent(inout) :: G      !< ocean grid structure
  type(verticalgrid_type),  intent(inout) :: GV     !< ocean vertical grid structure
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                            intent(inout) :: u      !< zonal velocity (m/s)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                            intent(inout) :: v      !< meridional velocity (m/s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                            intent(inout) :: h      !< layer thickness (m or kg/m2)
  type(thermo_var_ptrs),    intent(inout) :: tv     !< A structure pointing to various thermodynamic variables
  type(forcing),            intent(inout) :: fluxes !< pointers to forcing fields
  real,                     intent(in)    :: dtdia  !< The time interval over which to advance, in s

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: eta_predia, eta_preale
  integer :: i, j, k, is, ie, js, je, nz! , Isq, Ieq, Jsq, Jeq, n
  logical :: use_ice_shelf ! Needed for selecting the right ALE interface.
  logical :: showCallTree

  is   = g%isc  ; ie   = g%iec  ; js   = g%jsc  ; je   = g%jec ; nz = g%ke
  showcalltree = calltree_showquery()
  if (showcalltree) call calltree_enter("step_MOM_thermo(), MOM.F90")

  use_ice_shelf = .false.
  if (associated(fluxes%frac_shelf_h)) use_ice_shelf = .true.

  if (.not.cs%adiabatic) then

    if (cs%debug) then
      call uvchksum("Pre-diabatic [uv]", u, v, g%HI, haloshift=2)
      call hchksum(h,"Pre-diabatic h", g%HI, haloshift=1, scale=gv%H_to_m)
      call uvchksum("Pre-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
                    haloshift=0, scale=gv%H_to_m)
    ! call MOM_state_chksum("Pre-diabatic ",u, v, h, CS%uhtr, CS%vhtr, G, GV)
      call mom_thermo_chksum("Pre-diabatic ", cs%tv, g,haloshift=0)
      call check_redundant("Pre-diabatic ", u, v, g)
      call mom_forcing_chksum("Pre-diabatic", fluxes, g, haloshift=0)
    endif

    if (cs%id_u_predia > 0) call post_data(cs%id_u_predia, u, cs%diag)
    if (cs%id_v_predia > 0) call post_data(cs%id_v_predia, v, cs%diag)
    if (cs%id_h_predia > 0) call post_data(cs%id_h_predia, h, cs%diag)
    if (cs%id_T_predia > 0) call post_data(cs%id_T_predia, cs%tv%T, cs%diag)
    if (cs%id_S_predia > 0) call post_data(cs%id_S_predia, cs%tv%S, cs%diag)
    if (cs%id_e_predia > 0) then
      call find_eta(h, cs%tv, gv%g_Earth, g, gv, eta_predia)
      call post_data(cs%id_e_predia, eta_predia, cs%diag)
    endif

    call cpu_clock_begin(id_clock_diabatic)
    call diabatic(u, v, h, tv, cs%Hml, fluxes, cs%visc, cs%ADp, cs%CDp, &
                  dtdia, g, gv, cs%diabatic_CSp)
    fluxes%fluxes_used = .true.
    call cpu_clock_end(id_clock_diabatic)

    if (cs%id_u_preale > 0) call post_data(cs%id_u_preale, u,    cs%diag)
    if (cs%id_v_preale > 0) call post_data(cs%id_v_preale, v,    cs%diag)
    if (cs%id_h_preale > 0) call post_data(cs%id_h_preale, h,    cs%diag)
    if (cs%id_T_preale > 0) call post_data(cs%id_T_preale, tv%T, cs%diag)
    if (cs%id_S_preale > 0) call post_data(cs%id_S_preale, tv%S, cs%diag)
    if (cs%id_e_preale > 0) then
      call find_eta(h, tv, gv%g_Earth, g, gv, eta_preale)
      call post_data(cs%id_e_preale, eta_preale, cs%diag)
    endif

    if (showcalltree) call calltree_waypoint("finished diabatic (step_MOM_thermo)")

    ! Regridding/remapping is done here, at end of thermodynamics time step
    ! (that may comprise several dynamical time steps)
    ! The routine 'ALE_main' can be found in 'MOM_ALE.F90'.
    if ( cs%use_ALE_algorithm ) then
!         call pass_vector(u, v, G%Domain)
      call do_group_pass(cs%pass_T_S_h, g%Domain)

      ! update squared quantities
      if (associated(cs%S_squared)) then ; do k=1,nz ; do j=js,je ; do i=is,ie
        cs%S_squared(i,j,k) = tv%S(i,j,k)**2
      enddo ; enddo ; enddo ; endif
      if (associated(cs%T_squared)) then ; do k=1,nz ; do j=js,je ; do i=is,ie
        cs%T_squared(i,j,k) = tv%T(i,j,k)**2
      enddo ; enddo ; enddo ; endif

      if (cs%debug) then
        call mom_state_chksum("Pre-ALE ", u, v, h, cs%uh, cs%vh, g, gv)
        call hchksum(tv%T,"Pre-ALE T", g%HI, haloshift=1)
        call hchksum(tv%S,"Pre-ALE S", g%HI, haloshift=1)
        call check_redundant("Pre-ALE ", u, v, g)
      endif
      call cpu_clock_begin(id_clock_ale)
      if (use_ice_shelf) then
        call ale_main(g, gv, h, u, v, tv, cs%tracer_Reg, cs%ALE_CSp, dtdia, &
                      fluxes%frac_shelf_h)
      else
        call ale_main(g, gv, h, u, v, tv, cs%tracer_Reg, cs%ALE_CSp, dtdia)
      endif

      if (showcalltree) call calltree_waypoint("finished ALE_main (step_MOM_thermo)")
      call cpu_clock_end(id_clock_ale)
    endif   ! endif for the block "if ( CS%use_ALE_algorithm )"

    call cpu_clock_begin(id_clock_pass)
    call do_group_pass(cs%pass_uv_T_S_h, g%Domain)
    call cpu_clock_end(id_clock_pass)

    if (cs%debug .and. cs%use_ALE_algorithm) then
      call mom_state_chksum("Post-ALE ", u, v, h, cs%uh, cs%vh, g, gv)
      call hchksum(tv%T, "Post-ALE T", g%HI, haloshift=1)
      call hchksum(tv%S, "Post-ALE S", g%HI, haloshift=1)
      call check_redundant("Post-ALE ", u, v, g)
    endif

    ! Whenever thickness changes let the diag manager know, target grids
    ! for vertical remapping may need to be regenerated. This needs to
    ! happen after the H update and before the next post_data.
    call diag_update_remap_grids(cs%diag)

    call post_diags_ts_vardec(g, cs, dtdia)

    if (cs%debug) then
      call uvchksum("Post-diabatic u", u, v, g%HI, haloshift=2)
      call hchksum(h, "Post-diabatic h", g%HI, haloshift=1, scale=gv%H_to_m)
      call uvchksum("Post-diabatic [uv]h", cs%uhtr, cs%vhtr, g%HI, &
                    haloshift=0, scale=gv%H_to_m)
    ! call MOM_state_chksum("Post-diabatic ", u, v, &
    !                       h, CS%uhtr, CS%vhtr, G, GV, haloshift=1)
      if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", g%HI, haloshift=1)
      if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", g%HI, haloshift=1)
      if (associated(tv%frazil)) call hchksum(tv%frazil, &
                               "Post-diabatic frazil", g%HI, haloshift=0)
      if (associated(tv%salt_deficit)) call hchksum(tv%salt_deficit, &
                               "Post-diabatic salt deficit", g%HI, haloshift=0)
    ! call MOM_thermo_chksum("Post-diabatic ", tv, G)
      call check_redundant("Post-diabatic ", u, v, g)
    endif

  else   ! complement of "if (.not.CS%adiabatic)"

    call cpu_clock_begin(id_clock_diabatic)
    call adiabatic(h, tv, fluxes, dtdia, g, gv, cs%diabatic_CSp)
    fluxes%fluxes_used = .true.
    call cpu_clock_end(id_clock_diabatic)

    if (cs%use_temperature) then
      call cpu_clock_begin(id_clock_pass)
      call do_group_pass(cs%pass_T_S, g%Domain)
      call cpu_clock_end(id_clock_pass)
      if (cs%debug) then
        if (associated(tv%T)) call hchksum(tv%T, "Post-diabatic T", g%HI, haloshift=1)
        if (associated(tv%S)) call hchksum(tv%S, "Post-diabatic S", g%HI, haloshift=1)
      endif
    endif

  endif   ! endif for the block "if (.not.CS%adiabatic)"

  if (showcalltree) call calltree_leave("step_MOM_thermo(), MOM.F90")

end subroutine step_mom_thermo


!> step_offline is the main driver for running tracers offline in MOM6. This has been primarily
!! developed with ALE configurations in mind. Some work has been done in isopycnal configuration, but
!! the work is very preliminary. Some more detail about this capability along with some of the subroutines
!! called here can be found in tracers/MOM_offline_control.F90
subroutine step_offline(fluxes, state, Time_start, time_interval, CS)
  type(forcing),    intent(inout)    :: fluxes        !< pointers to forcing fields
  type(surface),    intent(inout)    :: state         !< surface ocean state
  type(time_type),  intent(in)       :: Time_start    !< starting time of a segment, as a time type
  real,             intent(in)       :: time_interval !< time interval
  type(mom_control_struct), pointer  :: CS            !< control structure from initialize_MOM

  ! Local pointers
  type(ocean_grid_type),      pointer :: G  => null() ! Pointer to a structure containing
                                                      ! metrics and related information
  type(verticalgrid_type),    pointer :: GV => null() ! Pointer to structure containing information
                                                      ! about the vertical grid

  logical :: first_iter    !< True if this is the first time step_offline has been called in a given interval
  logical :: last_iter     !< True if this is the last time step_tracer is to be called in an offline interval
  logical :: do_vertical   !< If enough time has elapsed, do the diabatic tracer sources/sinks
  logical :: adv_converged !< True if all the horizontal fluxes have been used

  integer :: dt_offline, dt_offline_vertical
  logical :: skip_diffusion
  integer :: id_eta_diff_end

  integer, pointer :: accumulated_time
  integer :: i,j,k
  integer :: is, ie, js, je, isd, ied, jsd, jed

  ! 3D pointers
  real, dimension(:,:,:), pointer   :: &
    uhtr, vhtr, &
    eatr, ebtr, &
    h_end

  ! 2D Array for diagnostics
  real, dimension(SZI_(CS%G),SZJ_(CS%G)) :: eta_pre, eta_end
  type(time_type) :: Time_end    ! End time of a segment, as a time type

  ! Grid-related pointer assignments
  g => cs%G
  gv => cs%GV

  is  = g%isc  ; ie  = g%iec  ; js  = g%jsc  ; je  = g%jec
  isd = g%isd  ; ied = g%ied  ; jsd = g%jsd  ; jed = g%jed

  call cpu_clock_begin(id_clock_offline_tracer)
  call extract_offline_main(cs%offline_CSp, uhtr, vhtr, eatr, ebtr, h_end, accumulated_time, &
                            dt_offline, dt_offline_vertical, skip_diffusion)
  time_end = increment_date(time_start, seconds=floor(time_interval+0.001))
  call enable_averaging(time_interval, time_end, cs%diag)

  ! Check to see if this is the first iteration of the offline interval
  if(accumulated_time==0) then
    first_iter = .true.
  else ! This is probably unnecessary but is used to guard against unwanted behavior
    first_iter = .false.
  endif

  ! Check to see if vertical tracer functions should  be done
  if ( mod(accumulated_time, dt_offline_vertical) == 0 ) then
    do_vertical = .true.
  else
    do_vertical = .false.
  endif

  ! Increment the amount of time elapsed since last read and check if it's time to roll around
  accumulated_time = mod(accumulated_time + int(time_interval), dt_offline)
  if(accumulated_time==0) then
    last_iter = .true.
  else
    last_iter = .false.
  endif

  if(cs%use_ALE_algorithm) then
    ! If this is the first iteration in the offline timestep, then we need to read in fields and
    ! perform the main advection.
    if (first_iter) then
      if(is_root_pe()) print *, "Reading in new offline fields"
      ! Read in new transport and other fields
      ! call update_transport_from_files(G, GV, CS%offline_CSp, h_end, eatr, ebtr, uhtr, vhtr, &
      !     CS%tv%T, CS%tv%S, fluxes, CS%use_ALE_algorithm)
      ! call update_transport_from_arrays(CS%offline_CSp)
      call update_offline_fields(cs%offline_CSp, cs%h, fluxes, cs%use_ALE_algorithm)

      ! Apply any fluxes into the ocean
      call offline_fw_fluxes_into_ocean(g, gv, cs%offline_CSp, fluxes, cs%h)

      if (.not.cs%diabatic_first) then
        call offline_advection_ale(fluxes, time_start, time_interval, cs%offline_CSp, id_clock_ale, &
            cs%h, uhtr, vhtr, converged=adv_converged)

        ! Redistribute any remaining transport
        call offline_redistribute_residual(cs%offline_CSp, cs%h, uhtr, vhtr, adv_converged)

        ! Perform offline diffusion if requested
        if (.not. skip_diffusion) then
          if (associated(cs%VarMix)) then
            call pass_var(cs%h,g%Domain)
            call calc_resoln_function(cs%h, cs%tv, g, gv, cs%VarMix)
            call calc_slope_functions(cs%h, cs%tv, REAL(dt_offline), G, GV, CS%VarMix)
          endif
          call tracer_hordiff(cs%h, REAL(dt_offline), CS%MEKE, CS%VarMix, G, GV, &
              CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv)
        endif
      endif
    endif
    ! The functions related to column physics of tracers is performed separately in ALE mode
    if (do_vertical) then
      call offline_diabatic_ale(fluxes, time_start, time_end, cs%offline_CSp, cs%h, eatr, ebtr)
    endif

    ! Last thing that needs to be done is the final ALE remapping
    if(last_iter) then
      if (cs%diabatic_first) then
        call offline_advection_ale(fluxes, time_start, time_interval, cs%offline_CSp, id_clock_ale, &
            cs%h, uhtr, vhtr, converged=adv_converged)

        ! Redistribute any remaining transport and perform the remaining advection
        call offline_redistribute_residual(cs%offline_CSp, cs%h, uhtr, vhtr, adv_converged)
                ! Perform offline diffusion if requested
        if (.not. skip_diffusion) then
          if (associated(cs%VarMix)) then
            call pass_var(cs%h,g%Domain)
            call calc_resoln_function(cs%h, cs%tv, g, gv, cs%VarMix)
            call calc_slope_functions(cs%h, cs%tv, REAL(dt_offline), G, GV, CS%VarMix)
          endif
          call tracer_hordiff(cs%h, REAL(dt_offline), CS%MEKE, CS%VarMix, G, GV, &
              CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv)
        endif
      endif

      if(is_root_pe()) print *, "Last iteration of offline interval"

      ! Apply freshwater fluxes out of the ocean
      call offline_fw_fluxes_out_ocean(g, gv, cs%offline_CSp, fluxes, cs%h)
      ! These diagnostic can be used to identify which grid points did not converge within
      ! the specified number of advection sub iterations
      call post_offline_convergence_diags(cs%offline_CSp, cs%h, h_end, uhtr, vhtr)

      ! Call ALE one last time to make sure that tracers are remapped onto the layer thicknesses
      ! stored from the forward run
      call cpu_clock_begin(id_clock_ale)
      call ale_offline_tracer_final( g, gv, cs%h, cs%tv, h_end, cs%tracer_Reg, cs%ALE_CSp)
      call cpu_clock_end(id_clock_ale)
      call pass_var(cs%h, g%Domain)
    endif
  else ! NON-ALE MODE...NOT WELL TESTED
    call mom_error(warning, &
        "Offline tracer mode in non-ALE configuration has not been thoroughly tested")
    ! Note that for the layer mode case, the calls to tracer sources and sinks is embedded in
    ! main_offline_advection_layer. Warning: this may not be appropriate for tracers that
    ! exchange with the atmosphere
    if(time_interval .NE. dt_offline) then
      call mom_error(fatal, &
          "For offline tracer mode in a non-ALE configuration, dt_offline must equal time_interval")
    endif
    call update_offline_fields(cs%offline_CSp, cs%h, fluxes, cs%use_ALE_algorithm)
    call offline_advection_layer(fluxes, time_start, time_interval, cs%offline_CSp, &
        cs%h, eatr, ebtr, uhtr, vhtr)
    ! Perform offline diffusion if requested
    if (.not. skip_diffusion) then
      call tracer_hordiff(h_end, REAL(dt_offline), CS%MEKE, CS%VarMix, G, GV, &
        CS%tracer_diff_CSp, CS%tracer_Reg, CS%tv)
    endif

    cs%h = h_end

    call pass_var(cs%tv%T, g%Domain)
    call pass_var(cs%tv%S, g%Domain)
    call pass_var(cs%h, g%Domain)

  endif

  call adjust_ssh_for_p_atm(cs, g, gv, cs%ave_ssh, fluxes%p_surf_SSH)
  call calculate_surface_state(state, cs%u, cs%v, cs%h, cs%ave_ssh, g, gv, cs)

  call disable_averaging(cs%diag)
  call pass_var(cs%tv%T,g%Domain)
  call pass_var(cs%tv%S,g%Domain)
  call pass_var(cs%h,g%Domain)

  fluxes%fluxes_used = .true.

  call cpu_clock_end(id_clock_offline_tracer)

end subroutine step_offline

!> This subroutine initializes MOM.
subroutine initialize_mom(Time, param_file, dirs, CS, Time_in, offline_tracer_mode)
  type(time_type), target,   intent(inout) :: Time        !< model time, set in this routine
  type(param_file_type),     intent(out)   :: param_file  !< structure indicating paramater file to parse
  type(directories),         intent(out)   :: dirs        !< structure with directory paths
  type(mom_control_struct),  pointer       :: CS          !< pointer set in this routine to MOM control structure
  type(time_type), optional, intent(in)    :: Time_in     !< time passed to MOM_initialize_state when
                                                          !! model is not being started from a restart file
  logical,         optional, intent(out)   :: offline_tracer_mode !< True if tracers are being run offline

  ! local
  type(ocean_grid_type),  pointer :: G => null() ! A pointer to a structure with metrics and related
  type(hor_index_type)            :: HI  !  A hor_index_type for array extents
  type(verticalgrid_type), pointer :: GV => null()
  type(dyn_horgrid_type), pointer :: dG => null()
  type(diag_ctrl),        pointer :: diag

  character(len=4), parameter :: vers_num = 'v2.0'

! This include declares and sets the variable "version".
#include "version_variable.h"

  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
  integer :: IsdB, IedB, JsdB, JedB
  real    :: dtbt
  real    :: Z_diag_int  ! minimum interval between calc depth-space diagnostics (sec)

  real, allocatable, dimension(:,:,:) :: e   ! interface heights (meter)
  real, allocatable, dimension(:,:)   :: eta ! free surface height (m) or bottom press (Pa)
  real, allocatable, dimension(:,:)   :: area_shelf_h ! area occupied by ice shelf
  real, dimension(:,:), allocatable, target  :: frac_shelf_h ! fraction of total area occupied by ice shelf
  real, dimension(:,:), pointer :: shelf_area
  type(mom_restart_cs),  pointer      :: restart_CSp_tmp => null()
  type(group_pass_type) :: tmp_pass_uv_T_S_h

  real    :: default_val       ! default value for a parameter
  logical :: write_geom_files  ! If true, write out the grid geometry files.
  logical :: new_sim
  logical :: use_geothermal    ! If true, apply geothermal heating.
  logical :: use_EOS           ! If true, density calculated from T & S using an equation of state.
  logical :: symmetric         ! If true, use symmetric memory allocation.
  logical :: save_IC           ! If true, save the initial conditions.
  logical :: do_unit_tests     ! If true, call unit tests.
  logical :: test_grid_copy = .false.
  logical :: use_ice_shelf     ! Needed for ALE
  logical :: global_indexing   ! If true use global horizontal index values instead
                               ! of having the data domain on each processor start at 1.
  logical :: bathy_at_vel      ! If true, also define bathymetric fields at the
                               ! the velocity points.
  integer :: first_direction   ! An integer that indicates which direction is to be
                               ! updated first in directionally split parts of the
                               ! calculation.  This can be altered during the course
                               ! of the run via calls to set_first_direction.
  integer :: nkml, nkbl, verbosity, write_geom
  integer :: dynamics_stencil  ! The computational stencil for the calculations
                               ! in the dynamic core.

  type(time_type)                 :: Start_time
  type(ocean_internal_state)      :: MOM_internal_state
  character(len=200) :: area_varname, ice_shelf_file, inputdir, filename

  if (associated(cs)) then
    call mom_error(warning, "initialize_MOM called with an associated "// &
                            "control structure.")
    return
  endif
  allocate(cs)

  if (test_grid_copy) then ; allocate(g)
  else ; g => cs%G ; endif

  cs%Time => time

  id_clock_init = cpu_clock_id('Ocean Initialization', grain=clock_subcomponent)
  call cpu_clock_begin(id_clock_init)

  start_time = time ; if (present(time_in)) start_time = time_in

  call get_mom_input(param_file, dirs)

  verbosity = 2 ; call read_param(param_file, "VERBOSITY", verbosity)
  call mom_set_verbosity(verbosity)
  call calltree_enter("initialize_MOM(), MOM.F90")

  call find_obsolete_params(param_file)

  ! Read relevant parameters and write them to the model log.
  call log_version(param_file, "MOM", version, "")
  call get_param(param_file, "MOM", "VERBOSITY", verbosity,  &
                 "Integer controlling level of messaging\n" // &
                 "\t0 = Only FATAL messages\n" // &
                 "\t2 = Only FATAL, WARNING, NOTE [default]\n" // &
                 "\t9 = All)", default=2)
  call get_param(param_file, "MOM", "DO_UNIT_TESTS", do_unit_tests, &
                 "If True, exercises unit tests at model start up.", &
                 default=.false.)
  if (do_unit_tests) then
    call unit_tests(verbosity)
  endif

  call get_param(param_file, "MOM", "SPLIT", cs%split, &
                 "Use the split time stepping if true.", default=.true.)
  if (cs%split) then
    call get_param(param_file, "MOM", "USE_LEGACY_SPLIT", cs%legacy_split, &
                 "If true, use the full range of options available from \n"//&
                 "the older GOLD-derived split time stepping code.", &
                 default=.false.)
    cs%use_RK2 = .false.
  else
    call get_param(param_file, "MOM", "USE_RK2", cs%use_RK2, &
                 "If true, use RK2 instead of RK3 in the unsplit time stepping.", &
                 default=.false.)
    cs%legacy_split = .false.
  endif

  call get_param(param_file, "MOM", "CALC_RHO_FOR_SEA_LEVEL", cs%calc_rho_for_sea_lev, &
                 "If true, the in-situ density is used to calculate the\n"//&
                 "effective sea level that is returned to the coupler. If false,\n"//&
                 "the Boussinesq parameter RHO_0 is used.", default=.false.)
  call get_param(param_file, "MOM", "ENABLE_THERMODYNAMICS", cs%use_temperature, &
                 "If true, Temperature and salinity are used as state \n"//&
                 "variables.", default=.true.)
  call get_param(param_file, "MOM", "USE_EOS", use_eos, &
                 "If true,  density is calculated from temperature and \n"//&
                 "salinity with an equation of state.  If USE_EOS is \n"//&
                 "true, ENABLE_THERMODYNAMICS must be true as well.", &
                 default=cs%use_temperature)
  call get_param(param_file, "MOM", "DIABATIC_FIRST", cs%diabatic_first, &
                 "If true, apply diabatic and thermodynamic processes, \n"//&
                 "including buoyancy forcing and mass gain or loss, \n"//&
                 "before stepping the dynamics forward.", default=.false.)
  call get_param(param_file, "MOM", "USE_CONTEMP_ABSSAL", cs%use_conT_absS, &
                 "If true, , the prognostics T&S are the conservative temperature \n"//&
                 "and absolute salinity. Care should be taken to convert them \n"//&
                 "to potential temperature and practical salinity before  \n"//&
                 "exchanging them with the coupler and/or reporting T&S diagnostics. \n"&
                 , default=.false.)
  call get_param(param_file, "MOM", "ADIABATIC", cs%adiabatic, &
                 "There are no diapycnal mass fluxes if ADIABATIC is \n"//&
                 "true. This assumes that KD = KDML = 0.0 and that \n"//&
                 "there is no buoyancy forcing, but makes the model \n"//&
                 "faster by eliminating subroutine calls.", default=.false.)
  call get_param(param_file, "MOM", "DO_DYNAMICS", cs%do_dynamics, &
                 "If False, skips the dynamics calls that update u & v, as well as\n"//&
                 "the gravity wave adjustment to h. This is a fragile feature and\n"//&
                 "thus undocumented.", default=.true., do_not_log=.true. )
  call get_param(param_file, "MOM", "ADVECT_TS", cs%advect_TS , &
                 "If True, advect temperature and salinity horizontally\n"//&
                 "If False, T/S are registered for advection.\n"//&
                 "This is intended only to be used in offline tracer mode.", &
                 "and is by default false in that case", &
                 do_not_log = .true., default=.true. )
  if (present(offline_tracer_mode)) then ! Only read this parameter in solo mode
    call get_param(param_file, "MOM", "OFFLINE_TRACER_MODE", cs%offline_tracer_mode, &
                 "If true, barotropic and baroclinic dynamics, thermodynamics\n"//&
                 "are all bypassed with all the fields necessary to integrate\n"//&
                 "the tracer advection and diffusion equation are read in from\n"//&
                 "files stored from a previous integration of the prognostic model.\n"//&
                 "NOTE: This option only used in the ocean_solo_driver.", default=.false.)
    if(cs%offline_tracer_mode) then
      call get_param(param_file, "MOM", "ADVECT_TS", cs%advect_TS , &
                   "If True, advect temperature and salinity horizontally\n"//&
                   "If False, T/S are registered for advection.\n"//&
                   "This is intended only to be used in offline tracer mode."//&
                   "and is by default false in that case", &
                   default=.false. )
    endif
  endif
  call get_param(param_file, "MOM", "USE_REGRIDDING", cs%use_ALE_algorithm , &
                 "If True, use the ALE algorithm (regridding/remapping).\n"//&
                 "If False, use the layered isopycnal algorithm.", default=.false. )
  call get_param(param_file, "MOM", "BULKMIXEDLAYER", cs%bulkmixedlayer, &
                 "If true, use a Kraus-Turner-like bulk mixed layer \n"//&
                 "with transitional buffer layers.  Layers 1 through  \n"//&
                 "NKML+NKBL have variable densities. There must be at \n"//&
                 "least NKML+NKBL+1 layers if BULKMIXEDLAYER is true. \n"//&
                 "BULKMIXEDLAYER can not be used with USE_REGRIDDING. \n"//&
                 "The default is influenced by ENABLE_THERMODYNAMICS.", &
                 default=cs%use_temperature .and. .not.cs%use_ALE_algorithm)
  call get_param(param_file, "MOM", "THICKNESSDIFFUSE", cs%thickness_diffuse, &
                 "If true, interface heights are diffused with a \n"//&
                 "coefficient of KHTH.", default=.false.)
  call get_param(param_file, "MOM",  "THICKNESSDIFFUSE_FIRST", &
                                      cs%thickness_diffuse_first, &
                 "If true, do thickness diffusion before dynamics.\n"//&
                 "This is only used if THICKNESSDIFFUSE is true.", &
                 default=.false.)
  call get_param(param_file, "MOM", "BATHYMETRY_AT_VEL", bathy_at_vel, &
                 "If true, there are separate values for the basin depths \n"//&
                 "at velocity points.  Otherwise the effects of topography \n"//&
                 "are entirely determined from thickness points.", &
                 default=.false.)

  call get_param(param_file, "MOM", "DEBUG", cs%debug, &
                 "If true, write out verbose debugging data.", default=.false.)
  call get_param(param_file, "MOM", "DEBUG_TRUNCATIONS", cs%debug_truncations, &
                 "If true, calculate all diagnostics that are useful for \n"//&
                 "debugging truncations.", default=.false.)

  call get_param(param_file, "MOM", "DT", cs%dt, &
                 "The (baroclinic) dynamics time step.  The time-step that \n"//&
                 "is actually used will be an integer fraction of the \n"//&
                 "forcing time-step (DT_FORCING in ocean-only mode or the \n"//&
                 "coupling timestep in coupled mode.)", units="s", &
                 fail_if_missing=.true.)
  call get_param(param_file, "MOM", "DT_THERM", cs%dt_therm, &
                 "The thermodynamic and tracer advection time step. \n"//&
                 "Ideally DT_THERM should be an integer multiple of DT \n"//&
                 "and less than the forcing or coupling time-step, unless \n"//&
                 "THERMO_SPANS_COUPLING is true, in which case DT_THERM \n"//&
                 "can be an integer multiple of the coupling timestep.  By \n"//&
                 "default DT_THERM is set to DT.", units="s", default=cs%dt)
  call get_param(param_file, "MOM", "THERMO_SPANS_COUPLING", cs%thermo_spans_coupling, &
                 "If true, the MOM will take thermodynamic and tracer \n"//&
                 "timesteps that can be longer than the coupling timestep. \n"//&
                 "The actual thermodynamic timestep that is used in this \n"//&
                 "case is the largest integer multiple of the coupling \n"//&
                 "timestep that is less than or equal to DT_THERM.", default=.false.)

  if (.not.cs%bulkmixedlayer) then
    call get_param(param_file, "MOM", "HMIX_SFC_PROP", cs%Hmix, &
                 "If BULKMIXEDLAYER is false, HMIX_SFC_PROP is the depth \n"//&
                 "over which to average to find surface properties like \n"//&
                 "SST and SSS or density (but not surface velocities).", &
                 units="m", default=1.0)
    call get_param(param_file, "MOM", "HMIX_UV_SFC_PROP", cs%Hmix_UV, &
                 "If BULKMIXEDLAYER is false, HMIX_UV_SFC_PROP is the depth\n"//&
                 "over which to average to find surface flow properties,\n"//&
                 "SSU, SSV. A non-positive value indicates no averaging.", &
                 units="m", default=0.)
  endif
  call get_param(param_file, "MOM", "MIN_Z_DIAG_INTERVAL", z_diag_int, &
                 "The minimum amount of time in seconds between \n"//&
                 "calculations of depth-space diagnostics. Making this \n"//&
                 "larger than DT_THERM reduces the  performance penalty \n"//&
                 "of regridding to depth online.", units="s", default=0.0)
  call get_param(param_file, "MOM", "INTERPOLATE_P_SURF", cs%interp_p_surf, &
                 "If true, linearly interpolate the surface pressure \n"//&
                 "over the coupling time step, using the specified value \n"//&
                 "at the end of the step.", default=.false.)

  if (cs%split) then
    call get_param(param_file, "MOM", "DTBT", dtbt, default=-0.98)
    default_val = cs%dt_therm ; if (dtbt > 0.0) default_val = -1.0
    cs%dtbt_reset_period = -1.0
    call get_param(param_file, "MOM", "DTBT_RESET_PERIOD", cs%dtbt_reset_period, &
                 "The period between recalculations of DTBT (if DTBT <= 0). \n"//&
                 "If DTBT_RESET_PERIOD is negative, DTBT is set based \n"//&
                 "only on information available at initialization.  If \n"//&
                 "dynamic, DTBT will be set at least every forcing time \n"//&
                 "step, and if 0, every dynamics time step.  The default is \n"//&
                 "set by DT_THERM.  This is only used if SPLIT is true.", &
                 units="s", default=default_val, do_not_read=(dtbt > 0.0))
  endif

  ! This is here in case these values are used inappropriately.
  cs%use_frazil = .false. ; cs%bound_salinity = .false. ; cs%tv%P_Ref = 2.0e7
  if (cs%use_temperature) then
    call get_param(param_file, "MOM", "FRAZIL", cs%use_frazil, &
                 "If true, water freezes if it gets too cold, and the \n"//&
                 "the accumulated heat deficit is returned in the \n"//&
                 "surface state.  FRAZIL is only used if \n"//&
                 "ENABLE_THERMODYNAMICS is true.", default=.false.)
    call get_param(param_file, "MOM", "DO_GEOTHERMAL", use_geothermal, &
                 "If true, apply geothermal heating.", default=.false.)
    call get_param(param_file, "MOM", "BOUND_SALINITY", cs%bound_salinity, &
                 "If true, limit salinity to being positive. (The sea-ice \n"//&
                 "model may ask for more salt than is available and \n"//&
                 "drive the salinity negative otherwise.)", default=.false.)
    call get_param(param_file, "MOM", "C_P", cs%tv%C_p, &
                 "The heat capacity of sea water, approximated as a \n"//&
                 "constant. This is only used if ENABLE_THERMODYNAMICS is \n"//&
                 "true. The default value is from the TEOS-10 definition \n"//&
                 "of conservative temperature.", units="J kg-1 K-1", &
                 default=3991.86795711963)
  endif
  if (use_eos) call get_param(param_file, "MOM", "P_REF", cs%tv%P_Ref, &
                 "The pressure that is used for calculating the coordinate \n"//&
                 "density.  (1 Pa = 1e4 dbar, so 2e7 is commonly used.) \n"//&
                 "This is only used if USE_EOS and ENABLE_THERMODYNAMICS \n"//&
                 "are true.", units="Pa", default=2.0e7)

  if (cs%bulkmixedlayer) then
    call get_param(param_file, "MOM", "NKML", nkml, &
                 "The number of sublayers within the mixed layer if \n"//&
                 "BULKMIXEDLAYER is true.", units="nondim", default=2)
    call get_param(param_file, "MOM", "NKBL", nkbl, &
                 "The number of layers that are used as variable density \n"//&
                 "buffer layers if BULKMIXEDLAYER is true.", units="nondim", &
                 default=2)
  endif

  call get_param(param_file, "MOM", "GLOBAL_INDEXING", global_indexing, &
                 "If true, use a global lateral indexing convention, so \n"//&
                 "that corresponding points on different processors have \n"//&
                 "the same index. This does not work with static memory.", &
                 default=.false., layoutparam=.true.)
#ifdef STATIC_MEMORY_
  if (global_indexing) call mom_error(fatal, "initialize_MOM: "//&
       "GLOBAL_INDEXING can not be true with STATIC_MEMORY.")
#endif
  call get_param(param_file, "MOM", "FIRST_DIRECTION", first_direction, &
                 "An integer that indicates which direction goes first \n"//&
                 "in parts of the code that use directionally split \n"//&
                 "updates, with even numbers (or 0) used for x- first \n"//&
                 "and odd numbers used for y-first.", default=0)

  call get_param(param_file, "MOM", "CHECK_BAD_SURFACE_VALS", &
                                     cs%check_bad_surface_vals, &
                 "If true, check the surface state for ridiculous values.", &
                 default=.false.)
  if (cs%check_bad_surface_vals) then
    call get_param(param_file, "MOM", "BAD_VAL_SSH_MAX", cs%bad_val_ssh_max, &
                 "The value of SSH above which a bad value message is \n"//&
                 "triggered, if CHECK_BAD_SURFACE_VALS is true.", units="m", &
                 default=20.0)
    call get_param(param_file, "MOM", "BAD_VAL_SSS_MAX", cs%bad_val_sss_max, &
                 "The value of SSS above which a bad value message is \n"//&
                 "triggered, if CHECK_BAD_SURFACE_VALS is true.", units="PPT", &
                 default=45.0)
    call get_param(param_file, "MOM", "BAD_VAL_SST_MAX", cs%bad_val_sst_max, &
                 "The value of SST above which a bad value message is \n"//&
                 "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
                 units="deg C", default=45.0)
    call get_param(param_file, "MOM", "BAD_VAL_SST_MIN", cs%bad_val_sst_min, &
                 "The value of SST below which a bad value message is \n"//&
                 "triggered, if CHECK_BAD_SURFACE_VALS is true.", &
                 units="deg C", default=-2.1)
    call get_param(param_file, "MOM", "BAD_VAL_COLUMN_THICKNESS", cs%bad_val_column_thickness, &
         "The value of column thickness below which a bad value message is \n"//&
         "triggered, if CHECK_BAD_SURFACE_VALS is true.", units="m", &
                          default=0.0)
  endif

  call get_param(param_file, "MOM", "SAVE_INITIAL_CONDS", save_ic, &
                 "If true, write the initial conditions to a file given \n"//&
                 "by IC_OUTPUT_FILE.", default=.false.)
  call get_param(param_file, "MOM", "IC_OUTPUT_FILE", cs%IC_file, &
                 "The file into which to write the initial conditions.", &
                 default="MOM_IC")
  call get_param(param_file, "MOM", "WRITE_GEOM", write_geom, &
                 "If =0, never write the geometry and vertical grid files.\n"//&
                 "If =1, write the geometry and vertical grid files only for\n"//&
                 "a new simulation. If =2, always write the geometry and\n"//&
                 "vertical grid files. Other values are invalid.", default=1)
  if (write_geom<0 .or. write_geom>2) call mom_error(fatal,"MOM: "//&
         "WRITE_GEOM must be equal to 0, 1 or 2.")
  write_geom_files = ((write_geom==2) .or. ((write_geom==1) .and. &
     ((dirs%input_filename(1:1)=='n') .and. (len_trim(dirs%input_filename)==1))))

  ! Check for inconsistent parameter settings.
  if (cs%use_ALE_algorithm .and. cs%bulkmixedlayer) call mom_error(fatal, &
    "MOM: BULKMIXEDLAYER can not currently be used with the ALE algorithm.")
  if (cs%use_ALE_algorithm .and. .not.cs%use_temperature) call mom_error(fatal, &
     "MOM: At this time, USE_EOS should be True when using the ALE algorithm.")
  if (cs%adiabatic .and. cs%use_temperature) call mom_error(warning, &
    "MOM: ADIABATIC and ENABLE_THERMODYNAMICS both defined is usually unwise.")
  if (use_eos .and. .not.cs%use_temperature) call mom_error(fatal, &
    "MOM: ENABLE_THERMODYNAMICS must be defined to use USE_EOS.")
  if (cs%adiabatic .and. cs%bulkmixedlayer) call mom_error(fatal, &
    "MOM: ADIABATIC and BULKMIXEDLAYER can not both be defined.")
  if (cs%bulkmixedlayer .and. .not.use_eos) call mom_error(fatal, &
      "initialize_MOM: A bulk mixed layer can only be used with T & S as "//&
      "state variables. Add USE_EOS = True to MOM_input.")

  call get_param(param_file, 'MOM', "ICE_SHELF", use_ice_shelf, default=.false., do_not_log=.true.)
  if (use_ice_shelf) then
     inputdir = "." ;  call get_param(param_file, 'MOM', "INPUTDIR", inputdir)
     inputdir = slasher(inputdir)
     call get_param(param_file, 'MOM', "ICE_THICKNESS_FILE", ice_shelf_file, &
                    "The file from which the ice bathymetry and area are read.", &
                    fail_if_missing=.true.)
     call get_param(param_file, 'MOM', "ICE_AREA_VARNAME", area_varname, &
                    "The name of the area variable in ICE_THICKNESS_FILE.", &
                    fail_if_missing=.true.)
  endif


  call calltree_waypoint("MOM parameters read (initialize_MOM)")

  ! Set up the model domain and grids.
#ifdef SYMMETRIC_MEMORY_
  symmetric = .true.
#else
  symmetric = .false.
#endif
#ifdef STATIC_MEMORY_
  call mom_domains_init(g%domain, param_file, symmetric=symmetric, &
            static_memory=.true., nihalo=nihalo_, njhalo=njhalo_, &
            niglobal=niglobal_, njglobal=njglobal_, niproc=niproc_, &
            njproc=njproc_)
#else
  call mom_domains_init(g%domain, param_file, symmetric=symmetric)
#endif
  call calltree_waypoint("domains initialized (initialize_MOM)")

  call mom_debugging_init(param_file)
  call diag_mediator_infrastructure_init()
  call mom_io_init(param_file)

  call hor_index_init(g%Domain, hi, param_file, &
                      local_indexing=.not.global_indexing)

  call create_dyn_horgrid(dg, hi, bathymetry_at_vel=bathy_at_vel)
  call clone_mom_domain(g%Domain, dg%Domain)

  call verticalgridinit( param_file, cs%GV )
  gv => cs%GV
!  dG%g_Earth = GV%g_Earth

  ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
  if (cs%debug .or. dg%symmetric) &
    call clone_mom_domain(dg%Domain, dg%Domain_aux, symmetric=.false.)

  call calltree_waypoint("grids initialized (initialize_MOM)")


  call mom_timing_init(cs)

  call tracer_registry_init(param_file, cs%tracer_Reg)

  is   = dg%isc   ; ie   = dg%iec  ; js   = dg%jsc  ; je   = dg%jec ; nz = gv%ke
  isd  = dg%isd   ; ied  = dg%ied  ; jsd  = dg%jsd  ; jed  = dg%jed
  isdb = dg%IsdB  ; iedb = dg%IedB ; jsdb = dg%JsdB ; jedb = dg%JedB

  ! Allocate and initialize space for primary MOM variables.
  alloc_(cs%u(isdb:iedb,jsd:jed,nz))   ; cs%u(:,:,:) = 0.0
  alloc_(cs%v(isd:ied,jsdb:jedb,nz))   ; cs%v(:,:,:) = 0.0
  alloc_(cs%h(isd:ied,jsd:jed,nz))     ; cs%h(:,:,:) = gv%Angstrom
  alloc_(cs%uh(isdb:iedb,jsd:jed,nz))  ; cs%uh(:,:,:) = 0.0
  alloc_(cs%vh(isd:ied,jsdb:jedb,nz))  ; cs%vh(:,:,:) = 0.0
  if (cs%use_temperature) then
    alloc_(cs%T(isd:ied,jsd:jed,nz))   ; cs%T(:,:,:) = 0.0
    alloc_(cs%S(isd:ied,jsd:jed,nz))   ; cs%S(:,:,:) = 0.0
    cs%tv%T => cs%T ; cs%tv%S => cs%S
    cs%vd_T = var_desc(name="T",units="degC",longname="Potential Temperature", &
                       cmor_field_name="thetao",cmor_units="C",                &
                       conversion=cs%tv%C_p)
    cs%vd_S = var_desc(name="S",units="PPT",longname="Salinity",&
                       cmor_field_name="so",cmor_units="ppt",   &
                       conversion=0.001)
    if(cs%advect_TS) then
      call register_tracer(cs%tv%T, cs%vd_T, param_file, dg%HI, gv, cs%tracer_Reg, cs%vd_T)
      call register_tracer(cs%tv%S, cs%vd_S, param_file, dg%HI, gv, cs%tracer_Reg, cs%vd_S)
    endif
  endif
  if (cs%use_frazil) then
    allocate(cs%tv%frazil(isd:ied,jsd:jed)) ; cs%tv%frazil(:,:) = 0.0
  endif
  if (cs%bound_salinity) then
    allocate(cs%tv%salt_deficit(isd:ied,jsd:jed)) ; cs%tv%salt_deficit(:,:)=0.0
  endif

  if (cs%bulkmixedlayer .or. cs%use_temperature) then
     allocate(cs%Hml(isd:ied,jsd:jed)) ; cs%Hml(:,:) = 0.0
  endif

  if (cs%bulkmixedlayer) then
    gv%nkml = nkml ; gv%nk_rho_varies = nkml + nkbl
  else
    gv%nkml = 0 ; gv%nk_rho_varies = 0
  endif
  if (cs%use_ALE_algorithm) then
    call get_param(param_file, "MOM", "NK_RHO_VARIES", gv%nk_rho_varies, default=0) ! Will default to nz later... -AJA
  endif

  alloc_(cs%uhtr(isdb:iedb,jsd:jed,nz)) ; cs%uhtr(:,:,:) = 0.0
  alloc_(cs%vhtr(isd:ied,jsdb:jedb,nz)) ; cs%vhtr(:,:,:) = 0.0
  cs%t_dyn_rel_adv = 0.0 ; cs%t_dyn_rel_thermo = 0.0 ; cs%t_dyn_rel_diag = 0.0

  if (cs%debug_truncations) then
    allocate(cs%u_prev(isdb:iedb,jsd:jed,nz)) ; cs%u_prev(:,:,:) = 0.0
    allocate(cs%v_prev(isd:ied,jsdb:jedb,nz)) ; cs%u_prev(:,:,:) = 0.0
  endif

  mom_internal_state%u => cs%u ; mom_internal_state%v => cs%v
  mom_internal_state%h => cs%h
  mom_internal_state%uh => cs%uh ; mom_internal_state%vh => cs%vh
  if (cs%use_temperature) then
    mom_internal_state%T => cs%T ; mom_internal_state%S => cs%S
  endif

  cs%CDp%uh => cs%uh ; cs%CDp%vh => cs%vh

  if (cs%interp_p_surf) then
    allocate(cs%p_surf_prev(isd:ied,jsd:jed)) ; cs%p_surf_prev(:,:) = 0.0
  endif

  alloc_(cs%ave_ssh(isd:ied,jsd:jed)) ; cs%ave_ssh(:,:) = 0.0

  ! Use the Wright equation of state by default, unless otherwise specified
  ! Note: this line and the following block ought to be in a separate
  ! initialization routine for tv.
  if (use_eos) call eos_init(param_file, cs%tv%eqn_of_state)
  if (cs%use_temperature) then
    allocate(cs%tv%TempxPmE(isd:ied,jsd:jed))
    cs%tv%TempxPmE(:,:) = 0.0
    if (use_geothermal) then
      allocate(cs%tv%internal_heat(isd:ied,jsd:jed))
      cs%tv%internal_heat(:,:) = 0.0
    endif
  endif
  call calltree_waypoint("state variables allocated (initialize_MOM)")

  ! Set the fields that are needed for bitwise identical restarting
  ! the time stepping scheme.
  call restart_init(param_file, cs%restart_CSp)
  call set_restart_fields(gv, param_file, cs)
  if (cs%split) then
    if (cs%legacy_split) then
      call register_restarts_dyn_legacy_split(dg%HI, gv, param_file, &
               cs%dyn_legacy_split_CSp, cs%restart_CSp, cs%uh, cs%vh)
    else
      call register_restarts_dyn_split_rk2(dg%HI, gv, param_file, &
               cs%dyn_split_RK2_CSp, cs%restart_CSp, cs%uh, cs%vh)
    endif
  else
    if (cs%use_RK2) then
      call register_restarts_dyn_unsplit_rk2(dg%HI, gv, param_file, &
             cs%dyn_unsplit_RK2_CSp, cs%restart_CSp)
    else
      call register_restarts_dyn_unsplit(dg%HI, gv, param_file, &
             cs%dyn_unsplit_CSp, cs%restart_CSp)
    endif
  endif

  ! This subroutine calls user-specified tracer registration routines.
  ! Additional calls can be added to MOM_tracer_flow_control.F90.
  call call_tracer_register(dg%HI, gv, param_file, cs%tracer_flow_CSp, &
                            cs%tracer_Reg, cs%restart_CSp)

  call meke_alloc_register_restart(dg%HI, param_file, cs%MEKE, cs%restart_CSp)
  call set_visc_register_restarts(dg%HI, gv, param_file, cs%visc, cs%restart_CSp)
  call mixedlayer_restrat_register_restarts(dg%HI, param_file, cs%mixedlayer_restrat_CSp, cs%restart_CSp)

  ! Initialize fields
  call calltree_waypoint("restart registration complete (initialize_MOM)")

  call cpu_clock_begin(id_clock_mom_init)
  call mom_initialize_fixed(dg, cs%OBC, param_file, write_geom_files, dirs%output_directory)
  call calltree_waypoint("returned from MOM_initialize_fixed() (initialize_MOM)")

  if (associated(cs%OBC)) call call_obc_register(param_file, cs%update_OBC_CSp, cs%OBC)

  call mom_initialize_coord(gv, param_file, write_geom_files, &
                            dirs%output_directory, cs%tv, dg%max_depth)
  call calltree_waypoint("returned from MOM_initialize_coord() (initialize_MOM)")

  if (cs%use_ALE_algorithm) then
    call ale_init(param_file, gv, dg%max_depth, cs%ALE_CSp)
    call calltree_waypoint("returned from ALE_init() (initialize_MOM)")
  endif

  !   Shift from using the temporary dynamic grid type to using the final
  ! (potentially static) ocean-specific grid type.
  !   The next line would be needed if G%Domain had not already been init'd above:
  !     call clone_MOM_domain(dG%Domain, G%Domain)
  call mom_grid_init(g, param_file, hi, bathymetry_at_vel=bathy_at_vel)
  call copy_dyngrid_to_mom_grid(dg, g)
  call destroy_dyn_horgrid(dg)

  ! Set a few remaining fields that are specific to the ocean grid type.
  call set_first_direction(g, first_direction)
  ! Allocate the auxiliary non-symmetric domain for debugging or I/O purposes.
  if (cs%debug .or. g%symmetric) &
    call clone_mom_domain(g%Domain, g%Domain_aux, symmetric=.false.)
  ! Copy common variables from the vertical grid to the horizontal grid.
  ! Consider removing this later?
  g%ke = gv%ke ; g%g_Earth = gv%g_Earth

  call mom_initialize_state(cs%u, cs%v, cs%h, cs%tv, time, g, gv, param_file, &
                            dirs, cs%restart_CSp, cs%ALE_CSp, cs%tracer_Reg, &
                            cs%sponge_CSp, cs%ALE_sponge_CSp, cs%OBC, time_in)
  call cpu_clock_end(id_clock_mom_init)
  call calltree_waypoint("returned from MOM_initialize_state() (initialize_MOM)")

  ! From this point, there may be pointers being set, so the final grid type
  ! that will persist throughout the run has to be used.

  if (test_grid_copy) then
    !  Copy the data from the temporary grid to the dyn_hor_grid to CS%G.
    call create_dyn_horgrid(dg, g%HI)
    call clone_mom_domain(g%Domain, dg%Domain)

    call clone_mom_domain(g%Domain, cs%G%Domain)
    call mom_grid_init(cs%G, param_file)

    call copy_mom_grid_to_dyngrid(g, dg)
    call copy_dyngrid_to_mom_grid(dg, cs%G)

    call destroy_dyn_horgrid(dg)
    call mom_grid_end(g) ; deallocate(g)

    g => cs%G
    if (cs%debug .or. cs%G%symmetric) &
      call clone_mom_domain(cs%G%Domain, cs%G%Domain_aux, symmetric=.false.)
    g%ke = gv%ke ; g%g_Earth = gv%g_Earth
  endif


  ! At this point, all user-modified initialization code has been called.  The
  ! remainder of this subroutine is controlled by the parameters that have
  ! have already been set.

  if (ale_remap_init_conds(cs%ALE_CSp) .and. .not. query_initialized(cs%h,"h",cs%restart_CSp)) then
    ! This block is controlled by the ALE parameter REMAP_AFTER_INITIALIZATION.
    ! \todo This block exists for legacy reasons and we should phase it out of
    ! all examples.
    if (cs%debug) then
      call uvchksum("Pre ALE adjust init cond [uv]", &
                    cs%u, cs%v, g%HI, haloshift=1)
      call hchksum(cs%h,"Pre ALE adjust init cond h", g%HI, haloshift=1, scale=gv%H_to_m)
    endif
    call calltree_waypoint("Calling adjustGridForIntegrity() to remap initial conditions (initialize_MOM)")
    call adjustgridforintegrity(cs%ALE_CSp, g, gv, cs%h )
    call calltree_waypoint("Calling ALE_main() to remap initial conditions (initialize_MOM)")
    if (use_ice_shelf) then
      filename = trim(inputdir)//trim(ice_shelf_file)
      if (.not.file_exists(filename, g%Domain)) call mom_error(fatal, &
        "MOM: Unable to open "//trim(filename))

      allocate(area_shelf_h(isd:ied,jsd:jed))
      allocate(frac_shelf_h(isd:ied,jsd:jed))
      call read_data(filename,trim(area_varname),area_shelf_h,domain=g%Domain%mpp_domain)
      ! initialize frac_shelf_h with zeros (open water everywhere)
      frac_shelf_h(:,:) = 0.0
      ! compute fractional ice shelf coverage of h
      do j=jsd,jed ; do i=isd,ied
        if (g%areaT(i,j) > 0.0) &
          frac_shelf_h(i,j) = area_shelf_h(i,j) / g%areaT(i,j)
      enddo ; enddo
      ! pass to the pointer
      shelf_area => frac_shelf_h
      call ale_main(g, gv, cs%h, cs%u, cs%v, cs%tv, cs%tracer_Reg, cs%ALE_CSp, &
                    frac_shelf_h = shelf_area)
    else
      call ale_main( g, gv, cs%h, cs%u, cs%v, cs%tv, cs%tracer_Reg, cs%ALE_CSp )
    endif
    call cpu_clock_begin(id_clock_pass_init)
    call create_group_pass(tmp_pass_uv_t_s_h, cs%u, cs%v, g%Domain)
    if (cs%use_temperature) then
      call create_group_pass(tmp_pass_uv_t_s_h, cs%tv%T, g%Domain, halo=1)
      call create_group_pass(tmp_pass_uv_t_s_h, cs%tv%S, g%Domain, halo=1)
    endif
    call create_group_pass(tmp_pass_uv_t_s_h, cs%h, g%Domain, halo=1)
    call do_group_pass(tmp_pass_uv_t_s_h, g%Domain)
    call cpu_clock_end(id_clock_pass_init)

    if (cs%debug) then
      call uvchksum("Post ALE adjust init cond [uv]", cs%u, cs%v, g%HI, haloshift=1)
      call hchksum(cs%h, "Post ALE adjust init cond h", g%HI, haloshift=1, scale=gv%H_to_m)
    endif
  endif
  if ( cs%use_ALE_algorithm ) call ale_updateverticalgridtype( cs%ALE_CSp, gv )

   diag    => cs%diag
  ! Initialize the diag mediator.
  call diag_mediator_init(g, gv%ke, param_file, diag, doc_file_dir=dirs%output_directory)

  ! Initialize the diagnostics mask arrays.
  ! This step has to be done after call to MOM_initialize_state
  ! and before MOM_diagnostics_init
  call diag_masks_set(g, gv%ke, cs%missing, diag)

  ! Set up pointers within diag mediator control structure,
  ! this needs to occur _after_ CS%h etc. have been allocated.
  call diag_set_state_ptrs(cs%h, cs%T, cs%S, cs%tv%eqn_of_state, diag)

  ! This call sets up the diagnostic axes. These are needed,
  ! e.g. to generate the target grids below.
  call set_axes_info(g, gv, param_file, diag)

  ! Whenever thickness/T/S changes let the diag manager know, target grids
  ! for vertical remapping may need to be regenerated.
  ! FIXME: are h, T, S updated at the same time? Review these for T, S updates.
  call diag_update_remap_grids(diag)

  ! Diagnose static fields AND associate areas/volumes with axes
  call write_static_fields(g, cs%diag)
  call calltree_waypoint("static fields written (initialize_MOM)")

  call cpu_clock_begin(id_clock_mom_init)
  if (cs%use_ALE_algorithm) then
    call ale_writecoordinatefile( cs%ALE_CSp, gv, dirs%output_directory )
  endif
  call cpu_clock_end(id_clock_mom_init)
  call calltree_waypoint("ALE initialized (initialize_MOM)")

  cs%useMEKE = meke_init(time, g, param_file, diag, cs%MEKE_CSp, cs%MEKE, cs%restart_CSp)

  call varmix_init(time, g, param_file, diag, cs%VarMix)
  call set_visc_init(time, g, gv, param_file, diag, cs%visc, cs%set_visc_CSp,cs%OBC)
  if (cs%split) then
    allocate(eta(szi_(g),szj_(g))) ; eta(:,:) = 0.0
    if (cs%legacy_split) then
      call initialize_dyn_legacy_split(cs%u, cs%v, cs%h, cs%uh, cs%vh, eta, time,   &
                  g, gv, param_file, diag, cs%dyn_legacy_split_CSp, cs%restart_CSp, &
                  cs%dt, cs%ADp, cs%CDp, mom_internal_state, cs%VarMix, cs%MEKE,    &
                  cs%OBC, cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp, cs%visc,  &
                  dirs, cs%ntrunc)
    else
      call initialize_dyn_split_rk2(cs%u, cs%v, cs%h, cs%uh, cs%vh, eta, time,   &
                  g, gv, param_file, diag, cs%dyn_split_RK2_CSp, cs%restart_CSp, &
                  cs%dt, cs%ADp, cs%CDp, mom_internal_state, cs%VarMix, cs%MEKE, &
                  cs%OBC, cs%update_OBC_CSp, cs%ALE_CSp, cs%set_visc_CSp,        &
                  cs%visc, dirs, cs%ntrunc)
    endif
  else
    if (cs%use_RK2) then
      call initialize_dyn_unsplit_rk2(cs%u, cs%v, cs%h, time, g, gv,         &
              param_file, diag, cs%dyn_unsplit_RK2_CSp, cs%restart_CSp,      &
              cs%ADp, cs%CDp, mom_internal_state, cs%OBC, cs%update_OBC_CSp, &
              cs%ALE_CSp, cs%set_visc_CSp, cs%visc, dirs, cs%ntrunc)
    else
      call initialize_dyn_unsplit(cs%u, cs%v, cs%h, time, g, gv,             &
              param_file, diag, cs%dyn_unsplit_CSp, cs%restart_CSp,          &
              cs%ADp, cs%CDp, mom_internal_state, cs%OBC, cs%update_OBC_CSp, &
              cs%ALE_CSp, cs%set_visc_CSp, cs%visc, dirs, cs%ntrunc)
    endif
  endif
  call calltree_waypoint("dynamics initialized (initialize_MOM)")

  call thickness_diffuse_init(time, g, gv, param_file, diag, cs%CDp, cs%thickness_diffuse_CSp)
  cs%mixedlayer_restrat = mixedlayer_restrat_init(time, g, gv, param_file, diag, &
                                                  cs%mixedlayer_restrat_CSp)
  if (cs%mixedlayer_restrat) then
    if (.not.(cs%bulkmixedlayer .or. cs%use_ALE_algorithm)) &
      call mom_error(fatal, "MOM: MIXEDLAYER_RESTRAT true requires a boundary layer scheme.")
    ! When DIABATIC_FIRST=False and using CS%visc%ML in mixedlayer_restrat we need to update after a restart
    if (.not. cs%diabatic_first .and. associated(cs%visc%MLD)) &
      call pass_var(cs%visc%MLD, g%domain, halo=1)
  endif

  call mom_diagnostics_init(mom_internal_state, cs%ADp, cs%CDp, time, g, gv, &
                            param_file, diag, cs%diagnostics_CSp)

  cs%Z_diag_interval = set_time(int((cs%dt_therm) * &
       max(1,floor(0.01 + z_diag_int/(cs%dt_therm)))))
  call mom_diag_to_z_init(time, g, gv, param_file, diag, cs%diag_to_Z_CSp)
  cs%Z_diag_time = start_time + cs%Z_diag_interval * (1 + &
    ((time + set_time(int(cs%dt_therm))) - start_time) / cs%Z_diag_interval)

  if (associated(cs%sponge_CSp)) &
    call init_sponge_diags(time, g, diag, cs%sponge_CSp)

  if (associated(cs%ALE_sponge_CSp)) &
    call init_ale_sponge_diags(time, g, diag, cs%ALE_sponge_CSp)

  if (cs%adiabatic) then
    call adiabatic_driver_init(time, g, param_file, diag, cs%diabatic_CSp, &
                               cs%tracer_flow_CSp, cs%diag_to_Z_CSp)
  else
    call diabatic_driver_init(time, g, gv, param_file, cs%use_ALE_algorithm, diag,     &
                              cs%ADp, cs%CDp, cs%diabatic_CSp, cs%tracer_flow_CSp, &
                              cs%sponge_CSp, cs%ALE_sponge_CSp, cs%diag_to_Z_CSp)
  endif

  call tracer_advect_init(time, g, param_file, diag, cs%tracer_adv_CSp)
  call tracer_hor_diff_init(time, g, param_file, diag, cs%tracer_diff_CSp, cs%neutral_diffusion_CSp)

  if (cs%use_ALE_algorithm) &
    call register_diags_ts_vardec(time, g%HI, gv, param_file, cs)

  call lock_tracer_registry(cs%tracer_Reg)
  call calltree_waypoint("tracer registry now locked (initialize_MOM)")

  ! now register some diagnostics since tracer registry is locked
  call register_diags(time, g, gv, cs, cs%ADp)
  call register_diags_ts_tendency(time, g, cs)
  if (cs%use_ALE_algorithm) then
    call ale_register_diags(time, g, diag, cs%tv%C_p, cs%tracer_Reg, cs%ALE_CSp)
  endif



  ! If need a diagnostic field, then would have been allocated in register_diags.
  if (cs%use_temperature) then
    if(cs%advect_TS) then
      call add_tracer_diagnostics("T", cs%tracer_Reg, cs%T_adx, cs%T_ady, &
                        cs%T_diffx, cs%T_diffy, cs%T_adx_2d, cs%T_ady_2d, &
                        cs%T_diffx_2d, cs%T_diffy_2d, cs%T_advection_xy)
      call add_tracer_diagnostics("S", cs%tracer_Reg, cs%S_adx, cs%S_ady, &
                        cs%S_diffx, cs%S_diffy, cs%S_adx_2d, cs%S_ady_2d, &
                        cs%S_diffx_2d, cs%S_diffy_2d, cs%S_advection_xy)
    endif
    call register_z_tracer(cs%tv%T, "temp", "Potential Temperature", "degC", time,   &
                      g, cs%diag_to_Z_CSp, cmor_field_name="thetao", cmor_units="C", &
                      cmor_standard_name="sea_water_potential_temperature",          &
                      cmor_long_name ="Sea Water Potential Temperature")
    call register_z_tracer(cs%tv%S, "salt", "Salinity", "PPT", time,               &
                      g, cs%diag_to_Z_CSp, cmor_field_name="so", cmor_units="ppt", &
                      cmor_standard_name="sea_water_salinity",                     &
                      cmor_long_name ="Sea Water Salinity")
  endif

  ! This subroutine initializes any tracer packages.
  new_sim = ((dirs%input_filename(1:1) == 'n') .and. &
             (len_trim(dirs%input_filename) == 1))
  call tracer_flow_control_init(.not.new_sim, time, g, gv, cs%h, param_file, &
             cs%diag, cs%OBC, cs%tracer_flow_CSp, cs%sponge_CSp, &
             cs%ALE_sponge_CSp, cs%diag_to_Z_CSp, cs%tv)


  ! If running in offline tracer mode, initialize the necessary control structure and
  ! parameters
  if(present(offline_tracer_mode)) offline_tracer_mode=cs%offline_tracer_mode

  if(cs%offline_tracer_mode) then
    ! Setup some initial parameterizations and also assign some of the subtypes
    call offline_transport_init(param_file, cs%offline_CSp, cs%diabatic_CSp, g, gv)
    call insert_offline_main( cs=cs%offline_CSp, ale_csp=cs%ALE_CSp, diabatic_csp=cs%diabatic_CSp, &
                              diag=cs%diag, obc=cs%OBC, tracer_adv_csp=cs%tracer_adv_CSp,              &
                              tracer_flow_csp=cs%tracer_flow_CSp, tracer_reg=cs%tracer_Reg,            &
                              tv=cs%tv, x_before_y = (mod(first_direction,2)==0), debug=cs%debug )
    call register_diags_offline_transport(time, cs%diag, cs%offline_CSp)
  endif

  call cpu_clock_begin(id_clock_pass_init)
  !--- set up group pass for u,v,T,S and h. pass_uv_T_S_h also is used in step_MOM

  dynamics_stencil = min(3, g%Domain%nihalo, g%Domain%njhalo)
  call create_group_pass(cs%pass_uv_T_S_h, cs%u, cs%v, g%Domain, halo=dynamics_stencil)
  if (cs%use_temperature) then
    call create_group_pass(cs%pass_uv_T_S_h, cs%tv%T, g%Domain, halo=dynamics_stencil)
    call create_group_pass(cs%pass_uv_T_S_h, cs%tv%S, g%Domain, halo=dynamics_stencil)
  endif
  call create_group_pass(cs%pass_uv_T_S_h, cs%h, g%Domain, halo=dynamics_stencil)

  call do_group_pass(cs%pass_uv_T_S_h, g%Domain)
  call cpu_clock_end(id_clock_pass_init)

  call register_obsolete_diagnostics(param_file, cs%diag)
  call neutral_diffusion_diag_init(time, g, diag, cs%tv%C_p, cs%tracer_Reg, cs%neutral_diffusion_CSp)

  if (cs%use_frazil) then
    if (.not.query_initialized(cs%tv%frazil,"frazil",cs%restart_CSp)) &
      cs%tv%frazil(:,:) = 0.0
  endif

  if (cs%interp_p_surf) then
    cs%p_surf_prev_set = &
      query_initialized(cs%p_surf_prev,"p_surf_prev",cs%restart_CSp)

    if (cs%p_surf_prev_set) call pass_var(cs%p_surf_prev, g%domain)
  endif

  if (.not.query_initialized(cs%ave_ssh,"ave_ssh",cs%restart_CSp)) then
    if (cs%split) then
      call find_eta(cs%h, cs%tv, gv%g_Earth, g, gv, cs%ave_ssh, eta)
    else
      call find_eta(cs%h, cs%tv, gv%g_Earth, g, gv, cs%ave_ssh)
    endif
  endif
  if (cs%split) deallocate(eta)

  ! Flag whether to save initial conditions in finish_MOM_initialization() or not.
  cs%write_IC = save_ic .and. &
                .not.((dirs%input_filename(1:1) == 'r') .and. &
                      (len_trim(dirs%input_filename) == 1))

  call calltree_leave("initialize_MOM()")
  call cpu_clock_end(id_clock_init)

end subroutine initialize_mom

!> This subroutine finishes initializing MOM and writes out the initial conditions.
subroutine finish_mom_initialization(Time, dirs, CS, fluxes)
  type(time_type),           intent(in)    :: Time        !< model time, used in this routine
  type(directories),         intent(in)    :: dirs        !< structure with directory paths
  type(mom_control_struct),  pointer       :: CS          !< pointer set in this routine to MOM control structure
  type(forcing),             intent(inout) :: fluxes      !< pointers to forcing fields
  ! Local variables
  type(ocean_grid_type), pointer :: G => null()
  type(verticalgrid_type), pointer :: GV => null()
  type(mom_restart_cs), pointer :: restart_CSp_tmp => null()
  real, allocatable :: z_interface(:,:,:) ! Interface heights (meter)
  real, allocatable :: eta(:,:) ! Interface heights (meter)
  type(vardesc) :: vd

  call cpu_clock_begin(id_clock_init)
  call calltree_enter("finish_MOM_initialization()")

  ! Pointers for convenience
  g => cs%G ; gv => cs%GV

  ! Write initial conditions
  if (cs%write_IC) then
    allocate(restart_csp_tmp)
    restart_csp_tmp = cs%restart_CSp
    allocate(z_interface(szi_(g),szj_(g),szk_(g)+1))
    call find_eta(cs%h, cs%tv, gv%g_Earth, g, gv, z_interface)
    vd = var_desc("eta","meter","Interface heights",z_grid='i')
    call register_restart_field(z_interface, vd, .true., restart_csp_tmp)

    call save_restart(dirs%output_directory, time, g, &
                      restart_csp_tmp, filename=cs%IC_file, gv=gv)
    deallocate(z_interface)
    deallocate(restart_csp_tmp)
  endif

  call calltree_leave("finish_MOM_initialization()")
  call cpu_clock_end(id_clock_init)

end subroutine finish_mom_initialization

!> Register the diagnostics
subroutine register_diags(Time, G, GV, CS, ADp)
  type(time_type),           intent(in)    :: Time  !< current model time
  type(ocean_grid_type),     intent(inout) :: G     !< ocean grid structu
  type(verticalgrid_type),   intent(inout) :: GV    !< ocean vertical grid structure
  type(mom_control_struct),  pointer       :: CS    !< control structure set up by initialize_MOM
  type(accel_diag_ptrs),     intent(inout) :: ADp   !< structure pointing to accelerations in momentum equation

  character(len=48) :: thickness_units, flux_units, T_flux_units, S_flux_units
  type(diag_ctrl), pointer :: diag
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
  isd  = g%isd  ; ied  = g%ied  ; jsd  = g%jsd  ; jed  = g%jed ; nz = g%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  diag => cs%diag

  thickness_units = get_thickness_units(gv)
  flux_units      = get_flux_units(gv)
  t_flux_units    = get_tr_flux_units(gv, "Celsius")
  s_flux_units    = get_tr_flux_units(gv, "PPT")

  !Initialize the diagnostics mask arrays.
  !This has to be done after MOM_initialize_state call.
  !call diag_masks_set(G, CS%missing)

  cs%id_u = register_diag_field('ocean_model', 'u', diag%axesCuL, time,              &
      'Zonal velocity', 'meter  second-1', cmor_field_name='uo', cmor_units='m s-1', &
      cmor_standard_name='sea_water_x_velocity', cmor_long_name='Sea Water X Velocity')
  cs%id_v = register_diag_field('ocean_model', 'v', diag%axesCvL, time,                  &
      'Meridional velocity', 'meter second-1', cmor_field_name='vo', cmor_units='m s-1', &
      cmor_standard_name='sea_water_y_velocity', cmor_long_name='Sea Water Y Velocity')
  cs%id_h = register_diag_field('ocean_model', 'h', diag%axesTL, time, &
      'Layer Thickness', thickness_units, v_extensive=.true.)

  cs%id_volo = register_scalar_field('ocean_model', 'volo', time, diag,&
      long_name='Total volume of liquid ocean', units='m3',            &
      standard_name='sea_water_volume')
  cs%id_zos = register_diag_field('ocean_model', 'zos', diag%axesT1, time,&
      standard_name = 'sea_surface_height_above_geoid',                   &
      long_name= 'Sea surface height above geoid', units='meter', missing_value=cs%missing)
  cs%id_zossq = register_diag_field('ocean_model', 'zossq', diag%axesT1, time,&
      standard_name='square_of_sea_surface_height_above_geoid',             &
      long_name='Square of sea surface height above geoid', units='m2', missing_value=cs%missing)
  cs%id_ssh = register_diag_field('ocean_model', 'SSH', diag%axesT1, time, &
      'Sea Surface Height', 'meter', cs%missing)
  cs%id_ssh_ga = register_scalar_field('ocean_model', 'ssh_ga', time, diag,&
      long_name='Area averaged sea surface height', units='m',            &
      standard_name='area_averaged_sea_surface_height')
  cs%id_ssh_inst = register_diag_field('ocean_model', 'SSH_inst', diag%axesT1, time, &
      'Instantaneous Sea Surface Height', 'meter', cs%missing)
  cs%id_ssu = register_diag_field('ocean_model', 'SSU', diag%axesCu1, time, &
      'Sea Surface Zonal Velocity', 'meter second-1', cs%missing)
  cs%id_ssv = register_diag_field('ocean_model', 'SSV', diag%axesCv1, time, &
      'Sea Surface Meridional Velocity', 'meter second-1', cs%missing)
  cs%id_speed = register_diag_field('ocean_model', 'speed', diag%axesT1, time, &
      'Sea Surface Speed', 'meter second-1', cs%missing)

  if (cs%use_temperature) then
    cs%id_T = register_diag_field('ocean_model', 'temp', diag%axesTL, time, &
        'Potential Temperature', 'Celsius',                                 &
         cmor_field_name="thetao", cmor_units="C",                          &
         cmor_standard_name="sea_water_potential_temperature",              &
         cmor_long_name ="Sea Water Potential Temperature")
    cs%id_S = register_diag_field('ocean_model', 'salt', diag%axesTL, time, &
        long_name='Salinity', units='PPT', cmor_field_name='so',            &
        cmor_long_name='Sea Water Salinity', cmor_units='ppt',              &
        cmor_standard_name='sea_water_salinity')
    cs%id_tob = register_diag_field('ocean_model','tob', diag%axesT1, time,          &
        long_name='Sea Water Potential Temperature at Sea Floor',                    &
        standard_name='sea_water_potential_temperature_at_sea_floor', units='degC')
    cs%id_sob = register_diag_field('ocean_model','sob',diag%axesT1, time,           &
        long_name='Sea Water Salinity at Sea Floor',                                 &
        standard_name='sea_water_salinity_at_sea_floor', units='ppt')
    cs%id_sst = register_diag_field('ocean_model', 'SST', diag%axesT1, time,     &
        'Sea Surface Temperature', 'Celsius', cs%missing, cmor_field_name='tos', &
        cmor_long_name='Sea Surface Temperature', cmor_units='degC',             &
        cmor_standard_name='sea_surface_temperature')
    cs%id_sst_sq = register_diag_field('ocean_model', 'SST_sq', diag%axesT1, time, &
        'Sea Surface Temperature Squared', 'Celsius**2', cs%missing, cmor_field_name='tossq', &
        cmor_long_name='Square of Sea Surface Temperature ', cmor_units='degC^2', &
        cmor_standard_name='square_of_sea_surface_temperature')
    cs%id_sss = register_diag_field('ocean_model', 'SSS', diag%axesT1, time, &
        'Sea Surface Salinity', 'PPT', cs%missing, cmor_field_name='sos', &
        cmor_long_name='Sea Surface Salinity', cmor_units='ppt',          &
        cmor_standard_name='sea_surface_salinity')
    cs%id_sss_sq = register_diag_field('ocean_model', 'SSS_sq', diag%axesT1, time, &
        'Sea Surface Salinity Squared', 'ppt**2', cs%missing, cmor_field_name='sossq', &
        cmor_long_name='Square of Sea Surface Salinity ', cmor_units='ppt^2', &
        cmor_standard_name='square_of_sea_surface_salinity')
    cs%id_Tcon = register_diag_field('ocean_model', 'contemp', diag%axesTL, time, &
        'Conservative Temperature', 'Celsius')
    cs%id_Sabs = register_diag_field('ocean_model', 'abssalt', diag%axesTL, time, &
        long_name='Absolute Salinity', units='g/Kg')
    cs%id_sstcon = register_diag_field('ocean_model', 'conSST', diag%axesT1, time,     &
        'Sea Surface Conservative Temperature', 'Celsius', cs%missing)
    cs%id_sssabs = register_diag_field('ocean_model', 'absSSS', diag%axesT1, time,     &
        'Sea Surface Absolute Salinity', 'g/Kg', cs%missing)

  endif

  if (cs%use_temperature .and. cs%use_frazil) then
    cs%id_fraz = register_diag_field('ocean_model', 'frazil', diag%axesT1, time,                         &
          'Heat from frazil formation', 'Watt meter-2', cmor_field_name='hfsifrazil',                    &
          cmor_units='W m-2', cmor_standard_name='heat_flux_into_sea_water_due_to_frazil_ice_formation', &
          cmor_long_name='Heat Flux into Sea Water due to Frazil Ice Formation')
  endif

  cs%id_salt_deficit = register_diag_field('ocean_model', 'salt_deficit', diag%axesT1, time, &
         'Salt sink in ocean due to ice flux', 'g Salt meter-2 s-1')
  cs%id_Heat_PmE = register_diag_field('ocean_model', 'Heat_PmE', diag%axesT1, time, &
         'Heat flux into ocean from mass flux into ocean', 'Watt meter-2')
  cs%id_intern_heat = register_diag_field('ocean_model', 'internal_heat', diag%axesT1, time,&
         'Heat flux into ocean from geothermal or other internal sources', 'Watt meter-2')


  ! lateral heat advective and diffusive fluxes
  cs%id_Tadx = register_diag_field('ocean_model', 'T_adx', diag%axesCuL, time, &
      'Advective (by residual mean) Zonal Flux of Potential Temperature', t_flux_units)
  cs%id_Tady = register_diag_field('ocean_model', 'T_ady', diag%axesCvL, time, &
      'Advective (by residual mean) Meridional Flux of Potential Temperature', t_flux_units)
  cs%id_Tdiffx = register_diag_field('ocean_model', 'T_diffx', diag%axesCuL, time, &
      'Diffusive Zonal Flux of Potential Temperature', t_flux_units)
  cs%id_Tdiffy = register_diag_field('ocean_model', 'T_diffy', diag%axesCvL, time, &
      'Diffusive Meridional Flux of Potential Temperature', t_flux_units)
  if (cs%id_Tadx   > 0) call safe_alloc_ptr(cs%T_adx,isdb,iedb,jsd,jed,nz)
  if (cs%id_Tady   > 0) call safe_alloc_ptr(cs%T_ady,isd,ied,jsdb,jedb,nz)
  if (cs%id_Tdiffx > 0) call safe_alloc_ptr(cs%T_diffx,isdb,iedb,jsd,jed,nz)
  if (cs%id_Tdiffy > 0) call safe_alloc_ptr(cs%T_diffy,isd,ied,jsdb,jedb,nz)


  ! lateral salt advective and diffusive fluxes
  cs%id_Sadx = register_diag_field('ocean_model', 'S_adx', diag%axesCuL, time, &
      'Advective (by residual mean) Zonal Flux of Salinity', s_flux_units)
  cs%id_Sady = register_diag_field('ocean_model', 'S_ady', diag%axesCvL, time, &
      'Advective (by residual mean) Meridional Flux of Salinity', s_flux_units)
  cs%id_Sdiffx = register_diag_field('ocean_model', 'S_diffx', diag%axesCuL, time, &
      'Diffusive Zonal Flux of Salinity', s_flux_units)
  cs%id_Sdiffy = register_diag_field('ocean_model', 'S_diffy', diag%axesCvL, time, &
      'Diffusive Meridional Flux of Salinity', s_flux_units)
  if (cs%id_Sadx   > 0) call safe_alloc_ptr(cs%S_adx,isdb,iedb,jsd,jed,nz)
  if (cs%id_Sady   > 0) call safe_alloc_ptr(cs%S_ady,isd,ied,jsdb,jedb,nz)
  if (cs%id_Sdiffx > 0) call safe_alloc_ptr(cs%S_diffx,isdb,iedb,jsd,jed,nz)
  if (cs%id_Sdiffy > 0) call safe_alloc_ptr(cs%S_diffy,isd,ied,jsdb,jedb,nz)


  ! vertically integrated lateral heat advective and diffusive fluxes
  cs%id_Tadx_2d = register_diag_field('ocean_model', 'T_adx_2d', diag%axesCu1, time, &
      'Vertically Integrated Advective Zonal Flux of Potential Temperature', t_flux_units)
  cs%id_Tady_2d = register_diag_field('ocean_model', 'T_ady_2d', diag%axesCv1, time, &
      'Vertically Integrated Advective Meridional Flux of Potential Temperature', t_flux_units)
  cs%id_Tdiffx_2d = register_diag_field('ocean_model', 'T_diffx_2d', diag%axesCu1, time, &
      'Vertically Integrated Diffusive Zonal Flux of Potential Temperature', t_flux_units)
  cs%id_Tdiffy_2d = register_diag_field('ocean_model', 'T_diffy_2d', diag%axesCv1, time, &
      'Vertically Integrated Diffusive Meridional Flux of Potential Temperature', t_flux_units)
  if (cs%id_Tadx_2d   > 0) call safe_alloc_ptr(cs%T_adx_2d,isdb,iedb,jsd,jed)
  if (cs%id_Tady_2d   > 0) call safe_alloc_ptr(cs%T_ady_2d,isd,ied,jsdb,jedb)
  if (cs%id_Tdiffx_2d > 0) call safe_alloc_ptr(cs%T_diffx_2d,isdb,iedb,jsd,jed)
  if (cs%id_Tdiffy_2d > 0) call safe_alloc_ptr(cs%T_diffy_2d,isd,ied,jsdb,jedb)

  ! vertically integrated lateral salt advective and diffusive fluxes
  cs%id_Sadx_2d = register_diag_field('ocean_model', 'S_adx_2d', diag%axesCu1, time, &
      'Vertically Integrated Advective Zonal Flux of Salinity', s_flux_units)
  cs%id_Sady_2d = register_diag_field('ocean_model', 'S_ady_2d', diag%axesCv1, time, &
      'Vertically Integrated Advective Meridional Flux of Salinity', s_flux_units)
  cs%id_Sdiffx_2d = register_diag_field('ocean_model', 'S_diffx_2d', diag%axesCu1, time, &
      'Vertically Integrated Diffusive Zonal Flux of Salinity', s_flux_units)
  cs%id_Sdiffy_2d = register_diag_field('ocean_model', 'S_diffy_2d', diag%axesCv1, time, &
      'Vertically Integrated Diffusive Meridional Flux of Salinity', s_flux_units)
  if (cs%id_Sadx_2d   > 0) call safe_alloc_ptr(cs%S_adx_2d,isdb,iedb,jsd,jed)
  if (cs%id_Sady_2d   > 0) call safe_alloc_ptr(cs%S_ady_2d,isd,ied,jsdb,jedb)
  if (cs%id_Sdiffx_2d > 0) call safe_alloc_ptr(cs%S_diffx_2d,isdb,iedb,jsd,jed)
  if (cs%id_Sdiffy_2d > 0) call safe_alloc_ptr(cs%S_diffy_2d,isd,ied,jsdb,jedb)


  if (cs%debug_truncations) then
    call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)
    call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
    if (.not.cs%adiabatic) then
      call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)
      call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
    endif
  endif

  ! fields posted prior to dynamics step
  cs%id_h_pre_dyn = register_diag_field('ocean_model', 'h_pre_dyn', diag%axesTL, time, &
      'Layer Thickness before dynamics step', thickness_units)

  ! diagnostics for values prior to diabatic and prior to ALE
  cs%id_u_predia = register_diag_field('ocean_model', 'u_predia', diag%axesCuL, time, &
      'Zonal velocity before diabatic forcing', 'meter second-1')
  cs%id_v_predia = register_diag_field('ocean_model', 'v_predia', diag%axesCvL, time, &
      'Meridional velocity before diabatic forcing', 'meter second-1')
  cs%id_h_predia = register_diag_field('ocean_model', 'h_predia', diag%axesTL, time, &
      'Layer Thickness before diabatic forcing', thickness_units, v_extensive=.true.)
  cs%id_e_predia = register_diag_field('ocean_model', 'e_predia', diag%axesTi, time, &
      'Interface Heights before diabatic forcing', 'meter')
  if (.not. cs%adiabatic) then
    cs%id_u_preale = register_diag_field('ocean_model', 'u_preale', diag%axesCuL, time, &
        'Zonal velocity before remapping', 'meter second-1')
    cs%id_v_preale = register_diag_field('ocean_model', 'v_preale', diag%axesCvL, time, &
        'Meridional velocity before remapping', 'meter second-1')
    cs%id_h_preale = register_diag_field('ocean_model', 'h_preale', diag%axesTL, time, &
        'Layer Thickness before remapping', thickness_units, v_extensive=.true.)
    cs%id_T_preale = register_diag_field('ocean_model', 'T_preale', diag%axesTL, time, &
        'Temperature before remapping', 'degC')
    cs%id_S_preale = register_diag_field('ocean_model', 'S_preale', diag%axesTL, time, &
        'Salinity before remapping', 'ppt')
    cs%id_e_preale = register_diag_field('ocean_model', 'e_preale', diag%axesTi, time, &
        'Interface Heights before remapping', 'meter')
  endif

  if (cs%use_temperature) then
    cs%id_T_predia = register_diag_field('ocean_model', 'temp_predia', diag%axesTL, time, &
        'Potential Temperature', 'Celsius')
    cs%id_S_predia = register_diag_field('ocean_model', 'salt_predia', diag%axesTL, time, &
        'Salinity', 'PPT')
  endif

  ! Diagnostics related to tracer transport
  cs%id_uhtr = register_diag_field('ocean_model', 'uhtr', diag%axesCuL, time, &
      'Accumulated zonal thickness fluxes to advect tracers', 'kg', &
      y_cell_method='sum', v_extensive=.true.)
  cs%id_vhtr = register_diag_field('ocean_model', 'vhtr', diag%axesCvL, time, &
      'Accumulated meridional thickness fluxes to advect tracers', 'kg', &
      x_cell_method='sum', v_extensive=.true.)
  cs%id_umo = register_diag_field('ocean_model', 'umo', &
      diag%axesCuL, time, 'Ocean Mass X Transport', 'kg/s', &
      standard_name='ocean_mass_x_transport', y_cell_method='sum', v_extensive=.true.)
  cs%id_vmo = register_diag_field('ocean_model', 'vmo', &
      diag%axesCvL, time, 'Ocean Mass Y Transport', 'kg/s', &
      standard_name='ocean_mass_y_transport', x_cell_method='sum', v_extensive=.true.)
  cs%id_umo_2d = register_diag_field('ocean_model', 'umo_2d', &
      diag%axesCu1, time, 'Ocean Mass X Transport Vertical Sum', 'kg/s', &
      standard_name='ocean_mass_x_transport_vertical_sum', y_cell_method='sum')
  cs%id_vmo_2d = register_diag_field('ocean_model', 'vmo_2d', &
      diag%axesCv1, time, 'Ocean Mass Y Transport Vertical Sum', 'kg/s', &
      standard_name='ocean_mass_y_transport_vertical_sum', x_cell_method='sum')

end subroutine register_diags


!> Initialize diagnostics for temp/heat and saln/salt tendencies.
subroutine register_diags_ts_tendency(Time, G, CS)
  type(time_type),           intent(in)    :: Time  !< current model time
  type(ocean_grid_type),     intent(inout) :: G     !< ocean grid structure
  type(mom_control_struct),  pointer       :: CS    !< control structure set up by initialize_MOM

  type(diag_ctrl), pointer :: diag
  integer :: i, j, k
  integer :: isd, ied, jsd, jed, nz
  integer :: is, ie, js, je

  diag => cs%diag
  isd  = g%isd  ; ied  = g%ied  ; jsd  = g%jsd  ; jed  = g%jed ; nz = g%ke
  is   = g%isc  ; ie   = g%iec  ; js   = g%jsc  ; je   = g%jec


  ! heat tendencies from lateral advection
  cs%id_T_advection_xy = register_diag_field('ocean_model', 'T_advection_xy', diag%axesTL, time, &
      'Horizontal convergence of residual mean heat advective fluxes', 'W/m2',v_extensive=.true.)
  cs%id_T_advection_xy_2d = register_diag_field('ocean_model', 'T_advection_xy_2d', diag%axesT1, time,&
      'Vertical sum of horizontal convergence of residual mean heat advective fluxes', 'W/m2')
  if (cs%id_T_advection_xy > 0 .or. cs%id_T_advection_xy_2d > 0) then
    call safe_alloc_ptr(cs%T_advection_xy,isd,ied,jsd,jed,nz)
    cs%tendency_diagnostics = .true.
  endif

  ! net temperature and heat tendencies
  cs%id_T_tendency = register_diag_field('ocean_model', 'T_tendency', diag%axesTL, time, &
      'Net time tendency for temperature', 'degC/s')
  cs%id_Th_tendency = register_diag_field('ocean_model', 'Th_tendency', diag%axesTL, time,        &
      'Net time tendency for heat', 'W/m2',                                                       &
      cmor_field_name="opottemptend", cmor_units="W m-2",                                         &
      cmor_standard_name="tendency_of_sea_water_potential_temperature_expressed_as_heat_content", &
      cmor_long_name ="Tendency of Sea Water Potential Temperature Expressed as Heat Content",    &
      v_extensive=.true.)
  cs%id_Th_tendency_2d = register_diag_field('ocean_model', 'Th_tendency_2d', diag%axesT1, time,              &
      'Vertical sum of net time tendency for heat', 'W/m2',                                                   &
      cmor_field_name="opottemptend_2d", cmor_units="W m-2",                                                   &
      cmor_standard_name="tendency_of_sea_water_potential_temperature_expressed_as_heat_content_vertical_sum",&
      cmor_long_name ="Tendency of Sea Water Potential Temperature Expressed as Heat Content Vertical Sum")
  if (cs%id_T_tendency > 0) then
    cs%tendency_diagnostics = .true.
    call safe_alloc_ptr(cs%T_prev,isd,ied,jsd,jed,nz)
    do k=1,nz ; do j=js,je ; do i=is,ie
      cs%T_prev(i,j,k) = cs%tv%T(i,j,k)
    enddo ; enddo ; enddo
  endif
  if (cs%id_Th_tendency > 0 .or. cs%id_Th_tendency_2d > 0) then
    cs%tendency_diagnostics = .true.
    call safe_alloc_ptr(cs%Th_prev,isd,ied,jsd,jed,nz)
    do k=1,nz ; do j=js,je ; do i=is,ie
      cs%Th_prev(i,j,k) = cs%tv%T(i,j,k) * cs%h(i,j,k)
    enddo ; enddo ; enddo
  endif


  ! salt tendencies from lateral advection
  cs%id_S_advection_xy = register_diag_field('ocean_model', 'S_advection_xy', diag%axesTL, time, &
      'Horizontal convergence of residual mean salt advective fluxes', 'kg/(m2 * s)', v_extensive=.true.)
  cs%id_S_advection_xy_2d = register_diag_field('ocean_model', 'S_advection_xy_2d', diag%axesT1, time,&
      'Vertical sum of horizontal convergence of residual mean salt advective fluxes', 'kg/(m2 * s)')
  if (cs%id_S_advection_xy > 0 .or. cs%id_S_advection_xy_2d > 0) then
    call safe_alloc_ptr(cs%S_advection_xy,isd,ied,jsd,jed,nz)
    cs%tendency_diagnostics = .true.
  endif

  ! net salinity and salt tendencies
  cs%id_S_tendency = register_diag_field('ocean_model', 'S_tendency', diag%axesTL, time, &
      'Net time tendency for salinity', 'PPT/s')
  cs%id_Sh_tendency = register_diag_field('ocean_model', 'Sh_tendency', diag%axesTL, time,&
      'Net time tendency for salt', 'kg/(m2 * s)',                                        &
      cmor_field_name="osalttend", cmor_units="kg m-2 s-1",                               &
      cmor_standard_name="tendency_of_sea_water_salinity_expressed_as_salt_content",      &
      cmor_long_name ="Tendency of Sea Water Salinity Expressed as Salt Content",         &
      v_extensive=.true.)
  cs%id_Sh_tendency_2d = register_diag_field('ocean_model', 'Sh_tendency_2d', diag%axesT1, time, &
      'Vertical sum of net time tendency for salt', 'kg/(m2 * s)',                               &
      cmor_field_name="osalttend_2d", cmor_units="kg m-2 s-1",                                   &
      cmor_standard_name="tendency_of_sea_water_salinity_expressed_as_salt_content_vertical_sum",&
      cmor_long_name ="Tendency of Sea Water Salinity Expressed as Salt Content Vertical Sum")
  if (cs%id_S_tendency > 0) then
    cs%tendency_diagnostics = .true.
    call safe_alloc_ptr(cs%S_prev,isd,ied,jsd,jed,nz)
    do k=1,nz ; do j=js,je ; do i=is,ie
      cs%S_prev(i,j,k) = cs%tv%S(i,j,k)
    enddo ; enddo ; enddo
  endif
  if (cs%id_Sh_tendency > 0 .or. cs%id_Sh_tendency_2d > 0) then
    cs%tendency_diagnostics = .true.
    call safe_alloc_ptr(cs%Sh_prev,isd,ied,jsd,jed,nz)
    do k=1,nz ; do j=js,je ; do i=is,ie
      cs%Sh_prev(i,j,k) = cs%tv%S(i,j,k) * cs%h(i,j,k)
    enddo ; enddo ; enddo
  endif

end subroutine register_diags_ts_tendency


!> Initialize diagnostics for the variance decay of temp/salt
!! across regridding/remapping
subroutine register_diags_ts_vardec(Time, HI, GV, param_file, CS)
  type(time_type),         intent(in) :: Time     !< current model time
  type(hor_index_type),    intent(in) :: HI       !< horizontal index type
  type(verticalgrid_type), intent(in) :: GV       !< ocean vertical grid structure
  type(param_file_type),   intent(in) :: param_file !< parameter file
  type(mom_control_struct), pointer :: CS   !< control structure for MOM

  integer :: isd, ied, jsd, jed, nz
  type(vardesc) :: vd_tmp
  type(diag_ctrl), pointer :: diag

  diag => cs%diag
  isd  = hi%isd  ; ied  = hi%ied  ; jsd  = hi%jsd  ; jed  = hi%jed ; nz = gv%ke

  ! variancy decay through ALE operation
  cs%id_T_vardec = register_diag_field('ocean_model', 'T_vardec', diag%axesTL, time, &
      'ALE variance decay for temperature', 'degC2 s-1')
  if (cs%id_T_vardec > 0) then
    call safe_alloc_ptr(cs%T_squared,isd,ied,jsd,jed,nz)
    cs%T_squared(:,:,:) = 0.

    vd_tmp = var_desc(name="T2", units="degC2", longname="Squared Potential Temperature")
    call register_tracer(cs%T_squared, vd_tmp, param_file, hi, gv, cs%tracer_reg)
  endif

  cs%id_S_vardec = register_diag_field('ocean_model', 'S_vardec', diag%axesTL, time, &
      'ALE variance decay for salinity', 'PPT2 s-1')
  if (cs%id_S_vardec > 0) then
    call safe_alloc_ptr(cs%S_squared,isd,ied,jsd,jed,nz)
    cs%S_squared(:,:,:) = 0.

    vd_tmp = var_desc(name="S2", units="PPT2", longname="Squared Salinity")
    call register_tracer(cs%S_squared, vd_tmp, param_file, hi, gv, cs%tracer_reg)
  endif

end subroutine register_diags_ts_vardec

!> This subroutine sets up clock IDs for timing various subroutines.
subroutine mom_timing_init(CS)
  type(mom_control_struct), intent(in) :: CS  !< control structure set up by initialize_MOM.

 id_clock_ocean    = cpu_clock_id('Ocean', grain=clock_component)
 id_clock_dynamics = cpu_clock_id('Ocean dynamics', grain=clock_subcomponent)
 id_clock_thermo   = cpu_clock_id('Ocean thermodynamics and tracers', grain=clock_subcomponent)
 id_clock_other    = cpu_clock_id('Ocean Other', grain=clock_subcomponent)
 id_clock_tracer   = cpu_clock_id('(Ocean tracer advection)', grain=clock_module_driver)
 if (.not.cs%adiabatic) &
   id_clock_diabatic = cpu_clock_id('(Ocean diabatic driver)', grain=clock_module_driver)

 id_clock_continuity = cpu_clock_id('(Ocean continuity equation *)', grain=clock_module)
 id_clock_bbl_visc = cpu_clock_id('(Ocean set BBL viscosity)', grain=clock_module)
 id_clock_pass       = cpu_clock_id('(Ocean message passing *)', grain=clock_module)
 id_clock_mom_init   = cpu_clock_id('(Ocean MOM_initialize_state)', grain=clock_module)
 id_clock_pass_init  = cpu_clock_id('(Ocean init message passing *)', grain=clock_routine)
 if (cs%thickness_diffuse) &
   id_clock_thick_diff = cpu_clock_id('(Ocean thickness diffusion *)', grain=clock_module)
!if (CS%mixedlayer_restrat) &
   id_clock_ml_restrat = cpu_clock_id('(Ocean mixed layer restrat)', grain=clock_module)
 id_clock_diagnostics  = cpu_clock_id('(Ocean collective diagnostics)', grain=clock_module)
 id_clock_z_diag       = cpu_clock_id('(Ocean Z-space diagnostics)', grain=clock_module)
 id_clock_ale          = cpu_clock_id('(Ocean ALE)', grain=clock_module)
 if(cs%offline_tracer_mode) then
  id_clock_offline_tracer = cpu_clock_id('Ocean offline tracers', grain=clock_subcomponent)
 endif

end subroutine mom_timing_init

!> This routine posts diagnostics of the transports, including the subgridscale
!! contributions.
subroutine post_transport_diagnostics(G, GV, CS, diag, dt_trans, h, h_pre_dyn)
  type(ocean_grid_type),    intent(inout) :: G   !< ocean grid structure
  type(verticalgrid_type),  intent(in)    :: GV  !< ocean vertical grid structure
  type(mom_control_struct), intent(in)    :: CS  !< control structure
  type(diag_ctrl),          intent(inout) :: diag !< regulates diagnostic output
  real                    , intent(in)    :: dt_trans !< total time step associated with the transports, in s.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                            intent(in)    :: h   !< The updated layer thicknesses, in H
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                            intent(in)    :: h_pre_dyn !< The thickness before the transports, in H.

  real, dimension(SZIB_(G), SZJ_(G)) :: umo2d ! Diagnostics of integrated mass transport, in kg s-1
  real, dimension(SZI_(G), SZJB_(G)) :: vmo2d ! Diagnostics of integrated mass transport, in kg s-1
  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)) :: umo ! Diagnostics of layer mass transport, in kg s-1
  real, dimension(SZI_(G), SZJB_(G), SZK_(G)) :: vmo ! Diagnostics of layer mass transport, in kg s-1
  real :: H_to_kg_m2_dt   ! A conversion factor from accumulated transports to fluxes, in kg m-2 H-1 s-1.
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  call cpu_clock_begin(id_clock_z_diag)
  call calculate_z_transport(cs%uhtr, cs%vhtr, h, dt_trans, g, gv, &
                             cs%diag_to_Z_CSp)
  call cpu_clock_end(id_clock_z_diag)

  ! Post mass transports, including SGS
  ! Build the remap grids using the layer thicknesses from before the dynamics
  call diag_update_remap_grids(diag, alt_h = h_pre_dyn)

  h_to_kg_m2_dt = gv%H_to_kg_m2 / dt_trans
  if (cs%id_umo_2d > 0) then
    umo2d(:,:) = 0.0
    do k=1,nz ; do j=js,je ; do i=is-1,ie
      umo2d(i,j) = umo2d(i,j) + cs%uhtr(i,j,k) * h_to_kg_m2_dt
    enddo ; enddo ; enddo
    call post_data(cs%id_umo_2d, umo2d, diag)
  endif
  if (cs%id_umo > 0) then
    ! Convert to kg/s. Modifying the array for diagnostics is allowed here since it is set to zero immediately below
    do k=1,nz ; do j=js,je ; do i=is-1,ie
      umo(i,j,k) =  cs%uhtr(i,j,k) * h_to_kg_m2_dt
    enddo ; enddo ; enddo
    call post_data(cs%id_umo, umo, diag, alt_h = h_pre_dyn)
  endif
  if (cs%id_vmo_2d > 0) then
    vmo2d(:,:) = 0.0
    do k=1,nz ; do j=js-1,je ; do i=is,ie
      vmo2d(i,j) = vmo2d(i,j) + cs%vhtr(i,j,k) * h_to_kg_m2_dt
    enddo ; enddo ; enddo
    call post_data(cs%id_vmo_2d, vmo2d, diag)
  endif
  if (cs%id_vmo > 0) then
    ! Convert to kg/s. Modifying the array for diagnostics is allowed here since it is set to zero immediately below
    do k=1,nz ; do j=js-1,je ; do i=is,ie
      vmo(i,j,k) = cs%vhtr(i,j,k) * h_to_kg_m2_dt
    enddo ; enddo ; enddo
    call post_data(cs%id_vmo, vmo, diag, alt_h = h_pre_dyn)
  endif

  if (cs%id_uhtr > 0) call post_data(cs%id_uhtr, cs%uhtr, diag, alt_h = h_pre_dyn)
  if (cs%id_vhtr > 0) call post_data(cs%id_vhtr, cs%vhtr, diag, alt_h = h_pre_dyn)

end subroutine post_transport_diagnostics

!> Post diagnostics of temperatures and salinities, their fluxes, and tendencies.
subroutine post_ts_diagnostics(CS, G, GV, tv, diag, dt)
  type(mom_control_struct), intent(inout) :: CS  !< control structure
  type(ocean_grid_type),    intent(in)    :: G   !< ocean grid structure
  type(verticalgrid_type),  intent(in)    :: GV  !< ocean vertical grid structure
  type(thermo_var_ptrs),    intent(in)    :: tv  !< A structure pointing to various thermodynamic variables  type(forcing),                             intent(inout) :: fluxes !< pointers to forcing fields
  type(diag_ctrl),          intent(in)    :: diag !< regulates diagnostic output
  real,                     intent(in)    :: dt  !< total time step for T,S update

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: potTemp, pracSal !TEOS10 Diagnostics
  real    :: work3d(szi_(g),szj_(g),szk_(g))
  real    :: work2d(szi_(g),szj_(g))
  real    :: Idt, ppt2mks
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  if (.NOT. cs%use_conT_absS) then
    !Internal T&S variables are assumed to be potential&practical
    if (cs%id_T > 0) call post_data(cs%id_T, tv%T, diag)
    if (cs%id_S > 0) call post_data(cs%id_S, tv%S, diag)

    if (cs%id_tob > 0) call post_data(cs%id_tob, tv%T(:,:,g%ke), diag, mask=g%mask2dT)
    if (cs%id_sob > 0) call post_data(cs%id_sob, tv%S(:,:,g%ke), diag, mask=g%mask2dT)
  else
    !Internal T&S variables are assumed to be conservative&absolute
    if (cs%id_Tcon > 0) call post_data(cs%id_Tcon, tv%T, diag)
    if (cs%id_Sabs > 0) call post_data(cs%id_Sabs, tv%S, diag)
    !Using TEOS-10 function calls convert T&S diagnostics
    !from conservative temp to potential temp and
    !from absolute salinity to practical salinity
    do k=1,nz ; do j=js,je ; do i=is,ie
      pracsal(i,j,k) = gsw_sp_from_sr(tv%S(i,j,k))
      pottemp(i,j,k) = gsw_pt_from_ct(tv%S(i,j,k),tv%T(i,j,k))
    enddo; enddo ; enddo
    if (cs%id_T > 0) call post_data(cs%id_T, pottemp, diag)
    if (cs%id_S > 0) call post_data(cs%id_S, pracsal, diag)
    if (cs%id_tob > 0) call post_data(cs%id_tob, pottemp(:,:,g%ke), diag, mask=g%mask2dT)
    if (cs%id_sob > 0) call post_data(cs%id_sob, pracsal(:,:,g%ke), diag, mask=g%mask2dT)
  endif

  if (cs%id_Tadx   > 0) call post_data(cs%id_Tadx,   cs%T_adx,   diag)
  if (cs%id_Tady   > 0) call post_data(cs%id_Tady,   cs%T_ady,   diag)
  if (cs%id_Tdiffx > 0) call post_data(cs%id_Tdiffx, cs%T_diffx, diag)
  if (cs%id_Tdiffy > 0) call post_data(cs%id_Tdiffy, cs%T_diffy, diag)

  if (cs%id_Sadx   > 0) call post_data(cs%id_Sadx,   cs%S_adx,   diag)
  if (cs%id_Sady   > 0) call post_data(cs%id_Sady,   cs%S_ady,   diag)
  if (cs%id_Sdiffx > 0) call post_data(cs%id_Sdiffx, cs%S_diffx, diag)
  if (cs%id_Sdiffy > 0) call post_data(cs%id_Sdiffy, cs%S_diffy, diag)

  if (cs%id_Tadx_2d   > 0) call post_data(cs%id_Tadx_2d,   cs%T_adx_2d,   diag)
  if (cs%id_Tady_2d   > 0) call post_data(cs%id_Tady_2d,   cs%T_ady_2d,   diag)
  if (cs%id_Tdiffx_2d > 0) call post_data(cs%id_Tdiffx_2d, cs%T_diffx_2d, diag)
  if (cs%id_Tdiffy_2d > 0) call post_data(cs%id_Tdiffy_2d, cs%T_diffy_2d, diag)

  if (cs%id_Sadx_2d   > 0) call post_data(cs%id_Sadx_2d,   cs%S_adx_2d,   diag)
  if (cs%id_Sady_2d   > 0) call post_data(cs%id_Sady_2d,   cs%S_ady_2d,   diag)
  if (cs%id_Sdiffx_2d > 0) call post_data(cs%id_Sdiffx_2d, cs%S_diffx_2d, diag)
  if (cs%id_Sdiffy_2d > 0) call post_data(cs%id_Sdiffy_2d, cs%S_diffy_2d, diag)

  if(.not. cs%tendency_diagnostics) return

  idt = 0.; if (dt/=0.) idt = 1.0 / dt ! The "if" is in case the diagnostic is called for a zero length interval
  ppt2mks       = 0.001
  work3d(:,:,:) = 0.0
  work2d(:,:)   = 0.0

  ! Diagnose tendency of heat from convergence of lateral advective,
  ! fluxes, where advective transport arises from residual mean velocity.
  if (cs%id_T_advection_xy > 0 .or. cs%id_T_advection_xy_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work3d(i,j,k) = cs%T_advection_xy(i,j,k) * gv%H_to_kg_m2 * tv%C_p
    enddo ; enddo ; enddo
    if (cs%id_T_advection_xy    > 0) call post_data(cs%id_T_advection_xy, work3d, diag)
    if (cs%id_T_advection_xy_2d > 0) then
      do j=js,je ; do i=is,ie
        work2d(i,j) = 0.0
        do k=1,nz
          work2d(i,j) = work2d(i,j) + work3d(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_T_advection_xy_2d, work2d, diag)
    endif
  endif

  ! Diagnose tendency of salt from convergence of lateral advective
  ! fluxes, where advective transport arises from residual mean velocity.
  if (cs%id_S_advection_xy > 0 .or. cs%id_S_advection_xy_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work3d(i,j,k) = cs%S_advection_xy(i,j,k) * gv%H_to_kg_m2 * ppt2mks
    enddo ; enddo ; enddo
    if (cs%id_S_advection_xy    > 0) call post_data(cs%id_S_advection_xy, work3d, diag)
    if (cs%id_S_advection_xy_2d > 0) then
      do j=js,je ; do i=is,ie
        work2d(i,j) = 0.0
        do k=1,nz
          work2d(i,j) = work2d(i,j) + work3d(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_S_advection_xy_2d, work2d, diag)
    endif
  endif

  ! diagnose net tendency for temperature over a time step and update T_prev
  if (cs%id_T_tendency > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work3d(i,j,k)    = (tv%T(i,j,k) - cs%T_prev(i,j,k))*idt
      cs%T_prev(i,j,k) =  tv%T(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_T_tendency, work3d, diag)
  endif

  ! diagnose net tendency for salinity over a time step and update S_prev
  if (cs%id_S_tendency > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work3d(i,j,k)    = (tv%S(i,j,k) - cs%S_prev(i,j,k))*idt
      cs%S_prev(i,j,k) =  tv%S(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_S_tendency, work3d, diag)
  endif

  ! diagnose net tendency for heat content of a grid cell over a time step and update Th_prev
  if (cs%id_Th_tendency > 0 .or. cs%id_Th_tendency_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work3d(i,j,k)     = (tv%T(i,j,k)*cs%h(i,j,k) - cs%Th_prev(i,j,k)) * idt * gv%H_to_kg_m2 * tv%C_p
      cs%Th_prev(i,j,k) =  tv%T(i,j,k)*cs%h(i,j,k)
    enddo ; enddo ; enddo
    if (cs%id_Th_tendency    > 0) call post_data(cs%id_Th_tendency, work3d, diag)
    if (cs%id_Th_tendency_2d > 0) then
      do j=js,je ; do i=is,ie
        work2d(i,j) = 0.0
        do k=1,nz
          work2d(i,j) = work2d(i,j) + work3d(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_Th_tendency_2d, work2d, diag)
    endif
  endif

  ! diagnose net tendency for salt content of a grid cell over a time step and update Sh_prev
  if (cs%id_Sh_tendency > 0 .or. cs%id_Sh_tendency_2d > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work3d(i,j,k)     = (tv%S(i,j,k)*cs%h(i,j,k) - cs%Sh_prev(i,j,k)) * idt * gv%H_to_kg_m2 * ppt2mks
      cs%Sh_prev(i,j,k) =  tv%S(i,j,k)*cs%h(i,j,k)
    enddo ; enddo ; enddo
    if (cs%id_Sh_tendency    > 0) call post_data(cs%id_Sh_tendency, work3d, diag)
    if (cs%id_Sh_tendency_2d > 0) then
      do j=js,je ; do i=is,ie
        work2d(i,j) = 0.0
        do k=1,nz
          work2d(i,j) = work2d(i,j) + work3d(i,j,k)
        enddo
      enddo ; enddo
      call post_data(cs%id_Sh_tendency_2d, work2d, diag)
    endif
  endif

end subroutine post_ts_diagnostics

!> Calculate and post variance decay diagnostics for temp/salt
subroutine post_diags_ts_vardec(G, CS, dt)
  type(ocean_grid_type),    intent(in) :: G    !< ocean grid structure
  type(mom_control_struct), intent(in) :: CS   !< control structure
  real, intent(in) :: dt                       !< total time step

  real :: work(szi_(g),szj_(g),szk_(g))
  real :: Idt
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  idt = 0.; if (dt/=0.) idt = 1.0 / dt ! The "if" is in case the diagnostic is called for a zero length interval

  if (cs%id_T_vardec > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work(i,j,k) = (cs%T_squared(i,j,k) - cs%tv%T(i,j,k)**2) * idt
    enddo ; enddo ; enddo
    call post_data(cs%id_T_vardec, work, cs%diag)
  endif

  if (cs%id_S_vardec > 0) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      work(i,j,k) = (cs%S_squared(i,j,k) - cs%tv%S(i,j,k)**2) * idt
    enddo ; enddo ; enddo
    call post_data(cs%id_S_vardec, work, cs%diag)
  endif
end subroutine post_diags_ts_vardec

!> This routine posts diagnostics of various integrated quantities.
subroutine post_integrated_diagnostics(CS, G, GV, diag, dt_int, tv, fluxes)
  type(mom_control_struct), intent(in)    :: CS  !< control structure
  type(ocean_grid_type),    intent(in)    :: G   !< ocean grid structure
  type(verticalgrid_type),  intent(in)    :: GV  !< ocean vertical grid structure
  type(diag_ctrl),          intent(in)    :: diag  !< regulates diagnostic output
  real,                     intent(in)    :: dt_int  !< total time step associated with these diagnostics, in s.
  type(thermo_var_ptrs),    intent(in)    :: tv  !< A structure pointing to various thermodynamic variables  type(forcing),                             intent(inout) :: fluxes !< pointers to forcing fields
  type(forcing),            intent(in)    :: fluxes  !< pointers to forcing fields

  real, allocatable, dimension(:,:) :: &
    tmp,              & ! temporary 2d field
    zos,              & ! dynamic sea lev (zero area mean) from inverse-barometer adjusted ssh (meter)
    zossq,            & ! square of zos (m^2)
    sfc_speed,        & ! sea surface speed at h-points (m/s)
    frazil_ave,       & ! average frazil heat flux required to keep temp above freezing (W/m2)
    salt_deficit_ave, & ! average salt flux required to keep salinity above 0.01ppt (gSalt m-2 s-1)
    Heat_PmE_ave,     & ! average effective heat flux into the ocean due to
                        ! the exchange of water with other components, times the
                        ! heat capacity of water, in W m-2.
    intern_heat_ave     ! avg heat flux into ocean from geothermal or
                        ! other internal heat sources (W/m2)
  real :: I_time_int    ! The inverse of the time interval in s-1.
  real :: zos_area_mean, volo, ssh_ga
  integer :: i, j, k, is, ie, js, je, nz! , Isq, Ieq, Jsq, Jeq

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  ! area mean SSH
  if (cs%id_ssh_ga > 0) then
    ssh_ga = global_area_mean(cs%ave_ssh, g)
    call post_data(cs%id_ssh_ga, ssh_ga, diag)
  endif

  i_time_int = 1.0 / dt_int
  if (cs%id_ssh > 0) &
    call post_data(cs%id_ssh, cs%ave_ssh, diag, mask=g%mask2dT)

  ! post the dynamic sea level, zos, and zossq.
  ! zos is ave_ssh with sea ice inverse barometer removed,
  ! and with zero global area mean.
  if(cs%id_zos > 0 .or. cs%id_zossq > 0) then
     allocate(zos(g%isd:g%ied,g%jsd:g%jed))
     zos(:,:) = 0.0
     do j=js,je ; do i=is,ie
       zos(i,j) = cs%ave_ssh(i,j)
     enddo ; enddo
     if (ASSOCIATED(fluxes%p_surf)) then
       do j=js,je ; do i=is,ie
         zos(i,j) = zos(i,j) + g%mask2dT(i,j)*fluxes%p_surf(i,j) / &
                              (gv%Rho0 * gv%g_Earth)
       enddo ; enddo
     endif
     zos_area_mean = global_area_mean(zos, g)
     do j=js,je ; do i=is,ie
       zos(i,j) = zos(i,j) - g%mask2dT(i,j)*zos_area_mean
     enddo ; enddo
     if(cs%id_zos > 0) then
       call post_data(cs%id_zos, zos, diag, mask=g%mask2dT)
     endif
     if(cs%id_zossq > 0) then
       allocate(zossq(g%isd:g%ied,g%jsd:g%jed))
       zossq(:,:) = 0.0
       do j=js,je ; do i=is,ie
         zossq(i,j) = zos(i,j)*zos(i,j)
       enddo ; enddo
       call post_data(cs%id_zossq, zossq, diag, mask=g%mask2dT)
       deallocate(zossq)
     endif
     deallocate(zos)
  endif

  ! post total volume of the liquid ocean
  if(cs%id_volo > 0) then
    allocate(tmp(g%isd:g%ied,g%jsd:g%jed))
    do j=js,je ; do i=is,ie
      tmp(i,j) = g%mask2dT(i,j)*(cs%ave_ssh(i,j) + g%bathyT(i,j))
    enddo ; enddo
    volo = global_area_integral(tmp, g)
    call post_data(cs%id_volo, volo, diag)
    deallocate(tmp)
  endif

  ! post frazil
  if (ASSOCIATED(tv%frazil) .and. (cs%id_fraz > 0)) then
    allocate(frazil_ave(g%isd:g%ied,g%jsd:g%jed))
    do j=js,je ; do i=is,ie
      frazil_ave(i,j) = tv%frazil(i,j) * i_time_int
    enddo ; enddo
    call post_data(cs%id_fraz, frazil_ave, diag, mask=g%mask2dT)
    deallocate(frazil_ave)
  endif

  ! post the salt deficit
  if (ASSOCIATED(tv%salt_deficit) .and. (cs%id_salt_deficit > 0)) then
    allocate(salt_deficit_ave(g%isd:g%ied,g%jsd:g%jed))
    do j=js,je ; do i=is,ie
      salt_deficit_ave(i,j) = tv%salt_deficit(i,j) * i_time_int
    enddo ; enddo
    call post_data(cs%id_salt_deficit, salt_deficit_ave, diag, mask=g%mask2dT)
    deallocate(salt_deficit_ave)
  endif

  ! post temperature of P-E+R
  if (ASSOCIATED(tv%TempxPmE) .and. (cs%id_Heat_PmE > 0)) then
    allocate(heat_pme_ave(g%isd:g%ied,g%jsd:g%jed))
    do j=js,je ; do i=is,ie
      heat_pme_ave(i,j) = tv%TempxPmE(i,j) * (tv%C_p * i_time_int)
    enddo ; enddo
    call post_data(cs%id_Heat_PmE, heat_pme_ave, diag, mask=g%mask2dT)
    deallocate(heat_pme_ave)
  endif

  ! post geothermal heating or internal heat source/sinks
  if (ASSOCIATED(tv%internal_heat) .and. (cs%id_intern_heat > 0)) then
    allocate(intern_heat_ave(g%isd:g%ied,g%jsd:g%jed))
    do j=js,je ; do i=is,ie
      intern_heat_ave(i,j) = tv%internal_heat(i,j) * (tv%C_p * i_time_int)
    enddo ; enddo
    call post_data(cs%id_intern_heat, intern_heat_ave, diag, mask=g%mask2dT)
    deallocate(intern_heat_ave)
  endif

end subroutine post_integrated_diagnostics

!> This routine posts diagnostics of various ocean surface quantities.
subroutine post_surface_diagnostics(CS, G, diag, state)
  type(mom_control_struct), intent(in)    :: CS  !< control structure
  type(ocean_grid_type),    intent(in)    :: G   !< ocean grid structure
  type(diag_ctrl),          intent(in)    :: diag  !< regulates diagnostic output
  type(surface),            intent(in)    :: state !< ocean surface state

  real, dimension(SZI_(G),SZJ_(G)) :: &
    potTemp, &  ! TEOS10 potential temperature (deg C)
    pracSal, &  ! TEOS10 practical salinity
    SST_sq, &   ! Surface temperature squared, in degC^2
    SSS_sq, &   ! Surface salinity squared, in salnity units^2
    sfc_speed   ! sea surface speed at h-points (m/s)

  integer :: i, j, is, ie, js, je

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  if (.NOT.cs%use_conT_absS) then
    !Internal T&S variables are assumed to be potential&practical
    if (cs%id_sst > 0) call post_data(cs%id_sst, state%SST, diag, mask=g%mask2dT)
    if (cs%id_sss > 0) call post_data(cs%id_sss, state%SSS, diag, mask=g%mask2dT)
  else
    !Internal T&S variables are assumed to be conservative&absolute
    if (cs%id_sstcon > 0) call post_data(cs%id_sstcon, state%SST, diag, mask=g%mask2dT)
    if (cs%id_sssabs > 0) call post_data(cs%id_sssabs, state%SSS, diag, mask=g%mask2dT)
    !Using TEOS-10 function calls convert T&S diagnostics
    !from conservative temp to potential temp and
    !from absolute salinity to practical salinity
    do j=js,je ; do i=is,ie
      pracsal(i,j) = gsw_sp_from_sr(state%SSS(i,j))
      pottemp(i,j) = gsw_pt_from_ct(state%SSS(i,j),state%SST(i,j))
    enddo ; enddo
    if (cs%id_sst > 0) call post_data(cs%id_sst, pottemp, diag, mask=g%mask2dT)
    if (cs%id_sss > 0) call post_data(cs%id_sss, pracsal, diag, mask=g%mask2dT)
  endif

  if (cs%id_sst_sq > 0) then
    do j=js,je ; do i=is,ie
      sst_sq(i,j) = state%SST(i,j)*state%SST(i,j)
    enddo ; enddo
    call post_data(cs%id_sst_sq, sst_sq, diag, mask=g%mask2dT)
  endif
  if (cs%id_sss_sq > 0) then
    do j=js,je ; do i=is,ie
      sss_sq(i,j) = state%SSS(i,j)*state%SSS(i,j)
    enddo ; enddo
    call post_data(cs%id_sss_sq, sss_sq, diag, mask=g%mask2dT)
  endif

  if (cs%id_ssu > 0) &
    call post_data(cs%id_ssu, state%u, diag, mask=g%mask2dCu)
  if (cs%id_ssv > 0) &
    call post_data(cs%id_ssv, state%v, diag, mask=g%mask2dCv)

  if (cs%id_speed > 0) then
    do j=js,je ; do i=is,ie
      sfc_speed(i,j) = sqrt(0.5*(state%u(i-1,j)**2 + state%u(i,j)**2) + &
                            0.5*(state%v(i,j-1)**2 + state%v(i,j)**2))
    enddo ; enddo
    call post_data(cs%id_speed, sfc_speed, diag, mask=g%mask2dT)
  endif

end subroutine post_surface_diagnostics

!> Offers the static fields in the ocean grid type
!! for output via the diag_manager.
subroutine write_static_fields(G, diag)
  type(ocean_grid_type),   intent(in)    :: G      !< ocean grid structure
  type(diag_ctrl), target, intent(inout) :: diag   !< regulates diagnostic output

  ! The out_X arrays are needed because some of the elements of the grid
  ! type may be reduced rank macros.
  real    :: out_h(szi_(g),szj_(g))
  integer :: id, i, j, is, ie, js, je
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec

  out_h(:,:) = 0.0

  id = register_static_field('ocean_model', 'geolat', diag%axesT1, &
        'Latitude of tracer (T) points', 'degrees_N')
  if (id > 0) call post_data(id, g%geoLatT, diag, .true.)

  id = register_static_field('ocean_model', 'geolon', diag%axesT1, &
        'Longitude of tracer (T) points', 'degrees_E')
  if (id > 0) call post_data(id, g%geoLonT, diag, .true.)

  id = register_static_field('ocean_model', 'geolat_c', diag%axesB1, &
        'Latitude of corner (Bu) points', 'degrees_N', interp_method='none')
  if (id > 0) call post_data(id, g%geoLatBu, diag, .true.)

  id = register_static_field('ocean_model', 'geolon_c', diag%axesB1, &
        'Longitude of corner (Bu) points', 'degrees_E', interp_method='none')
  if (id > 0) call post_data(id, g%geoLonBu, diag, .true.)

  id = register_static_field('ocean_model', 'geolat_v', diag%axesCv1, &
        'Latitude of meridional velocity (Cv) points', 'degrees_N', interp_method='none')
  if (id > 0) call post_data(id, g%geoLatCv, diag, .true.)

  id = register_static_field('ocean_model', 'geolon_v', diag%axesCv1, &
        'Longitude of meridional velocity (Cv) points', 'degrees_E', interp_method='none')
  if (id > 0) call post_data(id, g%geoLonCv, diag, .true.)

  id = register_static_field('ocean_model', 'geolat_u', diag%axesCu1, &
        'Latitude of zonal velocity (Cu) points', 'degrees_N', interp_method='none')
  if (id > 0) call post_data(id, g%geoLatCu, diag, .true.)

  id = register_static_field('ocean_model', 'geolon_u', diag%axesCu1, &
        'Longitude of zonal velocity (Cu) points', 'degrees_E', interp_method='none')
  if (id > 0) call post_data(id, g%geoLonCu, diag, .true.)

  id = register_static_field('ocean_model', 'area_t', diag%axesT1,   &
        'Surface area of tracer (T) cells', 'm2',                    &
        cmor_field_name='areacello', cmor_standard_name='cell_area', &
        cmor_units='m2', cmor_long_name='Ocean Grid-Cell Area')
  if (id > 0) then
    do j=js,je ; do i=is,ie ; out_h(i,j) = g%areaT(i,j) ; enddo ; enddo
    call post_data(id, out_h, diag, .true.)
    call diag_register_area_ids(diag, id_area_t=id)
  endif

  id = register_static_field('ocean_model', 'depth_ocean', diag%axesT1,  &
        'Depth of the ocean at tracer points', 'm',                      &
        standard_name='sea_floor_depth_below_geoid',                     &
        cmor_field_name='deptho', cmor_long_name='Sea Floor Depth',      &
        cmor_units='m', cmor_standard_name='sea_floor_depth_below_geoid',&
        area=diag%axesT1%id_area)
  if (id > 0) call post_data(id, g%bathyT, diag, .true., mask=g%mask2dT)

  id = register_static_field('ocean_model', 'wet', diag%axesT1, &
        '0 if land, 1 if ocean at tracer points', 'none', area=diag%axesT1%id_area)
  if (id > 0) call post_data(id, g%mask2dT, diag, .true.)

  id = register_static_field('ocean_model', 'wet_c', diag%axesB1, &
        '0 if land, 1 if ocean at corner (Bu) points', 'none', interp_method='none')
  if (id > 0) call post_data(id, g%mask2dBu, diag, .true.)

  id = register_static_field('ocean_model', 'wet_u', diag%axesCu1, &
        '0 if land, 1 if ocean at zonal velocity (Cu) points', 'none', interp_method='none')
  if (id > 0) call post_data(id, g%mask2dCu, diag, .true.)

  id = register_static_field('ocean_model', 'wet_v', diag%axesCv1, &
        '0 if land, 1 if ocean at meridional velocity (Cv) points', 'none', interp_method='none')
  if (id > 0) call post_data(id, g%mask2dCv, diag, .true.)

  id = register_static_field('ocean_model', 'Coriolis', diag%axesB1, &
        'Coriolis parameter at corner (Bu) points', 's-1', interp_method='none')
  if (id > 0) call post_data(id, g%CoriolisBu, diag, .true.)

  id = register_static_field('ocean_model', 'dxt', diag%axesT1, &
        'Delta(x) at thickness/tracer points (meter)', 'm', interp_method='none')
  if (id > 0) call post_data(id, g%dxt, diag, .true.)

  id = register_static_field('ocean_model', 'dyt', diag%axesT1, &
        'Delta(y) at thickness/tracer points (meter)', 'm', interp_method='none')
  if (id > 0) call post_data(id, g%dyt, diag, .true.)

  id = register_static_field('ocean_model', 'dxCu', diag%axesCu1, &
        'Delta(x) at u points (meter)', 'm', interp_method='none')
  if (id > 0) call post_data(id, g%dxCu, diag, .true.)

  id = register_static_field('ocean_model', 'dyCu', diag%axesCu1, &
        'Delta(y) at u points (meter)', 'm', interp_method='none')
  if (id > 0) call post_data(id, g%dyCu, diag, .true.)

  id = register_static_field('ocean_model', 'dxCv', diag%axesCv1, &
        'Delta(x) at v points (meter)', 'm', interp_method='none')
  if (id > 0) call post_data(id, g%dxCv, diag, .true.)

  id = register_static_field('ocean_model', 'dyCv', diag%axesCv1, &
        'Delta(y) at v points (meter)', 'm', interp_method='none')
  if (id > 0) call post_data(id, g%dyCv, diag, .true.)

end subroutine write_static_fields


!> Set the fields that are needed for bitwise identical restarting
!! the time stepping scheme.  In addition to those specified here
!! directly, there may be fields related to the forcing or to the
!! barotropic solver that are needed; these are specified in sub-
!! routines that are called from this one.
!!
!! This routine should be altered if there are any changes to the
!! time stepping scheme.  The CHECK_RESTART facility may be used to
!! confirm that all needed restart fields have been included.
subroutine set_restart_fields(GV, param_file, CS)
  type(verticalgrid_type),  intent(inout) :: GV         !< ocean vertical grid structure
  type(param_file_type),    intent(in) :: param_file    !< opened file for parsing to get parameters
  type(mom_control_struct), intent(in) :: CS            !< control structure set up by inialize_MOM
  ! Local variables
  logical :: use_ice_shelf ! Needed to determine whether to add CS%Hml to restarts
  type(vardesc) :: vd
  character(len=48) :: thickness_units, flux_units

  call get_param(param_file, '', "ICE_SHELF", use_ice_shelf, default=.false., do_not_log=.true.)

  thickness_units = get_thickness_units(gv)
  flux_units = get_flux_units(gv)

  if (cs%use_temperature) then
    vd = var_desc("Temp","degC","Potential Temperature")
    call register_restart_field(cs%tv%T, vd, .true., cs%restart_CSp)

    vd = var_desc("Salt","PPT","Salinity")
    call register_restart_field(cs%tv%S, vd, .true., cs%restart_CSp)
  endif

  vd = var_desc("h",thickness_units,"Layer Thickness")
  call register_restart_field(cs%h, vd, .true., cs%restart_CSp)

  vd = var_desc("u","meter second-1","Zonal velocity",'u','L')
  call register_restart_field(cs%u, vd, .true., cs%restart_CSp)

  vd = var_desc("v","meter second-1","Meridional velocity",'v','L')
  call register_restart_field(cs%v, vd, .true., cs%restart_CSp)

  if (cs%use_frazil) then
    vd = var_desc("frazil","J m-2","Frazil heat flux into ocean",'h','1')
    call register_restart_field(cs%tv%frazil, vd, .false., cs%restart_CSp)
  endif

  if (cs%interp_p_surf) then
    vd = var_desc("p_surf_prev","Pa","Previous ocean surface pressure",'h','1')
    call register_restart_field(cs%p_surf_prev, vd, .false., cs%restart_CSp)
  endif

  vd = var_desc("ave_ssh","meter","Time average sea surface height",'h','1')
  call register_restart_field(cs%ave_ssh, vd, .false., cs%restart_CSp)

  ! hML is needed when using the ice shelf module
  if (use_ice_shelf .and. associated(cs%Hml)) then
     vd = var_desc("hML","meter","Mixed layer thickness",'h','1')
     call register_restart_field(cs%Hml, vd, .false., cs%restart_CSp)
  endif

end subroutine set_restart_fields

!> This subroutine applies a correction to the sea surface height to compensate
!! for the atmospheric pressure (the inverse barometer).
subroutine adjust_ssh_for_p_atm(CS, G, GV, ssh, p_atm)
  type(mom_control_struct),          intent(in)    :: CS     !< control structure
  type(ocean_grid_type),             intent(in)    :: G      !< ocean grid structure
  type(verticalgrid_type),           intent(in)    :: GV     !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G)),  intent(inout) :: ssh    !< time mean surface height (m)
  real, dimension(:,:),    optional, pointer       :: p_atm  !< atmospheric pressure (Pascal)

  real :: Rho_conv    ! The density used to convert surface pressure to
                      ! a corrected effective SSH, in kg m-3.
  real :: IgR0        ! The SSH conversion factor from Pa to m.
  integer :: i, j, is, ie, js, je

  is  = g%isc ; ie  = g%iec ; js  = g%jsc ; je  = g%jec
  if (ASSOCIATED(p_atm)) then
    ! Correct the output sea surface height for the contribution from the
    ! atmospheric pressure
    do j=js,je ; do i=is,ie
      if ((ASSOCIATED(cs%tv%eqn_of_state)) .and. (cs%calc_rho_for_sea_lev)) then
        call calculate_density(cs%tv%T(i,j,1), cs%tv%S(i,j,1), p_atm(i,j)/2.0, &
                               rho_conv, cs%tv%eqn_of_state)
      else
        rho_conv=gv%Rho0
      endif
      igr0 = 1.0 / (rho_conv * gv%g_Earth)
      ssh(i,j) = ssh(i,j) + p_atm(i,j) * igr0
    enddo ; enddo
  endif

end subroutine adjust_ssh_for_p_atm

!> This subroutine sets the surface (return) properties of the ocean
!! model by setting the appropriate fields in state.  Unused fields
!! are set to NULL or are unallocated.
subroutine calculate_surface_state(state, u, v, h, ssh, G, GV, CS)
  type(ocean_grid_type),                     intent(inout) :: G      !< ocean grid structure
  type(verticalgrid_type),                   intent(inout) :: GV     !< ocean vertical grid structure
  type(surface),                             intent(inout) :: state  !< ocean surface state
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u      !< zonal velocity (m/s)
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v      !< meridional velocity (m/s)
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h      !< layer thickness (m or kg/m2)
  real, dimension(SZI_(G),SZJ_(G)),          intent(in)    :: ssh    !< time mean surface height (m)
  type(mom_control_struct),                  intent(inout) :: CS     !< control structure

  ! local
  real :: depth(szi_(g))              ! distance from the surface (meter)
  real :: depth_ml                    ! depth over which to average to
                                      ! determine mixed layer properties (meter)
  real :: dh                          ! thickness of a layer within mixed layer (meter)
  real :: mass                        ! mass per unit area of a layer (kg/m2)

  real :: hu, hv
  integer :: i, j, k, is, ie, js, je, nz, numberOfErrors
  integer :: isd, ied, jsd, jed
  integer :: iscB, iecB, jscB, jecB, isdB, iedB, jsdB, jedB
  logical :: localError
  character(240) :: msg

  call calltree_enter("calculate_surface_state(), MOM.F90")
  is  = g%isc ; ie  = g%iec ; js  = g%jsc ; je  = g%jec ; nz = g%ke
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed
  iscb = g%iscB ; iecb = g%iecB; jscb = g%jscB ; jecb = g%jecB
  isdb = g%isdB ; iedb = g%iedB; jsdb = g%jsdB ; jedb = g%jedB

  if (.not.state%arrays_allocated) then
    if (cs%use_temperature) then
      allocate(state%SST(isd:ied,jsd:jed)) ; state%SST(:,:) = 0.0
      allocate(state%SSS(isd:ied,jsd:jed)) ; state%SSS(:,:) = 0.0
    else
      allocate(state%sfc_density(isd:ied,jsd:jed)) ; state%sfc_density(:,:) = 0.0
    endif
    allocate(state%sea_lev(isd:ied,jsd:jed)) ; state%sea_lev(:,:) = 0.0
    allocate(state%Hml(isd:ied,jsd:jed)) ; state%Hml(:,:) = 0.0
    allocate(state%u(isdb:iedb,jsd:jed)) ; state%u(:,:) = 0.0
    allocate(state%v(isd:ied,jsdb:jedb)) ; state%v(:,:) = 0.0

  ! Allocate structures for ocean_mass, ocean_heat, and ocean_salt.  This could
  ! be wrapped in a run-time flag to disable it for economy, since the 3-d
  ! sums are not negligible.
    allocate(state%ocean_mass(isd:ied,jsd:jed)) ; state%ocean_mass(:,:) = 0.0
    if (cs%use_temperature) then
      allocate(state%ocean_heat(isd:ied,jsd:jed)) ; state%ocean_heat(:,:) = 0.0
      allocate(state%ocean_salt(isd:ied,jsd:jed)) ; state%ocean_salt(:,:) = 0.0
    endif
    allocate(state%salt_deficit(isd:ied,jsd:jed)) ; state%salt_deficit(:,:) = 0.0

    state%arrays_allocated = .true.
  endif
  state%frazil => cs%tv%frazil
  state%TempxPmE => cs%tv%TempxPmE
  state%internal_heat => cs%tv%internal_heat
  if (associated(cs%visc%taux_shelf)) state%taux_shelf => cs%visc%taux_shelf
  if (associated(cs%visc%tauy_shelf)) state%tauy_shelf => cs%visc%tauy_shelf

  do j=js,je ; do i=is,ie
    state%sea_lev(i,j) = ssh(i,j)
  enddo ; enddo

  if (cs%bulkmixedlayer) then
    if (cs%use_temperature) then ; do j=js,je ; do i=is,ie
      state%SST(i,j) = cs%tv%T(i,j,1)
      state%SSS(i,j) = cs%tv%S(i,j,1)
    enddo ; enddo ; endif
    do j=js,je ; do i=iscb,iecb
      state%u(i,j) = u(i,j,1)
    enddo ; enddo
    do j=jscb,jecb ; do i=is,ie
      state%v(i,j) = v(i,j,1)
    enddo ; enddo

    if (associated(cs%Hml)) then ; do j=js,je ; do i=is,ie
      state%Hml(i,j) = cs%Hml(i,j)
    enddo ; enddo ; endif
  else

    depth_ml = cs%Hmix
  !   Determine the mean tracer properties of the uppermost depth_ml fluid.
    !$OMP parallel do default(shared) private(depth,dh)
    do j=js,je
      do i=is,ie
        depth(i) = 0.0
        if (cs%use_temperature) then
          state%SST(i,j) = 0.0 ; state%SSS(i,j) = 0.0
        else
          state%sfc_density(i,j) = 0.0
        endif
      enddo

      do k=1,nz ; do i=is,ie
        if (depth(i) + h(i,j,k)*gv%H_to_m < depth_ml) then
          dh = h(i,j,k)*gv%H_to_m
        elseif (depth(i) < depth_ml) then
          dh = depth_ml - depth(i)
        else
          dh = 0.0
        endif
        if (cs%use_temperature) then
          state%SST(i,j) = state%SST(i,j) + dh * cs%tv%T(i,j,k)
          state%SSS(i,j) = state%SSS(i,j) + dh * cs%tv%S(i,j,k)
        else
          state%sfc_density(i,j) = state%sfc_density(i,j) + dh * gv%Rlay(k)
        endif
        depth(i) = depth(i) + dh
      enddo ; enddo
  ! Calculate the average properties of the mixed layer depth.
      do i=is,ie
        if (depth(i) < gv%H_subroundoff*gv%H_to_m) &
            depth(i) = gv%H_subroundoff*gv%H_to_m
        if (cs%use_temperature) then
          state%SST(i,j) = state%SST(i,j) / depth(i)
          state%SSS(i,j) = state%SSS(i,j) / depth(i)
        else
          state%sfc_density(i,j) = state%sfc_density(i,j) / depth(i)
        endif
        state%Hml(i,j) = depth(i)
      enddo
    enddo ! end of j loop

!   Determine the mean velocities in the uppermost depth_ml fluid.
    if (cs%Hmix_UV>0.) then
      depth_ml = cs%Hmix_UV
      !$OMP parallel do default(shared) private(depth,dh,hv)
      do j=jscb,jecb
        do i=is,ie
          depth(i) = 0.0
          state%v(i,j) = 0.0
        enddo
        do k=1,nz ; do i=is,ie
          hv = 0.5 * (h(i,j,k) + h(i,j+1,k)) * gv%H_to_m
          if (depth(i) + hv < depth_ml) then
            dh = hv
          elseif (depth(i) < depth_ml) then
            dh = depth_ml - depth(i)
          else
            dh = 0.0
          endif
          state%v(i,j) = state%v(i,j) + dh * v(i,j,k)
          depth(i) = depth(i) + dh
        enddo ; enddo
        ! Calculate the average properties of the mixed layer depth.
        do i=is,ie
          if (depth(i) < gv%H_subroundoff*gv%H_to_m) &
              depth(i) = gv%H_subroundoff*gv%H_to_m
          state%v(i,j) = state%v(i,j) / depth(i)
        enddo
      enddo ! end of j loop

      !$OMP parallel do default(shared) private(depth,dh,hu)
      do j=js,je
        do i=iscb,iecb
          depth(i) = 0.0
          state%u(i,j) = 0.0
        enddo
        do k=1,nz ; do i=iscb,iecb
          hu = 0.5 * (h(i,j,k) + h(i+1,j,k)) * gv%H_to_m
          if (depth(i) + hu < depth_ml) then
            dh = hu
          elseif (depth(i) < depth_ml) then
            dh = depth_ml - depth(i)
          else
            dh = 0.0
          endif
          state%u(i,j) = state%u(i,j) + dh * u(i,j,k)
          depth(i) = depth(i) + dh
        enddo ; enddo
        ! Calculate the average properties of the mixed layer depth.
        do i=iscb,iecb
          if (depth(i) < gv%H_subroundoff*gv%H_to_m) &
              depth(i) = gv%H_subroundoff*gv%H_to_m
          state%u(i,j) = state%u(i,j) / depth(i)
        enddo
      enddo ! end of j loop
    else ! Hmix_UV<=0.
      do j=js,je ; do i=iscb,iecb
        state%u(i,j) = u(i,j,1)
      enddo ; enddo
      do j=jscb,jecb ; do i=is,ie
        state%v(i,j) = v(i,j,1)
      enddo ; enddo
    endif
  endif                                              ! end BULKMIXEDLAYER

  if (allocated(state%salt_deficit) .and. associated(cs%tv%salt_deficit)) then
    !$OMP parallel do default(shared)
    do j=js,je ; do i=is,ie
      ! Convert from gSalt to kgSalt
      state%salt_deficit(i,j) = 1000.0 * cs%tv%salt_deficit(i,j)
    enddo ; enddo
  endif

  if (allocated(state%ocean_mass) .and. allocated(state%ocean_heat) .and. &
      allocated(state%ocean_salt)) then
    !$OMP parallel do default(shared)
    do j=js,je ; do i=is,ie
      state%ocean_mass(i,j) = 0.0
      state%ocean_heat(i,j) = 0.0 ; state%ocean_salt(i,j) = 0.0
    enddo ; enddo
    !$OMP parallel do default(shared) private(mass)
    do j=js,je ; do k=1,nz; do i=is,ie
      mass = gv%H_to_kg_m2*h(i,j,k)
      state%ocean_mass(i,j) = state%ocean_mass(i,j) + mass
      state%ocean_heat(i,j) = state%ocean_heat(i,j) + mass*cs%tv%T(i,j,k)
      state%ocean_salt(i,j) = state%ocean_salt(i,j) + &
                              mass * (1.0e-3*cs%tv%S(i,j,k))
    enddo ; enddo ; enddo
  else
    if (allocated(state%ocean_mass)) then
      !$OMP parallel do default(shared)
      do j=js,je ; do i=is,ie ; state%ocean_mass(i,j) = 0.0 ; enddo ; enddo
      !$OMP parallel do default(shared)
      do j=js,je ; do k=1,nz ; do i=is,ie
        state%ocean_mass(i,j) = state%ocean_mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
      enddo ; enddo ; enddo
    endif
    if (allocated(state%ocean_heat)) then
      !$OMP parallel do default(shared)
      do j=js,je ; do i=is,ie ; state%ocean_heat(i,j) = 0.0 ; enddo ; enddo
      !$OMP parallel do default(shared) private(mass)
      do j=js,je ; do k=1,nz ; do i=is,ie
        mass = gv%H_to_kg_m2*h(i,j,k)
        state%ocean_heat(i,j) = state%ocean_heat(i,j) + mass*cs%tv%T(i,j,k)
      enddo ; enddo ; enddo
    endif
    if (allocated(state%ocean_salt)) then
      !$OMP parallel do default(shared)
      do j=js,je ; do i=is,ie ; state%ocean_salt(i,j) = 0.0 ; enddo ; enddo
      !$OMP parallel do default(shared) private(mass)
      do j=js,je ; do k=1,nz ; do i=is,ie
        mass = gv%H_to_kg_m2*h(i,j,k)
        state%ocean_salt(i,j) = state%ocean_salt(i,j) + &
                                mass * (1.0e-3*cs%tv%S(i,j,k))
      enddo ; enddo ; enddo
    endif
  endif

  if (associated(cs%tracer_flow_CSp)) then
    if (.not.associated(state%tr_fields)) allocate(state%tr_fields)
    call call_tracer_surface_state(state, h, g, cs%tracer_flow_CSp)
  endif

  if (cs%check_bad_surface_vals) then
    numberoferrors=0 ! count number of errors
    do j=js,je; do i=is,ie
      if (g%mask2dT(i,j)>0.) then
        localerror = state%sea_lev(i,j)<=-g%bathyT(i,j)       &
                .or. state%sea_lev(i,j)>= cs%bad_val_ssh_max  &
                .or. state%sea_lev(i,j)<=-cs%bad_val_ssh_max  &
                .or. state%sea_lev(i,j)+g%bathyT(i,j) < cs%bad_val_column_thickness
        if (cs%use_temperature) localerror = localerror &
                .or. state%SSS(i,j)<0.                        &
                .or. state%SSS(i,j)>=cs%bad_val_sss_max       &
                .or. state%SST(i,j)< cs%bad_val_sst_min       &
                .or. state%SST(i,j)>=cs%bad_val_sst_max
        if (localerror) then
          numberoferrors=numberoferrors+1
          if (numberoferrors<9) then ! Only report details for the first few errors
            if (cs%use_temperature) then
              write(msg(1:240),'(2(a,i4,x),2(a,f8.3,x),8(a,es11.4,x))') &
                'Extreme surface state detected: i=',i,'j=',j, &
                'x=',g%geoLonT(i,j),'y=',g%geoLatT(i,j), &
                'D=',g%bathyT(i,j),                      &
                'SSH=',state%sea_lev(i,j),               &
                'SST=',state%SST(i,j),                   &
                'SSS=',state%SSS(i,j),                   &
                'U-=',state%u(i-1,j),                    &
                'U+=',state%u(i,j),                      &
                'V-=',state%v(i,j-1),                    &
                'V+=',state%v(i,j)
            else
              write(msg(1:240),'(2(a,i4,x),2(a,f8.3,x),6(a,es11.4))') &
                'Extreme surface state detected: i=',i,'j=',j, &
                'x=',g%geoLonT(i,j),'y=',g%geoLatT(i,j), &
                'D=',g%bathyT(i,j),                      &
                'SSH=',state%sea_lev(i,j),               &
                'U-=',state%u(i-1,j),                    &
                'U+=',state%u(i,j),                      &
                'V-=',state%v(i,j-1),                    &
                'V+=',state%v(i,j)
            endif
            call mom_error(warning, trim(msg), all_print=.true.)
          elseif (numberoferrors==9) then ! Indicate once that there are more errors
            call mom_error(warning, 'There were more unreported extreme events!', all_print=.true.)
          endif ! numberOfErrors
        endif ! localError
      endif ! mask2dT
    enddo; enddo
    call sum_across_pes(numberoferrors)
    if (numberoferrors>0) then
      write(msg(1:240),'(3(a,i9,x))') 'There were a total of ',numberoferrors, &
          'locations detected with extreme surface values!'
      call mom_error(fatal, trim(msg))
    endif
  endif

  call calltree_leave("calculate_surface_state()")
end subroutine calculate_surface_state


!> End of model
subroutine mom_end(CS)
  type(mom_control_struct), pointer :: CS   !< MOM control structure

  if (cs%use_ALE_algorithm) then
    call ale_end(cs%ALE_CSp)
  endif

  dealloc_(cs%u) ; dealloc_(cs%v) ; dealloc_(cs%h)
  dealloc_(cs%uh) ; dealloc_(cs%vh)

  if (cs%use_temperature) then
    dealloc_(cs%T) ; cs%tv%T => null() ; dealloc_(cs%S) ; cs%tv%S => null()
  endif
  if (associated(cs%tv%frazil)) deallocate(cs%tv%frazil)
  if (associated(cs%tv%salt_deficit)) deallocate(cs%tv%salt_deficit)
  if (associated(cs%Hml)) deallocate(cs%Hml)

  call tracer_advect_end(cs%tracer_adv_CSp)
  call tracer_hor_diff_end(cs%tracer_diff_CSp)
  call tracer_registry_end(cs%tracer_Reg)
  call tracer_flow_control_end(cs%tracer_flow_CSp)

  if(cs%offline_tracer_mode) then
    call offline_transport_end(cs%offline_CSp)
  endif

  dealloc_(cs%uhtr) ; dealloc_(cs%vhtr)
  if (cs%split) then ; if (cs%legacy_split) then
      call end_dyn_legacy_split(cs%dyn_legacy_split_CSp)
    else
      call end_dyn_split_rk2(cs%dyn_split_RK2_CSp)
  endif ; else ; if (cs%use_RK2) then
      call end_dyn_unsplit_rk2(cs%dyn_unsplit_RK2_CSp)
    else
      call end_dyn_unsplit(cs%dyn_unsplit_CSp)
  endif ; endif
  dealloc_(cs%ave_ssh)
  if (associated(cs%update_OBC_CSp)) call obc_register_end(cs%update_OBC_CSp)

  call verticalgridend(cs%GV)
  call mom_grid_end(cs%G)

  deallocate(cs)

end subroutine mom_end

!> \namespace mom
!!
!! Modular Ocean Model (MOM) Version 6.0 (MOM6)
!!
!! \authors Alistair Adcroft, Robert Hallberg, and Stephen Griffies
!!
!!  Additional contributions from:
!!    * Whit Anderson
!!    * Brian Arbic
!!    * Will Cooke
!!    * Anand Gnanadesikan
!!    * Matthew Harrison
!!    * Mehmet Ilicak
!!    * Laura Jackson
!!    * Jasmine John
!!    * John Krasting
!!    * Zhi Liang
!!    * Bonnie Samuels
!!    * Harper Simmons
!!    * Laurent White
!!    * Niki Zadeh
!!
!!  MOM ice-shelf code was developed by
!!  * Daniel Goldberg
!!  * Robert Hallberg
!!  * Chris Little
!!  * Olga Sergienko
!!
!!  \section section_overview Overview of MOM
!!
!!  This program (MOM) simulates the ocean by numerically solving
!!  the hydrostatic primitive equations in generalized Lagrangian
!!  vertical coordinates, typically tracking stretched pressure (p*)
!!  surfaces or following isopycnals in the ocean's interior, and
!!  general orthogonal horizontal coordinates. Unlike earlier versions
!!  of MOM, in MOM6 these equations are horizontally discretized on an
!!  Arakawa C-grid.  (It remains to be seen whether a B-grid dynamic
!!  core will be revived in MOM6 at a later date; for now applications
!!  requiring a B-grid discretization should use MOM5.1.)  MOM6 offers
!!  a range of options for the physical parameterizations, from those
!!  most appropriate to highly idealized models for geophysical fluid
!!  dynamics studies to a rich suite of processes appropriate for
!!  realistic ocean simulations.  The thermodynamic options typically
!!  use conservative temperature and preformed salinity as conservative
!!  state variables and a full nonlinear equation of state, but there
!!  are also idealized adiabatic configurations of the model that use
!!  fixed density layers.  Version 6.0 of MOM continues in the long
!!  tradition of a commitment to climate-quality ocean simulations
!!  embodied in previous versions of MOM, even as it draws extensively
!!  on the lessons learned in the development of the Generalized Ocean
!!  Layered Dynamics (GOLD) ocean model, which was also primarily
!!  developed at NOAA/GFDL.  MOM has also benefited tremendously from
!!  the FMS infrastructure, which it utilizes and shares with other
!!  component models developed at NOAA/GFDL.
!!
!!    When run is isopycnal-coordinate mode, the uppermost few layers
!!  are often 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.  MOM can run with an
!!  absurdly thin minimum layer thickness. A variety of non-isopycnal
!!  vertical coordinate options are under development, but all exploit
!!  the advantages of a Lagrangian vertical coordinate, as discussed
!!  in detail by Adcroft and Hallberg (Ocean Modelling, 2006).
!!
!!    Details of the numerics and physical parameterizations are
!!  provided in the appropriate source files.  All of the available
!!  options are selected at run-time by parsing the input files,
!!  usually MOM_input and MOM_override, and the options choices are
!!  then documented for each run in MOM_param_docs.
!!
!!    MOM6 integrates the equations forward in time in three distinct
!!  phases.  In one phase, the dynamic equations for the velocities
!!  and layer thicknesses are advanced, capturing the propagation of
!!  external and internal inertia-gravity waves, Rossby waves, and
!!  other strictly adiabatic processes, including lateral stresses,
!!  vertical viscosity and momentum forcing, and interface height
!!  diffusion (commonly called Gent-McWilliams diffusion in depth-
!!  coordinate models).  In the second phase, all tracers are advected
!!  and diffused along the layers.  The third phase applies diabatic
!!  processes, vertical mixing of water properties, and perhaps
!!  vertical remapping to cause the layers to track the desired
!!  vertical coordinate.
!!
!!    The present file (MOM.F90) orchestrates 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.  As
!!  described in Hallberg and Adcroft (2009), a barotropic correction
!!  is applied to the time-mean layer velocities to ensure that the
!!  sum of the layer transports agrees with the time-mean barotropic
!!  transport, thereby ensuring that the estimates of the free surface
!!  from the sum of the layer thicknesses agrees with the final free
!!  surface height as calculated by the barotropic solver.  The
!!  barotropic and baroclinic velocities are kept consistent by
!!  recalculating the barotropic velocities from the baroclinic
!!  transports each time step. This scheme is described in Hallberg,
!!  1997, J. Comp. Phys. 135, 54-65 and in Hallberg and Adcroft, 2009,
!!  Ocean Modelling, 29, 15-26.
!!
!!    The other time integration options use non-split time stepping
!!  schemes based on the 3-step third order Runge-Kutta scheme
!!  described in Matsuno, 1966, J. Met. Soc. Japan, 44, 85-88, or on
!!  a two-step quasi-2nd order Runge-Kutta scheme.  These are much
!!  slower than the split time-stepping scheme, but they are useful
!!  for providing a more robust solution for debugging cases where the
!!  more complicated split time-stepping scheme may be giving suspect
!!  solutions.
!!
!!    There are a range of closure options available.  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, like that described in Jackson et al.
!!  (JPO, 2008).  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.  A refined bulk-mixed layer
!!  is often used to describe the planetary boundary layer in realistic
!!  ocean simulations.
!!
!!    MOM has a number of noteworthy debugging capabilities.
!!  Excessively large velocities are truncated and MOM will stop
!!  itself after a number of such instances to keep the model from
!!  crashing altogether.  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 using the DEBUG option,
!!  which causes extensive checksums to be written out along with
!!  comments indicating where in the algorithm the sums originate and
!!  what variable is being summed.  The point where these checksums
!!  differ between runs is usually a good indication of where in the
!!  code the problem lies.  All of the test cases provided with MOM
!!  are routinely tested to ensure that they give bitwise identical
!!  results regardless of the domain decomposition, or whether they
!!  use static or dynamic memory allocation.
!!
!!  \section section_structure Structure of MOM
!!
!!  About 115 other files of source code and 4 header files comprise
!!  the MOM code, although there are several hundred more files that
!!  make up the FMS infrastructure upon which MOM is built.  Each of
!!  the MOM files contains comments documenting what it does, and
!!  most of the file names are fairly self-evident. In addition, all
!!  subroutines and data types are referenced via a module use, only
!!  statement, and the module names are consistent with the file names,
!!  so it is not too hard to find the source file for a subroutine.
!!
!!    The typical MOM directory tree is as follows:
!!
!! \verbatim
!!        ../MOM
!!        |-- config_src
!!        |   |-- coupled_driver
!!        |   |-- dynamic
!!        |   `-- solo_driver
!!        |-- examples
!!        |   |-- CM2G
!!        |   |-- ...
!!        |   `-- torus_advection_test
!!        `-- src
!!            |-- core
!!            |-- diagnostics
!!            |-- equation_of_state
!!            |-- framework
!!            |-- ice_shelf
!!            |-- initialization
!!            |-- parameterizations
!!            |   |-- lateral
!!            |   `-- vertical
!!            |-- tracer
!!            `-- user
!! \endverbatim
!!
!!  Rather than describing each file here, each directory contents
!!  will be described to give a broad overview of the MOM code
!!  structure.
!!
!!    The directories under config_src contain files that are used for
!!  configuring the code, for instance for coupled or ocean-only runs.
!!  Only one or two of these directories are used in compiling any,
!!  particular run.
!!
!!  * config_src/coupled_driver:
!!    The files here are used to couple MOM as a component in a larger
!!    run driven by the FMS coupler.  This includes code that converts
!!    various forcing fields into the code structures and flux and unit
!!    conventions used by MOM, and converts the MOM surface fields
!!    back to the forms used by other FMS components.
!!
!!  * config_src/dynamic:
!!    The only file here is the version of MOM_memory.h that is used
!!    for dynamic memory configurations of MOM.
!!
!!  * config_src/solo_driver:
!!    The files here are include the _main driver that is used when
!!    MOM is configured as an ocean-only model, as well as the files
!!    that specify the surface forcing in this configuration.
!!
!!    The directories under examples provide a large number of working
!!  configurations of MOM, along with reference solutions for several
!!  different compilers on GFDL's latest large computer.  The versions
!!  of MOM_memory.h in these directories need not be used if dynamic
!!  memory allocation is desired, and the answers should be unchanged.
!!
!!    The directories under src contain most of the MOM files.  These
!!  files are used in every configuration using MOM.
!!
!!  * src/core:
!!    The files here constitute the MOM dynamic core.  This directory
!!    also includes files with the types that describe the model's
!!    lateral grid and have defined types that are shared across
!!    various MOM modules to allow for more succinct and flexible
!!    subroutine argument lists.
!!
!!  * src/diagnostics:
!!    The files here calculate various diagnostics that are anciliary
!!    to the model itself.  While most of these diagnostics do not
!!    directly affect the model's solution, there are some, like the
!!    calculation of the deformation radius, that are used in some
!!    of the process parameterizations.
!!
!!  * src/equation_of_state:
!!    These files describe the physical properties of sea-water,
!!    including both the equation of state and when it freezes.
!!
!!  * src/framework:
!!    These files provide infrastructure utilities for MOM.  Many are
!!    simply wrappers for capabilities provided by FMS, although others
!!    provide capabilities (like the file_parser) that are unique to
!!    MOM. When MOM is adapted to use a modeling infrastructure
!!    distinct from FMS, most of the required changes are in this
!!    directory.
!!
!!  * src/initialization:
!!    These are the files that are used to initialize the MOM grid
!!    or provide the initial physical state for MOM.  These files are
!!    not intended to be modified, but provide a means for calling
!!    user-specific initialization code like the examples in src/user.
!!
!!  * src/parameterizations/lateral:
!!    These files implement a number of quasi-lateral (along-layer)
!!    process parameterizations, including lateral viscosities,
!!    parameterizations of eddy effects, and the calculation of tidal
!!    forcing.
!!
!!  * src/parameterizations/vertical:
!!    These files implement a number of vertical mixing or diabatic
!!    processes, including the effects of vertical viscosity and
!!    code to parameterize the planetary boundary layer.  There is a
!!    separate driver that orchestrates this portion of the algorithm,
!!    and there is a diversity of parameterizations to be found here.
!!
!!  * src/tracer:
!!    These files handle the lateral transport and diffusion of
!!    tracers, or are the code to implement various passive tracer
!!    packages.  Additional tracer packages are readily accomodated.
!!
!!  * src/user:
!!    These are either stub routines that a user could use to change
!!    the model's initial conditions or forcing, or are examples that
!!    implement specific test cases.  These files can easily  be hand
!!    edited to create new analytically specified configurations.
!!
!!
!!  Most simulations can be set up by modifying only the files
!!  MOM_input, and possibly one or two of the files in src/user.
!!  In addition, the diag_table (MOM_diag_table) will commonly be
!!  modified to tailor the output to the needs of the question at
!!  hand.  The FMS utility mkmf works with a file called path_names
!!  to build an appropriate makefile, and path_names should be edited
!!  to reflect the actual location of the desired source code.
!!
!!
!!  There are 3 publicly visible subroutines in this file (MOM.F90).
!!  * step_MOM steps MOM over a specified interval of time.
!!  * MOM_initialize calls initialize and does other initialization
!!    that does not warrant user modification.
!!  * calculate_surface_state determines the surface (bulk mixed layer
!!    if traditional isoycnal vertical coordinate) properties of the
!!    current model state and packages pointers to these fields into an
!!    exported structure.
!!
!!    The remaining subroutines in this file (src/core/MOM.F90) are:
!!  * find_total_transport determines the barotropic mass transport.
!!  * register_diags registers many diagnostic fields for the dynamic
!!    solver, or of the main model variables.
!!  * MOM_timing_init initializes various CPU time clocks.
!!  * write_static_fields writes out various time-invariant fields.
!!  * set_restart_fields is used to specify those fields that are
!!    written to and read from the restart file.
!!
!!  \section section_heat_budget Diagnosing MOM heat budget
!!
!!  Here are some example heat budgets for the ALE version of MOM6.
!!
!!  \subsection subsection_2d_heat_budget Depth integrated heat budget
!!
!!  Depth integrated heat budget diagnostic for MOM.
!!
!! * OPOTTEMPTEND_2d = T_ADVECTION_XY_2d + OPOTTEMPPMDIFF_2d + HFDS + HFGEOU
!!
!! * T_ADVECTION_XY_2d = horizontal advection
!! * OPOTTEMPPMDIFF_2d = neutral diffusion
!! * HFDS              = net surface boundary heat flux
!! * HFGEOU            = geothermal heat flux
!!
!! * HFDS = net surface boundary heat flux entering the ocean
!!        = rsntds + rlntds + hfls + hfss + heat_pme + hfsifrazil
!!
!! * More heat flux cross-checks
!!   * hfds     = net_heat_coupler + hfsifrazil + heat_pme
!!   * heat_pme = heat_content_surfwater
!!              = heat_content_massin + heat_content_massout
!!              = heat_content_fprec + heat_content_cond + heat_content_vprec
!!               + hfrunoffds + hfevapds + hfrainds
!!
!!  \subsection subsection_3d_heat_budget Depth integrated heat budget
!!
!!  Here is an example 3d heat budget diagnostic for MOM.
!!
!! * OPOTTEMPTEND = T_ADVECTION_XY + TH_TENDENCY_VERT_REMAP + OPOTTEMPDIFF + OPOTTEMPPMDIFF
!!                + BOUNDARY_FORCING_HEAT_TENDENCY + FRAZIL_HEAT_TENDENCY
!!
!! * OPOTTEMPTEND                   = net tendency of heat as diagnosed in MOM.F90
!! * T_ADVECTION_XY                 = heating of a cell from lateral advection
!! * TH_TENDENCY_VERT_REMAP         = heating of a cell from vertical remapping
!! * OPOTTEMPDIFF                   = heating of a cell from diabatic diffusion
!! * OPOTTEMPPMDIFF                 = heating of a cell from neutral diffusion
!! * BOUNDARY_FORCING_HEAT_TENDENCY = heating of cell from boundary fluxes
!! * FRAZIL_HEAT_TENDENCY           = heating of cell from frazil
!!
!! * TH_TENDENCY_VERT_REMAP has zero vertical sum, as it redistributes heat in vertical.
!!
!! * OPOTTEMPDIFF has zero vertical sum, as it redistributes heat in the vertical.
!!
!! * BOUNDARY_FORCING_HEAT_TENDENCY generally has 3d structure, with k > 1 contributions from
!!   penetrative shortwave, and from other fluxes for the case when layers are tiny, in which
!!   case MOM6 partitions tendencies into k > 1 layers.
!!
!! * FRAZIL_HEAT_TENDENCY generally has 3d structure, since MOM6 frazil calculation checks the
!!   full ocean column.
!!
!! * FRAZIL_HEAT_TENDENCY[k=\@sum] = HFSIFRAZIL = column integrated frazil heating.
!!
!! * HFDS = FRAZIL_HEAT_TENDENCY[k=\@sum] + BOUNDARY_FORCING_HEAT_TENDENCY[k=\@sum]
!!
!!  Here is an example 2d heat budget (depth summed) diagnostic for MOM.
!!
!! * OPOTTEMPTEND_2d = T_ADVECTION_XY_2d + OPOTTEMPPMDIFF_2d + HFDS
!!
!!
!!  Here is an example 3d salt budget diagnostic for MOM.
!!
!! * OSALTTEND = S_ADVECTION_XY + SH_TENDENCY_VERT_REMAP + OSALTDIFF + OSALTPMDIFF
!!                + BOUNDARY_FORCING_SALT_TENDENCY
!!
!! * OSALTTEND                      = net tendency of salt as diagnosed in MOM.F90
!! * S_ADVECTION_XY                 = salt convergence to cell from lateral advection
!! * SH_TENDENCY_VERT_REMAP         = salt convergence to cell from vertical remapping
!! * OSALTDIFF                      = salt convergence to cell from diabatic diffusion
!! * OSALTPMDIFF                    = salt convergence to cell from neutral diffusion
!! * BOUNDARY_FORCING_SALT_TENDENCY = salt convergence to cell from boundary fluxes
!!
!! * SH_TENDENCY_VERT_REMAP has zero vertical sum, as it redistributes salt in vertical.
!!
!! * OSALTDIFF has zero vertical sum, as it redistributes salt in the vertical.
!!
!! * BOUNDARY_FORCING_SALT_TENDENCY generally has 3d structure, with k > 1 contributions from
!!   the case when layers are tiny, in which case MOM6 partitions tendencies into k > 1 layers.
!!
!! * SFDSI = BOUNDARY_FORCING_SALT_TENDENCY[k=\@sum]
!!
!!  Here is an example 2d salt budget (depth summed) diagnostic for MOM.
!!
!! * OSALTTEND_2d = S_ADVECTION_XY_2d + OSALTPMDIFF_2d + SFDSI (+ SALT_FLUX_RESTORE)
!!
!!
!!


end module mom