MOM_diag_to_Z.F90

Go to the documentation of this file.

module mom_diag_to_z

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

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  By Robert Hallberg, July 2006                                      *
!*                                                                     *
!*    This subroutine maps tracers and velocities into depth space     *
!*  for output as diagnostic quantities.  Currently, a piecewise       *
!*  linear subgrid structure is used for tracers, while velocities can *
!*  use either piecewise constant or piecewise linear structures.      *
!*                                                                     *
!*     A small fragment of the grid is shown below:                    *
!*                                                                     *
!*    j+1  x ^ x ^ x   At x:  q, CoriolisBu                            *
!*    j+1  > o > o >   At ^:  v                                        *
!*    j    x ^ x ^ x   At >:  u                                        *
!*    j    > o > o >   At o:  h, bathyT                                *
!*    j-1  x ^ x ^ x                                                   *
!*        i-1  i  i+1  At x & ^:                                       *
!*           i  i+1    At > & o:                                       *
!*                                                                     *
!*  The boundaries always run through q grid points (x).               *
!*                                                                     *
!********+*********+*********+*********+*********+*********+*********+**
use mom_domains,       only : pass_var
use mom_coms,          only : reproducing_sum
use mom_diag_mediator, only : post_data, post_data_1d_k, register_diag_field, safe_alloc_ptr
use mom_diag_mediator, only : diag_ctrl, time_type, diag_axis_init
use mom_diag_mediator, only : axes_grp, define_axes_group
use mom_diag_mediator, only : ocean_register_diag
use mom_error_handler, only : mom_error, fatal, warning
use mom_file_parser,   only : get_param, log_param, log_version, param_file_type
use mom_grid,          only : ocean_grid_type
use mom_io,            only : slasher, vardesc, query_vardesc, modify_vardesc
use mom_spatial_means, only : global_layer_mean
use mom_variables,     only : p3d, p2d
use mom_verticalgrid,  only : verticalgrid_type

use netcdf

implicit none ; private

#include <MOM_memory.h>

public calculate_z_diag_fields
public register_z_tracer
public mom_diag_to_z_init
public calculate_z_transport
public mom_diag_to_z_end
public ocean_register_diag_with_z
public find_overlap
public find_limited_slope
public register_zint_diag
public calc_zint_diags

type, public :: diag_to_z_cs ; private
  ! The following arrays are used to store diagnostics calculated in this
  ! module and unavailable outside of it.

  real, pointer, dimension(:,:,:) :: &
    u_z  => null(), &   ! zonal velocity remapped to depth space (m/s)
    v_z  => null(), &   ! meridional velocity remapped to depth space (m/s)
    uh_z => null(), &   ! zonal transport remapped to depth space (m3/s or kg/s)
    vh_z => null()      ! meridional transport remapped to depth space (m3/s or kg/s)

  type(p3d) :: tr_z(max_fields_)     ! array of tracers, remapped to depth space
  type(p3d) :: tr_model(max_fields_) ! pointers to an array of tracers

  real    :: missing_vel             = -1.0e34
  real    :: missing_trans           = -1.0e34
  real    :: missing_value           = -1.0e34
  real    :: missing_tr(max_fields_) = -1.0e34

  integer :: id_u_z  = -1
  integer :: id_v_z  = -1
  integer :: id_uh_z = -1
  integer :: id_vh_z = -1
  integer :: id_tr(max_fields_) = -1
  integer :: id_tr_xyave(max_fields_) = -1
  integer :: num_tr_used = 0
  integer :: nk_zspace = -1

  real, pointer :: z_int(:) => null()  ! interface depths of the z-space file (meter)

  type(axes_grp) :: axesbz,  axestz,  axescuz,  axescvz
  type(axes_grp) :: axesbzi, axestzi, axescuzi, axescvzi
  type(axes_grp) :: axesz
  integer, dimension(1) :: axesz_out

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

end type diag_to_z_cs

integer, parameter :: no_zspace = -1

contains

function global_z_mean(var,G,CS,tracer)
  type(ocean_grid_type),                           intent(in)  :: G    !< The ocean's grid structure
  type(diag_to_z_cs),                              intent(in)  :: CS
  real, dimension(SZI_(G), SZJ_(G), CS%nk_zspace), intent(in)  :: var
  real, dimension(SZI_(G), SZJ_(G), CS%nk_zspace)  :: tmpForSumming, weight
  real, dimension(SZI_(G), SZJ_(G), CS%nk_zspace)  :: localVar, valid_point, depth_weight
  real, dimension(CS%nk_zspace)                    :: global_z_mean, scalarij, weightij
  real, dimension(CS%nk_zspace)                    :: global_temp_scalar, global_weight_scalar
  integer :: i, j, k, is, ie, js, je, nz, tracer
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ;
  nz = cs%nk_zspace

  ! Initialize local arrays
  valid_point = 1. ; depth_weight = 0.
  tmpforsumming(:,:,:) = 0. ; weight(:,:,:) = 0.

  ! Local array copy of tracer field pointer
  localvar = var

  do k=1, nz ; do j=js,je ; do i=is, ie
    ! Weight factor for partial bottom cells
    depth_weight(i,j,k) = min( max( (-1.*g%bathyT(i,j)), cs%Z_int(k+1) ) - cs%Z_int(k), 0.)

    ! Flag the point as invalid if it contains missing data, or is below the bathymetry
    if (var(i,j,k) == cs%missing_tr(tracer)) valid_point(i,j,k) = 0.
    if (depth_weight(i,j,k) == 0.) valid_point(i,j,k) = 0.

    ! If the point is flagged, set the variable itsef to zero to avoid NaNs
    if (valid_point(i,j,k) == 0.) localvar(i,j,k) = 0.

    weight(i,j,k)  =  depth_weight(i,j,k) * ( (valid_point(i,j,k) * (g%areaT(i,j) * g%mask2dT(i,j))) )
    tmpforsumming(i,j,k) =  localvar(i,j,k) * weight(i,j,k)
  enddo ; enddo ; enddo

  global_temp_scalar   = reproducing_sum(tmpforsumming,sums=scalarij)
  global_weight_scalar = reproducing_sum(weight,sums=weightij)

  do k=1, nz
    if (scalarij(k) == 0) then
        global_z_mean(k) = 0.0
    else
        global_z_mean(k) = scalarij(k) / weightij(k)
    endif
  enddo

end function global_z_mean

!> This subroutine maps tracers and velocities into depth space for diagnostics.
subroutine calculate_z_diag_fields(u, v, h, ssh_in, frac_shelf_h, G, GV, CS)
  type(ocean_grid_type),                     intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),                   intent(in)    :: GV   !< The ocean's vertical grid
                                                                   !! structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: u    !< The zonal velocity, in m s-1.
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: v    !< The meridional velocity,
                                                                   !! in m s-1.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h    !< Layer thicknesses, in H
                                                                   !! (usually m or kg m-2).
  real, dimension(SZI_(G),SZJ_(G)),          intent(in)    :: ssh_in !< Sea surface height
                                                                   !! (meter or kg/m2).
  real, dimension(:,:),                      pointer       :: frac_shelf_h
  type(diag_to_z_cs),                        pointer       :: CS   !< Control structure returned by
                                                                   !! previous call to
                                                                   !! diagnostics_init.

! This subroutine maps tracers and velocities into depth space for diagnostics.

