MOM_tracer_registry.F90

!> This module contains the tracer_registry_type and the subroutines
!! that handle registration of tracers and related subroutines.
!! The primary subroutine, register_tracer, is called to indicate the
!! tracers advected and diffused.
module mom_tracer_registry

! This file is part of MOM6. See LICENSE.md for the license.

! use MOM_diag_mediator, only : diag_ctrl
use mom_coms,          only : reproducing_sum
use mom_debugging,     only : hchksum
use mom_diag_mediator, only : diag_ctrl, register_diag_field, post_data, safe_alloc_ptr
use mom_diag_mediator, only : diag_grid_storage
use mom_diag_mediator, only : diag_copy_storage_to_diag, diag_save_grids, diag_restore_grids
use mom_error_handler, only : mom_error, fatal, warning, mom_mesg, is_root_pe
use mom_file_parser,   only : get_param, log_version, param_file_type
use mom_hor_index,     only : hor_index_type
use mom_grid,          only : ocean_grid_type
use mom_io,            only : vardesc, query_vardesc, cmor_long_std
use mom_restart,       only : register_restart_field, mom_restart_cs
use mom_string_functions, only : lowercase
use mom_time_manager,  only : time_type
use mom_unit_scaling,  only : unit_scale_type
use mom_verticalgrid,  only : verticalgrid_type

implicit none ; private

#include <MOM_memory.h>

public register_tracer
public mom_tracer_chksum, mom_tracer_chkinv
public register_tracer_diagnostics, post_tracer_diagnostics_at_sync, post_tracer_transport_diagnostics
public preale_tracer_diagnostics, postale_tracer_diagnostics
public tracer_registry_init, lock_tracer_registry, tracer_registry_end
public tracer_name_lookup

!> The tracer type
type, public :: tracer_type

  real, dimension(:,:,:), pointer :: t              => null() !< tracer concentration array [conc]
! real                            :: OBC_inflow_conc=  0.0    !< tracer concentration for generic inflows [conc]
! real, dimension(:,:,:), pointer :: OBC_in_u       => NULL() !< structured values for flow into the domain
!                                                             !! specified in OBCs through u-face of cell
! real, dimension(:,:,:), pointer :: OBC_in_v       => NULL() !< structured values for flow into the domain
!                                                             !! specified in OBCs through v-face of cell

  real, dimension(:,:,:), pointer :: ad_x           => null() !< diagnostic array for x-advective tracer flux
                                                              !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
  real, dimension(:,:,:), pointer :: ad_y           => null() !< diagnostic array for y-advective tracer flux
                                                              !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
  real, dimension(:,:),   pointer :: ad2d_x         => null() !< diagnostic vertical sum x-advective tracer flux
                                                              !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
  real, dimension(:,:),   pointer :: ad2d_y         => null() !< diagnostic vertical sum y-advective tracer flux
                                                              !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1]

  real, dimension(:,:,:), pointer :: df_x           => null() !< diagnostic array for x-diffusive tracer flux
                                                              !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
  real, dimension(:,:,:), pointer :: df_y           => null() !< diagnostic array for y-diffusive tracer flux
                                                              !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
  real, dimension(:,:,:), pointer :: lbd_dfx       => null()  !< diagnostic array for x-diffusive tracer flux
                                                              !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1]
  real, dimension(:,:,:), pointer :: lbd_dfy       => null()  !< diagnostic array for y-diffusive tracer flux
                                                              !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1]
  real, dimension(:,:), pointer :: lbd_dfx_2d       => null() !< diagnostic array for x-diffusive tracer flux
                                                              !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1]
  real, dimension(:,:), pointer :: lbd_dfy_2d       => null() !< diagnostic array for y-diffusive tracer flux
                                                              !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1]
  real, dimension(:,:), pointer :: lbd_bulk_df_x       => null() !< diagnostic array for x-diffusive tracer flux
                                                              !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1]
  real, dimension(:,:), pointer :: lbd_bulk_df_y       => null() !< diagnostic array for y-diffusive tracer flux
                                                              !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1]
  real, dimension(:,:),   pointer :: df2d_x         => null() !< diagnostic vertical sum x-diffusive flux
                                                              !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
  real, dimension(:,:),   pointer :: df2d_y         => null() !< diagnostic vertical sum y-diffusive flux
                                                              !! [conc H L2 T-1 ~> conc m3 s-1 or conc kg s-1]
!  real, dimension(:,:),   pointer :: df2d_conc_x    => NULL() !< diagnostic vertical sum x-diffusive content flux
!                                                              !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1]
!  real, dimension(:,:),   pointer :: df2d_conc_y    => NULL() !< diagnostic vertical sum y-diffusive content flux
!                                                              !! [conc H m2 s-1 ~> conc m3 s-1 or conc kg s-1]

  real, dimension(:,:,:), pointer :: advection_xy   => null() !< convergence of lateral advective tracer fluxes
                                                              !! [conc H T-1 ~> conc m s-1 or conc kg m-2 s-1]
!  real, dimension(:,:,:), pointer :: diff_cont_xy   => NULL() !< convergence of lateral diffusive tracer fluxes
!                                                              !! [conc H s-1 ~> conc m s-1 or conc kg m-2 s-1]
!  real, dimension(:,:,:), pointer :: diff_conc_xy   => NULL() !< convergence of lateral diffusive tracer fluxes
!                                                              !! expressed as a change in concentration [conc s-1]
  real, dimension(:,:,:), pointer :: t_prev         => null() !< tracer concentration array at a previous
                                                              !! timestep used for diagnostics [conc]
  real, dimension(:,:,:), pointer :: trxh_prev      => null() !< layer integrated tracer concentration array
                                                              !! at a previous timestep used for diagnostics

  character(len=32)               :: name                     !< tracer name used for diagnostics and error messages
  character(len=64)               :: units                    !< Physical dimensions of the tracer concentration
  character(len=240)              :: longname                 !< Long name of the variable
