pseudo_salt_tracer.F90

!> A tracer package that mimics salinity
module pseudo_salt_tracer

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

use mom_debugging,     only : hchksum
use mom_diag_mediator, only : post_data, register_diag_field, safe_alloc_ptr
use mom_diag_mediator, only : diag_ctrl
use mom_error_handler, only : mom_error, fatal, warning
use mom_file_parser, only : get_param, log_param, log_version, param_file_type
use mom_forcing_type, only : forcing
use mom_grid, only : ocean_grid_type
use mom_hor_index, only : hor_index_type
use mom_io, only : file_exists, read_data, slasher, vardesc, var_desc, query_vardesc
use mom_open_boundary, only : ocean_obc_type
use mom_restart, only : query_initialized, mom_restart_cs
use mom_sponge, only : set_up_sponge_field, sponge_cs
use mom_time_manager, only : time_type
use mom_tracer_registry, only : register_tracer, tracer_registry_type
use mom_tracer_diabatic, only : tracer_vertdiff, applytracerboundaryfluxesinout
use mom_tracer_z_init, only : tracer_z_init
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : surface
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type

use coupler_types_mod, only : coupler_type_set_data, ind_csurf
use atmos_ocean_fluxes_mod, only : aof_set_coupler_flux

implicit none ; private

#include <MOM_memory.h>

public register_pseudo_salt_tracer, initialize_pseudo_salt_tracer
public pseudo_salt_tracer_column_physics, pseudo_salt_tracer_surface_state
public pseudo_salt_stock, pseudo_salt_tracer_end

!> The control structure for the pseudo-salt tracer
type, public :: pseudo_salt_tracer_cs ; private
  type(time_type), pointer :: time => null() !< A pointer to the ocean model's clock.
  type(tracer_registry_type), pointer :: tr_reg => null() !< A pointer to the MOM tracer registry
  real, pointer :: ps(:,:,:) => null()   !< The array of pseudo-salt tracer used in this
                                         !! subroutine [ppt}
  real, pointer :: diff(:,:,:) => null() !< The difference between the pseudo-salt
                                         !! tracer and the real salt [ppt].
  logical :: pseudo_salt_may_reinit = .true. !< Hard coding since this should not matter

  integer :: id_psd = -1   !< A diagnostic ID

  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() !< A pointer to the restart control structure

  type(vardesc) :: tr_desc !< A description and metadata for the pseudo-salt tracer
end type pseudo_salt_tracer_cs

contains

!> Register the pseudo-salt tracer with MOM6
function register_pseudo_salt_tracer(HI, GV, param_file, CS, tr_Reg, restart_CS)
  type(hor_index_type),       intent(in) :: hi   !< A horizontal index type structure
  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(pseudo_salt_tracer_cs),  pointer  :: cs !< The control structure returned by a previous
                                               !! call to register_pseudo_salt_tracer.
  type(tracer_registry_type), pointer    :: tr_reg !< A pointer that is set to point to the control
                                                  !! structure for the tracer advection and
                                                  !! diffusion module
  type(mom_restart_cs),       pointer    :: restart_cs !< A pointer to the restart control structure
! This subroutine is used to register tracer fields and subroutines
! to be used with MOM.

  ! Local variables
  character(len=40)  :: mdl = "pseudo_salt_tracer" ! This module's name.
  character(len=200) :: inputdir ! The directory where the input files are.
  character(len=48)  :: var_name ! The variable's name.
  character(len=3)   :: name_tag ! String for creating identifying pseudo_salt
! This include declares and sets the variable "version".
#include "version_variable.h"
  real, pointer :: tr_ptr(:,:,:) => null()
  logical :: register_pseudo_salt_tracer
  integer :: isd, ied, jsd, jed, nz, i, j
  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke

  if (associated(cs)) then
    call mom_error(warning, "register_pseudo_salt_tracer called with an "// &
                             "associated control structure.")
    return
  endif
  allocate(cs)

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

  allocate(cs%ps(isd:ied,jsd:jed,nz)) ; cs%ps(:,:,:) = 0.0
  allocate(cs%diff(isd:ied,jsd:jed,nz)) ; cs%diff(:,:,:) = 0.0

  cs%tr_desc = var_desc(trim("pseudo_salt"), "psu", &
                     "Pseudo salt passive tracer", caller=mdl)

  tr_ptr => cs%ps(:,:,:)
  call query_vardesc(cs%tr_desc, name=var_name, caller="register_pseudo_salt_tracer")
  ! Register the tracer for horizontal advection, diffusion, and restarts.
  call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, name="pseudo_salt", &
                       longname="Pseudo salt passive tracer", units="psu", &
                       registry_diags=.true., restart_cs=restart_cs, &
                       mandatory=.not.cs%pseudo_salt_may_reinit)

  cs%tr_Reg => tr_reg
  cs%restart_CSp => restart_cs
  register_pseudo_salt_tracer = .true.