! Arguments:
!  (in)  u   - zonal velocity component (m/s)
!  (in)  v   - meridional velocity component (m/s)
!  (in)  h   - layer thickness (meter or kg/m2)
!  (in)  ssh_in - sea surface height (meter or kg/m2)
!  (in)  G   - ocean grid structure
!  (in)  GV  - The ocean's vertical grid structure.
!  (in)  CS  - control structure returned by previous call to diagnostics_init

  ! Note the deliberately reversed axes in h_f, u_f, v_f, and tr_f.
  real :: ssh(szi_(g),szj_(g))   ! copy of ssh_in (meter or kg/m2)
  real :: e(szk_(g)+2)           ! z-star interface heights (meter or kg/m2)
  real :: h_f(szk_(g)+1,szi_(g)) ! thicknesses of massive layers (meter or kg/m2)
  real :: u_f(szk_(g)+1,szib_(g))! zonal velocity component in any massive layer
  real :: v_f(szk_(g)+1,szi_(g)) ! meridional velocity component in any massive layer

  real :: tr_f(szk_(g),max(cs%num_tr_used,1),szi_(g)) ! tracer concentration in massive layers
  integer :: nk_valid(szib_(g))  ! number of massive layers in a column

  real :: D_pt(szib_(g))        ! bottom depth (meter or kg/m2)
  real :: shelf_depth(szib_(g)) ! ice shelf depth (meter or kg/m2)
  real :: htot           ! summed layer thicknesses (meter or kg/m2)
  real :: dilate         ! proportion by which to dilate every layer
  real :: wt(szk_(g)+1)  ! fractional weight for each layer in the
                         ! range between k_top and k_bot (nondim)
  real :: z1(szk_(g)+1)  ! z1 and z2 are the depths of the top and bottom
  real :: z2(szk_(g)+1)  ! limits of the part of a layer that contributes
                         ! to a depth level, relative to the cell center
                         ! and normalized by the cell thickness (nondim)
                         ! Note that -1/2 <= z1 < z2 <= 1/2.
  real :: sl_tr(max(cs%num_tr_used,1)) ! normalized slope of the tracer
                                       ! within the cell, in tracer units
  real :: Angstrom ! A minimal layer thickness, in H.
  real :: slope ! normalized slope of a variable within the cell

  real :: layer_ave(cs%nk_zspace)

  logical :: linear_velocity_profiles, ice_shelf

  integer :: k_top, k_bot, k_bot_prev
  integer :: i, j, k, k2, kz, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nk, m, nkml
  integer :: IsgB, IegB, JsgB, JegB
  is   = g%isc  ; ie  = g%iec  ; js  = g%jsc  ; je  = g%jec ; nk = g%ke
  isq  = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB
  isgb = g%IsgB ; iegb = g%IegB ; jsgb = g%JsgB ; jegb = g%JegB
  nkml = max(gv%nkml, 1)
  angstrom = gv%Angstrom
  ssh(:,:) = ssh_in
  linear_velocity_profiles = .true.
  ! Update the halos
  call pass_var(ssh, g%Domain)

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

  ice_shelf = associated(frac_shelf_h)

  ! If no fields are needed, return
  if ((cs%id_u_z <= 0) .and. (cs%id_v_z <= 0) .and. (cs%num_tr_used < 1)) return

  ! zonal velocity component
  if (cs%id_u_z > 0) then

    do kz=1,cs%nk_zspace ; do j=js,je ; do i=isq,ieq
      cs%u_z(i,j,kz) = cs%missing_vel
    enddo ; enddo ; enddo


    do j=js,je
      shelf_depth(:) = 0. ! initially all is open ocean
      ! Remove all massless layers.
      do i=isq,ieq
        nk_valid(i) = 0
        d_pt(i) = 0.5*(g%bathyT(i+1,j)+g%bathyT(i,j))
        if (ice_shelf) then
          if (frac_shelf_h(i,j)+frac_shelf_h(i+1,j) > 0.) then ! under shelf
            shelf_depth(i) = abs(0.5*(ssh(i+1,j)+ssh(i,j)))
          endif
        endif
      enddo
      do k=1,nk ; do i=isq,ieq
        if ((g%mask2dCu(i,j) > 0.5) .and. (h(i,j,k)+h(i+1,j,k) > 4.0*angstrom)) then
          nk_valid(i) = nk_valid(i) + 1 ; k2 = nk_valid(i)
          h_f(k2,i) = 0.5*(h(i,j,k)+h(i+1,j,k)) ; u_f(k2,i) = u(i,j,k)
        endif
      enddo ; enddo
      do i=isq,ieq ; if (g%mask2dCu(i,j) > 0.5) then
        ! Add an Angstrom thick layer at the bottom with 0 velocity to impose a
        ! no-slip BBC in the output, if anything but piecewise constant is used.
        nk_valid(i) = nk_valid(i) + 1 ; k2 = nk_valid(i)
        h_f(k2,i) = angstrom ; u_f(k2,i) = 0.0
        ! GM: D_pt is always slightly larger (by 1E-6 or so) than shelf_depth, so
        ! I consider that the ice shelf is grounded when
        ! shelf_depth(I) + 1.0E-3 > D_pt(i)
        if (ice_shelf .and. shelf_depth(i) + 1.0e-3 > d_pt(i)) nk_valid(i)=0
      endif ; enddo


      do i=isq,ieq ; if (nk_valid(i) > 0) then
      ! Calculate the z* interface heights for tracers.
        htot = 0.0 ; do k=1,nk_valid(i) ; htot = htot + h_f(k,i) ; enddo
        dilate = 0.0
        if (htot*gv%H_to_m > 2.0*angstrom) then
           dilate = max((d_pt(i) - shelf_depth(i)),angstrom)/htot
        endif
        e(nk_valid(i)+1) = -d_pt(i)
        do k=nk_valid(i),1,-1 ; e(k) = e(k+1) + h_f(k,i)*dilate ; enddo

      ! Interpolate each variable into depth space.
        k_bot = 1 ; k_bot_prev = -1
        do kz=1,cs%nk_zspace
          call find_overlap(e, cs%Z_int(kz), cs%Z_int(kz+1), nk_valid(i), &
                            k_bot, k_top, k_bot, wt, z1, z2)
          if (k_top>nk_valid(i)) exit

          !GM if top range that is being map is below the shelf, interpolate
          ! otherwise keep missing_vel
          if (cs%Z_int(kz)<=-shelf_depth(i)) then

            if (linear_velocity_profiles) then
              k = k_top
              if (k /= k_bot_prev) then
                ! Calculate the intra-cell profile.
                slope = 0.0 ! ; curv = 0.0
                if ((k < nk_valid(i)) .and. (k > nkml)) call &
                  find_limited_slope(u_f(:,i), e, slope, k)
              endif
              ! This is the piecewise linear form.
              cs%u_z(i,j,kz) = wt(k) * (u_f(k,i) + 0.5*slope*(z2(k) + z1(k)))
              ! For the piecewise parabolic form add the following...
              !     + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2))
              do k=k_top+1,k_bot-1
                cs%u_z(i,j,kz) = cs%u_z(i,j,kz) + wt(k)*u_f(k,i)
              enddo
              if (k_bot > k_top) then ; k = k_bot
                ! Calculate the intra-cell profile.
                slope = 0.0 ! ; curv = 0.0
                if ((k < nk_valid(i)) .and. (k > nkml)) call &
                  find_limited_slope(u_f(:,i), e, slope, k)
                 ! This is the piecewise linear form.
                cs%u_z(i,j,kz) = cs%u_z(i,j,kz) + wt(k) * &
                    (u_f(k,i) + 0.5*slope*(z2(k) + z1(k)))
                ! For the piecewise parabolic form add the following...
                !     + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2))
              endif
              k_bot_prev = k_bot
            else ! Use piecewise constant profiles.
              cs%u_z(i,j,kz) = wt(k_top)*u_f(k_top,i)
              do k=k_top+1,k_bot
                cs%u_z(i,j,kz) = cs%u_z(i,j,kz) + wt(k)*u_f(k,i)
              enddo
            endif ! linear profiles
          endif ! below shelf
        enddo ! kz-loop
      endif ; enddo ! I-loop and mask
    enddo ! j-loop

    call post_data(cs%id_u_z, cs%u_z, cs%diag)
  endif

  ! meridional velocity component
  if (cs%id_v_z > 0) then
    do kz=1,cs%nk_zspace ; do j=jsq,jeq ; do i=is,ie
      cs%v_z(i,j,kz) = cs%missing_vel
    enddo ; enddo ; enddo

    do j=jsq,jeq
      shelf_depth(:) = 0.0 ! initially all is open ocean
      ! Remove all massless layers.
      do i=is,ie
        nk_valid(i) = 0 ; d_pt(i) = 0.5*(g%bathyT(i,j)+g%bathyT(i,j+1))
        if (ice_shelf) then
          if (frac_shelf_h(i,j)+frac_shelf_h(i,j+1) > 0.) then ! under shelf
            shelf_depth(i) = abs(0.5*(ssh(i,j)+ssh(i,j+1)))
          endif
        endif
      enddo
      do k=1,nk ; do i=is,ie
        if ((g%mask2dCv(i,j) > 0.5) .and. (h(i,j,k)+h(i,j+1,k) > 4.0*angstrom)) then
          nk_valid(i) = nk_valid(i) + 1 ; k2 = nk_valid(i)
          h_f(k2,i) = 0.5*(h(i,j,k)+h(i,j+1,k)) ; v_f(k2,i) = v(i,j,k)
        endif
      enddo ; enddo
      do i=is,ie ; if (g%mask2dCv(i,j) > 0.5) then
        ! Add an Angstrom thick layer at the bottom with 0 velocity to impose a
        ! no-slip BBC in the output, if anything but piecewise constant is used.
        nk_valid(i) = nk_valid(i) + 1 ; k2 = nk_valid(i)
        h_f(k2,i) = angstrom ; v_f(k2,i) = 0.0
        if (ice_shelf .and. shelf_depth(i) + 1.0e-3 > d_pt(i)) nk_valid(i)=0
      endif ; enddo

      do i=is,ie ; if (nk_valid(i) > 0) then
      ! Calculate the z* interface heights for tracers.
        htot = 0.0 ; do k=1,nk_valid(i) ; htot = htot + h_f(k,i) ; enddo
        dilate = 0.0
        if (htot > 2.0*angstrom) then
           dilate = max((d_pt(i) - shelf_depth(i)),angstrom)/htot
        endif
        e(nk_valid(i)+1) = -d_pt(i)
        do k=nk_valid(i),1,-1 ; e(k) = e(k+1) + h_f(k,i)*dilate ; enddo

      ! Interpolate each variable into depth space.
        k_bot = 1 ; k_bot_prev = -1
        do kz=1,cs%nk_zspace
          call find_overlap(e, cs%Z_int(kz), cs%Z_int(kz+1), nk_valid(i), &
                            k_bot, k_top, k_bot, wt, z1, z2)
          if (k_top>nk_valid(i)) exit
          !GM if top range that is being map is below the shelf, interpolate
          ! otherwise keep missing_vel
          if (cs%Z_int(kz)<=-shelf_depth(i)) then
            if (linear_velocity_profiles) then
              k = k_top
              if (k /= k_bot_prev) then
                ! Calculate the intra-cell profile.
                slope = 0.0 ! ; curv = 0.0
                if ((k < nk_valid(i)) .and. (k > nkml)) call &
                  find_limited_slope(v_f(:,i), e, slope, k)
              endif
              ! This is the piecewise linear form.
              cs%v_z(i,j,kz) = wt(k) * (v_f(k,i) + 0.5*slope*(z2(k) + z1(k)))
              ! For the piecewise parabolic form add the following...
              !     + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2))
              do k=k_top+1,k_bot-1
                cs%v_z(i,j,kz) = cs%v_z(i,j,kz) + wt(k)*v_f(k,i)
              enddo
              if (k_bot > k_top) then ; k = k_bot
                ! Calculate the intra-cell profile.
                slope = 0.0 ! ; curv = 0.0
                if ((k < nk_valid(i)) .and. (k > nkml)) call &
                  find_limited_slope(v_f(:,i), e, slope, k)
                 ! This is the piecewise linear form.
                cs%v_z(i,j,kz) = cs%v_z(i,j,kz) + wt(k) * &
                    (v_f(k,i) + 0.5*slope*(z2(k) + z1(k)))
                ! For the piecewise parabolic form add the following...
                !     + C1_3*curv*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2))
              endif
              k_bot_prev = k_bot
            else ! Use piecewise constant profiles.
              cs%v_z(i,j,kz) = wt(k_top)*v_f(k_top,i)
              do k=k_top+1,k_bot
                cs%v_z(i,j,kz) = cs%v_z(i,j,kz) + wt(k)*v_f(k,i)
              enddo
            endif ! linear profiles
          endif ! below shelf
          enddo ! kz-loop
      endif ; enddo ! i-loop and mask
    enddo ! J-loop

    call post_data(cs%id_v_z, cs%v_z, cs%diag)
  endif

  ! tracer concentrations
  if (cs%num_tr_used > 0) then

    do m=1,cs%num_tr_used ; do kz=1,cs%nk_zspace ; do j=js,je ; do i=is,ie
      cs%tr_z(m)%p(i,j,kz) = cs%missing_tr(m)
    enddo ; enddo ; enddo ; enddo

    do j=js,je
      shelf_depth(:) = 0.0 ! initially all is open ocean
      ! Remove all massless layers.
      do i=is,ie
        nk_valid(i) = 0 ; d_pt(i) = g%bathyT(i,j)
        if (ice_shelf) then
          if (frac_shelf_h(i,j) > 0.) then ! under shelf
            shelf_depth(i) = abs(ssh(i,j))
          endif
        endif
      enddo
      do k=1,nk ; do i=is,ie
        if ((g%mask2dT(i,j) > 0.5) .and. (h(i,j,k) > 2.0*angstrom)) then
          nk_valid(i) = nk_valid(i) + 1 ; k2 = nk_valid(i)
          h_f(k2,i) = h(i,j,k)
          if (ice_shelf .and. shelf_depth(i) + 1.0e-3 > d_pt(i)) nk_valid(i)=0
          do m=1,cs%num_tr_used ; tr_f(k2,m,i) = cs%tr_model(m)%p(i,j,k) ; enddo
        endif
      enddo ; enddo

      do i=is,ie ; if (nk_valid(i) > 0) then
      ! Calculate the z* interface heights for tracers.
        htot = 0.0 ;  do k=1,nk_valid(i) ; htot = htot + h_f(k,i) ; enddo
        dilate = 0.0
        if (htot > 2.0*angstrom) then
           dilate = max((d_pt(i) - shelf_depth(i)),angstrom)/htot
        endif
        e(nk_valid(i)+1) = -d_pt(i)
        do k=nk_valid(i),1,-1 ; e(k) = e(k+1) + h_f(k,i)*dilate ; enddo

      ! Interpolate each variable into depth space.
        k_bot = 1 ; k_bot_prev = -1
        do kz=1,cs%nk_zspace
          call find_overlap(e, cs%Z_int(kz), cs%Z_int(kz+1), nk_valid(i), &
                            k_bot, k_top, k_bot, wt, z1, z2)
          if (k_top>nk_valid(i)) exit
          if (cs%Z_int(kz)<=-shelf_depth(i)) then
            do m=1,cs%num_tr_used
              k = k_top
              if (k /= k_bot_prev) then
                ! Calculate the intra-cell profile.
                sl_tr(m) = 0.0 ! ; cur_tr(m) = 0.0
                if ((k < nk_valid(i)) .and. (k > nkml)) call &
                  find_limited_slope(tr_f(:,m,i), e, sl_tr(m), k)
              endif
              ! This is the piecewise linear form.
              cs%tr_z(m)%p(i,j,kz) = wt(k) * &
                  (tr_f(k,m,i) + 0.5*sl_tr(m)*(z2(k) + z1(k)))
              ! For the piecewise parabolic form add the following...
              !     + C1_3*cur_tr(m)*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2))
              do k=k_top+1,k_bot-1
                cs%tr_z(m)%p(i,j,kz) = cs%tr_z(m)%p(i,j,kz) + wt(k)*tr_f(k,m,i)
              enddo
              if (k_bot > k_top) then
                k = k_bot
                ! Calculate the intra-cell profile.
                sl_tr(m) = 0.0 ! ; cur_tr(m) = 0.0
                if ((k < nk_valid(i)) .and. (k > nkml)) call &
                  find_limited_slope(tr_f(:,m,i), e, sl_tr(m), k)
                ! This is the piecewise linear form.
                cs%tr_z(m)%p(i,j,kz) = cs%tr_z(m)%p(i,j,kz) + wt(k) * &
                    (tr_f(k,m,i) + 0.5*sl_tr(m)*(z2(k) + z1(k)))
                ! For the piecewise parabolic form add the following...
                !     + C1_3*cur_tr(m)*(z2(k)**2 + z2(k)*z1(k) + z1(k)**2))
              endif
            enddo
            k_bot_prev = k_bot
          endif ! below shelf
        enddo ! kz-loop
      endif ; enddo ! i-loop and mask

    enddo ! j-loop

    do m=1,cs%num_tr_used
      if (cs%id_tr(m) > 0) call post_data(cs%id_tr(m), cs%tr_z(m)%p, cs%diag)
      if (cs%id_tr_xyave(m) > 0) then
        layer_ave = global_z_mean(cs%tr_z(m)%p,g,cs,m)
        call post_data_1d_k(cs%id_tr_xyave(m), layer_ave, cs%diag)
      endif
    enddo
  endif

