MOM_diagnostics.F90

!> Calculates any requested diagnostic quantities
!! that are not calculated in the various subroutines.
!! Diagnostic quantities are requested by allocating them memory.
module mom_diagnostics
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_coms,              only : reproducing_sum
use mom_density_integrals, only : int_density_dz
use mom_diag_mediator,     only : post_data, get_diag_time_end
use mom_diag_mediator,     only : register_diag_field, register_scalar_field
use mom_diag_mediator,     only : register_static_field, diag_register_area_ids
use mom_diag_mediator,     only : diag_ctrl, time_type, safe_alloc_ptr
use mom_diag_mediator,     only : diag_get_volume_cell_measure_dm_id
use mom_diag_mediator,     only : diag_grid_storage
use mom_diag_mediator,     only : diag_save_grids, diag_restore_grids, diag_copy_storage_to_diag
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, calculate_density_derivs, eos_domain
use mom_eos,               only : gsw_sp_from_sr, gsw_pt_from_ct
use mom_error_handler,     only : mom_error, fatal, warning
use mom_file_parser,       only : get_param, log_version, param_file_type
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
use mom_spatial_means,     only : global_volume_mean, global_area_integral
use mom_tracer_registry,   only : tracer_registry_type, post_tracer_transport_diagnostics
use mom_unit_scaling,      only : unit_scale_type
use mom_variables,         only : thermo_var_ptrs, ocean_internal_state, p3d
use mom_variables,         only : accel_diag_ptrs, cont_diag_ptrs, surface
use mom_verticalgrid,      only : verticalgrid_type, get_thickness_units
use mom_wave_speed,        only : wave_speed, wave_speed_cs, wave_speed_init
use coupler_types_mod,     only : coupler_type_send_data
 
implicit none ; private
 
#include <MOM_memory.h>
 
public calculate_diagnostic_fields, register_time_deriv, write_static_fields
public find_eta
public register_surface_diags, post_surface_dyn_diags, post_surface_thermo_diags
public register_transport_diags, post_transport_diagnostics
public mom_diagnostics_init, mom_diagnostics_end
 
! A note on unit descriptions in comments: MOM6 uses units that can be rescaled for dimensional
! consistency testing. These are noted in comments with units like Z, H, L, and T, along with
! their mks counterparts with notation like "a velocity [Z T-1 ~> m s-1]".  If the units
! vary with the Boussinesq approximation, the Boussinesq variant is given first.
 
!> The control structure for the MOM_diagnostics module
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 [Z ~> m].
 
  type(diag_ctrl), pointer :: diag => null() !< A structure that is used to
                                       !! regulate the timing of diagnostic output.
 
  ! following arrays store diagnostics calculated here and unavailable outside.
 
  ! following fields have nz+1 levels.
  real, pointer, dimension(:,:,:) :: &
    e => null(), &   !< interface height [Z ~> m]
    e_d => null()    !< interface height above bottom [Z ~> m]
 
  ! following fields have nz layers.
  real, pointer, dimension(:,:,:) :: &
    du_dt => null(), & !< net i-acceleration [L T-2 ~> m s-2]
    dv_dt => null(), & !< net j-acceleration [L T-2 ~> m s-2]
    dh_dt => null(), & !< thickness rate of change [H T-1 ~> m s-1 or kg m-2 s-1]
    p_ebt => null()    !< Equivalent barotropic modal structure [nondim]
    ! hf_du_dt => NULL(), hf_dv_dt => NULL() !< du_dt, dv_dt x fract. thickness [L T-2 ~> m s-2].
    ! 3D diagnostics hf_du(dv)_dt are commented because there is no clarity on proper remapping grid option.
    ! The code is retained for degugging purposes in the future.
 
  real, pointer, dimension(:,:,:) :: h_rlay => null() !< Layer thicknesses in potential density
                                              !! coordinates [H ~> m or kg m-2]
  real, pointer, dimension(:,:,:) :: uh_rlay => null() !< Zonal transports in potential density
                                              !! coordinates [H m2 s-1 ~> m3 s-1 or kg s-1]
  real, pointer, dimension(:,:,:) :: vh_rlay => null() !< Meridional transports in potential density
                                              !! coordinates [H m2 s-1 ~> m3 s-1 or kg s-1]
  real, pointer, dimension(:,:,:) :: uhgm_rlay => null() !< Zonal Gent-McWilliams transports in potential density
                                              !! coordinates [H m2 s-1 ~> m3 s-1 or kg s-1]
  real, pointer, dimension(:,:,:) :: vhgm_rlay => null() !< Meridional Gent-McWilliams transports in potential density
                                              !! coordinates [H m2 s-1 ~> m3 s-1 or kg s-1]
 
  ! following fields are 2-D.
  real, pointer, dimension(:,:) :: &
    cg1 => null(),       & !< First baroclinic gravity wave speed [L T-1 ~> m s-1]
    rd1 => null(),       & !< First baroclinic deformation radius [L ~> 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]
 
  ! The following arrays hold diagnostics in the layer-integrated energy budget.
  real, pointer, dimension(:,:,:) :: &
    ke        => null(), &  !< KE per unit mass [L2 T-2 ~> m2 s-2]
    dke_dt    => null(), &  !< time derivative of the layer KE [H L2 T-3 ~> m3 s-3]
    pe_to_ke  => null(), &  !< potential energy to KE term [m3 s-3]
    ke_coradv => null(), &  !< KE source from the combined Coriolis and
                            !! advection terms [H L2 T-3 ~> m3 s-3].
                            !! 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 [H L2 T-3 ~> m3 s-3]
    ke_visc    => null(), & !< KE source from vertical viscosity [H L2 T-3 ~> m3 s-3]
    ke_horvisc => null(), & !< KE source from horizontal viscosity [H L2 T-3 ~> m3 s-3]
    ke_dia     => null()    !< KE source from diapycnal diffusion [H L2 T-3 ~> m3 s-3]
 
  !>@{ Diagnostic IDs
  integer :: id_u = -1,   id_v = -1, id_h = -1
  integer :: id_e              = -1, id_e_d            = -1
  integer :: id_du_dt          = -1, id_dv_dt          = -1
  ! integer :: id_hf_du_dt       = -1, id_hf_dv_dt       = -1
  integer :: id_hf_du_dt_2d    = -1, id_hf_dv_dt_2d    = -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_volcello       = -1
  integer :: id_tpot           = -1, id_sprac          = -1
  integer :: id_tob            = -1, id_sob            = -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
  integer :: id_drho_dt        = -1, id_drho_ds        = -1
  integer :: id_h_pre_sync     = -1
  !>@}
  !> The control structure for calculating wave speed.
  type(wave_speed_cs), pointer :: wave_speed_csp => null()
 
  type(p3d) :: var_ptr(max_fields_)  !< pointers to variables used in the calculation
                                     !! of time derivatives
  type(p3d) :: deriv(max_fields_)    !< Time derivatives of various fields
  type(p3d) :: prev_val(max_fields_) !< Previous values of variables used in the calculation
                                     !! of time derivatives
  !< previous values of variables used in calculation of time derivatives
  integer   :: nlay(max_fields_) !< The number of layers in each diagnostics
  integer   :: num_time_deriv = 0 !< The number of time derivative diagnostics
 
  type(group_pass_type) :: pass_ke_uv !< A handle used for group halo passes
 
end type diagnostics_cs
 
 
!> A structure with diagnostic IDs of the surface and integrated variables
type, public :: surface_diag_ids ; private
  !>@{ Diagnostic IDs for 2-d surface and bottom flux and state fields
  !Diagnostic IDs for 2-d surface and bottom fields
  integer :: id_zos  = -1, id_zossq  = -1
  integer :: id_volo = -1, id_speed  = -1
  integer :: id_ssh  = -1, id_ssh_ga = -1
  integer :: id_sst  = -1, id_sst_sq = -1, id_sstcon = -1
  integer :: id_sss  = -1, id_sss_sq = -1, id_sssabs = -1
  integer :: id_ssu  = -1, id_ssv    = -1
 
  ! Diagnostic IDs for  heat and salt flux fields
  integer :: id_fraz         = -1
  integer :: id_salt_deficit = -1
  integer :: id_heat_pme     = -1
  integer :: id_intern_heat  = -1
  !>@}
end type surface_diag_ids
 
 
!> A structure with diagnostic IDs of mass transport related diagnostics
type, public :: transport_diag_ids ; private
  !>@{  Diagnostics for tracer horizontal transport
  integer :: id_uhtr = -1, id_umo = -1, id_umo_2d = -1
  integer :: id_vhtr = -1, id_vmo = -1, id_vmo_2d = -1
  integer :: id_dynamics_h = -1, id_dynamics_h_tendency = -1
  !>@}
