advection_test_tracer.F90

!> This tracer package is used to test advection schemes
module advection_test_tracer

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

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_unit_scaling, only : unit_scale_type
use mom_variables, only : surface
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_advection_test_tracer, initialize_advection_test_tracer
public advection_test_tracer_surface_state, advection_test_tracer_end
public advection_test_tracer_column_physics, advection_test_stock

integer, parameter :: ntr = 11  !< The number of tracers in this module.

!> The control structure for the advect_test_tracer module
type, public :: advection_test_tracer_cs ; private
  integer :: ntr = ntr                 !< Number of tracers in this module
  logical :: coupled_tracers = .false. !< These tracers are not offered to the coupler.
  character(len=200) :: tracer_ic_file !< The full path to the IC file, or " " to initialize internally.
  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 :: tr(:,:,:,:) => null() !< The array of tracers used in this subroutine, in g m-3?
  real :: land_val(ntr) = -1.0 !< The value of tr used where land is masked out.
  logical :: use_sponge    !< If true, sponges may be applied somewhere in the domain.
  logical :: tracers_may_reinit !< If true, the tracers may be set up via the initialization code if
                           !! they are not found in the restart files.  Otherwise it is a fatal error
                           !! if the tracers are not found in the restart files of a restarted run.
  real :: x_origin !< Parameters describing the test functions
  real :: x_width  !< Parameters describing the test functions
  real :: y_origin !< Parameters describing the test functions
  real :: y_width  !< Parameters describing the test functions

  integer, dimension(NTR) :: ind_tr !< Indices returned by aof_set_coupler_flux if it is used and
                   !! the surface tracer concentrations are to be provided to the coupler.

  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(ntr) !< Descriptions and metadata for the tracers
end type advection_test_tracer_cs

contains

