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 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_io, only : file_exists, mom_read_data, slasher
  use mom_restart, only : register_restart_field, query_initialized, mom_restart_cs
  use mom_spatial_means, only : global_area_mean
  use mom_sponge, only : set_up_sponge_field, sponge_cs
  use mom_ale_sponge, only : set_up_ale_sponge_field, ale_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_open_boundary, only : ocean_obc_type
  use mom_verticalgrid, only : verticalgrid_type


  implicit none ; private

  !> An 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

  !> 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.
     real :: tracer_land_val = -1.0 !< The values of tracers used where  land is masked out.
     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

     !> Pointer to the first element of the linked list of generic tracers.
     type(g_tracer_type), pointer :: g_tracer_list => null()

     integer :: h_to_m !< Auxiliary to access GV%H_to_m in routines that do not have access to GV

  end type mom_generic_tracer_cs

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

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),       pointer      :: restart_cs !< Pointer to the restart control structure.

! Local variables
    logical :: register_mom_generic_tracer

    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, k,i,j,axes(3)
    type(g_tracer_type), pointer      :: g_tracer,g_tracer_next
    character(len=fm_string_len)      :: g_tracer_name,longname,units
    real, dimension(:,:,:,:), pointer   :: tr_field
    real, dimension(:,:,:), pointer     :: tr_ptr
    real, dimension(HI%isd:HI%ied, HI%jsd:HI%jed,GV%ke)         :: grid_tmask
    integer, dimension(HI%isd:HI%ied, HI%jsd:HI%jed)           :: grid_kmt

    register_mom_generic_tracer = .false.
    if (associated(cs)) then
      call mom_error(warning, "register_MOM_generic_tracer called with an "// &
            "associated control structure.")
      return
    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 trcer 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.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 prognastic 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

  !>  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, 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_(G)), intent(in) :: h    !< Layer thicknesses [H ~> m or kg m-2]
    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
    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
    real, dimension(:,:,:), pointer     :: tr_ptr
    real,    dimension(G%isd:G%ied, G%jsd:G%jed,1:G%ke) :: grid_tmask
    integer, dimension(G%isd:G%ied, G%jsd:G%jed)        :: grid_kmt

    !! 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 = g%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

    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(h, 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               )

         !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, us)
            if (.not.ok) then
              ok = tracer_z_init(tr_ptr, h, cs%IC_file, trim(g_tracer_name), g, 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
      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 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) then
        grid_tmask(i,j,:) = 1.0
        grid_kmt(i,j) = g%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

    cs%H_to_m = gv%H_to_m

  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, 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_(G)), &
                             intent(in) :: h_old !< Layer thickness before entrainment [H ~> m or kg m-2].
    real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                             intent(in) :: h_new !< Layer thickness after entrainment [H ~> m or kg m-2].
    real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                             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_(G)), &
                             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 [s]
    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 !< Limits how much water can be fluxed out of
                                                 !! the top layer 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,trunoff_array,runoff_tracer_flux_array

    real :: surface_field(szi_(g),szj_(g))
    real :: dz_ml(szi_(g),szj_(g))  ! The mixed layer depth in the MKS units used for generic tracers [m]
    real :: sosga

    real, dimension(G%isd:G%ied,G%jsd:G%jed,G%ke) :: rho_dzt, dzt
    real, dimension(SZI_(G),SZJ_(G),SZK_(G))      :: h_work
    integer :: i, j, k, isc, iec, jsc, jec, nk

    isc = g%isc ; iec = g%iec ; jsc = g%jsc ; jec = g%jec ; nk = g%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)) 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(:,:) * &
                   g%US%RZ_T_to_kg_m2s*fluxes%lrunoff(:,:)
          stf_array = stf_array + runoff_tracer_flux_array
       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
    do k = 1, nk ; do j = jsc, jec ; do i = isc, iec  !{
      dzt(i,j,k) = gv%H_to_m * h_old(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) = g%US%Z_to_m * hml(i,j)
    enddo ; enddo
    sosga = global_area_mean(surface_field, g)

    !
    !Calculate tendencies (i.e., field changes at dt) from the sources / sinks
    !
    if ((g%US%L_to_m == 1.0) .and. (g%US%RZ_to_kg_m2 == 1.0) .and. (g%US%s_to_T == 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
      call generic_tracer_source(tv%T, 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, optics%sw_pen_band, optics%opacity_band, &
               internal_heat=g%US%RZ_to_kg_m2*tv%internal_heat(:,:), &
               frunoff=g%US%RZ_T_to_kg_m2s*fluxes%frunoff(:,:), sosga=sosga)
    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), g%US%s_to_T*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, 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, 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(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_(G)), intent(in) :: h !< Layer thicknesses [H ~> m or kg m-2]
    real, 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
    real, dimension(:,:,:), pointer     :: tr_ptr
    character(len=128), parameter :: sub_name = 'MOM_generic_tracer_stock'

    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

    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)

      stocks(m) = 0.0
      tr_ptr => tr_field(:,:,:,1)
      do k=1,nz ; do j=js,je ; do i=is,ie
        stocks(m) = stocks(m) + tr_ptr(i,j,k) * &
                               (g%mask2dT(i,j) * g%US%L_to_m**2*g%areaT(i,j) * h(i,j,k))
      enddo ; enddo ; enddo
      stocks(m) = gv%H_to_kg_m2 * stocks(m)

      !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 find 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 gone through.
  function mom_generic_tracer_min_max(ind_start, got_minmax, gmin, gmax, xgmin, ygmin, zgmin, &
                                      xgmax, ygmax, zgmax , G, CS, names, units)
    use mpp_utilities_mod, only: mpp_array_global_min_max
    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, in kg
                                                            !! times concentration units.
    real, dimension(:),             intent(out)   :: gmax   !< Global maximum of each tracer, in kg
                                                            !! times concentration units.
    real, dimension(:),             intent(out)   :: xgmin  !< The x-position of the global minimum
    real, dimension(:),             intent(out)   :: ygmin  !< The y-position of the global minimum
    real, dimension(:),             intent(out)   :: zgmin  !< The z-position of the global minimum
    real, dimension(:),             intent(out)   :: xgmax  !< The x-position of the global maximum
    real, dimension(:),             intent(out)   :: ygmax  !< The y-position of the global maximum
    real, dimension(:),             intent(out)   :: zgmax  !< The z-position of the global maximum
    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.
    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
    real, dimension(:,:,:), pointer     :: tr_ptr
    character(len=128), parameter :: sub_name = 'MOM_generic_tracer_min_max'

    real, dimension(:,:,:),pointer :: grid_tmask
    integer :: isc,iec,jsc,jec,isd,ied,jsd,jed,nk,ntau

    integer :: i, j, k, is, ie, js, je, nz, m
    real, allocatable, dimension(:) :: geo_z

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

    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,grid_tmask=grid_tmask)

    !  Because the use of a simple z-coordinate can not be assumed, simply
    ! use the layer index as the vertical label.
    allocate(geo_z(nk))
    do k=1,nk ; geo_z(k) = real(k) ; enddo

    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))
      units(m) = trim(units(m))//" kg"
      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)

      call mpp_array_global_min_max(tr_ptr, grid_tmask,isd,jsd,isc,iec,jsc,jec,nk , gmin(m), gmax(m), &
                                    g%geoLonT,g%geoLatT,geo_z,xgmin(m), ygmin(m), zgmin(m), &
                                    xgmax(m), ygmax(m), zgmax(m))

      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, CS)
    type(ocean_grid_type),                 intent(in)    :: g    !< The ocean's 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_(G)), 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

    character(len=128), parameter :: sub_name = 'MOM_generic_tracer_surface_state'
    real, dimension(G%isd:G%ied,G%jsd:G%jed,1:G%ke,1) :: rho0
    real, dimension(G%isd:G%ied,G%jsd:G%jed,1:G%ke) ::  dzt
    type(g_tracer_type), pointer :: g_tracer

    !Set coupler values
    !nnz: fake rho0
    rho0=1.0

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

    sosga = global_area_mean(sfc_state%SSS, g)

    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))

    !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 subrotuine 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.

    integer :: ind
    character(len=fm_string_len)   :: g_tracer_name,longname, package,units,old_package,file_in,file_out
    real :: const_init_value
    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) !### Add this after ocean shared is updated.

      ! 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

    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.
    type(mom_generic_tracer_cs), pointer :: cs   !< Pointer to the control structure for this module.

    real, dimension(:,:,:),   pointer :: array_ptr
    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