MOM_diagnostics.F90

Go to the documentation of this file.

module mom_diagnostics

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

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, February 2001                                  *
!*                                                                     *
!*    This subroutine calculates any requested diagnostic quantities   *
!*  that are not calculated in the various subroutines.  Diagnostic    *
!*  quantities are requested by allocating them memory.                *
!*                                                                     *
!*  Macros written all in capital letters are defined in MOM_memory.h. *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q, CoriolisBu                            *
!*    j+1  > o > o >   At ^:  v                                        *
!*    j    x ^ x ^ x   At >:  u                                        *
!*    j    > o > o >   At o:  h, bathyT                                *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**

use mom_coms,              only : reproducing_sum
use mom_diag_mediator,     only : post_data, post_data_1d_k
use mom_diag_mediator,     only : register_diag_field, register_scalar_field
use mom_diag_mediator,     only : diag_ctrl, time_type, safe_alloc_ptr
use mom_domains,           only : create_group_pass, do_group_pass, group_pass_type
use mom_domains,           only : to_north, to_east
use mom_eos,               only : calculate_density, int_density_dz
use mom_error_handler,     only : mom_error, fatal, warning
use mom_file_parser,       only : get_param, log_version, param_file_type
use mom_forcing_type,      only : forcing
use mom_grid,              only : ocean_grid_type
use mom_interface_heights, only : find_eta
use mom_spatial_means,     only : global_area_mean, global_layer_mean, global_volume_mean
use mom_variables,         only : thermo_var_ptrs, ocean_internal_state, p3d
use mom_variables,         only : accel_diag_ptrs, cont_diag_ptrs
use mom_verticalgrid,      only : verticalgrid_type
use mom_wave_speed,        only : wave_speed, wave_speed_cs, wave_speed_init

implicit none ; private

#include <MOM_memory.h>

public calculate_diagnostic_fields
public register_time_deriv
public find_eta
public mom_diagnostics_init
public mom_diagnostics_end


type, public :: diagnostics_cs ; private
  real :: mono_n2_column_fraction = 0. !< The lower fraction of water column over which N2 is limited as
                                       !! monotonic for the purposes of calculating the equivalent barotropic wave speed.
  real :: mono_n2_depth = -1.          !< The depth below which N2 is limited as monotonic for the purposes of
                                       !! calculating the equivalent barotropic wave speed. (m)

  type(diag_ctrl), pointer :: diag ! structure to regulate diagnostics timing

  ! following arrays store diagnostics calculated here and unavailable outside.

  ! following fields have nz+1 levels.
  real, pointer, dimension(:,:,:) :: &
    e => null(), &   ! interface height (metre)
    e_d => null()    ! interface height above bottom (metre)

  ! following fields have nz layers.
  real, pointer, dimension(:,:,:) :: &
    du_dt => null(), &    ! net i-acceleration in m/s2
    dv_dt => null(), &    ! net j-acceleration in m/s2
    dh_dt => null(), &    ! thickness rate of change in (m/s) or kg/(m2*s)

    h_rlay => null(),    & ! layer thicknesses in layered potential density
                           ! coordinates, in m (Bouss) or kg/m2 (non-Bouss)
    uh_rlay   => null(), & ! zonal and meridional transports in layered
    vh_rlay   => null(), & ! potential rho coordinates: m3/s(Bouss) kg/s(non-Bouss)
    uhgm_rlay => null(), & ! zonal and meridional Gent-McWilliams transports in layered
    vhgm_rlay => null(), & ! potential density coordinates, m3/s (Bouss) kg/s(non-Bouss)
    p_ebt     => null()    ! Equivalent barotropic modal structure

  ! following fields are 2-D.
  real, pointer, dimension(:,:) :: &
    cg1 => null(),       & ! first baroclinic gravity wave speed, in m s-1
    rd1 => null(),       & ! first baroclinic deformation radius, in m
    cfl_cg1 => null(),   & ! CFL for first baroclinic gravity wave speed, nondim
    cfl_cg1_x => null(), & ! i-component of CFL for first baroclinic gravity wave speed, nondim
    cfl_cg1_y => null()    ! j-component of CFL for first baroclinic gravity wave speed, nondim

  ! arrays to hold diagnostics in the layer-integrated energy budget.
  ! all except KE have units of m3 s-3 (when Boussinesq).
  real, pointer, dimension(:,:,:) :: &
    ke        => null(), &  ! KE per unit mass, in m2 s-2
    dke_dt    => null(), &  ! time derivative of the layer KE
    pe_to_ke  => null(), &  ! potential energy to KE term
    ke_coradv => null(), &  ! KE source from the combined Coriolis and
                            ! advection terms.  The Coriolis source should be
                            ! zero, but is not due to truncation errors.  There
                            ! should be near-cancellation of the global integral
                            ! of this spurious Coriolis source.
    ke_adv     => null(),&  ! KE source from along-layer advection
    ke_visc    => null(),&  ! KE source from vertical viscosity
    ke_horvisc => null(),&  ! KE source from horizontal viscosity
    ke_dia     => null(),&  ! KE source from diapycnal diffusion
    diag_tmp3d => null()    ! 3D re-usable arrays for diagnostics

  ! diagnostic IDs
  integer :: id_e              = -1, id_e_d            = -1
  integer :: id_du_dt          = -1, id_dv_dt          = -1
  integer :: id_col_ht         = -1, id_dh_dt          = -1
  integer :: id_ke             = -1, id_dkedt          = -1
  integer :: id_pe_to_ke       = -1, id_ke_coradv      = -1
  integer :: id_ke_adv         = -1, id_ke_visc        = -1
  integer :: id_ke_horvisc     = -1, id_ke_dia         = -1
  integer :: id_uh_rlay        = -1, id_vh_rlay        = -1
  integer :: id_uhgm_rlay      = -1, id_vhgm_rlay      = -1
  integer :: id_h_rlay         = -1, id_rd1            = -1
  integer :: id_rml            = -1, id_rcv            = -1
  integer :: id_cg1            = -1, id_cfl_cg1        = -1
  integer :: id_cfl_cg1_x      = -1, id_cfl_cg1_y      = -1
  integer :: id_cg_ebt         = -1, id_rd_ebt         = -1
  integer :: id_p_ebt          = -1
  integer :: id_temp_int       = -1, id_salt_int       = -1
  integer :: id_mass_wt        = -1, id_col_mass       = -1
  integer :: id_masscello      = -1, id_masso          = -1
  integer :: id_thetaoga       = -1, id_soga           = -1
  integer :: id_sosga          = -1, id_tosga          = -1
  integer :: id_temp_layer_ave = -1, id_salt_layer_ave = -1
  integer :: id_pbo            = -1
  integer :: id_thkcello       = -1, id_rhoinsitu      = -1
  integer :: id_rhopot0        = -1, id_rhopot2        = -1

  type(wave_speed_cs), pointer :: wave_speed_csp => null()

  ! pointers used in calculation of time derivatives
  type(p3d) :: var_ptr(max_fields_)
  type(p3d) :: deriv(max_fields_)
  type(p3d) :: prev_val(max_fields_)
  integer   :: nlay(max_fields_)
  integer   :: num_time_deriv = 0

  ! for group halo pass
  type(group_pass_type) :: pass_ke_uv

end type diagnostics_cs

contains
!> Diagnostics not more naturally calculated elsewhere are computed here.
subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, fluxes, &
                                       dt, G, GV, CS, eta_bt)
  type(ocean_grid_type),                     intent(inout) :: 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).
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: uh   !< Transport through zonal faces
                                                                   !! = u*h*dy, m3/s(Bouss)
                                                                   !! kg/s(non-Bouss).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: vh   !< transport through meridional
                                                                   !! faces = v*h*dx, m3/s(Bouss)
                                                                   !! kg/s(non-Bouss).
  type(thermo_var_ptrs),                     intent(in)    :: tv   !< A structure pointing to
                                                                   !! various thermodynamic
                                                                   !! variables.
  type(accel_diag_ptrs),                     intent(in)    :: ADp  !< structure with pointers to
                                                                   !! accelerations in momentum
                                                                   !! equation.
  type(cont_diag_ptrs),                      intent(in)    :: CDp  !< structure with pointers to
                                                                   !! terms in continuity equation.
  type(forcing),                             intent(in)    :: fluxes !< A structure containing the
                                                                   !! surface fluxes.
  real,                                      intent(in)    :: dt   !< The time difference in s since
                                                                   !! the last call to this
                                                                   !! subroutine.
  type(diagnostics_cs),                      intent(inout) :: CS   !< Control structure returned by
                                                                   !! a previous call to
                                                                   !! diagnostics_init.
  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in)   :: eta_bt !< An optional barotropic
    !! variable that gives the "correct" free surface height (Boussinesq) or total water column
    !! mass per unit area (non-Boussinesq).  This is used to dilate the layer thicknesses when
    !! calculating interface heights, in m or kg m-2.

! Diagnostics not more naturally calculated elsewhere are computed here.

