MOM_sum_output.F90

Go to the documentation of this file.

module mom_sum_output
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  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, the day is increased to Huge_time  *
!*  so that the model will gracefully halt itself.                     *
!*                                                                     *
!*    This file also contains a few miscelaneous initialization        *
!*  calls to FMS-related modules.                                      *
!*                                                                     *
!*  Macros written all in capital letters are defined in MOM_memory.h. *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q                                        *
!*    j+1  > o > o >   At ^:  v                                        *
!*    j    x ^ x ^ x   At >:  u                                        *
!*    j    > o > o >   At o:  h, bathyT                                *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use mom_coms, only : sum_across_pes, pe_here, root_pe, num_pes, max_across_pes
use mom_coms, only : reproducing_sum
use mom_coms, only : efp_type, operator(+), operator(-), assignment(=), efp_to_real, real_to_efp
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, get_filename_appendix
use mom_io, only : file_exists, slasher, vardesc, var_desc, write_field
use mom_io, only : append_file, ascii_file, single_file, writeonly_file
use mom_time_manager, only : time_type, get_time, get_date, set_time, operator(>), operator(-)
use mom_time_manager, only : get_calendar_type, no_calendar
use mom_tracer_flow_control, only : tracer_flow_control_cs, call_tracer_stocks
use mom_variables, only : surface, thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type

use netcdf

implicit none ; private

#include <MOM_memory.h>

public write_energy, accumulate_net_input, mom_sum_output_init

!-----------------------------------------------------------------------

integer, parameter :: num_fields = 17

type :: depth_list
  real :: depth       ! A depth, in m.
  real :: area        ! The cross-sectional area of the ocean at that depth, in m2.
  real :: vol_below   ! The ocean volume below that depth, in m3.
end type depth_list

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

  integer allocable_, dimension(NKMEM_) :: 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, in m, between the
                                ! depths of the entries in the depth-list file,
                                ! 0 by default.
  logical :: use_temperature    !   If true, temperature and salinity are state
                                ! variables.
  real    :: fresh_water_input  !   The total mass of fresh water added by
                                ! surface fluxes since the last time that
  real    :: mass_prev          !   The total ocean mass the last time that
                                ! write_energy was called, in kg.
  real    :: salt_prev          !   The total amount of salt in the ocean the last
                                ! time that write_energy was called, in PSU kg.
  real    :: net_salt_input     !   The total salt added by surface fluxes since
                                ! the last time that write_energy was called,
                                ! in PSU kg.
  real    :: heat_prev          !   The total amount of heat in the ocean the last
                                ! time that write_energy was called, in Joules.
  real    :: net_heat_input     !   The total heat added by surface fluxes since
                                ! the last time that write_energy was called,
                                ! in Joules.
  type(efp_type) :: &
    fresh_water_in_efp, &       ! These are extended fixed point versions of the
    net_salt_in_efp, &          ! correspondingly named variables above.
    net_heat_in_efp, heat_prev_efp, salt_prev_efp, mass_prev_efp
  real    :: dt                 ! The baroclinic dynamics time step, in s.
  real    :: timeunit           !   The length of the units for the time
                                ! axis, in 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
  type(time_type) :: huge_time  ! A large time, which is used to indicate
                                ! that an error has been encountered
                                ! and the run should be terminated with
                                ! an error code.
  integer, pointer :: ntrunc    ! 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, in 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

! #@# This subroutine needs a doxygen description.
subroutine mom_sum_output_init(G, param_file, directory, ntrnc, &
                                Input_start_time, CS)
  type(ocean_grid_type),  intent(in)    :: G          !< The ocean's grid structure.
  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.
! Arguments: G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in)      directory - The directory where the energy file goes.
!  (in/out)  ntrnc - The integer that stores the number of times the velocity
!                     has been truncated since the last call to write_energy.
!  (in)      Input_start_time - The start time of the simulation.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                 for this module
  real :: Rho_0, maxvel
! 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 \n"//&
                 "the interfaces.  Setting this to false reduces the \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "variables.", default=.true.)
  call get_param(param_file, mdl, "DT", cs%dt, &
                 "The (baroclinic) dynamics time step.", units="s", &
                 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 \n"//&
                 "large value if the velocity is truncated more than \n"//&
                 "MAXTRUNC times between energy saves.  Set MAXTRUNC to 0 \n"//&
                 "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 \n"//&
                 "model will be stopped if there is more energy than \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "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)
  end if

  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 \n"//&
                   "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 \n"//&
                   "entries in the depth-list file.", units="m", &
                   default=1.0e-10)
    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)
    endif

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

  cs%Huge_time = set_time(int(1e9),0)
  cs%Start_time = input_start_time
  cs%ntrunc => ntrnc

