mom_sum_output Module Reference

Detailed Description

Reports integrated quantities for monitoring the model state.

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.

Data Types
type	depth_list
	A list of depths and corresponding globally integrated ocean area at each depth and the ocean volume below each depth. More...
type	sum_output_cs
	The control structure for the MOM_sum_output module. More...

Functions/Subroutines
subroutine, public	mom_sum_output_init (G, US, param_file, directory, ntrnc, Input_start_time, CS)
	MOM_sum_output_init initializes the parameters and settings for the MOM_sum_output module. More...
subroutine	mom_sum_output_end (CS)
	MOM_sum_output_end deallocates memory used by the MOM_sum_output module. More...
subroutine, public	write_energy (u, v, h, tv, day, n, G, GV, US, CS, tracer_CSp, OBC, dt_forcing)
	This subroutine calculates and writes the total model energy, the energy and mass of each layer, and other globally integrated physical quantities. More...
subroutine, public	accumulate_net_input (fluxes, sfc_state, dt, G, CS)
	This subroutine accumates the net input of volume, salt and heat, through the ocean surface for use in diagnosing conservation. More...
subroutine	depth_list_setup (G, US, 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. More...
subroutine	create_depth_list (G, CS)
	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. More...
subroutine	write_depth_list (G, US, CS, filename, list_size)
	This subroutine writes out the depth list to the specified file. More...
subroutine	read_depth_list (G, US, CS, filename)
	This subroutine reads in the depth list to the specified file and allocates and sets up CSDL and CSlist_size . More...
subroutine	get_depth_list_checksums (G, depth_chksum, area_chksum)
	Return the checksums required to verify DEPTH_LIST_FILE contents. More...

Variables
integer, parameter	num_fields = 17
	Number of diagnostic fields.
character(*), parameter	depth_chksum_attr = "bathyT_checksum"
	Checksum attribute name of GbathyT over the compute domain.
character(*), parameter	area_chksum_attr = "mask2dT_areaT_checksum"
	Checksum attribute of name of Gmask2dT * GareaT over the compute domain.

Function/Subroutine Documentation

accumulate_net_input()

subroutine, public mom_sum_output::accumulate_net_input	(	type(forcing), intent(in)	fluxes,
		type(surface), intent(in)	sfc_state,
		real, intent(in)	dt,
		type(ocean_grid_type), intent(in)	G,
		type(sum_output_cs), pointer	CS
	)

This subroutine accumates the net input of volume, salt and heat, through the ocean surface for use in diagnosing conservation.

Parameters

[in]	fluxes	A structure containing pointers to any possible forcing fields. Unused fields are unallocated.
[in]	sfc_state	A structure containing fields that describe the surface state of the ocean.
[in]	dt	The amount of time over which to average [s].
[in]	g	The ocean's grid structure.
	cs	The control structure returned by a previous call to MOM_sum_output_init.

Definition at line 937 of file MOM_sum_output.F90.

  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.
  real,                  intent(in) :: dt     !< The amount of time over which to average [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.
  ! 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 :: C_p        ! The heat capacity of seawater [J degC-1 kg-1].
 
  type(EFP_type) :: &
    FW_in_EFP,   & ! Extended fixed point version of FW_input [kg]
    salt_in_EFP, & ! Extended fixed point version of salt_input [ppt kg]
    heat_in_EFP    ! Extended fixed point version of heat_input [J]
 
  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
 
  if (associated(fluxes%seaice_melt)) then ; do j=js,je ; do i=is,ie
    fw_in(i,j) = fw_in(i,j) + 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*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*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*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(sfc_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)) * sfc_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 * sfc_state%SST(i,j)) * fw_in(i,j)
      enddo ; enddo
    endif
 
 
    ! The following heat sources may or may not be used.
    if (associated(sfc_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)) * &
                     sfc_state%internal_heat(i,j)
      enddo ; enddo
    endif
    if (associated(sfc_state%frazil)) then ; do j=js,je ; do i=is,ie
      heat_in(i,j) = heat_in(i,j) + g%areaT(i,j) * sfc_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(sfc_state%sw_lost)) then ; do j=js,je ; do i=is,ie
!      heat_in(i,j) = heat_in(i,j) - 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) = 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
 

create_depth_list()

subroutine mom_sum_output::create_depth_list ( type(ocean_grid_type), intent(in)  G, type(sum_output_cs), pointer  CS  )

private

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.

Parameters

[in]	g	The ocean's grid structure.
	cs	The control structure set up in MOM_sum_output_init, in which the ordered depth list is stored.