! Arguments:
!  (in)      u   - zonal velocity component (m/s).
!  (in)      v   - meridional velocity component (m/s).
!  (in)      h   - layer thickness,  meter(Bouss)  kg/m^2(non-Bouss).
!  (in)      uh  - transport through zonal faces = u*h*dy, m3/s(Bouss) kg/s(non-Bouss).
!  (in)      vh  - transport through meridional faces = v*h*dx, m3/s(Bouss) kg/s(non-Bouss).
!  (in)      tv  - structure pointing to various thermodynamic variables.
!  (in)      ADp - structure with pointers to accelerations in momentum equation.
!  (in)      CDp - structure with pointers to terms in continuity equation.
!  (in)      dt  - time difference in s since the last call to this subroutine.
!  (in)      G   - ocean grid structure.
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS  - control structure returned by a previous call to diagnostics_init.
!  (in,opt)  eta_bt - An optional barotropic variable that gives the "correct"
!                     free surface height (Boussinesq) or total water column
!                     mass per unit area (non-Boussinesq).  This is used to
!                     dilate the layer thicknesses when calculating interface
!                     heights, in m or kg m-2.

  integer i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz, nkmb

  ! coordinate variable potential density, in kg m-3.
  real :: Rcv(szi_(g),szj_(g),szk_(g))

  ! tmp array for surface properties
  real :: surface_field(szi_(g),szj_(g))
  real :: pressure_1d(szi_(g)) ! Temporary array for pressure when calling EOS
  real :: wt, wt_p

  ! squared Coriolis parameter at to h-points (1/s2)
  real :: f2_h

  ! magnitude of the gradient of f (1/(m*s))
  real :: mag_beta

  ! frequency squared used to avoid division by 0 (1/s2)
  ! value is roughly (pi / (the age of the universe) )^2.
  real, parameter :: absurdly_small_freq2 = 1e-34

  integer :: k_list

  real, dimension(SZK_(G)) :: temp_layer_ave, salt_layer_ave
  real :: thetaoga, soga, masso, tosga, sosga

  is  = g%isc  ; ie   = g%iec  ; js  = g%jsc  ; je  = g%jec
  isq = g%IscB ; ieq  = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  nz  = g%ke   ; nkmb = gv%nk_rho_varies

  ! smg: is the following robust to ALE? It seems a bit opaque.
  ! If the model is NOT in isopycnal mode then nkmb=0. But we need all the
  ! following diagnostics to treat all layers as variable density, so we set
  ! nkmb = nz, on the expectation that loops nkmb+1,nz will not iterate.
  ! This behavior is ANSI F77 but some compiler options can force at least
  ! one iteration that would break the following one-line workaround!
  if (nkmb==0) nkmb = nz

  if (loc(cs)==0) call mom_error(fatal, &
         "calculate_diagnostic_fields: Module must be initialized before used.")

  call calculate_derivs(dt, g, cs)

  if (ASSOCIATED(cs%e)) then
    call find_eta(h, tv, gv%g_Earth, g, gv, cs%e, eta_bt)
    if (cs%id_e > 0) call post_data(cs%id_e, cs%e, cs%diag)
  endif

  if (ASSOCIATED(cs%e_D)) then
    if (ASSOCIATED(cs%e)) then
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        cs%e_D(i,j,k) = cs%e(i,j,k) + g%bathyT(i,j)
      enddo ; enddo ; enddo
    else
      call find_eta(h, tv, gv%g_Earth, g, gv, cs%e_D, eta_bt)
      do k=1,nz+1 ; do j=js,je ; do i=is,ie
        cs%e_D(i,j,k) = cs%e_D(i,j,k) + g%bathyT(i,j)
      enddo ; enddo ; enddo
    endif

    if (cs%id_e_D > 0) call post_data(cs%id_e_D, cs%e_D, cs%diag)
  endif

  ! mass per area of grid cell (for Bouss, use Rho0)
  if (cs%id_masscello > 0) then
    do k=1,nz; do j=js,je ; do i=is,ie
       cs%diag_tmp3d(i,j,k) = gv%H_to_kg_m2*h(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_masscello, cs%diag_tmp3d, cs%diag)
  endif

  ! mass of liquid ocean (for Bouss, use Rho0)
  if (cs%id_masso > 0) then
    do k=1,nz; do j=js,je ; do i=is,ie
       cs%diag_tmp3d(i,j,k) = gv%H_to_kg_m2*h(i,j,k)*g%areaT(i,j)
    enddo ; enddo ; enddo
    masso = (reproducing_sum(sum(cs%diag_tmp3d,3)))
    call post_data(cs%id_masso, masso, cs%diag)
  endif

  ! diagnose thickness of grid cells (meter)
  if (cs%id_thkcello > 0) then

    ! thkcello = h for Boussinesq
    if (gv%Boussinesq) then
      call post_data(cs%id_thkcello, gv%H_to_m*h, cs%diag)

    ! thkcello = dp/(rho*g) for non-Boussinesq
    else
      do j=js,je

        ! pressure loading at top of surface layer (Pa)
        if(ASSOCIATED(fluxes%p_surf)) then
          do i=is,ie
            pressure_1d(i) = fluxes%p_surf(i,j)
          enddo
        else
          do i=is,ie
            pressure_1d(i) = 0.0
          enddo
        endif

        do k=1,nz
          ! pressure for EOS at the layer center (Pa)
          do i=is,ie
            pressure_1d(i) = pressure_1d(i) + 0.5*(gv%g_Earth*gv%H_to_kg_m2)*h(i,j,k)
          enddo
          ! store in-situ density (kg/m3) in diag_tmp3d
          call calculate_density(tv%T(:,j,k),tv%S(:,j,k), pressure_1d, &
                                 cs%diag_tmp3d(:,j,k), is, ie-is+1, tv%eqn_of_state)
          ! cell thickness = dz = dp/(g*rho) (meter); store in diag_tmp3d
          do i=is,ie
            cs%diag_tmp3d(i,j,k) = (gv%H_to_kg_m2*h(i,j,k))/cs%diag_tmp3d(i,j,k)
          enddo
          ! pressure for EOS at the bottom interface (Pa)
          do i=is,ie
            pressure_1d(i) = pressure_1d(i) + 0.5*(gv%g_Earth*gv%H_to_kg_m2)*h(i,j,k)
          enddo
        enddo ! k

      enddo ! j
      call post_data(cs%id_thkcello, cs%diag_tmp3d, cs%diag)
    endif
  endif

  ! volume mean potential temperature
  if (cs%id_thetaoga>0) then
    thetaoga = global_volume_mean(tv%T, h, g, gv)
    call post_data(cs%id_thetaoga, thetaoga, cs%diag)
  endif

  ! area mean SST
  if (cs%id_tosga > 0) then
    do j=js,je ; do i=is,ie
       surface_field(i,j) = tv%T(i,j,1)
    enddo ; enddo
    tosga = global_area_mean(surface_field, g)
    call post_data(cs%id_tosga, tosga, cs%diag)
  endif

  ! volume mean salinity
  if (cs%id_soga>0) then
    soga = global_volume_mean(tv%S, h, g, gv)
    call post_data(cs%id_soga, soga, cs%diag)
  endif

  ! area mean SSS
  if (cs%id_sosga > 0) then
    do j=js,je ; do i=is,ie
       surface_field(i,j) = tv%S(i,j,1)
    enddo ; enddo
    sosga = global_area_mean(surface_field, g)
    call post_data(cs%id_sosga, sosga, cs%diag)
  endif

  ! layer mean potential temperature
  if (cs%id_temp_layer_ave>0) then
    temp_layer_ave = global_layer_mean(tv%T, h, g, gv)
    call post_data_1d_k(cs%id_temp_layer_ave, temp_layer_ave, cs%diag)
  endif

  ! layer mean salinity
  if (cs%id_salt_layer_ave>0) then
    salt_layer_ave = global_layer_mean(tv%S, h, g, gv)
    call post_data_1d_k(cs%id_salt_layer_ave, salt_layer_ave, cs%diag)
  endif

  call calculate_vertical_integrals(h, tv, fluxes, g, gv, cs)

  if ((cs%id_Rml > 0) .or. (cs%id_Rcv > 0) .or. ASSOCIATED(cs%h_Rlay) .or. &
      ASSOCIATED(cs%uh_Rlay) .or. ASSOCIATED(cs%vh_Rlay) .or. &
      ASSOCIATED(cs%uhGM_Rlay) .or. ASSOCIATED(cs%vhGM_Rlay)) then

    if (associated(tv%eqn_of_state)) then
      pressure_1d(:) = tv%P_Ref
!$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d)
      do k=1,nz ; do j=js-1,je+1
        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, &
                               rcv(:,j,k), is-1, ie-is+3, tv%eqn_of_state)
      enddo ; enddo
    else ! Rcv should not be used much in this case, so fill in sensible values.
      do k=1,nz ; do j=js-1,je+1 ; do i=is-1,ie+1
        rcv(i,j,k) = gv%Rlay(k)
      enddo ; enddo ; enddo
    endif
    if (cs%id_Rml > 0) call post_data(cs%id_Rml, rcv, cs%diag)
    if (cs%id_Rcv > 0) call post_data(cs%id_Rcv, rcv, cs%diag)

    if (ASSOCIATED(cs%h_Rlay)) then
      k_list = nz/2