!  type(vardesc), pointer          :: vd             => NULL() !< metadata describing the tracer
  logical                         :: registry_diags = .false. !< If true, use the registry to set up the
                                                              !! diagnostics associated with this tracer.
  character(len=64)               :: cmor_name                !< CMOR name of this tracer
  character(len=64)               :: cmor_units               !< CMOR physical dimensions of the tracer
  character(len=240)              :: cmor_longname            !< CMOR long name of the tracer
  character(len=32)               :: flux_nameroot = ""       !< Short tracer name snippet used construct the
                                                              !! names of flux diagnostics.
  character(len=64)               :: flux_longname = ""       !< A word or phrase used construct the long
                                                              !! names of flux diagnostics.
  real                            :: flux_scale= 1.0          !< A scaling factor used to convert the fluxes
                                                              !! of this tracer to its desired units.
  character(len=48)               :: flux_units = ""          !< The units for fluxes of this variable.
  character(len=48)               :: conv_units = ""          !< The units for the flux convergence of this tracer.
  real                            :: conv_scale = 1.0         !< A scaling factor used to convert the flux
                                                              !! convergence of this tracer to its desired units.
  character(len=48)               :: cmor_tendprefix = ""     !< The CMOR variable prefix for tendencies of this
                                                              !! tracer, required because CMOR does not follow any
                                                              !! discernable pattern for these names.
  integer :: ind_tr_squared = -1 !< The tracer registry index for the square of this tracer

  !### THESE CAPABILITIES HAVE NOT YET BEEN IMPLEMENTED.
  ! logical :: advect_tr = .true.     !< If true, this tracer should be advected
  ! logical :: hordiff_tr = .true.    !< If true, this tracer should experience epineutral diffusion
  logical :: remap_tr = .true.      !< If true, this tracer should be vertically remapped

  integer :: diag_form = 1  !< An integer indicating which template is to be used to label diagnostics.
  !>@{ Diagnostic IDs
  integer :: id_tr = -1, id_tr_post_horzn = -1
  integer :: id_adx = -1, id_ady = -1, id_dfx = -1, id_dfy = -1
  integer :: id_lbd_bulk_dfx = -1, id_lbd_bulk_dfy = -1, id_lbd_dfx = -1, id_lbd_dfy = -1
  integer :: id_lbd_dfx_2d = -1  , id_lbd_dfy_2d = -1
  integer :: id_adx_2d = -1, id_ady_2d = -1, id_dfx_2d = -1, id_dfy_2d = -1
  integer :: id_adv_xy = -1, id_adv_xy_2d = -1
  integer :: id_dfxy_cont = -1, id_dfxy_cont_2d = -1, id_dfxy_conc = -1
  integer :: id_lbdxy_cont = -1, id_lbdxy_cont_2d = -1, id_lbdxy_conc = -1
  integer :: id_remap_conc = -1, id_remap_cont = -1, id_remap_cont_2d = -1
  integer :: id_tendency = -1, id_trxh_tendency = -1, id_trxh_tendency_2d = -1
  integer :: id_tr_vardec = -1
  !>@}
end type tracer_type

!> Type to carry basic tracer information
type, public :: tracer_registry_type
  integer                  :: ntr = 0           !< number of registered tracers
  type(tracer_type)        :: tr(max_fields_)   !< array of registered tracers
! type(diag_ctrl), pointer :: diag              !< structure to regulate timing of diagnostics
  logical                  :: locked = .false.  !< New tracers may be registered if locked=.false.
                                                !! When locked=.true., no more tracers can be registered,
                                                !! at which point common diagnostics can be set up
                                                !! for the registered tracers
end type tracer_registry_type

contains

