MOM_tracer_diabatic.F90

!> This module contains routines that implement physical fluxes of tracers (e.g. due
!! to surface fluxes or mixing). These are intended to be called from call_tracer_column_fns
!! in the MOM_tracer_flow_control module.
module mom_tracer_diabatic
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_grid,             only : ocean_grid_type
use mom_verticalgrid,     only : verticalgrid_type
use mom_forcing_type,     only : forcing
use mom_error_handler,    only : mom_error, fatal, warning
 
implicit none ; private
 
#include <MOM_memory.h>
public tracer_vertdiff
public applytracerboundaryfluxesinout
 
contains
 
!> This subroutine solves a tridiagonal equation for the final tracer
!! concentrations after the dual-entrainments, and possibly sinking or surface
!! and bottom sources, are applied.  The sinking is implemented with an
!! fully implicit upwind advection scheme.
subroutine tracer_vertdiff(h_old, ea, eb, dt, tr, G, GV, &
                           sfc_flux, btm_flux, btm_reservoir, sink_rate, convert_flux_in)
  type(ocean_grid_type),                     intent(in)    :: g      !< ocean grid structure
  type(verticalgrid_type),                   intent(in)    :: gv     !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: h_old  !< layer thickness before entrainment
                                                                     !! [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: ea     !< amount of fluid entrained from the layer
                                                                     !! above [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(in)    :: eb     !< amount of fluid entrained from the layer
                                                                     !! below [H ~> m or kg m-2]
  real, dimension(SZI_(G),SZJ_(G),SZK_(GV)), intent(inout) :: tr     !< tracer concentration in concentration units [CU]
  real,                                      intent(in)    :: dt     !< amount of time covered by this call [s]
  real, dimension(SZI_(G),SZJ_(G)), optional,intent(in)    :: sfc_flux !< surface flux of the tracer [CU kg m-2 s-1]
  real, dimension(SZI_(G),SZJ_(G)), optional,intent(in)    :: btm_flux !< The (negative upward) bottom flux of the
                                                                     !! tracer [CU kg m-2 s-1]
  real, dimension(SZI_(G),SZJ_(G)), optional,intent(inout) :: btm_reservoir !< amount of tracer in a bottom reservoir
                                                                     !! [CU kg m-2]; formerly [CU m]
  real,                             optional,intent(in)    :: sink_rate !< rate at which the tracer sinks [m s-1]
  logical,                          optional,intent(in)    :: convert_flux_in !< True if the specified sfc_flux needs
                                                                     !! to be integrated in time
 
  ! local variables
  real :: sink_dist !< The distance the tracer sinks in a time step [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G)) :: &
    sfc_src, &      !< The time-integrated surface source of the tracer [CU H ~> CU m or CU kg m-2].
    btm_src         !< The time-integrated bottom source of the tracer [CU H ~> CU m or CU kg m-2].
  real, dimension(SZI_(G)) :: &
    b1, &           !< b1 is used by the tridiagonal solver [H-1 ~> m-1 or m2 kg-1].
    d1              !! d1=1-c1 is used by the tridiagonal solver, nondimensional.
  real :: c1(szi_(g),szk_(gv))    !< c1 is used by the tridiagonal solver [nondim].
  real :: h_minus_dsink(szi_(g),szk_(gv)) !< The layer thickness minus the
                    !! difference in sinking rates across the layer [H ~> m or kg m-2].
                    !! By construction, 0 <= h_minus_dsink < h_work.
  real :: sink(szi_(g),szk_(gv)+1) !< The tracer's sinking distances at the
                    !! interfaces, limited to prevent characteristics from
                    !! crossing within a single timestep [H ~> m or kg m-2].
  real :: b_denom_1 !< The first term in the denominator of b1 [H ~> m or kg m-2].
  real :: h_tr      !< h_tr is h at tracer points with a h_neglect added to
                    !! ensure positive definiteness [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 :: convert_flux = .true.
 
 
  integer :: i, j, k, is, ie, js, je, nz
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke
 
  if (nz == 1) then
    call mom_error(warning, "MOM_tracer_diabatic.F90, tracer_vertdiff called "//&
                            "with only one vertical level")
    return
  endif
 
  if (present(convert_flux_in)) convert_flux = convert_flux_in
  h_neglect = gv%H_subroundoff
  sink_dist = 0.0
  if (present(sink_rate)) sink_dist = (dt*sink_rate) * gv%m_to_H
!$OMP parallel default(none) shared(is,ie,js,je,sfc_src,btm_src,sfc_flux,dt,G,GV,btm_flux, &
!$OMP                               sink_rate,btm_reservoir,nz,sink_dist,ea,      &
!$OMP                               h_old,convert_flux,h_neglect,eb,tr) &
!$OMP                       private(sink,h_minus_dsink,b_denom_1,b1,d1,h_tr,c1)
!$OMP do
  do j=js,je; do i=is,ie ; sfc_src(i,j) = 0.0 ; btm_src(i,j) = 0.0 ; enddo ; enddo
  if (present(sfc_flux)) then
    if (convert_flux) then
!$OMP do
      do j = js, je; do i = is,ie
        sfc_src(i,j) = (sfc_flux(i,j)*dt) * gv%kg_m2_to_H
      enddo ; enddo
    else
!$OMP do
      do j = js, je; do i = is,ie
        sfc_src(i,j) = sfc_flux(i,j)
      enddo ; enddo
    endif
  endif
  if (present(btm_flux)) then
    if (convert_flux) then
!$OMP do
      do j = js, je; do i = is,ie
        btm_src(i,j) = (btm_flux(i,j)*dt) * gv%kg_m2_to_H
      enddo ; enddo
    else
!$OMP do
      do j = js, je; do i = is,ie
        btm_src(i,j) = btm_flux(i,j)
      enddo ; enddo
    endif
  endif
 
  if (present(sink_rate)) then
!$OMP do
    do j=js,je
      ! Find the sinking rates at all interfaces, limiting them if necesary
      ! so that the characteristics do not cross within a timestep.
      !   If a non-constant sinking rate were used, that would be incorprated
      ! here.
      if (present(btm_reservoir)) then
        do i=is,ie ; sink(i,nz+1) = sink_dist ; enddo
        do k=2,nz ; do i=is,ie
          sink(i,k) = sink_dist ; h_minus_dsink(i,k) = h_old(i,j,k)
        enddo ; enddo
      else
        do i=is,ie ; sink(i,nz+1) = 0.0 ; enddo
        ! Find the limited sinking distance at the interfaces.
        do k=nz,2,-1 ; do i=is,ie
          if (sink(i,k+1) >= sink_dist) then
            sink(i,k) = sink_dist
            h_minus_dsink(i,k) = h_old(i,j,k) + (sink(i,k+1) - sink(i,k))
          elseif (sink(i,k+1) + h_old(i,j,k) < sink_dist) then
            sink(i,k) = sink(i,k+1) + h_old(i,j,k)
            h_minus_dsink(i,k) = 0.0
          else
            sink(i,k) = sink_dist
            h_minus_dsink(i,k) = (h_old(i,j,k) + sink(i,k+1)) - sink(i,k)
          endif
        enddo ; enddo
      endif
      do i=is,ie
        sink(i,1) = 0.0 ; h_minus_dsink(i,1) = (h_old(i,j,1) + sink(i,2))
      enddo
 
      ! Now solve the tridiagonal equation for the tracer concentrations.
      do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
        b_denom_1 = h_minus_dsink(i,1) + ea(i,j,1) + h_neglect
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
        d1(i) = b_denom_1 * b1(i)
        h_tr = h_old(i,j,1) + h_neglect
        tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j)
      endif ; enddo
      do k=2,nz-1 ; do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
        c1(i,k) = eb(i,j,k-1) * b1(i)
        b_denom_1 = h_minus_dsink(i,k) + d1(i) * (ea(i,j,k) + sink(i,k)) + &
                    h_neglect
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
        d1(i) = b_denom_1 * b1(i)
        h_tr = h_old(i,j,k) + h_neglect
        tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + &
                             (ea(i,j,k) + sink(i,k)) * tr(i,j,k-1))
      endif ; enddo ; enddo
      do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
        c1(i,nz) = eb(i,j,nz-1) * b1(i)
        b_denom_1 = h_minus_dsink(i,nz) + d1(i) * (ea(i,j,nz) + sink(i,nz)) + &
                    h_neglect
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,nz))
        h_tr = h_old(i,j,nz) + h_neglect
        tr(i,j,nz) = b1(i) * ((h_tr * tr(i,j,nz) + btm_src(i,j)) + &
                              (ea(i,j,nz) + sink(i,nz)) * tr(i,j,nz-1))
      endif ; enddo
      if (present(btm_reservoir)) then ; do i=is,ie ; if (g%mask2dT(i,j)>0.5) then
        btm_reservoir(i,j) = btm_reservoir(i,j) + &
                             (sink(i,nz+1)*tr(i,j,nz)) * gv%H_to_kg_m2
      endif ; enddo ; endif
 
      do k=nz-1,1,-1 ; do i=is,ie ; if (g%mask2dT(i,j) > 0.5) then
        tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
      endif ; enddo ; enddo
    enddo
  else
