MOM_generic_tracer.F90

!> Drives the generic version of tracers TOPAZ and CFC and other GFDL BGC components
module MOM_generic_tracer

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

#include <MOM_memory.h>

! The following macro is usually defined in <fms_platform.h> but since MOM6 should not directly
! include files from FMS we replicate the macro lines here:
#ifdef NO_F2000
#define _ALLOCATED associated
#else
#define _ALLOCATED allocated
#endif

  ! ### These imports should not reach into FMS directly ###
  use field_manager_mod, only: fm_string_len

  use generic_tracer, only: generic_tracer_register, generic_tracer_get_diag_list
  use generic_tracer, only: generic_tracer_init, generic_tracer_source, generic_tracer_register_diag
  use generic_tracer, only: generic_tracer_coupler_get, generic_tracer_coupler_set
  use generic_tracer, only: generic_tracer_end, generic_tracer_get_list, do_generic_tracer
  use generic_tracer, only: generic_tracer_update_from_bottom,generic_tracer_vertdiff_G
  use generic_tracer, only: generic_tracer_coupler_accumulate

  use g_tracer_utils,   only: g_tracer_get_name,g_tracer_set_values,g_tracer_set_common,g_tracer_get_common
  use g_tracer_utils,   only: g_tracer_get_next,g_tracer_type,g_tracer_is_prog,g_tracer_flux_init
  use g_tracer_utils,   only: g_tracer_send_diag,g_tracer_get_values
  use g_tracer_utils,   only: g_tracer_get_pointer,g_tracer_get_alias,g_tracer_set_csdiag
  use g_tracer_utils,   only: g_tracer_get_obc_segment_props

  use MOM_ALE_sponge, only : set_up_ALE_sponge_field, ALE_sponge_CS
  use MOM_coms, only : EFP_type, max_across_PEs, min_across_PEs, PE_here
  use MOM_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
  use MOM_diag_mediator, only : diag_ctrl, get_diag_time_end
  use MOM_error_handler, only : MOM_error, FATAL, WARNING, NOTE, is_root_pe
  use MOM_file_parser, only : get_param, log_param, log_version, param_file_type
  use MOM_forcing_type, only : forcing, optics_type
  use MOM_grid, only : ocean_grid_type
  use MOM_hor_index, only : hor_index_type
  use MOM_interface_heights, only : thickness_to_dz
  use MOM_io, only : file_exists, MOM_read_data, slasher
  use MOM_open_boundary, only : ocean_OBC_type
  use MOM_open_boundary, only : register_obgc_segments, fill_obgc_segments
  use MOM_open_boundary, only : set_obgc_segments_props
  use MOM_restart, only : register_restart_field, query_initialized, set_initialized, MOM_restart_CS
  use MOM_spatial_means, only : global_area_mean, global_mass_int_EFP, array_global_min_max
  use MOM_sponge, only : set_up_sponge_field, sponge_CS
  use MOM_time_manager, only : time_type, set_time
  use MOM_tracer_diabatic, only : tracer_vertdiff, applyTracerBoundaryFluxesInOut
  use MOM_tracer_registry, only : register_tracer, tracer_registry_type
  use MOM_tracer_Z_init, only : tracer_Z_init
  use MOM_tracer_initialization_from_Z, only : MOM_initialize_tracer_from_Z
  use MOM_unit_scaling, only : unit_scale_type
  use MOM_variables, only : surface, thermo_var_ptrs
  use MOM_verticalGrid, only : verticalGrid_type


  implicit none ; private

  !> A state hidden in module data that is very much not allowed in MOM6
  ! ### This needs to be fixed
  logical :: g_registered = .false.

  public register_MOM_generic_tracer, initialize_MOM_generic_tracer
  public MOM_generic_tracer_column_physics, MOM_generic_tracer_surface_state
  public end_MOM_generic_tracer, MOM_generic_tracer_get
  public MOM_generic_tracer_stock
  public MOM_generic_flux_init
  public MOM_generic_tracer_min_max
  public MOM_generic_tracer_fluxes_accumulate
  public register_MOM_generic_tracer_segments

  !> Control structure for generic tracers
  type, public :: MOM_generic_tracer_CS ; private
    character(len = 200) :: IC_file !< The file in which the generic tracer initial values can
                                    !! be found, or an empty string for internal initialization.
    logical :: Z_IC_file !< If true, the generic_tracer IC_file is in Z-space.  The default is false.
    real :: tracer_IC_val = 0.0    !< The initial value assigned to tracers, in
                                   !! concentration units [conc]
    real :: tracer_land_val = -1.0 !< The values of tracers used where land is masked out, in
                                   !! concentration units [conc]
    logical :: tracers_may_reinit  !< If true, tracers may go through the
                                   !! initialization code if they are not found in the restart files.

    type(diag_ctrl), pointer :: diag => NULL() !< A structure that is used to
                                               !! regulate the timing of diagnostic output.
    type(MOM_restart_CS), pointer :: restart_CSp => NULL() !< Restart control structure
    type(ocean_OBC_type), pointer :: OBC => NULL() !<open boundary condition type
    !> Pointer to the first element of the linked list of generic tracers.
    type(g_tracer_type), pointer :: g_tracer_list => NULL()

  end type MOM_generic_tracer_CS

contains

  !> Initializes the generic tracer packages and adds their tracers to the list
  !! Adds the tracers in the list of generic tracers to the set of MOM tracers (i.e., MOM-register them)
  !! Register these tracers for restart
  function register_MOM_generic_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
    type(hor_index_type),       intent(in)   :: HI         !< Horizontal index ranges
    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(MOM_generic_tracer_CS), pointer     :: CS         !< Pointer to the control structure for this module
    type(tracer_registry_type), pointer      :: tr_Reg     !< Pointer to the control structure for the tracer
                                                           !! advection and diffusion module.
    type(MOM_restart_CS), target, intent(inout)  :: restart_CS !< MOM restart control struct

    ! Local variables
    logical :: register_MOM_generic_tracer
    logical :: obc_has
    ! This include declares and sets the variable "version".
