MOM_sum_output.F90

!> Reports integrated quantities for monitoring the model state
module mom_sum_output

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

use iso_fortran_env, only : int64
use mom_coms, only : sum_across_pes, pe_here, root_pe, num_pes, max_across_pes
use mom_coms, only : reproducing_sum, reproducing_sum_efp, efp_to_real, real_to_efp
use mom_coms, only : efp_type, operator(+), operator(-), assignment(=), efp_sum_across_pes
use mom_error_handler, only : mom_error, fatal, warning, is_root_pe, mom_mesg
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_forcing_type, only : forcing
use mom_grid, only : ocean_grid_type
use mom_interface_heights, only : find_eta
use mom_io, only : create_file, fieldtype, flush_file, open_file, reopen_file
use mom_io, only : file_exists, slasher, vardesc, var_desc, write_field, get_filename_appendix
use mom_io, only : append_file, ascii_file, single_file, writeonly_file
use mom_open_boundary, only : ocean_obc_type, obc_segment_type
use mom_open_boundary, only : obc_direction_e, obc_direction_w, obc_direction_n, obc_direction_s
use mom_time_manager, only : time_type, get_time, get_date, set_time, operator(>)
use mom_time_manager, only : operator(+), operator(-), operator(*), operator(/)
use mom_time_manager, only : operator(/=), operator(<=), operator(>=), operator(<)
use mom_time_manager, only : get_calendar_type, time_type_to_real, no_calendar
use mom_tracer_flow_control, only : tracer_flow_control_cs, call_tracer_stocks
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : surface, thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mpp_mod, only : mpp_chksum

use netcdf

implicit none ; private

#include <MOM_memory.h>

public write_energy, accumulate_net_input, mom_sum_output_init

! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.

integer, parameter :: num_fields = 17 !< Number of diagnostic fields
character (*), parameter :: depth_chksum_attr = "bathyT_checksum"
                                      !< Checksum attribute name of G%bathyT
                                      !! over the compute domain
character (*), parameter :: area_chksum_attr = "mask2dT_areaT_checksum"
                                      !< Checksum attribute of name of
                                      !! G%mask2dT * G%areaT over the compute
                                      !! domain

!> A list of depths and corresponding globally integrated ocean area at each
!! depth and the ocean volume below each depth.
type :: depth_list
  real :: depth       !< A depth [Z ~> m].
  real :: area        !< The cross-sectional area of the ocean at that depth [L2 ~> m2].
  real :: vol_below   !< The ocean volume below that depth [Z m2 ~> m3].
end type depth_list

!> The control structure for the MOM_sum_output module
type, public :: sum_output_cs ; private
  type(depth_list), pointer, dimension(:) :: dl => null() !< The sorted depth list.
  integer :: list_size          !< length of sorting vector <= niglobal*njglobal

  integer, allocatable, dimension(:) :: lh
                                !< This saves the entry in DL with a volume just
                                !! less than the volume of fluid below the interface.
  logical :: do_ape_calc        !<   If true, calculate the available potential energy of the
                                !! interfaces.  Disabling this reduces the memory footprint of
                                !! high-PE-count models dramatically.
  logical :: read_depth_list    !<   Read the depth list from a file if it exists
                                !! and write it if it doesn't.
  character(len=200) :: depth_list_file  !< The name of the depth list file.
  real    :: d_list_min_inc     !<  The minimum increment [Z ~> m], between the depths of the
                                !! entries in the depth-list file, 0 by default.
  logical :: require_depth_list_chksum
                                !< Require matching checksums in Depth_list.nc when reading
                                !! the file.
  logical :: update_depth_list_chksum
                                !< Automatically update the Depth_list.nc file if the
                                !! checksums are missing or do not match current values.
  logical :: use_temperature    !<   If true, temperature and salinity are state variables.
  type(efp_type) :: fresh_water_in_efp !< The total mass of fresh water added by surface fluxes on
                                  !! this PE since the last time that write_energy was called [kg].
  type(efp_type) :: net_salt_in_efp !< The total salt added by surface fluxes on this PE since
                                  !! the last time that write_energy was called [ppt kg].
  type(efp_type) :: net_heat_in_efp !<  The total heat added by surface fluxes on this PE since
                                  !! the last time that write_energy was called [J].
  type(efp_type) :: heat_prev_efp !<  The total amount of heat in the ocean the last
                                  !! time that write_energy was called [J].
  type(efp_type) :: salt_prev_efp !< The total amount of salt in the ocean the last
                                  !! time that write_energy was called [ppt kg].
  type(efp_type) :: mass_prev_efp !< The total ocean mass the last time that
                                  !! write_energy was called [kg].
  real    :: dt_in_t            !< The baroclinic dynamics time step [T ~> s].

  type(time_type) :: energysavedays            !< The interval between writing the energies
                                               !! and other integral quantities of the run.
  type(time_type) :: energysavedays_geometric  !< The starting interval for computing a geometric
                                               !! progression of time deltas between calls to
                                               !! write_energy. This interval will increase by a factor of 2.
                                               !! after each call to write_energy.
  logical         :: energysave_geometric      !< Logical to control whether calls to write_energy should
                                               !! follow a geometric progression
  type(time_type) :: write_energy_time         !< The next time to write to the energy file.
  type(time_type) :: geometric_end_time        !< Time at which to stop the geometric progression
                                               !! of calls to write_energy and revert to the standard
                                               !! energysavedays interval

  real    :: timeunit           !<  The length of the units for the time axis [s].
  logical :: date_stamped_output !< If true, use dates (not times) in messages to stdout.
  type(time_type) :: start_time !< The start time of the simulation.
                                ! Start_time is set in MOM_initialization.F90
  integer, pointer :: ntrunc => null() !< The number of times the velocity has been
                                !! truncated since the last call to write_energy.
  real    :: max_energy         !< The maximum permitted energy per unit mass.  If there is
                                !! more energy than this, the model should stop [m2 s-2].
  integer :: maxtrunc           !< The number of truncations per energy save
                                !! interval at which the run is stopped.
  logical :: write_stocks       !< If true, write the integrated tracer amounts
                                !! to stdout when the energy files are written.
  integer :: previous_calls = 0 !< The number of times write_energy has been called.
  integer :: prev_n = 0         !< The value of n from the last call.
  integer :: fileenergy_nc      !< NetCDF id of the energy file.
  integer :: fileenergy_ascii   !< The unit number of the ascii version of the energy file.
  type(fieldtype), dimension(NUM_FIELDS+MAX_FIELDS_) :: &
             fields             !< fieldtype variables for the output fields.
  character(len=200) :: energyfile  !< The name of the energy file with path.
end type sum_output_cs

contains

!> MOM_sum_output_init initializes the parameters and settings for the MOM_sum_output module.
subroutine mom_sum_output_init(G, US, param_file, directory, ntrnc, &
                               Input_start_time, CS)
  type(ocean_grid_type),  intent(in)    :: g          !< The ocean's grid structure.
  type(unit_scale_type),  intent(in)    :: us         !< A dimensional unit scaling type
  type(param_file_type),  intent(in)    :: param_file !< A structure to parse for run-time
                                                      !! parameters.
  character(len=*),       intent(in)    :: directory  !< The directory where the energy file goes.
  integer, target,        intent(inout) :: ntrnc      !< The integer that stores the number of times
                                                      !! the velocity has been truncated since the
                                                      !! last call to write_energy.
  type(time_type),        intent(in)    :: input_start_time !< The start time of the simulation.
  type(sum_output_cs),    pointer       :: cs         !< A pointer that is set to point to the
                                                      !! control structure for this module.
  ! Local variables
  real :: time_unit ! The time unit in seconds for ENERGYSAVEDAYS.
  real :: rho_0     ! A reference density [kg m-3]
  real :: maxvel    ! The maximum permitted velocity [m s-1]
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_sum_output" ! This module's name.
  character(len=200) :: energyfile  ! The name of the energy file.
  character(len=32) :: filename_appendix = '' !fms appendix to filename for ensemble runs

  if (associated(cs)) then
    call mom_error(warning, "MOM_sum_output_init called with associated control structure.")
    return
  endif
  allocate(cs)

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "CALCULATE_APE", cs%do_APE_calc, &
                 "If true, calculate the available potential energy of "//&
                 "the interfaces.  Setting this to false reduces the "//&
                 "memory footprint of high-PE-count models dramatically.", &
                 default=.true.)
  call get_param(param_file, mdl, "WRITE_STOCKS", cs%write_stocks, &
                 "If true, write the integrated tracer amounts to stdout "//&
                 "when the energy files are written.", default=.true.)
  call get_param(param_file, mdl, "ENABLE_THERMODYNAMICS", cs%use_temperature, &
                 "If true, Temperature and salinity are used as state "//&
                 "variables.", default=.true.)
  call get_param(param_file, mdl, "DT", cs%dt_in_T, &
                 "The (baroclinic) dynamics time step.", &
                 units="s", scale=us%s_to_T, fail_if_missing=.true.)
  call get_param(param_file, mdl, "MAXTRUNC", cs%maxtrunc, &
                 "The run will be stopped, and the day set to a very "//&
                 "large value if the velocity is truncated more than "//&
                 "MAXTRUNC times between energy saves.  Set MAXTRUNC to 0 "//&
                 "to stop if there is any truncation of velocities.", &
                 units="truncations save_interval-1", default=0)

  call get_param(param_file, mdl, "MAX_ENERGY", cs%max_Energy, &
                 "The maximum permitted average energy per unit mass; the "//&
                 "model will be stopped if there is more energy than "//&
                 "this.  If zero or negative, this is set to 10*MAXVEL^2.", &
                 units="m2 s-2", default=0.0)
  if (cs%max_Energy <= 0.0) then
    call get_param(param_file, mdl, "MAXVEL", maxvel, &
                 "The maximum velocity allowed before the velocity "//&
                 "components are truncated.", units="m s-1", default=3.0e8)
    cs%max_Energy = 10.0 * maxvel**2
    call log_param(param_file, mdl, "MAX_ENERGY as used", cs%max_Energy)
  endif

  call get_param(param_file, mdl, "ENERGYFILE", energyfile, &
                 "The file to use to write the energies and globally "//&
                 "summed diagnostics.", default="ocean.stats")

  !query fms_io if there is a filename_appendix (for ensemble runs)
  call get_filename_appendix(filename_appendix)
  if (len_trim(filename_appendix) > 0) then
     energyfile = trim(energyfile) //'.'//trim(filename_appendix)
  endif

  cs%energyfile = trim(slasher(directory))//trim(energyfile)
  call log_param(param_file, mdl, "output_path/ENERGYFILE", cs%energyfile)