!$OMP do
    do j=js,je
      do i=is,ie ; if (g%mask2dT(i,j) > -0.5) then
        h_tr = h_old(i,j,1) + h_neglect
        b_denom_1 = h_tr + ea(i,j,1)
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,1))
        d1(i) = h_tr * b1(i)
        tr(i,j,1) = (b1(i)*h_tr)*tr(i,j,1) + sfc_src(i,j)
       endif
      enddo
      do k=2,nz-1 ; do i=is,ie ; if (g%mask2dT(i,j) > -0.5) then
        c1(i,k) = eb(i,j,k-1) * b1(i)
        h_tr = h_old(i,j,k) + h_neglect
        b_denom_1 = h_tr + d1(i) * ea(i,j,k)
        b1(i) = 1.0 / (b_denom_1 + eb(i,j,k))
        d1(i) = b_denom_1 * b1(i)
        tr(i,j,k) = b1(i) * (h_tr * tr(i,j,k) + ea(i,j,k) * tr(i,j,k-1))
      endif ; enddo ; enddo
      do i=is,ie ; if (g%mask2dT(i,j) > -0.5) then
        c1(i,nz) = eb(i,j,nz-1) * b1(i)
        h_tr = h_old(i,j,nz) + h_neglect
        b_denom_1 = h_tr + d1(i)*ea(i,j,nz)
        b1(i) = 1.0 / ( b_denom_1 + eb(i,j,nz))
        tr(i,j,nz) = b1(i) * (( h_tr * tr(i,j,nz) + btm_src(i,j)) + &
                              ea(i,j,nz) * tr(i,j,nz-1))
      endif ; enddo
      do k=nz-1,1,-1 ; do i=is,ie ; if (g%mask2dT(i,j) > -0.5) then
        tr(i,j,k) = tr(i,j,k) + c1(i,k+1)*tr(i,j,k+1)
      endif ; enddo ; enddo
    enddo
  endif
 
