MOM_tracer_advect.F90

!>  This module contains the subroutines that advect tracers along coordinate surfaces.
module mom_tracer_advect

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

use mom_cpu_clock,       only : cpu_clock_id, cpu_clock_begin, cpu_clock_end
use mom_cpu_clock,       only : clock_module, clock_routine
use mom_diag_mediator,   only : post_data, query_averaging_enabled, diag_ctrl
use mom_diag_mediator,   only : register_diag_field, safe_alloc_ptr, time_type
use mom_domains,         only : sum_across_pes, max_across_pes
use mom_domains,         only : create_group_pass, do_group_pass, group_pass_type, pass_var
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_grid,            only : ocean_grid_type
use mom_open_boundary,   only : ocean_obc_type, obc_none, obc_direction_e
use mom_open_boundary,   only : obc_direction_w, obc_direction_n, obc_direction_s
use mom_open_boundary,   only : obc_segment_type
use mom_tracer_registry, only : tracer_registry_type, tracer_type
use mom_unit_scaling,    only : unit_scale_type
use mom_verticalgrid,    only : verticalgrid_type
implicit none ; private

#include <MOM_memory.h>

public advect_tracer
public tracer_advect_init
public tracer_advect_end

!> Control structure for this module
type, public :: tracer_advect_cs ; private
  real    :: dt                    !< The baroclinic dynamics time step [T ~> s].
  type(diag_ctrl), pointer :: diag !< A structure that is used to regulate the
                                   !< timing of diagnostic output.
  logical :: debug                 !< If true, write verbose checksums for debugging purposes.
  logical :: useppm                !< If true, use PPM instead of PLM
  logical :: usehuynh              !< If true, use the Huynh scheme for PPM interface values
  type(group_pass_type) :: pass_uhr_vhr_t_hprev !< A structred used for group passes
end type tracer_advect_cs

!>@{ CPU time clocks
integer :: id_clock_advect
integer :: id_clock_pass
integer :: id_clock_sync
!>@}

contains

!> This routine time steps the tracer concentration using a
!! monotonic, conservative, weakly diffusive scheme.
subroutine advect_tracer(h_end, uhtr, vhtr, OBC, dt, G, GV, US, CS, Reg, &
      h_prev_opt, max_iter_in, x_first_in, uhr_out, vhr_out, h_out)
  type(ocean_grid_type),   intent(inout) :: g     !< ocean grid structure
  type(verticalgrid_type), intent(in)    :: gv    !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: h_end !< layer thickness after advection [H ~> m or kg m-2]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: uhtr  !< accumulated volume/mass flux through zonal face [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: vhtr  !< accumulated volume/mass flux through merid face [H L2 ~> m3 or kg]
  type(ocean_obc_type),    pointer       :: obc   !< specifies whether, where, and what OBCs are used
  real,                    intent(in)    :: dt    !< time increment [T ~> s]
  type(unit_scale_type),   intent(in)    :: us    !< A dimensional unit scaling type
  type(tracer_advect_cs),  pointer       :: cs    !< control structure for module
  type(tracer_registry_type), pointer    :: reg   !< pointer to tracer registry
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                 optional, intent(in)    :: h_prev_opt !< layer thickness before advection [H ~> m or kg m-2]
  integer,       optional, intent(in)    :: max_iter_in !< The maximum number of iterations
  logical,       optional, intent(in)    :: x_first_in !< If present, indicate whether to update
                                                  !! first in the x- or y-direction.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(out)    :: uhr_out  !< accumulated volume/mass flux through zonal face
                                                  !! [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                 optional, intent(out)    :: vhr_out  !< accumulated volume/mass flux through merid face
                                                  !! [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                 optional, intent(out)    :: h_out !< layer thickness before advection [H ~> m or kg m-2]

  type(tracer_type) :: tr(max_fields_) ! The array of registered tracers
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: &
    hprev           ! cell volume at the end of previous tracer change [H L2 ~> m3 or kg]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)) :: &
    uhr             ! The remaining zonal thickness flux [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)) :: &
    vhr             ! The remaining meridional thickness fluxes [H L2 ~> m3 or kg]
  real :: uh_neglect(szib_(g),szj_(g)) ! uh_neglect and vh_neglect are the
  real :: vh_neglect(szi_(g),szjb_(g)) ! magnitude of remaining transports that
                                       ! can be simply discarded [H L2 ~> m3 or kg].

  real :: landvolfill                   ! An arbitrary? nonzero cell volume [H L2 ~> m3 or kg].
  real :: idt                           ! 1/dt [T-1 ~> s-1].
  logical :: domore_u(szj_(g),szk_(g))  ! domore__ indicate whether there is more
  logical :: domore_v(szjb_(g),szk_(g)) ! advection to be done in the corresponding
                                        ! row or column.
  logical :: x_first            ! If true, advect in the x-direction first.
  integer :: max_iter           ! maximum number of iterations in each layer
  integer :: domore_k(szk_(g))
  integer :: stencil            ! stencil of the advection scheme
  integer :: nsten_halo         ! number of stencils that fit in the halos
  integer :: i, j, k, m, is, ie, js, je, isd, ied, jsd, jed, nz, itt, ntr, do_any
  integer :: isv, iev, jsv, jev ! The valid range of the indices.
  integer :: isdb, iedb, jsdb, jedb

  domore_u(:,:) = .false.
  domore_v(:,:) = .false.
  is  = g%isc ; ie  = g%iec ; js  = g%jsc ; je  = g%jec ; nz = gv%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
  landvolfill = 1.0e-20         ! This is arbitrary, but must be positive.
  stencil = 2                   ! The scheme's stencil; 2 for PLM and PPM:H3

  if (.not. associated(cs)) call mom_error(fatal, "MOM_tracer_advect: "// &
       "tracer_advect_init must be called before advect_tracer.")
  if (.not. associated(reg)) call mom_error(fatal, "MOM_tracer_advect: "// &
       "register_tracer must be called before advect_tracer.")
  if (reg%ntr==0) return
  call cpu_clock_begin(id_clock_advect)
  x_first = (mod(g%first_direction,2) == 0)

  ! increase stencil size for Colella & Woodward PPM
  if (cs%usePPM .and. .not. cs%useHuynh) stencil = 3

  ntr = reg%ntr
  do m=1,ntr ; tr(m) = reg%Tr(m) ; enddo
  idt = 1.0 / dt

  max_iter = 2*int(ceiling(dt/cs%dt)) + 1

  if (present(max_iter_in)) max_iter = max_iter_in
  if (present(x_first_in))  x_first = x_first_in
  call cpu_clock_begin(id_clock_pass)
  call create_group_pass(cs%pass_uhr_vhr_t_hprev, uhr, vhr, g%Domain)
  call create_group_pass(cs%pass_uhr_vhr_t_hprev, hprev, g%Domain)
  do m=1,ntr
    call create_group_pass(cs%pass_uhr_vhr_t_hprev, tr(m)%t, g%Domain)
  enddo
  call cpu_clock_end(id_clock_pass)