end subroutine calculate_z_diag_fields

!> This subroutine maps horizontal transport into depth space for diagnostic output.
subroutine calculate_z_transport(uh_int, vh_int, h, dt, G, GV, CS)
  type(ocean_grid_type),                     intent(inout) :: G    !< The ocean's grid structure.
  type(verticalgrid_type),                   intent(in)    :: GV   !< The ocean's vertical grid
                                                                   !! structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(in)    :: uh_int !< Time integrated zonal
                                                                   !! transport (m3 or kg).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(in)    :: vh_int !< Time integrated meridional
                                                                   !! transport (m3 or kg).
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(in)    :: h    !< Layer thicknesses, in H
                                                                   !! (usually m or kg m-2).
  real,                                      intent(in)    :: dt   !< The time difference in s since
                                                                   !! the last call to this
                                                                   !! subroutine.
  type(diag_to_z_cs),                        pointer       :: CS   !< Control structure returned by
                                                                   !! previous call to
                                                                   !! diagnostics_init.

!   This subroutine maps horizontal transport into depth space for diagnostic output.

! Arguments:
!  (in)      uh_int - time integrated zonal transport (m3 or kg)
!  (in)      vh_int - time integrated meridional transport (m3 or kg)
!  (in)      h      - layer thickness (meter or kg/m2)
!  (in)      dt     - time difference (sec) since last call to this routine
!  (in)      G      - ocean grid structure
!  (in)      GV     - The ocean's vertical grid structure
!  (in)      CS     - control structure returned by previous call to diagnostics_init

  real, dimension(SZI_(G), SZJ_(G)) :: &
    htot, &        ! total layer thickness (meter or kg/m2)
    dilate         ! nondimensional factor by which to dilate layers to
                   ! convert them into z* space.  (-G%D < z* < 0)

  real, dimension(SZI_(G), max(CS%nk_zspace,1)) :: &
    uh_Z           ! uh_int interpolated into depth space (m3 or kg)
  real, dimension(SZIB_(G), max(CS%nk_zspace,1)) :: &
    vh_Z           ! vh_int interpolated into depth space (m3 or kg)

  real :: h_rem    ! dilated thickness of a layer that has yet to be mapped
                   ! into depth space (meter or kg/m2)
  real :: uh_rem   ! integrated zonal transport of a layer that has yet to be
                   ! mapped into depth space (m3 or kg)
  real :: vh_rem   ! integrated meridional transport of a layer that has yet
                   ! to be mapped into depth space (m3 or kg)
  real :: h_here   ! thickness of a layer that is within the range of the
                   ! current depth level (meter or kg/m2)
  real :: h_above  ! thickness of a layer that is above the current depth
                   ! level (meter or kg.m2)
  real :: uh_here  ! zonal transport of a layer that is attributed to the
                   ! current depth level (m3 or kg)
  real :: vh_here  ! meridional transport of a layer that is attributed to
                   ! the current depth level (m3 or kg)
  real :: Idt      ! inverse of the time step (sec)

  real :: Z_int_above(szib_(g)) ! height of the interface atop a layer (meter or kg/m2)

  integer :: kz(szib_(g)) ! index of depth level that is being contributed to

  integer :: i, j, k, is, ie, js, je, Isq, Ieq, Jsq, Jeq, nk, nk_z
  is  = g%isc  ; ie  = g%iec  ; js  = g%jsc  ; je  = g%jec ; nk = g%ke
  isq = g%IscB ; ieq = g%IecB ; jsq = g%JscB ; jeq = g%JecB

  if (.not.associated(cs)) call mom_error(fatal, &
         "calculate_Z_transport: Module must be initialized before it is used.")
  if ((cs%id_uh_Z <= 0) .and. (cs%id_vh_Z <= 0)) return

  idt = 1.0 ; if (dt > 0.0) idt = 1.0 / dt
  nk_z = cs%nk_zspace

  if (nk_z <= 0) return

  ! Determine how much the layers will be dilated in recasting them into z*
  ! coordiantes.  (-G%D < z* < 0).
  do j=jsq,jeq+1 ; do i=isq,ieq+1
    htot(i,j) = gv%H_subroundoff
  enddo ; enddo
  do k=1,nk ; do j=jsq,jeq+1 ; do i=isq,ieq+1
    htot(i,j) = htot(i,j) + h(i,j,k)
  enddo ; enddo ; enddo
  do j=jsq,jeq+1 ; do i=isq,ieq+1
    dilate(i,j) = g%bathyT(i,j) / htot(i,j)
  enddo ; enddo

  ! zonal transport
  if (cs%id_uh_Z > 0) then ; do j=js,je
    do i=isq,ieq
      kz(i) = nk_z ; z_int_above(i) = -0.5*(g%bathyT(i,j)+g%bathyT(i+1,j))
    enddo
    do k=nk_z,1,-1 ; do i=isq,ieq
      uh_z(i,k) = 0.0
      if (cs%Z_int(k) < z_int_above(i)) kz(i) = k-1
    enddo ; enddo
    do k=nk,1,-1 ; do i=isq,ieq
      h_rem = 0.5*(dilate(i,j)*h(i,j,k) + dilate(i+1,j)*h(i+1,j,k))
      uh_rem = uh_int(i,j,k)
      z_int_above(i) = z_int_above(i) + h_rem

      do ! Distribute this layer's transport into the depth-levels.
        h_above = z_int_above(i) - cs%Z_int(kz(i))
        if ((kz(i) == 1) .or. (h_above <= 0.0) .or. (h_rem <= 0.0)) then
          ! The entire remaining transport is on this level.
          uh_z(i,kz(i)) = uh_z(i,kz(i)) + uh_rem ; exit
        else
          h_here = h_rem - h_above
          uh_here = uh_rem * (h_here / h_rem)

          h_rem = h_rem - h_here ! = h_above
          uh_z(i,kz(i)) = uh_z(i,kz(i)) + uh_here
          uh_rem = uh_rem - uh_here
          kz(i) = kz(i) - 1
        endif
      enddo ! End of loop through the target depth-space levels.
    enddo ; enddo
    do k=1,nk_z ; do i=isq,ieq
      cs%uh_z(i,j,k) = uh_z(i,k)*idt
    enddo ; enddo
  enddo ; endif

  ! meridional transport
  if (cs%id_vh_Z > 0) then ; do j=jsq,jeq
    do i=is,ie
      kz(i) = nk_z ; z_int_above(i) = -0.5*(g%bathyT(i,j)+g%bathyT(i,j+1))
    enddo
    do k=nk_z,1,-1 ; do i=is,ie
      vh_z(i,k) = 0.0
      if (cs%Z_int(k) < z_int_above(i)) kz(i) = k-1
    enddo ; enddo
    do k=nk,1,-1 ; do i=is,ie
      h_rem = 0.5*(dilate(i,j)*h(i,j,k) + dilate(i,j+1)*h(i,j+1,k))
      vh_rem = vh_int(i,j,k)
      z_int_above(i) = z_int_above(i) + h_rem

      do ! Distribute this layer's transport into the depth-levels.
        h_above = z_int_above(i) - cs%Z_int(kz(i))
        if ((kz(i) == 1) .or. (h_above <= 0.0) .or. (h_rem <= 0.0)) then
          ! The entire remaining transport is on this level.
          vh_z(i,kz(i)) = vh_z(i,kz(i)) + vh_rem ; exit
        else
          h_here = h_rem - h_above
          vh_here = vh_rem * (h_here / h_rem)

          h_rem = h_rem - h_here ! = h_above
          vh_z(i,kz(i)) = vh_z(i,kz(i)) + vh_here
          vh_rem = vh_rem - vh_here
          kz(i) = kz(i) - 1
        endif
      enddo ! End of loop through the target depth-space levels.
    enddo ; enddo
    do k=1,nk_z ; do i=is,ie
      cs%vh_z(i,j,k) = vh_z(i,k)*idt
    enddo ; enddo
  enddo ; endif

  if (cs%id_uh_Z > 0) then
    do k=1,nk_z ; do j=js,je ; do i=isq,ieq
      cs%uh_z(i,j,k) = cs%uh_z(i,j,k)*gv%H_to_kg_m2
    enddo ; enddo ; enddo
    call post_data(cs%id_uh_Z, cs%uh_z, cs%diag)
  endif

  if (cs%id_vh_Z > 0) then
    do k=1,nk_z ; do j=jsq,jeq ; do i=is,ie
      cs%vh_z(i,j,k) = cs%vh_z(i,j,k)*gv%H_to_kg_m2
    enddo ; enddo ; enddo
    call post_data(cs%id_vh_Z, cs%vh_z, cs%diag)
  endif