#ifdef STATSLABEL
  cs%energyfile = trim(cs%energyfile)//"."//trim(adjustl(statslabel))
#endif

  call get_param(param_file, mdl, "DATE_STAMPED_STDOUT", cs%date_stamped_output, &
                 "If true, use dates (not times) in messages to stdout", &
                 default=.true.)
  call get_param(param_file, mdl, "TIMEUNIT", cs%Timeunit, &
                 "The time unit in seconds a number of input fields", &
                 units="s", default=86400.0)
  if (cs%Timeunit < 0.0) cs%Timeunit = 86400.0



  if (cs%do_APE_calc) then
    call get_param(param_file, mdl, "READ_DEPTH_LIST", cs%read_depth_list, &
                   "Read the depth list from a file if it exists or "//&
                   "create that file otherwise.", default=.false.)
    call get_param(param_file, mdl, "DEPTH_LIST_MIN_INC", cs%D_list_min_inc, &
                   "The minimum increment between the depths of the "//&
                   "entries in the depth-list file.", &
                   units="m", default=1.0e-10, scale=us%m_to_Z)
    if (cs%read_depth_list) then
      call get_param(param_file, mdl, "DEPTH_LIST_FILE", cs%depth_list_file, &
                   "The name of the depth list file.", default="Depth_list.nc")
      if (scan(cs%depth_list_file,'/') == 0) &
        cs%depth_list_file = trim(slasher(directory))//trim(cs%depth_list_file)

      call get_param(param_file, mdl, "REQUIRE_DEPTH_LIST_CHECKSUMS", &
                     cs%require_depth_list_chksum, &
                 "Require that matching checksums be in Depth_list.nc "//&
                 "when reading the file.", default=.true.)
      if (.not. cs%require_depth_list_chksum) &
        call get_param(param_file, mdl, "UPDATE_DEPTH_LIST_CHECKSUMS", &
                     cs%update_depth_list_chksum, &
                 "Automatically update the Depth_list.nc file if the "//&
                 "checksums are missing or do not match current values.", &
                 default=.false.)
    endif

    allocate(cs%lH(g%ke))
    call depth_list_setup(g, us, cs)
  else
    cs%list_size = 0
  endif

  call get_param(param_file, mdl, "TIMEUNIT", time_unit, &
                 "The time unit for ENERGYSAVEDAYS.", &
                 units="s", default=86400.0)
  call get_param(param_file, mdl, "ENERGYSAVEDAYS",cs%energysavedays, &
                 "The interval in units of TIMEUNIT between saves of the "//&
                 "energies of the run and other globally summed diagnostics.",&
                 default=set_time(0,days=1), timeunit=time_unit)
  call get_param(param_file, mdl, "ENERGYSAVEDAYS_GEOMETRIC",cs%energysavedays_geometric, &
                 "The starting interval in units of TIMEUNIT for the first call "//&
                 "to save the energies of the run and other globally summed diagnostics. "//&
                 "The interval increases by a factor of 2. after each call to write_energy.",&
                 default=set_time(seconds=0), timeunit=time_unit)

  if ((time_type_to_real(cs%energysavedays_geometric) > 0.) .and. &
     (cs%energysavedays_geometric < cs%energysavedays)) then
         cs%energysave_geometric = .true.
  else
         cs%energysave_geometric = .false.
  endif

  cs%Start_time = input_start_time
  cs%ntrunc => ntrnc

end subroutine mom_sum_output_init

!> MOM_sum_output_end deallocates memory used by the MOM_sum_output module.
subroutine mom_sum_output_end(CS)
  type(Sum_output_CS), pointer :: CS  !< The control structure returned by a
                                      !! previous call to MOM_sum_output_init.
  if (associated(cs)) then
    if (cs%do_APE_calc) then
      deallocate(cs%lH, cs%DL)
    endif

    deallocate(cs)
  endif
end subroutine mom_sum_output_end