!$OMP parallel default(none) shared(nz,jsd,jed,IsdB,IedB,uhr,jsdB,jedB,Isd,Ied,vhr, &
!$OMP                               hprev,domore_k,js,je,is,ie,uhtr,vhtr,G,GV,h_end,&
!$OMP                               uh_neglect,vh_neglect,ntr,Tr,h_prev_opt)

  ! This initializes the halos of uhr and vhr because pass_vector might do
  ! calculations on them, even though they are never used.
  !$OMP do
  do k=1,nz
    do j=jsd,jed ; do i=isdb,iedb ; uhr(i,j,k) = 0.0 ; enddo ; enddo
    do j=jsdb,jedb ; do i=isd,ied ; vhr(i,j,k) = 0.0 ; enddo ; enddo
    do j=jsd,jed ; do i=isd,ied ; hprev(i,j,k) = 0.0 ; enddo ; enddo
    domore_k(k)=1
    !  Put the remaining (total) thickness fluxes into uhr and vhr.
    do j=js,je ; do i=is-1,ie ; uhr(i,j,k) = uhtr(i,j,k) ; enddo ; enddo
    do j=js-1,je ; do i=is,ie ; vhr(i,j,k) = vhtr(i,j,k) ; enddo ; enddo
    if (.not. present(h_prev_opt)) then
    !   This loop reconstructs the thickness field the last time that the
    ! tracers were updated, probably just after the diabatic forcing.  A useful
    ! diagnostic could be to compare this reconstruction with that older value.
      do i=is,ie ; do j=js,je
        hprev(i,j,k) = max(0.0, g%areaT(i,j)*h_end(i,j,k) + &
             ((uhr(i,j,k) - uhr(i-1,j,k)) + (vhr(i,j,k) - vhr(i,j-1,k))))
    ! In the case that the layer is now dramatically thinner than it was previously,
    ! add a bit of mass to avoid truncation errors.  This will lead to
    ! non-conservation of tracers
        hprev(i,j,k) = hprev(i,j,k) + &
                       max(0.0, 1.0e-13*hprev(i,j,k) - g%areaT(i,j)*h_end(i,j,k))
      enddo ; enddo
    else
      do i=is,ie ; do j=js,je
        hprev(i,j,k) = h_prev_opt(i,j,k)
      enddo ; enddo
    endif
  enddo


  !$OMP do
  do j=jsd,jed ; do i=isd,ied-1
    uh_neglect(i,j) = gv%H_subroundoff * min(g%areaT(i,j), g%areaT(i+1,j))
  enddo ; enddo
  !$OMP do
  do j=jsd,jed-1 ; do i=isd,ied
    vh_neglect(i,j) = gv%H_subroundoff * min(g%areaT(i,j), g%areaT(i,j+1))
  enddo ; enddo

  ! initialize diagnostic fluxes and tendencies
  !$OMP do
  do m=1,ntr
    if (associated(tr(m)%ad_x)) then
      do k=1,nz ; do j=jsd,jed ; do i=isd,ied
        tr(m)%ad_x(i,j,k) = 0.0
      enddo ; enddo ; enddo
    endif
    if (associated(tr(m)%ad_y)) then
      do k=1,nz ; do j=jsd,jed ; do i=isd,ied
        tr(m)%ad_y(i,j,k) = 0.0
      enddo ; enddo ; enddo
    endif
    if (associated(tr(m)%advection_xy)) then
      do k=1,nz ; do j=jsd,jed ; do i=isd,ied
        tr(m)%advection_xy(i,j,k) = 0.0
      enddo ; enddo ; enddo
    endif
    if (associated(tr(m)%ad2d_x)) then
      do j=jsd,jed ; do i=isd,ied ; tr(m)%ad2d_x(i,j) = 0.0 ; enddo ; enddo
    endif
    if (associated(tr(m)%ad2d_y)) then
      do j=jsd,jed ; do i=isd,ied ; tr(m)%ad2d_y(i,j) = 0.0 ; enddo ; enddo
    endif
  enddo
  !$OMP end parallel

  isv = is ; iev = ie ; jsv = js ; jev = je

  do itt=1,max_iter

    if (isv > is-stencil) then
      call do_group_pass(cs%pass_uhr_vhr_t_hprev, g%Domain, clock=id_clock_pass)

      nsten_halo = min(is-isd,ied-ie,js-jsd,jed-je)/stencil
      isv = is-nsten_halo*stencil ; jsv = js-nsten_halo*stencil
      iev = ie+nsten_halo*stencil ; jev = je+nsten_halo*stencil
      ! Reevaluate domore_u & domore_v unless the valid range is the same size as
      ! before.  Also, do this if there is Strang splitting.
      if ((nsten_halo > 1) .or. (itt==1)) then
        !$OMP parallel do default(shared)
        do k=1,nz ; if (domore_k(k) > 0) then
          do j=jsv,jev ; if (.not.domore_u(j,k)) then
            do i=isv+stencil-1,iev-stencil; if (uhr(i,j,k) /= 0.0) then
              domore_u(j,k) = .true. ; exit
            endif ; enddo ! i-loop
          endif ; enddo
          do j=jsv+stencil-1,jev-stencil ; if (.not.domore_v(j,k)) then
            do i=isv+stencil,iev-stencil; if (vhr(i,j,k) /= 0.0) then
              domore_v(j,k) = .true. ; exit
            endif ; enddo ! i-loop
          endif ; enddo

          !   At this point, domore_k is global.  Change it so that it indicates
          ! whether any work is needed on a layer on this processor.
          domore_k(k) = 0
          do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo
          do j=jsv+stencil-1,jev-stencil ; if (domore_v(j,k)) domore_k(k) = 1 ; enddo

        endif ; enddo ! k-loop
      endif
    endif

    ! Set the range of valid points after this iteration.
    isv = isv + stencil ; iev = iev - stencil
    jsv = jsv + stencil ; jev = jev - stencil

    !  To ensure positive definiteness of the thickness at each iteration, the
    !  mass fluxes out of each layer are checked each step, and limited to keep
    !  the thicknesses positive.  This means that several iterations may be required
    !  for all the transport to happen.  The sum over domore_k keeps the processors
    !  synchronized.  This may not be very efficient, but it should be reliable.

    !$OMP parallel default(shared)

    if (x_first) then

      !$OMP do ordered
      do k=1,nz ; if (domore_k(k) > 0) then
        ! First, advect zonally.
        call advect_x(tr, hprev, uhr, uh_neglect, obc, domore_u, ntr, idt, &
                      isv, iev, jsv-stencil, jev+stencil, k, g, gv, us, cs%usePPM, cs%useHuynh)
      endif ; enddo

      !$OMP do ordered
      do k=1,nz ; if (domore_k(k) > 0) then
        !  Next, advect meridionally.
        call advect_y(tr, hprev, vhr, vh_neglect, obc, domore_v, ntr, idt, &
                      isv, iev, jsv, jev, k, g, gv, us, cs%usePPM, cs%useHuynh)

        ! Update domore_k(k) for the next iteration
        domore_k(k) = 0
        do j=jsv-stencil,jev+stencil ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo
        do j=jsv-1,jev ; if (domore_v(j,k)) domore_k(k) = 1 ; enddo

      endif ; enddo

    else

      !$OMP do ordered
      do k=1,nz ; if (domore_k(k) > 0) then
        ! First, advect meridionally.
        call advect_y(tr, hprev, vhr, vh_neglect, obc, domore_v, ntr, idt, &
                      isv-stencil, iev+stencil, jsv, jev, k, g, gv, us, cs%usePPM, cs%useHuynh)
      endif ; enddo

      !$OMP do ordered
      do k=1,nz ; if (domore_k(k) > 0) then
        ! Next, advect zonally.
        call advect_x(tr, hprev, uhr, uh_neglect, obc, domore_u, ntr, idt, &
                      isv, iev, jsv, jev, k, g, gv, us, cs%usePPM, cs%useHuynh)

        ! Update domore_k(k) for the next iteration
        domore_k(k) = 0
        do j=jsv,jev ; if (domore_u(j,k)) domore_k(k) = 1 ; enddo
        do j=jsv-1,jev ; if (domore_v(j,k)) domore_k(k) = 1 ; enddo
      endif ; enddo

    endif ! x_first

    !$OMP end parallel

    ! If the advection just isn't finishing after max_iter, move on.
    if (itt >= max_iter) then
      exit
    endif

    ! Exit if there are no layers that need more iterations.
    if (isv > is-stencil) then
      do_any = 0
      call cpu_clock_begin(id_clock_sync)
      call sum_across_pes(domore_k(:), nz)
      call cpu_clock_end(id_clock_sync)
      do k=1,nz ; do_any = do_any + domore_k(k) ; enddo
      if (do_any == 0) then
        exit
      endif

    endif

  enddo ! Iterations loop

  if (present(uhr_out)) uhr_out(:,:,:) = uhr(:,:,:)
  if (present(vhr_out)) vhr_out(:,:,:) = vhr(:,:,:)
  if (present(h_out)) h_out(:,:,:) = hprev(:,:,:)

  call cpu_clock_end(id_clock_advect)