end subroutine calculate_z_transport

!>   This subroutine determines the layers bounded by interfaces e that overlap
!! with the depth range between Z_top and Z_bot, and the fractional weights
!! of each layer. It also calculates the normalized relative depths of the range
!! of each layer that overlaps that depth range.
subroutine find_overlap(e, Z_top, Z_bot, k_max, k_start, k_top, k_bot, wt, z1, z2)
  real, dimension(:), intent(in)    :: e      !< Column interface heights (meter or kg/m2).
  real,               intent(in)    :: Z_top  !< Top of range being mapped to (meter or kg/m2).
  real,               intent(in)    :: Z_bot  !< Bottom of range being mapped to (meter or kg/m2).
  integer,            intent(in)    :: k_max  !< Number of valid layers.
  integer,            intent(in)    :: k_start !< Layer at which to start searching.
  integer,            intent(inout) :: k_top  !< Indices of top layers that overlap with the depth
                                              !! range.
  integer,            intent(inout) :: k_bot  !< Indices of bottom layers that overlap with the
                                              !! depth range.
  real, dimension(:), intent(out)   :: wt     !< Relative weights of each layer from k_top to k_bot.
  real, dimension(:), intent(out)   :: z1, z2 !< Depths of the top and bottom limits of the part of
       !! a layer that contributes to a depth level, relative to the cell center and normalized
       !! by the cell thickness (nondim).  Note that -1/2 <= z1 < z2 <= 1/2.