!$OMP parallel do default(none) shared(is,ie,js,je,nz,nkmb,CS,Rcv,h,GV) &
!$OMP                          private(wt,wt_p) firstprivate(k_list)
      do j=js,je
        do k=1,nkmb ; do i=is,ie
          cs%h_Rlay(i,j,k) = 0.0
        enddo ; enddo
        do k=nkmb+1,nz ; do i=is,ie
          cs%h_Rlay(i,j,k) = h(i,j,k)
        enddo ; enddo
        do k=1,nkmb ; do i=is,ie
          call find_weights(gv%Rlay, rcv(i,j,k), k_list, nz, wt, wt_p)
          cs%h_Rlay(i,j,k_list)   = cs%h_Rlay(i,j,k_list)   + h(i,j,k)*wt
          cs%h_Rlay(i,j,k_list+1) = cs%h_Rlay(i,j,k_list+1) + h(i,j,k)*wt_p
        enddo ; enddo
      enddo

      if (cs%id_h_Rlay > 0) call post_data(cs%id_h_Rlay, cs%h_Rlay, cs%diag)
    endif

    if (ASSOCIATED(cs%uh_Rlay)) then
      k_list = nz/2
!$OMP parallel do default(none) shared(Isq,Ieq,js,je,nz,nkmb,Rcv,CS,GV,uh) &
!$OMP                          private(wt,wt_p) firstprivate(k_list)
      do j=js,je
        do k=1,nkmb ; do i=isq,ieq
          cs%uh_Rlay(i,j,k) = 0.0
        enddo ; enddo
        do k=nkmb+1,nz ; do i=isq,ieq
          cs%uh_Rlay(i,j,k) = uh(i,j,k)
        enddo ; enddo
        k_list = nz/2
        do k=1,nkmb ; do i=isq,ieq
          call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i+1,j,k)), k_list, nz, wt, wt_p)
          cs%uh_Rlay(i,j,k_list)   = cs%uh_Rlay(i,j,k_list)   + uh(i,j,k)*wt
          cs%uh_Rlay(i,j,k_list+1) = cs%uh_Rlay(i,j,k_list+1) + uh(i,j,k)*wt_p
        enddo ; enddo
      enddo

      if (cs%id_uh_Rlay > 0) call post_data(cs%id_uh_Rlay, cs%uh_Rlay, cs%diag)
    endif

    if (ASSOCIATED(cs%vh_Rlay)) then
      k_list = nz/2
!$OMP parallel do default(none)  shared(Jsq,Jeq,is,ie,nz,nkmb,Rcv,CS,GV,vh) &
!$OMP                          private(wt,wt_p) firstprivate(k_list)
      do j=jsq,jeq
        do k=1,nkmb ; do i=is,ie
          cs%vh_Rlay(i,j,k) = 0.0
        enddo ; enddo
        do k=nkmb+1,nz ; do i=is,ie
          cs%vh_Rlay(i,j,k) = vh(i,j,k)
        enddo ; enddo
        do k=1,nkmb ; do i=is,ie
          call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i,j+1,k)), k_list, nz, wt, wt_p)
          cs%vh_Rlay(i,j,k_list)   = cs%vh_Rlay(i,j,k_list)   + vh(i,j,k)*wt
          cs%vh_Rlay(i,j,k_list+1) = cs%vh_Rlay(i,j,k_list+1) + vh(i,j,k)*wt_p
        enddo ; enddo
      enddo

      if (cs%id_vh_Rlay > 0) call post_data(cs%id_vh_Rlay, cs%vh_Rlay, cs%diag)
    endif

    if (ASSOCIATED(cs%uhGM_Rlay) .and. ASSOCIATED(cdp%uhGM)) then
      k_list = nz/2
!$OMP parallel do default(none) shared(Isq,Ieq,js,je,nz,nkmb,Rcv,CDP,CS,GV) &
!$OMP                          private(wt,wt_p) firstprivate(k_list)
      do j=js,je
        do k=1,nkmb ; do i=isq,ieq
          cs%uhGM_Rlay(i,j,k) = 0.0
        enddo ; enddo
        do k=nkmb+1,nz ; do i=isq,ieq
          cs%uhGM_Rlay(i,j,k) = cdp%uhGM(i,j,k)
        enddo ; enddo
        do k=1,nkmb ; do i=isq,ieq
          call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i+1,j,k)), k_list, nz, wt, wt_p)
          cs%uhGM_Rlay(i,j,k_list)   = cs%uhGM_Rlay(i,j,k_list)   + cdp%uhGM(i,j,k)*wt
          cs%uhGM_Rlay(i,j,k_list+1) = cs%uhGM_Rlay(i,j,k_list+1) + cdp%uhGM(i,j,k)*wt_p
        enddo ; enddo
      enddo

      if (cs%id_uh_Rlay > 0) call post_data(cs%id_uhGM_Rlay, cs%uhGM_Rlay, cs%diag)
    endif

    if (ASSOCIATED(cs%vhGM_Rlay) .and. ASSOCIATED(cdp%vhGM)) then
      k_list = nz/2
!$OMP parallel do default(none) shared(is,ie,Jsq,Jeq,nz,nkmb,CS,CDp,Rcv,GV) &
!$OMP                          private(wt,wt_p) firstprivate(k_list)
      do j=jsq,jeq
        do k=1,nkmb ; do i=is,ie
          cs%vhGM_Rlay(i,j,k) = 0.0
        enddo ; enddo
        do k=nkmb+1,nz ; do i=is,ie
          cs%vhGM_Rlay(i,j,k) = cdp%vhGM(i,j,k)
        enddo ; enddo
        do k=1,nkmb ; do i=is,ie
          call find_weights(gv%Rlay, 0.5*(rcv(i,j,k)+rcv(i,j+1,k)), k_list, nz, wt, wt_p)
          cs%vhGM_Rlay(i,j,k_list)   = cs%vhGM_Rlay(i,j,k_list)   + cdp%vhGM(i,j,k)*wt
          cs%vhGM_Rlay(i,j,k_list+1) = cs%vhGM_Rlay(i,j,k_list+1) + cdp%vhGM(i,j,k)*wt_p
        enddo ; enddo
      enddo

      if (cs%id_vhGM_Rlay > 0) call post_data(cs%id_vhGM_Rlay, cs%vhGM_Rlay, cs%diag)
    endif
  endif

  if (associated(tv%eqn_of_state)) then
    if (cs%id_rhopot0 > 0) then
      pressure_1d(:) = 0.
!$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d)
      do k=1,nz ; do j=js,je
        call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
                               rcv(:,j,k),is,ie-is+1, tv%eqn_of_state)
      enddo ; enddo
      if (cs%id_rhopot0 > 0) call post_data(cs%id_rhopot0, rcv, cs%diag)
    endif
    if (cs%id_rhopot2 > 0) then
      pressure_1d(:) = 2.e7 ! 2000 dbars
!$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d)
      do k=1,nz ; do j=js,je
        call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
                               rcv(:,j,k),is,ie-is+1, tv%eqn_of_state)
      enddo ; enddo
      if (cs%id_rhopot2 > 0) call post_data(cs%id_rhopot2, rcv, cs%diag)
    endif
    if (cs%id_rhoinsitu > 0) then
!$OMP parallel do default(none) shared(tv,Rcv,is,ie,js,je,nz,pressure_1d,h,GV)
      do j=js,je
        pressure_1d(:) = 0. ! Start at p=0 Pa at surface
        do k=1,nz
          pressure_1d(:) =  pressure_1d(:) + 0.5 * h(:,j,k) * gv%H_to_Pa ! Pressure in middle of layer k
          call calculate_density(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
                                 rcv(:,j,k),is,ie-is+1, tv%eqn_of_state)
          pressure_1d(:) =  pressure_1d(:) + 0.5 * h(:,j,k) * gv%H_to_Pa ! Pressure at bottom of layer k
        enddo
      enddo
      if (cs%id_rhoinsitu > 0) call post_data(cs%id_rhoinsitu, rcv, cs%diag)
    endif
  endif

  if ((cs%id_cg1>0) .or. (cs%id_Rd1>0) .or. (cs%id_cfl_cg1>0) .or. &
      (cs%id_cfl_cg1_x>0) .or. (cs%id_cfl_cg1_y>0)) then
    call wave_speed(h, tv, g, gv, cs%cg1, cs%wave_speed_CSp)
    if (cs%id_cg1>0) call post_data(cs%id_cg1, cs%cg1, cs%diag)
    if (cs%id_Rd1>0) then