end subroutine mom_sum_output_init

subroutine mom_sum_output_end(CS)
  type(sum_output_cs), pointer :: CS
!   This subroutine deallocates the memory owned by this module.
! Argument: CS - The control structure returned by a previous call to
!                MOM_sum_output_init.

  if (associated(cs)) then
    if (cs%do_APE_calc) then
      dealloc_(cs%lH)
      deallocate(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, CS, tracer_CSp)
  type(ocean_grid_type),                     intent(in)    :: G   !< The ocean's grid structure.
  type(verticalgrid_type),                   intent(in)    :: GV  !< The ocean's vertical grid
                                                                  !! structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u   !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v   !< The meridional velocity,
                                                                  !! in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h   !< Layer thicknesses, in H
                                                                  !! (usually m or kg m-2).
  type(thermo_var_ptrs),                     intent(in)    :: tv  !< A structure pointing to various
                                                                  !! thermodynamic variables.
  type(time_type),                           intent(inout) :: 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


!  This subroutine calculates and writes the total model energy, the
! energy and mass of each layer, and other globally integrated
! physical quantities.

! Arguments: u - Zonal velocity, in m s-1.
!  (in)      v - Meridional velocity, in m s-1.
!  (in)      h - Layer thickness, in m.
!  (in)      tv - A structure containing pointers to any available
!                 thermodynamic fields, including potential temperature and
!                 salinity or mixed layer density. Absent fields have NULL ptrs.
!  (in/out)  day - The current model time.
!  (in)      n - The time step number of the current execution.
!  (in)      G - The ocean's grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 MOM_sum_output_init.

  real :: eta(szi_(g),szj_(g),szk_(g)+1) ! The height of interfaces, in m.
  real :: areaTm(szi_(g),szj_(g)) ! A masked version of areaT, in m2.
  real :: KE(szk_(g))  ! The total kinetic energy of a layer, in J.
  real :: PE(szk_(g)+1)! The available potential energy of an interface, in J.
  real :: KE_tot       ! The total kinetic energy, in J.
  real :: PE_tot       ! The total available potential energy, in J.
  real :: H_0APE(szk_(g)+1) ! The uniform depth which overlies the same
                       ! volume as is below an interface, in m.
                       ! H is usually positive.
  real :: toten        ! The total kinetic & potential energies of
                       ! all layers, in Joules (i.e. kg m2 s-2).
  real :: En_mass      ! The total kinetic and potential energies divided by
                       ! the total mass of the ocean, in m2 s-2.
  real :: vol_lay(szk_(g))  ! The volume of fluid in a layer, in m3.
  real :: volbelow     ! The volume of all layers beneath an interface in m3.
  real :: mass_lay(szk_(g)) ! The mass of fluid in a layer, in kg.
  real :: mass_tot     ! The total mass of the ocean in kg.
  real :: vol_tot      ! The total ocean volume in m3.
  real :: mass_chg     ! The change in total ocean mass of fresh water since
                       ! the last call to this subroutine, in kg.
  real :: mass_anom    ! The change in fresh water that cannot be accounted for
                       ! by the surface fluxes, in kg.
  real :: Salt         ! The total amount of salt in the ocean, in PSU kg.
  real :: Salt_chg     ! The change in total ocean salt since the last call
                       ! to this subroutine, in PSU kg.
  real :: Salt_anom    ! The change in salt that cannot be accounted for by
                       ! the surface fluxes, in PSU kg.
  real :: salin        ! The mean salinity of the ocean, in PSU.
  real :: salin_chg    ! The change in total salt since the last call
                       ! to this subroutine divided by total mass, in PSU.
  real :: salin_anom   ! The change in total salt that cannot be accounted for by
                       ! the surface fluxes divided by total mass in PSU.
  real :: salin_mass_in ! The mass of salt input since the last call, kg.
  real :: Heat         ! The total amount of Heat in the ocean, in Joules.
  real :: Heat_chg     ! The change in total ocean heat since the last call
                       ! to this subroutine, in Joules.
  real :: Heat_anom    ! The change in heat that cannot be accounted for by
                       ! the surface fluxes, in Joules.
  real :: temp         ! The mean potential temperature of the ocean, in C.
  real :: temp_chg     ! The change in total heat divided by total heat capacity
                       ! of the ocean since the last call to this subroutine, C.
  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, in C.
  real :: hint         ! The deviation of an interface from H, in m.
  real :: hbot         ! 0 if the basin is deeper than H, or the
                       ! height of the basin depth over H otherwise,
                       ! in m. This makes PE only include real fluid.
  real :: hbelow       ! The depth of fluid in all layers beneath
                       ! an interface, in m.
  type(efp_type) :: &
    mass_EFP, &        ! Extended fixed point sums of total mass, etc.
    salt_EFP, heat_EFP, salt_chg_EFP, heat_chg_EFP, mass_chg_EFP, &
    mass_anom_EFP, salt_anom_EFP, heat_anom_EFP
  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
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    tmp1
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)+1) :: &
    PE_pt
  real, dimension(SZI_(G),SZJ_(G)) :: &
    Temp_int, Salt_int
  real :: H_to_m, H_to_kg_m2  ! Local copies of unit conversion factors.
  integer :: num_nc_fields  ! The number of fields that will actually go into
                            ! the NetCDF file.
  integer :: i, j, k, is, ie, js, je, nz, m, Isq, Ieq, Jsq, Jeq
  integer :: l, lbelow, labove   ! indices of deep_area_vol, used to find
                                 ! H.  lbelow & labove are lower & upper
                                 ! limits for l in the search for lH.
  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
  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
  real, allocatable :: toten_PE(:)
  integer :: pe_num
  integer :: iyear, imonth, iday, ihour, iminute, isecond, itick ! For call to get_date()

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

  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

  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
  h_to_m = gv%H_to_m ; h_to_kg_m2 = gv%H_to_kg_m2

  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) * (h_to_kg_m2*areatm(i,j))
    enddo ; enddo ; enddo
    mass_tot = reproducing_sum(tmp1, sums=mass_lay, efp_sum=mass_efp)
    do k=1,nz ; vol_lay(k) = (h_to_m/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) = h_to_kg_m2 * h(i,j,k) * areatm(i,j)
      enddo ; enddo ; enddo
      mass_tot = reproducing_sum(tmp1, sums=mass_lay, efp_sum=mass_efp)

      call find_eta(h, tv, gv%g_Earth, g, gv, eta)
      do k=1,nz ; do j=js,je ; do i=is,ie
        tmp1(i,j,k) = (eta(i,j,k)-eta(i,j,k+1)) * areatm(i,j)
      enddo ; enddo ; enddo
      vol_tot = h_to_m*reproducing_sum(tmp1, sums=vol_lay)
    else
      do k=1,nz ; do j=js,je ; do i=is,ie
        tmp1(i,j,k) = h_to_kg_m2 * h(i,j,k) * areatm(i,j)
      enddo ; enddo ; enddo
      mass_tot = reproducing_sum(tmp1, sums=mass_lay, efp_sum=mass_efp)
      do k=1,nz ; vol_lay(k) = 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 = mass_tot ; cs%fresh_water_input = 0.0

    cs%mass_prev_EFP = mass_efp
    cs%fresh_water_in_EFP = real_to_efp(0.0)

    !  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]     "
          else if ((cs%timeunit >= 3599.0) .and. (cs%timeunit < 3601.0)) then
            time_units = "            [hours]      "
          else if ((cs%timeunit >= 86399.0) .and. (cs%timeunit < 86401.0)) then
            time_units = "             [days]      "
          else if ((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
      h_0ape(k) = cs%DL(l)%depth - (volbelow - cs%DL(l)%vol_below) / cs%DL(l)%area
    enddo
    h_0ape(nz+1) = cs%DL(2)%depth
  else
    do k=1,nz+1 ; h_0ape(k) = 0.0 ; enddo
  endif

! Calculate the Kinetic Energy integrated over each layer.
  tmp1(:,:,:) = 0.0
  do k=1,nz ; do j=js,je ; do i=is,ie
    tmp1(i,j,k) = (0.25 * h_to_kg_m2 * (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

!   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.
  do k=1,nz+1 ; pe(k) = 0.0 ; enddo
  if (cs%do_APE_calc) then
    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) * h_to_m
          hint = (h_0ape(k) + hbelow - g%bathyT(i,j))
          hbot = h_0ape(k) - g%bathyT(i,j)
          hbot = (hbot + abs(hbot)) * 0.5
          pe_pt(i,j,k) = 0.5 * 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
        hbelow = 0.0
        do k=nz,1,-1
          hint = h_0ape(k) + eta(i,j,k)  ! eta and H_0 have opposite signs.
          hbot = max(h_0ape(k) - g%bathyT(i,j), 0.0)
          pe_pt(i,j,k) = 0.5 * (areatm(i,j) * (gv%Rho0*gv%g_prime(k))) * &
                  (hint * hint - hbot * hbot)
        enddo
      enddo ; enddo
    endif
  endif

  ke_tot = reproducing_sum(tmp1, sums=ke)
  pe_tot = 0.0
  if (cs%do_APE_calc) &
    pe_tot = reproducing_sum(pe_pt, sums=pe)
  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)*(h_to_kg_m2 * areatm(i,j)))
      temp_int(i,j) = temp_int(i,j) + (tv%C_p * tv%T(i,j,k)) * &
                      (h(i,j,k)*(h_to_kg_m2 * areatm(i,j)))
    enddo ; enddo ; enddo
    salt = reproducing_sum(salt_int, efp_sum=salt_efp)
    heat = reproducing_sum(temp_int, efp_sum=heat_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
    if (u(i,j,k) < 0.0) then
      cfl_trans = (-u(i,j,k) * cs%dt) * (g%dy_Cu(i,j) * g%IareaT(i+1,j))
    else
      cfl_trans = (u(i,j,k) * cs%dt) * (g%dy_Cu(i,j) * g%IareaT(i,j))
    endif
    cfl_lin = abs(u(i,j,k) * cs%dt) * 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
    if (v(i,j,k) < 0.0) then
      cfl_trans = (-v(i,j,k) * cs%dt) * (g%dx_Cv(i,j) * g%IareaT(i,j+1))
    else
      cfl_trans = (v(i,j,k) * cs%dt) * (g%dx_Cv(i,j) * g%IareaT(i,j))
    endif
    cfl_lin = abs(v(i,j,k) * cs%dt) * 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(1))
  call max_across_pes(max_cfl(2))
  if (cs%use_temperature .and. cs%previous_calls == 0) then
    cs%salt_prev = salt ; cs%net_salt_input = 0.0
    cs%heat_prev = heat ; cs%net_heat_input = 0.0

    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
  irho0 = 1.0/gv%Rho0

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

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

  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*tv%C_p) ; temp_anom = heat_anom / (mass_tot*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 ",ES18.12, &
                               &", 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 ",ES18.12, &
                                &", 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(1), 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
    day = cs%Huge_time
    write(mesg,'("Energy per unit mass of ",ES11.4," exceeds ",ES11.4)') &
                  en_mass, cs%max_Energy
    call mom_error(warning, "write_energy : "//trim(mesg))
    call mom_error(warning, &
      "write_energy : Time set to a large value to force model termination.")
  endif
  if (cs%ntrunc>cs%maxtrunc) then
    day = cs%Huge_time
    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 = mass_tot ; cs%fresh_water_input = 0.0
  if (cs%use_temperature) then
    cs%salt_prev = salt ; cs%net_salt_input = 0.0
    cs%heat_prev = heat ; cs%net_heat_input = 0.0
  endif

  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, and perhaps later salt and
!! heat, through the ocean surface for use in diagnosing conservation.
subroutine accumulate_net_input(fluxes, state, dt, G, CS)
  type(forcing),         intent(in) :: fluxes !< A structure containing pointers to any possible forcing fields.  Unused fields are unallocated.
  type(surface),         intent(in) :: state
  real,                  intent(in) :: dt     !< The amount of time over which to average, in s.
  type(ocean_grid_type), intent(in) :: G      !< The ocean's grid structure.
  type(sum_output_cs),   pointer    :: CS     !< The control structure returned by a previous call to MOM_sum_output_init.

! This subroutine accumates the net input of volume, and perhaps later salt and
! heat, through the ocean surface for use in diagnosing conservation.
! Arguments: fluxes - A structure containing pointers to any possible
!                     forcing fields.  Unused fields are unallocated.
!  (in)      dt - The amount of time over which to average.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 MOM_sum_output_init.
  real, dimension(SZI_(G),SZJ_(G)) :: &
    FW_in, &   ! The net fresh water input, integrated over a timestep in kg.
    salt_in, & ! The total salt added by surface fluxes, integrated
               ! over a time step in ppt*kg.
    heat_in    ! The total heat added by surface fluxes, integrated
               ! over a time step in Joules.
  real :: FW_input   ! The net fresh water input, integrated over a timestep
                     ! and summed over space, in kg.
  real :: salt_input ! The total salt added by surface fluxes, integrated
                     ! over a time step and summed over space, in ppt * kg.
  real :: heat_input ! The total heat added by boundary fluxes, integrated
                     ! over a time step and summed over space, in Joules.
  real :: C_p        ! The heat capacity of seawater, in J K-1 kg-1.

  type(efp_type) :: &
    FW_in_EFP,      &  ! Extended fixed point versions of FW_input, salt_input, and
    salt_in_EFP,    &  ! heat_input, in kg, ppt*kg, and Joules.
    heat_in_EFP

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

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

  fw_in(:,:) = 0.0 ; fw_input = 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) = 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

  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*g%areaT(i,j) * (fluxes%sw(i,j) + &
             (fluxes%lw(i,j) + (fluxes%latent(i,j) + fluxes%sens(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*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(state%TempxPmE)) then
      do j=js,je ; do i=is,ie
        heat_in(i,j) = heat_in(i,j) + (c_p * g%areaT(i,j)) * state%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) + (c_p * state%SST(i,j)) * fw_in(i,j)
      enddo ; enddo
    endif


    ! The following heat sources may or may not be used.
    if (ASSOCIATED(state%internal_heat)) then
      do j=js,je ; do i=is,ie
        heat_in(i,j) = heat_in(i,j) + (c_p * g%areaT(i,j)) * &
                     state%internal_heat(i,j)
      enddo ; enddo
    endif
    if (ASSOCIATED(state%frazil)) then ; do j=js,je ; do i=is,ie
      heat_in(i,j) = heat_in(i,j) + g%areaT(i,j) * state%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) + dt*g%areaT(i,j)*fluxes%heat_added(i,j)
    enddo ; enddo ; endif
!    if (ASSOCIATED(state%sw_lost)) then ; do j=js,je ; do i=is,ie
!      heat_in(i,j) = heat_in(i,j) - G%areaT(i,j) * 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).
      salt_in(i,j) = 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
    fw_input   = reproducing_sum(fw_in,   efp_sum=fw_in_efp)
    heat_input = reproducing_sum(heat_in, efp_sum=heat_in_efp)
    salt_input = reproducing_sum(salt_in, efp_sum=salt_in_efp)

    cs%fresh_water_input = cs%fresh_water_input + fw_input
    cs%net_salt_input    = cs%net_salt_input    + salt_input
    cs%net_heat_input    = cs%net_heat_input    + heat_input

    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, CS)
  type(ocean_grid_type), intent(in) :: G    !< The ocean's grid structure
  type(sum_output_cs),   pointer    :: CS