!   This subroutine determines the layers bounded by interfaces e that overlap
! with the depth range between Z_top and Z_bot, and the fractional weights
! of each layer. It also calculates the normalized relative depths of the range
! of each layer that overlaps that depth range.

!   Note that by convention, e decreases with increasing k and Z_top > Z_bot.
!
! Arguments:
!  (in)          e        - column interface heights (meter or kg/m2)
!  (in)      Z_top        - top of range being mapped to (meter or kg/m2)
!  (in)      Z_bot        - bottom of range being mapped to (meter or kg/m2)
!  (in)      k_max        - number of valid layers
!  (in)      k_start      - layer at which to start searching
!  (out)     k_top, k_bot - indices of top and bottom layers that
!                           overlap with the depth range
!  (out)     wt           - relative weights of each layer from k_top to k_bot
!  (out)     z1, z2       - depths of the top and bottom limits of
!                           the part of a layer that contributes to a depth level,
!                           relative to the cell center and normalized by the cell
!                           thickness (nondim).  Note that -1/2 <= z1 < z2 <= 1/2.

  real    :: Ih, e_c, tot_wt, I_totwt
  integer :: k

  do k=k_start,k_max ; if (e(k+1)<z_top) exit ; enddo
  k_top = k
  if (k>k_max) return

  ! Determine the fractional weights of each layer.
  ! Note that by convention, e and Z_int decrease with increasing k.
  if (e(k+1)<=z_bot) then
    wt(k) = 1.0 ; k_bot = k
    ih = 1.0 / (e(k)-e(k+1))
    e_c = 0.5*(e(k)+e(k+1))
    z1(k) = (e_c - min(e(k),z_top)) * ih
    z2(k) = (e_c - z_bot) * ih
  else
    wt(k) = min(e(k),z_top) - e(k+1) ; tot_wt = wt(k) ! These are always > 0.
    z1(k) = (0.5*(e(k)+e(k+1)) - min(e(k),z_top)) / (e(k)-e(k+1))
    z2(k) = 0.5
    k_bot = k_max
    do k=k_top+1,k_max
      if (e(k+1)<=z_bot) then
        k_bot = k
        wt(k) = e(k) - z_bot ; z1(k) = -0.5
        z2(k) = (0.5*(e(k)+e(k+1)) - z_bot) / (e(k)-e(k+1))
      else
        wt(k) = e(k) - e(k+1) ; z1(k) = -0.5 ; z2(k) = 0.5
      endif
      tot_wt = tot_wt + wt(k) ! wt(k) is always > 0.
      if (k>=k_bot) exit
    enddo

    i_totwt = 1.0 / tot_wt
    do k=k_top,k_bot ; wt(k) = i_totwt*wt(k) ; enddo
  endif

end subroutine find_overlap

!>   This subroutine determines a limited slope for val to be advected with
!! a piecewise limited scheme.
subroutine find_limited_slope(val, e, slope, k)
  real, dimension(:), intent(in)  :: val !< A column of values that are being interpolated.
  real, dimension(:), intent(in)  :: e   !< Column interface heights (meter or kg/m2).
  real,               intent(out) :: slope !< Normalized slope in the intracell distribution of val.
  integer,            intent(in)  :: k   !< Layer whose slope is being determined.

!   This subroutine determines a limited slope for val to be advected with
! a piecewise limited scheme.

! Arguments:
!  (in)      val   - a column of values that are being interpolated
!  (in)      e     - column interface heights (meter or kg/m2)
!  (in)      slope - normalized slope in the intracell distribution of val
!  (in)      k     - layer whose slope is being determined

  real :: d1, d2

  if ((val(k)-val(k-1)) * (val(k)-val(k+1)) >= 0.0) then
    slope = 0.0 ! ; curvature = 0.0
  else
    d1 = 0.5*(e(k-1)-e(k+1)) ; d2 = 0.5*(e(k)-e(k+2))
    slope = (d1**2*(val(k+1) - val(k)) + d2**2*(val(k) - val(k-1))) * &
            ((e(k) - e(k+1)) / (d1*d2*(d1+d2)))
    ! slope = 0.5*(val(k+1) - val(k-1))
    ! This is S.J. Lin's form of the PLM limiter.
    slope = sign(1.0,slope) * min(abs(slope), &
        2.0*(max(val(k-1),val(k),val(k+1)) - val(k)), &
        2.0*(val(k) - min(val(k-1),val(k),val(k+1))))
    ! curvature = 0.0
  endif

end subroutine find_limited_slope

! #@# This subroutine needs a doxygen description
subroutine calc_zint_diags(h, in_ptrs, ids, num_diags, G, GV, CS)
  type(ocean_grid_type),                    intent(in) :: G    !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h    !< Layer thicknesses, in H
                                                               !! (usually m or kg m-2).
  type(p3d), dimension(:),                  intent(in) :: in_ptrs
  integer,   dimension(:),                  intent(in) :: ids
  integer,                                  intent(in) :: num_diags
  type(verticalgrid_type),                  intent(in) :: GV   !< The ocean's vertical grid
                                                               !! structure.
  type(diag_to_z_cs),                       pointer    :: CS

  real, dimension(SZI_(G),SZJ_(G),max(CS%nk_zspace+1,1),max(num_diags,1)) :: &
    diag_on_Z  ! diagnostics interpolated to depth space
  real, dimension(SZI_(G),SZK_(G)+1) :: e
  real, dimension(max(num_diags,1),SZI_(G),SZK_(G)+1) :: diag2d

  real, dimension(SZI_(G)) :: &
    htot, &              ! summed layer thicknesses (meter or kg/m2)
    dilate               ! proportion by which to dilate every layer
  real :: wt             ! weighting of the interface above in the
                         ! interpolation to target depths
  integer :: kL(szi_(g)) ! layer-space index of shallowest interface
                         ! below the target depth

  integer :: i, j, k, k2, kz, is, ie, js, je, nk, m
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nk = g%ke

  if (num_diags < 1) return
  if (.not.associated(cs)) call mom_error(fatal, &
       "calc_Zint_diags: Module must be initialized before it is used.")

  do j=js,je
    ! Calculate the stretched z* interface depths.
    do i=is,ie ; htot(i) = 0.0 ; kl(i) = 1 ; enddo
    do k=1,nk ; do i=is,ie ; htot(i) = htot(i) + h(i,j,k) ; enddo ; enddo
    do i=is,ie
      dilate(i) = 0.0
      if (htot(i)*gv%H_to_m > 0.5) dilate(i) = (g%bathyT(i,j) - 0.0) / htot(i)
      e(i,nk+1) = -g%bathyT(i,j)
    enddo
    do k=nk,1,-1 ; do i=is,ie
      e(i,k) = e(i,k+1) + h(i,j,k) * dilate(i)
    enddo ; enddo
    ! e(i,1) should be 0 as a consistency check.

    do k=1,nk+1 ; do i=is,ie ; do m=1,num_diags
      diag2d(m,i,k) = in_ptrs(m)%p(i,j,k)
    enddo ; enddo ; enddo

    do kz=1,cs%nk_zspace+1 ; do i=is,ie
      ! Find the interface below the target Z-file depth, kL.
      if (cs%Z_int(kz) < e(i,nk+1)) then
        kl(i) = nk+2
      else
        do k=kl(i),nk+1 ; if (cs%Z_int(kz) > e(i,k)) exit ; enddo
        kl(i) = k
      endif
      if (kl(i)>1) then
        if (cs%Z_int(kz) > e(i,kl(i)-1)) call mom_error(fatal, &
        "calc_Zint_diags: Interface depth mapping is incorrect.")
      endif
      if ((kl(i)>1) .and. (kl(i)<=nk+1)) then
        if (e(i,kl(i)-1) == e(i,kl(i))) call mom_error(warning, &
          "calc_Zint_diags: Interface depths equal.", all_print=.true.)
        if (e(i,kl(i)-1) - e(i,kl(i)) < 0.0) call mom_error(fatal, &
          "calc_Zint_diags: Interface depths inverted.")
      endif

      if (kl(i) <= 1) then
        do m=1,num_diags
          diag_on_z(i,j,kz,m) = diag2d(m,i,1)
        enddo
      elseif (kl(i) > nk+1) then
        do m=1,num_diags
          diag_on_z(i,j,kz,m) = cs%missing_value
        enddo
      else
        wt = 0.0 ! This probably should not happen?
        if (e(i,kl(i)-1) - e(i,kl(i)) > 0.0) &
          wt = (cs%Z_int(kz) - e(i,kl(i))) / (e(i,kl(i)-1) - e(i,kl(i)))
        if ((wt < 0.0) .or. (wt > 1.0)) call mom_error(fatal, &
          "calc_Zint_diags: Bad value of wt found.")
        do m=1,num_diags
          diag_on_z(i,j,kz,m) = wt * diag2d(m,i,kl(i)-1) + &
                                (1.0-wt) * diag2d(m,i,kl(i))
        enddo
      endif
    enddo ; enddo

  enddo

  do m=1,num_diags
    if (ids(m) > 0) call post_data(ids(m), diag_on_z(:,:,:,m), cs%diag)
  enddo

