MOM_continuity.F90

!> Solve the layer continuity equation.
module mom_continuity
 
! This file is part of MOM6. See LICENSE.md for the license.
 
use mom_continuity_ppm, only : continuity_ppm, continuity_ppm_init
use mom_continuity_ppm, only : continuity_ppm_stencil
use mom_continuity_ppm, only : continuity_ppm_end, continuity_ppm_cs
use mom_diag_mediator, only : time_type, diag_ctrl
use mom_error_handler, only : mom_error, mom_mesg, fatal, warning, is_root_pe
use mom_file_parser, only : get_param, log_version, param_file_type
use mom_string_functions, only : uppercase
use mom_grid, only : ocean_grid_type
use mom_open_boundary, only : ocean_obc_type
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : bt_cont_type
use mom_verticalgrid, only : verticalgrid_type
 
implicit none ; private
 
#include <MOM_memory.h>
 
public continuity, continuity_init, continuity_end, continuity_stencil
 
!> Control structure for mom_continuity
type, public :: continuity_cs ; private
  integer :: continuity_scheme !< Selects the discretization for the continuity solver.
                               !! Valid values are:
                               !! - PPM - A directionally split piecewise parabolic reconstruction solver.
                               !! The default, PPM, seems most appropriate for use with our current
                               !! time-splitting strategies.
  type(continuity_ppm_cs), pointer :: ppm_csp => null() !< Control structure for mom_continuity_ppm
end type continuity_cs
 
integer, parameter :: ppm_scheme = 1 !< Enumerated constant to select PPM
character(len=20), parameter :: ppm_string = "PPM" !< String to select PPM
 
contains
 