end subroutine advect_tracer


!> This subroutine does 1-d flux-form advection in the zonal direction using
!! a monotonic piecewise linear scheme.
subroutine advect_x(Tr, hprev, uhr, uh_neglect, OBC, domore_u, ntr, Idt, &
                    is, ie, js, je, k, G, GV, US, usePPM, useHuynh)
  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(tracer_type), dimension(ntr),         intent(inout) :: Tr   !< The array of registered tracers to work on
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: hprev !< cell volume at the end of previous
                                                                  !! tracer change [H L2 ~> m3 or kg]
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), intent(inout) :: uhr !< accumulated volume/mass flux through
                                                                  !! the zonal face [H L2 ~> m3 or kg]
   real, dimension(SZIB_(G),SZJ_(G)),        intent(in)    :: uh_neglect !< A tiny zonal mass flux that can
                                                                  !! be neglected [H L2 ~> m3 or kg]
  type(ocean_OBC_type),                      pointer       :: OBC !< specifies whether, where, and what OBCs are used
  logical, dimension(SZJ_(G),SZK_(G)),       intent(inout) :: domore_u !< If true, there is more advection to be
                                                                  !! done in this u-row
  real,                                      intent(in)    :: Idt !< The inverse of dt [T-1 ~> s-1]
  integer,                                   intent(in)    :: ntr !< The number of tracers
  integer,                                   intent(in)    :: is  !< The starting tracer i-index to work on
  integer,                                   intent(in)    :: ie  !< The ending tracer i-index to work on
  integer,                                   intent(in)    :: js  !< The starting tracer j-index to work on
  integer,                                   intent(in)    :: je  !< The ending tracer j-index to work on
  integer,                                   intent(in)    :: k   !< The k-level to work on
  type(unit_scale_type),                     intent(in)    :: US  !< A dimensional unit scaling type
  logical,                                   intent(in)    :: usePPM !< If true, use PPM instead of PLM
  logical,                                   intent(in)    :: useHuynh !< If true, use the Huynh scheme
                                                                     !! for PPM interface values

  real, dimension(SZI_(G),ntr) :: &
    slope_x             ! The concentration slope per grid point [conc].
  real, dimension(SZIB_(G),SZJ_(G),ntr) :: &
    flux_x              ! The tracer flux across a boundary [H L2 conc ~> m3 conc or kg conc].
  real, dimension(SZI_(G),ntr) :: &
    T_tmp               ! The copy of the tracer concentration at constant i,k [H m2 conc ~> m3 conc or kg conc].

  real :: maxslope      ! The maximum concentration slope per grid point
                        ! consistent with monotonicity [conc].
  real :: hup, hlos     ! hup is the upwind volume, hlos is the
                        ! part of that volume that might be lost
                        ! due to advection out the other side of
                        ! the grid box, both in [H L2 ~> m3 or kg].
  real :: uhh(SZIB_(G)) ! The zonal flux that occurs during the
                        ! current iteration [H L2 ~> m3 or kg].
  real, dimension(SZIB_(G)) :: &
    hlst, &             ! Work variable [H L2 ~> m3 or kg].
    Ihnew, &            ! Work variable [H-1 L-2 ~> m-3 or kg-1].
    CFL                 ! The absolute value of the advective upwind-cell CFL number [nondim].
  real :: min_h         ! The minimum thickness that can be realized during
                        ! any of the passes [H ~> m or kg m-2].
  real :: tiny_h        ! The smallest numerically invertable thickness [H ~> m or kg m-2].
  real :: h_neglect     ! A thickness that is so small it is usually lost
                        ! in roundoff and can be neglected [H ~> m or kg m-2].
  logical :: do_i(SZIB_(G),SZJ_(G))     ! If true, work on given points.
  logical :: do_any_i
  integer :: i, j, m, n, i_up, stencil
  real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
  real :: fac1,u_L_in,u_L_out  ! terms used for time-stepping OBC reservoirs
  type(OBC_segment_type), pointer :: segment=>null()
  logical :: usePLMslope
  logical, dimension(SZJ_(G),SZK_(G)) :: domore_u_initial

  ! keep a local copy of the initial values of domore_u, which is to be used when computing ad2d_x
  ! diagnostic at the end of this subroutine.
  domore_u_initial = domore_u

  useplmslope = .not. (useppm .and. usehuynh)
  ! stencil for calculating slope values
  stencil = 1
  if (useppm .and. .not. usehuynh) stencil = 2

  min_h = 0.1*gv%Angstrom_H
  tiny_h = tiny(min_h)
  h_neglect = gv%H_subroundoff

  do i=is-1,ie ; cfl(i) = 0.0 ; enddo

  do j=js,je ; if (domore_u(j,k)) then
    domore_u(j,k) = .false.

    ! Calculate the i-direction profiles (slopes) of each tracer that is being advected.
    if (useplmslope) then
      do m=1,ntr ; do i=is-stencil,ie+stencil
       !if (ABS(Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k)) < &
       !    ABS(Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k))) then
       !  maxslope = 4.0*(Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k))
       !else
       !  maxslope = 4.0*(Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k))
       !endif
       !if ((Tr(m)%t(i+1,j,k)-Tr(m)%t(i,j,k)) * (Tr(m)%t(i,j,k)-Tr(m)%t(i-1,j,k)) < 0.0) then
       !  slope_x(i,m) = 0.0
       !elseif (ABS(Tr(m)%t(i+1,j,k)-Tr(m)%t(i-1,j,k))<ABS(maxslope)) then
       !  slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * &
       !                 0.5*(Tr(m)%t(i+1,j,k)-Tr(m)%t(i-1,j,k))
       !else
       !  slope_x(i,m) = G%mask2dCu(I,j)*G%mask2dCu(I-1,j) * 0.5*maxslope
       !endif
        tp = tr(m)%t(i+1,j,k) ; tc = tr(m)%t(i,j,k) ; tm = tr(m)%t(i-1,j,k)
        dmx = max( tp, tc, tm ) - tc
        dmn= tc - min( tp, tc, tm )
        slope_x(i,m) = g%mask2dCu(i,j)*g%mask2dCu(i-1,j) * &
            sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
      enddo ; enddo
    endif ! usePLMslope

    ! make a copy of the tracers in case values need to be overridden for OBCs
    do m = 1,ntr
      do i=g%isd,g%ied
        t_tmp(i,m) = tr(m)%t(i,j,k)
      enddo
    enddo
    ! loop through open boundaries and recalculate flux terms
    if (associated(obc)) then ; if (obc%OBC_pe) then
       do n=1,obc%number_of_segments
         segment=>obc%segment(n)
         if (.not. associated(segment%tr_Reg)) cycle
         if (segment%is_E_or_W) then
           if (j>=segment%HI%jsd .and. j<=segment%HI%jed) then
              i = segment%HI%IsdB
              do m = 1,ntr ! replace tracers with OBC values
                if (associated(segment%tr_Reg%Tr(m)%tres)) then
                   if (segment%direction == obc_direction_w) then
                      t_tmp(i,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
                   else
                      t_tmp(i+1,m) = segment%tr_Reg%Tr(m)%tres(i,j,k)
                   endif
                else
                   if (segment%direction == obc_direction_w) then
                      t_tmp(i,m) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
                   else
                      t_tmp(i+1,m) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
                   endif
                endif
              enddo
              do m = 1,ntr ! Apply update tracer values for slope calculation
                do i=segment%HI%IsdB-1,segment%HI%IsdB+1
                  tp = t_tmp(i+1,m) ; tc = t_tmp(i,m) ; tm = t_tmp(i-1,m)
                  dmx = max( tp, tc, tm ) - tc
                  dmn= tc - min( tp, tc, tm )
                  slope_x(i,m) = g%mask2dCu(i,j)*g%mask2dCu(i-1,j) * &
                       sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
                enddo
              enddo

           endif
         endif
       enddo
    endif; endif


    ! Calculate the i-direction fluxes of each tracer, using as much
    ! the minimum of the remaining mass flux (uhr) and the half the mass
    ! in the cell plus whatever part of its half of the mass flux that
    ! the flux through the other side does not require.
    do i=is-1,ie
      if ((uhr(i,j,k) == 0.0) .or. &
          ((uhr(i,j,k) < 0.0) .and. (hprev(i+1,j,k) <= tiny_h)) .or. &
          ((uhr(i,j,k) > 0.0) .and. (hprev(i,j,k) <= tiny_h)) ) then
        uhh(i) = 0.0
        cfl(i) = 0.0
      elseif (uhr(i,j,k) < 0.0) then
        hup = hprev(i+1,j,k) - g%areaT(i+1,j)*min_h
        hlos = max(0.0, uhr(i+1,j,k))
        if ((((hup - hlos) + uhr(i,j,k)) < 0.0) .and. &
            ((0.5*hup + uhr(i,j,k)) < 0.0)) then
          uhh(i) = min(-0.5*hup, -hup+hlos, 0.0)
          domore_u(j,k) = .true.
        else
          uhh(i) = uhr(i,j,k)
        endif
        cfl(i) = - uhh(i) / (hprev(i+1,j,k))  ! CFL is positive
      else
        hup = hprev(i,j,k) - g%areaT(i,j)*min_h
        hlos = max(0.0, -uhr(i-1,j,k))
        if ((((hup - hlos) - uhr(i,j,k)) < 0.0) .and. &
            ((0.5*hup - uhr(i,j,k)) < 0.0)) then
          uhh(i) = max(0.5*hup, hup-hlos, 0.0)
          domore_u(j,k) = .true.
        else
          uhh(i) = uhr(i,j,k)
        endif
        cfl(i) = uhh(i) / (hprev(i,j,k))  ! CFL is positive
      endif
    enddo


    if (useppm) then
      do m=1,ntr ; do i=is-1,ie
        ! centre cell depending on upstream direction
        if (uhh(i) >= 0.0) then
          i_up = i
        else
          i_up = i+1
        endif

        ! Implementation of PPM-H3
        tp = t_tmp(i_up+1,m) ; tc = t_tmp(i_up,m) ; tm = t_tmp(i_up-1,m)

        if (usehuynh) then
          al = ( 5.*tc + ( 2.*tm - tp ) )/6. ! H3 estimate
          al = max( min(tc,tm), al) ; al = min( max(tc,tm), al) ! Bound
          ar = ( 5.*tc + ( 2.*tp - tm ) )/6. ! H3 estimate
          ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar) ! Bound
        else
          al = 0.5 * ((tm + tc) + (slope_x(i_up-1,m) - slope_x(i_up,m)) / 3.)
          ar = 0.5 * ((tc + tp) + (slope_x(i_up,m) - slope_x(i_up+1,m)) / 3.)
        endif

        da = ar - al ; ma = 0.5*( ar + al )
        if (g%mask2dCu(i_up,j)*g%mask2dCu(i_up-1,j)*(tp-tc)*(tc-tm) <= 0.) then
          al = tc ; ar = tc ! PCM for local extremum and bounadry cells
        elseif ( da*(tc-ma) > (da*da)/6. ) then
          al = 3.*tc - 2.*ar
        elseif ( da*(tc-ma) < - (da*da)/6. ) then
          ar = 3.*tc - 2.*al
        endif

        a6 = 6.*tc - 3. * (ar + al) ! Curvature

        if (uhh(i) >= 0.0) then
          flux_x(i,j,m) = uhh(i)*( ar - 0.5 * cfl(i) * ( &
               ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
        else
          flux_x(i,j,m) = uhh(i)*( al + 0.5 * cfl(i) * ( &
               ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
        endif
      enddo ; enddo
    else ! PLM
      do m=1,ntr ; do i=is-1,ie
        if (uhh(i) >= 0.0) then
          ! Indirect implementation of PLM
         !aL = Tr(m)%t(i,j,k) - 0.5 * slope_x(i,m)
         !aR = Tr(m)%t(i,j,k) + 0.5 * slope_x(i,m)
         !flux_x(I,j,m) = uhh(I)*( aR - 0.5 * (aR-aL) * CFL(I) )
          ! Alternative implementation of PLM
          tc = t_tmp(i,m)
          flux_x(i,j,m) = uhh(i)*( tc + 0.5 * slope_x(i,m) * ( 1. - cfl(i) ) )
        else
          ! Indirect implementation of PLM
         !aL = Tr(m)%t(i+1,j,k) - 0.5 * slope_x(i+1,m)
         !aR = Tr(m)%t(i+1,j,k) + 0.5 * slope_x(i+1,m)
         !flux_x(I,j,m) = uhh(I)*( aL + 0.5 * (aR-aL) * CFL(I) )
          ! Alternative implementation of PLM
          tc = t_tmp(i+1,m)
          flux_x(i,j,m) = uhh(i)*( tc - 0.5 * slope_x(i+1,m) * ( 1. - cfl(i) ) )
        endif
      enddo ; enddo
    endif ! usePPM

    if (associated(obc)) then ; if (obc%OBC_pe) then
      if (obc%specified_u_BCs_exist_globally .or. obc%open_u_BCs_exist_globally) then
        do n=1,obc%number_of_segments
          segment=>obc%segment(n)
          if (.not. associated(segment%tr_Reg)) cycle
          if (segment%is_E_or_W) then
            if (j>=segment%HI%jsd .and. j<=segment%HI%jed) then
              i = segment%HI%IsdB
              ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
              ! Now changing to simply fixed inflows.
              if ((uhr(i,j,k) > 0.0) .and. (segment%direction == obc_direction_w) .or. &
                 (uhr(i,j,k) < 0.0) .and. (segment%direction == obc_direction_e)) then
                uhh(i) = uhr(i,j,k)
              ! should the reservoir evolve for this case Kate ?? - Nope
                do m=1,ntr
                  if (associated(segment%tr_Reg%Tr(m)%tres)) then
                    flux_x(i,j,m) = uhh(i)*segment%tr_Reg%Tr(m)%tres(i,j,k)
                  else ; flux_x(i,j,m) = uhh(i)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
                enddo
              endif
            endif
          endif
        enddo
      endif

      if (obc%open_u_BCs_exist_globally) then
        do n=1,obc%number_of_segments
          segment=>obc%segment(n)
          i = segment%HI%IsdB
          if (segment%is_E_or_W .and. (j >= segment%HI%jsd .and. j<= segment%HI%jed)) then
            if (segment%specified) cycle
            if (.not. associated(segment%tr_Reg)) cycle

            ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
            if ((uhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
               (uhr(i,j,k) < 0.0) .and. (g%mask2dT(i+1,j) < 0.5)) then
              uhh(i) = uhr(i,j,k)
              do m=1,ntr
                if (associated(segment%tr_Reg%Tr(m)%tres)) then
                  flux_x(i,j,m) = uhh(i)*segment%tr_Reg%Tr(m)%tres(i,j,k)
                else; flux_x(i,j,m) = uhh(i)*segment%tr_Reg%Tr(m)%OBC_inflow_conc; endif
              enddo
            endif
          endif
        enddo
      endif
    endif ; endif

    ! Calculate new tracer concentration in each cell after accounting
    ! for the i-direction fluxes.
    do i=is-1,ie
      uhr(i,j,k) = uhr(i,j,k) - uhh(i)
      if (abs(uhr(i,j,k)) < uh_neglect(i,j)) uhr(i,j,k) = 0.0
    enddo
    do i=is,ie
      if ((uhh(i) /= 0.0) .or. (uhh(i-1) /= 0.0)) then
        do_i(i,j) = .true.
        hlst(i) = hprev(i,j,k)
        hprev(i,j,k) = hprev(i,j,k) - (uhh(i) - uhh(i-1))
        if (hprev(i,j,k) <= 0.0) then ; do_i(i,j) = .false.
        elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j)) then
          hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
          ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
        else ;  ihnew(i) = 1.0 / hprev(i,j,k) ; endif
      else
        do_i(i,j) = .false.
      endif
    enddo

    ! update tracer concentration from i-flux and save some diagnostics
    do m=1,ntr

      ! update tracer
      do i=is,ie
        if (do_i(i,j)) then
          if (ihnew(i) > 0.0) then
            tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
                              (flux_x(i,j,m) - flux_x(i-1,j,m))) * ihnew(i)
          endif
        endif
      enddo

      ! diagnostics
      if (associated(tr(m)%ad_x)) then ; do i=is,ie ; if (do_i(i,j)) then
        tr(m)%ad_x(i,j,k) = tr(m)%ad_x(i,j,k) + flux_x(i,j,m)*idt
      endif ; enddo ; endif

      ! diagnose convergence of flux_x (do not use the Ihnew(i) part of the logic).
      ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt.
      if (associated(tr(m)%advection_xy)) then
        do i=is,ie ; if (do_i(i,j)) then
          tr(m)%advection_xy(i,j,k) = tr(m)%advection_xy(i,j,k) - (flux_x(i,j,m) - flux_x(i-1,j,m)) * &
                                          idt * g%IareaT(i,j)
        endif ; enddo
      endif

    enddo

  endif


  enddo ! End of j-loop.

  ! compute ad2d_x diagnostic outside above j-loop so as to make the summation ordered when OMP is active.

  !$OMP ordered
  do j=js,je ; if (domore_u_initial(j,k)) then
    do m=1,ntr
      if (associated(tr(m)%ad2d_x)) then ; do i=is,ie ; if (do_i(i,j)) then
        tr(m)%ad2d_x(i,j) = tr(m)%ad2d_x(i,j) + flux_x(i,j,m)*idt
      endif ; enddo ; endif
    enddo
  endif ; enddo ! End of j-loop.
  !$OMP end ordered

end subroutine advect_x

!> This subroutine does 1-d flux-form advection using a monotonic piecewise
!! linear scheme.
subroutine advect_y(Tr, hprev, vhr, vh_neglect, OBC, domore_v, ntr, Idt, &
                    is, ie, js, je, k, G, GV, US, usePPM, useHuynh)
  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(tracer_type), dimension(ntr),         intent(inout) :: Tr   !< The array of registered tracers to work on
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  intent(inout) :: hprev !< cell volume at the end of previous
                                                                  !! tracer change [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), intent(inout) :: vhr !< accumulated volume/mass flux through
                                                                  !! the meridional face [H L2 ~> m3 or kg]
  real, dimension(SZI_(G),SZJB_(G)),         intent(inout) :: vh_neglect !< A tiny meridional mass flux that can
                                                                  !! be neglected [H L2 ~> m3 or kg]
  type(ocean_OBC_type),                      pointer       :: OBC !< specifies whether, where, and what OBCs are used
  logical, dimension(SZJB_(G),SZK_(G)),      intent(inout) :: domore_v !< If true, there is more advection to be
                                                                  !! done in this v-row
  real,                                      intent(in)    :: Idt !< The inverse of dt [T-1 ~> s-1]
  integer,                                   intent(in)    :: ntr !< The number of tracers
  integer,                                   intent(in)    :: is  !< The starting tracer i-index to work on
  integer,                                   intent(in)    :: ie  !< The ending tracer i-index to work on
  integer,                                   intent(in)    :: js  !< The starting tracer j-index to work on
  integer,                                   intent(in)    :: je  !< The ending tracer j-index to work on
  integer,                                   intent(in)    :: k   !< The k-level to work on
  type(unit_scale_type),                     intent(in)    :: US  !< A dimensional unit scaling type
  logical,                                   intent(in)    :: usePPM !< If true, use PPM instead of PLM
  logical,                                   intent(in)    :: useHuynh !< If true, use the Huynh scheme
                                                                     !! for PPM interface values

  real, dimension(SZI_(G),ntr,SZJ_(G)) :: &
    slope_y                     ! The concentration slope per grid point [conc].
  real, dimension(SZI_(G),ntr,SZJB_(G)) :: &
       flux_y                      ! The tracer flux across a boundary [H m2 conc ~> m3 conc or kg conc].
  real, dimension(SZI_(G),ntr,SZJB_(G)) :: &
    T_tmp               ! The copy of the tracer concentration at constant i,k [H m2 conc ~> m3 conc or kg conc].
  real :: maxslope              ! The maximum concentration slope per grid point
                                ! consistent with monotonicity [conc].
  real :: vhh(SZI_(G),SZJB_(G)) ! The meridional flux that occurs during the
                                ! current iteration [H L2 ~> m3 or kg].
  real :: hup, hlos             ! hup is the upwind volume, hlos is the
                                ! part of that volume that might be lost
                                ! due to advection out the other side of
                                ! the grid box, both in  [H L2 ~> m3 or kg].
  real, dimension(SZIB_(G)) :: &
    hlst, &             ! Work variable [H L2 ~> m3 or kg].
    Ihnew, &            ! Work variable [H-1 L-2 ~> m-3 or kg-1].
    CFL                 ! The absolute value of the advective upwind-cell CFL number [nondim].
  real :: min_h         ! The minimum thickness that can be realized during
                        ! any of the passes [H ~> m or kg m-2].
  real :: tiny_h        ! The smallest numerically invertable thickness [H ~> m or kg m-2].
  real :: h_neglect     ! A thickness that is so small it is usually lost
                        ! in roundoff and can be neglected [H ~> m or kg m-2].
  logical :: do_j_tr(SZJ_(G))   ! If true, calculate the tracer profiles.
  logical :: do_i(SZIB_(G), SZJ_(G))     ! If true, work on given points.
  logical :: do_any_i
  integer :: i, j, j2, m, n, j_up, stencil
  real :: aR, aL, dMx, dMn, Tp, Tc, Tm, dA, mA, a6
  real :: fac1,v_L_in,v_L_out  ! terms used for time-stepping OBC reservoirs
  type(OBC_segment_type), pointer :: segment=>null()
  logical :: usePLMslope

  useplmslope = .not. (useppm .and. usehuynh)
  ! stencil for calculating slope values
  stencil = 1
  if (useppm .and. .not. usehuynh) stencil = 2

  min_h = 0.1*gv%Angstrom_H
  tiny_h = tiny(min_h)
  h_neglect = gv%H_subroundoff

  ! We conditionally perform work on tracer points: calculating the PLM slope,
  ! and updating tracer concentration within a cell
  ! this depends on whether there is a flux which would affect this tracer point,
  ! as indicated by domore_v. In the case of PPM reconstruction, a flux requires
  ! slope calculations at the two tracer points on either side (as indicated by
  ! the stencil variable), so we account for this with the do_j_tr flag array
  !
  ! Note: this does lead to unnecessary work in updating tracer concentrations,
  ! since that doesn't need a wider stencil with the PPM advection scheme, but
  ! this would require an additional loop, etc.
  do_j_tr(:) = .false.
  do j=js-1,je ; if (domore_v(j,k)) then ; do j2=1-stencil,stencil ; do_j_tr(j+j2) = .true. ; enddo ; endif ; enddo

  ! Calculate the j-direction profiles (slopes) of each tracer that
  ! is being advected.
  if (useplmslope) then
    do j=js-stencil,je+stencil ; if (do_j_tr(j)) then ; do m=1,ntr ; do i=is,ie
      !if (ABS(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k)) < &
      !    ABS(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k))) then
      !  maxslope = 4.0*(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k))
      !else
      !  maxslope = 4.0*(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k))
      !endif
      !if ((Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j,k))*(Tr(m)%t(i,j,k)-Tr(m)%t(i,j-1,k)) < 0.0) then
      !  slope_y(i,m,j) = 0.0
      !elseif (ABS(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j-1,k))<ABS(maxslope)) then
      !  slope_y(i,m,j) = G%mask2dCv(i,J) * G%mask2dCv(i,J-1) * &
      !                 0.5*(Tr(m)%t(i,j+1,k)-Tr(m)%t(i,j-1,k))
      !else
      !  slope_y(i,m,j) = G%mask2dCv(i,J) * G%mask2dCv(i,J-1) * 0.5*maxslope
      !endif
       tp = tr(m)%t(i,j+1,k) ; tc = tr(m)%t(i,j,k) ; tm = tr(m)%t(i,j-1,k)
       dmx = max( tp, tc, tm ) - tc
       dmn= tc - min( tp, tc, tm )
       slope_y(i,m,j) = g%mask2dCv(i,j)*g%mask2dCv(i,j-1) * &
           sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
    enddo ; enddo ; endif ; enddo ! End of i-, m-, & j- loops.
  endif ! usePLMslope


  ! make a copy of the tracers in case values need to be overridden for OBCs

  do j=g%jsd,g%jed ; do m=1,ntr ; do i=g%isd,g%ied
    t_tmp(i,m,j) = tr(m)%t(i,j,k)
  enddo ; enddo ; enddo

  ! loop through open boundaries and recalculate flux terms
  if (associated(obc)) then ; if (obc%OBC_pe) then
     do n=1,obc%number_of_segments
       segment=>obc%segment(n)
       if (.not. associated(segment%tr_Reg)) cycle
       do i=is,ie
         if (segment%is_N_or_S) then
           if (i>=segment%HI%isd .and. i<=segment%HI%ied) then
              j = segment%HI%JsdB
              do m = 1,ntr ! replace tracers with OBC values
                if (associated(segment%tr_Reg%Tr(m)%tres)) then
                   if (segment%direction == obc_direction_s) then
                      t_tmp(i,m,j) = segment%tr_Reg%Tr(m)%tres(i,j,k)
                   else
                      t_tmp(i,m,j+1) = segment%tr_Reg%Tr(m)%tres(i,j,k)
                   endif
                else
                   if (segment%direction == obc_direction_s) then
                      t_tmp(i,m,j) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
                   else
                      t_tmp(i,m,j+1) = segment%tr_Reg%Tr(m)%OBC_inflow_conc
                   endif
                endif
              enddo
              do m = 1,ntr ! Apply update tracer values for slope calculation
                do j=segment%HI%JsdB-1,segment%HI%JsdB+1
                  tp = t_tmp(i,m,j+1) ; tc = t_tmp(i,m,j) ; tm = t_tmp(i,m,j-1)
                  dmx = max( tp, tc, tm ) - tc
                  dmn= tc - min( tp, tc, tm )
                  slope_y(i,m,j) = g%mask2dCv(i,j)*g%mask2dCv(i,j-1) * &
                       sign( min(0.5*abs(tp-tm), 2.0*dmx, 2.0*dmn), tp-tm )
                enddo
              enddo
           endif
         endif ! is_N_S
       enddo ! i-loop
     enddo ! segment loop
  endif; endif

  ! Calculate the j-direction fluxes of each tracer, using as much
  ! the minimum of the remaining mass flux (vhr) and the half the mass
  ! in the cell plus whatever part of its half of the mass flux that
  ! the flux through the other side does not require.
  do j=js-1,je ; if (domore_v(j,k)) then
    domore_v(j,k) = .false.

    do i=is,ie
      if ((vhr(i,j,k) == 0.0) .or. &
          ((vhr(i,j,k) < 0.0) .and. (hprev(i,j+1,k) <= tiny_h)) .or. &
          ((vhr(i,j,k) > 0.0) .and. (hprev(i,j,k) <= tiny_h)) ) then
        vhh(i,j) = 0.0
        cfl(i) = 0.0
      elseif (vhr(i,j,k) < 0.0) then
        hup = hprev(i,j+1,k) - g%areaT(i,j+1)*min_h
        hlos = max(0.0, vhr(i,j+1,k))
        if ((((hup - hlos) + vhr(i,j,k)) < 0.0) .and. &
            ((0.5*hup + vhr(i,j,k)) < 0.0)) then
          vhh(i,j) = min(-0.5*hup, -hup+hlos, 0.0)
          domore_v(j,k) = .true.
        else
          vhh(i,j) = vhr(i,j,k)
        endif
        cfl(i) = - vhh(i,j) / hprev(i,j+1,k)  ! CFL is positive
      else
        hup = hprev(i,j,k) - g%areaT(i,j)*min_h
        hlos = max(0.0, -vhr(i,j-1,k))
        if ((((hup - hlos) - vhr(i,j,k)) < 0.0) .and. &
            ((0.5*hup - vhr(i,j,k)) < 0.0)) then
          vhh(i,j) = max(0.5*hup, hup-hlos, 0.0)
          domore_v(j,k) = .true.
        else
          vhh(i,j) = vhr(i,j,k)
        endif
        cfl(i) = vhh(i,j) / hprev(i,j,k)  ! CFL is positive
      endif
    enddo

    if (useppm) then
      do m=1,ntr ; do i=is,ie
        ! centre cell depending on upstream direction
        if (vhh(i,j) >= 0.0) then
          j_up = j
        else
          j_up = j + 1
        endif

        ! Implementation of PPM-H3
        tp = t_tmp(i,m,j_up+1) ; tc = t_tmp(i,m,j_up) ; tm = t_tmp(i,m,j_up-1)

        if (usehuynh) then
          al = ( 5.*tc + ( 2.*tm - tp ) )/6. ! H3 estimate
          al = max( min(tc,tm), al) ; al = min( max(tc,tm), al) ! Bound
          ar = ( 5.*tc + ( 2.*tp - tm ) )/6. ! H3 estimate
          ar = max( min(tc,tp), ar) ; ar = min( max(tc,tp), ar) ! Bound
        else
          al = 0.5 * ((tm + tc) + (slope_y(i,m,j_up-1) - slope_y(i,m,j_up)) / 3.)
          ar = 0.5 * ((tc + tp) + (slope_y(i,m,j_up) - slope_y(i,m,j_up+1)) / 3.)
        endif

        da = ar - al ; ma = 0.5*( ar + al )
        if (g%mask2dCv(i,j_up)*g%mask2dCv(i,j_up-1)*(tp-tc)*(tc-tm) <= 0.) then
          al = tc ; ar = tc ! PCM for local extremum and bounadry cells
        elseif ( da*(tc-ma) > (da*da)/6. ) then
          al = 3.*tc - 2.*ar
        elseif ( da*(tc-ma) < - (da*da)/6. ) then
          ar = 3.*tc - 2.*al
        endif

        a6 = 6.*tc - 3. * (ar + al) ! Curvature

        if (vhh(i,j) >= 0.0) then
          flux_y(i,m,j) = vhh(i,j)*( ar - 0.5 * cfl(i) * ( &
               ( ar - al ) - a6 * ( 1. - 2./3. * cfl(i) ) ) )
        else
          flux_y(i,m,j) = vhh(i,j)*( al + 0.5 * cfl(i) * ( &
               ( ar - al ) + a6 * ( 1. - 2./3. * cfl(i) ) ) )
        endif
      enddo ; enddo
    else ! PLM
      do m=1,ntr ; do i=is,ie
        if (vhh(i,j) >= 0.0) then
          ! Indirect implementation of PLM
         !aL = Tr(m)%t(i,j,k) - 0.5 * slope_y(i,m,j)
         !aR = Tr(m)%t(i,j,k) + 0.5 * slope_y(i,m,j)
         !flux_y(i,m,J) = vhh(i,J)*( aR - 0.5 * (aR-aL) * CFL(i) )
          ! Alternative implementation of PLM
          tc = t_tmp(i,m,j)
          flux_y(i,m,j) = vhh(i,j)*( tc + 0.5 * slope_y(i,m,j) * ( 1. - cfl(i) ) )
        else
          ! Indirect implementation of PLM
         !aL = Tr(m)%t(i,j+1,k) - 0.5 * slope_y(i,m,j+1)
         !aR = Tr(m)%t(i,j+1,k) + 0.5 * slope_y(i,m,j+1)
         !flux_y(i,m,J) = vhh(i,J)*( aL + 0.5 * (aR-aL) * CFL(i) )
          ! Alternative implementation of PLM
          tc = t_tmp(i,m,j+1)
          flux_y(i,m,j) = vhh(i,j)*( tc - 0.5 * slope_y(i,m,j+1) * ( 1. - cfl(i) ) )
        endif
      enddo ; enddo
    endif ! usePPM

    if (associated(obc)) then ; if (obc%OBC_pe) then
      if (obc%specified_v_BCs_exist_globally .or. obc%open_v_BCs_exist_globally) then
        do n=1,obc%number_of_segments
          segment=>obc%segment(n)
          if (.not. segment%specified) cycle
          if (.not. associated(segment%tr_Reg)) cycle
          if (obc%segment(n)%is_N_or_S) then
            if (j >= segment%HI%JsdB .and. j<= segment%HI%JedB) then
              do i=segment%HI%isd,segment%HI%ied
                ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
                ! Now changing to simply fixed inflows.
                if ((vhr(i,j,k) > 0.0) .and. (segment%direction == obc_direction_s) .or. &
                   (vhr(i,j,k) < 0.0) .and. (segment%direction == obc_direction_n)) then
                  vhh(i,j) = vhr(i,j,k)
                  do m=1,ntr
                    if (associated(segment%tr_Reg%Tr(m)%t)) then
                      flux_y(i,m,j) = vhh(i,j)*obc%segment(n)%tr_Reg%Tr(m)%tres(i,j,k)
                    else ; flux_y(i,m,j) = vhh(i,j)*obc%segment(n)%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
                  enddo
                endif
              enddo
            endif
          endif
        enddo
      endif

      if (obc%open_v_BCs_exist_globally) then
        do n=1,obc%number_of_segments
          segment=>obc%segment(n)
          if (segment%specified) cycle
          if (.not. associated(segment%tr_Reg)) cycle
          if (segment%is_N_or_S .and. (j >= segment%HI%JsdB .and. j<= segment%HI%JedB)) then
            do i=segment%HI%isd,segment%HI%ied
              ! Tracer fluxes are set to prescribed values only for inflows from masked areas.
              if ((vhr(i,j,k) > 0.0) .and. (g%mask2dT(i,j) < 0.5) .or. &
                  (vhr(i,j,k) < 0.0) .and. (g%mask2dT(i,j+1) < 0.5)) then
                vhh(i,j) = vhr(i,j,k)
                do m=1,ntr
                  if (associated(segment%tr_Reg%Tr(m)%t)) then
                    flux_y(i,m,j) = vhh(i,j)*segment%tr_Reg%Tr(m)%tres(i,j,k)
                  else ; flux_y(i,m,j) = vhh(i,j)*segment%tr_Reg%Tr(m)%OBC_inflow_conc ; endif
                enddo
              endif
            enddo
          endif
        enddo
      endif
    endif; endif

  else ! not domore_v.
    do i=is,ie ; vhh(i,j) = 0.0 ; enddo
    do m=1,ntr ; do i=is,ie ; flux_y(i,m,j) = 0.0 ; enddo ; enddo
  endif ; enddo ! End of j-loop

  do j=js-1,je ; do i=is,ie
    vhr(i,j,k) = vhr(i,j,k) - vhh(i,j)
    if (abs(vhr(i,j,k)) < vh_neglect(i,j)) vhr(i,j,k) = 0.0
  enddo ; enddo

  ! Calculate new tracer concentration in each cell after accounting
  ! for the j-direction fluxes.
  do j=js,je ; if (do_j_tr(j)) then
    do i=is,ie
      if ((vhh(i,j) /= 0.0) .or. (vhh(i,j-1) /= 0.0)) then
        do_i(i,j) = .true.
        hlst(i) = hprev(i,j,k)
        hprev(i,j,k) = max(hprev(i,j,k) - (vhh(i,j) - vhh(i,j-1)), 0.0)
        if (hprev(i,j,k) <= 0.0) then ; do_i(i,j) = .false.
        elseif (hprev(i,j,k) < h_neglect*g%areaT(i,j)) then
          hlst(i) = hlst(i) + (h_neglect*g%areaT(i,j) - hprev(i,j,k))
          ihnew(i) = 1.0 / (h_neglect*g%areaT(i,j))
        else ;  ihnew(i) = 1.0 / hprev(i,j,k) ; endif
      else ; do_i(i,j) = .false. ; endif
    enddo

    ! update tracer and save some diagnostics
    do m=1,ntr
      do i=is,ie ; if (do_i(i,j)) then
        tr(m)%t(i,j,k) = (tr(m)%t(i,j,k) * hlst(i) - &
                          (flux_y(i,m,j) - flux_y(i,m,j-1))) * ihnew(i)
      endif ; enddo

      ! diagnostics
      if (associated(tr(m)%ad_y)) then ; do i=is,ie ; if (do_i(i,j)) then
        tr(m)%ad_y(i,j,k) = tr(m)%ad_y(i,j,k) + flux_y(i,m,j)*idt
      endif ; enddo ; endif

      ! diagnose convergence of flux_y and add to convergence of flux_x.
      ! division by areaT to get into W/m2 for heat and kg/(s*m2) for salt.
      if (associated(tr(m)%advection_xy)) then
        do i=is,ie ; if (do_i(i,j)) then
          tr(m)%advection_xy(i,j,k) = tr(m)%advection_xy(i,j,k) - (flux_y(i,m,j) - flux_y(i,m,j-1))* idt * &
                                          g%IareaT(i,j)
        endif ; enddo
      endif

    enddo
  endif ; enddo ! End of j-loop.

  ! compute ad2d_y diagnostic outside above j-loop so as to make the summation ordered when OMP is active.

  !$OMP ordered
  do j=js,je ; if (do_j_tr(j)) then
    do m=1,ntr
      if (associated(tr(m)%ad2d_y)) then ; do i=is,ie ; if (do_i(i,j)) then
        tr(m)%ad2d_y(i,j) = tr(m)%ad2d_y(i,j) + flux_y(i,m,j)*idt
      endif ; enddo ; endif
    enddo
  endif ; enddo ! End of j-loop.
  !$OMP end ordered

