mom_sum_output Module Reference

Data Types
type	depth_list
type	sum_output_cs

Functions/Subroutines
subroutine, public	mom_sum_output_init (G, param_file, directory, ntrnc, Input_start_time, CS)
subroutine	mom_sum_output_end (CS)
subroutine, public	write_energy (u, v, h, tv, day, n, G, GV, CS, tracer_CSp)
	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, state, dt, G, CS)
	This subroutine accumates the net input of volume, and perhaps later salt and heat, through the ocean surface for use in diagnosing conservation. More...
subroutine	depth_list_setup (G, 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, CS, filename, list_size)
	This subroutine writes out the depth list to the specified file. More...
subroutine	read_depth_list (G, CS, filename)
	This subroutine reads in the depth list to the specified file and allocates and sets up CSDL and CSlist_size . More...

Variables
integer, parameter	num_fields = 17

Function/Subroutine Documentation

accumulate_net_input()

subroutine, public mom_sum_output::accumulate_net_input	(	type(forcing), intent(in)	fluxes,
		type(surface), intent(in)	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, and perhaps later 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]	dt	The amount of time over which to average, in 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 890 of file MOM_sum_output.F90.

References mom_error_handler::mom_error().

Referenced by mom_main().

  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

Here is the call graph for this function:

Here is the caller graph for this function:

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 1050 of file MOM_sum_output.F90.

Referenced by depth_list_setup().

  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

Here is the caller graph for this function:

depth_list_setup()

subroutine mom_sum_output::depth_list_setup ( type(ocean_grid_type), intent(in)  G, 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

Definition at line 1018 of file MOM_sum_output.F90.

References create_depth_list(), mom_error_handler::is_root_pe(), mom_error_handler::mom_error(), read_depth_list(), and write_depth_list().

Referenced by mom_sum_output_init().

  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

Here is the call graph for this function:

Here is the caller graph for this function:

mom_sum_output_end()

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

private

Definition at line 286 of file MOM_sum_output.F90.

  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

mom_sum_output_init()

subroutine, public mom_sum_output::mom_sum_output_init	(	type(ocean_grid_type), intent(in)	G,
		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
	)

Parameters

[in]	g	The ocean's grid structure.
[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 164 of file MOM_sum_output.F90.

References depth_list_setup(), and mom_error_handler::mom_error().

Referenced by mom_main(), and ocean_model_mod::ocean_model_init().

  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

Here is the call graph for this function:

Here is the caller graph for this function:

read_depth_list()

subroutine mom_sum_output::read_depth_list ( type(ocean_grid_type), intent(in)  G, 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

Definition at line 1257 of file MOM_sum_output.F90.

References mom_error_handler::mom_error().

Referenced by depth_list_setup().

  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)

Here is the call graph for this function:

Here is the caller graph for this function:

write_depth_list()

subroutine mom_sum_output::write_depth_list ( type(ocean_grid_type), intent(in)  G, 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.

Definition at line 1175 of file MOM_sum_output.F90.

References mom_error_handler::is_root_pe(), and mom_error_handler::mom_error().

Referenced by depth_list_setup().

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

Here is the call graph for this function:

Here is the caller graph for this function:

write_energy()

subroutine, public mom_sum_output::write_energy	(	real, dimension( g %isdb: g %iedb, g %jsd: g %jed, g %ke), intent(in)	u,
		real, dimension( g %isd: g %ied, g %jsdb: g %jedb, g %ke), intent(in)	v,
		real, dimension( g %isd: g %ied, g %jsd: g %jed, g %ke), intent(in)	h,
		type(thermo_var_ptrs), intent(in)	tv,
		type(time_type), intent(inout)	day,
		integer, intent(in)	n,
		type(ocean_grid_type), intent(in)	G,
		type(verticalgrid_type), intent(in)	GV,
		type(sum_output_cs), pointer	CS,
		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.

Parameters

[in]	g	The ocean's grid structure.
[in]	gv	The ocean's vertical grid structure.
[in]	u	The zonal velocity, in m s-1.
[in]	v	The meridional velocity, in m s-1.
[in]	h	Layer thicknesses, in H (usually m or kg m-2).
[in]	tv	A structure pointing to various thermodynamic variables.
[in,out]	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.

Definition at line 305 of file MOM_sum_output.F90.

References mom_tracer_flow_control::call_tracer_stocks(), mom_coms::efp_to_real(), mom_error_handler::is_root_pe(), mom_error_handler::mom_error(), and mom_coms::real_to_efp().

Referenced by mom_main().

  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

Here is the call graph for this function:

Here is the caller graph for this function:

Variable Documentation

num_fields

integer, parameter mom_sum_output::num_fields = 17

private

Definition at line 85 of file MOM_sum_output.F90.

integer, parameter :: num_fields = 17