!> Time steps the layer thicknesses, using a monotonically limited, directionally split PPM scheme,
!! based on Lin (1994).
subroutine continuity(u, v, hin, h, uh, vh, dt, G, GV, US, CS, uhbt, vhbt, OBC, &
                      visc_rem_u, visc_rem_v, u_cor, v_cor, &
                      uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, BT_cont)
  type(ocean_grid_type),   intent(inout) :: g   !< Ocean grid structure.
  type(verticalgrid_type), intent(in)    :: gv  !< Vertical grid structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(in)    :: u   !< Zonal velocity [m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(in)    :: v   !< Meridional velocity [m s-1].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(in)    :: hin !< Initial layer thickness [H ~> m or kg m-2].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)),  &
                           intent(inout) :: h   !< Final layer thickness [H ~> m or kg m-2].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(out)   :: uh  !< Volume flux through zonal faces =
                                                !! u*h*dy [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(out)   :: vh  !< Volume flux through meridional faces =
                                                !! v*h*dx [H m2 s-1 ~> m3 s-1 or kg s-1].
  real,                    intent(in)    :: dt  !< Time increment [s].
  type(unit_scale_type),   intent(in)    :: us  !< A dimensional unit scaling type
  type(continuity_cs),     pointer       :: cs  !< Control structure for mom_continuity.
  real, dimension(SZIB_(G),SZJ_(G)), &
                 optional, intent(in)    :: uhbt !< The vertically summed volume
                                                !! flux through zonal faces [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G)), &
                 optional, intent(in)    :: vhbt !< The vertically summed volume
                                                !! flux through meridional faces [H m2 s-1 ~> m3 s-1 or kg s-1].
  type(ocean_obc_type), &
                 optional, pointer       :: obc !< Open boundaries control structure.
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(in)    :: visc_rem_u !< Both the fraction of
          !! zonal momentum that remains after a time-step of viscosity, and the fraction of a time-step's
          !! worth of a barotropic acceleration that a layer experiences after viscosity is applied.
          !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                 optional, intent(in)    :: visc_rem_v !< Both the fraction of
          !! meridional momentum that remains after a time-step of viscosity, and the fraction of a time-step's
          !! worth of a barotropic acceleration that a layer experiences after viscosity is applied.
          !! Non-dimensional between 0 (at the bottom) and 1 (far above the bottom).
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(out)   :: u_cor !< The zonal velocities that
          !! give uhbt as the depth-integrated transport [m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                 optional, intent(out)   :: v_cor !< The meridional velocities that
          !! give vhbt as the depth-integrated transport [m s-1].
  real, dimension(SZIB_(G),SZJ_(G)), &
                 optional, intent(in)    :: uhbt_aux !< A second summed zonal
          !! volume flux [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZI_(G),SZJB_(G)), &
                 optional, intent(in)    :: vhbt_aux !< A second summed meridional
          !! volume flux [H m2 s-1 ~> m3 s-1 or kg s-1].
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(inout) :: u_cor_aux !< The zonal velocities
          !! that give uhbt_aux as the depth-integrated transport [m s-1].
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                 optional, intent(inout) :: v_cor_aux !< The meridional velocities
          !! that give vhbt_aux as the depth-integrated transport [m s-1].
  type(bt_cont_type), &
                 optional, pointer       :: bt_cont !< A structure with elements
          !! that describe the effective open face areas as a function of barotropic flow.
 
  if (present(visc_rem_u) .neqv. present(visc_rem_v)) call mom_error(fatal, &
      "MOM_continuity: Either both visc_rem_u and visc_rem_v or neither"// &
       " one must be present in call to continuity.")
  if (present(u_cor) .neqv. present(v_cor)) call mom_error(fatal, &
      "MOM_continuity: Either both u_cor and v_cor or neither"// &
       " one must be present in call to continuity.")
  if (present(uhbt_aux) .neqv. present(vhbt_aux)) call mom_error(fatal, &
      "MOM_continuity: Either both uhbt_aux and uhbt_aux or neither"// &
       " one must be present in call to continuity.")
  if (present(u_cor_aux) .neqv. present(v_cor_aux)) call mom_error(fatal, &
      "MOM_continuity: Either both u_cor_aux and v_cor_aux or neither"// &
       " one must be present in call to continuity.")
  if (present(u_cor_aux) .neqv. present(uhbt_aux)) call mom_error(fatal, &
      "MOM_continuity: u_cor_aux can only be calculated if uhbt_aux is"// &
      " provided, and uhbt_aux has no other purpose.  Include both arguments"//&
      " or neither.")
 
  if (cs%continuity_scheme == ppm_scheme) then
    call continuity_ppm(u, v, hin, h, uh, vh, dt, g, gv, us, cs%PPM_CSp, uhbt, vhbt, obc, &
                        visc_rem_u, visc_rem_v, u_cor, v_cor, &
                        uhbt_aux, vhbt_aux, u_cor_aux, v_cor_aux, bt_cont)
  else
    call mom_error(fatal, "continuity: Unrecognized value of continuity_scheme")
  endif
 
end subroutine continuity
 
!> Initializes continuity_cs
subroutine continuity_init(Time, G, GV, 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(verticalgrid_type), intent(in)    :: gv         !< Vertical grid structure.
  type(param_file_type),   intent(in)    :: param_file !< Parameter file handles.
  type(diag_ctrl), target, intent(inout) :: diag       !< Diagnostics control structure.
  type(continuity_cs),     pointer       :: cs         !< Control structure for mom_continuity.
! This include declares and sets the variable "version".
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_continuity" ! This module's name.
  character(len=20)  :: tmpstr
 
  if (associated(cs)) then
    call mom_error(warning, "continuity_init called with 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, "CONTINUITY_SCHEME", tmpstr, &
                 "CONTINUITY_SCHEME selects the discretization for the "//&
                 "continuity solver. The only valid value currently is: \n"//&
                 "\t PPM - use a positive-definite (or monotonic) \n"//&
                 "\t       piecewise parabolic reconstruction solver.", &
                 default=ppm_string)
 
  tmpstr = uppercase(tmpstr) ; cs%continuity_scheme = 0
  select case (trim(tmpstr))
    case (ppm_string) ; cs%continuity_scheme = ppm_scheme
    case default
      call mom_mesg('continuity_init: CONTINUITY_SCHEME ="'//trim(tmpstr)//'"', 0)
      call mom_mesg("continuity_init: The only valid value is currently "// &
                     trim(ppm_string), 0)
      call mom_error(fatal, "continuity_init: Unrecognized setting "// &
            "#define CONTINUITY_SCHEME "//trim(tmpstr)//" found in input file.")
  end select
 
  if (cs%continuity_scheme == ppm_scheme) then
    call continuity_ppm_init(time, g, gv, param_file, diag, cs%PPM_CSp)
  endif
 
end subroutine continuity_init
 
 
!> continuity_stencil returns the continuity solver stencil size
function continuity_stencil(CS) result(stencil)
  type(continuity_cs), pointer       :: cs  !< Module's control structure.
  integer ::  stencil !< The continuity solver stencil size with the current settings.
 
  stencil = 1
 
  if (cs%continuity_scheme == ppm_scheme) then
    stencil = continuity_ppm_stencil(cs%PPM_CSp)
  endif
 
end function continuity_stencil
 
!> Destructor for continuity_cs.
subroutine continuity_end(CS)
  type(continuity_cs), pointer :: cs !< Control structure for mom_continuity.
 
  if (cs%continuity_scheme == ppm_scheme) then
    call continuity_ppm_end(cs%PPM_CSp)
  endif
 
  deallocate(cs)
 
end subroutine continuity_end
 
end module mom_continuity