#   include "version_variable.h"

    character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer'
    character(len=200) :: inputdir ! The directory where NetCDF input files are.
    ! These can be overridden later in via the field manager?

    integer :: ntau, axes(3)
    type(g_tracer_type), pointer      :: g_tracer, g_tracer_next
    character(len=fm_string_len)      :: g_tracer_name, longname,units
    character(len=fm_string_len)      :: obc_src_file_name, obc_src_field_name
    real :: lfac_in   ! Multiplicative factor used in setting the tracer-specific inverse length
                      ! scales associated with inflowing tracer reservoirs at OBCs [nondim]
    real :: lfac_out  ! Multiplicative factor used in setting the tracer-specific inverse length
                      ! scales associated with outflowing tracer reservoirs at OBCs [nondim]
    real, dimension(:,:,:,:), pointer   :: tr_field ! A pointer to a generic tracer field, in concentration units [conc]
    real, dimension(:,:,:), pointer     :: tr_ptr   ! A pointer to a generic tracer field, in concentration units [conc]
    real,    dimension(SZI_(HI),SZJ_(HI),SZK_(GV)) :: grid_tmask ! A 3-d copy of G%mask2dT [nondim]
    integer, dimension(SZI_(HI),SZJ_(HI))          :: grid_kmt   ! A 2-d array of nk

    register_MOM_generic_tracer = .false.
    if (associated(CS)) then
      call MOM_error(FATAL, "register_MOM_generic_tracer called with an "// &
                            "associated control structure.")
    endif
    allocate(CS)


    !Register all the generic tracers used and create the list of them.
    !This can be called by ALL PE's. No array fields allocated.
    if (.not. g_registered) then
      call generic_tracer_register()
      g_registered = .true.
    endif


  ! Read all relevant parameters and write them to the model log.
    call log_version(param_file, sub_name, version, "")
    call get_param(param_file, sub_name, "GENERIC_TRACER_IC_FILE", CS%IC_file, &
                 "The file in which the generic tracer initial values can "//&
                 "be found, or an empty string for internal initialization.", &
                 default=" ")
    if ((len_trim(CS%IC_file) > 0) .and. (scan(CS%IC_file,'/') == 0)) then
      ! Add the directory if CS%IC_file is not already a complete path.
      call get_param(param_file, sub_name, "INPUTDIR", inputdir, default=".")
      CS%IC_file = trim(slasher(inputdir))//trim(CS%IC_file)
      call log_param(param_file, sub_name, "INPUTDIR/GENERIC_TRACER_IC_FILE", CS%IC_file)
    endif
    call get_param(param_file, sub_name, "GENERIC_TRACER_IC_FILE_IS_Z", CS%Z_IC_file, &
                 "If true, GENERIC_TRACER_IC_FILE is in depth space, not "//&
                 "layer space.",default=.false.)
    call get_param(param_file, sub_name, "TRACERS_MAY_REINIT", CS%tracers_may_reinit, &
                 "If true, tracers may go through the initialization code "//&
                 "if they are not found in the restart files.  Otherwise "//&
                 "it is a fatal error if tracers are not found in the "//&
                 "restart files of a restarted run.", default=.false.)

    CS%restart_CSp => restart_CS

    ntau=1 ! MOM needs the fields at only one time step


    !   At this point G%mask2dT and CS%diag%axesTL are not allocated.
    ! postpone diag_registeration to initialize_MOM_generic_tracer

    !Fields cannot be diag registered as they are allocated and have to registered later.
    grid_tmask(:,:,:) = 0.0
    grid_kmt(:,:) = 0
    axes(:) = -1

    !
    ! Initialize all generic tracers
    !
    call generic_tracer_init(HI%isc,HI%iec,HI%jsc,HI%jec,HI%isd,HI%ied,HI%jsd,HI%jed,&
         GV%ke,ntau,axes,grid_tmask,grid_kmt,set_time(0,0))


    !
    ! MOM-register the generic tracers
    !

    !Get the tracer list
    call generic_tracer_get_list(CS%g_tracer_list)
    if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//&
         ": No tracer in the list.")
    ! For each tracer name get its T_prog index and get its fields

    g_tracer=>CS%g_tracer_list
    do
      call g_tracer_get_alias(g_tracer,g_tracer_name)

      call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',tr_field)
      call g_tracer_get_values(g_tracer,g_tracer_name,'longname', longname)
      call g_tracer_get_values(g_tracer,g_tracer_name,'units',units )

      !!nnz: MOM field is 3D. Does this affect performance? Need it be override field?
      tr_ptr => tr_field(:,:,:,1)
      ! Register prognostic tracer for horizontal advection, diffusion, and restarts.
      if (g_tracer_is_prog(g_tracer)) then
        call register_tracer(tr_ptr, tr_Reg, param_file, HI, GV, &
                             name=g_tracer_name, longname=longname, units=units, &
                             registry_diags=.false., &   !### CHANGE TO TRUE?
                             restart_CS=restart_CS, mandatory=.not.CS%tracers_may_reinit)
      else
        call register_restart_field(tr_ptr, g_tracer_name, .not.CS%tracers_may_reinit, &
                                    restart_CS, longname=longname, units=units)
      endif

      !traverse the linked list till hit NULL
      call g_tracer_get_next(g_tracer, g_tracer_next)
      if (.NOT. associated(g_tracer_next)) exit
      g_tracer=>g_tracer_next

    enddo

    register_MOM_generic_tracer = .true.
  end function register_MOM_generic_tracer

  !> Register OBC segments for generic tracers
  subroutine register_MOM_generic_tracer_segments(CS, GV, OBC, tr_Reg, param_file)
    type(MOM_generic_tracer_CS), pointer    :: CS         !< Pointer to the control structure for this module.
    type(verticalGrid_type),     intent(in) :: GV         !< The ocean's vertical grid structure
    type(ocean_OBC_type),        pointer    :: OBC        !< This open boundary condition type specifies whether,
                                                          !! where, and what open boundary conditions are used.
    type(tracer_registry_type),  pointer    :: tr_Reg     !< Pointer to the control structure for the tracer
                                                          !! advection and diffusion module.
    type(param_file_type),       intent(in) :: param_file !< A structure to parse for run-time parameters

    ! Local variables
    logical :: obc_has
    ! This include declares and sets the variable "version".