!> This subroutine registers a tracer to be advected and laterally diffused.
subroutine register_tracer(tr_ptr, Reg, param_file, HI, GV, name, longname, units, &
                           cmor_name, cmor_units, cmor_longname, tr_desc, &
                           OBC_inflow, OBC_in_u, OBC_in_v, ad_x, ad_y, df_x, df_y, &
                           ad_2d_x, ad_2d_y, df_2d_x, df_2d_y, advection_xy, registry_diags, &
                           flux_nameroot, flux_longname, flux_units, flux_scale, &
                           convergence_units, convergence_scale, cmor_tendprefix, diag_form, &
                           restart_CS, mandatory)
  type(hor_index_type),           intent(in)    :: hi           !< horizontal index type
  type(verticalgrid_type),        intent(in)    :: gv           !< ocean vertical grid structure
  type(tracer_registry_type),     pointer       :: reg          !< pointer to the tracer registry
  real, dimension(SZI_(HI),SZJ_(HI),SZK_(GV)), &
                                  target        :: tr_ptr       !< target or pointer to the tracer array
  type(param_file_type), intent(in)             :: param_file   !< file to parse for model parameter values
  character(len=*),     optional, intent(in)    :: name         !< Short tracer name
  character(len=*),     optional, intent(in)    :: longname     !< The long tracer name
  character(len=*),     optional, intent(in)    :: units        !< The units of this tracer
  character(len=*),     optional, intent(in)    :: cmor_name    !< CMOR name
  character(len=*),     optional, intent(in)    :: cmor_units   !< CMOR physical dimensions of variable
  character(len=*),     optional, intent(in)    :: cmor_longname !< CMOR long name
  type(vardesc),        optional, intent(in)    :: tr_desc      !< A structure with metadata about the tracer

  real,                 optional, intent(in)    :: obc_inflow   !< the tracer for all inflows via OBC for which OBC_in_u
                                                                !! or OBC_in_v are not specified (units of tracer CONC)
  real, dimension(:,:,:), optional, pointer     :: obc_in_u     !< tracer at inflows through u-faces of
                                                                !! tracer cells (units of tracer CONC)
  real, dimension(:,:,:), optional, pointer     :: obc_in_v     !< tracer at inflows through v-faces of
                                                                !! tracer cells (units of tracer CONC)

  ! The following are probably not necessary if registry_diags is present and true.
  real, dimension(:,:,:), optional, pointer     :: ad_x         !< diagnostic x-advective flux
                                                                !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1]
  real, dimension(:,:,:), optional, pointer     :: ad_y         !< diagnostic y-advective flux
                                                                !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1]
  real, dimension(:,:,:), optional, pointer     :: df_x         !< diagnostic x-diffusive flux
                                                                !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1]
  real, dimension(:,:,:), optional, pointer     :: df_y         !< diagnostic y-diffusive flux
                                                                !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1]
  real, dimension(:,:),   optional, pointer     :: ad_2d_x      !< vert sum of diagnostic x-advect flux
                                                                !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1]
  real, dimension(:,:),   optional, pointer     :: ad_2d_y      !< vert sum of diagnostic y-advect flux
                                                                !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1]
  real, dimension(:,:),   optional, pointer     :: df_2d_x      !< vert sum of diagnostic x-diffuse flux
                                                                !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1]
  real, dimension(:,:),   optional, pointer     :: df_2d_y      !< vert sum of diagnostic y-diffuse flux
                                                                !! [conc H L2 T-1 ~> CONC m3 s-1 or CONC kg s-1]

  real, dimension(:,:,:), optional, pointer     :: advection_xy !< convergence of lateral advective tracer fluxes
  logical,              optional, intent(in)    :: registry_diags !< If present and true, use the registry for
                                                                !! the diagnostics of this tracer.
  character(len=*),     optional, intent(in)    :: flux_nameroot !< Short tracer name snippet used construct the
                                                                !! names of flux diagnostics.
  character(len=*),     optional, intent(in)    :: flux_longname !< A word or phrase used construct the long
                                                                !! names of flux diagnostics.
  character(len=*),     optional, intent(in)    :: flux_units   !< The units for the fluxes of this tracer.
  real,                 optional, intent(in)    :: flux_scale   !< A scaling factor used to convert the fluxes
                                                                !! of this tracer to its desired units.
  character(len=*),     optional, intent(in)    :: convergence_units !< The units for the flux convergence of
                                                                !! this tracer.
  real,                 optional, intent(in)    :: convergence_scale !< A scaling factor used to convert the flux
                                                                !! convergence of this tracer to its desired units.
  character(len=*),     optional, intent(in)    :: cmor_tendprefix !< The CMOR name for the layer-integrated
                                                                !! tendencies of this tracer.
  integer,              optional, intent(in)    :: diag_form    !< An integer (1 or 2, 1 by default) indicating the
                                                                !! character string template to use in
                                                                !! labeling diagnostics
  type(mom_restart_cs), optional, pointer       :: restart_cs   !< A pointer to the restart control structure
                                                                !! this tracer will be registered for
                                                                !! restarts if this argument is present
  logical,              optional, intent(in)    :: mandatory    !< If true, this tracer must be read
                                                                !! from a restart file.

  logical :: mand
  type(tracer_type), pointer :: tr=>null()
  character(len=256) :: mesg    ! Message for error messages.

  if (.not. associated(reg)) call tracer_registry_init(param_file, reg)

  if (reg%ntr>=max_fields_) then
    write(mesg,'("Increase MAX_FIELDS_ in MOM_memory.h to at least ",I3," to allow for &
        &all the tracers being registered via register_tracer.")') reg%ntr+1
    call mom_error(fatal,"MOM register_tracer: "//mesg)
  endif
  reg%ntr = reg%ntr + 1

  tr => reg%Tr(reg%ntr)

  if (present(name)) then
    tr%name = name
    tr%longname = name ; if (present(longname)) tr%longname = longname
    tr%units = "Conc" ; if (present(units)) tr%units = units

    tr%cmor_name = ""
    if (present(cmor_name)) tr%cmor_name = cmor_name

    tr%cmor_units = tr%units
    if (present(cmor_units)) tr%cmor_units = cmor_units

    tr%cmor_longname = ""
    if (present(cmor_longname)) tr%cmor_longname = cmor_longname

    if (present(tr_desc)) call mom_error(warning, "MOM register_tracer: "//&
      "It is a bad idea to use both name and tr_desc when registring "//trim(name))
  elseif (present(tr_desc)) then
    call query_vardesc(tr_desc, name=tr%name, units=tr%units, &
                       longname=tr%longname, cmor_field_name=tr%cmor_name, &
                       cmor_longname=tr%cmor_longname, caller="register_tracer")
    tr%cmor_units = tr%units
  else
    call mom_error(fatal,"MOM register_tracer: Either name or "//&
                   "tr_desc must be present when registering a tracer.")
  endif

  if (reg%locked) call mom_error(fatal, &
      "MOM register_tracer was called for variable "//trim(tr%name)//&
      " with a locked tracer registry.")

  tr%flux_nameroot = tr%name
  if (present(flux_nameroot)) then
    if (len_trim(flux_nameroot) > 0) tr%flux_nameroot = flux_nameroot
  endif

  tr%flux_longname = tr%longname
  if (present(flux_longname)) then
    if (len_trim(flux_longname) > 0) tr%flux_longname = flux_longname
  endif

  tr%flux_units = ""
  if (present(flux_units)) tr%flux_units = flux_units

  tr%flux_scale = 1.0
  if (present(flux_scale)) tr%flux_scale = flux_scale

  tr%conv_units = ""
  if (present(convergence_units)) tr%conv_units = convergence_units

  tr%cmor_tendprefix = ""
  if (present(cmor_tendprefix)) tr%cmor_tendprefix = cmor_tendprefix

  tr%conv_scale = 1.0
  if (present(convergence_scale)) then
    tr%conv_scale = convergence_scale
  elseif (present(flux_scale)) then
    tr%conv_scale = flux_scale
  endif

  tr%diag_form = 1
  if (present(diag_form)) tr%diag_form = diag_form

  tr%t => tr_ptr

  if (present(registry_diags)) tr%registry_diags = registry_diags

  if (present(ad_x)) then ; if (associated(ad_x)) tr%ad_x => ad_x ; endif
  if (present(ad_y)) then ; if (associated(ad_y)) tr%ad_y => ad_y ; endif
  if (present(df_x)) then ; if (associated(df_x)) tr%df_x => df_x ; endif
  if (present(df_y)) then ; if (associated(df_y)) tr%df_y => df_y ; endif
! if (present(OBC_inflow)) Tr%OBC_inflow_conc = OBC_inflow
! if (present(OBC_in_u)) then ; if (associated(OBC_in_u)) &
!                                   Tr%OBC_in_u => OBC_in_u ; endif
! if (present(OBC_in_v)) then ; if (associated(OBC_in_v)) &
!                                   Tr%OBC_in_v => OBC_in_v ; endif
  if (present(ad_2d_x)) then ; if (associated(ad_2d_x)) tr%ad2d_x => ad_2d_x ; endif
  if (present(ad_2d_y)) then ; if (associated(ad_2d_y)) tr%ad2d_y => ad_2d_y ; endif
  if (present(df_2d_x)) then ; if (associated(df_2d_x)) tr%df2d_x => df_2d_x ; endif

  if (present(advection_xy)) then ; if (associated(advection_xy)) tr%advection_xy => advection_xy ; endif

  if (present(restart_cs)) then ; if (associated(restart_cs)) then
    ! Register this tracer to be read from and written to restart files.
    mand = .true. ; if (present(mandatory)) mand = mandatory

    call register_restart_field(tr_ptr, tr%name, mand, restart_cs, &
                                longname=tr%longname, units=tr%units)
  endif ; endif

end subroutine register_tracer


!> This subroutine locks the tracer registry to prevent the addition of more
!! tracers.  After locked=.true., can then register common diagnostics.
subroutine lock_tracer_registry(Reg)
  type(tracer_registry_type), pointer    :: reg    !< pointer to the tracer registry

  if (.not. associated(reg)) call mom_error(warning, &
    "lock_tracer_registry called with an unassociated registry.")

  reg%locked = .true.

end subroutine lock_tracer_registry

!> register_tracer_diagnostics does a set of register_diag_field calls for any previously
!! registered in a tracer registry with a value of registry_diags set to .true.
subroutine register_tracer_diagnostics(Reg, h, Time, diag, G, GV, US, use_ALE)
  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(tracer_registry_type), pointer    :: reg  !< pointer to the tracer registry
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                              intent(in) :: h    !< Layer thicknesses
  type(time_type),            intent(in) :: time !< current model time
  type(diag_ctrl),            intent(in) :: diag !< structure to regulate diagnostic output
  logical,                    intent(in) :: use_ale !< If true active diagnostics that only
                                                 !! apply to ALE configurations

  character(len=24) :: name     ! A variable's name in a NetCDF file.
  character(len=24) :: shortnm  ! A shortened version of a variable's name for
                                ! creating additional diagnostics.
  character(len=72) :: longname ! The long name of that tracer variable.
  character(len=72) :: flux_longname ! The tracer name in the long names of fluxes.
  character(len=48) :: units    ! The dimensions of the tracer.
  character(len=48) :: flux_units ! The units for fluxes, either
                                ! [units] m3 s-1 or [units] kg s-1.
  character(len=48) :: conv_units ! The units for flux convergences, either
                                ! [units] m2 s-1 or [units] kg s-1.
  character(len=48) :: unit2    ! The dimensions of the tracer squared
  character(len=72)  :: cmorname ! The CMOR name of this tracer.
  character(len=120) :: cmor_longname ! The CMOR long name of that variable.
  character(len=120) :: var_lname      ! A temporary longname for a diagnostic.
  character(len=120) :: cmor_var_lname ! The temporary CMOR long name for a diagnostic
  character(len=72)  :: cmor_varname ! The temporary CMOR name for a diagnostic
  type(tracer_type), pointer :: tr=>null()
  integer :: i, j, k, is, ie, js, je, nz, m, m2, ntr_in
  integer :: isd, ied, jsd, jed, isdb, iedb, jsdb, jedb
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  isd  = g%isd  ; ied  = g%ied  ; jsd  = g%jsd  ; jed  = g%jed
  isdb = g%IsdB ; iedb = g%IedB ; jsdb = g%JsdB ; jedb = g%JedB

  if (.not. associated(reg)) call mom_error(fatal, "register_tracer_diagnostics: "//&
       "register_tracer must be called before register_tracer_diagnostics")

  ntr_in = reg%ntr

  do m=1,ntr_in ; if (reg%Tr(m)%registry_diags) then
    tr => reg%Tr(m)
!    call query_vardesc(Tr%vd, name, units=units, longname=longname, &
!                       cmor_field_name=cmorname, cmor_longname=cmor_longname, &
!                       caller="register_tracer_diagnostics")
    name = tr%name ; units=adjustl(tr%units) ; longname = tr%longname
    cmorname = tr%cmor_name ; cmor_longname = tr%cmor_longname
    shortnm = tr%flux_nameroot
    flux_longname = tr%flux_longname
    if (len_trim(cmor_longname) == 0) cmor_longname = longname

    if (len_trim(tr%flux_units) > 0) then ; flux_units = tr%flux_units
    elseif (gv%Boussinesq) then ; flux_units = trim(units)//" m3 s-1"
    else ; flux_units = trim(units)//" kg s-1" ; endif

    if (len_trim(tr%conv_units) > 0) then ; conv_units = tr%conv_units
    elseif (gv%Boussinesq) then ; conv_units = trim(units)//" m s-1"
    else ; conv_units = trim(units)//" kg m-2 s-1" ; endif

    if (len_trim(cmorname) == 0) then
      tr%id_tr = register_diag_field("ocean_model", trim(name), diag%axesTL, &
        time, trim(longname), trim(units))
    else
      tr%id_tr = register_diag_field("ocean_model", trim(name), diag%axesTL, &
        time, trim(longname), trim(units), cmor_field_name=cmorname, &
        cmor_long_name=cmor_longname, cmor_units=tr%cmor_units, &
        cmor_standard_name=cmor_long_std(cmor_longname))
    endif
    tr%id_tr_post_horzn = register_diag_field("ocean_model",                &
      trim(name)//"_post_horzn", diag%axesTL, time,                         &
      trim(longname)//" after horizontal transport (advection/diffusion) "//&
      "has occurred", trim(units))
    if (tr%diag_form == 1) then
      tr%id_adx = register_diag_field("ocean_model", trim(shortnm)//"_adx", &
          diag%axesCuL, time, trim(flux_longname)//" advective zonal flux" , &
          trim(flux_units), v_extensive = .true., y_cell_method = 'sum', &
          conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T)
      tr%id_ady = register_diag_field("ocean_model", trim(shortnm)//"_ady", &
          diag%axesCvL, time, trim(flux_longname)//" advective meridional flux" , &
          trim(flux_units), v_extensive = .true., x_cell_method = 'sum', &
          conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T)
      tr%id_dfx = register_diag_field("ocean_model", trim(shortnm)//"_dfx", &
          diag%axesCuL, time, trim(flux_longname)//" diffusive zonal flux" , &
          trim(flux_units), v_extensive = .true., y_cell_method = 'sum', &
          conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T)
      tr%id_dfy = register_diag_field("ocean_model", trim(shortnm)//"_dfy", &
          diag%axesCvL, time, trim(flux_longname)//" diffusive zonal flux" , &
          trim(flux_units), v_extensive = .true., x_cell_method = 'sum', &
          conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T)
      tr%id_lbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffx", &
          diag%axesCuL, time, trim(flux_longname)//" diffusive zonal flux from the lateral boundary diffusion "//&
          "scheme", trim(flux_units), v_extensive = .true., y_cell_method = 'sum', &
          conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T)
      tr%id_lbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffy", &
          diag%axesCvL, time, trim(flux_longname)//" diffusive meridional flux from the lateral boundary diffusion "//&
          "scheme", trim(flux_units), v_extensive = .true., x_cell_method = 'sum', &
          conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T)
    else
      tr%id_adx = register_diag_field("ocean_model", trim(shortnm)//"_adx", &
          diag%axesCuL, time, "Advective (by residual mean) Zonal Flux of "//trim(flux_longname), &
          flux_units, v_extensive=.true., conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T, y_cell_method = 'sum')
      tr%id_ady = register_diag_field("ocean_model", trim(shortnm)//"_ady", &
          diag%axesCvL, time, "Advective (by residual mean) Meridional Flux of "//trim(flux_longname), &
          flux_units, v_extensive=.true., conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T, x_cell_method = 'sum')
      tr%id_dfx = register_diag_field("ocean_model", trim(shortnm)//"_diffx", &
          diag%axesCuL, time, "Diffusive Zonal Flux of "//trim(flux_longname), &
          flux_units, v_extensive=.true., conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
          y_cell_method='sum')
      tr%id_dfy = register_diag_field("ocean_model", trim(shortnm)//"_diffy", &
          diag%axesCvL, time, "Diffusive Meridional Flux of "//trim(flux_longname), &
          flux_units, v_extensive=.true., conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
          x_cell_method='sum')
      tr%id_lbd_dfx = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffx", &
          diag%axesCuL, time, "Lateral Boundary Diffusive Zonal Flux of "//trim(flux_longname), &
          flux_units, v_extensive=.true., conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
          y_cell_method='sum')
      tr%id_lbd_dfy = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffy", &
          diag%axesCvL, time, "Lateral Boundary Diffusive Meridional Flux of "//trim(flux_longname), &
          flux_units, v_extensive=.true., conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
          x_cell_method='sum')
    endif
    if (tr%id_adx > 0) call safe_alloc_ptr(tr%ad_x,isdb,iedb,jsd,jed,nz)
    if (tr%id_ady > 0) call safe_alloc_ptr(tr%ad_y,isd,ied,jsdb,jedb,nz)
    if (tr%id_dfx > 0) call safe_alloc_ptr(tr%df_x,isdb,iedb,jsd,jed,nz)
    if (tr%id_dfy > 0) call safe_alloc_ptr(tr%df_y,isd,ied,jsdb,jedb,nz)
    if (tr%id_lbd_dfx > 0) call safe_alloc_ptr(tr%lbd_dfx,isdb,iedb,jsd,jed,nz)
    if (tr%id_lbd_dfy > 0) call safe_alloc_ptr(tr%lbd_dfy,isd,ied,jsdb,jedb,nz)

    tr%id_adx_2d = register_diag_field("ocean_model", trim(shortnm)//"_adx_2d", &
        diag%axesCu1, time, &
        "Vertically Integrated Advective Zonal Flux of "//trim(flux_longname), &
        flux_units, conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T, y_cell_method = 'sum')
    tr%id_ady_2d = register_diag_field("ocean_model", trim(shortnm)//"_ady_2d", &
        diag%axesCv1, time, &
        "Vertically Integrated Advective Meridional Flux of "//trim(flux_longname), &
        flux_units, conversion=tr%flux_scale*(us%L_to_m**2)*us%s_to_T, x_cell_method = 'sum')
    tr%id_dfx_2d = register_diag_field("ocean_model", trim(shortnm)//"_diffx_2d", &
        diag%axesCu1, time, &
        "Vertically Integrated Diffusive Zonal Flux of "//trim(flux_longname), &
        flux_units, conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
        y_cell_method='sum')
    tr%id_dfy_2d = register_diag_field("ocean_model", trim(shortnm)//"_diffy_2d", &
        diag%axesCv1, time, &
        "Vertically Integrated Diffusive Meridional Flux of "//trim(flux_longname), &
        flux_units, conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
        x_cell_method='sum')
    tr%id_lbd_bulk_dfx = register_diag_field("ocean_model", trim(shortnm)//"_lbd_bulk_diffx", &
        diag%axesCu1, time, &
        "Total Bulk Diffusive Zonal Flux of "//trim(flux_longname), &
        flux_units, conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
        y_cell_method='sum')
    tr%id_lbd_bulk_dfy = register_diag_field("ocean_model", trim(shortnm)//"_lbd_bulk_diffy", &
        diag%axesCv1, time, &
        "Total Bulk Diffusive Meridional Flux of "//trim(flux_longname), &
        flux_units, conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
        x_cell_method='sum')
    tr%id_lbd_dfx_2d = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffx_2d", &
          diag%axesCu1, time, "Vertically-integrated zonal diffusive flux from the lateral boundary diffusion "//&
          "scheme for "//trim(flux_longname), flux_units, conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
          y_cell_method = 'sum')
    tr%id_lbd_dfy_2d = register_diag_field("ocean_model", trim(shortnm)//"_lbd_diffy_2d", &
          diag%axesCv1, time, "Vertically-integrated meridional diffusive flux from the lateral boundary diffusion "//&
          "scheme for "//trim(flux_longname), flux_units, conversion=(us%L_to_m**2)*tr%flux_scale*us%s_to_T, &
           x_cell_method = 'sum')

    if (tr%id_adx_2d > 0) call safe_alloc_ptr(tr%ad2d_x,isdb,iedb,jsd,jed)
    if (tr%id_ady_2d > 0) call safe_alloc_ptr(tr%ad2d_y,isd,ied,jsdb,jedb)
    if (tr%id_dfx_2d > 0) call safe_alloc_ptr(tr%df2d_x,isdb,iedb,jsd,jed)
    if (tr%id_dfy_2d > 0) call safe_alloc_ptr(tr%df2d_y,isd,ied,jsdb,jedb)
    if (tr%id_lbd_bulk_dfx > 0) call safe_alloc_ptr(tr%lbd_bulk_df_x,isdb,iedb,jsd,jed)
    if (tr%id_lbd_bulk_dfy > 0) call safe_alloc_ptr(tr%lbd_bulk_df_y,isd,ied,jsdb,jedb)
    if (tr%id_lbd_dfx_2d > 0) call safe_alloc_ptr(tr%lbd_dfx_2d,isdb,iedb,jsd,jed)
    if (tr%id_lbd_dfy_2d > 0) call safe_alloc_ptr(tr%lbd_dfy_2d,isd,ied,jsdb,jedb)

    tr%id_adv_xy = register_diag_field('ocean_model', trim(shortnm)//"_advection_xy", &
        diag%axesTL, time, &
        'Horizontal convergence of residual mean advective fluxes of '//&
        trim(lowercase(flux_longname)), conv_units, v_extensive=.true., &
        conversion=tr%conv_scale*us%s_to_T)
    tr%id_adv_xy_2d = register_diag_field('ocean_model', trim(shortnm)//"_advection_xy_2d", &
        diag%axesT1, time, &
        'Vertical sum of horizontal convergence of residual mean advective fluxes of '//&
        trim(lowercase(flux_longname)), conv_units, conversion=tr%conv_scale*us%s_to_T)
    if ((tr%id_adv_xy > 0) .or. (tr%id_adv_xy_2d > 0)) &
      call safe_alloc_ptr(tr%advection_xy,isd,ied,jsd,jed,nz)

    tr%id_tendency = register_diag_field('ocean_model', trim(shortnm)//'_tendency', &
        diag%axesTL, time, &
        'Net time tendency for '//trim(lowercase(longname)), trim(units)//' s-1', conversion=us%s_to_T)

    if (tr%id_tendency > 0) then
      call safe_alloc_ptr(tr%t_prev,isd,ied,jsd,jed,nz)
      do k=1,nz ; do j=js,je ; do i=is,ie
        tr%t_prev(i,j,k) = tr%t(i,j,k)
      enddo ; enddo ; enddo
    endif

    ! Neutral/Lateral diffusion convergence tendencies
    if (tr%diag_form == 1) then
      tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', &
          diag%axesTL, time, "Neutral diffusion tracer content tendency for "//trim(shortnm), &
          conv_units, conversion=tr%conv_scale*us%s_to_T, x_cell_method='sum', y_cell_method='sum', v_extensive=.true.)

      tr%id_dfxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency_2d', &
          diag%axesT1, time, "Depth integrated neutral diffusion tracer content "//&
          "tendency for "//trim(shortnm), conv_units, conversion=tr%conv_scale*us%s_to_T, &
          x_cell_method='sum', y_cell_method= 'sum')

      tr%id_lbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency', &
          diag%axesTL, time, "Lateral diffusion tracer content tendency for "//trim(shortnm), &
          conv_units, conversion=tr%conv_scale*us%s_to_T, x_cell_method='sum', y_cell_method='sum', v_extensive=.true.)

      tr%id_lbdxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency_2d', &
          diag%axesT1, time, "Depth integrated lateral diffusion tracer content "//&
          "tendency for "//trim(shortnm), conv_units, conversion=tr%conv_scale*us%s_to_T, &
          x_cell_method='sum', y_cell_method= 'sum')
    else
      cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//&
           ' expressed as '//trim(lowercase(flux_longname))//&
           ' content due to parameterized mesoscale neutral diffusion'
      tr%id_dfxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency', &
          diag%axesTL, time, "Neutral diffusion tracer content tendency for "//trim(shortnm), &
          conv_units, conversion=tr%conv_scale*us%s_to_T, cmor_field_name = trim(tr%cmor_tendprefix)//'pmdiff', &
          cmor_long_name = trim(cmor_var_lname), cmor_standard_name = trim(cmor_long_std(cmor_var_lname)), &
          x_cell_method = 'sum', y_cell_method = 'sum', v_extensive = .true.)

      cmor_var_lname = 'Tendency of '//trim(lowercase(cmor_longname))//' expressed as '//&
                       trim(lowercase(flux_longname))//' content due to parameterized mesoscale neutral diffusion'
      tr%id_dfxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_cont_tendency_2d', &
          diag%axesT1, time, "Depth integrated neutral diffusion tracer "//&
          "content tendency for "//trim(shortnm), conv_units, &
          conversion=tr%conv_scale*us%s_to_T, cmor_field_name=trim(tr%cmor_tendprefix)//'pmdiff_2d', &
          cmor_long_name=trim(cmor_var_lname), cmor_standard_name=trim(cmor_long_std(cmor_var_lname)), &
          x_cell_method='sum', y_cell_method='sum')

      tr%id_lbdxy_cont = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency', &
          diag%axesTL, time, "Lateral diffusion tracer content tendency for "//trim(shortnm), &
          conv_units, conversion=tr%conv_scale*us%s_to_T, &
          x_cell_method = 'sum', y_cell_method = 'sum', v_extensive = .true.)

      tr%id_lbdxy_cont_2d = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_cont_tendency_2d', &
          diag%axesT1, time, "Depth integrated lateral diffusion tracer "//&
          "content tendency for "//trim(shortnm), conv_units, &
          conversion=tr%conv_scale*us%s_to_T, x_cell_method='sum', y_cell_method='sum')
    endif
    tr%id_dfxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_dfxy_conc_tendency', &
        diag%axesTL, time, "Neutral diffusion tracer concentration tendency for "//trim(shortnm), &
        trim(units)//' s-1', conversion=us%s_to_T)

    tr%id_lbdxy_conc = register_diag_field("ocean_model", trim(shortnm)//'_lbdxy_conc_tendency', &
        diag%axesTL, time, "Lateral diffusion tracer concentration tendency for "//trim(shortnm), &
        trim(units)//' s-1', conversion=us%s_to_T)

    var_lname = "Net time tendency for "//lowercase(flux_longname)
    if (len_trim(tr%cmor_tendprefix) == 0) then
      tr%id_trxh_tendency = register_diag_field('ocean_model', trim(shortnm)//'h_tendency', &
          diag%axesTL, time, var_lname, conv_units, v_extensive=.true., &
          conversion=tr%conv_scale*us%s_to_T)
      tr%id_trxh_tendency_2d = register_diag_field('ocean_model', trim(shortnm)//'h_tendency_2d', &
          diag%axesT1, time, "Vertical sum of "//trim(lowercase(var_lname)), conv_units, &
          conversion=tr%conv_scale*us%s_to_T)
    else
      cmor_var_lname = "Tendency of "//trim(cmor_longname)//" Expressed as "//&
                        trim(flux_longname)//" Content"
      tr%id_trxh_tendency = register_diag_field('ocean_model', trim(shortnm)//'h_tendency', &
          diag%axesTL, time, var_lname, conv_units, &
          cmor_field_name=trim(tr%cmor_tendprefix)//"tend", &
          cmor_standard_name=cmor_long_std(cmor_var_lname), cmor_long_name=cmor_var_lname, &
          v_extensive=.true., conversion=tr%conv_scale*us%s_to_T)
      cmor_var_lname = trim(cmor_var_lname)//" Vertical Sum"
      tr%id_trxh_tendency_2d = register_diag_field('ocean_model', trim(shortnm)//'h_tendency_2d', &
          diag%axesT1, time, "Vertical sum of "//trim(lowercase(var_lname)), conv_units, &
          cmor_field_name=trim(tr%cmor_tendprefix)//"tend_2d", &
          cmor_standard_name=cmor_long_std(cmor_var_lname), cmor_long_name=cmor_var_lname, &
          conversion=tr%conv_scale*us%s_to_T)
    endif
    if ((tr%id_trxh_tendency > 0) .or. (tr%id_trxh_tendency_2d > 0)) then
      call safe_alloc_ptr(tr%Trxh_prev,isd,ied,jsd,jed,nz)
      do k=1,nz ; do j=js,je ; do i=is,ie
        tr%Trxh_prev(i,j,k) = tr%t(i,j,k) * h(i,j,k)
      enddo ; enddo ; enddo
    endif

    ! Vertical regridding/remapping tendencies
    if (use_ale .and. tr%remap_tr) then
      var_lname = "Vertical remapping tracer concentration tendency for "//trim(reg%Tr(m)%name)
      tr%id_remap_conc= register_diag_field('ocean_model',                          &
        trim(tr%flux_nameroot)//'_tendency_vert_remap', diag%axesTL, time, var_lname, &
        trim(units)//' s-1', conversion=us%s_to_T)

      var_lname = "Vertical remapping tracer content tendency for "//trim(reg%Tr(m)%flux_longname)
      tr%id_remap_cont = register_diag_field('ocean_model', &
        trim(tr%flux_nameroot)//'h_tendency_vert_remap',         &
        diag%axesTL, time, var_lname, flux_units, v_extensive=.true., conversion=tr%conv_scale*us%s_to_T)

      var_lname = "Vertical sum of vertical remapping tracer content tendency for "//&
                  trim(reg%Tr(m)%flux_longname)
      tr%id_remap_cont_2d = register_diag_field('ocean_model', &
        trim(tr%flux_nameroot)//'h_tendency_vert_remap_2d',         &
        diag%axesT1, time, var_lname, flux_units, conversion=tr%conv_scale*us%s_to_T)

    endif

    if (use_ale .and. (reg%ntr<max_fields_) .and. tr%remap_tr) then
      unit2 = trim(units)//"2"
      if (index(units(1:len_trim(units))," ") > 0) unit2 = "("//trim(units)//")2"
      tr%id_tr_vardec = register_diag_field('ocean_model', trim(shortnm)//"_vardec", diag%axesTL, &
        time, "ALE variance decay for "//lowercase(longname), trim(unit2)//" s-1", conversion=us%s_to_T)
      if (tr%id_tr_vardec > 0) then
        ! Set up a new tracer for this tracer squared
        m2 = reg%ntr+1
        tr%ind_tr_squared = m2
        call safe_alloc_ptr(reg%Tr(m2)%t,isd,ied,jsd,jed,nz) ; reg%Tr(m2)%t(:,:,:) = 0.0
        reg%Tr(m2)%name = trim(shortnm)//"2"
        reg%Tr(m2)%longname = "Squared "//trim(longname)
        reg%Tr(m2)%units = unit2
        reg%Tr(m2)%registry_diags = .false.
        reg%Tr(m2)%ind_tr_squared = -1
        ! Augment the total number of tracers, including the squared tracers.
        reg%ntr = reg%ntr + 1
      endif
    endif

  endif ; enddo

end subroutine register_tracer_diagnostics

subroutine preale_tracer_diagnostics(Reg, G, GV)
  type(tracer_registry_type), pointer    :: reg  !< pointer to the tracer registry
  type(ocean_grid_type),      intent(in) :: g    !< The ocean's grid structure
  type(verticalgrid_type),    intent(in) :: gv   !< ocean vertical grid structure

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

  do m=1,reg%ntr ; if (reg%Tr(m)%ind_tr_squared > 0) then
    m2 = reg%Tr(m)%ind_tr_squared
  ! Update squared quantities
    do k=1,nz ; do j=js,je ; do i=is,ie
      reg%Tr(m2)%T(i,j,k) = reg%Tr(m)%T(i,j,k)**2
    enddo ; enddo ; enddo
  endif ; enddo

end subroutine preale_tracer_diagnostics

subroutine postale_tracer_diagnostics(Reg, G, GV, diag, dt)
  type(tracer_registry_type), pointer    :: reg  !< pointer to the tracer registry
  type(ocean_grid_type),      intent(in) :: g    !< The ocean's grid structure
  type(verticalgrid_type),    intent(in) :: gv   !< ocean vertical grid structure
  type(diag_ctrl),            intent(in) :: diag !< regulates diagnostic output
  real,                       intent(in) :: dt   !< total time interval for these diagnostics [T ~> s]

  real    :: work(szi_(g),szj_(g),szk_(g))
  real    :: idt ! The inverse of the time step [T-1 ~> s-1]
  integer :: i, j, k, is, ie, js, je, nz, m, m2
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  ! The "if" is to avoid NaNs if the diagnostic is called for a zero length interval
  idt = 0.0 ; if (dt /= 0.0) idt = 1.0 / dt

  do m=1,reg%ntr ; if (reg%Tr(m)%id_tr_vardec > 0) then
    m2 = reg%Tr(m)%ind_tr_squared
    if (m2 < 1) call mom_error(fatal, "Bad value of Tr%ind_tr_squared for "//trim(reg%Tr(m)%name))
  ! Update squared quantities
    do k=1,nz ; do j=js,je ; do i=is,ie
      work(i,j,k) = (reg%Tr(m2)%T(i,j,k) - reg%Tr(m)%T(i,j,k)**2) * idt
    enddo ; enddo ; enddo
    call post_data(reg%Tr(m)%id_tr_vardec, work, diag)
  endif ; enddo

end subroutine postale_tracer_diagnostics

!> Post tracer diganostics when that should only be posted when MOM's state
!! is self-consistent (also referred to as 'synchronized')
subroutine post_tracer_diagnostics_at_sync(Reg, h, diag_prev, diag, G, GV, dt)
  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(tracer_registry_type), pointer    :: reg  !< pointer to the tracer registry
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                              intent(in) :: h    !< Layer thicknesses
  type(diag_grid_storage),    intent(in) :: diag_prev !< Contains diagnostic grids from previous timestep
  type(diag_ctrl),            intent(inout) :: diag !< structure to regulate diagnostic output
  real,                       intent(in) :: dt   !< total time step for tracer updates [T ~> s]

  real    :: work3d(szi_(g),szj_(g),szk_(gv))
  real    :: work2d(szi_(g),szj_(g))
  real    :: idt ! The inverse of the time step [T-1 ~> s-1]
  type(tracer_type), pointer :: tr=>null()
  integer :: i, j, k, is, ie, js, je, nz, m
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  idt = 0.; if (dt/=0.) idt = 1.0 / dt ! The "if" is in case the diagnostic is called for a zero length interval

  ! Tendency diagnostics need to be posted on the grid from the last call to this routine
  call diag_save_grids(diag)
  call diag_copy_storage_to_diag(diag, diag_prev)
  do m=1,reg%ntr ; if (reg%Tr(m)%registry_diags) then
    tr => reg%Tr(m)
    if (tr%id_tr > 0) call post_data(tr%id_tr, tr%t, diag)
    if (tr%id_tendency > 0) then
      work3d(:,:,:) = 0.0
      do k=1,nz ; do j=js,je ; do i=is,ie
        work3d(i,j,k)    = (tr%t(i,j,k) - tr%t_prev(i,j,k))*idt
        tr%t_prev(i,j,k) =  tr%t(i,j,k)
      enddo ; enddo ; enddo
      call post_data(tr%id_tendency, work3d, diag, alt_h=diag_prev%h_state)
    endif
    if ((tr%id_trxh_tendency > 0) .or. (tr%id_trxh_tendency_2d > 0)) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        work3d(i,j,k)     = (tr%t(i,j,k)*h(i,j,k) - tr%Trxh_prev(i,j,k)) * idt
        tr%Trxh_prev(i,j,k) =  tr%t(i,j,k) * h(i,j,k)
      enddo ; enddo ; enddo
      if (tr%id_trxh_tendency > 0) call post_data(tr%id_trxh_tendency, work3d, diag, alt_h=diag_prev%h_state)
      if (tr%id_trxh_tendency_2d > 0) then
        work2d(:,:) = 0.0
        do k=1,nz ; do j=js,je ; do i=is,ie
          work2d(i,j) = work2d(i,j) + work3d(i,j,k)
        enddo ; enddo ; enddo
        call post_data(tr%id_trxh_tendency_2d, work2d, diag)
      endif
    endif
  endif ; enddo
  call diag_restore_grids(diag)

end subroutine post_tracer_diagnostics_at_sync

!> Post the advective and diffusive tendencies
subroutine post_tracer_transport_diagnostics(G, GV, Reg, h_diag, diag)
  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(tracer_registry_type), pointer    :: reg  !< pointer to the tracer registry
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                              intent(in) :: h_diag !< Layer thicknesses on which to post fields
  type(diag_ctrl),            intent(in) :: diag !< structure to regulate diagnostic output

  integer :: i, j, k, is, ie, js, je, nz, m
  real    :: work2d(szi_(g),szj_(g))
  type(tracer_type), pointer :: tr=>null()

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke

  do m=1,reg%ntr ; if (reg%Tr(m)%registry_diags) then
    tr => reg%Tr(m)
    if (tr%id_tr_post_horzn> 0) call post_data(tr%id_tr_post_horzn, tr%t, diag)
    if (tr%id_adx > 0) call post_data(tr%id_adx, tr%ad_x, diag, alt_h=h_diag)
    if (tr%id_ady > 0) call post_data(tr%id_ady, tr%ad_y, diag, alt_h=h_diag)
    if (tr%id_dfx > 0) call post_data(tr%id_dfx, tr%df_x, diag, alt_h=h_diag)
    if (tr%id_dfy > 0) call post_data(tr%id_dfy, tr%df_y, diag, alt_h=h_diag)
    if (tr%id_adx_2d > 0) call post_data(tr%id_adx_2d, tr%ad2d_x, diag)
    if (tr%id_ady_2d > 0) call post_data(tr%id_ady_2d, tr%ad2d_y, diag)
    if (tr%id_dfx_2d > 0) call post_data(tr%id_dfx_2d, tr%df2d_x, diag)
    if (tr%id_dfy_2d > 0) call post_data(tr%id_dfy_2d, tr%df2d_y, diag)
    if (tr%id_adv_xy > 0) call post_data(tr%id_adv_xy, tr%advection_xy, diag, alt_h=h_diag)
    if (tr%id_adv_xy_2d > 0) then
      work2d(:,:) = 0.0
      do k=1,nz ; do j=js,je ; do i=is,ie
        work2d(i,j) = work2d(i,j) + tr%advection_xy(i,j,k)
      enddo ; enddo ; enddo
      call post_data(tr%id_adv_xy_2d, work2d, diag)
    endif
  endif ; enddo

end subroutine post_tracer_transport_diagnostics

!> This subroutine writes out chksums for tracers.
subroutine mom_tracer_chksum(mesg, Tr, ntr, G)
  character(len=*),         intent(in) :: mesg   !< message that appears on the chksum lines
  type(tracer_type),        intent(in) :: tr(:)  !< array of all of registered tracers
  integer,                  intent(in) :: ntr    !< number of registered tracers
  type(ocean_grid_type),    intent(in) :: g      !< ocean grid structure

  integer :: is, ie, js, je, nz
  integer :: i, j, k, m

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  do m=1,ntr
    call hchksum(tr(m)%t, mesg//trim(tr(m)%name), g%HI)
  enddo

end subroutine mom_tracer_chksum

!> Calculates and prints the global inventory of all tracers in the registry.
subroutine mom_tracer_chkinv(mesg, G, h, Tr, ntr)
  character(len=*),                         intent(in) :: mesg !< message that appears on the chksum lines
  type(ocean_grid_type),                    intent(in) :: g    !< ocean grid structure
  type(tracer_type), dimension(:),          intent(in) :: tr   !< array of all of registered tracers
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), intent(in) :: h    !< Layer thicknesses
  integer,                                  intent(in) :: ntr  !< number of registered tracers

  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: tr_inv !< Tracer inventory
  real :: total_inv
  integer :: is, ie, js, je, nz
  integer :: i, j, k, m

  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
  do m=1,ntr
    do k=1,nz ; do j=js,je ; do i=is,ie
      tr_inv(i,j,k) = tr(m)%t(i,j,k)*h(i,j,k)*g%US%L_to_m**2*g%areaT(i,j)*g%mask2dT(i,j)
    enddo ; enddo ; enddo
    total_inv = reproducing_sum(tr_inv, is+(1-g%isd), ie+(1-g%isd), js+(1-g%jsd), je+(1-g%jsd))
    if (is_root_pe()) write(0,'(A,1X,A5,1X,ES25.16,1X,A)') "h-point: inventory", tr(m)%name, total_inv, mesg
  enddo

end subroutine mom_tracer_chkinv

!> Find a tracer in the tracer registry by name.
subroutine tracer_name_lookup(Reg, tr_ptr, name)
  type(tracer_registry_type), pointer    :: reg     !< pointer to tracer registry
  type(tracer_type), pointer             :: tr_ptr  !< target or pointer to the tracer array
  character(len=32), intent(in)          :: name    !< tracer name

  integer n
  do n=1,reg%ntr
    if (lowercase(reg%Tr(n)%name) == lowercase(name)) tr_ptr => reg%Tr(n)
  enddo

end subroutine tracer_name_lookup

!> Initialize the tracer registry.
subroutine tracer_registry_init(param_file, Reg)
  type(param_file_type),      intent(in) :: param_file !< open file to parse for model parameters
  type(tracer_registry_type), pointer    :: reg        !< pointer to tracer registry

  integer, save :: init_calls = 0

! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_tracer_registry" ! This module's name.
  character(len=256) :: mesg    ! Message for error messages.

  if (.not.associated(reg)) then ; allocate(reg)
  else ; return ; endif

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "", all_default=.true.)

  init_calls = init_calls + 1
  if (init_calls > 1) then
    write(mesg,'("tracer_registry_init called ",I3, &
      &" times with different registry pointers.")') init_calls
    if (is_root_pe()) call mom_error(warning,"MOM_tracer"//mesg)
  endif

end subroutine tracer_registry_init


!> This routine closes the tracer registry module.
subroutine tracer_registry_end(Reg)
  type(tracer_registry_type), pointer :: reg !< The tracer registry that will be deallocated
  if (associated(reg)) deallocate(reg)
end subroutine tracer_registry_end

end module mom_tracer_registry