Definition at line 1098 of file MOM_sum_output.F90.

  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 [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 m2 ~> m3].
  real    :: area  !< The open area at the current depth [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
 

depth_list_setup()

subroutine mom_sum_output::depth_list_setup ( type(ocean_grid_type), intent(in)  G, type(unit_scale_type), intent(in)  US, type(sum_output_cs), pointer  CS  )

private

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.

Parameters

[in]	g	The ocean's grid structure
[in]	us	A dimensional unit scaling type
	cs	The control structure returned by a previous call to MOM_sum_output_init.

Definition at line 1068 of file MOM_sum_output.F90.

  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
 

get_depth_list_checksums()

subroutine mom_sum_output::get_depth_list_checksums ( type(ocean_grid_type), intent(in)  G, character(len=16), intent(out)  depth_chksum, character(len=16), intent(out)  area_chksum  )

private

Return the checksums required to verify DEPTH_LIST_FILE contents.

This function computes checksums for the bathymetry (GbathyT) 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.

Parameters

[in]	g	Ocean grid structure
[out]	depth_chksum	Depth checksum hexstring
[out]	area_chksum	Area checksum hexstring

Definition at line 1474 of file MOM_sum_output.F90.

  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%areaT(i,j)
  enddo ; enddo
  write(area_chksum, '(Z16)') mpp_chksum(field(:,:))
 
  deallocate(field)

mom_sum_output_end()

subroutine mom_sum_output::mom_sum_output_end ( type(sum_output_cs), pointer  CS )

private

MOM_sum_output_end deallocates memory used by the MOM_sum_output module.

Parameters

cs	The control structure returned by a previous call to MOM_sum_output_init.

Definition at line 290 of file MOM_sum_output.F90.

  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

mom_sum_output_init()

subroutine, public mom_sum_output::mom_sum_output_init	(	type(ocean_grid_type), intent(in)	G,
		type(unit_scale_type), intent(in)	US,
		type(param_file_type), intent(in)	param_file,
		character(len=*), intent(in)	directory,
		integer, intent(inout), target	ntrnc,
		type(time_type), intent(in)	Input_start_time,
		type(sum_output_cs), pointer	CS
	)

MOM_sum_output_init initializes the parameters and settings for the MOM_sum_output module.

Parameters

[in]	g	The ocean's grid structure.
[in]	us	A dimensional unit scaling type
[in]	param_file	A structure to parse for run-time parameters.
[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.
	cs	A pointer that is set to point to the control structure for this module.

Definition at line 142 of file MOM_sum_output.F90.

  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, &
                 "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 "//&
                 "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
 

read_depth_list()

subroutine mom_sum_output::read_depth_list ( type(ocean_grid_type), intent(in)  G, type(unit_scale_type), intent(in)  US, type(sum_output_cs), pointer  CS, character(len=*), intent(in)  filename  )

private

This subroutine reads in the depth list to the specified file and allocates and sets up CSDL and CSlist_size .

Parameters

[in]	g	The ocean's grid structure
[in]	us	A dimensional unit scaling type
	cs	The control structure returned by a previous call to MOM_sum_output_init.
[in]	filename	The path to the depth list file to read.

Definition at line 1318 of file MOM_sum_output.F90.

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

write_depth_list()

subroutine mom_sum_output::write_depth_list ( type(ocean_grid_type), intent(in)  G, type(unit_scale_type), intent(in)  US, type(sum_output_cs), pointer  CS, character(len=*), intent(in)  filename, integer, intent(in)  list_size  )

private

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

Parameters

[in]	g	The ocean's grid structure.
[in]	us	A dimensional unit scaling type
	cs	The control structure returned by a previous call to MOM_sum_output_init.
[in]	filename	The path to the depth list file to write.
[in]	list_size	The size of the depth list.

Definition at line 1223 of file MOM_sum_output.F90.

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

write_energy()

subroutine, public mom_sum_output::write_energy	(	real, dimension(szib_(g),szj_(g),szk_(g)), intent(in)	u,
		real, dimension(szi_(g),szjb_(g),szk_(g)), intent(in)	v,
		real, dimension(szi_(g),szj_(g),szk_(g)), intent(in)	h,
		type(thermo_var_ptrs), intent(in)	tv,
		type(time_type), intent(in)	day,
		integer, intent(in)	n,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(unit_scale_type), intent(in)	US,
		type(sum_output_cs), pointer	CS,
		type(tracer_flow_control_cs), optional, pointer	tracer_CSp,
		type(ocean_obc_type), optional, pointer	OBC,
		type(time_type), intent(in), optional	dt_forcing
	)

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

Parameters

[in]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	us	A dimensional unit scaling type
[in]	u	The zonal velocity [m s-1].
[in]	v	The meridional velocity [m s-1].
[in]	h	Layer thicknesses [H ~> m or kg m-2].
[in]	tv	A structure pointing to various thermodynamic variables.
[in]	day	The current model time.
[in]	n	The time step number of the current execution.
	cs	The control structure returned by a previous call to MOM_sum_output_init.
	tracer_csp	tracer control structure.
	obc	Open boundaries control structure.
[in]	dt_forcing	The forcing time step

Definition at line 304 of file MOM_sum_output.F90.

  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 [m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: v   !< The meridional velocity [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 [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 m2 ~> m3].
  real :: volbelow     ! The volume of all layers beneath an interface [Z m2 ~> 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 :: salin_mass_in ! The mass of salt input since the last call [kg].
  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, &        ! 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        ! 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 :: H_to_kg_m2   ! Local copy of a unit conversion factor.
  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
  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
  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
 
    ! 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, sums=mass_lay, efp_sum=mass_efp)
    do k=1,nz ; vol_lay(k) = (gv%H_to_Z/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, g, gv, us, 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 = us%Z_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) = us%m_to_Z * (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]     "
          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.
    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 * areatm(i,j) * us%Z_to_m*(gv%Rho0*us%L_to_m**2*us%s_to_T**2*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 * (areatm(i,j) * us%Z_to_m*(gv%Rho0*us%L_to_m**2*us%s_to_T**2*gv%g_prime(k))) * &
                  (hint * hint - hbot * hbot)
        enddo
      enddo ; enddo
    endif
 
    pe_tot = reproducing_sum(pe_pt, 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.
  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
  ke_tot = reproducing_sum(tmp1, 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)*(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 ppt 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 ",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(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
    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 = 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