!$OMP parallel do default(none) shared(is,ie,js,je,G,CS) &
!$OMP                          private(f2_h,mag_beta)
      do j=js,je ; do i=is,ie
        ! Blend the equatorial deformation radius with the standard one.
        f2_h = absurdly_small_freq2 + 0.25 * &
            ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
             (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
        mag_beta = sqrt(0.5 * ( &
            (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
             ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
            (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
             ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ))
        cs%Rd1(i,j) = cs%cg1(i,j) / sqrt(f2_h + cs%cg1(i,j) * mag_beta)

      enddo ; enddo
      call post_data(cs%id_Rd1, cs%Rd1, cs%diag)
    endif
    if (cs%id_cfl_cg1>0) then
      do j=js,je ; do i=is,ie
        cs%cfl_cg1(i,j) = (dt*cs%cg1(i,j)) * (g%IdxT(i,j) + g%IdyT(i,j))
      enddo ; enddo
      call post_data(cs%id_cfl_cg1, cs%cfl_cg1, cs%diag)
    endif
    if (cs%id_cfl_cg1_x>0) then
      do j=js,je ; do i=is,ie
        cs%cfl_cg1_x(i,j) = (dt*cs%cg1(i,j)) * g%IdxT(i,j)
      enddo ; enddo
      call post_data(cs%id_cfl_cg1_x, cs%cfl_cg1_x, cs%diag)
    endif
    if (cs%id_cfl_cg1_y>0) then
      do j=js,je ; do i=is,ie
        cs%cfl_cg1_y(i,j) = (dt*cs%cg1(i,j)) * g%IdyT(i,j)
      enddo ; enddo
      call post_data(cs%id_cfl_cg1_y, cs%cfl_cg1_y, cs%diag)
    endif
  endif
  if ((cs%id_cg_ebt>0) .or. (cs%id_Rd_ebt>0) .or. (cs%id_p_ebt>0)) then
    if (cs%id_p_ebt>0) then
      call wave_speed(h, tv, g, gv, cs%cg1, cs%wave_speed_CSp, use_ebt_mode=.true., &
                      mono_n2_column_fraction=cs%mono_N2_column_fraction, &
                      mono_n2_depth=cs%mono_N2_depth, modal_structure=cs%p_ebt)
      call post_data(cs%id_p_ebt, cs%p_ebt, cs%diag)
    else
      call wave_speed(h, tv, g, gv, cs%cg1, cs%wave_speed_CSp, use_ebt_mode=.true., &
                      mono_n2_column_fraction=cs%mono_N2_column_fraction, &
                      mono_n2_depth=cs%mono_N2_depth)
    endif
    if (cs%id_cg_ebt>0) call post_data(cs%id_cg_ebt, cs%cg1, cs%diag)
    if (cs%id_Rd_ebt>0) then
!$OMP parallel do default(none) shared(is,ie,js,je,G,CS) &
!$OMP                          private(f2_h,mag_beta)
      do j=js,je ; do i=is,ie
        ! Blend the equatorial deformation radius with the standard one.
        f2_h = absurdly_small_freq2 + 0.25 * &
            ((g%CoriolisBu(i,j)**2 + g%CoriolisBu(i-1,j-1)**2) + &
             (g%CoriolisBu(i-1,j)**2 + g%CoriolisBu(i,j-1)**2))
        mag_beta = sqrt(0.5 * ( &
            (((g%CoriolisBu(i,j)-g%CoriolisBu(i-1,j)) * g%IdxCv(i,j))**2 + &
             ((g%CoriolisBu(i,j-1)-g%CoriolisBu(i-1,j-1)) * g%IdxCv(i,j-1))**2) + &
            (((g%CoriolisBu(i,j)-g%CoriolisBu(i,j-1)) * g%IdyCu(i,j))**2 + &
             ((g%CoriolisBu(i-1,j)-g%CoriolisBu(i-1,j-1)) * g%IdyCu(i-1,j))**2) ))
        cs%Rd1(i,j) = cs%cg1(i,j) / sqrt(f2_h + cs%cg1(i,j) * mag_beta)

      enddo ; enddo
      call post_data(cs%id_Rd_ebt, cs%Rd1, cs%diag)
    endif
  endif

  if (dt > 0.0) then
    if (cs%id_du_dt>0) call post_data(cs%id_du_dt, cs%du_dt, cs%diag)

    if (cs%id_dv_dt>0) call post_data(cs%id_dv_dt, cs%dv_dt, cs%diag)

    if (cs%id_dh_dt>0) call post_data(cs%id_dh_dt, cs%dh_dt, cs%diag)

    call calculate_energy_diagnostics(u, v, h, uh, vh, adp, cdp, g, cs)
  endif

end subroutine calculate_diagnostic_fields

!> This subroutine finds location of R_in in an increasing ordered
!! list, Rlist, returning as k the element such that
!! Rlist(k) <= R_in < Rlist(k+1), and where wt and wt_p are the linear
!! weights that should be assigned to elements k and k+1.
subroutine find_weights(Rlist, R_in, k, nz, wt, wt_p)
  real,     intent(in)    :: Rlist(:), R_in
  integer,  intent(inout) :: k
  integer,  intent(in)    :: nz
  real,     intent(out)   :: wt, wt_p

  ! This subroutine finds location of R_in in an increasing ordered
  ! list, Rlist, returning as k the element such that
  !  Rlist(k) <= R_in < Rlist(k+1), and where wt and wt_p are the linear
  ! weights that should be assigned to elements k and k+1.

  integer :: k_upper, k_lower, k_new, inc

  ! First, bracket the desired point.
  if ((k < 1) .or. (k > nz)) k = nz/2

  k_upper = k ; k_lower = k ; inc = 1
  if (r_in < rlist(k)) then
    do
      k_lower = max(k_lower-inc, 1)
      if ((k_lower == 1) .or. (r_in >= rlist(k_lower))) exit
      k_upper = k_lower
      inc = inc*2
    end do
  else
    do
      k_upper = min(k_upper+inc, nz)
      if ((k_upper == nz) .or. (r_in < rlist(k_upper))) exit
      k_lower = k_upper
      inc = inc*2
    end do
  endif

  if ((k_lower == 1) .and. (r_in <= rlist(k_lower))) then
    k = 1 ; wt = 1.0 ; wt_p = 0.0
  else if ((k_upper == nz) .and. (r_in >= rlist(k_upper))) then
    k = nz-1 ; wt = 0.0 ; wt_p = 1.0
  else
    do
      if (k_upper <= k_lower+1) exit
      k_new = (k_upper + k_lower) / 2
      if (r_in < rlist(k_new)) then
        k_upper = k_new
      else
        k_lower = k_new
      endif
    end do

!   Uncomment this as a code check
!    if ((R_in < Rlist(k_lower)) .or. (R_in >= Rlist(k_upper)) .or. (k_upper-k_lower /= 1)) &
!      write (*,*) "Error: ",R_in," is not between R(",k_lower,") = ", &
!        Rlist(k_lower)," and R(",k_upper,") = ",Rlist(k_upper),"."
    k = k_lower
    wt = (rlist(k_upper) - r_in) / (rlist(k_upper) - rlist(k_lower))
    wt_p = 1.0 - wt

  endif

end subroutine find_weights

!> Subroutine calculates vertical integrals of several tracers, along
!! with the mass-weight of these tracers, the total column mass, and the
!! carefully calculated column height.
subroutine calculate_vertical_integrals(h, tv, fluxes, G, GV, CS)
  type(ocean_grid_type),                    intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),                  intent(in)    :: GV   !< The ocean's vertical grid
                                                                  !! structure.
  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(forcing),                            intent(in)    :: fluxes !< A structure containing the
                                                                  !! surface fluxes.
  type(diagnostics_cs),                     intent(inout) :: CS   !< A control structure returned
                                                                  !! by a previous call to
                                                                  !! diagnostics_init.

! Subroutine calculates vertical integrals of several tracers, along
! with the mass-weight of these tracers, the total column mass, and the
! carefully calculated column height.

! Arguments:
!  (in)      h  - layer thickness: metre (Bouss) or kg/ m2 (non-Bouss)
!  (in)      tv - structure pointing to thermodynamic variables
!  (in)      fluxes - a structure containing the surface fluxes.
!  (in)      G  - ocean grid structure
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      CS - control structure returned by a previous call to diagnostics_init

  real, dimension(SZI_(G), SZJ_(G)) :: &
    z_top, &  ! Height of the top of a layer or the ocean, in m.
    z_bot, &  ! Height of the bottom of a layer (for id_mass) or the
              ! (positive) depth of the ocean (for id_col_ht), in m.
    mass, &   ! integrated mass of the water column, in kg m-2.  For
              ! non-Boussinesq models this is rho*dz. For Boussiensq
              ! models, this is either the integral of in-situ density
              ! (rho*dz for col_mass) or reference dens (Rho_0*dz for mass_wt).
    btm_pres,&! The pressure at the ocean bottom, or CMIP variable 'pbo'.
              ! This is the column mass multiplied by gravity plus the pressure
              ! at the ocean surface.
    dpress, &    ! Change in hydrostatic pressure across a layer, in Pa.
    tr_int    ! vertical integral of a tracer times density,
              ! (Rho_0 in a Boussinesq model) in TR kg m-2.
  real    :: IG_Earth  ! Inverse of gravitational acceleration, in s2 m-1.

  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  if (cs%id_mass_wt > 0) then
    do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo
    do k=1,nz ; do j=js,je ; do i=is,ie
      mass(i,j) = mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_mass_wt, mass, cs%diag)
  endif

  if (cs%id_temp_int > 0) then
    do j=js,je ; do i=is,ie ; tr_int(i,j) = 0.0 ; enddo ; enddo
    do k=1,nz ; do j=js,je ; do i=is,ie
      tr_int(i,j) = tr_int(i,j) + (gv%H_to_kg_m2*h(i,j,k))*tv%T(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_temp_int, tr_int, cs%diag)
  endif

  if (cs%id_salt_int > 0) then
    do j=js,je ; do i=is,ie ; tr_int(i,j) = 0.0 ; enddo ; enddo
    do k=1,nz ; do j=js,je ; do i=is,ie
      tr_int(i,j) = tr_int(i,j) + (gv%H_to_kg_m2*h(i,j,k))*tv%S(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_salt_int, tr_int, cs%diag)
  endif

  if (cs%id_col_ht > 0) then
    call find_eta(h, tv, gv%g_Earth, g, gv, z_top)
    do j=js,je ; do i=is,ie
      z_bot(i,j) = z_top(i,j) + g%bathyT(i,j)
    enddo ; enddo
    call post_data(cs%id_col_ht, z_bot, cs%diag)
  endif

  if (cs%id_col_mass > 0 .or. cs%id_pbo > 0) then
    do j=js,je ; do i=is,ie ; mass(i,j) = 0.0 ; enddo ; enddo
    if (gv%Boussinesq) then
      if (associated(tv%eqn_of_state)) then
        ig_earth = 1.0 / gv%g_Earth
!       do j=js,je ; do i=is,ie ; z_bot(i,j) = -P_SURF(i,j)/GV%H_to_Pa ; enddo ; enddo
        do j=js,je ; do i=is,ie ; z_bot(i,j) = 0.0 ; enddo ; enddo
        do k=1,nz
          do j=js,je ; do i=is,ie
            z_top(i,j) = z_bot(i,j)
            z_bot(i,j) = z_top(i,j) - gv%H_to_m*h(i,j,k)
          enddo ; enddo
          call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), &
                              z_top, z_bot, 0.0, gv%H_to_kg_m2, gv%g_Earth, &
                              g%HI, g%HI, tv%eqn_of_state, dpress)
          do j=js,je ; do i=is,ie
            mass(i,j) = mass(i,j) + dpress(i,j) * ig_earth
          enddo ; enddo
        enddo
      else
        do k=1,nz ; do j=js,je ; do i=is,ie
          mass(i,j) = mass(i,j) + (gv%H_to_m*gv%Rlay(k))*h(i,j,k)
        enddo ; enddo ; enddo
      endif
    else
      do k=1,nz ; do j=js,je ; do i=is,ie
        mass(i,j) = mass(i,j) + gv%H_to_kg_m2*h(i,j,k)
      enddo ; enddo ; enddo
    endif
    if (cs%id_col_mass > 0) then
      call post_data(cs%id_col_mass, mass, cs%diag)
    endif
    if (cs%id_pbo > 0) then
      do j=js,je ; do i=is,ie ; btm_pres(i,j) = 0.0 ; enddo ; enddo
      ! 'pbo' is defined as the sea water pressure at the sea floor
      !     pbo = (mass * g) + pso
      ! where pso is the sea water pressure at sea water surface
      ! note that pso is equivalent to fluxes%p_surf
      do j=js,je ; do i=is,ie
        btm_pres(i,j) = mass(i,j) * gv%g_Earth
        if (ASSOCIATED(fluxes%p_surf)) then
          btm_pres(i,j) = btm_pres(i,j) + fluxes%p_surf(i,j)
        endif
      enddo ; enddo
      call post_data(cs%id_pbo, btm_pres, cs%diag)
    endif
  endif