!> Register tracer fields and subroutines to be used with MOM.
function register_advection_test_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(advection_test_tracer_cs), pointer :: cs !< The control structure returned by a previous
                                                !! call to register_advection_test_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

  ! Local variables
  character(len=80)  :: name, longname
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "advection_test_tracer" ! This module's name.
  character(len=200) :: inputdir
  character(len=48)  :: flux_units ! The units for tracer fluxes, usually
                            ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1.
  real, pointer :: tr_ptr(:,:,:) => null()
  logical :: register_advection_test_tracer
  integer :: isd, ied, jsd, jed, nz, m
  isd = hi%isd ; ied = hi%ied ; jsd = hi%jsd ; jed = hi%jed ; nz = gv%ke

  if (associated(cs)) then
    call mom_error(warning, "register_advection_test_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, "")

  call get_param(param_file, mdl, "ADVECTION_TEST_X_ORIGIN", cs%x_origin, &
        "The x-coordinate of the center of the test-functions.", units="same as geoLon", default=0.)
  call get_param(param_file, mdl, "ADVECTION_TEST_Y_ORIGIN", cs%y_origin, &
        "The y-coordinate of the center of the test-functions.", units="same as geoLat", default=0.)
  call get_param(param_file, mdl, "ADVECTION_TEST_X_WIDTH", cs%x_width, &
        "The x-width of the test-functions.", units="same as geoLon", default=0.)
  call get_param(param_file, mdl, "ADVECTION_TEST_Y_WIDTH", cs%y_width, &
        "The y-width of the test-functions.", units="same as geoLat", default=0.)
  call get_param(param_file, mdl, "ADVECTION_TEST_TRACER_IC_FILE", cs%tracer_IC_file, &
                 "The name of a file from which to read the initial "//&
                 "conditions for the tracers, or blank to initialize "//&
                 "them internally.", default=" ")

  if (len_trim(cs%tracer_IC_file) >= 1) then
    call get_param(param_file, mdl, "INPUTDIR", inputdir, default=".")
    cs%tracer_IC_file = trim(slasher(inputdir))//trim(cs%tracer_IC_file)
    call log_param(param_file, mdl, "INPUTDIR/ADVECTION_TEST_TRACER_IC_FILE", &
                   cs%tracer_IC_file)
  endif
  call get_param(param_file, mdl, "SPONGE", cs%use_sponge, &
                 "If true, sponges may be applied anywhere in the domain. "//&
                 "The exact location and properties of those sponges are "//&
                 "specified from MOM_initialization.F90.", default=.false.)

  call get_param(param_file, mdl, "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 the tracers are not found in the "//&
                 "restart files of a restarted run.", default=.false.)


  allocate(cs%tr(isd:ied,jsd:jed,nz,ntr)) ; cs%tr(:,:,:,:) = 0.0

  do m=1,ntr
    if (m < 10) then ; write(name,'("tr",I1.1)') m
    else ; write(name,'("tr",I2.2)') m ; endif
    write(longname,'("Concentration of Tracer ",I2.2)') m
    cs%tr_desc(m) = var_desc(name, units="kg kg-1", longname=longname, caller=mdl)
    if (gv%Boussinesq) then ; flux_units = "kg kg-1 m3 s-1"
    else ; flux_units = "kg s-1" ; endif


    ! This is needed to force the compiler not to do a copy in the registration
    ! calls.  Curses on the designers and implementers of Fortran90.
    tr_ptr => cs%tr(:,:,:,m)
    ! Register the tracer for horizontal advection, diffusion, and restarts.
    call register_tracer(tr_ptr, tr_reg, param_file, hi, gv, &
                         name=name, longname=longname, units="kg kg-1", &
                         registry_diags=.true., flux_units=flux_units, &
                         restart_cs=restart_cs, mandatory=.not.cs%tracers_may_reinit)

    !   Set coupled_tracers to be true (hard-coded above) to provide the surface
    ! values to the coupler (if any).  This is meta-code and its arguments will
    ! currently (deliberately) give fatal errors if it is used.
    if (cs%coupled_tracers) &
      cs%ind_tr(m) = aof_set_coupler_flux(trim(name)//'_flux', &
          flux_type=' ', implementation=' ', caller="register_advection_test_tracer")
  enddo

  cs%tr_Reg => tr_reg
  cs%restart_CSp => restart_cs
  register_advection_test_tracer = .true.
end function register_advection_test_tracer

!>   Initializes the NTR tracer fields in tr(:,:,:,:) and it sets up the tracer output.
subroutine initialize_advection_test_tracer(restart, day, G, GV, h,diag, OBC, CS, &
                                            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(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(advection_test_tracer_cs),     pointer    :: cs !< The control structure returned by a previous
                                                       !! call to register_advection_test_tracer.
  type(sponge_cs),                    pointer    :: sponge_csp !< Pointer to the control structure for the sponges.

  ! Local variables
  real, allocatable :: temp(:,:,:)
  real, pointer, dimension(:,:,:) :: &
    obc_tr1_u => null(), & ! These arrays should be allocated and set to
    obc_tr1_v => null()    ! specify the values of tracer 1 that should come
                           ! in through u- and v- points through the open
                           ! boundary conditions, in the same units as tr.
  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 tracer fluxes, usually
                            ! kg(tracer) kg(water)-1 m3 s-1 or kg(tracer) s-1.
  real, pointer :: tr_ptr(:,:,:) => null()
  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].
  integer :: i, j, k, is, ie, js, je, isd, ied, jsd, jed, nz, m
  integer :: isdb, iedb, jsdb, jedb
  real :: tmpx, tmpy, locx, locy

  if (.not.associated(cs)) 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
  h_neglect = gv%H_subroundoff

  cs%diag => diag
  cs%ntr = ntr
  do m=1,ntr
    call query_vardesc(cs%tr_desc(m), name=name, &
                       caller="initialize_advection_test_tracer")
    if ((.not.restart) .or. (cs%tracers_may_reinit .and. .not. &
        query_initialized(cs%tr(:,:,:,m), name, cs%restart_CSp))) then
      do k=1,nz ; do j=js,je ; do i=is,ie
        cs%tr(i,j,k,m) = 0.0
      enddo ; enddo ; enddo
      k=1 ! Square wave
      do j=js,je ; do i=is,ie
        if (abs(g%geoLonT(i,j)-cs%x_origin)<0.5*cs%x_width .and. &
            abs(g%geoLatT(i,j)-cs%y_origin)<0.5*cs%y_width) cs%tr(i,j,k,m) = 1.0
      enddo ; enddo
      k=2 ! Triangle wave
      do j=js,je ; do i=is,ie
        locx=abs(g%geoLonT(i,j)-cs%x_origin)/cs%x_width
        locy=abs(g%geoLatT(i,j)-cs%y_origin)/cs%y_width
        cs%tr(i,j,k,m) = max(0.0, 1.0-locx)*max(0.0, 1.0-locy)
      enddo ; enddo
      k=3 ! Cosine bell
      do j=js,je ; do i=is,ie
        locx=min(1.0, abs(g%geoLonT(i,j)-cs%x_origin)/cs%x_width)*(acos(0.0)*2.)
        locy=min(1.0, abs(g%geoLatT(i,j)-cs%y_origin)/cs%y_width)*(acos(0.0)*2.)
        cs%tr(i,j,k,m) = (1.0+cos(locx))*(1.0+cos(locy))*0.25
      enddo ; enddo
      k=4 ! Cylinder
      do j=js,je ; do i=is,ie
        locx=abs(g%geoLonT(i,j)-cs%x_origin)/cs%x_width
        locy=abs(g%geoLatT(i,j)-cs%y_origin)/cs%y_width
        if (locx**2+locy**2<=1.0) cs%tr(i,j,k,m) = 1.0
      enddo ; enddo
      k=5 ! Cut cylinder
      do j=js,je ; do i=is,ie
        locx=(g%geoLonT(i,j)-cs%x_origin)/cs%x_width
        locy=(g%geoLatT(i,j)-cs%y_origin)/cs%y_width
        if (locx**2+locy**2<=1.0) cs%tr(i,j,k,m) = 1.0
        if (locx>0.0.and.abs(locy)<0.2) cs%tr(i,j,k,m) = 0.0
      enddo ; enddo
    endif ! restart
  enddo


end subroutine initialize_advection_test_tracer


!>   Applies diapycnal diffusion and any other column tracer physics or chemistry to the tracers
!! from this package.  This is a simple example of a set of advected passive tracers.
subroutine advection_test_tracer_column_physics(h_old, h_new,  ea,  eb, fluxes, dt, G, GV, US, CS, &
              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(advection_test_tracer_cs), pointer :: cs !< The control structure returned by a previous
                                              !! call to register_advection_test_tracer.
  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.
! This is a simple example of a set of advected passive tracers.

! 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, dimension(SZI_(G),SZJ_(G),SZK_(G)) :: h_work ! Used so that h can be modified
  real :: b1(szi_(g))          ! b1 and c1 are variables used by the
  real :: c1(szi_(g),szk_(g))  ! tridiagonal solver.
  integer :: i, j, k, is, ie, js, je, nz, m
  is = g%isc ; ie = g%iec ; js = g%jsc ; je = g%jec ; nz = gv%ke

  if (.not.associated(cs)) return

  if (present(evap_cfl_limit) .and. present(minimum_forcing_depth)) then
    do m=1,cs%ntr
      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%tr(:,:,:,m), dt, fluxes, h_work, &
                                          evap_cfl_limit, minimum_forcing_depth)
      call tracer_vertdiff(h_work, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
    enddo
  else
    do m=1,ntr
      call tracer_vertdiff(h_old, ea, eb, dt, cs%tr(:,:,:,m), g, gv)
    enddo
  endif

end subroutine advection_test_tracer_column_physics

!> 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 advection_test_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(advection_test_tracer_cs), pointer :: cs !< The control structure returned by a previous
                                              !! call to register_advection_test_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

  if (cs%coupled_tracers) then
    do m=1,cs%ntr
      !   This call loads the surface values into the appropriate array in the
      ! coupler-type structure.
      call coupler_type_set_data(cs%tr(:,:,1,m), cs%ind_tr(m), ind_csurf, &
                   sfc_state%tr_fields, idim=(/isd, is, ie, ied/), &
                   jdim=(/jsd, js, je, jed/) )
    enddo
  endif

end subroutine advection_test_tracer_surface_state

!> Calculate 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 advection_test_stock(h, stocks, G, GV, CS, names, units, stock_index)
  type(ocean_grid_type),              intent(in)    :: g      !< The ocean's 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(verticalgrid_type),            intent(in)    :: gv     !< The ocean's vertical grid structure
  type(advection_test_tracer_cs),     pointer       :: cs     !< The control structure returned by a previous
                                                              !! call to register_advection_test_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                                           :: advection_test_stock !< the number of stocks calculated here.

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

  advection_test_stock = 0
  if (.not.associated(cs)) return
  if (cs%ntr < 1) 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

  do m=1,cs%ntr
    call query_vardesc(cs%tr_desc(m), name=names(m), units=units(m), caller="advection_test_stock")
    stocks(m) = 0.0
    do k=1,nz ; do j=js,je ; do i=is,ie
      stocks(m) = stocks(m) + cs%tr(i,j,k,m) * &
                             (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)
  enddo
  advection_test_stock = cs%ntr

end function advection_test_stock

!> Deallocate memory associated with this module
subroutine advection_test_tracer_end(CS)
  type(advection_test_tracer_cs), pointer :: cs !< The control structure returned by a previous
                                              !! call to register_advection_test_tracer.
  integer :: m

  if (associated(cs)) then
    if (associated(cs%tr)) deallocate(cs%tr)
    deallocate(cs)
  endif
end subroutine advection_test_tracer_end

end module advection_test_tracer