end subroutine advect_y

!> Initialize lateral tracer advection module
subroutine tracer_advect_init(Time, G, US, param_file, diag, CS)
  type(time_type), target, intent(in)    :: time        !< current model time
  type(ocean_grid_type),   intent(in)    :: g           !< ocean grid structure
  type(unit_scale_type),   intent(in)    :: us          !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file  !< open file to parse for model parameters
  type(diag_ctrl), target, intent(inout) :: diag        !< regulates diagnostic output
  type(tracer_advect_cs),  pointer       :: cs          !< module control structure

  integer, save :: init_calls = 0

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

  if (associated(cs)) then
    call mom_error(warning, "tracer_advect_init called with associated control structure.")
    return
  endif
  allocate(cs)

  cs%diag => diag

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "DT", cs%dt, fail_if_missing=.true., &
          desc="The (baroclinic) dynamics time step.", units="s", scale=us%s_to_T)
  call get_param(param_file, mdl, "DEBUG", cs%debug, default=.false.)
  call get_param(param_file, mdl, "TRACER_ADVECTION_SCHEME", mesg, &
          desc="The horizontal transport scheme for tracers:\n"//&
          "  PLM    - Piecewise Linear Method\n"//&
          "  PPM:H3 - Piecewise Parabolic Method (Huyhn 3rd order)\n"// &
          "  PPM    - Piecewise Parabolic Method (Colella-Woodward)" &
          , default='PLM')
  select case (trim(mesg))
    case ("PLM")
      cs%usePPM = .false.
    case ("PPM:H3")
      cs%usePPM = .true.
      cs%useHuynh = .true.
    case ("PPM")
      cs%usePPM = .true.
      cs%useHuynh = .false.
    case default
      call mom_error(fatal, "MOM_tracer_advect, tracer_advect_init: "//&
           "Unknown TRACER_ADVECTION_SCHEME = "//trim(mesg))
  end select

  id_clock_advect = cpu_clock_id('(Ocean advect tracer)', grain=clock_module)
  id_clock_pass = cpu_clock_id('(Ocean tracer halo updates)', grain=clock_routine)
  id_clock_sync = cpu_clock_id('(Ocean tracer global synch)', grain=clock_routine)

end subroutine tracer_advect_init

!> Close the tracer advection module
subroutine tracer_advect_end(CS)
  type(tracer_advect_cs), pointer :: cs  !< module control structure

  if (associated(cs)) deallocate(cs)

end subroutine tracer_advect_end


!> \namespace mom_tracer_advect
!!
!!    This program contains the subroutines that advect tracers
!!  horizontally (i.e. along layers).
!!
!! \section section_mom_advect_intro
!!
!!  * advect_tracer advects tracer concentrations using a combination
!!  of the modified flux advection scheme from Easter (Mon. Wea. Rev.,
!!  1993) with tracer distributions given by the monotonic modified
!!  van Leer scheme proposed by Lin et al. (Mon. Wea. Rev., 1994).
!!  This scheme conserves the total amount of tracer while avoiding
!!  spurious maxima and minima of the tracer concentration.  If a
!!  higher order accuracy scheme is needed, suggest monotonic
!!  piecewise parabolic method, as described in Carpenter et al.
!!  (MWR, 1990).
!!
!!  * advect_tracer has 4 arguments, described below. This
!!  subroutine determines the volume of a layer in a grid cell at the
!!  previous instance when the tracer concentration was changed, so
!!  it is essential that the volume fluxes should be correct.  It is
!!  also important that the tracer advection occurs before each
!!  calculation of the diabatic forcing.

end module mom_tracer_advect