end subroutine calculate_vertical_integrals

!> This subroutine calculates terms in the mechanical energy budget.
subroutine calculate_energy_diagnostics(u, v, h, uh, vh, ADp, CDp, G, CS)
  type(ocean_grid_type),                     intent(inout) :: G    !< The ocean's 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).
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: uh   !< Transport through zonal
                                                                   !! faces=u*h*dy: m3/s (Bouss)
                                                                   !! kg/s(non-Bouss).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: vh   !< Transport through merid
                                                                   !! faces=v*h*dx: m3/s (Bouss)
                                                                   !! kg/s(non-Bouss).
  type(accel_diag_ptrs),                     intent(in)    :: ADp  !< Structure pointing to
                                                                   !! accelerations in momentum
                                                                   !! equation.
  type(cont_diag_ptrs),                      intent(in)    :: CDp  !< Structure pointing to terms
                                                                   !! in continuity equations.
  type(diagnostics_cs),                      intent(inout) :: CS   !< Control structure returned by
                                                                   !! a previous call to
                                                                   !! diagnostics_init.

! This subroutine calculates terms in the mechanical energy budget.

! Arguments:
!  (in)      u   - zonal velocity component (m/s)
!  (in)      v   - meridional velocity componnent (m/s)
!  (in)      h   - layer thickness: metre(Bouss) of kg/m2(non-Bouss)
!  (in)      uh  - transport through zonal faces=u*h*dy: m3/s (Bouss) kg/s(non-Bouss)
!  (in)      vh  - transport through merid faces=v*h*dx: m3/s (Bouss) kg/s(non-Bouss)
!  (in)      ADp - structure pointing to accelerations in momentum equation
!  (in)      CDp - structure pointing to terms in continuity equations
!  (in)      G   - ocean grid structure
!  (in)      CS  - control structure returned by a previous call to diagnostics_init

  real :: KE_u(szib_(g),szj_(g))
  real :: KE_v(szi_(g),szjb_(g))
  real :: KE_h(szi_(g),szj_(g))

  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nz
  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

  do j=js-1,je ; do i=is-1,ie
    ke_u(i,j) = 0.0 ; ke_v(i,j) = 0.0
  enddo ; enddo

  if (ASSOCIATED(cs%KE)) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      cs%KE(i,j,k) = ((u(i,j,k)*u(i,j,k) + u(i-1,j,k)*u(i-1,j,k)) + &
          (v(i,j,k)*v(i,j,k) + v(i,j-1,k)*v(i,j-1,k)))*0.25
      ! DELETE THE FOLLOWING...  Make this 0 to test the momentum balance,
      ! or a huge number to test the continuity balance.
      ! CS%KE(i,j,k) *= 1e20
    enddo ; enddo ; enddo
    if (cs%id_KE > 0) call post_data(cs%id_KE, cs%KE, cs%diag)
  endif

  if(.not.g%symmetric) then
    if(ASSOCIATED(cs%dKE_dt) .OR. ASSOCIATED(cs%PE_to_KE) .OR. ASSOCIATED(cs%KE_CorAdv) .OR. &
       ASSOCIATED(cs%KE_adv) .OR. ASSOCIATED(cs%KE_visc)  .OR. ASSOCIATED(cs%KE_horvisc).OR. &
       ASSOCIATED(cs%KE_dia) ) then
        call create_group_pass(cs%pass_KE_uv, ke_u, ke_v, g%Domain, to_north+to_east)
    endif
  endif

  if (ASSOCIATED(cs%dKE_dt)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*cs%du_dt(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*cs%dv_dt(i,j,k)
      enddo ; enddo
      do j=js,je ; do i=is,ie
        ke_h(i,j) = cs%KE(i,j,k)*cs%dh_dt(i,j,k)
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%dKE_dt(i,j,k) = ke_h(i,j) + 0.5 * g%IareaT(i,j) * &
            (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_dKEdt > 0) call post_data(cs%id_dKEdt, cs%dKE_dt, cs%diag)
  endif

  if (ASSOCIATED(cs%PE_to_KE)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%PFu(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%PFv(i,j,k)
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%PE_to_KE(i,j,k) = 0.5 * g%IareaT(i,j) * &
            (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_PE_to_KE > 0) call post_data(cs%id_PE_to_KE, cs%PE_to_KE, cs%diag)
  endif

  if (ASSOCIATED(cs%KE_CorAdv)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%CAu(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%CAv(i,j,k)
      enddo ; enddo
      do j=js,je ; do i=is,ie
        ke_h(i,j) = -cs%KE(i,j,k) * g%IareaT(i,j) * &
            (uh(i,j,k) - uh(i-1,j,k) + vh(i,j,k) - vh(i,j-1,k))
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%KE_CorAdv(i,j,k) = ke_h(i,j) + 0.5 * g%IareaT(i,j) * &
            (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_KE_Coradv > 0) call post_data(cs%id_KE_Coradv, cs%KE_Coradv, cs%diag)
  endif

  if (ASSOCIATED(cs%KE_adv)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%gradKEu(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%gradKEv(i,j,k)
      enddo ; enddo
      do j=js,je ; do i=is,ie
        ke_h(i,j) = -cs%KE(i,j,k) * g%IareaT(i,j) * &
            (uh(i,j,k) - uh(i-1,j,k) + vh(i,j,k) - vh(i,j-1,k))
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%KE_adv(i,j,k) = ke_h(i,j) + 0.5 * g%IareaT(i,j) * &
            (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_KE_adv > 0) call post_data(cs%id_KE_adv, cs%KE_adv, cs%diag)
  endif

  if (ASSOCIATED(cs%KE_visc)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%du_dt_visc(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%dv_dt_visc(i,j,k)
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%KE_visc(i,j,k) = 0.5 * g%IareaT(i,j) * &
            (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_KE_visc > 0) call post_data(cs%id_KE_visc, cs%KE_visc, cs%diag)
  endif

  if (ASSOCIATED(cs%KE_horvisc)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%diffu(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%diffv(i,j,k)
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%KE_horvisc(i,j,k) = 0.5 * g%IareaT(i,j) * &
            (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_KE_horvisc > 0) call post_data(cs%id_KE_horvisc, cs%KE_horvisc, cs%diag)
  endif

  if (ASSOCIATED(cs%KE_dia)) then
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        ke_u(i,j) = uh(i,j,k)*g%dxCu(i,j)*adp%du_dt_dia(i,j,k)
      enddo ; enddo
      do j=jsq,jeq ; do i=is,ie
        ke_v(i,j) = vh(i,j,k)*g%dyCv(i,j)*adp%dv_dt_dia(i,j,k)
      enddo ; enddo
      do j=js,je ; do i=is,ie
        ke_h(i,j) = cs%KE(i,j,k) * &
            (cdp%diapyc_vel(i,j,k) - cdp%diapyc_vel(i,j,k+1))
      enddo ; enddo
      if (.not.g%symmetric) &
         call do_group_pass(cs%pass_KE_uv, g%domain)
      do j=js,je ; do i=is,ie
        cs%KE_dia(i,j,k) = ke_h(i,j) + 0.5 * g%IareaT(i,j) * &
            (ke_u(i,j) + ke_u(i-1,j) + ke_v(i,j) + ke_v(i,j-1))
      enddo ; enddo
    enddo
    if (cs%id_KE_dia > 0) call post_data(cs%id_KE_dia, cs%KE_dia, cs%diag)
  endif

end subroutine calculate_energy_diagnostics

!> This subroutine registers fields to calculate a diagnostic time derivative.
subroutine register_time_deriv(f_ptr, deriv_ptr, CS)
  real, dimension(:,:,:), target :: f_ptr     !< Field whose derivative is taken.
  real, dimension(:,:,:), target :: deriv_ptr !< Field in which the calculated time derivatives
                                              !! placed.
  type(diagnostics_cs),  pointer :: CS        !< Control structure returned by previous call to
                                              !! diagnostics_init.

! This subroutine registers fields to calculate a diagnostic time derivative.
! Arguments:
!  (target)  f_ptr     - field whose derivative is taken
!  (in)      deriv_ptr - field in which the calculated time derivatives placed
!  (in)      num_lay   - number of layers in this field
!  (in)      CS        - control structure returned by previous call to diagnostics_init

  integer :: m

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

  if (cs%num_time_deriv >= max_fields_) then
    call mom_error(warning,"MOM_diagnostics:  Attempted to register more than " // &
                   "MAX_FIELDS_ diagnostic time derivatives via register_time_deriv.")
    return
  endif

  m = cs%num_time_deriv+1 ; cs%num_time_deriv = m

  cs%nlay(m) = size(f_ptr(:,:,:),3)
  cs%deriv(m)%p => deriv_ptr
  allocate(cs%prev_val(m)%p(size(f_ptr(:,:,:),1), size(f_ptr(:,:,:),2), cs%nlay(m)) )

  cs%var_ptr(m)%p => f_ptr
  cs%prev_val(m)%p(:,:,:) = f_ptr(:,:,:)

end subroutine register_time_deriv

!> This subroutine calculates all registered time derivatives.
subroutine calculate_derivs(dt, G, CS)
  real,                  intent(in)    :: dt   !< The time interval over which differences occur,
                                               !! in s.
  type(ocean_grid_type), intent(inout) :: G    !< The ocean's grid structure.
  type(diagnostics_cs),  intent(inout) :: CS   !< Control structure returned by previous call to
                                               !! diagnostics_init.

! This subroutine calculates all registered time derivatives.
! Arguments:
!  (in)      dt - time interval in s over which differences occur
!  (in)      G  - ocean grid structure.
!  (in)      CS - control structure returned by previous call to diagnostics_init

  integer i, j, k, m
  real Idt

  if (dt > 0.0) then ; idt = 1.0/dt
  else ; return ; endif

  do m=1,cs%num_time_deriv
    do k=1,cs%nlay(m) ; do j=g%jsc,g%jec ; do i=g%isc,g%iec
      cs%deriv(m)%p(i,j,k) = (cs%var_ptr(m)%p(i,j,k) - cs%prev_val(m)%p(i,j,k)) * idt
      cs%prev_val(m)%p(i,j,k) = cs%var_ptr(m)%p(i,j,k)
    enddo ; enddo ; enddo
  enddo

end subroutine calculate_derivs

! #@# This subroutine needs a doxygen description
subroutine mom_diagnostics_init(MIS, ADp, CDp, Time, G, GV, param_file, diag, CS)
  type(ocean_internal_state), intent(in)    :: MIS  !< For "MOM Internal State" a set of pointers to
                                                    !! the fields and accelerations that make up the
                                                    !! ocean's internal physical state.
  type(accel_diag_ptrs),      intent(inout) :: ADp  !< Structure with pointers to momentum equation
                                                    !! terms.
  type(cont_diag_ptrs),       intent(inout) :: CDp  !< Structure with pointers to continuity
                                                    !! equation terms.
  type(time_type),            intent(in)    :: Time !< Current model time.
  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(param_file_type),      intent(in)    :: param_file !< A structure to parse for run-time
                                                    !! parameters.
  type(diag_ctrl), target,    intent(inout) :: diag !< Structure to regulate diagnostic output.
  type(diagnostics_cs),       pointer       :: CS   !< Pointer set to point to control structure
                                                    !! for this module.

! Arguments
!  (in)     MIS    - For "MOM Internal State" a set of pointers to the fields and
!                    accelerations that make up the ocean's internal physical
!                    state.
!  (inout)  ADp    - structure with pointers to momentum equation terms
!  (inout)  CDp    - structure with pointers to continuity equation terms
!  (in)     Time   - current model time
!  (in)     G      - ocean grid structure
!  (in)     GV     - The ocean's vertical grid structure.
!  (in) param_file - structure indicating the open file to parse for
!                     model parameter values
!  (in)     diag   - structure to regulate diagnostic output
!  (in/out) CS     - pointer set to point to control structure for this module

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

  character(len=40)  :: mdl = "MOM_diagnostics" ! This module's name.
  real :: omega, f2_min
  character(len=48) :: thickness_units, flux_units
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz, nkml, nkbl
  integer :: is, ie, js, je, Isq, Ieq, Jsq, Jeq, i, j

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

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

  cs%diag => diag

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version)
  call get_param(param_file, mdl, "DIAG_EBT_MONO_N2_COLUMN_FRACTION", cs%mono_N2_column_fraction, &
                 "The lower fraction of water column over which N2 is limited as monotonic\n"// &
                 "for the purposes of calculating the equivalent barotropic wave speed.", &
                 units='nondim', default=0.)
  call get_param(param_file, mdl, "DIAG_EBT_MONO_N2_DEPTH", cs%mono_N2_depth, &
                 "The depth below which N2 is limited as monotonic for the\n"// &
                 "purposes of calculating the equivalent barotropic wave speed.", &
                 units='m', default=-1.)

  if (gv%Boussinesq) then
    thickness_units = "meter" ; flux_units = "meter3 second-1"
  else
    thickness_units = "kilogram meter-2" ; flux_units = "kilogram second-1"
  endif

  cs%id_temp_layer_ave = register_diag_field('ocean_model', 'temp_layer_ave', diag%axesZL, time, &
      'Layer Average Ocean Temperature', 'Celsius')

  cs%id_salt_layer_ave = register_diag_field('ocean_model', 'salt_layer_ave', diag%axesZL, time, &
      'Layer Average Ocean Salinity', 'ppt')

  cs%id_masscello = register_diag_field('ocean_model', 'masscello', diag%axesTL,&
      time, 'Mass per unit area of liquid ocean grid cell', 'kg m-2',           &
      standard_name='sea_water_mass_per_unit_area', v_extensive=.true.)

  cs%id_masso = register_scalar_field('ocean_model', 'masso', time,  &
      diag, 'Mass of liquid ocean', 'kg', standard_name='sea_water_mass')

  cs%id_thkcello = register_diag_field('ocean_model', 'thkcello', diag%axesTL, time, &
      long_name = 'Cell Thickness', standard_name='cell_thickness', units='m', v_extensive=.true.)

  if (((cs%id_masscello>0) .or. (cs%id_masso>0) .or. (cs%id_thkcello>0.and..not.gv%Boussinesq)) &
      .and. .not.ASSOCIATED(cs%diag_tmp3d)) then
    call safe_alloc_ptr(cs%diag_tmp3d,isd,ied,jsd,jed,nz)
  endif

  cs%id_thetaoga = register_scalar_field('ocean_model', 'thetaoga',    &
      time, diag, 'Global Mean Ocean Potential Temperature', 'Celsius',&
      standard_name='sea_water_potential_temperature')

  cs%id_soga = register_scalar_field('ocean_model', 'soga', &
      time, diag, 'Global Mean Ocean Salinity', 'ppt',      &
      standard_name='sea_water_salinity')

  cs%id_tosga = register_scalar_field('ocean_model', 'sst_global', time, diag,&
      long_name='Global Area Average Sea Surface Temperature',                &
      units='degC', standard_name='sea_surface_temperature',                  &
      cmor_field_name='tosga', cmor_standard_name='sea_surface_temperature',  &
      cmor_units='degC', cmor_long_name='Sea Surface Temperature')

  cs%id_sosga = register_scalar_field('ocean_model', 'sss_global', time, diag,&
      long_name='Global Area Average Sea Surface Salinity',                   &
      units='ppt', standard_name='sea_surface_salinity',                      &
      cmor_field_name='sosga', cmor_standard_name='sea_surface_salinity',     &
      cmor_units='ppt', cmor_long_name='Sea Surface Salinity')

  cs%id_e = register_diag_field('ocean_model', 'e', diag%axesTi, time, &
      'Interface Height Relative to Mean Sea Level', 'meter')
  if (cs%id_e>0) call safe_alloc_ptr(cs%e,isd,ied,jsd,jed,nz+1)

  cs%id_e_D = register_diag_field('ocean_model', 'e_D', diag%axesTi, time, &
      'Interface Height above the Seafloor', 'meter')
  if (cs%id_e_D>0) call safe_alloc_ptr(cs%e_D,isd,ied,jsd,jed,nz+1)

  cs%id_Rml = register_diag_field('ocean_model', 'Rml', diag%axesTL, time, &
      'Mixed Layer Coordinate Potential Density', 'kg meter-3')

  cs%id_Rcv = register_diag_field('ocean_model', 'Rho_cv', diag%axesTL, time, &
      'Coordinate Potential Density', 'kg meter-3')

  cs%id_rhopot0 = register_diag_field('ocean_model', 'rhopot0', diag%axesTL, time, &
      'Potential density referenced to surface', 'kg meter-3')
  cs%id_rhopot2 = register_diag_field('ocean_model', 'rhopot2', diag%axesTL, time, &
      'Potential density referenced to 2000 dbar', 'kg meter-3')
  cs%id_rhoinsitu = register_diag_field('ocean_model', 'rhoinsitu', diag%axesTL, time, &
      'In situ density', 'kg meter-3')

  cs%id_du_dt = register_diag_field('ocean_model', 'dudt', diag%axesCuL, time, &
      'Zonal Acceleration', 'meter second-2')
  if ((cs%id_du_dt>0) .and. .not.ASSOCIATED(cs%du_dt)) then
    call safe_alloc_ptr(cs%du_dt,isdb,iedb,jsd,jed,nz)
    call register_time_deriv(mis%u, cs%du_dt, cs)
  endif

  cs%id_dv_dt = register_diag_field('ocean_model', 'dvdt', diag%axesCvL, time, &
      'Meridional Acceleration', 'meter second-2')
  if ((cs%id_dv_dt>0) .and. .not.ASSOCIATED(cs%dv_dt)) then
    call safe_alloc_ptr(cs%dv_dt,isd,ied,jsdb,jedb,nz)
    call register_time_deriv(mis%v, cs%dv_dt, cs)
  endif

  cs%id_dh_dt = register_diag_field('ocean_model', 'dhdt', diag%axesTL, time, &
      'Thickness tendency', trim(thickness_units)//" second-1")
  if ((cs%id_dh_dt>0) .and. .not.ASSOCIATED(cs%dh_dt)) then
    call safe_alloc_ptr(cs%dh_dt,isd,ied,jsd,jed,nz)
    call register_time_deriv(mis%h, cs%dh_dt, cs)
  endif

  ! layer thickness variables
  !if (GV%nk_rho_varies > 0) then
    cs%id_h_Rlay = register_diag_field('ocean_model', 'h_rho', diag%axesTL, time, &
        'Layer thicknesses in pure potential density coordinates', thickness_units)
    if (cs%id_h_Rlay>0) call safe_alloc_ptr(cs%h_Rlay,isd,ied,jsd,jed,nz)

    cs%id_uh_Rlay = register_diag_field('ocean_model', 'uh_rho', diag%axesCuL, time, &
        'Zonal volume transport in pure potential density coordinates', flux_units)
    if (cs%id_uh_Rlay>0) call safe_alloc_ptr(cs%uh_Rlay,isdb,iedb,jsd,jed,nz)

    cs%id_vh_Rlay = register_diag_field('ocean_model', 'vh_rho', diag%axesCvL, time, &
        'Meridional volume transport in pure potential density coordinates', flux_units)
    if (cs%id_vh_Rlay>0) call safe_alloc_ptr(cs%vh_Rlay,isd,ied,jsdb,jedb,nz)

    cs%id_uhGM_Rlay = register_diag_field('ocean_model', 'uhGM_rho', diag%axesCuL, time, &
        'Zonal volume transport due to interface height diffusion in pure potential &
        &density coordinates', flux_units)
    if (cs%id_uhGM_Rlay>0) call safe_alloc_ptr(cs%uhGM_Rlay,isdb,iedb,jsd,jed,nz)

    cs%id_vhGM_Rlay = register_diag_field('ocean_model', 'vhGM_rho', diag%axesCvL, time, &
        'Meridional volume transport due to interface height diffusion in pure &
        &potential density coordinates', flux_units)
    if (cs%id_vhGM_Rlay>0) call safe_alloc_ptr(cs%vhGM_Rlay,isd,ied,jsdb,jedb,nz)
  !endif


  ! terms in the kinetic energy budget
  cs%id_KE = register_diag_field('ocean_model', 'KE', diag%axesTL, time, &
      'Layer kinetic energy per unit mass', 'meter2 second-2')
  if (cs%id_KE>0) call safe_alloc_ptr(cs%KE,isd,ied,jsd,jed,nz)

  cs%id_dKEdt = register_diag_field('ocean_model', 'dKE_dt', diag%axesTL, time, &
      'Kinetic Energy Tendency of Layer', 'meter3 second-3')
  if (cs%id_dKEdt>0) call safe_alloc_ptr(cs%dKE_dt,isd,ied,jsd,jed,nz)

  cs%id_PE_to_KE = register_diag_field('ocean_model', 'PE_to_KE', diag%axesTL, time, &
      'Potential to Kinetic Energy Conversion of Layer', 'meter3 second-3')
  if (cs%id_PE_to_KE>0) call safe_alloc_ptr(cs%PE_to_KE,isd,ied,jsd,jed,nz)

  cs%id_KE_Coradv = register_diag_field('ocean_model', 'KE_Coradv', diag%axesTL, time, &
      'Kinetic Energy Source from Coriolis and Advection', 'meter3 second-3')
  if (cs%id_KE_Coradv>0) call safe_alloc_ptr(cs%KE_Coradv,isd,ied,jsd,jed,nz)

  cs%id_KE_adv = register_diag_field('ocean_model', 'KE_adv', diag%axesTL, time, &
      'Kinetic Energy Source from Advection', 'meter3 second-3')
  if (cs%id_KE_adv>0) call safe_alloc_ptr(cs%KE_adv,isd,ied,jsd,jed,nz)

  cs%id_KE_visc = register_diag_field('ocean_model', 'KE_visc', diag%axesTL, time, &
      'Kinetic Energy Source from Vertical Viscosity and Stresses', 'meter3 second-3')
  if (cs%id_KE_visc>0) call safe_alloc_ptr(cs%KE_visc,isd,ied,jsd,jed,nz)

  cs%id_KE_horvisc = register_diag_field('ocean_model', 'KE_horvisc', diag%axesTL, time, &
      'Kinetic Energy Source from Horizontal Viscosity', 'meter3 second-3')
  if (cs%id_KE_horvisc>0) call safe_alloc_ptr(cs%KE_horvisc,isd,ied,jsd,jed,nz)

  cs%id_KE_dia = register_diag_field('ocean_model', 'KE_dia', diag%axesTL, time, &
      'Kinetic Energy Source from Diapycnal Diffusion', 'meter3 second-3')
  if (cs%id_KE_dia>0) call safe_alloc_ptr(cs%KE_dia,isd,ied,jsd,jed,nz)


  ! gravity wave CFLs
  cs%id_cg1 = register_diag_field('ocean_model', 'cg1', diag%axesT1, time, &
      'First baroclinic gravity wave speed', 'meter second-1')
  cs%id_Rd1 = register_diag_field('ocean_model', 'Rd1', diag%axesT1, time, &
      'First baroclinic deformation radius', 'meter')
  cs%id_cfl_cg1 = register_diag_field('ocean_model', 'CFL_cg1', diag%axesT1, time, &
      'CFL of first baroclinic gravity wave = dt*cg1*(1/dx+1/dy)', 'nondim')
  cs%id_cfl_cg1_x = register_diag_field('ocean_model', 'CFL_cg1_x', diag%axesT1, time, &
      'i-component of CFL of first baroclinic gravity wave = dt*cg1*/dx', 'nondim')
  cs%id_cfl_cg1_y = register_diag_field('ocean_model', 'CFL_cg1_y', diag%axesT1, time, &
      'j-component of CFL of first baroclinic gravity wave = dt*cg1*/dy', 'nondim')
  cs%id_cg_ebt = register_diag_field('ocean_model', 'cg_ebt', diag%axesT1, time, &
      'Equivalent barotropic gravity wave speed', 'meter second-1')
  cs%id_Rd_ebt = register_diag_field('ocean_model', 'Rd_ebt', diag%axesT1, time, &
      'Equivalent barotropic deformation radius', 'meter')
  cs%id_p_ebt = register_diag_field('ocean_model', 'p_ebt', diag%axesTL, time, &
      'Equivalent barotropic modal strcuture', 'nondim')

  if ((cs%id_cg1>0) .or. (cs%id_Rd1>0) .or. (cs%id_cfl_cg1>0) .or. &
      (cs%id_cfl_cg1_x>0) .or. (cs%id_cfl_cg1_y>0) .or. &
      (cs%id_cg_ebt>0) .or. (cs%id_Rd_ebt>0) .or. (cs%id_p_ebt>0)) then
    call wave_speed_init(cs%wave_speed_CSp)
    call safe_alloc_ptr(cs%cg1,isd,ied,jsd,jed)
    if (cs%id_Rd1>0)       call safe_alloc_ptr(cs%Rd1,isd,ied,jsd,jed)
    if (cs%id_Rd_ebt>0)    call safe_alloc_ptr(cs%Rd1,isd,ied,jsd,jed)
    if (cs%id_cfl_cg1>0)   call safe_alloc_ptr(cs%cfl_cg1,isd,ied,jsd,jed)
    if (cs%id_cfl_cg1_x>0) call safe_alloc_ptr(cs%cfl_cg1_x,isd,ied,jsd,jed)
    if (cs%id_cfl_cg1_y>0) call safe_alloc_ptr(cs%cfl_cg1_y,isd,ied,jsd,jed)
    if (cs%id_p_ebt>0) call safe_alloc_ptr(cs%p_ebt,isd,ied,jsd,jed,nz)
  endif

  cs%id_mass_wt = register_diag_field('ocean_model', 'mass_wt', diag%axesT1, time,  &
      'The column mass for calculating mass-weighted average properties', 'kg m-2')

  cs%id_temp_int = register_diag_field('ocean_model', 'temp_int', diag%axesT1, time,                &
      'Density weighted column integrated potential temperature', 'degC kg m-2',                    &
      cmor_field_name='opottempmint',                                                               &
      cmor_long_name='integral_wrt_depth_of_product_of_sea_water_density_and_potential_temperature',&
      cmor_units='degC kg m-2', cmor_standard_name='Depth integrated density times potential temperature')

  cs%id_salt_int = register_diag_field('ocean_model', 'salt_int', diag%axesT1, time,   &
      'Density weighted column integrated salinity', 'ppt kg m-2',                     &
      cmor_field_name='somint',                                                        &
      cmor_long_name='integral_wrt_depth_of_product_of_sea_water_density_and_salinity',&
      cmor_units='ppt kg m-2', cmor_standard_name='Depth integrated density times salinity')

  cs%id_col_mass = register_diag_field('ocean_model', 'col_mass', diag%axesT1, time, &
      'The column integrated in situ density', 'kg m-2')

  cs%id_col_ht = register_diag_field('ocean_model', 'col_height', diag%axesT1, time, &
      'The height of the water column', 'm')
  cs%id_pbo = register_diag_field('ocean_model', 'pbo', diag%axesT1, time, &
      long_name='Sea Water Pressure at Sea Floor', standard_name='sea_water_pressure_at_sea_floor', &
      units='Pa')

  call set_dependent_diagnostics(mis, adp, cdp, g, cs)

end subroutine mom_diagnostics_init

!> This subroutine sets up diagnostics upon which other diagnostics depend.
subroutine set_dependent_diagnostics(MIS, ADp, CDp, G, CS)
  type(ocean_internal_state), intent(in)    :: MIS !< For "MOM Internal State" a set of pointers to
                                                   !! the fields and accelerations making up ocean
                                                   !! internal physical state.
  type(accel_diag_ptrs),      intent(inout) :: ADp !< Structure pointing to accelerations in
                                                   !! momentum equation.
  type(cont_diag_ptrs),       intent(inout) :: CDp !< Structure pointing to terms in continuity
                                                   !! equation.
  type(ocean_grid_type),      intent(in)    :: G   !< The ocean's grid structure.
  type(diagnostics_cs),       pointer       :: CS  !< Pointer to the control structure for this
                                                   !! module.

! This subroutine sets up diagnostics upon which other diagnostics depend.
! Arguments:
!  (in)      MIS - For "MOM Internal State" a set of pointers to the fields and
!                  accelerations making up ocean internal physical state.
!  (inout)   ADp - structure pointing to accelerations in momentum equation
!  (inout)   CDp - structure pointing to terms in continuity equation
!  (in)      G   - ocean grid structure
!  (in)      CS -  pointer to the control structure for this module

  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nz
  isd  = g%isd  ; ied  = g%ied  ; jsd  = g%jsd  ; jed  = g%jed ; nz = g%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  if (ASSOCIATED(cs%dKE_dt) .or. ASSOCIATED(cs%PE_to_KE) .or. &
      ASSOCIATED(cs%KE_CorAdv) .or. ASSOCIATED(cs%KE_adv) .or. &
      ASSOCIATED(cs%KE_visc) .or. ASSOCIATED(cs%KE_horvisc) .or. &
      ASSOCIATED(cs%KE_dia)) &
    call safe_alloc_ptr(cs%KE,isd,ied,jsd,jed,nz)

  if (ASSOCIATED(cs%dKE_dt)) then
    if (.not.ASSOCIATED(cs%du_dt)) then
      call safe_alloc_ptr(cs%du_dt,isdb,iedb,jsd,jed,nz)
      call register_time_deriv(mis%u, cs%du_dt, cs)
    endif
    if (.not.ASSOCIATED(cs%dv_dt)) then
      call safe_alloc_ptr(cs%dv_dt,isd,ied,jsdb,jedb,nz)
      call register_time_deriv(mis%v, cs%dv_dt, cs)
    endif
    if (.not.ASSOCIATED(cs%dh_dt)) then
      call safe_alloc_ptr(cs%dh_dt,isd,ied,jsd,jed,nz)
      call register_time_deriv(mis%h, cs%dh_dt, cs)
    endif
  endif

  if (ASSOCIATED(cs%KE_adv)) then
    call safe_alloc_ptr(adp%gradKEu,isdb,iedb,jsd,jed,nz)
    call safe_alloc_ptr(adp%gradKEv,isd,ied,jsdb,jedb,nz)
  endif

  if (ASSOCIATED(cs%KE_visc)) then
    call safe_alloc_ptr(adp%du_dt_visc,isdb,iedb,jsd,jed,nz)
    call safe_alloc_ptr(adp%dv_dt_visc,isd,ied,jsdb,jedb,nz)
  endif

  if (ASSOCIATED(cs%KE_dia)) then
    call safe_alloc_ptr(adp%du_dt_dia,isdb,iedb,jsd,jed,nz)
    call safe_alloc_ptr(adp%dv_dt_dia,isd,ied,jsdb,jedb,nz)
  endif

  if (ASSOCIATED(cs%uhGM_Rlay)) call safe_alloc_ptr(cdp%uhGM,isdb,iedb,jsd,jed,nz)
  if (ASSOCIATED(cs%vhGM_Rlay)) call safe_alloc_ptr(cdp%vhGM,isd,ied,jsdb,jedb,nz)

end subroutine set_dependent_diagnostics


subroutine mom_diagnostics_end(CS, ADp)
  type(diagnostics_cs),   pointer       :: CS
  type(accel_diag_ptrs),  intent(inout) :: ADp
  integer :: m

  if (ASSOCIATED(cs%e))          deallocate(cs%e)
  if (ASSOCIATED(cs%e_D))        deallocate(cs%e_D)
  if (ASSOCIATED(cs%KE))         deallocate(cs%KE)
  if (ASSOCIATED(cs%dKE_dt))     deallocate(cs%dKE_dt)
  if (ASSOCIATED(cs%PE_to_KE))   deallocate(cs%PE_to_KE)
  if (ASSOCIATED(cs%KE_Coradv))  deallocate(cs%KE_Coradv)
  if (ASSOCIATED(cs%KE_adv))     deallocate(cs%KE_adv)
  if (ASSOCIATED(cs%KE_visc))    deallocate(cs%KE_visc)
  if (ASSOCIATED(cs%KE_horvisc)) deallocate(cs%KE_horvisc)
  if (ASSOCIATED(cs%KE_dia))     deallocate(cs%KE_dia)
  if (ASSOCIATED(cs%dv_dt))      deallocate(cs%dv_dt)
  if (ASSOCIATED(cs%dh_dt))      deallocate(cs%dh_dt)
  if (ASSOCIATED(cs%du_dt))      deallocate(cs%du_dt)
  if (ASSOCIATED(cs%h_Rlay))     deallocate(cs%h_Rlay)
  if (ASSOCIATED(cs%uh_Rlay))    deallocate(cs%uh_Rlay)
  if (ASSOCIATED(cs%vh_Rlay))    deallocate(cs%vh_Rlay)
  if (ASSOCIATED(cs%uhGM_Rlay))  deallocate(cs%uhGM_Rlay)
  if (ASSOCIATED(cs%vhGM_Rlay))  deallocate(cs%vhGM_Rlay)
  if (ASSOCIATED(cs%diag_tmp3d)) deallocate(cs%diag_tmp3d)

  if (ASSOCIATED(adp%gradKEu))    deallocate(adp%gradKEu)
  if (ASSOCIATED(adp%gradKEu))    deallocate(adp%gradKEu)
  if (ASSOCIATED(adp%du_dt_visc)) deallocate(adp%du_dt_visc)
  if (ASSOCIATED(adp%dv_dt_visc)) deallocate(adp%dv_dt_visc)
  if (ASSOCIATED(adp%du_dt_dia))  deallocate(adp%du_dt_dia)
  if (ASSOCIATED(adp%dv_dt_dia))  deallocate(adp%dv_dt_dia)
  if (ASSOCIATED(adp%du_other))   deallocate(adp%du_other)
  if (ASSOCIATED(adp%dv_other))   deallocate(adp%dv_other)

  do m=1,cs%num_time_deriv ; deallocate(cs%prev_val(m)%p) ; enddo

  deallocate(cs)

end subroutine mom_diagnostics_end

end module mom_diagnostics