!$OMP end parallel
 
end subroutine tracer_vertdiff
 
!> This routine is modeled after applyBoundaryFluxesInOut in MOM_diabatic_aux.F90
!! NOTE: Please note that in this routine sfc_flux gets set to zero to ensure that the surface
!! flux of the tracer does not get applied again during a subsequent call to tracer_vertdif
subroutine applytracerboundaryfluxesinout(G, GV, Tr, dt, fluxes, h, evap_CFL_limit, minimum_forcing_depth, &
               in_flux_optional, out_flux_optional, update_h_opt)
 
  type(ocean_grid_type),                      intent(in   ) :: g  !< Grid structure
  type(verticalgrid_type),                    intent(in   ) :: gv !< ocean vertical grid structure
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(inout) :: tr !< Tracer concentration on T-cell
  real,                                       intent(in   ) :: dt !< Time-step over which forcing is applied [s]
  type(forcing),                              intent(in   ) :: fluxes !< Surface fluxes container
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),   intent(inout) :: h  !< Layer thickness [H ~> m or kg m-2]
  real,                                       intent(in   ) :: evap_cfl_limit !< Limit on the fraction of the
                                                                  !! water that can be fluxed out of the top
                                                                  !! layer in a timestep [nondim]
  real,                                       intent(in   ) :: minimum_forcing_depth !< The smallest depth over
                                                                  !! which fluxes can be applied [m]
  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in   ) :: in_flux_optional !< The total time-integrated
                                                                  !! amount of tracer that enters with freshwater
  real, dimension(SZI_(G),SZJ_(G)), optional, intent(in) :: out_flux_optional !< The total time-integrated
                                                                  !! amount of tracer that leaves with freshwater
  logical,                          optional, intent(in) :: update_h_opt  !< Optional flag to determine whether
                                                                  !! h should be updated
 
  integer, parameter :: maxgroundings = 5
  integer :: numberofgroundings, iground(maxgroundings), jground(maxgroundings)
  real :: h_limit_fluxes, iforcingdepthscale, idt
  real :: dthickness, dtracer
  real :: fractionofforcing, hold, ithickness
  real :: rivermixconst  ! A constant used in implementing river mixing [Pa s].
  real, dimension(SZI_(G)) :: &
    netmassinout, &  ! surface water fluxes [H ~> m or kg m-2] over time step
    netmassin,    &  ! mass entering ocean surface [H ~> m or kg m-2] over a time step
    netmassout       ! mass leaving ocean surface [H ~> m or kg m-2] over a time step
 
  real, dimension(SZI_(G), SZK_(G)) :: h2d, tr2d
  real, dimension(SZI_(G),SZJ_(G))  :: in_flux  ! The total time-integrated amount of tracer!
                                                   ! that enters with freshwater
  real, dimension(SZI_(G),SZJ_(G))  :: out_flux ! The total time-integrated amount of tracer!
                                                    ! that leaves with freshwater
  real, dimension(SZI_(G))          :: in_flux_1d, out_flux_1d
  real                              :: hgrounding(maxgroundings)
  real    :: tr_in
  logical :: update_h
  integer :: i, j, is, ie, js, je, k, nz, n, nsw
  character(len=45) :: mesg
 
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = g%ke
 
  ! If no freshwater fluxes, nothing needs to be done in this routine
  if ( (.not. associated(fluxes%netMassIn)) .or. (.not. associated(fluxes%netMassOut)) ) return
 
  in_flux(:,:) = 0.0 ; out_flux(:,:) = 0.0
  if (present(in_flux_optional)) then
    do j=js,je ; do i=is,ie
      in_flux(i,j) = in_flux_optional(i,j)
    enddo ; enddo
  endif
  if (present(out_flux_optional)) then
    do j=js,je ; do i=is,ie
      out_flux(i,j) = out_flux_optional(i,j)
    enddo ; enddo
  endif
 
  if (present(update_h_opt)) then
    update_h = update_h_opt
  else
    update_h = .true.
  endif
 
  idt = 1.0/dt
  numberofgroundings = 0
 