end function register_pseudo_salt_tracer

!> Initialize the pseudo-salt tracer
subroutine initialize_pseudo_salt_tracer(restart, day, G, GV, h, diag, OBC, CS, &
                                  sponge_CSp, tv)
  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(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]
  type(diag_ctrl),            target, intent(in) :: diag !< A structure that is used to regulate
                                                         !! diagnostic output.
  type(ocean_obc_type),               pointer    :: obc  !< This open boundary condition type specifies
                                                         !! whether, where, and what open boundary
                                                         !! conditions are used.
  type(pseudo_salt_tracer_cs),        pointer    :: cs !< The control structure returned by a previous
                                                       !! call to register_pseudo_salt_tracer.
  type(sponge_cs),                    pointer    :: sponge_csp !< Pointer to the control structure for the sponges.
  type(thermo_var_ptrs),              intent(in) :: tv   !< A structure pointing to various thermodynamic variables
!   This subroutine initializes the tracer fields in CS%ps(:,:,:).

  ! Local variables
  character(len=16) :: name     ! A variable's name in a NetCDF file.
  character(len=72) :: longname ! The long name of that variable.
  character(len=48) :: units    ! The dimensions of the variable.
  character(len=48) :: flux_units ! The units for age tracer fluxes, either
                                ! years m3 s-1 or years kg s-1.
  logical :: ok
  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz
  integer :: isdb, iedb, jsdb, jedb

  if (.not.associated(cs)) return
  if (.not.associated(cs%diff)) return

  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

  cs%Time => day
  cs%diag => diag
  name = "pseudo_salt"

  call query_vardesc(cs%tr_desc, name=name, caller="initialize_pseudo_salt_tracer")
  if ((.not.restart) .or. (.not.query_initialized(cs%ps, name, cs%restart_CSp))) then
    do k=1,nz ; do j=jsd,jed ; do i=isd,ied
      cs%ps(i,j,k) = tv%S(i,j,k)
    enddo ; enddo ; enddo
  endif

  if (associated(obc)) then
  ! Steal from updated DOME in the fullness of time.
  endif

  cs%id_psd = register_diag_field("ocean_model", "pseudo_salt_diff", cs%diag%axesTL, &
        day, "Difference between pseudo salt passive tracer and salt tracer", "psu")

end subroutine initialize_pseudo_salt_tracer