#   include "version_variable.h"
    character(len=128), parameter :: sub_name = 'register_MOM_generic_tracer_segments'
    type(g_tracer_type), pointer      :: g_tracer,g_tracer_next
    character(len=fm_string_len)      :: g_tracer_name
    character(len=fm_string_len)      :: obc_src_file_name, obc_src_field_name
    real :: lfac_in   ! Multiplicative factor used in setting the tracer-specific inverse length
                      ! scales associated with inflowing tracer reservoirs at OBCs [nondim]
    real :: lfac_out  ! Multiplicative factor used in setting the tracer-specific inverse length
                      ! scales associated with outflowing tracer reservoirs at OBCs [nondim]

    if (.NOT. associated(OBC)) return
    !Get the tracer list
    call generic_tracer_get_list(CS%g_tracer_list)
    if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//&
         ": No tracer in the list.")

    g_tracer=>CS%g_tracer_list
    do
      call g_tracer_get_alias(g_tracer,g_tracer_name)
      if (g_tracer_is_prog(g_tracer)) then
          call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has ,&
                                   obc_src_file_name,obc_src_field_name,lfac_in,lfac_out)
          if (obc_has) then
            call set_obgc_segments_props(OBC,g_tracer_name,obc_src_file_name,obc_src_field_name,lfac_in,lfac_out)
            call register_obgc_segments(GV, OBC, tr_Reg, param_file, g_tracer_name)
          endif
      endif

      !traverse the linked list till hit NULL
      call g_tracer_get_next(g_tracer, g_tracer_next)
      if (.NOT. associated(g_tracer_next)) exit
      g_tracer=>g_tracer_next

    enddo

  end subroutine register_MOM_generic_tracer_segments

  !>  Initialize phase II:  Initialize required variables for generic tracers
  !!  There are some steps of initialization that cannot be done in register_MOM_generic_tracer
  !!  This is the place and time to do them:
  !!      Set the grid mask and initial time for all generic tracers.
  !!      Diag_register them.
  !!      Z_diag_register them.
  !!
  !!   This subroutine initializes the NTR tracer fields in tr(:,:,:,:)
  !! and it sets up the tracer output.
  subroutine initialize_MOM_generic_tracer(restart, day, G, GV, US, h, tv, param_file, diag, OBC, &
                                           CS, sponge_CSp, ALE_sponge_CSp)
    logical,                               intent(in) :: restart !< .true. if the fields have already been
                                                                 !! read from a restart file.
    type(time_type), target,               intent(in) :: day     !< Time of the start of the run.
    type(ocean_grid_type),                 intent(inout) :: G    !< The ocean's grid structure
    type(verticalGrid_type),               intent(in)    :: GV   !< The ocean's vertical grid structure
    type(unit_scale_type),                 intent(in)    :: US   !< A dimensional unit scaling type
    real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h   !< Layer thicknesses [H ~> m or kg m-2]
    type(thermo_var_ptrs),                 intent(in) :: tv      !< A structure pointing to various thermodynamic
                                                                 !! variables
    type(param_file_type),                 intent(in) :: param_file !< A structure to parse for run-time parameters
    type(diag_ctrl),               target, intent(in) :: diag    !< Regulates diagnostic output.
    type(ocean_OBC_type),                  pointer    :: OBC     !< This open boundary condition type specifies whether,
                                                                 !! where, and what open boundary conditions are used.
    type(MOM_generic_tracer_CS),           pointer    :: CS      !< Pointer to the control structure for this module.
    type(sponge_CS),                       pointer    :: sponge_CSp !< Pointer to the control structure for the sponges.
    type(ALE_sponge_CS),                   pointer    :: ALE_sponge_CSp !< Pointer  to the control structure for the
                                                                 !! ALE sponges.

    character(len=128), parameter :: sub_name = 'initialize_MOM_generic_tracer'
    logical :: OK,obc_has
    integer :: i, j, k, isc, iec, jsc, jec, nk
    type(g_tracer_type), pointer    :: g_tracer,g_tracer_next
    character(len=fm_string_len)      :: g_tracer_name
    real, dimension(:,:,:,:), pointer   :: tr_field ! A pointer to a generic tracer field, in concentration units [conc]
    real, dimension(:,:,:), pointer     :: tr_ptr   ! A pointer to a generic tracer field, in concentration units [conc]
    real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dz ! Layer vertical extent [Z ~> m]
    real,    dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: grid_tmask ! A 3-d copy of G%mask2dT [nondim]
    integer, dimension(SZI_(G),SZJ_(G))          :: grid_kmt   ! A 2-d array of nk

    !! 2010/02/04  Add code to re-initialize Generic Tracers if needed during a model simulation
    !! By default, restart cpio should not contain a Generic Tracer IC file and step below will be skipped.
    !! Ideally, the generic tracer IC file should have the tracers on Z levels.

    isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke

    CS%diag=>diag
    !Get the tracer list
    if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//&
         ": No tracer in the list.")
    !For each tracer name get its  fields
    g_tracer=>CS%g_tracer_list

    call thickness_to_dz(h, tv, dz, G, GV, US)

    do
      if (INDEX(CS%IC_file, '_NULL_') /= 0) then
        call MOM_error(WARNING, "The name of the IC_file "//trim(CS%IC_file)//&
                              " indicates no MOM initialization was asked for the generic tracers."//&
                              "Bypassing the MOM initialization of ALL generic tracers!")
        exit
      endif
      call g_tracer_get_alias(g_tracer,g_tracer_name)
      call g_tracer_get_pointer(g_tracer,g_tracer_name,'field',tr_field)
      tr_ptr => tr_field(:,:,:,1)

      if (.not.restart .or. (CS%tracers_may_reinit .and. &
          .not.query_initialized(tr_ptr, g_tracer_name, CS%restart_CSp))) then

        if (g_tracer%requires_src_info ) then
          call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//&
                              "initializing generic tracer "//trim(g_tracer_name)//&
                              " using MOM_initialize_tracer_from_Z ")

          call MOM_initialize_tracer_from_Z(dz, tr_ptr, G, GV, US, param_file, &
                                 src_file=g_tracer%src_file, src_var_nam=g_tracer%src_var_name, &
                                 src_var_unit_conversion=g_tracer%src_var_unit_conversion, &
                                 src_var_record=g_tracer%src_var_record, src_var_gridspec=g_tracer%src_var_gridspec, &
                                 h_in_Z_units=.true.)

          !Check/apply the bounds for each g_tracer
          do k=1,nk ; do j=jsc,jec ; do i=isc,iec
            if (tr_ptr(i,j,k) /= CS%tracer_land_val) then
              if (tr_ptr(i,j,k) < g_tracer%src_var_valid_min) tr_ptr(i,j,k) = g_tracer%src_var_valid_min
              !Jasmin does not want to apply the maximum for now
              !if (tr_ptr(i,j,k) > g_tracer%src_var_valid_max) tr_ptr(i,j,k) = g_tracer%src_var_valid_max
            endif
          enddo ; enddo ; enddo

          !jgj: Reset CASED to 0 below K=1
          if ( (trim(g_tracer_name) == 'cased') .or. (trim(g_tracer_name) == 'ca13csed') ) then
            do k=2,nk ; do j=jsc,jec ; do i=isc,iec
              if (tr_ptr(i,j,k) /= CS%tracer_land_val) then
                tr_ptr(i,j,k) = 0.0
              endif
            enddo ; enddo ; enddo
          endif
        elseif(.not. g_tracer%requires_restart) then
         !Do nothing for this tracer, it is initialized by the tracer package
          call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//&
                            "skip initialization of generic tracer "//trim(g_tracer_name))
        else !Do it old way if the tracer is not registered to start from a specific source file.
            !This path should be deprecated if all generic tracers are required to start from specified sources.
          if (len_trim(CS%IC_file) > 0) then
          !  Read the tracer concentrations from a netcdf file.
            if (.not.file_exists(CS%IC_file)) call MOM_error(FATAL, &
                    "initialize_MOM_Generic_tracer: Unable to open "//CS%IC_file)
            if (CS%Z_IC_file) then
              OK = tracer_Z_init(tr_ptr, h, CS%IC_file, g_tracer_name, G, GV, US)
              if (.not.OK) then
                OK = tracer_Z_init(tr_ptr, h, CS%IC_file, trim(g_tracer_name), G, GV, US)
                if (.not.OK) call MOM_error(FATAL,"initialize_MOM_Generic_tracer: "//&
                        "Unable to read "//trim(g_tracer_name)//" from "//&
                        trim(CS%IC_file)//".")
              endif
              call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//&
                              "initialized generic tracer "//trim(g_tracer_name)//&
                              " using Generic Tracer File on Z: "//CS%IC_file)
            else
              ! native grid
              call MOM_error(NOTE,"initialize_MOM_generic_tracer: "//&
                    "Using Generic Tracer IC file on native grid "//trim(CS%IC_file)//&
                    " for tracer "//trim(g_tracer_name))
              call MOM_read_data(CS%IC_file, trim(g_tracer_name), tr_ptr, G%Domain)
            endif
          else
            call MOM_error(FATAL,"initialize_MOM_generic_tracer: "//&
                    "check Generic Tracer IC filename "//trim(CS%IC_file)//&
                    " for tracer "//trim(g_tracer_name))
          endif

        endif

        call set_initialized(tr_ptr, g_tracer_name, CS%restart_CSp)
      endif

      call g_tracer_get_obc_segment_props(g_tracer,g_tracer_name,obc_has )
      if(obc_has .and. g_tracer_is_prog(g_tracer) .and. .not.restart) &
             call fill_obgc_segments(G, GV, OBC, tr_ptr, g_tracer_name)
      !traverse the linked list till hit NULL
      call g_tracer_get_next(g_tracer, g_tracer_next)
      if (.NOT. associated(g_tracer_next)) exit
      g_tracer=>g_tracer_next
    enddo
    !! end section to re-initialize generic tracers


    !Now we can reset the grid mask, axes and time to their true values
    !Note that grid_tmask must be set correctly on the data domain boundary
    !so that coast mask can be deduced from it.
    grid_tmask(:,:,:) = 0.0
    grid_kmt(:,:) = 0
    do j = G%jsd, G%jed ; do i = G%isd, G%ied
      if (G%mask2dT(i,j) > 0.0) then
        grid_tmask(i,j,:) = 1.0
        grid_kmt(i,j) = GV%ke ! Tell the code that a layer thicker than 1m is the bottom layer.
      endif
    enddo ; enddo
    call g_tracer_set_common(G%isc,G%iec,G%jsc,G%jec,G%isd,G%ied,G%jsd,G%jed,&
                             GV%ke,1,CS%diag%axesTL%handles,grid_tmask,grid_kmt,day)

    ! Register generic tracer modules diagnostics

#ifdef _USE_MOM6_DIAG
    call g_tracer_set_csdiag(CS%diag)
#endif
    call generic_tracer_register_diag()
#ifdef _USE_MOM6_DIAG
    call g_tracer_set_csdiag(CS%diag)
#endif

  end subroutine initialize_MOM_generic_tracer

  !>  Column physics for generic tracers.
  !!      Get the coupler values for generic tracers that exchange with atmosphere
  !!      Update generic tracer concentration fields from sources and sinks.
  !!      Vertically diffuse generic tracer concentration fields.
  !!      Update generic tracers from bottom and their bottom reservoir.
  !!
  !!   This subroutine applies diapycnal diffusion and any other column
  !! tracer physics or chemistry to the tracers from this file.
  !! CFCs are relatively simple, as they are passive tracers. with only a surface
  !! flux as a source.
  subroutine MOM_generic_tracer_column_physics(h_old, h_new, ea, eb, fluxes, Hml, dt, G, GV, US, CS, tv, optics, &
        evap_CFL_limit, minimum_forcing_depth)
    type(ocean_grid_type),   intent(in) :: G     !< The ocean's grid structure
    type(verticalGrid_type), intent(in) :: GV    !< The ocean's vertical grid structure
    real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                             intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
    real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                             intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
    real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                             intent(in) :: ea    !< The amount of fluid entrained from the layer
                                                 !! above during this call [H ~> m or kg m-2].
    real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), &
                             intent(in) :: eb    !< The amount of fluid entrained from the layer
                                                 !! below during this call [H ~> m or kg m-2].
    type(forcing),           intent(in) :: fluxes !< A structure containing pointers to thermodynamic
                                                 !! and tracer forcing fields.
    real, dimension(SZI_(G),SZJ_(G)), intent(in) :: Hml  !< Mixed layer depth [Z ~> m]
    real,                    intent(in) :: dt    !< The amount of time covered by this call [T ~> s]
    type(unit_scale_type),   intent(in) :: US    !< A dimensional unit scaling type
    type(MOM_generic_tracer_CS), pointer :: CS   !< Pointer to the control structure for this module.
    type(thermo_var_ptrs),   intent(in) :: tv    !< A structure pointing to various thermodynamic variables
    type(optics_type),       intent(in) :: optics !< The structure containing optical properties.
    real,          optional, intent(in) :: evap_CFL_limit !< Limit on the fraction of the water that can
                                                 !! be fluxed out of the top layer in a timestep [nondim]
                                                 !   Stored previously in diabatic CS.
    real,          optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which fluxes
                                                 !!  can be applied [H ~> m or kg m-2]
                                                 !   Stored previously in diabatic CS.
    ! The arguments to this subroutine are redundant in that
    !     h_new(k) = h_old(k) + ea(k) - eb(k-1) + eb(k) - ea(k+1)

    ! Local variables
    character(len=128), parameter :: sub_name = 'MOM_generic_tracer_column_physics'

    type(g_tracer_type), pointer  :: g_tracer, g_tracer_next
    character(len=fm_string_len)  :: g_tracer_name
    real, dimension(:,:), pointer :: stf_array   ! The surface flux of the tracer [conc kg m-2 s-1]
    real, dimension(:,:), pointer :: trunoff_array  ! The tracer concentration in the river runoff [conc]
    real, dimension(:,:), pointer :: runoff_tracer_flux_array ! The runoff tracer flux [conc kg m-2 s-1]

    real :: surface_field(SZI_(G),SZJ_(G))  ! The surface value of some field, here only used for salinity [S ~> ppt]
    real :: dz_ml(SZI_(G),SZJ_(G))  ! The mixed layer depth in the MKS units used for generic tracers [m]
    real :: sosga ! The global mean surface salinity [ppt]

    real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: rho_dzt ! Layer mass per unit area [kg m-2]
    real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: dzt     ! Layer vertical extents [m]
    real, dimension(SZI_(G),SZJ_(G),SZK_(GV)) :: h_work  ! A work array of thicknesses [H ~> m or kg m-2]
    integer :: i, j, k, isc, iec, jsc, jec, nk

    isc = G%isc ; iec = G%iec ; jsc = G%jsc ; jec = G%jec ; nk = GV%ke

    !Get the tracer list
    if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL,&
         trim(sub_name)//": No tracer in the list.")

#ifdef _USE_MOM6_DIAG
    call g_tracer_set_csdiag(CS%diag)
#endif

    !
    !Extract the tracer surface fields from coupler and update tracer fields from sources
    !
    !call generic_tracer_coupler_get(fluxes%tr_fluxes)
    !Niki: This is moved out to ocean_model_MOM.F90 because if dt_therm>dt_cpld we need to average
    !      the fluxes without coming into this subroutine.
    !      MOM5 has to modified to conform.

    !
    !Add contribution of river to surface flux
    !
    g_tracer=>CS%g_tracer_list
    do
      if (_ALLOCATED(g_tracer%trunoff) .and. (.NOT. g_tracer%runoff_added_to_stf)) then
        call g_tracer_get_alias(g_tracer,g_tracer_name)
        call g_tracer_get_pointer(g_tracer,g_tracer_name,'stf',   stf_array)
        call g_tracer_get_pointer(g_tracer,g_tracer_name,'trunoff',trunoff_array)
        call g_tracer_get_pointer(g_tracer,g_tracer_name,'runoff_tracer_flux',runoff_tracer_flux_array)
        !nnz: Why is fluxes%river = 0?
        runoff_tracer_flux_array(:,:) = trunoff_array(:,:) * &
                 US%RZ_T_to_kg_m2s*fluxes%lrunoff(:,:)
        stf_array = stf_array + runoff_tracer_flux_array
        g_tracer%runoff_added_to_stf = .true.
      endif

      !traverse the linked list till hit NULL
      call g_tracer_get_next(g_tracer, g_tracer_next)
      if (.NOT. associated(g_tracer_next)) exit
      g_tracer => g_tracer_next

    enddo

    !
    !Prepare input arrays for source update
    !

    rho_dzt(:,:,:) = GV%H_to_kg_m2 * GV%Angstrom_H
    do k=1,nk ; do j=jsc,jec ; do i=isc,iec
      rho_dzt(i,j,k) = GV%H_to_kg_m2 * h_old(i,j,k)
    enddo ; enddo ; enddo

    dzt(:,:,:) = 1.0
    call thickness_to_dz(h_old, tv, dzt, G, GV, US)
    do k=1,nk ; do j=jsc,jec ; do i=isc,iec
      dzt(i,j,k) = US%Z_to_m * dzt(i,j,k)
    enddo ; enddo ; enddo
    dz_ml(:,:) = 0.0
    do j=jsc,jec ; do i=isc,iec
      surface_field(i,j) = tv%S(i,j,1)
      dz_ml(i,j) = US%Z_to_m * Hml(i,j)
    enddo ; enddo
    sosga = global_area_mean(surface_field, G, unscale=US%S_to_ppt)

    !
    !Calculate tendencies (i.e., field changes at dt) from the sources / sinks
    !
    if ((G%US%L_to_m == 1.0) .and. (G%US%s_to_T == 1.0) .and. (G%US%Z_to_m == 1.0) .and. &
        (G%US%Q_to_J_kg == 1.0) .and. (G%US%RZ_to_kg_m2 == 1.0) .and. &
        (US%C_to_degC == 1.0) .and. (US%S_to_ppt == 1.0)) then
      ! Avoid unnecessary copies when no unit conversion is needed.
      call generic_tracer_source(tv%T, tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, &
               G%areaT, get_diag_time_end(CS%diag), &
               optics%nbands, optics%max_wavelength_band, optics%sw_pen_band, optics%opacity_band, &
               internal_heat=tv%internal_heat, frunoff=fluxes%frunoff, sosga=sosga)
    else
      ! tv%internal_heat is a null pointer unless DO_GEOTHERMAL = True,
      ! so we have to check and only do the scaling if it is associated.
      if(associated(tv%internal_heat)) then
        call generic_tracer_source(US%C_to_degC*tv%T, US%S_to_ppt*tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, &
               G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), &
               optics%nbands, optics%max_wavelength_band, &
               sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), &
               opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), &
               internal_heat=G%US%RZ_to_kg_m2*US%C_to_degC*tv%internal_heat(:,:), &
               frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga)
      else
        call generic_tracer_source(US%C_to_degC*tv%T, US%S_to_ppt*tv%S, rho_dzt, dzt, dz_ml, G%isd, G%jsd, 1, dt, &
               G%US%L_to_m**2*G%areaT(:,:), get_diag_time_end(CS%diag), &
               optics%nbands, optics%max_wavelength_band, &
               sw_pen_band=G%US%QRZ_T_to_W_m2*optics%sw_pen_band(:,:,:), &
               opacity_band=G%US%m_to_Z*optics%opacity_band(:,:,:,:), &
               frunoff=G%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga)
      endif
    endif

    ! This uses applyTracerBoundaryFluxesInOut to handle the change in tracer due to freshwater fluxes
    ! usually in ALE mode
    if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
      g_tracer=>CS%g_tracer_list
      do
        if (g_tracer_is_prog(g_tracer)) then
          do k=1,nk ;do j=jsc,jec ; do i=isc,iec
            h_work(i,j,k) = h_old(i,j,k)
          enddo ; enddo ; enddo
          call applyTracerBoundaryFluxesInOut(G, GV, g_tracer%field(:,:,:,1), dt, &
                            fluxes, h_work, evap_CFL_limit, minimum_forcing_depth)
        endif

         !traverse the linked list till hit NULL
         call g_tracer_get_next(g_tracer, g_tracer_next)
        if (.NOT. associated(g_tracer_next)) exit
        g_tracer=>g_tracer_next
      enddo
    endif

    !
    !Update Tr(n)%field from explicit vertical diffusion
    !
    ! Use a tridiagonal solver to determine the concentrations after the
    ! surface source is applied and diapycnal advection and diffusion occurs.
    if (present(evap_CFL_limit) .and. present(minimum_forcing_depth)) then
      ! Last arg is tau which is always 1 for MOM6
      call generic_tracer_vertdiff_G(h_work, ea, eb, US%T_to_s*dt, GV%kg_m2_to_H, GV%m_to_H, 1)
    else
      ! Last arg is tau which is always 1 for MOM6
      call generic_tracer_vertdiff_G(h_old, ea, eb, US%T_to_s*dt, GV%kg_m2_to_H, GV%m_to_H, 1)
    endif

    ! Update bottom fields after vertical processes

    ! Second arg is tau which is always 1 for MOM6
    call generic_tracer_update_from_bottom(US%T_to_s*dt, 1, get_diag_time_end(CS%diag))

    !Output diagnostics via diag_manager for all generic tracers and their fluxes
    call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1)