end subroutine calc_zint_diags

!> This subroutine registers a tracer to be output in depth space.
subroutine register_z_tracer(tr_ptr, name, long_name, units, Time, G, CS, standard_name,   &
     cmor_field_name, cmor_long_name, cmor_units, cmor_standard_name)
  type(ocean_grid_type),            intent(in) :: G      !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                    target,         intent(in) :: tr_ptr !< Tracer for translation to Z-space.
  character(len=*),                 intent(in) :: name   !< name for the output tracer.
  character(len=*),                 intent(in) :: long_name !< Long name for the output tracer.
  character(len=*),                 intent(in) :: units  !< Units of output tracer.
  character(len=*), optional,       intent(in) :: standard_name
  type(time_type),                  intent(in) :: Time   !< Current model time.
  type(diag_to_z_cs),               pointer    :: CS     !< Control struct returned by previous
                                                         !! call to diagnostics_init.
  character(len=*), optional,       intent(in) :: cmor_field_name !< cmor name of a field.
  character(len=*), optional,       intent(in) :: cmor_long_name  !< cmor long name of a field.
  character(len=*), optional,       intent(in) :: cmor_units      !< cmor units of a field.
  character(len=*), optional,       intent(in) :: cmor_standard_name !< cmor standardized name
                                                                  !! associated with a field.

!   This subroutine registers a tracer to be output in depth space.
! Arguments:
!  (in)      tr_ptr    - tracer for translation to Z-space
!  (in)      name      - name for the output tracer
!  (in)      long_name - long name for the output tracer
!  (in)      units     - units of output tracer
!  (in)      Time      - current model time
!  (in)      G         - ocean grid structure
!  (in)      CS        - control struct returned by previous call to diagnostics_init
!  (in,opt)  cmor_field_name    - cmor name of a field
!  (in,opt)  cmor_long_name     - cmor long name of a field
!  (in,opt)  cmor_units         - cmor units of a field
!  (in,opt)  cmor_standard_name - cmor standardized name associated with a field

  character(len=256) :: posted_standard_name
  character(len=256) :: posted_cmor_units
  character(len=256) :: posted_cmor_standard_name
  character(len=256) :: posted_cmor_long_name

  if (cs%nk_zspace<1) return

  if (present(standard_name)) then
    posted_standard_name = standard_name
  else
    posted_standard_name = 'not provided'
  endif

  call register_z_tracer_low(tr_ptr, name, long_name, units, trim(posted_standard_name), time, g, cs)

  if (present(cmor_field_name)) then
    ! Fallback values for strings set to "NULL"
    posted_cmor_units         = "not provided"   !
    posted_cmor_standard_name = "not provided"   ! values might be replaced with a CS%missing field?
    posted_cmor_long_name     = "not provided"   !

    ! If attributes are present for MOM variable names, use them first for the register_diag_field
    ! call for CMOR verison of the variable
    posted_cmor_units         = units
    posted_cmor_long_name     = long_name
    posted_cmor_standard_name = posted_standard_name

    ! If specified in the call to register_diag_field, override attributes with the CMOR versions
    if (present(cmor_units)) posted_cmor_units                 = cmor_units
    if (present(cmor_standard_name)) posted_cmor_standard_name = cmor_standard_name
    if (present(cmor_long_name)) posted_cmor_long_name         = cmor_long_name

    call register_z_tracer_low(tr_ptr, trim(cmor_field_name), trim(posted_cmor_long_name),&
         trim(posted_cmor_units), trim(posted_cmor_standard_name), time, g, cs)

  endif

end subroutine register_z_tracer

!> This subroutine registers a tracer to be output in depth space.
subroutine register_z_tracer_low(tr_ptr, name, long_name, units, standard_name, Time, G, CS)
  type(ocean_grid_type),      intent(in) :: G      !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                      target, intent(in) :: tr_ptr !< Tracer for translation to Z-space.
  character(len=*),           intent(in) :: name   !< Name for the output tracer.
  character(len=*),           intent(in) :: long_name !< Long name for output tracer.
  character(len=*),           intent(in) :: units  !< Units of output tracer.
  character(len=*),           intent(in) :: standard_name
  type(time_type),            intent(in) :: Time   !< Current model time.
  type(diag_to_z_cs),         pointer    :: CS     !< Control struct returned by previous call to
                                                   !! diagnostics_init.

!   This subroutine registers a tracer to be output in depth space.

! Arguments:
!  (in)      tr_ptr    - tracer for translation to Z-space
!  (in)      name      - name for the output tracer
!  (in)      long_name - long name for output tracer
!  (in)      units     - units of output tracer
!  (in)      Time      - current model time
!  (in)      G         - ocean grid structure
!  (in)      CS        - control struct returned by previous call to diagnostics_init

  character(len=256) :: posted_standard_name
  integer :: isd, ied, jsd, jed, nk, m, id_test
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nk = g%ke

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

  if (cs%num_tr_used >= max_fields_) then
    call mom_error(warning,"MOM_diag_to_Z:  Attempted to register and use "//&
                   "more than MAX_FIELDS_ z-space tracers via register_Z_tracer.")
    return
  endif

  m = cs%num_tr_used + 1

  cs%missing_tr(m) = cs%missing_value ! This could be changed later, if desired.
  if (cs%nk_zspace > 0) then
    cs%id_tr(m) = register_diag_field('ocean_model_zold', name, cs%axesTz, time,           &
                                      long_name, units, missing_value=cs%missing_tr(m), &
                                      standard_name=standard_name)
    cs%id_tr_xyave(m) = register_diag_field('ocean_model_zold', name//'_xyave', cs%axesZ, time,           &
                                      long_name, units, missing_value=cs%missing_tr(m), &
                                      standard_name=standard_name)
  else
    id_test = register_diag_field('ocean_model_zold', name, cs%diag%axesT1, time,      &
                                  long_name, units, missing_value=cs%missing_tr(m), &
                                  standard_name=standard_name)
    if (id_test>0) call mom_error(warning, &
        "MOM_diag_to_Z_init: "//trim(name)// &
        " cannot be output without an appropriate depth-space target file.")
  endif

  if (cs%id_tr(m) <= 0) cs%id_tr(m) = -1
  if (cs%id_tr_xyave(m) <= 0) cs%id_tr_xyave(m) = -1
  if (cs%id_tr(m) > 0 .or. cs%id_tr_xyave(m) > 0) then
    cs%num_tr_used = m
    call safe_alloc_ptr(cs%tr_z(m)%p,isd,ied,jsd,jed,cs%nk_zspace)
    cs%tr_model(m)%p => tr_ptr
  endif

end subroutine register_z_tracer_low

! #@# This subroutine needs a doxygen comment.
subroutine mom_diag_to_z_init(Time, G, GV, param_file, diag, CS)
  type(time_type),         intent(in)    :: Time !< Current model time.
  type(ocean_grid_type),   intent(in)    :: G    !< The ocean's grid structure.
  type(verticalgrid_type), intent(in)    :: GV   !< The ocean's vertical grid structure.
  type(param_file_type),   intent(in)    :: param_file !< A structure to parse for run-time
                                                 !! parameters.
  type(diag_ctrl), target, intent(inout) :: diag !< Struct to regulate diagnostic output.
  type(diag_to_z_cs),      pointer       :: CS   !< Pointer to point to control structure for
                                                 !! this module.