!$OMP parallel do default(none) shared(is,ie,js,je,nz,h,Tr,G,GV,fluxes,dt,    &
!$OMP                                  IforcingDepthScale,minimum_forcing_depth, &
!$OMP                                  numberOfGroundings,iGround,jGround,update_h, &
!$OMP                                  in_flux,out_flux,hGrounding,Idt,evap_CFL_limit) &
!$OMP                          private(h2d,Tr2d,netMassInOut,netMassOut,      &
!$OMP                                  in_flux_1d,out_flux_1d,fractionOfForcing,     &
!$OMP                                  dThickness,dTracer,hOld,Ithickness,           &
!$OMP                                  netMassIn, Tr_in)
 
  ! Work in vertical slices for efficiency
  do j=js,je
 
    ! Copy state into 2D-slice arrays
    do k=1,nz ; do i=is,ie
      h2d(i,k) = h(i,j,k)
      tr2d(i,k) = tr(i,j,k)
    enddo ; enddo
 
    do i = is,ie
      in_flux_1d(i) = in_flux(i,j)
      out_flux_1d(i) = out_flux(i,j)
    enddo
    ! The surface forcing is contained in the fluxes type.
    ! We aggregate the thermodynamic forcing for a time step into the following:
    ! These should have been set and stored during a call to applyBoundaryFluxesInOut
    ! netMassIn    = net mass entering at ocean surface over a timestep
    ! netMassOut   = net mass leaving ocean surface [H ~> m or kg m-2] over a time step.
    !                netMassOut < 0 means mass leaves ocean.
 
    ! Note here that the aggregateFW flag has already been taken care of in the call to
    ! applyBoundaryFluxesInOut
    do i=is,ie
        netmassout(i) = fluxes%netMassOut(i,j)
        netmassin(i)  = fluxes%netMassIn(i,j)
    enddo
 
    ! Apply the surface boundary fluxes in three steps:
    ! A/ update concentration from mass entering the ocean
    ! B/ update concentration from mass leaving ocean.
    do i=is,ie
      if (g%mask2dT(i,j)>0.) then
 
        ! A/ Update tracer due to incoming mass flux.
        do k=1,1
 
          ! Change in state due to forcing
          dthickness = netmassin(i) ! Since we are adding mass, we can use all of it
          dtracer = 0.
 
          ! Update the forcing by the part to be consumed within the present k-layer.
          ! If fractionOfForcing = 1, then updated netMassIn, netHeat, and netSalt vanish.
          netmassin(i) = netmassin(i) - dthickness
          dtracer = dtracer + in_flux_1d(i)
          in_flux_1d(i) = 0.0
 
          ! Update state
          hold     = h2d(i,k)               ! Keep original thickness in hand
          h2d(i,k) = h2d(i,k) + dthickness  ! New thickness
          if (h2d(i,k) > 0.0) then
            ithickness  = 1.0/h2d(i,k)      ! Inverse new thickness
            ! The "if"s below avoid changing T/S by roundoff unnecessarily
            if (dthickness /= 0. .or. dtracer /= 0.) tr2d(i,k) = (hold*tr2d(i,k)+ dtracer)*ithickness
          endif
 
        enddo ! k=1,1
 
        ! B/ Update tracer from mass leaving ocean
        do k=1,nz
 
          ! Place forcing into this layer if this layer has nontrivial thickness.
          ! For layers thin relative to 1/IforcingDepthScale, then distribute
          ! forcing into deeper layers.
          iforcingdepthscale = 1. / max(gv%H_subroundoff, minimum_forcing_depth*gv%m_to_H - netmassout(i) )
          ! fractionOfForcing = 1.0, unless h2d is less than IforcingDepthScale.
          fractionofforcing = min(1.0, h2d(i,k)*iforcingdepthscale)
 
          ! In the case with (-1)*netMassOut*fractionOfForcing greater than cfl*h, we
          ! limit the forcing applied to this cell, leaving the remaining forcing to
          ! be distributed downwards.
          if (-fractionofforcing*netmassout(i) > evap_cfl_limit*h2d(i,k)) then
            fractionofforcing = -evap_cfl_limit*h2d(i,k)/netmassout(i)
          endif
 
          ! Change in state due to forcing
          dthickness = max( fractionofforcing*netmassout(i), -h2d(i,k) )
          ! Note this is slightly different to how salt is currently treated
          dtracer =  fractionofforcing*out_flux_1d(i)
 
          ! Update the forcing by the part to be consumed within the present k-layer.
          ! If fractionOfForcing = 1, then new netMassOut vanishes.
          netmassout(i) = netmassout(i) - dthickness
          out_flux_1d(i) = out_flux_1d(i) - dtracer
 
          ! Update state by the appropriate increment.
          hold     = h2d(i,k)               ! Keep original thickness in hand
          h2d(i,k) = h2d(i,k) + dthickness  ! New thickness
          if (h2d(i,k) > 0.) then
            ithickness  = 1.0/h2d(i,k) ! Inverse of new thickness
            tr2d(i,k)    = (hold*tr2d(i,k) + dtracer)*ithickness
          endif
 
        enddo ! k
 
      endif
 
      ! If anything remains after the k-loop, then we have grounded out, which is a problem.
      if (abs(in_flux_1d(i))+abs(out_flux_1d(i)) /= 0.0) then