#ifdef _USE_MOM6_DIAG
    call g_tracer_set_csdiag(CS%diag)
#endif

  end subroutine MOM_generic_tracer_column_physics

  !> This subroutine calculates mass-weighted integral on the PE either
  !! of all available tracer concentrations, or of a tracer that is
  !! being requested specifically, returning the number of stocks it has
  !! calculated. If the stock_index is present, only the stock corresponding
  !! to that coded index is returned.
  function MOM_generic_tracer_stock(h, stocks, G, GV, CS, names, units, stock_index)
    type(ocean_grid_type),              intent(in)    :: G    !< The ocean's grid structure
    type(verticalGrid_type),            intent(in)    :: GV   !< The ocean's vertical grid structure
    real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
    type(EFP_type), dimension(:),       intent(out)   :: stocks !< The mass-weighted integrated amount of each
                                                                !! tracer, in kg times concentration units [kg conc]
    type(MOM_generic_tracer_CS),        pointer       :: CS     !< Pointer to the control structure for this module.
    character(len=*), dimension(:),     intent(out)   :: names  !< The names of the stocks calculated.
    character(len=*), dimension(:),     intent(out)   :: units  !< The units of the stocks calculated.
    integer, optional,                  intent(in)    :: stock_index !< The coded index of a specific stock
                                                                     !! being sought.
    integer                                           :: MOM_generic_tracer_stock !< Return value, the
                                                                     !! number of stocks calculated here.

    ! Local variables
    type(g_tracer_type), pointer  :: g_tracer, g_tracer_next
    real, dimension(:,:,:,:), pointer   :: tr_field ! A pointer to a generic tracer field, in concentration units [conc]
    real, dimension(:,:,:), pointer     :: tr_ptr   ! A pointer to a generic tracer field, in concentration units [conc]
    character(len=128), parameter :: sub_name = 'MOM_generic_tracer_stock'

    integer :: m

    MOM_generic_tracer_stock = 0
    if (.not.associated(CS)) return

    if (present(stock_index)) then ; if (stock_index > 0) then
      ! Check whether this stock is available from this routine.

      ! No stocks from this routine are being checked yet.  Return 0.
      return
    endif ; endif

    if (.NOT. associated(CS%g_tracer_list)) return ! No stocks.

    m=1 ; g_tracer=>CS%g_tracer_list
    do
      call g_tracer_get_alias(g_tracer,names(m))
      call g_tracer_get_values(g_tracer,names(m),'units',units(m))
      units(m) = trim(units(m))//" kg"
      call g_tracer_get_pointer(g_tracer,names(m),'field',tr_field)

      tr_ptr => tr_field(:,:,:,1)
      stocks(m) = global_mass_int_EFP(h, G, GV, tr_ptr, on_PE_only=.true.)

      !traverse the linked list till hit NULL
      call g_tracer_get_next(g_tracer, g_tracer_next)
      if (.NOT. associated(g_tracer_next)) exit
      g_tracer=>g_tracer_next
      m = m+1
    enddo

    MOM_generic_tracer_stock = m

  end function MOM_generic_tracer_stock

  !> This subroutine finds the global min and max of either of all available
  !! tracer concentrations, or of a tracer that is being requested specifically,
  !! returning the number of tracers it has evaluated.
  !! It also optionally returns the locations of the extrema.
  function MOM_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, G, CS, names, units, &
                                      xgmin, ygmin, zgmin, xgmax, ygmax, zgmax)
    integer,                        intent(in)    :: ind_start !< The index of the tracer to start with
    logical, dimension(:),          intent(out)   :: got_minmax !< Indicates whether the global min and
                                                            !! max are found for each tracer
    real, dimension(:),             intent(out)   :: gmin   !< Global minimum of each tracer [conc]
    real, dimension(:),             intent(out)   :: gmax   !< Global maximum of each tracer [conc]
    type(ocean_grid_type),          intent(in)    :: G      !< The ocean's grid structure
    type(MOM_generic_tracer_CS),    pointer       :: CS     !< Pointer to the control structure for this module.
    character(len=*), dimension(:), intent(out)   :: names  !< The names of the stocks calculated.
    character(len=*), dimension(:), intent(out)   :: units  !< The units of the stocks calculated.
    real, dimension(:),   optional, intent(out)   :: xgmin  !< The x-position of the global minimum in the
                                                            !! units of G%geoLonT, often [degrees_E] or [km] or [m]
    real, dimension(:),   optional, intent(out)   :: ygmin  !< The y-position of the global minimum in the
                                                            !! units of G%geoLatT, often [degrees_N] or [km] or [m]
    real, dimension(:),   optional, intent(out)   :: zgmin  !< The z-position of the global minimum [layer]
    real, dimension(:),   optional, intent(out)   :: xgmax  !< The x-position of the global maximum in the
                                                            !! units of G%geoLonT, often [degrees_E] or [km] or [m]
    real, dimension(:),   optional, intent(out)   :: ygmax  !< The y-position of the global maximum in the
                                                            !! units of G%geoLatT, often [degrees_N] or [km] or [m]
    real, dimension(:),   optional, intent(out)   :: zgmax  !< The z-position of the global maximum [layer]
    integer                                       :: MOM_generic_tracer_min_max !< Return value, the
                                                            !! number of tracers done here.

    ! Local variables
    type(g_tracer_type), pointer  :: g_tracer, g_tracer_next
    real, dimension(:,:,:,:), pointer   :: tr_field ! The tracer array whose extrema are being sought [conc]
    real, dimension(:,:,:), pointer     :: tr_ptr   ! The tracer array whose extrema are being sought [conc]
    real :: x_min  ! The x-position of the global minimum in the units of G%geoLonT, often [degrees_E] or [km] or [m]
    real :: y_min  ! The y-position of the global minimum in the units of G%geoLatT, often [degrees_N] or [km] or [m]
    real :: z_min  ! The z-position of the global minimum [layer]
    real :: x_max  ! The x-position of the global maximum in the units of G%geoLonT, often [degrees_E] or [km] or [m]
    real :: y_max  ! The y-position of the global maximum in the units of G%geoLatT, often [degrees_N] or [km] or [m]
    real :: z_max  ! The z-position of the global maximum [layer]
    character(len=128), parameter :: sub_name = 'MOM_generic_tracer_min_max'

    logical :: find_location
    integer :: isc, iec, jsc, jec, isd, ied, jsd, jed, nk, ntau
    integer :: k, is, ie, js, je, m

    is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec

    MOM_generic_tracer_min_max = 0
    if (.not.associated(CS)) return

    if (.NOT. associated(CS%g_tracer_list)) return ! No stocks.

    call g_tracer_get_common(isc, iec, jsc, jec, isd, ied, jsd, jed, nk, ntau)
    find_location = present(xgmin) .or. present(ygmin) .or. present(zgmin) .or. &
                    present(xgmax) .or. present(ygmax) .or. present(zgmax)

    m=ind_start ; g_tracer=>CS%g_tracer_list
    do
      call g_tracer_get_alias(g_tracer,names(m))
      call g_tracer_get_values(g_tracer,names(m),'units',units(m))
      call g_tracer_get_pointer(g_tracer,names(m),'field',tr_field)

      gmin(m) = -1.0
      gmax(m) = -1.0

      tr_ptr => tr_field(:,:,:,1)

      if (find_location) then
        call array_global_min_max(tr_ptr, G, nk, gmin(m), gmax(m), &
                                  x_min, y_min, z_min, x_max, y_max, z_max)
        if (present(xgmin)) xgmin(m) = x_min
        if (present(ygmin)) ygmin(m) = y_min
        if (present(zgmin)) zgmin(m) = z_min
        if (present(xgmax)) xgmax(m) = x_max
        if (present(ygmax)) ygmax(m) = y_max
        if (present(zgmax)) zgmax(m) = z_max
      else
        call array_global_min_max(tr_ptr, G, nk, gmin(m), gmax(m))
      endif

      got_minmax(m) = .true.

      !traverse the linked list till hit NULL
      call g_tracer_get_next(g_tracer, g_tracer_next)
      if (.NOT. associated(g_tracer_next)) exit
      g_tracer=>g_tracer_next
      m = m+1
    enddo

    MOM_generic_tracer_min_max = m

  end function MOM_generic_tracer_min_max

  !> This subroutine calculates the surface state and sets coupler values for
  !! those generic tracers that have flux exchange with atmosphere.
  !!
  !! This subroutine sets up the fields that the coupler needs to calculate the
  !! CFC fluxes between the ocean and atmosphere.
  subroutine MOM_generic_tracer_surface_state(sfc_state, h, G, GV, CS)
    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(surface),                         intent(inout) :: sfc_state !< A structure containing fields that
                                                                 !! describe the surface state of the ocean.
    real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in) :: h    !< Layer thicknesses [H ~> m or kg m-2]
    type(MOM_generic_tracer_CS),           pointer       :: CS   !< Pointer to the control structure for this module.

    ! Local variables
    real :: sosga ! The global mean surface salinity [ppt]
    real, dimension(SZI_(G),SZJ_(G),SZK_(GV),1) :: rho0 ! An unused array of densities [kg m-3]
    real, dimension(SZI_(G),SZJ_(G),SZK_(GV))   :: dzt  ! Layer vertical extents [m]

    character(len=128), parameter :: sub_name = 'MOM_generic_tracer_surface_state'

    !Set coupler values
    !nnz: fake rho0
    rho0(:,:,:,:) = 1.0

    dzt(:,:,:) = GV%H_to_m * h(:,:,:)

    sosga = global_area_mean(sfc_state%SSS, G, unscale=G%US%S_to_ppt)

    if ((G%US%C_to_degC == 1.0) .and. (G%US%S_to_ppt == 1.0)) then
      call generic_tracer_coupler_set(sfc_state%tr_fields, &
              ST=sfc_state%SST, SS=sfc_state%SSS, &
              rho=rho0, & !nnz: required for MOM5 and previous versions.
              ilb=G%isd, jlb=G%jsd, &
              dzt=dzt,& !This is needed for the Mocsy method of carbonate system vars
              tau=1, sosga=sosga, model_time=get_diag_time_end(CS%diag))
    else
      call generic_tracer_coupler_set(sfc_state%tr_fields, &
              ST=G%US%C_to_degC*sfc_state%SST, SS=G%US%S_to_ppt*sfc_state%SSS, &
              rho=rho0, & !nnz: required for MOM5 and previous versions.
              ilb=G%isd, jlb=G%jsd, &
              dzt=dzt,& !This is needed for the Mocsy method of carbonate system vars
              tau=1, sosga=sosga, model_time=get_diag_time_end(CS%diag))
    endif

    !Output diagnostics via diag_manager for all tracers in this module