! Arguments:
!  (in)      Time       - current model time
!  (in)      G          - ocean grid structure
!  (in)      GV - The ocean's vertical grid structure.
!  (in)      param_file - struct indicating open file to parse for model param values
!  (in)      diag       - struct to regulate diagnostic output
!  (in/out)  CS         - pointer to point to control structure for this module

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

  character(len=40)  :: mdl = "MOM_diag_to_Z" ! module name
  character(len=200) :: in_dir, zgrid_file    ! strings for directory/file
  character(len=48)  :: flux_units, string
  integer :: z_axis, zint_axis
  integer :: isd, ied, jsd, jed, IsdB, IedB, JsdB, JedB, nk, id_test
  isd  = g%isd   ; ied = g%ied  ; jsd  = g%jsd  ; jed  = g%jed ; nk = g%ke
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

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

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

  cs%diag => diag

  ! Read parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  ! Read in z-space info from a NetCDF file.
  call get_param(param_file, mdl, "Z_OUTPUT_GRID_FILE", zgrid_file, &
                 "The file that specifies the vertical grid for \n"//&
                 "depth-space diagnostics, or blank to disable \n"//&
                 "depth-space output.", default="")

  if (len_trim(zgrid_file) > 0) then
    call get_param(param_file, mdl, "INPUTDIR", in_dir, &
                 "The directory in which input files are found.", default=".")
    in_dir = slasher(in_dir)
    call get_z_depths(trim(in_dir)//trim(zgrid_file), "zw", cs%Z_int, "zt", &
                      z_axis, zint_axis, cs%nk_zspace)
    call log_param(param_file, mdl, "!INPUTDIR/Z_OUTPUT_GRID_FILE", &
                   trim(in_dir)//trim(zgrid_file))
    call log_param(param_file, mdl, "!NK_ZSPACE (from file)", cs%nk_zspace, &
                 "The number of depth-space levels.  This is determined \n"//&
                 "from the size of the variable zw in the output grid file.", &
                 units="nondim")
  else
    cs%nk_zspace = -1
  endif

  if (cs%nk_zspace > 0) then

    call define_axes_group(diag, (/ diag%axesB1%handles(1),  diag%axesB1%handles(2),  z_axis /),    cs%axesBz)
    call define_axes_group(diag, (/ diag%axesT1%handles(1),  diag%axesT1%handles(2),  z_axis /),    cs%axesTz)
    call define_axes_group(diag, (/ diag%axesCu1%handles(1), diag%axesCu1%handles(2), z_axis /),    cs%axesCuz)
    call define_axes_group(diag, (/ diag%axesCv1%handles(1), diag%axesCv1%handles(2), z_axis /),    cs%axesCvz)
    call define_axes_group(diag, (/ diag%axesB1%handles(1),  diag%axesB1%handles(2),  zint_axis /), cs%axesBzi)
    call define_axes_group(diag, (/ diag%axesT1%handles(1),  diag%axesT1%handles(2),  zint_axis /), cs%axesTzi)
    call define_axes_group(diag, (/ diag%axesCu1%handles(1), diag%axesCu1%handles(2), zint_axis /), cs%axesCuzi)
    call define_axes_group(diag, (/ diag%axesCv1%handles(1), diag%axesCv1%handles(2), zint_axis /), cs%axesCvzi)
    call define_axes_group(diag, (/ z_axis /),    cs%axesZ)

    cs%id_u_z = register_diag_field('ocean_model_zold', 'u', cs%axesCuz, time,    &
        'Zonal Velocity in Depth Space', 'meter second-1',                     &
        missing_value=cs%missing_vel, cmor_field_name='uo', cmor_units='m s-1',&
        cmor_standard_name='sea_water_x_velocity', cmor_long_name='Sea Water X Velocity')
    if (cs%id_u_z>0) call safe_alloc_ptr(cs%u_z,isdb,iedb,jsd,jed,cs%nk_zspace)

    cs%id_v_z = register_diag_field('ocean_model_zold', 'v', cs%axesCvz, time,    &
        'Meridional Velocity in Depth Space', 'meter second-1',                &
        missing_value=cs%missing_vel, cmor_field_name='vo', cmor_units='m s-1',&
        cmor_standard_name='sea_water_y_velocity', cmor_long_name='Sea Water Y Velocity')
    if (cs%id_v_z>0) call safe_alloc_ptr(cs%v_z,isd,ied,jsdb,jedb,cs%nk_zspace)

    cs%id_uh_z = register_diag_field('ocean_model_zold', 'uh', cs%axesCuz, time,    &
        'Zonal Mass Transport (including SGS param) in Depth Space', flux_units, &
        missing_value=cs%missing_trans)
    if (cs%id_uh_z>0) call safe_alloc_ptr(cs%uh_z,isdb,iedb,jsd,jed,cs%nk_zspace)

    cs%id_vh_z = register_diag_field('ocean_model_zold', 'vh', cs%axesCvz, time,        &
        'Meridional Mass Transport (including SGS param) in Depth Space', flux_units,&
        missing_value=cs%missing_trans)
    if (cs%id_vh_z>0) call safe_alloc_ptr(cs%vh_z,isd,ied,jsdb,jedb,cs%nk_zspace)

  endif

end subroutine mom_diag_to_z_init

!> This subroutine reads the depths of the interfaces bounding the intended
!! layers from a NetCDF file.  If no appropriate file is found, -1 is returned
!! as the number of layers in the output file.  Also, a diag_manager axis is set
!! up with the same information as this axis.
subroutine get_z_depths(depth_file, int_depth_name, int_depth, cell_depth_name, &
                        z_axis_index, edge_index, nk_out)
  character(len=*),   intent(in)  :: depth_file
  character(len=*),   intent(in)  :: int_depth_name
  real, dimension(:), pointer     :: int_depth
  character(len=*),   intent(in)  :: cell_depth_name
  integer,            intent(out) :: z_axis_index
  integer,            intent(out) :: edge_index
  integer,            intent(out) :: nk_out

!   This subroutine reads the depths of the interfaces bounding the intended
! layers from a NetCDF file.  If no appropriate file is found, -1 is returned
! as the number of layers in the output file.  Also, a diag_manager axis is set
! up with the same information as this axis.

  real, allocatable   :: cell_depth(:)
  character (len=200) :: units, long_name
  integer :: ncid, status, intid, intvid, layid, layvid, k, ni

  nk_out = -1

  status = nf90_open(depth_file, nf90_nowrite, ncid);
  if (status /= nf90_noerr) then
    call mom_error(warning,"MOM_diag_to_Z get_Z_depths: "//&
        " Difficulties opening "//trim(depth_file)//" - "//&
        trim(nf90_strerror(status)))
    nk_out = -1 ; return
  endif

  status = nf90_inq_dimid(ncid, int_depth_name, intid)
  if (status /= nf90_noerr) then
    call mom_error(warning,"MOM_diag_to_Z get_Z_depths: "//&
      trim(nf90_strerror(status))//" Getting ID of dimension "//&
      trim(int_depth_name)//" in "//trim(depth_file))
    nk_out = -1 ; return
  endif

  status = nf90_inquire_dimension(ncid, intid, len=ni)
  if (status /= nf90_noerr) then
    call mom_error(warning,"MOM_diag_to_Z get_Z_depths: "//&
      trim(nf90_strerror(status))//" Getting number of interfaces of "//&
      trim(int_depth_name)//" in "//trim(depth_file))
    nk_out = -1 ; return
  endif

  if (ni < 2) then
    call mom_error(warning,"MOM_diag_to_Z get_Z_depths: "//&
        "At least two interface depths must be specified in "//trim(depth_file))
    nk_out = -1 ; return
  endif

  status = nf90_inq_dimid(ncid, cell_depth_name, layid)
  if (status /= nf90_noerr) call mom_error(warning,"MOM_diag_to_Z get_Z_depths: "//&
      trim(nf90_strerror(status))//" Getting ID of dimension "//&
      trim(cell_depth_name)//" in "//trim(depth_file))

  status = nf90_inquire_dimension(ncid, layid, len=nk_out)
  if (status /= nf90_noerr) call mom_error(warning,"MOM_diag_to_Z get_Z_depths: "//&
      trim(nf90_strerror(status))//" Getting number of interfaces of "//&
      trim(cell_depth_name)//" in "//trim(depth_file))

  if (ni /= nk_out+1) then
    call mom_error(warning,"MOM_diag_to_Z get_Z_depths: "//&
        "The interface depths must have one more point than cell centers in "//&
        trim(depth_file))
    nk_out = -1 ; return
  endif

  allocate(int_depth(nk_out+1))
  allocate(cell_depth(nk_out))

  status = nf90_inq_varid(ncid, int_depth_name, intvid)
  if (status /= nf90_noerr) call mom_error(fatal,"MOM_diag_to_Z get_Z_depths: "//&
      trim(nf90_strerror(status))//" Getting ID of variable "//&
      trim(int_depth_name)//" in "//trim(depth_file))
  status = nf90_get_var(ncid, intvid, int_depth)
  if (status /= nf90_noerr) call mom_error(fatal,"MOM_diag_to_Z get_Z_depths: "//&
      trim(nf90_strerror(status))//" Reading variable "//&
      trim(int_depth_name)//" in "//trim(depth_file))
  status = nf90_get_att(ncid, intvid, "units", units)
  if (status /= nf90_noerr) units = "meters"
  status = nf90_get_att(ncid, intvid, "long_name", long_name)
  if (status /= nf90_noerr) long_name = "Depth of edges"
  edge_index = diag_axis_init(int_depth_name, int_depth, units, 'z', &
                              long_name, direction=-1)

! Create an fms axis with the same data as the cell_depth array in the input file.
  status = nf90_inq_varid(ncid, cell_depth_name, layvid)
  if (status /= nf90_noerr) call mom_error(fatal,"MOM_diag_to_Z get_Z_depths: "//&
      trim(nf90_strerror(status))//" Getting ID of variable "//&
      trim(cell_depth_name)//" in "//trim(depth_file))
  status = nf90_get_var(ncid, layvid, cell_depth)
  if (status /= nf90_noerr) call mom_error(fatal,"MOM_diag_to_Z get_Z_depths: "//&
      trim(nf90_strerror(status))//" Reading variable "//&
      trim(cell_depth_name)//" in "//trim(depth_file))
  status = nf90_get_att(ncid, layvid, "units", units)
  if (status /= nf90_noerr) units = "meters"
  status = nf90_get_att(ncid, layvid, "long_name", long_name)
  if (status /= nf90_noerr) long_name = "Depth of cell center"

  z_axis_index = diag_axis_init(cell_depth_name, cell_depth, units, 'z',&
                                long_name, edges = edge_index, direction=-1)

  deallocate(cell_depth)

  status = nf90_close(ncid)

  ! Check the sign convention and change to the MOM "height" convention.
  if (int_depth(1) < int_depth(2)) then
    do k=1,nk_out+1 ; int_depth(k) = -1*int_depth(k) ; enddo
  endif

  ! Check for inversions in grid.
  do k=1,nk_out ; if (int_depth(k) < int_depth(k+1)) then
    call mom_error(warning,"MOM_diag_to_Z get_Z_depths: "//&
        "Inverted interface depths in output grid in "//depth_file)
    nk_out = -1 ; deallocate(int_depth) ; return
  endif ; enddo

end subroutine get_z_depths

subroutine mom_diag_to_z_end(CS)
  type(diag_to_z_cs), pointer :: CS
  integer :: m

  if (ASSOCIATED(cs%u_z))   deallocate(cs%u_z)
  if (ASSOCIATED(cs%v_z))   deallocate(cs%v_z)
  if (ASSOCIATED(cs%Z_int)) deallocate(cs%Z_int)
  do m=1,cs%num_tr_used ;   deallocate(cs%tr_z(m)%p) ; enddo

  deallocate(cs)

end subroutine mom_diag_to_z_end

!> This subroutine registers a tracer to be output in depth space.
function ocean_register_diag_with_z(tr_ptr, vardesc_tr, G, Time, CS)
  type(ocean_grid_type),             intent(in) :: G      !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                             target, intent(in) :: tr_ptr !< Tracer for translation to Z-space.
  type(vardesc),                     intent(in) :: vardesc_tr !< Variable descriptor.
  type(time_type),                   intent(in) :: Time   !< Current model time.
  type(diag_to_z_cs),                pointer    :: CS     !< Control struct returned by a previous
                                                          !! call to diagnostics_init.
  integer                                       :: ocean_register_diag_with_z

!   This subroutine registers a tracer to be output in depth space.
! Arguments:
!  (in)      tr_ptr     - tracer for translation to Z-space
!  (in)      vardesc_tr - variable descriptor
!  (in)      Time       - current model time
!  (in)      G          - ocean grid structure
!  (in)      CS         - control struct returned by a previous call to diagnostics_init

  type(vardesc) :: vardesc_z
  character(len=64) :: var_name         ! A variable's name.
  integer :: isd, ied, jsd, jed, nk, m, id_test
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed ; nk = g%ke
  if (.not.associated(cs)) call mom_error(fatal, &
         "register_Z_tracer: Module must be initialized before it is used.")
  if (cs%nk_zspace<1) return

  if (cs%num_tr_used >= max_fields_) then
    call mom_error(warning,"ocean_register_diag_with_z:  Attempted to register and use "//&
                   "more than MAX_FIELDS_ z-space tracers via ocean_register_diag_with_z.")
    return
  endif

  ! register the layer tracer
  ocean_register_diag_with_z = ocean_register_diag(vardesc_tr, g, cs%diag, time)

  ! copy layer tracer variable descriptor to a z-tracer descriptor;
  ! change the name and layer information.
  vardesc_z = vardesc_tr
  call modify_vardesc(vardesc_z, z_grid="z", caller="ocean_register_diag_with_z")
  m = cs%num_tr_used + 1
  cs%missing_tr(m) = cs%missing_value ! This could be changed later, if desired.
  cs%id_tr(m) =  register_z_diag(vardesc_z, cs, time, cs%missing_tr(m))

  if (cs%nk_zspace > no_zspace) then
! There is a depth-space target file.
    if (cs%id_tr(m)>0) then
! Only allocate the tr_z field id there is a diag_table entry looking
! for it.
      cs%num_tr_used = m
      call safe_alloc_ptr(cs%tr_z(m)%p,isd,ied,jsd,jed,cs%nk_zspace)
!Can we do the following at this point?
! tr_ptr might not be allocated yet
      cs%tr_model(m)%p => tr_ptr
    endif
  else
! There is no depth-space target file but warn if a diag_table entry is
! present.
    call query_vardesc(vardesc_z, name=var_name, caller="ocean_register_diag_with_z")
    if (cs%id_tr(m)>0) call mom_error(warning, &
        "ocean_register_diag_with_z: "//trim(var_name)// &
        " cannot be output without an appropriate depth-space target file.")
  endif

end function ocean_register_diag_with_z

function register_z_diag(var_desc, CS, day, missing)
  integer                         :: register_Z_diag
  type(vardesc),       intent(in) :: var_desc
  type(diag_to_z_cs),  pointer    :: CS
  type(time_type),     intent(in) :: day
  real,                intent(in) :: missing

  character(len=64) :: var_name         ! A variable's name.
  character(len=48) :: units            ! A variable's units.
  character(len=240) :: longname        ! A variable's longname.
  character(len=8) :: hor_grid, z_grid  ! Variable grid info.
  type(axes_grp), pointer :: axes

  call query_vardesc(var_desc, name=var_name, units=units, longname=longname, &
                     hor_grid=hor_grid, z_grid=z_grid, caller="register_Zint_diag")

  ! Use the hor_grid and z_grid components of vardesc to determine the
  ! desired axes to register the diagnostic field for.
  select case (z_grid)

    case ("z")
      select case (hor_grid)
        case ("q")
          axes => cs%axesBz
        case ("h")
          axes => cs%axesTz
        case ("u")
          axes => cs%axesCuz
        case ("v")
          axes => cs%axesCvz
        case ("Bu")
          axes => cs%axesBz
        case ("T")
          axes => cs%axesTz
        case ("Cu")
          axes => cs%axesCuz
        case ("Cv")
          axes => cs%axesCvz
        case default
          call mom_error(fatal,&
            "register_Z_diag: unknown hor_grid component "//trim(hor_grid))
      end select

    case default
        call mom_error(fatal,&
          "register_Z_diag: unknown z_grid component "//trim(z_grid))
  end select

  register_z_diag = register_diag_field("ocean_model_zold", trim(var_name), axes, &
        day, trim(longname), trim(units), missing_value=missing)

end function register_z_diag

function register_zint_diag(var_desc, CS, day)
  integer                         :: register_Zint_diag
  type(vardesc),       intent(in) :: var_desc
  type(diag_to_z_cs),  pointer    :: CS
  type(time_type),     intent(in) :: day

  character(len=64) :: var_name         ! A variable's name.
  character(len=48) :: units            ! A variable's units.
  character(len=240) :: longname        ! A variable's longname.
  character(len=8) :: hor_grid          ! Variable grid info.
  type(axes_grp), pointer :: axes

  call query_vardesc(var_desc, name=var_name, units=units, longname=longname, &
                     hor_grid=hor_grid, caller="register_Zint_diag")

  if (cs%nk_zspace < 0) then
    register_zint_diag = -1 ; return
  endif

  ! Use the hor_grid and z_grid components of vardesc to determine the
  ! desired axes to register the diagnostic field for.
  select case (hor_grid)
    case ("h")
      axes => cs%axesTzi
    case ("q")
      axes => cs%axesBzi
    case ("u")
      axes => cs%axesCuzi
    case ("v")
      axes => cs%axesCvzi
    case default
      call mom_error(fatal,&
        "register_Z_diag: unknown hor_grid component "//trim(hor_grid))
  end select

  register_zint_diag = register_diag_field("ocean_model_zold", trim(var_name),&
        axes, day, trim(longname), trim(units), missing_value=cs%missing_value)

end function register_zint_diag


end module mom_diag_to_z