end type transport_diag_ids
 
 
contains
!> Diagnostics not more naturally calculated elsewhere are computed here.
subroutine calculate_diagnostic_fields(u, v, h, uh, vh, tv, ADp, CDp, p_surf, &
                                       dt, diag_pre_sync, G, GV, US, 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.
  type(unit_scale_type),   intent(in)    :: us   !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: u    !< The zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: v    !< The meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: uh   !< Transport through zonal faces = u*h*dy,
                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: vh   !< Transport through meridional faces = v*h*dx,
                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  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.
  real, dimension(:,:),    pointer       :: p_surf !< A pointer to the surface pressure [R L2 T-2 ~> Pa].
                                                 !! If p_surf is not associated, it is the same
                                                 !! as setting the surface pressure to 0.
  real,                    intent(in)    :: dt   !< The time difference since the last
                                                 !! call to this subroutine [T ~> s].
  type(diag_grid_storage), intent(in)    :: diag_pre_sync !< Target grids from previous timestep
  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 [H ~> m or kg m-2].
 
  ! Local variables
  integer, dimension(2) :: eosdom ! The i-computational domain for the equation of state
  integer i, j, k, is, ie, js, je, isq, ieq, jsq, jeq, nz, nkmb
 
  real :: rcv(szi_(g),szj_(g),szk_(g))  ! Coordinate variable potential density [R ~> kg m-3].
  real :: work_3d(szi_(g),szj_(g),szk_(g)) ! A 3-d temporary work array.
  real :: work_2d(szi_(g),szj_(g))         ! A 2-d temporary work array.
  real :: rho_in_situ(szi_(g))             ! In situ density [R ~> kg m-3]
 
  real, allocatable, dimension(:,:) :: &
    hf_du_dt_2d, hf_dv_dt_2d ! z integeral of hf_du_dt, hf_dv_dt [L T-2 ~> m s-2].
 
  ! tmp array for surface properties
  real :: surface_field(szi_(g),szj_(g))
  real :: pressure_1d(szi_(g)) ! Temporary array for pressure when calling EOS [R L2 T-2 ~> Pa]
  real :: wt, wt_p
 
  real :: f2_h     ! Squared Coriolis parameter at to h-points [T-2 ~> s-2]
  real :: mag_beta ! Magnitude of the gradient of f [T-1 L-1 ~> s-1 m-1]
  real :: absurdly_small_freq2 ! Frequency squared used to avoid division by 0 [T-2 ~> s-2]
 
  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
 
  ! This value is roughly (pi / (the age of the universe) )^2.
  absurdly_small_freq2 = 1e-34*us%T_to_s**2
 
  if (loc(cs)==0) call mom_error(fatal, &
         "calculate_diagnostic_fields: Module must be initialized before used.")
 
  call calculate_derivs(dt, g, cs)
 
  if (dt > 0.0) then
    call diag_save_grids(cs%diag)
    call diag_copy_storage_to_diag(cs%diag, diag_pre_sync)
 
    if (cs%id_h_pre_sync > 0) &
        call post_data(cs%id_h_pre_sync, diag_pre_sync%h_state, cs%diag, alt_h = diag_pre_sync%h_state)
 
    if (cs%id_du_dt>0) call post_data(cs%id_du_dt, cs%du_dt, cs%diag, alt_h = diag_pre_sync%h_state)
 
    if (cs%id_dv_dt>0) call post_data(cs%id_dv_dt, cs%dv_dt, cs%diag, alt_h = diag_pre_sync%h_state)
 
    if (cs%id_dh_dt>0) call post_data(cs%id_dh_dt, cs%dh_dt, cs%diag, alt_h = diag_pre_sync%h_state)
 
    !! Diagnostics for terms multiplied by fractional thicknesses
 
    ! 3D diagnostics hf_du(dv)_dt are commented because there is no clarity on proper remapping grid option.
    ! The code is retained for degugging purposes in the future.
    !if (CS%id_hf_du_dt > 0) then
    !  do k=1,nz ; do j=js,je ; do I=Isq,Ieq
    !    CS%hf_du_dt(I,j,k) = CS%du_dt(I,j,k) * ADp%diag_hfrac_u(I,j,k)
    !  enddo ; enddo ; enddo
    !  call post_data(CS%id_hf_du_dt, CS%hf_du_dt, CS%diag, alt_h = diag_pre_sync%h_state)
    !endif
 
    !if (CS%id_hf_dv_dt > 0) then
    !  do k=1,nz ; do J=Jsq,Jeq ; do i=is,ie
    !    CS%hf_dv_dt(i,J,k) = CS%dv_dt(i,J,k) * ADp%diag_hfrac_v(i,J,k)
    !  enddo ; enddo ; enddo
    !  call post_data(CS%id_hf_dv_dt, CS%hf_dv_dt, CS%diag, alt_h = diag_pre_sync%h_state)
    !endif
 
    if (cs%id_hf_du_dt_2d > 0) then
      allocate(hf_du_dt_2d(g%IsdB:g%IedB,g%jsd:g%jed))
      hf_du_dt_2d(:,:) = 0.0
      do k=1,nz ; do j=js,je ; do i=isq,ieq
        hf_du_dt_2d(i,j) = hf_du_dt_2d(i,j) + cs%du_dt(i,j,k) * adp%diag_hfrac_u(i,j,k)
      enddo ; enddo ; enddo
      call post_data(cs%id_hf_du_dt_2d, hf_du_dt_2d, cs%diag)
      deallocate(hf_du_dt_2d)
    endif
 
    if (cs%id_hf_dv_dt_2d > 0) then
      allocate(hf_dv_dt_2d(g%isd:g%ied,g%JsdB:g%JedB))
      hf_dv_dt_2d(:,:) = 0.0
      do k=1,nz ; do j=jsq,jeq ; do i=is,ie
        hf_dv_dt_2d(i,j) = hf_dv_dt_2d(i,j) + cs%dv_dt(i,j,k) * adp%diag_hfrac_v(i,j,k)
      enddo ; enddo ; enddo
      call post_data(cs%id_hf_dv_dt_2d, hf_dv_dt_2d, cs%diag)
      deallocate(hf_dv_dt_2d)
    endif
 
    call diag_restore_grids(cs%diag)
 
    call calculate_energy_diagnostics(u, v, h, uh, vh, adp, cdp, g, gv, us, cs)
  endif
 
  ! 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 .and. nz > 1) nkmb = nz
 
  if (cs%id_u > 0) call post_data(cs%id_u, u, cs%diag)
 
  if (cs%id_v > 0) call post_data(cs%id_v, v, cs%diag)
 
  if (cs%id_h > 0) call post_data(cs%id_h, h, cs%diag)
 
  if (associated(cs%e)) then
    call find_eta(h, tv, g, gv, us, 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, g, gv, us, 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
      work_3d(i,j,k) = gv%H_to_kg_m2*h(i,j,k)
    enddo ; enddo ; enddo
    call post_data(cs%id_masscello, work_3d, cs%diag)
  endif
 
  ! mass of liquid ocean (for Bouss, use Rho0). The reproducing sum requires the use of MKS units.
  if (cs%id_masso > 0) then
    work_2d(:,:) = 0.0
    do k=1,nz ; do j=js,je ; do i=is,ie
      work_2d(i,j) = work_2d(i,j) + (gv%H_to_kg_m2*h(i,j,k)) * us%L_to_m**2*g%areaT(i,j)
    enddo ; enddo ; enddo
    masso = reproducing_sum(work_2d)
    call post_data(cs%id_masso, masso, cs%diag)
  endif
 
  ! diagnose thickness/volumes of grid cells [m]
  if (cs%id_thkcello>0 .or. cs%id_volcello>0) then
    if (gv%Boussinesq) then ! thkcello = h for Boussinesq
      if (cs%id_thkcello > 0) then ; if (gv%H_to_Z == 1.0) then
        call post_data(cs%id_thkcello, h, cs%diag)
      else
        do k=1,nz; do j=js,je ; do i=is,ie
          work_3d(i,j,k) = gv%H_to_Z*h(i,j,k)
        enddo ; enddo ; enddo
        call post_data(cs%id_thkcello, work_3d, cs%diag)
      endif ; endif
      if (cs%id_volcello > 0) then ! volcello = h*area for Boussinesq
        do k=1,nz; do j=js,je ; do i=is,ie
          work_3d(i,j,k) = ( gv%H_to_Z*h(i,j,k) ) * us%Z_to_m*us%L_to_m**2*g%areaT(i,j)
        enddo ; enddo ; enddo
        call post_data(cs%id_volcello, work_3d, cs%diag)
      endif
    else ! thkcello = dp/(rho*g) for non-Boussinesq
      eosdom(:) = eos_domain(g%HI)
      do j=js,je
        if (associated(p_surf)) then ! Pressure loading at top of surface layer [R L2 T-2 ~> Pa]
          do i=is,ie
            pressure_1d(i) = p_surf(i,j)
          enddo
        else
          do i=is,ie
            pressure_1d(i) = 0.0
          enddo
        endif
        do k=1,nz ! Integrate vertically downward for pressure
          do i=is,ie ! Pressure for EOS at the layer center [R L2 T-2 ~> Pa]
            pressure_1d(i) = pressure_1d(i) + 0.5*(gv%H_to_RZ*gv%g_Earth)*h(i,j,k)
          enddo
          ! Store in-situ density [R ~> kg m-3] in work_3d
          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, rho_in_situ, &
                                 tv%eqn_of_state, eosdom)
          do i=is,ie ! Cell thickness = dz = dp/(g*rho) (meter); store in work_3d
            work_3d(i,j,k) = (gv%H_to_RZ*h(i,j,k)) / rho_in_situ(i)
          enddo
          do i=is,ie ! Pressure for EOS at the bottom interface [R L2 T-2 ~> Pa]
            pressure_1d(i) = pressure_1d(i) + 0.5*(gv%H_to_RZ*gv%g_Earth)*h(i,j,k)
          enddo
        enddo ! k
      enddo ! j
      if (cs%id_thkcello > 0) call post_data(cs%id_thkcello, work_3d, cs%diag)
      if (cs%id_volcello > 0) then
        do k=1,nz; do j=js,je ; do i=is,ie ! volcello = dp/(rho*g)*area for non-Boussinesq
          work_3d(i,j,k) = us%Z_to_m*us%L_to_m**2*g%areaT(i,j) * work_3d(i,j,k)
        enddo ; enddo ; enddo
        call post_data(cs%id_volcello, work_3d, cs%diag)
      endif
    endif
  endif
 
  ! Calculate additional, potentially derived temperature diagnostics
  if (tv%T_is_conT) then
    ! Internal T&S variables are conservative temperature & absolute salinity,
    ! so they need to converted to potential temperature and practical salinity
    ! for some diagnostics using TEOS-10 function calls.
    if ((cs%id_Tpot > 0) .or. (cs%id_tob > 0)) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        work_3d(i,j,k) = gsw_pt_from_ct(tv%S(i,j,k),tv%T(i,j,k))
      enddo ; enddo ; enddo
      if (cs%id_Tpot > 0) call post_data(cs%id_Tpot, work_3d, cs%diag)
      if (cs%id_tob > 0) call post_data(cs%id_tob, work_3d(:,:,nz), cs%diag, mask=g%mask2dT)
    endif
  else
    ! Internal T&S variables are potential temperature & practical salinity
    if (cs%id_tob > 0) call post_data(cs%id_tob, tv%T(:,:,nz), cs%diag, mask=g%mask2dT)
  endif
 
  ! Calculate additional, potentially derived salinity diagnostics
  if (tv%S_is_absS) then
    ! Internal T&S variables are conservative temperature & absolute salinity,
    ! so they need to converted to potential temperature and practical salinity
    ! for some diagnostics using TEOS-10 function calls.
    if ((cs%id_Sprac > 0) .or. (cs%id_sob > 0)) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        work_3d(i,j,k) = gsw_sp_from_sr(tv%S(i,j,k))
      enddo ; enddo ; enddo
      if (cs%id_Sprac > 0) call post_data(cs%id_Sprac, work_3d, cs%diag)
      if (cs%id_sob > 0) call post_data(cs%id_sob, work_3d(:,:,nz), cs%diag, mask=g%mask2dT)
    endif
  else
    ! Internal T&S variables are potential temperature & practical salinity
    if (cs%id_sob > 0) call post_data(cs%id_sob, tv%S(:,:,nz), cs%diag, mask=g%mask2dT)
  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(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(cs%id_salt_layer_ave, salt_layer_ave, cs%diag)
  endif
 
  call calculate_vertical_integrals(h, tv, p_surf, g, gv, us, 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
      eosdom(:) = eos_domain(g%HI, halo=1)
      pressure_1d(:) = tv%P_Ref
      !$OMP parallel do default(shared)
      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), tv%eqn_of_state, &
                               eosdom)
      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_uhGM_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
    eosdom(:) = eos_domain(g%HI)
    if (cs%id_rhopot0 > 0) then
      pressure_1d(:) = 0.
      !$OMP parallel do default(shared)
      do k=1,nz ; do j=js,je
        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, rcv(:,j,k), &
                                tv%eqn_of_state, eosdom)
      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.0e7*us%kg_m3_to_R*us%m_s_to_L_T**2 ! 2000 dbars
      !$OMP parallel do default(shared)
      do k=1,nz ; do j=js,je
        call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, rcv(:,j,k), &
                                tv%eqn_of_state, eosdom)
      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(shared) private(pressure_1d)
      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_RZ*gv%g_Earth) ! Pressure in middle of layer k
          call calculate_density(tv%T(:,j,k), tv%S(:,j,k), pressure_1d, rcv(:,j,k), &
                                tv%eqn_of_state, eosdom)
          pressure_1d(:) =  pressure_1d(:) + 0.5 * h(:,j,k) * (gv%H_to_RZ*gv%g_Earth) ! Pressure at bottom of layer k
        enddo
      enddo
      if (cs%id_rhoinsitu > 0) call post_data(cs%id_rhoinsitu, rcv, cs%diag)
    endif
 
    if (cs%id_drho_dT > 0 .or. cs%id_drho_dS > 0) then
      !$OMP parallel do default(shared) private(pressure_1d)
      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
          ! To avoid storing more arrays, put drho_dT into Rcv, and drho_dS into work3d
          call calculate_density_derivs(tv%T(:,j,k),tv%S(:,j,k),pressure_1d, &
                                 rcv(:,j,k),work_3d(:,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_drho_dT > 0) call post_data(cs%id_drho_dT, rcv, cs%diag)
      if (cs%id_drho_dS > 0) call post_data(cs%id_drho_dS, work_3d, 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, us, 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(shared) 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, us, 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, us, 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(shared) 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
 
end subroutine calculate_diagnostic_fields
 
!> This subroutine finds the 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, dimension(:), &
            intent(in)    :: Rlist !< The list of target densities [R ~> kg m-3]
  real,     intent(in)    :: R_in !< The density being inserted into Rlist [R ~> kg m-3]
  integer,  intent(inout) :: k    !< The value of k such that Rlist(k) <= R_in < Rlist(k+1)
                                  !! The input value is a first guess
  integer,  intent(in)    :: nz   !< The number of layers in Rlist
  real,     intent(out)   :: wt   !< The weight of layer k for interpolation, nondim
  real,     intent(out)   :: wt_p !< The weight of layer k+1 for interpolation, nondim
 
  ! 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
    enddo
  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
    enddo
  endif
 
  if ((k_lower == 1) .and. (r_in <= rlist(k_lower))) then
    k = 1 ; wt = 1.0 ; wt_p = 0.0
  elseif ((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
    enddo
 
!   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
 
!> This 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, p_surf, G, GV, US, 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.
  type(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2].
  type(thermo_var_ptrs),   intent(in)    :: tv   !< A structure pointing to various
                                                 !! thermodynamic variables.
  real, dimension(:,:),    pointer       :: p_surf !< A pointer to the surface pressure [R L2 T-2 ~> Pa].
                                                 !! If p_surf is not associated, it is the same
                                                 !! as setting the surface pressure to 0.
  type(diagnostics_cs),    intent(inout) :: 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 [Z ~> m].
    z_bot, &  ! Height of the bottom of a layer (for id_mass) or the
              ! (positive) depth of the ocean (for id_col_ht) [Z ~> m].
    mass, &   ! integrated mass of the water column [R Z ~> kg m-2].  For
              ! non-Boussinesq models this is rho*dz. For Boussinesq
              ! models, this is either the integral of in-situ density
              ! (rho*dz for col_mass) or reference density (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 [R L2 T-2 ~> Pa].
    dpress, & ! Change in hydrostatic pressure across a layer [R L2 T-2 ~> Pa].
    tr_int    ! vertical integral of a tracer times density,
              ! (Rho_0 in a Boussinesq model) [TR kg m-2].
  real    :: IG_Earth  ! Inverse of gravitational acceleration [T2 Z L-2 ~> 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_RZ*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_RZ*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_RZ*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, g, gv, us, 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
 
  ! NOTE: int_density_z expects z_top and z_btm values from [ij]sq to [ij]eq+1
  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=g%jscB,g%jecB+1 ; do i=g%iscB,g%iecB+1
          z_bot(i,j) = 0.0
        enddo ; enddo
        do k=1,nz
          do j=g%jscB,g%jecB+1 ; do i=g%iscB,g%iecB+1
            z_top(i,j) = z_bot(i,j)
            z_bot(i,j) = z_top(i,j) - gv%H_to_Z*h(i,j,k)
          enddo ; enddo
          call int_density_dz(tv%T(:,:,k), tv%S(:,:,k), z_top, z_bot, 0.0, gv%Rho0, gv%g_Earth, &
                              g%HI, tv%eqn_of_state, us, 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_Z*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_RZ*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) + p_surf
      ! where p_surf is the sea water pressure at sea water surface.
      do j=js,je ; do i=is,ie
        btm_pres(i,j) = gv%g_Earth * mass(i,j)
        if (associated(p_surf)) then
          btm_pres(i,j) = btm_pres(i,j) + 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, GV, US, 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(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: u    !< The zonal velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: v    !< The meridional velocity [L T-1 ~> m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: h    !< Layer thicknesses [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: uh   !< Transport through zonal faces=u*h*dy,
                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: vh   !< Transport through merid faces=v*h*dx,
                                                 !! [H L2 T-1 ~> m3 s-1 or kg s-1].
  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(unit_scale_type),   intent(in)    :: US   !< A dimensional unit scaling type
  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.
 
  ! Local variables
  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
    ! NOTE: All terms in KE_adv are multipled by -1, which can easily produce
    ! negative zeros and may signal a reproducibility issue over land.
    ! We resolve this by re-initializing and only evaluating over water points.
    ke_u(:,:) = 0. ; ke_v(:,:) = 0.
    do k=1,nz
      do j=js,je ; do i=isq,ieq
        if (g%mask2dCu(i,j) /= 0.) &
          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
        if (g%mask2dCv(i,j) /= 0.) &
          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) &
            * (us%T_to_s * (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(lb, f_ptr, deriv_ptr, CS)
  integer, intent(in), dimension(3) :: lb     !< Lower index bound of f_ptr
  real, dimension(lb(1):,lb(2):,:), target :: f_ptr
                                              !< Time derivative operand
  real, dimension(lb(1):,lb(2):,:), target :: deriv_ptr
                                              !< Time derivative of f_ptr
  type(diagnostics_cs),  pointer :: cs        !< Control structure returned by previous call to
                                              !! diagnostics_init.
 
  ! This subroutine registers fields to calculate a diagnostic time derivative.
  ! NOTE: Lower bound is required for grid indexing in calculate_derivs().
  !       We assume that the vertical axis is 1-indexed.
 
  integer :: m      !< New index of deriv_ptr in CS%deriv
  integer :: ub(3)  !< Upper index bound of f_ptr, based on shape.
 
  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
 
  ub(:) = lb(:) + shape(f_ptr) - 1
 
  cs%nlay(m) = size(f_ptr, 3)
  cs%deriv(m)%p => deriv_ptr
  allocate(cs%prev_val(m)%p(lb(1):ub(1), lb(2):ub(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 [T ~> 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.
  real :: Idt  ! The inverse timestep [T-1 ~> s-1]
  integer :: i, j, k, m
 
  if (dt > 0.0) then ; idt = 1.0/dt
  else ; return ; endif
 
  ! Because the field is unknown, its grid index bounds are also unknown.
  ! Additionally, two of the fields (dudt, dvdt) require calculation of spatial
  ! derivatives when computing d(KE)/dt.  This raises issues in non-symmetric
  ! mode, where the symmetric boundaries (west, south) may not be updated.
 
  ! For this reason, we explicitly loop from isc-1:iec and jsc-1:jec, in order
  ! to force boundary value updates, even though it may not be strictly valid
  ! for all fields.  Note this assumes a halo, and that it has been updated.
 
  do m=1,cs%num_time_deriv
    do k=1,cs%nlay(m) ; do j=g%jsc-1,g%jec ; do i=g%isc-1,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 routine posts diagnostics of various dynamic ocean surface quantities,
!! including velocities, speed and sea surface height, at the time the ocean
!! state is reported back to the caller
subroutine post_surface_dyn_diags(IDs, G, diag, sfc_state, ssh)
  type(surface_diag_ids),   intent(in) :: ids !< A structure with the diagnostic IDs.
  type(ocean_grid_type),    intent(in) :: g   !< ocean grid structure
  type(diag_ctrl),          intent(in) :: diag !< regulates diagnostic output
  type(surface),            intent(in) :: sfc_state !< structure describing the ocean surface state
  real, dimension(SZI_(G),SZJ_(G)), &
                            intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m]
 
  ! Local variables
  real, dimension(SZI_(G),SZJ_(G)) :: speed  ! The surface speed [L T-1 ~> m s-1]
  integer :: i, j, is, ie, js, je
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
 
  if (ids%id_ssh > 0) &
    call post_data(ids%id_ssh, ssh, diag, mask=g%mask2dT)
 
  if (ids%id_ssu > 0) &
    call post_data(ids%id_ssu, sfc_state%u, diag, mask=g%mask2dCu)
 
  if (ids%id_ssv > 0) &
    call post_data(ids%id_ssv, sfc_state%v, diag, mask=g%mask2dCv)
 
  if (ids%id_speed > 0) then
    do j=js,je ; do i=is,ie
      speed(i,j) = sqrt(0.5*(sfc_state%u(i-1,j)**2 + sfc_state%u(i,j)**2) + &
                        0.5*(sfc_state%v(i,j-1)**2 + sfc_state%v(i,j)**2))
    enddo ; enddo
    call post_data(ids%id_speed, speed, diag, mask=g%mask2dT)
  endif
 
end subroutine post_surface_dyn_diags
 
 
!> This routine posts diagnostics of various ocean surface and integrated
!! quantities at the time the ocean state is reported back to the caller
subroutine post_surface_thermo_diags(IDs, G, GV, US, diag, dt_int, sfc_state, tv, &
                                    ssh, ssh_ibc)
  type(surface_diag_ids),   intent(in) :: ids !< A structure with the diagnostic IDs.
  type(ocean_grid_type),    intent(in) :: g   !< ocean grid structure
  type(verticalgrid_type),  intent(in) :: gv  !< ocean vertical grid structure
  type(unit_scale_type),    intent(in) :: us  !< A dimensional unit scaling type
  type(diag_ctrl),          intent(in) :: diag  !< regulates diagnostic output
  real,                     intent(in) :: dt_int !< total time step associated with these diagnostics [T ~> s].
  type(surface),            intent(in) :: sfc_state !< structure describing the ocean surface state
  type(thermo_var_ptrs),    intent(in) :: tv  !< A structure pointing to various thermodynamic variables
  real, dimension(SZI_(G),SZJ_(G)), &
                            intent(in) :: ssh !< Time mean surface height without corrections for ice displacement [m]
  real, dimension(SZI_(G),SZJ_(G)), intent(in) :: ssh_ibc !< Time mean surface height with corrections
                                                  !! for ice displacement and the inverse barometer [m]
 
  real, dimension(SZI_(G),SZJ_(G)) :: work_2d  ! A 2-d work array
  real, dimension(SZI_(G),SZJ_(G)) :: &
    zos  ! dynamic sea lev (zero area mean) from inverse-barometer adjusted ssh [m]
  real :: i_time_int    ! The inverse of the time interval [T-1 ~> s-1].
  real :: zos_area_mean, volo, ssh_ga
  integer :: i, j, is, ie, js, je
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
 
  ! area mean SSH
  if (ids%id_ssh_ga > 0) then
    ssh_ga = global_area_mean(ssh, g)
    call post_data(ids%id_ssh_ga, ssh_ga, diag)
  endif
 
  ! post the dynamic sea level, zos, and zossq.
  ! zos is ave_ssh with sea ice inverse barometer removed,
  ! and with zero global area mean.
  if (ids%id_zos > 0 .or. ids%id_zossq > 0) then
    zos(:,:) = 0.0
    do j=js,je ; do i=is,ie
      zos(i,j) = ssh_ibc(i,j)
    enddo ; enddo
    zos_area_mean = global_area_mean(zos, g)
    do j=js,je ; do i=is,ie
      zos(i,j) = zos(i,j) - g%mask2dT(i,j)*zos_area_mean
    enddo ; enddo
    if (ids%id_zos > 0) call post_data(ids%id_zos, zos, diag, mask=g%mask2dT)
    if (ids%id_zossq > 0) then
      do j=js,je ; do i=is,ie
        work_2d(i,j) = zos(i,j)*zos(i,j)
      enddo ; enddo
      call post_data(ids%id_zossq, work_2d, diag, mask=g%mask2dT)
    endif
  endif
 
  ! post total volume of the liquid ocean
  if (ids%id_volo > 0) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = g%mask2dT(i,j)*(ssh(i,j) + us%Z_to_m*g%bathyT(i,j))
    enddo ; enddo
    volo = global_area_integral(work_2d, g)
    call post_data(ids%id_volo, volo, diag)
  endif
 
  ! Use Adcroft's rule of reciprocals; it does the right thing here.
  i_time_int = 0.0 ; if (dt_int > 0.0) i_time_int = 1.0 / dt_int
 
  ! post time-averaged rate of frazil formation
  if (associated(tv%frazil) .and. (ids%id_fraz > 0)) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = tv%frazil(i,j) * i_time_int
    enddo ; enddo
    call post_data(ids%id_fraz, work_2d, diag, mask=g%mask2dT)
  endif
 
  ! post time-averaged salt deficit
  if (associated(tv%salt_deficit) .and. (ids%id_salt_deficit > 0)) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = tv%salt_deficit(i,j) * i_time_int
    enddo ; enddo
    call post_data(ids%id_salt_deficit, work_2d, diag, mask=g%mask2dT)
  endif
 
  ! post temperature of P-E+R
  if (associated(tv%TempxPmE) .and. (ids%id_Heat_PmE > 0)) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = tv%TempxPmE(i,j) * (tv%C_p * i_time_int)
    enddo ; enddo
    call post_data(ids%id_Heat_PmE, work_2d, diag, mask=g%mask2dT)
  endif
 
  ! post geothermal heating or internal heat source/sinks
  if (associated(tv%internal_heat) .and. (ids%id_intern_heat > 0)) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = tv%internal_heat(i,j) * (tv%C_p * i_time_int)
    enddo ; enddo
    call post_data(ids%id_intern_heat, work_2d, diag, mask=g%mask2dT)
  endif
 
  if (tv%T_is_conT) then
    ! Internal T&S variables are conservative temperature & absolute salinity
    if (ids%id_sstcon > 0) call post_data(ids%id_sstcon, sfc_state%SST, diag, mask=g%mask2dT)
    ! Use TEOS-10 function calls convert T&S diagnostics from conservative temp
    ! to potential temperature.
    do j=js,je ; do i=is,ie
      work_2d(i,j) = gsw_pt_from_ct(sfc_state%SSS(i,j), sfc_state%SST(i,j))
    enddo ; enddo
    if (ids%id_sst > 0) call post_data(ids%id_sst, work_2d, diag, mask=g%mask2dT)
  else
    ! Internal T&S variables are potential temperature & practical salinity
    if (ids%id_sst > 0) call post_data(ids%id_sst, sfc_state%SST, diag, mask=g%mask2dT)
  endif
 
  if (tv%S_is_absS) then
    ! Internal T&S variables are conservative temperature & absolute salinity
    if (ids%id_sssabs > 0) call post_data(ids%id_sssabs, sfc_state%SSS, diag, mask=g%mask2dT)
    ! Use TEOS-10 function calls convert T&S diagnostics from absolute salinity
    ! to practical salinity.
    do j=js,je ; do i=is,ie
      work_2d(i,j) = gsw_sp_from_sr(sfc_state%SSS(i,j))
    enddo ; enddo
    if (ids%id_sss > 0) call post_data(ids%id_sss, work_2d, diag, mask=g%mask2dT)
  else
    ! Internal T&S variables are potential temperature & practical salinity
    if (ids%id_sss > 0) call post_data(ids%id_sss, sfc_state%SSS, diag, mask=g%mask2dT)
  endif
 
  if (ids%id_sst_sq > 0) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = sfc_state%SST(i,j)*sfc_state%SST(i,j)
    enddo ; enddo
    call post_data(ids%id_sst_sq, work_2d, diag, mask=g%mask2dT)
  endif
  if (ids%id_sss_sq > 0) then
    do j=js,je ; do i=is,ie
      work_2d(i,j) = sfc_state%SSS(i,j)*sfc_state%SSS(i,j)
    enddo ; enddo
    call post_data(ids%id_sss_sq, work_2d, diag, mask=g%mask2dT)
  endif
 
  call coupler_type_send_data(sfc_state%tr_fields, get_diag_time_end(diag))
 
end subroutine post_surface_thermo_diags
 
 
!> This routine posts diagnostics of the transports, including the subgridscale
!! contributions.
subroutine post_transport_diagnostics(G, GV, US, uhtr, vhtr, h, IDs, diag_pre_dyn, &
                                      diag, dt_trans, Reg)
  type(ocean_grid_type),    intent(inout) :: g   !< ocean grid structure
  type(verticalgrid_type),  intent(in)    :: gv  !< ocean vertical grid structure
  type(unit_scale_type),    intent(in)    :: us  !< A dimensional unit scaling type
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in) :: uhtr !< Accumulated zonal thickness fluxes
                                                 !! used to advect tracers [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in) :: vhtr !< Accumulated meridional thickness fluxes
                                                 !! used to advect tracers [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                            intent(in)    :: h   !< The updated layer thicknesses [H ~> m or kg m-2]
  type(transport_diag_ids), intent(in)    :: ids !< A structure with the diagnostic IDs.
  type(diag_grid_storage),  intent(inout) :: diag_pre_dyn !< Stored grids from before dynamics
  type(diag_ctrl),          intent(inout) :: diag !< regulates diagnostic output
  real,                     intent(in)    :: dt_trans !< total time step associated with the transports [T ~> s].
  type(tracer_registry_type), pointer     :: reg !< Pointer to the tracer registry
 
  ! Local variables
  real, dimension(SZIB_(G), SZJ_(G)) :: umo2d ! Diagnostics of integrated mass transport [R Z L2 T-1 ~> kg s-1]
  real, dimension(SZI_(G), SZJB_(G)) :: vmo2d ! Diagnostics of integrated mass transport [R Z L2 T-1 ~> kg s-1]
  real, dimension(SZIB_(G), SZJ_(G), SZK_(G)) :: umo ! Diagnostics of layer mass transport [R Z L2 T-1 ~> kg s-1]
  real, dimension(SZI_(G), SZJB_(G), SZK_(G)) :: vmo ! Diagnostics of layer mass transport [R Z L2 T-1 ~> kg s-1]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G))    :: h_tend ! Change in layer thickness due to dynamics
                          ! [H s-1 ~> m s-1 or kg m-2 s-1].
  real :: idt             ! The inverse of the time interval [T-1 ~> s-1]
  real :: h_to_rz_dt   ! A conversion factor from accumulated transports to fluxes
                          ! [R Z H-1 T-1 ~> kg m-3 s-1 or s-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
 
  idt = 1. / dt_trans
  h_to_rz_dt = gv%H_to_RZ * idt
 
  call diag_save_grids(diag)
  call diag_copy_storage_to_diag(diag, diag_pre_dyn)
 
  if (ids%id_umo_2d > 0) then
    umo2d(:,:) = 0.0
    do k=1,nz ; do j=js,je ; do i=is-1,ie
      umo2d(i,j) = umo2d(i,j) + uhtr(i,j,k) * h_to_rz_dt
    enddo ; enddo ; enddo
    call post_data(ids%id_umo_2d, umo2d, diag)
  endif
  if (ids%id_umo > 0) then
    ! Convert to kg/s.
    do k=1,nz ; do j=js,je ; do i=is-1,ie
      umo(i,j,k) = uhtr(i,j,k) * h_to_rz_dt
    enddo ; enddo ; enddo
    call post_data(ids%id_umo, umo, diag, alt_h = diag_pre_dyn%h_state)
  endif
  if (ids%id_vmo_2d > 0) then
    vmo2d(:,:) = 0.0
    do k=1,nz ; do j=js-1,je ; do i=is,ie
      vmo2d(i,j) = vmo2d(i,j) + vhtr(i,j,k) * h_to_rz_dt
    enddo ; enddo ; enddo
    call post_data(ids%id_vmo_2d, vmo2d, diag)
  endif
  if (ids%id_vmo > 0) then
    ! Convert to kg/s.
    do k=1,nz ; do j=js-1,je ; do i=is,ie
      vmo(i,j,k) = vhtr(i,j,k) * h_to_rz_dt
    enddo ; enddo ; enddo
    call post_data(ids%id_vmo, vmo, diag, alt_h = diag_pre_dyn%h_state)
  endif
 
  if (ids%id_uhtr > 0) call post_data(ids%id_uhtr, uhtr, diag, alt_h = diag_pre_dyn%h_state)
  if (ids%id_vhtr > 0) call post_data(ids%id_vhtr, vhtr, diag, alt_h = diag_pre_dyn%h_state)
  if (ids%id_dynamics_h > 0) call post_data(ids%id_dynamics_h, diag_pre_dyn%h_state, diag, &
                                            alt_h = diag_pre_dyn%h_state)
  ! Post the change in thicknesses
  if (ids%id_dynamics_h_tendency > 0) then
    h_tend(:,:,:) = 0.
    do k=1,nz ; do j=js,je ; do i=is,ie
      h_tend(i,j,k) = (h(i,j,k) - diag_pre_dyn%h_state(i,j,k))*idt
    enddo ; enddo ; enddo
    call post_data(ids%id_dynamics_h_tendency, h_tend, diag, alt_h = diag_pre_dyn%h_state)
  endif
 
  call post_tracer_transport_diagnostics(g, gv, reg, diag_pre_dyn%h_state, diag)
 
  call diag_restore_grids(diag)
 
end subroutine post_transport_diagnostics
 
!> This subroutine registers various diagnostics and allocates space for fields
!! that other diagnostis depend upon.
subroutine mom_diagnostics_init(MIS, ADp, CDp, Time, G, GV, US, param_file, diag, CS, tv)
  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(unit_scale_type),      intent(in)    :: us   !< A dimensional unit scaling type
  type(param_file_type),      intent(in)    :: param_file !< A structure to parse for run-time
                                                    !! parameters.
  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.
  type(thermo_var_ptrs),      intent(in)    :: tv   !< A structure pointing to various
                                                    !! thermodynamic variables.
 
  ! Local variables
  real :: omega, f2_min, convert_h
  ! This include declares and sets the variable "version".
# include "version_variable.h"
  character(len=40)  :: mdl = "MOM_diagnostics" ! This module's name.
  character(len=48) :: thickness_units, flux_units
  real :: wave_speed_min      ! A floor in the first mode speed below which 0 is returned [L T-1 ~> m s-1]
  real :: wave_speed_tol      ! The fractional tolerance for finding the wave speeds [nondim]
  logical :: better_speed_est ! If true, use a more robust estimate of the first
                              ! mode wave speed as the starting point for iterations.
  logical :: use_temperature, adiabatic
  logical :: default_2018_answers, remap_answers_2018
  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
  use_temperature = associated(tv%T)
  call get_param(param_file, mdl, "ADIABATIC", adiabatic, default=.false., &
                 do_not_log=.true.)
 
  ! 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 "// &
                 "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 "// &
                 "purposes of calculating the equivalent barotropic wave speed.", &
                 units='m', scale=us%m_to_Z, default=-1.)
  call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_TOL", wave_speed_tol, &
                 "The fractional tolerance for finding the wave speeds.", &
                 units="nondim", default=0.001)
  !### Set defaults so that wave_speed_min*wave_speed_tol >= 1e-9 m s-1
  call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_MIN", wave_speed_min, &
                 "A floor in the first mode speed below which 0 used instead.", &
                 units="m s-1", default=0.0, scale=us%m_s_to_L_T)
  call get_param(param_file, mdl, "INTERNAL_WAVE_SPEED_BETTER_EST", better_speed_est, &
                 "If true, use a more robust estimate of the first mode wave speed as the "//&
                 "starting point for iterations.", default=.false.) !### Change the default.
  call get_param(param_file, mdl, "DEFAULT_2018_ANSWERS", default_2018_answers, &
                 "This sets the default value for the various _2018_ANSWERS parameters.", &
                 default=.false.)
  call get_param(param_file, mdl, "REMAPPING_2018_ANSWERS", remap_answers_2018, &
                 "If true, use the order of arithmetic and expressions that recover the "//&
                 "answers from the end of 2018.  Otherwise, use updated and more robust "//&
                 "forms of the same expressions.", default=default_2018_answers)
 
  if (gv%Boussinesq) then
    thickness_units = "m" ; flux_units = "m3 s-1" ; convert_h = gv%H_to_m
  else
    thickness_units = "kg m-2" ; flux_units = "kg s-1" ; convert_h = gv%H_to_kg_m2
  endif
 
  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', conversion=us%Z_to_m, v_extensive=.true.)
  cs%id_h_pre_sync = register_diag_field('ocean_model', 'h_pre_sync', diag%axesTL, time, &
      long_name = 'Cell thickness from the previous timestep', &
      units='m', conversion=gv%H_to_m, v_extensive=.true.)
 
  ! Note that CS%id_volcello would normally be registered here but because it is a "cell measure" and
  ! must be registered first. We earlier stored the handle of volcello but need it here for posting
  ! by this module.
  cs%id_volcello = diag_get_volume_cell_measure_dm_id(diag)
 
  if (use_temperature) then
    if (tv%T_is_conT) then
      cs%id_Tpot = register_diag_field('ocean_model', 'temp', diag%axesTL, &
          time, 'Potential Temperature', 'degC')
    endif
    if (tv%S_is_absS) then
      cs%id_Sprac = register_diag_field('ocean_model', 'salt', diag%axesTL, &
          time, 'Salinity', 'psu')
    endif
 
    cs%id_tob = register_diag_field('ocean_model','tob', diag%axesT1, time, &
        long_name='Sea Water Potential Temperature at Sea Floor', &
        standard_name='sea_water_potential_temperature_at_sea_floor', units='degC')
    cs%id_sob = register_diag_field('ocean_model','sob',diag%axesT1, time, &
        long_name='Sea Water Salinity at Sea Floor', &
        standard_name='sea_water_salinity_at_sea_floor', units='psu')
 
    cs%id_temp_layer_ave = register_diag_field('ocean_model', 'temp_layer_ave', &
        diag%axesZL, time, 'Layer Average Ocean Temperature', 'degC')
    cs%id_salt_layer_ave = register_diag_field('ocean_model', 'salt_layer_ave', &
        diag%axesZL, time, 'Layer Average Ocean Salinity', 'psu')
 
    cs%id_thetaoga = register_scalar_field('ocean_model', 'thetaoga', &
        time, diag, 'Global Mean Ocean Potential Temperature', 'degC',&
        standard_name='sea_water_potential_temperature')
    cs%id_soga = register_scalar_field('ocean_model', 'soga', &
        time, diag, 'Global Mean Ocean Salinity', 'psu', &
        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_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='psu', standard_name='sea_surface_salinity',                      &
        cmor_field_name='sosga', cmor_standard_name='sea_surface_salinity',     &
        cmor_long_name='Sea Surface Salinity')
  endif
 
  cs%id_u = register_diag_field('ocean_model', 'u', diag%axesCuL, time,              &
      'Zonal velocity', 'm s-1', conversion=us%L_T_to_m_s, cmor_field_name='uo', &
      cmor_standard_name='sea_water_x_velocity', cmor_long_name='Sea Water X Velocity')
  cs%id_v = register_diag_field('ocean_model', 'v', diag%axesCvL, time,                  &
      'Meridional velocity', 'm s-1', conversion=us%L_T_to_m_s, cmor_field_name='vo', &
      cmor_standard_name='sea_water_y_velocity', cmor_long_name='Sea Water Y Velocity')
  cs%id_h = register_diag_field('ocean_model', 'h', diag%axesTL, time, &
      'Layer Thickness', thickness_units, v_extensive=.true., conversion=convert_h)
 
  cs%id_e = register_diag_field('ocean_model', 'e', diag%axesTi, time, &
      'Interface Height Relative to Mean Sea Level', 'm', conversion=us%Z_to_m)
  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', 'm', conversion=us%Z_to_m)
  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 m-3', conversion=us%R_to_kg_m3)
 
  cs%id_Rcv = register_diag_field('ocean_model', 'Rho_cv', diag%axesTL, time, &
      'Coordinate Potential Density', 'kg m-3', conversion=us%R_to_kg_m3)
 
  cs%id_rhopot0 = register_diag_field('ocean_model', 'rhopot0', diag%axesTL, time, &
      'Potential density referenced to surface', 'kg m-3', conversion=us%R_to_kg_m3)
  cs%id_rhopot2 = register_diag_field('ocean_model', 'rhopot2', diag%axesTL, time, &
      'Potential density referenced to 2000 dbar', 'kg m-3', conversion=us%R_to_kg_m3)
  cs%id_rhoinsitu = register_diag_field('ocean_model', 'rhoinsitu', diag%axesTL, time, &
      'In situ density', 'kg m-3', conversion=us%R_to_kg_m3)
  cs%id_drho_dT = register_diag_field('ocean_model', 'drho_dT', diag%axesTL, time, &
      'Partial derivative of rhoinsitu with respect to temperature (alpha)', 'kg m-3 degC-1')
  cs%id_drho_dS = register_diag_field('ocean_model', 'drho_dS', diag%axesTL, time, &
      'Partial derivative of rhoinsitu with respect to salinity (beta)', 'kg^2 g-1 m-3')
 
  cs%id_du_dt = register_diag_field('ocean_model', 'dudt', diag%axesCuL, time, &
      'Zonal Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
  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(lbound(mis%u), mis%u, cs%du_dt, cs)
  endif
 
  cs%id_dv_dt = register_diag_field('ocean_model', 'dvdt', diag%axesCvL, time, &
      'Meridional Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
  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(lbound(mis%v), 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)//" s-1", conversion=convert_h*us%s_to_T, v_extensive=.true.)
  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(lbound(mis%h), mis%h, cs%dh_dt, cs)
  endif
 
  !CS%id_hf_du_dt = register_diag_field('ocean_model', 'hf_dudt', diag%axesCuL, Time, &
  !    'Fractional Thickness-weighted Zonal Acceleration', 'm s-2', v_extensive=.true., &
  !    conversion=US%L_T2_to_m_s2)
  !if (CS%id_hf_du_dt > 0) then
  !  call safe_alloc_ptr(CS%hf_du_dt,IsdB,IedB,jsd,jed,nz)
  !  if (.not.associated(CS%du_dt)) then
  !    call safe_alloc_ptr(CS%du_dt,IsdB,IedB,jsd,jed,nz)
  !    call register_time_deriv(lbound(MIS%u), MIS%u, CS%du_dt, CS)
  !  endif
  !  call safe_alloc_ptr(ADp%diag_hfrac_u,IsdB,IedB,jsd,jed,nz)
  !endif
 
  !CS%id_hf_dv_dt = register_diag_field('ocean_model', 'hf_dvdt', diag%axesCvL, Time, &
  !    'Fractional Thickness-weighted Meridional Acceleration', 'm s-2', v_extensive=.true., &
  !    conversion=US%L_T2_to_m_s2)
  !if (CS%id_hf_dv_dt > 0) then
  !  call safe_alloc_ptr(CS%hf_dv_dt,isd,ied,JsdB,JedB,nz)
  !  if (.not.associated(CS%dv_dt)) then
  !    call safe_alloc_ptr(CS%dv_dt,isd,ied,JsdB,JedB,nz)
  !    call register_time_deriv(lbound(MIS%v), MIS%v, CS%dv_dt, CS)
  !  endif
  !  call safe_alloc_ptr(ADp%diag_hfrac_v,isd,ied,Jsd,JedB,nz)
  !endif
 
  cs%id_hf_du_dt_2d = register_diag_field('ocean_model', 'hf_dudt_2d', diag%axesCu1, time, &
      'Depth-sum Fractional Thickness-weighted Zonal Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
  if (cs%id_hf_du_dt_2d > 0) then
    if (.not.associated(cs%du_dt)) then
      call safe_alloc_ptr(cs%du_dt,isdb,iedb,jsd,jed,nz)
      call register_time_deriv(lbound(mis%u), mis%u, cs%du_dt, cs)
    endif
    call safe_alloc_ptr(adp%diag_hfrac_u,isdb,iedb,jsd,jed,nz)
  endif
 
  cs%id_hf_dv_dt_2d = register_diag_field('ocean_model', 'hf_dvdt_2d', diag%axesCv1, time, &
      'Depth-sum Fractional Thickness-weighted Meridional Acceleration', 'm s-2', conversion=us%L_T2_to_m_s2)
  if (cs%id_hf_dv_dt_2d > 0) then
    if (.not.associated(cs%dv_dt)) then
      call safe_alloc_ptr(cs%dv_dt,isd,ied,jsdb,jedb,nz)
      call register_time_deriv(lbound(mis%v), mis%v, cs%dv_dt, cs)
    endif
    call safe_alloc_ptr(adp%diag_hfrac_v,isd,ied,jsd,jedb,nz)
  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, conversion=convert_h)
    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, conversion=us%L_to_m**2*us%s_to_T*convert_h)
    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, conversion=us%L_to_m**2*us%s_to_T*convert_h)
    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, conversion=us%L_to_m**2*us%s_to_T*convert_h)
    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, conversion=us%L_to_m**2*us%s_to_T*convert_h)
    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', &
      'm2 s-2', conversion=us%L_T_to_m_s**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', &
      'm3 s-3', conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
  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', &
      'm3 s-3', conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
  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', &
      'm3 s-3', conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
  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', &
      'm3 s-3', conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
  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', &
      'm3 s-3', conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
  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', &
      'm3 s-3', conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
  if (cs%id_KE_horvisc>0) call safe_alloc_ptr(cs%KE_horvisc,isd,ied,jsd,jed,nz)
 
  if (.not. adiabatic) then
    cs%id_KE_dia = register_diag_field('ocean_model', 'KE_dia', diag%axesTL, time, &
        'Kinetic Energy Source from Diapycnal Diffusion', &
        'm3 s-3', conversion=gv%H_to_m*(us%L_T_to_m_s**2)*us%s_to_T)
    if (cs%id_KE_dia>0) call safe_alloc_ptr(cs%KE_dia,isd,ied,jsd,jed,nz)
  endif
 
  ! gravity wave CFLs
  cs%id_cg1 = register_diag_field('ocean_model', 'cg1', diag%axesT1, time, &
      'First baroclinic gravity wave speed', 'm s-1', conversion=us%L_T_to_m_s)
  cs%id_Rd1 = register_diag_field('ocean_model', 'Rd1', diag%axesT1, time, &
      'First baroclinic deformation radius', 'm', conversion=us%L_to_m)
  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', 'm s-1', conversion=us%L_T_to_m_s)
  cs%id_Rd_ebt = register_diag_field('ocean_model', 'Rd_ebt', diag%axesT1, time, &
      'Equivalent barotropic deformation radius', 'm', conversion=us%L_to_m)
  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, remap_answers_2018=remap_answers_2018, &
                         better_speed_est=better_speed_est, min_speed=wave_speed_min, &
                         wave_speed_tol=wave_speed_tol)
    call wave_speed_init(cs%wave_speed_CSp, remap_answers_2018=remap_answers_2018)
    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', conversion=us%RZ_to_kg_m2)
 
  if (use_temperature) then
    cs%id_temp_int = register_diag_field('ocean_model', 'temp_int', diag%axesT1, time,                &
        'Density weighted column integrated potential temperature', 'degC kg m-2', conversion=us%RZ_to_kg_m2, &
        cmor_field_name='opottempmint',                                                               &
        cmor_long_name='integral_wrt_depth_of_product_of_sea_water_density_and_potential_temperature',&
        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', 'psu kg m-2', conversion=us%RZ_to_kg_m2, &
        cmor_field_name='somint',                                                        &
        cmor_long_name='integral_wrt_depth_of_product_of_sea_water_density_and_salinity',&
        cmor_standard_name='Depth integrated density times salinity')
  endif
 
  cs%id_col_mass = register_diag_field('ocean_model', 'col_mass', diag%axesT1, time, &
      'The column integrated in situ density', 'kg m-2', conversion=us%RZ_to_kg_m2)
 
  cs%id_col_ht = register_diag_field('ocean_model', 'col_height', diag%axesT1, time, &
      'The height of the water column', 'm', conversion=us%Z_to_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', conversion=us%RL2_T2_to_Pa)
 
  call set_dependent_diagnostics(mis, adp, cdp, g, cs)
 
end subroutine mom_diagnostics_init
 
 
!> Register diagnostics of the surface state and integrated quantities
subroutine register_surface_diags(Time, G, US, IDs, diag, tv)
  type(time_type),         intent(in)    :: time  !< current model time
  type(ocean_grid_type),   intent(in)    :: g     !< ocean grid structure
  type(unit_scale_type),   intent(in)    :: us   !< A dimensional unit scaling type
  type(surface_diag_ids),  intent(inout) :: ids   !< A structure with the diagnostic IDs.
  type(diag_ctrl),         intent(inout) :: diag  !< regulates diagnostic output
  type(thermo_var_ptrs),   intent(in)    :: tv    !< A structure pointing to various thermodynamic variables
 
  ! Vertically integrated, budget, and surface state diagnostics
  ids%id_volo = register_scalar_field('ocean_model', 'volo', time, diag,&
      long_name='Total volume of liquid ocean', units='m3',            &
      standard_name='sea_water_volume')
  ids%id_zos = register_diag_field('ocean_model', 'zos', diag%axesT1, time,&
      standard_name = 'sea_surface_height_above_geoid',                   &
      long_name= 'Sea surface height above geoid', units='m')
  ids%id_zossq = register_diag_field('ocean_model', 'zossq', diag%axesT1, time,&
      standard_name='square_of_sea_surface_height_above_geoid',             &
      long_name='Square of sea surface height above geoid', units='m2')
  ids%id_ssh = register_diag_field('ocean_model', 'SSH', diag%axesT1, time, &
      'Sea Surface Height', 'm')
  ids%id_ssh_ga = register_scalar_field('ocean_model', 'ssh_ga', time, diag,&
      long_name='Area averaged sea surface height', units='m',            &
      standard_name='area_averaged_sea_surface_height')
  ids%id_ssu = register_diag_field('ocean_model', 'SSU', diag%axesCu1, time, &
      'Sea Surface Zonal Velocity', 'm s-1', conversion=us%L_T_to_m_s)
  ids%id_ssv = register_diag_field('ocean_model', 'SSV', diag%axesCv1, time, &
      'Sea Surface Meridional Velocity', 'm s-1', conversion=us%L_T_to_m_s)
  ids%id_speed = register_diag_field('ocean_model', 'speed', diag%axesT1, time, &
      'Sea Surface Speed', 'm s-1', conversion=us%L_T_to_m_s)
 
  if (associated(tv%T)) then
    ids%id_sst = register_diag_field('ocean_model', 'SST', diag%axesT1, time,     &
        'Sea Surface Temperature', 'degC', cmor_field_name='tos', &
        cmor_long_name='Sea Surface Temperature',                                &
        cmor_standard_name='sea_surface_temperature')
    ids%id_sst_sq = register_diag_field('ocean_model', 'SST_sq', diag%axesT1, time, &
        'Sea Surface Temperature Squared', 'degC2', cmor_field_name='tossq', &
        cmor_long_name='Square of Sea Surface Temperature ',                      &
        cmor_standard_name='square_of_sea_surface_temperature')
    ids%id_sss = register_diag_field('ocean_model', 'SSS', diag%axesT1, time, &
        'Sea Surface Salinity', 'psu', cmor_field_name='sos', &
        cmor_long_name='Sea Surface Salinity',                            &
        cmor_standard_name='sea_surface_salinity')
    ids%id_sss_sq = register_diag_field('ocean_model', 'SSS_sq', diag%axesT1, time, &
        'Sea Surface Salinity Squared', 'psu', cmor_field_name='sossq', &
        cmor_long_name='Square of Sea Surface Salinity ',                     &
        cmor_standard_name='square_of_sea_surface_salinity')
    if (tv%T_is_conT) then
      ids%id_sstcon = register_diag_field('ocean_model', 'conSST', diag%axesT1, time,     &
          'Sea Surface Conservative Temperature', 'Celsius')
    endif
    if (tv%S_is_absS) then
      ids%id_sssabs = register_diag_field('ocean_model', 'absSSS', diag%axesT1, time,     &
          'Sea Surface Absolute Salinity', 'g kg-1')
    endif
    if (associated(tv%frazil)) then
      ids%id_fraz = register_diag_field('ocean_model', 'frazil', diag%axesT1, time, &
            'Heat from frazil formation', 'W m-2', conversion=us%QRZ_T_to_W_m2, &
            cmor_field_name='hfsifrazil', &
            cmor_standard_name='heat_flux_into_sea_water_due_to_frazil_ice_formation', &
            cmor_long_name='Heat Flux into Sea Water due to Frazil Ice Formation')
    endif
  endif
 
  ids%id_salt_deficit = register_diag_field('ocean_model', 'salt_deficit', diag%axesT1, time, &
         'Salt sink in ocean due to ice flux', &
         'psu m-2 s-1', conversion=us%RZ_T_to_kg_m2s)
  ids%id_Heat_PmE = register_diag_field('ocean_model', 'Heat_PmE', diag%axesT1, time, &
         'Heat flux into ocean from mass flux into ocean', &
         'W m-2', conversion=us%QRZ_T_to_W_m2)
  ids%id_intern_heat = register_diag_field('ocean_model', 'internal_heat', diag%axesT1, time,&
         'Heat flux into ocean from geothermal or other internal sources', &
         'W m-2', conversion=us%QRZ_T_to_W_m2)
 
end subroutine register_surface_diags
 
!> Register certain diagnostics related to transports
subroutine register_transport_diags(Time, G, GV, US, IDs, diag)
  type(time_type),          intent(in)    :: time  !< current model time
  type(ocean_grid_type),    intent(in)    :: g     !< ocean grid structure
  type(verticalgrid_type),  intent(in)    :: gv    !< ocean vertical grid structure
  type(unit_scale_type),    intent(in)    :: us   !< A dimensional unit scaling type
  type(transport_diag_ids), intent(inout) :: ids   !< A structure with the diagnostic IDs.
  type(diag_ctrl),          intent(inout) :: diag  !< regulates diagnostic output
 
  real :: h_convert
  character(len=48) :: thickness_units
 
  thickness_units = get_thickness_units(gv)
  if (gv%Boussinesq) then
    h_convert = gv%H_to_m
  else
    h_convert = gv%H_to_kg_m2
  endif
 
  ! Diagnostics related to tracer and mass transport
  ids%id_uhtr = register_diag_field('ocean_model', 'uhtr', diag%axesCuL, time, &
      'Accumulated zonal thickness fluxes to advect tracers', 'kg', &
      y_cell_method='sum', v_extensive=.true., conversion=h_convert*us%L_to_m**2)
  ids%id_vhtr = register_diag_field('ocean_model', 'vhtr', diag%axesCvL, time, &
      'Accumulated meridional thickness fluxes to advect tracers', 'kg', &
      x_cell_method='sum', v_extensive=.true., conversion=h_convert*us%L_to_m**2)
  ids%id_umo = register_diag_field('ocean_model', 'umo', &
      diag%axesCuL, time, 'Ocean Mass X Transport', &
      'kg s-1', conversion=us%RZ_T_to_kg_m2s*us%L_to_m**2, &
      standard_name='ocean_mass_x_transport', y_cell_method='sum', v_extensive=.true.)
  ids%id_vmo = register_diag_field('ocean_model', 'vmo', &
      diag%axesCvL, time, 'Ocean Mass Y Transport', &
      'kg s-1', conversion=us%RZ_T_to_kg_m2s*us%L_to_m**2, &
      standard_name='ocean_mass_y_transport', x_cell_method='sum', v_extensive=.true.)
  ids%id_umo_2d = register_diag_field('ocean_model', 'umo_2d', &
      diag%axesCu1, time, 'Ocean Mass X Transport Vertical Sum', &
      'kg s-1', conversion=us%RZ_T_to_kg_m2s*us%L_to_m**2, &
      standard_name='ocean_mass_x_transport_vertical_sum', y_cell_method='sum')
  ids%id_vmo_2d = register_diag_field('ocean_model', 'vmo_2d', &
      diag%axesCv1, time, 'Ocean Mass Y Transport Vertical Sum', &
      'kg s-1', conversion=us%RZ_T_to_kg_m2s*us%L_to_m**2, &
      standard_name='ocean_mass_y_transport_vertical_sum', x_cell_method='sum')
  ids%id_dynamics_h = register_diag_field('ocean_model','dynamics_h',  &
      diag%axesTl, time, 'Change in layer thicknesses due to horizontal dynamics', &
      'm s-1', v_extensive=.true., conversion=gv%H_to_m)
  ids%id_dynamics_h_tendency = register_diag_field('ocean_model','dynamics_h_tendency',  &
      diag%axesTl, time, 'Change in layer thicknesses due to horizontal dynamics', &
      'm s-1', v_extensive=.true., conversion=gv%H_to_m*us%s_to_T)
 
end subroutine register_transport_diags
 
!> Offers the static fields in the ocean grid type for output via the diag_manager.
subroutine write_static_fields(G, GV, US, tv, diag)
  type(ocean_grid_type),   intent(in)    :: g    !< ocean grid structure
  type(verticalgrid_type), intent(in)    :: gv   !< ocean vertical grid structure
  type(unit_scale_type),   intent(in)    :: us   !< A dimensional unit scaling type
  type(thermo_var_ptrs),   intent(in)    :: tv   !< A structure pointing to various thermodynamic variables
  type(diag_ctrl), target, intent(inout) :: diag !< regulates diagnostic output
 
  ! Local variables
  integer :: id
  logical :: use_temperature
 
  id = register_static_field('ocean_model', 'geolat', diag%axesT1, &
        'Latitude of tracer (T) points', 'degrees_north')
  if (id > 0) call post_data(id, g%geoLatT, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolon', diag%axesT1, &
        'Longitude of tracer (T) points', 'degrees_east')
  if (id > 0) call post_data(id, g%geoLonT, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolat_c', diag%axesB1, &
        'Latitude of corner (Bu) points', 'degrees_north', interp_method='none')
  if (id > 0) call post_data(id, g%geoLatBu, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolon_c', diag%axesB1, &
        'Longitude of corner (Bu) points', 'degrees_east', interp_method='none')
  if (id > 0) call post_data(id, g%geoLonBu, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolat_v', diag%axesCv1, &
        'Latitude of meridional velocity (Cv) points', 'degrees_north', interp_method='none')
  if (id > 0) call post_data(id, g%geoLatCv, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolon_v', diag%axesCv1, &
        'Longitude of meridional velocity (Cv) points', 'degrees_east', interp_method='none')
  if (id > 0) call post_data(id, g%geoLonCv, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolat_u', diag%axesCu1, &
        'Latitude of zonal velocity (Cu) points', 'degrees_north', interp_method='none')
  if (id > 0) call post_data(id, g%geoLatCu, diag, .true.)
 
  id = register_static_field('ocean_model', 'geolon_u', diag%axesCu1, &
        'Longitude of zonal velocity (Cu) points', 'degrees_east', interp_method='none')
  if (id > 0) call post_data(id, g%geoLonCu, diag, .true.)
 
  id = register_static_field('ocean_model', 'area_t', diag%axesT1, &
        'Surface area of tracer (T) cells', 'm2', conversion=us%L_to_m**2, &
        cmor_field_name='areacello', cmor_standard_name='cell_area', &
        cmor_long_name='Ocean Grid-Cell Area', &
        x_cell_method='sum', y_cell_method='sum', area_cell_method='sum')
  if (id > 0) then
    call post_data(id, g%areaT, diag, .true.)
    call diag_register_area_ids(diag, id_area_t=id)
  endif
 
  id = register_static_field('ocean_model', 'area_u', diag%axesCu1, &
        'Surface area of x-direction flow (U) cells', 'm2', conversion=us%L_to_m**2, &
        cmor_field_name='areacello_cu', cmor_standard_name='cell_area', &
        cmor_long_name='Ocean Grid-Cell Area', &
        x_cell_method='sum', y_cell_method='sum', area_cell_method='sum')
  if (id > 0) call post_data(id, g%areaCu, diag, .true.)
 
  id = register_static_field('ocean_model', 'area_v', diag%axesCv1, &
        'Surface area of y-direction flow (V) cells', 'm2', conversion=us%L_to_m**2, &
        cmor_field_name='areacello_cv', cmor_standard_name='cell_area', &
        cmor_long_name='Ocean Grid-Cell Area', &
        x_cell_method='sum', y_cell_method='sum', area_cell_method='sum')
  if (id > 0) call post_data(id, g%areaCv, diag, .true.)
 
  id = register_static_field('ocean_model', 'area_q', diag%axesB1, &
        'Surface area of B-grid flow (Q) cells', 'm2', conversion=us%L_to_m**2, &
        cmor_field_name='areacello_bu', cmor_standard_name='cell_area', &
        cmor_long_name='Ocean Grid-Cell Area', &
        x_cell_method='sum', y_cell_method='sum', area_cell_method='sum')
  if (id > 0) call post_data(id, g%areaBu, diag, .true.)
 
  id = register_static_field('ocean_model', 'depth_ocean', diag%axesT1,  &
        'Depth of the ocean at tracer points', 'm', conversion=us%Z_to_m, &
        standard_name='sea_floor_depth_below_geoid',                     &
        cmor_field_name='deptho', cmor_long_name='Sea Floor Depth',      &
        cmor_standard_name='sea_floor_depth_below_geoid', area=diag%axesT1%id_area, &
        x_cell_method='mean', y_cell_method='mean', area_cell_method='mean')
  if (id > 0) call post_data(id, g%bathyT, diag, .true., mask=g%mask2dT)
 
  id = register_static_field('ocean_model', 'wet', diag%axesT1, &
        '0 if land, 1 if ocean at tracer points', 'none', area=diag%axesT1%id_area)
  if (id > 0) call post_data(id, g%mask2dT, diag, .true.)
 
  id = register_static_field('ocean_model', 'wet_c', diag%axesB1, &
        '0 if land, 1 if ocean at corner (Bu) points', 'none', interp_method='none')
  if (id > 0) call post_data(id, g%mask2dBu, diag, .true.)
 
  id = register_static_field('ocean_model', 'wet_u', diag%axesCu1, &
        '0 if land, 1 if ocean at zonal velocity (Cu) points', 'none', interp_method='none')
  if (id > 0) call post_data(id, g%mask2dCu, diag, .true.)
 
  id = register_static_field('ocean_model', 'wet_v', diag%axesCv1, &
        '0 if land, 1 if ocean at meridional velocity (Cv) points', 'none', interp_method='none')
  if (id > 0) call post_data(id, g%mask2dCv, diag, .true.)
 
  id = register_static_field('ocean_model', 'Coriolis', diag%axesB1, &
        'Coriolis parameter at corner (Bu) points', 's-1', interp_method='none', conversion=us%s_to_T)
  if (id > 0) call post_data(id, g%CoriolisBu, diag, .true.)
 
  id = register_static_field('ocean_model', 'dxt', diag%axesT1, &
        'Delta(x) at thickness/tracer points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dxT, diag, .true.)
 
  id = register_static_field('ocean_model', 'dyt', diag%axesT1, &
        'Delta(y) at thickness/tracer points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dyT, diag, .true.)
 
  id = register_static_field('ocean_model', 'dxCu', diag%axesCu1, &
        'Delta(x) at u points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dxCu, diag, .true.)
 
  id = register_static_field('ocean_model', 'dyCu', diag%axesCu1, &
        'Delta(y) at u points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dyCu, diag, .true.)
 
  id = register_static_field('ocean_model', 'dxCv', diag%axesCv1, &
        'Delta(x) at v points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dxCv, diag, .true.)
 
  id = register_static_field('ocean_model', 'dyCv', diag%axesCv1, &
        'Delta(y) at v points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dyCv, diag, .true.)
 
  id = register_static_field('ocean_model', 'dyCuo', diag%axesCu1, &
        'Open meridional grid spacing at u points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dy_Cu, diag, .true.)
 
  id = register_static_field('ocean_model', 'dxCvo', diag%axesCv1, &
        'Open zonal grid spacing at v points (meter)', 'm', interp_method='none', conversion=us%L_to_m)
  if (id > 0) call post_data(id, g%dx_Cv, diag, .true.)
 
  id = register_static_field('ocean_model', 'sin_rot', diag%axesT1, &
        'sine of the clockwise angle of the ocean grid north to true north', 'none')
  if (id > 0) call post_data(id, g%sin_rot, diag, .true.)
 
  id = register_static_field('ocean_model', 'cos_rot', diag%axesT1, &
        'cosine of the clockwise angle of the ocean grid north to true north', 'none')
  if (id > 0) call post_data(id, g%cos_rot, diag, .true.)
 
 
  ! This static diagnostic is from CF 1.8, and is the fraction of a cell
  ! covered by ocean, given as a percentage (poorly named).
  id = register_static_field('ocean_model', 'area_t_percent', diag%axesT1, &
        'Percentage of cell area covered by ocean', '%', conversion=100.0, &
        cmor_field_name='sftof', cmor_standard_name='SeaAreaFraction', &
        cmor_long_name='Sea Area Fraction', &
        x_cell_method='mean', y_cell_method='mean', area_cell_method='mean')
  if (id > 0) call post_data(id, g%mask2dT, diag, .true.)
 
 
  id = register_static_field('ocean_model','Rho_0', diag%axesNull, &
       'mean ocean density used with the Boussinesq approximation', &
       'kg m-3', conversion=us%R_to_kg_m3, cmor_field_name='rhozero', &
       cmor_standard_name='reference_sea_water_density_for_boussinesq_approximation', &
       cmor_long_name='reference sea water density for boussinesq approximation')
  if (id > 0) call post_data(id, gv%Rho0, diag, .true.)
 
  use_temperature = associated(tv%T)
  if (use_temperature) then
    id = register_static_field('ocean_model','C_p', diag%axesNull, &
         'heat capacity of sea water', 'J kg-1 K-1', conversion=us%Q_to_J_kg, &
         cmor_field_name='cpocean', &
         cmor_standard_name='specific_heat_capacity_of_sea_water', &
         cmor_long_name='specific_heat_capacity_of_sea_water')
    if (id > 0) call post_data(id, tv%C_p, diag, .true.)
  endif
 
end subroutine write_static_fields
 
 
!> 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.
  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(lbound(mis%u), 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(lbound(mis%v), 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(lbound(mis%h), 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
 
!> Deallocate memory associated with the diagnostics module
subroutine mom_diagnostics_end(CS, ADp)
  type(diagnostics_cs),   pointer       :: cs  !< Control structure returned by a
                                               !! previous call to diagnostics_init.
  type(accel_diag_ptrs),  intent(inout) :: adp !< structure with pointers to
                                               !! accelerations in momentum equation.
  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(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)
 
  if (associated(adp%diag_hfrac_u)) deallocate(adp%diag_hfrac_u)
  if (associated(adp%diag_hfrac_v)) deallocate(adp%diag_hfrac_v)
 
  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