!    if (.NOT. associated(CS%g_tracer_list)) call MOM_error(FATAL, trim(sub_name)//&
!         "No tracer in the list.")
!    call g_tracer_send_diag(CS%g_tracer_list, get_diag_time_end(CS%diag), tau=1)
    !Niki: The problem with calling diagnostic outputs here is that this subroutine is called every dt_cpld
    !      hence if dt_therm > dt_cpld we get output (and contribution to the mean) at times that tracers
    !      had not been updated.
    !      Moving this to the end of column physics subroutine fixes this issue.

  end subroutine MOM_generic_tracer_surface_state

!ALL PE subroutine on Ocean!  Due to otpm design the fluxes should be initialized like this on ALL PE's!
  subroutine MOM_generic_flux_init(verbosity)
    integer, optional, intent(in) :: verbosity  !< A 0-9 integer indicating a level of verbosity.

    character(len=128), parameter :: sub_name = 'MOM_generic_flux_init'
    type(g_tracer_type), pointer :: g_tracer_list,g_tracer,g_tracer_next

    if (.not. g_registered) then
      call generic_tracer_register()
      g_registered = .true.
    endif

    call generic_tracer_get_list(g_tracer_list)
    if (.NOT. associated(g_tracer_list)) then
      call MOM_error(WARNING, trim(sub_name)// ": No generic tracer in the list.")
      return
    endif

    g_tracer=>g_tracer_list
    do

      call g_tracer_flux_init(g_tracer, verbosity=verbosity)

      ! traverse the linked list till hit NULL
      call g_tracer_get_next(g_tracer, g_tracer_next)
      if (.NOT. associated(g_tracer_next)) exit
      g_tracer=>g_tracer_next

    enddo

  end subroutine MOM_generic_flux_init

  subroutine MOM_generic_tracer_fluxes_accumulate(flux_tmp, weight)
    type(forcing), intent(in)    :: flux_tmp  !< A structure containing pointers to
                                              !! thermodynamic and tracer forcing fields.
    real,          intent(in)    :: weight    !< A weight for accumulating this flux [nondim]

    call generic_tracer_coupler_accumulate(flux_tmp%tr_fluxes, weight)

  end subroutine MOM_generic_tracer_fluxes_accumulate

  !> Copy the requested tracer into an array.
  subroutine MOM_generic_tracer_get(name,member,array, CS)
    character(len=*),         intent(in)  :: name   !< Name of requested tracer.
    character(len=*),         intent(in)  :: member !< The tracer element to return.
    real, dimension(:,:,:),   intent(out) :: array  !< Array filled by this routine, in arbitrary units [A]
    type(MOM_generic_tracer_CS), pointer  :: CS     !< Pointer to the control structure for this module.

    ! Local variables
    real, dimension(:,:,:),   pointer :: array_ptr  ! The tracer in the generic tracer structures, in
                                                    ! arbitrary units [A]
    character(len=128), parameter :: sub_name = 'MOM_generic_tracer_get'

    call g_tracer_get_pointer(CS%g_tracer_list,name,member,array_ptr)
    array(:,:,:) = array_ptr(:,:,:)

  end subroutine MOM_generic_tracer_get

  !> This subroutine deallocates the memory owned by this module.
  subroutine end_MOM_generic_tracer(CS)
    type(MOM_generic_tracer_CS), pointer :: CS   !< Pointer to the control structure for this module.

    call generic_tracer_end()

    if (associated(CS)) then
      deallocate(CS)
    endif
  end subroutine end_MOM_generic_tracer

!----------------------------------------------------------------
! <CONTACT EMAIL="Niki.Zadeh@noaa.gov"> Niki Zadeh
! </CONTACT>
!
! <REVIEWER EMAIL="William.Cooke@noaa.gov"> William Cooke
! </REVIEWER>
!
! <OVERVIEW>
!  This module drives the generic version of tracers TOPAZ and CFC
! </OVERVIEW>
!----------------------------------------------------------------

end module MOM_generic_tracer