!$OMP critical
        numberofgroundings = numberofgroundings +1
        if (numberofgroundings<=maxgroundings) then
          iground(numberofgroundings) = i ! Record i,j location of event for
          jground(numberofgroundings) = j ! warning message
          hgrounding(numberofgroundings) = abs(in_flux_1d(i))+abs(out_flux_1d(i))
        endif
!$OMP end critical
      endif
 
    enddo ! i
 
    ! Step C/ copy updated tracer concentration from the 2d slice now back into model state.
    do k=1,nz ; do i=is,ie
      tr(i,j,k) = tr2d(i,k)
    enddo ; enddo
 
    if (update_h) then
      do k=1,nz ; do i=is,ie
        h(i,j,k) = h2d(i,k)
      enddo ; enddo
    endif
 
  enddo ! j-loop finish
 
  if (numberofgroundings>0) then
    do i = 1, min(numberofgroundings, maxgroundings)
      write(mesg(1:45),'(3es15.3)') g%geoLonT( iground(i), jground(i) ), &
                             g%geoLatT( iground(i), jground(i)) , hgrounding(i)
      call mom_error(warning, "MOM_tracer_diabatic.F90, applyTracerBoundaryFluxesInOut(): "//&
                              "Tracer created. x,y,dh= "//trim(mesg), all_print=.true.)
    enddo
 
    if (numberofgroundings - maxgroundings > 0) then
      write(mesg, '(i4)') numberofgroundings - maxgroundings
      call mom_error(warning, "MOM_tracer_vertical.F90, applyTracerBoundaryFluxesInOut(): "//&
                              trim(mesg) // " groundings remaining")
    endif
  endif
 
end subroutine applytracerboundaryfluxesinout
end module mom_tracer_diabatic