!  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.

  integer :: k

  if (cs%read_depth_list) then
    if (file_exists(cs%depth_list_file)) then
      call read_depth_list(g, 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, 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.

  real, dimension(G%Domain%niglobal*G%Domain%njglobal + 1) :: &
    Dlist, &  !< The global list of bottom depths, in m.
    AreaList  !< The global list of cell areas, in 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, in m.
  real    :: Dprev !< The most recent depth that was considered, in m.
  real    :: vol   !< The running sum of open volume below a deptn, in m3.
  real    :: area  !< The open area at the current depth, in m2.
  real    :: D_list_prev !< The most recent depth added to the list, in 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, CS, filename, list_size)
  type(ocean_grid_type), intent(in) :: G    !< The ocean's grid structure.
  type(sum_output_cs),   pointer    :: CS
  character(len=*),      intent(in) :: filename
  integer,               intent(in) :: list_size

! This subroutine writes out the depth list to the specified file.

  real, allocatable :: tmp(:)
  integer :: ncid, dimid(1), Did, Aid, Vid, status, k

  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)))

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

  do k=1,list_size ; tmp(k) = 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) = 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) = 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, CS, filename)
  type(ocean_grid_type), intent(in) :: G    !< The ocean's grid structure
  type(sum_output_cs),   pointer    :: CS
  character(len=*),      intent(in) :: filename

! This subroutine reads in the depth list to the specified file
! and allocates and sets up CS%DL and CS%list_size .
  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)

  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

  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 = 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 = 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 = 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

end module mom_sum_output