!> Apply sources, sinks and diapycnal diffusion to the tracers in this package.
subroutine pseudo_salt_tracer_column_physics(h_old, h_new, ea, eb, fluxes, dt, G, GV, US, CS, tv, debug, &
              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   !< an array to which the amount of fluid entrained
                                              !! from the layer above during this call will be
                                              !! added [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in) :: eb   !< an array to which the amount of fluid entrained
                                              !! from the layer below during this call will be
                                              !! added [H ~> m or kg m-2].
  type(forcing),           intent(in) :: fluxes !< A structure containing pointers to thermodynamic
                                              !! and tracer forcing fields.  Unused fields have NULL ptrs.
  real,                    intent(in) :: dt   !< The amount of time covered by this call [T ~> s]
  type(unit_scale_type),   intent(in) :: us   !< A dimensional unit scaling type
  type(pseudo_salt_tracer_cs), pointer :: cs  !< The control structure returned by a previous
                                              !! call to register_pseudo_salt_tracer.
  type(thermo_var_ptrs),   intent(in) :: tv   !< A structure pointing to various thermodynamic variables
  logical,                 intent(in) :: debug !< If true calculate checksums
  real,          optional, intent(in) :: evap_cfl_limit !< Limit on the fraction of the water that can
                                              !! be fluxed out of the top layer in a timestep [nondim]
  real,          optional, intent(in) :: minimum_forcing_depth !< The smallest depth over which
                                              !! fluxes can be applied [H ~> m or kg m-2]

!   This subroutine applies diapycnal diffusion and any other column
! tracer physics or chemistry to the tracers from this file.

! 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
  real :: year, h_total, scale, htot, ih_limit
  integer :: secs, days
  integer :: i, j, k, is, ie, js, je, nz, k_max
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified

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

  if (.not.associated(cs)) return
  if (.not.associated(cs%diff)) return

  if (debug) then
    call hchksum(tv%S,"salt pre pseudo-salt vertdiff", g%HI)
    call hchksum(cs%ps,"pseudo_salt pre pseudo-salt vertdiff", g%HI)
  endif

  ! This uses applyTracerBoundaryFluxesInOut, usually in ALE mode
  if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
    do k=1,nz ; do j=js,je ; do i=is,ie
      h_work(i,j,k) = h_old(i,j,k)
    enddo ; enddo ; enddo
    call applytracerboundaryfluxesinout(g, gv, cs%ps, dt, fluxes, h_work, &
                                        evap_cfl_limit, minimum_forcing_depth, out_flux_optional=fluxes%netSalt)
    call tracer_vertdiff(h_work, ea, eb, dt, cs%ps, g, gv)
  else
    call tracer_vertdiff(h_old, ea, eb, dt, cs%ps, g, gv)
  endif

  do k=1,nz ; do j=js,je ; do i=is,ie
    cs%diff(i,j,k) = cs%ps(i,j,k)-tv%S(i,j,k)
  enddo ; enddo ; enddo

  if (debug) then
    call hchksum(tv%S,"salt post pseudo-salt vertdiff", g%HI)
    call hchksum(cs%ps,"pseudo_salt post pseudo-salt vertdiff", g%HI)
  endif

  if (cs%id_psd>0) call post_data(cs%id_psd, cs%diff, cs%diag)

end subroutine pseudo_salt_tracer_column_physics


!> Calculates the mass-weighted integral of all tracer stocks, 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 pseudo_salt_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(pseudo_salt_tracer_cs),        pointer       :: cs !< The control structure returned by a previous
                                                          !! call to register_pseudo_salt_tracer.
  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                                           :: pseudo_salt_stock !< Return value: the number of
                                                              !! stocks calculated here.

! This function calculates the mass-weighted integral of all tracer stocks,
! returning the number of stocks it has calculated.  If the stock_index
! is present, only the stock corresponding to that coded index is returned.

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

  pseudo_salt_stock = 0
  if (.not.associated(cs)) return
  if (.not.associated(cs%diff)) 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

  call query_vardesc(cs%tr_desc, name=names(1), units=units(1), caller="pseudo_salt_stock")
  units(1) = trim(units(1))//" kg"
  stocks(1) = 0.0
  do k=1,nz ; do j=js,je ; do i=is,ie
    stocks(1) = stocks(1) + cs%diff(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(1) = gv%H_to_kg_m2 * stocks(1)

  pseudo_salt_stock = 1

end function pseudo_salt_stock

!> This subroutine extracts the surface fields from this tracer package that
!! are to be shared with the atmosphere in coupled configurations.
!! This particular tracer package does not report anything back to the coupler.
subroutine pseudo_salt_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 thickness [H ~> m or kg m-2].
  type(pseudo_salt_tracer_cs), pointer  :: cs !< The control structure returned by a previous
                                              !! call to register_pseudo_salt_tracer.

  ! This particular tracer package does not report anything back to the coupler.
  ! The code that is here is just a rough guide for packages that would.

  integer :: m, is, ie, js, je, isd, ied, jsd, jed
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec
  isd = g%isd ; ied = g%ied ; jsd = g%jsd ; jed = g%jed

  if (.not.associated(cs)) return

  ! By design, this tracer package does not return any surface states.

end subroutine pseudo_salt_tracer_surface_state

!> Deallocate memory associated with this tracer package
subroutine pseudo_salt_tracer_end(CS)
  type(pseudo_salt_tracer_cs), pointer :: cs !< The control structure returned by a previous
                                              !! call to register_pseudo_salt_tracer.
  integer :: m

  if (associated(cs)) then
    if (associated(cs%ps)) deallocate(cs%ps)
    if (associated(cs%diff)) deallocate(cs%diff)
    deallocate(cs)
  endif
end subroutine pseudo_salt_tracer_end

!> \namespace pseudo_salt_tracer
!!
!!  By Andrew Shao, 2016
!!
!!    This file contains the routines necessary to model a passive
!!  tracer that uses the same boundary fluxes as salinity. At the
!!  beginning of the run, salt is set to the same as tv%S. Any
!!  deviations between this salt-like tracer and tv%S signifies a
!!  difference between how active and passive tracers are treated.

end module pseudo_salt_tracer