!>  This subroutine calculates and writes the total model energy, the energy and
!! mass of each layer, and other globally integrated  physical quantities.
subroutine write_energy(u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_forcing)
  type(ocean_grid_type),   intent(in)    :: g   !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)    :: gv  !< The ocean's vertical grid structure.
  type(unit_scale_type),   intent(in)    :: us  !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: u   !< The zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: v   !< The meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(in)    :: h   !< Layer thicknesses [H ~> m or kg m-2].
  type(thermo_var_ptrs),   intent(in)    :: tv  !< A structure pointing to various
                                                !! thermodynamic variables.
  type(time_type),         intent(in)    :: day !< The current model time.
  integer,                 intent(in)    :: n   !< The time step number of the
                                                !! current execution.
  type(sum_output_cs),     pointer       :: cs  !< The control structure returned by a
                                                !! previous call to MOM_sum_output_init.
  type(tracer_flow_control_cs), &
                    optional, pointer    :: tracer_csp !< tracer control structure.
  type(ocean_obc_type),         &
                    optional, pointer    :: obc !< Open boundaries control structure.
  type(time_type),  optional, intent(in) :: dt_forcing !< The forcing time step
  ! Local variables
  real :: eta(szi_(g),szj_(g),szk_(g)+1) ! The height of interfaces [Z ~> m].
  real :: areatm(szi_(g),szj_(g)) ! A masked version of areaT [L2 ~> m2].
  real :: ke(szk_(g))  ! The total kinetic energy of a layer [J].
  real :: pe(szk_(g)+1)! The available potential energy of an interface [J].
  real :: ke_tot       ! The total kinetic energy [J].
  real :: pe_tot       ! The total available potential energy [J].
  real :: z_0ape(szk_(g)+1) ! The uniform depth which overlies the same
                       ! volume as is below an interface [Z ~> m].
  real :: h_0ape(szk_(g)+1) ! A version of Z_0APE, converted to m, usually positive.
  real :: toten        ! The total kinetic & potential energies of
                       ! all layers [J] (i.e. kg m2 s-2).
  real :: en_mass      ! The total kinetic and potential energies divided by
                       ! the total mass of the ocean [m2 s-2].
  real :: vol_lay(szk_(g))  ! The volume of fluid in a layer [Z L2 ~> m3].
  real :: volbelow     ! The volume of all layers beneath an interface [Z L2 ~> m3].
  real :: mass_lay(szk_(g)) ! The mass of fluid in a layer [kg].
  real :: mass_tot     ! The total mass of the ocean [kg].
  real :: vol_tot      ! The total ocean volume [m3].
  real :: mass_chg     ! The change in total ocean mass of fresh water since
                       ! the last call to this subroutine [kg].
  real :: mass_anom    ! The change in fresh water that cannot be accounted for
                       ! by the surface fluxes [kg].
  real :: salt         ! The total amount of salt in the ocean [ppt kg].
  real :: salt_chg     ! The change in total ocean salt since the last call
                       ! to this subroutine [ppt kg].
  real :: salt_anom    ! The change in salt that cannot be accounted for by
                       ! the surface fluxes [ppt kg].
  real :: salin        ! The mean salinity of the ocean [ppt].
  real :: salin_chg    ! The change in total salt since the last call
                       ! to this subroutine divided by total mass [ppt].
  real :: salin_anom   ! The change in total salt that cannot be accounted for by
                       ! the surface fluxes divided by total mass [ppt].
  real :: heat         ! The total amount of Heat in the ocean [J].
  real :: heat_chg     ! The change in total ocean heat since the last call to this subroutine [J].
  real :: heat_anom    ! The change in heat that cannot be accounted for by the surface fluxes [J].
  real :: temp         ! The mean potential temperature of the ocean [degC].
  real :: temp_chg     ! The change in total heat divided by total heat capacity
                       ! of the ocean since the last call to this subroutine, degC.
  real :: temp_anom    ! The change in total heat that cannot be accounted for
                       ! by the surface fluxes, divided by the total heat
                       ! capacity of the ocean [degC].
  real :: hint         ! The deviation of an interface from H [Z ~> m].
  real :: hbot         ! 0 if the basin is deeper than H, or the
                       ! height of the basin depth over H otherwise [Z ~> m].
                       ! This makes PE only include real fluid.
  real :: hbelow       ! The depth of fluid in all layers beneath an interface [Z ~> m].
  type(efp_type) :: &
    mass_efp, &        ! The total mass of the ocean in extended fixed point form [kg].
    salt_efp, &        ! The total amount of salt in the ocean in extended fixed point form [ppt kg].
    heat_efp, &        ! The total amount of heat in the ocean in extended fixed point form [J].
    salt_chg_efp, &    ! The change in total ocean salt since the last call to this subroutine [ppt kg].
    heat_chg_efp, &    ! The change in total ocean heat since the last call to this subroutine [J].
    mass_chg_efp, &    ! The change in total ocean mass of fresh water since
                       ! the last call to this subroutine [kg].
    salt_anom_efp, &   ! The change in salt that cannot be accounted for by the surface
                       ! fluxes [ppt kg].
    heat_anom_efp, &   ! The change in heat that cannot be accounted for by the surface fluxes [J].
    mass_anom_efp      ! The change in fresh water that cannot be accounted for by the surface
                       ! fluxes [kg].
  type(efp_type), dimension(5) :: efp_list ! An array of EFP types for joint global sums.
  real :: cfl_iarea    ! Direction-based inverse area used in CFL test [L-2].
  real :: cfl_trans    ! A transport-based definition of the CFL number [nondim].
  real :: cfl_lin      ! A simpler definition of the CFL number [nondim].
  real :: max_cfl(2)   ! The maxima of the CFL numbers [nondim].
  real :: irho0        ! The inverse of the reference density [m3 kg-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    tmp1               ! A temporary array
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
    pe_pt              ! The potential energy at each point [J].
  real, dimension(SZI_(G),SZJ_(G)) :: &
    temp_int, salt_int ! Layer and cell integrated heat and salt [J] and [g Salt].
  real :: hl2_to_kg    ! A conversion factor from a thickness-volume to mass [kg H-1 L-2 ~> kg m-3 or 1]
  real :: ke_scale_factor   ! The combination of unit rescaling factors in the kinetic energy
                            ! calculation [kg T2 H-1 L-2 s-2 ~> kg m-3 or nondim]
  real :: pe_scale_factor   ! The combination of unit rescaling factors in the potential energy
                            ! calculation [kg T2 R-1 Z-1 L-2 s-2 ~> nondim]
  integer :: num_nc_fields  ! The number of fields that will actually go into
                            ! the NetCDF file.
  integer :: i, j, k, is, ie, js, je, ns, nz, m, isq, ieq, jsq, jeq, isr, ier, jsr, jer
  integer :: l, lbelow, labove   ! indices of deep_area_vol, used to find Z_0APE.
                                 ! lbelow & labove are lower & upper limits for l
                                 ! in the search for the entry in lH to use.
  integer :: start_of_day, num_days
  real    :: reday, var
  character(len=240) :: energypath_nc
  character(len=200) :: mesg
  character(len=32)  :: mesg_intro, time_units, day_str, n_str, date_str
  logical :: date_stamped
  type(time_type) :: dt_force ! A time_type version of the forcing timestep.
  real :: tr_stocks(max_fields_)
  real :: tr_min(max_fields_), tr_max(max_fields_)
  real :: tr_min_x(max_fields_), tr_min_y(max_fields_), tr_min_z(max_fields_)
  real :: tr_max_x(max_fields_), tr_max_y(max_fields_), tr_max_z(max_fields_)
  logical :: tr_minmax_got(max_fields_) = .false.
  character(len=40), dimension(MAX_FIELDS_) :: &
    tr_names, tr_units
  integer :: ntr_stocks
  integer :: iyear, imonth, iday, ihour, iminute, isecond, itick ! For call to get_date()
  logical :: local_open_bc
  type(obc_segment_type), pointer :: segment => null()

 ! A description for output of each of the fields.
  type(vardesc) :: vars(num_fields+max_fields_)

  ! write_energy_time is the next integral multiple of energysavedays.
  dt_force = set_time(seconds=2) ; if (present(dt_forcing)) dt_force = dt_forcing
  if (cs%previous_calls == 0) then
    if (cs%energysave_geometric) then
      if (cs%energysavedays_geometric < cs%energysavedays) then
        cs%write_energy_time = day + cs%energysavedays_geometric
        cs%geometric_end_time = cs%Start_time + cs%energysavedays * &
          (1 + (day - cs%Start_time) / cs%energysavedays)
      else
        cs%write_energy_time = cs%Start_time + cs%energysavedays * &
          (1 + (day - cs%Start_time) / cs%energysavedays)
      endif
    else
      cs%write_energy_time = cs%Start_time + cs%energysavedays * &
        (1 + (day - cs%Start_time) / cs%energysavedays)
    endif
  elseif (day + (dt_force/2) <= cs%write_energy_time) then
    return  ! Do not write this step
  else ! Determine the next write time before proceeding
    if (cs%energysave_geometric) then
      if (cs%write_energy_time + cs%energysavedays_geometric >= &
          cs%geometric_end_time) then
        cs%write_energy_time = cs%geometric_end_time
        cs%energysave_geometric = .false.  ! stop geometric progression
      else
        cs%write_energy_time = cs%write_energy_time + cs%energysavedays_geometric
      endif
      cs%energysavedays_geometric = cs%energysavedays_geometric*2
    else
      cs%write_energy_time = cs%write_energy_time + cs%energysavedays
    endif
  endif

  num_nc_fields = 17
  if (.not.cs%use_temperature) num_nc_fields = 11
  vars(1) = var_desc("Ntrunc","Nondim","Number of Velocity Truncations",'1','1')
  vars(2) = var_desc("En","Joules","Total Energy",'1','1')
  vars(3) = var_desc("APE","Joules","Total Interface APE",'1','i')
  vars(4) = var_desc("KE","Joules","Total Layer KE",'1','L')
  vars(5) = var_desc("H0","meter","Zero APE Depth of Interface",'1','i')
  vars(6) = var_desc("Mass_lay","kg","Total Layer Mass",'1','L')
  vars(7) = var_desc("Mass","kg","Total Mass",'1','1')
  vars(8) = var_desc("Mass_chg","kg","Total Mass Change between Entries",'1','1')
  vars(9) = var_desc("Mass_anom","kg","Anomalous Total Mass Change",'1','1')
  vars(10) = var_desc("max_CFL_trans","Nondim","Maximum finite-volume CFL",'1','1')
  vars(11) = var_desc("max_CFL_lin","Nondim","Maximum finite-difference CFL",'1','1')
  if (cs%use_temperature) then
    vars(12) = var_desc("Salt","kg","Total Salt",'1','1')
    vars(13) = var_desc("Salt_chg","kg","Total Salt Change between Entries",'1','1')
    vars(14) = var_desc("Salt_anom","kg","Anomalous Total Salt Change",'1','1')
    vars(15) = var_desc("Heat","Joules","Total Heat",'1','1')
    vars(16) = var_desc("Heat_chg","Joules","Total Heat Change between Entries",'1','1')
    vars(17) = var_desc("Heat_anom","Joules","Anomalous Total Heat Change",'1','1')
  endif

  local_open_bc = .false.
  if (present(obc)) then ; if (associated(obc)) then
    local_open_bc = (obc%open_u_BCs_exist_globally .or. obc%open_v_BCs_exist_globally)
  endif ; endif

  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
  isr = is - (g%isd-1) ; ier = ie - (g%isd-1) ; jsr = js - (g%jsd-1) ; jer = je - (g%jsd-1)


  hl2_to_kg = gv%H_to_kg_m2*us%L_to_m**2

  if (.not.associated(cs)) call mom_error(fatal, &
         "write_energy: Module must be initialized before it is used.")

  do j=js,je ; do i=is,ie
    areatm(i,j) = g%mask2dT(i,j)*g%areaT(i,j)
  enddo ; enddo

  if (gv%Boussinesq) then
    tmp1(:,:,:) = 0.0
    do k=1,nz ; do j=js,je ; do i=is,ie
      tmp1(i,j,k) = h(i,j,k) * (hl2_to_kg*areatm(i,j))
    enddo ; enddo ; enddo

    ! This block avoids using the points beyond an open boundary condition
    ! in the accumulation of mass, but perhaps it would be unnecessary if there
    ! were a more judicious use of masks in the loops 4 or 7 lines above.
    if (local_open_bc) then
      do ns=1, obc%number_of_segments
        segment => obc%segment(ns)
        if (.not. segment%on_pe .or. segment%specified) cycle
        i=segment%HI%IsdB ; j=segment%HI%JsdB
        if (segment%direction == obc_direction_e) then
          do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
            tmp1(i+1,j,k) = 0.0
          enddo ; enddo
        elseif (segment%direction == obc_direction_w) then
          do k=1,nz ; do j=segment%HI%jsd,segment%HI%jed
            tmp1(i,j,k) = 0.0
          enddo ; enddo
        elseif (segment%direction == obc_direction_n) then
          do k=1,nz ; do i=segment%HI%isd,segment%HI%ied
            tmp1(i,j+1,k) = 0.0
          enddo ; enddo
        elseif (segment%direction == obc_direction_s) then
          do k=1,nz ; do i=segment%HI%isd,segment%HI%ied
            tmp1(i,j,k) = 0.0
          enddo ; enddo
        endif
      enddo
    endif

    mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, efp_sum=mass_efp)
    do k=1,nz ; vol_lay(k) = (us%m_to_L**2*gv%H_to_Z/gv%H_to_kg_m2)*mass_lay(k) ; enddo
  else
    tmp1(:,:,:) = 0.0
    if (cs%do_APE_calc) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        tmp1(i,j,k) = hl2_to_kg * h(i,j,k) * areatm(i,j)
      enddo ; enddo ; enddo
      mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, efp_sum=mass_efp)

      call find_eta(h, tv, g, gv, us, eta)
      do k=1,nz ; do j=js,je ; do i=is,ie
        tmp1(i,j,k) = us%Z_to_m*us%L_to_m**2*(eta(i,j,k)-eta(i,j,k+1)) * areatm(i,j)
      enddo ; enddo ; enddo
      vol_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=vol_lay)
      do k=1,nz ; vol_lay(k) = us%m_to_Z*us%m_to_L**2 * vol_lay(k) ; enddo
    else
      do k=1,nz ; do j=js,je ; do i=is,ie
        tmp1(i,j,k) = hl2_to_kg * h(i,j,k) * areatm(i,j)
      enddo ; enddo ; enddo
      mass_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=mass_lay, efp_sum=mass_efp)
      do k=1,nz ; vol_lay(k) = us%m_to_Z*us%m_to_L**2*us%kg_m3_to_R * (mass_lay(k) / gv%Rho0) ; enddo
    endif
  endif ! Boussinesq

  ntr_stocks = 0
  if (present(tracer_csp)) then
    call call_tracer_stocks(h, tr_stocks, g, gv, tracer_csp, stock_names=tr_names, &
                            stock_units=tr_units, num_stocks=ntr_stocks,&
                            got_min_max=tr_minmax_got, global_min=tr_min, global_max=tr_max, &
                            xgmin=tr_min_x, ygmin=tr_min_y, zgmin=tr_min_z,&
                            xgmax=tr_max_x, ygmax=tr_max_y, zgmax=tr_max_z)
    if (ntr_stocks > 0) then
      do m=1,ntr_stocks
        vars(num_nc_fields+m) = var_desc(tr_names(m), units=tr_units(m), &
                      longname=tr_names(m), hor_grid='1', z_grid='1')
      enddo
      num_nc_fields = num_nc_fields + ntr_stocks
    endif
  endif

  if (cs%previous_calls == 0) then

    cs%mass_prev_EFP = mass_efp
    cs%fresh_water_in_EFP = real_to_efp(0.0)
    if (cs%use_temperature) then
      cs%net_salt_in_EFP = real_to_efp(0.0)  ; cs%net_heat_in_EFP = real_to_efp(0.0)
    endif

    !  Reopen or create a text output file, with an explanatory header line.
    if (is_root_pe()) then
      if (day > cs%Start_time) then
        call open_file(cs%fileenergy_ascii, trim(cs%energyfile), &
                       action=append_file, form=ascii_file, nohdrs=.true.)
      else
        call open_file(cs%fileenergy_ascii, trim(cs%energyfile), &
                       action=writeonly_file, form=ascii_file, nohdrs=.true.)
        if (abs(cs%timeunit - 86400.0) < 1.0) then
          if (cs%use_temperature) then
            write(cs%fileenergy_ascii,'("  Step,",7x,"Day,  Truncs,      &
                &Energy/Mass,      Maximum CFL,  Mean Sea Level,  Total Mass,  Mean Salin, &
                &Mean Temp, Frac Mass Err,   Salin Err,    Temp Err")')
            write(cs%fileenergy_ascii,'(12x,"[days]",17x,"[m2 s-2]",11x,"[Nondim]",7x,"[m]",13x,&
                &"[kg]",9x,"[PSU]",6x,"[degC]",7x,"[Nondim]",8x,"[PSU]",8x,"[degC]")')
          else
            write(cs%fileenergy_ascii,'("  Step,",7x,"Day,  Truncs,      &
                &Energy/Mass,      Maximum CFL,  Mean sea level,   Total Mass,    Frac Mass Err")')
            write(cs%fileenergy_ascii,'(12x,"[days]",17x,"[m2 s-2]",11x,"[Nondim]",8x,"[m]",13x,&
                &"[kg]",11x,"[Nondim]")')
          endif
        else
          if ((cs%timeunit >= 0.99) .and. (cs%timeunit < 1.01)) then
            time_units = "           [seconds]     "
          elseif ((cs%timeunit >= 3599.0) .and. (cs%timeunit < 3601.0)) then
            time_units = "            [hours]      "
          elseif ((cs%timeunit >= 86399.0) .and. (cs%timeunit < 86401.0)) then
            time_units = "             [days]      "
          elseif ((cs%timeunit >= 3.0e7) .and. (cs%timeunit < 3.2e7)) then
            time_units = "            [years]      "
          else
            write(time_units,'(9x,"[",es8.2," s]    ")') cs%timeunit
          endif

          if (cs%use_temperature) then
            write(cs%fileenergy_ascii,'("  Step,",7x,"Time, Truncs,      &
                &Energy/Mass,      Maximum CFL,  Mean Sea Level,  Total Mass,  Mean Salin, &
                &Mean Temp, Frac Mass Err,   Salin Err,    Temp Err")')
            write(cs%fileenergy_ascii,'(A25,10x,"[m2 s-2]",11x,"[Nondim]",7x,"[m]",13x,&
                &"[kg]",9x,"[PSU]",6x,"[degC]",7x,"[Nondim]",8x,"[PSU]",6x,&
                &"[degC]")') time_units
          else
            write(cs%fileenergy_ascii,'("  Step,",7x,"Time, Truncs,      &
                &Energy/Mass,      Maximum CFL,  Mean sea level,   Total Mass,    Frac Mass Err")')
            write(cs%fileenergy_ascii,'(A25,10x,"[m2 s-2]",11x,"[Nondim]",8x,"[m]",13x,&
                &"[kg]",11x,"[Nondim]")') time_units
          endif
        endif
      endif
    endif

    energypath_nc = trim(cs%energyfile) // ".nc"
    if (day > cs%Start_time) then
      call reopen_file(cs%fileenergy_nc, trim(energypath_nc), vars, &
                       num_nc_fields, cs%fields, single_file, cs%timeunit, &
                       g=g, gv=gv)
    else
      call create_file(cs%fileenergy_nc, trim(energypath_nc), vars, &
                       num_nc_fields, cs%fields, single_file, cs%timeunit, &
                       g=g, gv=gv)
    endif
  endif

  if (cs%do_APE_calc) then
    lbelow = 1 ; volbelow = 0.0
    do k=nz,1,-1
      volbelow = volbelow + vol_lay(k)
      if ((volbelow >= cs%DL(cs%lH(k))%vol_below) .and. &
          (volbelow < cs%DL(cs%lH(k)+1)%vol_below)) then
        l = cs%lH(k)
      else
        labove=cs%list_size+1
        l = (labove + lbelow) / 2
        do while (l > lbelow)
          if (volbelow < cs%DL(l)%vol_below) then ; labove = l
          else ; lbelow = l ; endif
          l = (labove + lbelow) / 2
        enddo
        cs%lH(k) = l
      endif
      lbelow = l
      z_0ape(k) = cs%DL(l)%depth - (volbelow - cs%DL(l)%vol_below) / cs%DL(l)%area
    enddo
    z_0ape(nz+1) = cs%DL(2)%depth

    !   Calculate the Available Potential Energy integrated over each interface.  With a nonlinear
    ! equation of state or with a bulk mixed layer this calculation is only approximate.
    ! With an ALE model this does not make sense and should be revisited.
    pe_scale_factor = us%RZ_to_kg_m2*us%L_to_m**2*us%L_T_to_m_s**2
    pe_pt(:,:,:) = 0.0
    if (gv%Boussinesq) then
      do j=js,je ; do i=is,ie
        hbelow = 0.0
        do k=nz,1,-1
          hbelow = hbelow + h(i,j,k) * gv%H_to_Z
          hint = z_0ape(k) + (hbelow - g%bathyT(i,j))
          hbot = z_0ape(k) - g%bathyT(i,j)
          hbot = (hbot + abs(hbot)) * 0.5
          pe_pt(i,j,k) = (0.5 * pe_scale_factor * areatm(i,j)) * (gv%Rho0*gv%g_prime(k)) * &
                  (hint * hint - hbot * hbot)
        enddo
      enddo ; enddo
    else
      do j=js,je ; do i=is,ie
        do k=nz,1,-1
          hint = z_0ape(k) + eta(i,j,k)  ! eta and H_0 have opposite signs.
          hbot = max(z_0ape(k) - g%bathyT(i,j), 0.0)
          pe_pt(i,j,k) = (0.5 * pe_scale_factor * areatm(i,j) * (gv%Rho0*gv%g_prime(k))) * &
                  (hint * hint - hbot * hbot)
        enddo
      enddo ; enddo
    endif

    pe_tot = reproducing_sum(pe_pt, isr, ier, jsr, jer, sums=pe)
    do k=1,nz+1 ; h_0ape(k) = us%Z_to_m*z_0ape(k) ; enddo
  else
    pe_tot = 0.0
    do k=1,nz+1 ; pe(k) = 0.0 ; h_0ape(k) = 0.0 ; enddo
  endif

! Calculate the Kinetic Energy integrated over each layer.
  ke_scale_factor = hl2_to_kg*us%L_T_to_m_s**2
  tmp1(:,:,:) = 0.0
  do k=1,nz ; do j=js,je ; do i=is,ie
    tmp1(i,j,k) = (0.25 * ke_scale_factor * (areatm(i,j) * h(i,j,k))) * &
            ((u(i-1,j,k)**2 + u(i,j,k)**2) + (v(i,j-1,k)**2 + v(i,j,k)**2))
  enddo ; enddo ; enddo

  ke_tot = reproducing_sum(tmp1, isr, ier, jsr, jer, sums=ke)

  toten = ke_tot + pe_tot

  salt = 0.0 ; heat = 0.0
  if (cs%use_temperature) then
    temp_int(:,:) = 0.0 ; salt_int(:,:) = 0.0
    do k=1,nz ; do j=js,je ; do i=is,ie
      salt_int(i,j) = salt_int(i,j) + tv%S(i,j,k) * &
                      (h(i,j,k)*(hl2_to_kg * areatm(i,j)))
      temp_int(i,j) = temp_int(i,j) + (us%Q_to_J_kg*tv%C_p * tv%T(i,j,k)) * &
                      (h(i,j,k)*(hl2_to_kg * areatm(i,j)))
    enddo ; enddo ; enddo
    salt_efp = reproducing_sum_efp(salt_int, isr, ier, jsr, jer, only_on_pe=.true.)
    heat_efp = reproducing_sum_efp(temp_int, isr, ier, jsr, jer, only_on_pe=.true.)

    ! Combining the sums avoids multiple blocking all-PE updates.
    efp_list(1) = salt_efp ;  efp_list(2) = heat_efp ; efp_list(3) = cs%fresh_water_in_EFP
    efp_list(4) = cs%net_salt_in_EFP ; efp_list(5) = cs%net_heat_in_EFP
    call efp_sum_across_pes(efp_list, 5)
    ! Return the globally summed values to the original variables.
    salt_efp = efp_list(1) ; heat_efp = efp_list(2) ; cs%fresh_water_in_EFP = efp_list(3)
    cs%net_salt_in_EFP = efp_list(4) ; cs%net_heat_in_EFP = efp_list(5)

    salt = efp_to_real(salt_efp)
    heat = efp_to_real(heat_efp)
  else
    call efp_sum_across_pes(cs%fresh_water_in_EFP)
  endif

! Calculate the maximum CFL numbers.
  max_cfl(1:2) = 0.0
  do k=1,nz ; do j=js,je ; do i=isq,ieq
    cfl_iarea = g%IareaT(i,j)
    if (u(i,j,k) < 0.0) &
      cfl_iarea = g%IareaT(i+1,j)

    cfl_trans = abs(u(i,j,k) * cs%dt_in_T) * (g%dy_Cu(i,j) * cfl_iarea)
    cfl_lin = abs(u(i,j,k) * cs%dt_in_T) * g%IdxCu(i,j)
    max_cfl(1) = max(max_cfl(1), cfl_trans)
    max_cfl(2) = max(max_cfl(2), cfl_lin)
  enddo ; enddo ; enddo
  do k=1,nz ; do j=jsq,jeq ; do i=is,ie
    cfl_iarea = g%IareaT(i,j)
    if (v(i,j,k) < 0.0) &
      cfl_iarea = g%IareaT(i,j+1)

    cfl_trans = abs(v(i,j,k) * cs%dt_in_T) * (g%dx_Cv(i,j) * cfl_iarea)
    cfl_lin = abs(v(i,j,k) * cs%dt_in_T) * g%IdyCv(i,j)
    max_cfl(1) = max(max_cfl(1), cfl_trans)
    max_cfl(2) = max(max_cfl(2), cfl_lin)
  enddo ; enddo ; enddo

  call sum_across_pes(cs%ntrunc)
  !   Sum the various quantities across all the processors.  This sum is NOT
  ! guaranteed to be bitwise reproducible, even on the same decomposition.
  !   The sum of Tr_stocks should be reimplemented using the reproducing sums.
  if (ntr_stocks > 0) call sum_across_pes(tr_stocks,ntr_stocks)

  call max_across_pes(max_cfl, 2)
  irho0 = 1.0 / (us%R_to_kg_m3*gv%Rho0)

  if (cs%use_temperature) then
    if (cs%previous_calls == 0) then
      cs%salt_prev_EFP = salt_efp ; cs%heat_prev_EFP = heat_efp
    endif
    salt_chg_efp = salt_efp - cs%salt_prev_EFP
    salt_chg = efp_to_real(salt_chg_efp)
    salt_anom_efp = salt_chg_efp - cs%net_salt_in_EFP
    salt_anom = efp_to_real(salt_anom_efp)
    heat_chg_efp = heat_efp - cs%heat_prev_EFP
    heat_chg = efp_to_real(heat_chg_efp)
    heat_anom_efp = heat_chg_efp - cs%net_heat_in_EFP
    heat_anom = efp_to_real(heat_anom_efp)
  endif

  mass_chg_efp = mass_efp - cs%mass_prev_EFP
  mass_anom_efp = mass_chg_efp - cs%fresh_water_in_EFP
  mass_anom = efp_to_real(mass_anom_efp)
  if (cs%use_temperature .and. .not.gv%Boussinesq) then
    ! net_salt_input needs to be converted from ppt m s-1 to kg m-2 s-1.
    mass_anom = mass_anom - 0.001*efp_to_real(cs%net_salt_in_EFP)
  endif
  mass_chg = efp_to_real(mass_chg_efp)

  if (cs%use_temperature) then
    salin = salt / mass_tot ; salin_anom = salt_anom / mass_tot
   ! salin_chg = Salt_chg / mass_tot
    temp = heat / (mass_tot*us%Q_to_J_kg*tv%C_p) ; temp_anom = heat_anom / (mass_tot*us%Q_to_J_kg*tv%C_p)
  endif
  en_mass = toten / mass_tot

  call get_time(day, start_of_day, num_days)
  date_stamped = (cs%date_stamped_output .and. (get_calendar_type() /= no_calendar))
  if (date_stamped) &
    call get_date(day, iyear, imonth, iday, ihour, iminute, isecond, itick)
  if (abs(cs%timeunit - 86400.0) < 1.0) then
    reday = REAL(num_days)+ (REAL(start_of_day)/86400.0)
    mesg_intro = "MOM Day "
  else
    reday = REAL(num_days)*(86400.0/cs%timeunit) + &
            REAL(start_of_day)/abs(cs%timeunit)
    mesg_intro = "MOM Time "
  endif
  if (reday < 1.0e8) then ;      write(day_str, '(F12.3)') reday
  elseif (reday < 1.0e11) then ; write(day_str, '(F15.3)') reday
  else ;                         write(day_str, '(ES15.9)') reday ; endif

  if     (n < 1000000)   then ; write(n_str, '(I6)')  n
  elseif (n < 10000000)  then ; write(n_str, '(I7)')  n
  elseif (n < 100000000) then ; write(n_str, '(I8)')  n
  else                        ; write(n_str, '(I10)') n ; endif

  if (date_stamped) then
    write(date_str,'("MOM Date",i7,2("/",i2.2)," ",i2.2,2(":",i2.2))') &
       iyear, imonth, iday, ihour, iminute, isecond
  else
    date_str = trim(mesg_intro)//trim(day_str)
  endif

  if (is_root_pe()) then
    if (cs%use_temperature) then
        write(*,'(A," ",A,": En ",ES12.6, ", MaxCFL ", F8.5, ", Mass ", &
                & ES18.12, ", Salt ", F15.11,", Temp ", F15.11)') &
          trim(date_str), trim(n_str), en_mass, max_cfl(1), mass_tot, salin, temp
    else
        write(*,'(A," ",A,": En ",ES12.6, ", MaxCFL ", F8.5, ", Mass ", &
                & ES18.12)') &
          trim(date_str), trim(n_str), en_mass, max_cfl(1), mass_tot
    endif

    if (cs%use_temperature) then
      write(cs%fileenergy_ascii,'(A,",",A,",", I6,", En ",ES22.16, &
                               &", CFL ", F8.5, ", SL ",&
                               &es11.4,", M ",ES11.5,", S",f8.4,", T",f8.4,&
                               &", Me ",ES9.2,", Se ",ES9.2,", Te ",ES9.2)') &
            trim(n_str), trim(day_str), cs%ntrunc, en_mass, max_cfl(1), &
            -h_0ape(1), mass_tot, salin, temp, mass_anom/mass_tot, salin_anom, &
            temp_anom
    else
      write(cs%fileenergy_ascii,'(A,",",A,",", I6,", En ",ES22.16, &
                                &", CFL ", F8.5, ", SL ",&
                                  &ES11.4,", Mass ",ES11.5,", Me ",ES9.2)') &
            trim(n_str), trim(day_str), cs%ntrunc, en_mass, max_cfl(1), &
            -h_0ape(1), mass_tot, mass_anom/mass_tot
    endif

    if (cs%ntrunc > 0) then
      write(*,'(A," Energy/Mass:",ES12.5," Truncations ",I0)') &
        trim(date_str), en_mass, cs%ntrunc
    endif

    if (cs%write_stocks) then
      write(*,'("    Total Energy: ",Z16.16,ES24.16)') toten, toten
      write(*,'("    Total Mass: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') &
            mass_tot, mass_chg, mass_anom, mass_anom/mass_tot
      if (cs%use_temperature) then
        if (salt == 0.) then
          write(*,'("    Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') &
              salt*0.001, salt_chg*0.001, salt_anom*0.001
        else
          write(*,'("    Total Salt: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') &
              salt*0.001, salt_chg*0.001, salt_anom*0.001, salt_anom/salt
        endif
        if (heat == 0.) then
          write(*,'("    Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5)') &
              heat, heat_chg, heat_anom
        else
          write(*,'("    Total Heat: ",ES24.16,", Change: ",ES24.16," Error: ",ES12.5," (",ES8.1,")")') &
              heat, heat_chg, heat_anom, heat_anom/heat
        endif
      endif
      do m=1,ntr_stocks

         write(*,'("      Total ",a,": ",ES24.16,X,a)') &
              trim(tr_names(m)), tr_stocks(m), trim(tr_units(m))

         if (tr_minmax_got(m)) then
           write(*,'(64X,"Global Min:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")"  )') &
                tr_min(m),tr_min_x(m),tr_min_y(m),tr_min_z(m)
           write(*,'(64X,"Global Max:",ES24.16,X,"at: (", f7.2,","f7.2,","f8.2,")"  )') &
                tr_max(m),tr_max_x(m),tr_max_y(m),tr_max_z(m)
        endif

      enddo
    endif
  endif

  var = real(cs%ntrunc)
  call write_field(cs%fileenergy_nc, cs%fields(1), var, reday)
  call write_field(cs%fileenergy_nc, cs%fields(2), toten, reday)
  call write_field(cs%fileenergy_nc, cs%fields(3), pe, reday)
  call write_field(cs%fileenergy_nc, cs%fields(4), ke, reday)
  call write_field(cs%fileenergy_nc, cs%fields(5), h_0ape, reday)
  call write_field(cs%fileenergy_nc, cs%fields(6), mass_lay, reday)

  call write_field(cs%fileenergy_nc, cs%fields(7), mass_tot, reday)
  call write_field(cs%fileenergy_nc, cs%fields(8), mass_chg, reday)
  call write_field(cs%fileenergy_nc, cs%fields(9), mass_anom, reday)
  call write_field(cs%fileenergy_nc, cs%fields(10), max_cfl(1), reday)
  call write_field(cs%fileenergy_nc, cs%fields(11), max_cfl(2), reday)
  if (cs%use_temperature) then
    call write_field(cs%fileenergy_nc, cs%fields(12), 0.001*salt, reday)
    call write_field(cs%fileenergy_nc, cs%fields(13), 0.001*salt_chg, reday)
    call write_field(cs%fileenergy_nc, cs%fields(14), 0.001*salt_anom, reday)
    call write_field(cs%fileenergy_nc, cs%fields(15), heat, reday)
    call write_field(cs%fileenergy_nc, cs%fields(16), heat_chg, reday)
    call write_field(cs%fileenergy_nc, cs%fields(17), heat_anom, reday)
    do m=1,ntr_stocks
      call write_field(cs%fileenergy_nc, cs%fields(17+m), tr_stocks(m), reday)
    enddo
  else
    do m=1,ntr_stocks
      call write_field(cs%fileenergy_nc, cs%fields(11+m), tr_stocks(m), reday)
    enddo
  endif

  call flush_file(cs%fileenergy_nc)

  ! The second (impossible-looking) test looks for a NaN in En_mass.
  if ((en_mass>cs%max_Energy) .or. &
     ((en_mass>cs%max_Energy) .and. (en_mass<cs%max_Energy))) then
    write(mesg,'("Energy per unit mass of ",ES11.4," exceeds ",ES11.4)') &
                  en_mass, cs%max_Energy
    call mom_error(fatal, &
      "write_energy : Excessive energy per unit mass or NaNs forced model termination.")
  endif
  if (cs%ntrunc>cs%maxtrunc) then
    call mom_error(fatal, "write_energy : Ocean velocity has been truncated too many times.")
  endif
  cs%ntrunc = 0
  cs%previous_calls = cs%previous_calls + 1

  cs%mass_prev_EFP = mass_efp ; cs%fresh_water_in_EFP = real_to_efp(0.0)
  if (cs%use_temperature) then
    cs%salt_prev_EFP = salt_efp ; cs%net_salt_in_EFP = real_to_efp(0.0)
    cs%heat_prev_EFP = heat_efp ; cs%net_heat_in_EFP = real_to_efp(0.0)
  endif

end subroutine write_energy

!> This subroutine accumates the net input of volume, salt and heat, through
!! the ocean surface for use in diagnosing conservation.
subroutine accumulate_net_input(fluxes, sfc_state, tv, dt, G, US, CS)
  type(forcing),         intent(in) :: fluxes !< A structure containing pointers to any possible
                                              !! forcing fields.  Unused fields are unallocated.
  type(surface),         intent(in) :: sfc_state !< A structure containing fields that
                                              !! describe the surface state of the ocean.
  type(thermo_var_ptrs), intent(in) :: tv     !< A structure pointing to various
                                              !! thermodynamic variables.
  real,                  intent(in) :: dt     !< The amount of time over which to average [T ~> s].
  type(ocean_grid_type), intent(in) :: g      !< The ocean's grid structure.
  type(unit_scale_type), intent(in) :: us     !< A dimensional unit scaling type
  type(sum_output_cs),   pointer    :: cs     !< The control structure returned by a previous call
                                              !! to MOM_sum_output_init.
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G)) :: &
    fw_in, &   ! The net fresh water input, integrated over a timestep [kg].
    salt_in, & ! The total salt added by surface fluxes, integrated
               ! over a time step [ppt kg].
    heat_in    ! The total heat added by surface fluxes, integrated
               ! over a time step [J].
  real :: fw_input   ! The net fresh water input, integrated over a timestep
                     ! and summed over space [kg].
  real :: salt_input ! The total salt added by surface fluxes, integrated
                     ! over a time step and summed over space [ppt kg].
  real :: heat_input ! The total heat added by boundary fluxes, integrated
                     ! over a time step and summed over space [J].
  real :: rzl2_to_kg ! A combination of scaling factors for mass [kg R-1 Z-1 L-2 ~> 1]
  real :: qrzl2_to_j ! A combination of scaling factors for heat [J Q-1 R-1 Z-1 L-2 ~> 1]

  type(efp_type) :: &
    fw_in_efp,   &   ! The net fresh water input, integrated over a timestep
                     ! and summed over space [kg].
    salt_in_efp, &   ! The total salt added by surface fluxes, integrated
                     ! over a time step and summed over space [ppt kg].
    heat_in_efp      ! The total heat added by boundary fluxes, integrated
                     ! over a time step and summed over space [J].

  real :: inputs(3)   ! A mixed array for combining the sums
  integer :: i, j, is, ie, js, je, isr, ier, jsr, jer

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

  rzl2_to_kg = us%L_to_m**2*us%RZ_to_kg_m2
  qrzl2_to_j = rzl2_to_kg*us%Q_to_J_kg

  fw_in(:,:) = 0.0
  if (associated(fluxes%evap)) then
    if (associated(fluxes%lprec) .and. associated(fluxes%fprec)) then
      do j=js,je ; do i=is,ie
        fw_in(i,j) = rzl2_to_kg * dt*g%areaT(i,j)*(fluxes%evap(i,j) + &
            (((fluxes%lprec(i,j) + fluxes%vprec(i,j)) + fluxes%lrunoff(i,j)) + &
              (fluxes%fprec(i,j) + fluxes%frunoff(i,j))))
      enddo ; enddo
    else
      call mom_error(warning, &
        "accumulate_net_input called with associated evap field, but no precip field.")
    endif
  endif

  if (associated(fluxes%seaice_melt)) then ; do j=js,je ; do i=is,ie
    fw_in(i,j) = fw_in(i,j) + rzl2_to_kg*dt * &
                 g%areaT(i,j) * fluxes%seaice_melt(i,j)
  enddo ; enddo ; endif

  salt_in(:,:) = 0.0 ; heat_in(:,:) = 0.0
  if (cs%use_temperature) then

    if (associated(fluxes%sw)) then ; do j=js,je ; do i=is,ie
      heat_in(i,j) = heat_in(i,j) + dt * qrzl2_to_j * g%areaT(i,j) * (fluxes%sw(i,j) + &
             (fluxes%lw(i,j) + (fluxes%latent(i,j) + fluxes%sens(i,j))))
    enddo ; enddo ; endif

    if (associated(fluxes%seaice_melt_heat)) then ; do j=js,je ; do i=is,ie
       heat_in(i,j) = heat_in(i,j) + dt * qrzl2_to_j * g%areaT(i,j) * &
                                     fluxes%seaice_melt_heat(i,j)
    enddo ; enddo ; endif

    ! smg: new code
    ! include heat content from water transport across ocean surface
!    if (associated(fluxes%heat_content_lprec)) then ; do j=js,je ; do i=is,ie
!      heat_in(i,j) = heat_in(i,j) + dt * QRZL2_to_J * G%areaT(i,j) * &
!         (fluxes%heat_content_lprec(i,j)   + (fluxes%heat_content_fprec(i,j)   &
!       + (fluxes%heat_content_lrunoff(i,j) + (fluxes%heat_content_frunoff(i,j) &
!       + (fluxes%heat_content_cond(i,j)    + (fluxes%heat_content_vprec(i,j)   &
!       +  fluxes%heat_content_massout(i,j)))))))
!    enddo ; enddo ; endif

    ! smg: old code
    if (associated(tv%TempxPmE)) then
      do j=js,je ; do i=is,ie
        heat_in(i,j) = heat_in(i,j) + (fluxes%C_p * qrzl2_to_j*g%areaT(i,j)) * tv%TempxPmE(i,j)
      enddo ; enddo
    elseif (associated(fluxes%evap)) then
      do j=js,je ; do i=is,ie
        heat_in(i,j) = heat_in(i,j) + (us%Q_to_J_kg*fluxes%C_p * sfc_state%SST(i,j)) * fw_in(i,j)
      enddo ; enddo
    endif

    ! The following heat sources may or may not be used.
    if (associated(tv%internal_heat)) then
      do j=js,je ; do i=is,ie
        heat_in(i,j) = heat_in(i,j) + (fluxes%C_p * qrzl2_to_j*g%areaT(i,j)) * tv%internal_heat(i,j)
      enddo ; enddo
    endif
    if (associated(tv%frazil)) then ; do j=js,je ; do i=is,ie
      heat_in(i,j) = heat_in(i,j) + qrzl2_to_j * g%areaT(i,j) * tv%frazil(i,j)
    enddo ; enddo ; endif
    if (associated(fluxes%heat_added)) then ; do j=js,je ; do i=is,ie
      heat_in(i,j) = heat_in(i,j) + qrzl2_to_j * dt*g%areaT(i,j) * fluxes%heat_added(i,j)
    enddo ; enddo ; endif
!    if (associated(sfc_state%sw_lost)) then ; do j=js,je ; do i=is,ie
!      heat_in(i,j) = heat_in(i,j) - US%L_to_m**2*G%areaT(i,j) * sfc_state%sw_lost(i,j)
!    enddo ; enddo ; endif

    if (associated(fluxes%salt_flux)) then ; do j=js,je ; do i=is,ie
      ! convert salt_flux from kg (salt)/(m^2 s) to ppt * [m s-1].
      salt_in(i,j) = rzl2_to_kg * dt * &
                     g%areaT(i,j)*(1000.0*fluxes%salt_flux(i,j))
    enddo ; enddo ; endif
  endif

  if ((cs%use_temperature) .or. associated(fluxes%lprec) .or. &
      associated(fluxes%evap)) then
    ! The on-PE sums are stored here, but the sums across PEs are deferred to
    ! the next call to write_energy to avoid extra barriers.
    isr = is - (g%isd-1) ; ier = ie - (g%isd-1) ; jsr = js - (g%jsd-1) ; jer = je - (g%jsd-1)
    fw_in_efp   = reproducing_sum_efp(fw_in,   isr, ier, jsr, jer, only_on_pe=.true.)
    heat_in_efp = reproducing_sum_efp(heat_in, isr, ier, jsr, jer, only_on_pe=.true.)
    salt_in_efp = reproducing_sum_efp(salt_in, isr, ier, jsr, jer, only_on_pe=.true.)

    cs%fresh_water_in_EFP = cs%fresh_water_in_EFP + fw_in_efp
    cs%net_salt_in_EFP    = cs%net_salt_in_EFP    + salt_in_efp
    cs%net_heat_in_EFP    = cs%net_heat_in_EFP    + heat_in_efp
  endif

end subroutine accumulate_net_input

!>  This subroutine sets up an ordered list of depths, along with the
!! cross sectional areas at each depth and the volume of fluid deeper
!! than each depth.  This might be read from a previously created file
!! or it might be created anew.  (For now only new creation occurs.
subroutine depth_list_setup(G, US, CS)
  type(ocean_grid_type), intent(in) :: G   !< The ocean's grid structure
  type(unit_scale_type), intent(in) :: US  !< A dimensional unit scaling type
  type(Sum_output_CS),   pointer    :: CS  !< The control structure returned by a
                                           !! previous call to MOM_sum_output_init.
  ! Local variables
  integer :: k

  if (cs%read_depth_list) then
    if (file_exists(cs%depth_list_file)) then
      call read_depth_list(g, us, cs, cs%depth_list_file)
    else
      if (is_root_pe()) call mom_error(warning, "depth_list_setup: "// &
        trim(cs%depth_list_file)//" does not exist.  Creating a new file.")
      call create_depth_list(g, cs)

      call write_depth_list(g, us, cs, cs%depth_list_file, cs%list_size+1)
    endif
  else
    call create_depth_list(g, cs)
  endif

  do k=1,g%ke
    cs%lH(k) = cs%list_size
  enddo

end subroutine depth_list_setup

!>  create_depth_list makes an ordered list of depths, along with the cross
!! sectional areas at each depth and the volume of fluid deeper than each depth.
subroutine create_depth_list(G, CS)
  type(ocean_grid_type), intent(in) :: G  !< The ocean's grid structure.
  type(Sum_output_CS),   pointer    :: CS !< The control structure set up in MOM_sum_output_init,
                                          !! in which the ordered depth list is stored.
  ! Local variables
  real, dimension(G%Domain%niglobal*G%Domain%njglobal + 1) :: &
    Dlist, &  !< The global list of bottom depths [Z ~> m].
    AreaList  !< The global list of cell areas [L2 ~> m2].
  integer, dimension(G%Domain%niglobal*G%Domain%njglobal+1) :: &
    indx2     !< The position of an element in the original unsorted list.
  real    :: Dnow  !< The depth now being considered for sorting [Z ~> m].
  real    :: Dprev !< The most recent depth that was considered [Z ~> m].
  real    :: vol   !< The running sum of open volume below a deptn [Z L2 ~> m3].
  real    :: area  !< The open area at the current depth [L2 ~> m2].
  real    :: D_list_prev !< The most recent depth added to the list [Z ~> m].
  logical :: add_to_list !< This depth should be included as an entry on the list.

  integer :: ir, indxt
  integer :: mls, list_size
  integer :: list_pos, i_global, j_global
  integer :: i, j, k, kl

  mls = g%Domain%niglobal*g%Domain%njglobal

! Need to collect the global data from compute domains to a 1D array for sorting.
  dlist(:) = 0.0
  arealist(:) = 0.0
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    ! Set global indices that start the global domain at 1 (Fortran convention).
    j_global = j + g%jdg_offset - (g%jsg-1)
    i_global = i + g%idg_offset - (g%isg-1)

    list_pos = (j_global-1)*g%Domain%niglobal + i_global
    dlist(list_pos) = g%bathyT(i,j)
    arealist(list_pos) = g%mask2dT(i,j) * g%areaT(i,j)
  enddo ; enddo

  ! These sums reproduce across PEs because the arrays are only nonzero on one PE.
  call sum_across_pes(dlist, mls+1)
  call sum_across_pes(arealist, mls+1)

  do j=1,mls+1 ; indx2(j) = j ; enddo
  k = mls / 2  + 1 ; ir = mls
  do
    if (k > 1) then
      k = k - 1
      indxt = indx2(k)
      dnow = dlist(indxt)
    else
      indxt = indx2(ir)
      dnow = dlist(indxt)
      indx2(ir) = indx2(1)
      ir = ir - 1
      if (ir == 1) then ; indx2(1) = indxt ; exit ; endif
    endif
    i=k ; j=k*2
    do ; if (j > ir) exit
      if (j < ir .AND. dlist(indx2(j)) < dlist(indx2(j+1))) j = j + 1
      if (dnow < dlist(indx2(j))) then ; indx2(i) = indx2(j) ; i = j ; j = j + i
      else ; j = ir+1 ; endif
    enddo
    indx2(i) = indxt
  enddo

!  At this point, the lists should perhaps be culled to save memory.
! Culling: (1) identical depths (e.g. land) - take the last one.
!          (2) the topmost and bottommost depths are always saved.
!          (3) very close depths
!          (4) equal volume changes.

  ! Count the unique elements in the list.
  d_list_prev = dlist(indx2(mls))
  list_size = 2
  do k=mls-1,1,-1
    if (dlist(indx2(k)) < d_list_prev-cs%D_list_min_inc) then
      list_size = list_size + 1
      d_list_prev = dlist(indx2(k))
    endif
  enddo

  cs%list_size = list_size
  allocate(cs%DL(cs%list_size+1))

  vol = 0.0 ; area = 0.0
  dprev = dlist(indx2(mls))
  d_list_prev = dprev

  kl = 0
  do k=mls,1,-1
    i = indx2(k)
    vol = vol + area * (dprev - dlist(i))
    area = area + arealist(i)

    add_to_list = .false.
    if ((kl == 0) .or. (k==1)) then
      add_to_list = .true.
    elseif (dlist(indx2(k-1)) < d_list_prev-cs%D_list_min_inc) then
      add_to_list = .true.
      d_list_prev = dlist(indx2(k-1))
    endif

    if (add_to_list) then
      kl = kl+1
      cs%DL(kl)%depth = dlist(i)
      cs%DL(kl)%area = area
      cs%DL(kl)%vol_below = vol
    endif
    dprev = dlist(i)
  enddo

  do while (kl < list_size)
    ! I don't understand why this is needed... RWH
    kl = kl+1
    cs%DL(kl)%vol_below = cs%DL(kl-1)%vol_below * 1.000001
    cs%DL(kl)%area = cs%DL(kl-1)%area
    cs%DL(kl)%depth = cs%DL(kl-1)%depth
  enddo

  cs%DL(cs%list_size+1)%vol_below = cs%DL(cs%list_size)%vol_below * 1000.0
  cs%DL(cs%list_size+1)%area = cs%DL(cs%list_size)%area
  cs%DL(cs%list_size+1)%depth = cs%DL(cs%list_size)%depth

end subroutine create_depth_list

!> This subroutine writes out the depth list to the specified file.
subroutine write_depth_list(G, US, CS, filename, list_size)
  type(ocean_grid_type), intent(in) :: G   !< The ocean's grid structure.
  type(unit_scale_type), intent(in) :: US  !< A dimensional unit scaling type
  type(Sum_output_CS),   pointer    :: CS  !< The control structure returned by a
                                           !! previous call to MOM_sum_output_init.
  character(len=*),      intent(in) :: filename !< The path to the depth list file to write.
  integer,               intent(in) :: list_size !< The size of the depth list.
  ! Local variables
  real, allocatable :: tmp(:)
  integer :: ncid, dimid(1), Did, Aid, Vid, status, k
  character(len=16) :: depth_chksum, area_chksum

  ! All ranks are required to compute the global checksum
  call get_depth_list_checksums(g, depth_chksum, area_chksum)

  if (.not.is_root_pe()) return

  allocate(tmp(list_size)) ; tmp(:) = 0.0

  status = nf90_create(filename, 0, ncid)
  if (status /= nf90_noerr) then
    call mom_error(warning, filename//trim(nf90_strerror(status)))
    return
  endif

  status = nf90_def_dim(ncid, "list", list_size, dimid(1))
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//trim(nf90_strerror(status)))

  status = nf90_def_var(ncid, "depth", nf90_double, dimid, did)
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" depth "//trim(nf90_strerror(status)))
  status = nf90_put_att(ncid, did, "long_name", "Sorted depth")
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" depth "//trim(nf90_strerror(status)))
  status = nf90_put_att(ncid, did, "units", "m")
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" depth "//trim(nf90_strerror(status)))

  status = nf90_def_var(ncid, "area", nf90_double, dimid, aid)
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" area "//trim(nf90_strerror(status)))
  status = nf90_put_att(ncid, aid, "long_name", "Open area at depth")
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" area "//trim(nf90_strerror(status)))
  status = nf90_put_att(ncid, aid, "units", "m2")
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" area "//trim(nf90_strerror(status)))

  status = nf90_def_var(ncid, "vol_below", nf90_double, dimid, vid)
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" vol_below "//trim(nf90_strerror(status)))
  status = nf90_put_att(ncid, vid, "long_name", "Open volume below depth")
   if (status /= nf90_noerr) call mom_error(warning, &
      filename//" vol_below "//trim(nf90_strerror(status)))
  status = nf90_put_att(ncid, vid, "units", "m3")
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" vol_below "//trim(nf90_strerror(status)))

  ! Dependency checksums
  status = nf90_put_att(ncid, nf90_global, depth_chksum_attr, depth_chksum)
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" "//depth_chksum_attr//" "//trim(nf90_strerror(status)))

  status = nf90_put_att(ncid, nf90_global, area_chksum_attr, area_chksum)
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" "//area_chksum_attr//" "//trim(nf90_strerror(status)))

  status = nf90_enddef(ncid)
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//trim(nf90_strerror(status)))

  do k=1,list_size ; tmp(k) = us%Z_to_m*cs%DL(k)%depth ; enddo
  status = nf90_put_var(ncid, did, tmp)
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" depth "//trim(nf90_strerror(status)))

  do k=1,list_size ; tmp(k) = us%L_to_m**2*cs%DL(k)%area ; enddo
  status = nf90_put_var(ncid, aid, tmp)
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" area "//trim(nf90_strerror(status)))

  do k=1,list_size ; tmp(k) = us%Z_to_m*us%L_to_m**2*cs%DL(k)%vol_below ; enddo
  status = nf90_put_var(ncid, vid, tmp)
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//" vol_below "//trim(nf90_strerror(status)))

  status = nf90_close(ncid)
  if (status /= nf90_noerr) call mom_error(warning, &
      filename//trim(nf90_strerror(status)))

end subroutine write_depth_list

!> This subroutine reads in the depth list to the specified file
!! and allocates and sets up CS%DL and CS%list_size .
subroutine read_depth_list(G, US, CS, filename)
  type(ocean_grid_type), intent(in) :: G   !< The ocean's grid structure
  type(unit_scale_type), intent(in) :: US  !< A dimensional unit scaling type
  type(Sum_output_CS),   pointer    :: CS  !< The control structure returned by a
                                           !! previous call to MOM_sum_output_init.
  character(len=*),      intent(in) :: filename !< The path to the depth list file to read.
  ! Local variables
  character(len=32) :: mdl
  character(len=240) :: var_name, var_msg
  real, allocatable :: tmp(:)
  integer :: ncid, status, varid, list_size, k
  integer :: ndim, len, var_dim_ids(NF90_MAX_VAR_DIMS)
  character(len=16) :: depth_file_chksum, depth_grid_chksum
  character(len=16) :: area_file_chksum, area_grid_chksum
  integer :: depth_attr_status, area_attr_status

  mdl = "MOM_sum_output read_depth_list:"

  status = nf90_open(filename, nf90_nowrite, ncid)
  if (status /= nf90_noerr) then
    call mom_error(fatal,mdl//" Difficulties opening "//trim(filename)// &
        " - "//trim(nf90_strerror(status)))
  endif

  ! Check bathymetric consistency
  depth_attr_status = nf90_get_att(ncid, nf90_global, depth_chksum_attr, &
                                   depth_file_chksum)
  area_attr_status = nf90_get_att(ncid, nf90_global, area_chksum_attr, &
                                  area_file_chksum)

  if (any([depth_attr_status, area_attr_status] == nf90_enotatt)) then
    var_msg = trim(cs%depth_list_file) // " checksums are missing;"
    if (cs%require_depth_list_chksum) then
      call mom_error(fatal, trim(var_msg) // " aborting.")
    elseif (cs%update_depth_list_chksum) then
      call mom_error(warning, trim(var_msg) // " updating file.")
      call create_depth_list(g, cs)
      call write_depth_list(g, us, cs, cs%depth_list_file, cs%list_size+1)
      return
    else
      call mom_error(warning, &
        trim(var_msg) // " some diagnostics may not be reproducible.")
    endif
  else
    ! Validate netCDF call
    if (depth_attr_status /= nf90_noerr) then
      var_msg = mdl // "Failed to read " // trim(filename) // ":" &
                // depth_chksum_attr
      call mom_error(fatal, &
        trim(var_msg) // " - " // nf90_strerror(depth_attr_status))
    endif

    if (area_attr_status /= nf90_noerr) then
      var_msg = mdl // "Failed to read " // trim(filename) // ":" &
                // area_chksum_attr
      call mom_error(fatal, &
        trim(var_msg) // " - " // nf90_strerror(area_attr_status))
    endif

    call get_depth_list_checksums(g, depth_grid_chksum, area_grid_chksum)

    if (depth_grid_chksum /= depth_file_chksum &
            .or. area_grid_chksum /= area_file_chksum) then
      var_msg = trim(cs%depth_list_file) // " checksums do not match;"
      if (cs%require_depth_list_chksum) then
        call mom_error(fatal, trim(var_msg) // " aborting.")
      elseif (cs%update_depth_list_chksum) then
        call mom_error(warning, trim(var_msg) // " updating file.")
        call create_depth_list(g, cs)
        call write_depth_list(g, us, cs, cs%depth_list_file, cs%list_size+1)
        return
      else
        call mom_error(warning, &
          trim(var_msg) // " some diagnostics may not be reproducible.")
      endif
    endif
  endif

  var_name = "depth"
  var_msg = trim(var_name)//" in "//trim(filename)//" - "
  status = nf90_inq_varid(ncid, var_name, varid)
  if (status /= nf90_noerr) call mom_error(fatal,mdl// &
        " Difficulties finding variable "//trim(var_msg)//&
        trim(nf90_strerror(status)))

  status = nf90_inquire_variable(ncid, varid, ndims=ndim, dimids=var_dim_ids)
  if (status /= nf90_noerr) then
    call mom_error(fatal,mdl//" cannot inquire about "//trim(var_msg)//&
        trim(nf90_strerror(status)))
  elseif (ndim > 1) then
    call mom_error(fatal,mdl//" "//trim(var_msg)//&
         " has too many or too few dimensions.")
  endif

  ! Get the length of the list.
  status = nf90_inquire_dimension(ncid, var_dim_ids(1), len=list_size)
  if (status /= nf90_noerr) call mom_error(fatal,mdl// &
        " cannot inquire about dimension(1) of "//trim(var_msg)//&
        trim(nf90_strerror(status)))

  cs%list_size = list_size-1
  allocate(cs%DL(list_size))
  allocate(tmp(list_size))

  status = nf90_get_var(ncid, varid, tmp)
  if (status /= nf90_noerr) call mom_error(fatal,mdl// &
        " Difficulties reading variable "//trim(var_msg)//&
        trim(nf90_strerror(status)))

  do k=1,list_size ; cs%DL(k)%depth = us%m_to_Z*tmp(k) ; enddo

  var_name = "area"
  var_msg = trim(var_name)//" in "//trim(filename)//" - "
  status = nf90_inq_varid(ncid, var_name, varid)
  if (status /= nf90_noerr) call mom_error(fatal,mdl// &
        " Difficulties finding variable "//trim(var_msg)//&
        trim(nf90_strerror(status)))
  status = nf90_get_var(ncid, varid, tmp)
  if (status /= nf90_noerr) call mom_error(fatal,mdl// &
        " Difficulties reading variable "//trim(var_msg)//&
        trim(nf90_strerror(status)))

  do k=1,list_size ; cs%DL(k)%area = us%m_to_L**2*tmp(k) ; enddo

  var_name = "vol_below"
  var_msg = trim(var_name)//" in "//trim(filename)
  status = nf90_inq_varid(ncid, var_name, varid)
  if (status /= nf90_noerr) call mom_error(fatal,mdl// &
        " Difficulties finding variable "//trim(var_msg)//&
        trim(nf90_strerror(status)))
  status = nf90_get_var(ncid, varid, tmp)
  if (status /= nf90_noerr) call mom_error(fatal,mdl// &
        " Difficulties reading variable "//trim(var_msg)//&
        trim(nf90_strerror(status)))

  do k=1,list_size ; cs%DL(k)%vol_below = us%m_to_Z*us%m_to_L**2*tmp(k) ; enddo

  status = nf90_close(ncid)
  if (status /= nf90_noerr) call mom_error(warning, mdl// &
    " Difficulties closing "//trim(filename)//" - "//trim(nf90_strerror(status)))

  deallocate(tmp)

end subroutine read_depth_list


!> Return the checksums required to verify DEPTH_LIST_FILE contents.
!!
!! This function computes checksums for the bathymetry (G%bathyT) and masked
!! area (mask2dT * areaT) fields of the model grid G, which are used to compute
!! the depth list.  A difference in checksum indicates that a different method
!! was used to compute the grid data, and that any results using the depth
!! list, such as APE, will not be reproducible.
!!
!! Checksums are saved as hexadecimal strings, in order to avoid potential
!! datatype issues with netCDF attributes.
subroutine get_depth_list_checksums(G, depth_chksum, area_chksum)
  type(ocean_grid_type), intent(in) :: G          !< Ocean grid structure
  character(len=16), intent(out) :: depth_chksum  !< Depth checksum hexstring
  character(len=16), intent(out) :: area_chksum   !< Area checksum hexstring

  integer :: i, j
  real, allocatable :: field(:,:)

  allocate(field(g%isc:g%iec, g%jsc:g%jec))

  ! Depth checksum
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    field(i,j) = g%bathyT(i,j)
  enddo ; enddo
  write(depth_chksum, '(Z16)') mpp_chksum(field(:,:))

  ! Area checksum
  do j=g%jsc,g%jec ; do i=g%isc,g%iec
    field(i,j) = g%mask2dT(i,j) * g%US%L_to_m**2*g%areaT(i,j)
  enddo ; enddo
  write(area_chksum, '(Z16)') mpp_chksum(field(:,:))

  deallocate(field)
end subroutine get_depth_list_checksums

!> \namespace mom_sum_output
!!
!! By Robert Hallberg, April 1994 - June 2002
!!
!!   This file contains the subroutine (write_energy) that writes
!! horizontally integrated quantities, such as energies and layer
!! volumes, and other summary information to an output file.  Some
!! of these quantities (APE or resting interface height) are defined
!! relative to the global histogram of topography.  The subroutine
!! that compiles that histogram (depth_list_setup) is also included
!! in this file.
!!
!!   In addition, if the number of velocity truncations since the
!! previous call to write_energy exceeds maxtrunc or the total energy
!! exceeds a very large threshold, a fatal termination is triggered.

end module mom_sum_output