MOM_PressureForce.F90

!> A thin wrapper for Boussinesq/non-Boussinesq forms of the pressure force calculation.
module mom_pressureforce

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

use mom_diag_mediator, only : diag_ctrl, time_type
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_grid, only : ocean_grid_type
use mom_pressureforce_fv, only : pressureforce_fv_bouss, pressureforce_fv_nonbouss
use mom_pressureforce_fv, only : pressureforce_fv_init, pressureforce_fv_end
use mom_pressureforce_fv, only : pressureforce_fv_cs
use mom_pressureforce_mont, only : pressureforce_mont_bouss, pressureforce_mont_nonbouss
use mom_pressureforce_mont, only : pressureforce_mont_init, pressureforce_mont_end
use mom_pressureforce_mont, only : pressureforce_mont_cs
use mom_tidal_forcing, only : tidal_forcing_cs
use mom_unit_scaling, only : unit_scale_type
use mom_variables, only : thermo_var_ptrs
use mom_verticalgrid, only : verticalgrid_type
use mom_ale, only: ale_cs
implicit none ; private

#include <MOM_memory.h>

public pressureforce, pressureforce_init, pressureforce_end

!> Pressure force control structure
type, public :: pressureforce_cs ; private
  logical :: analytic_fv_pgf !< If true, use the analytic finite volume form
                             !! (Adcroft et al., Ocean Mod. 2008) of the PGF.
  !> Control structure for the analytically integrated finite volume pressure force
  type(pressureforce_fv_cs), pointer :: pressureforce_fv_csp => null()
  !> Control structure for the Montgomery potential form of pressure force
  type(pressureforce_mont_cs), pointer :: pressureforce_mont_csp => null()
end type pressureforce_cs

contains

!> A thin layer between the model and the Boussinesq and non-Boussinesq pressure force routines.
subroutine pressureforce(h, tv, PFu, PFv, G, GV, US, CS, ALE_CSp, p_atm, pbce, eta)
  type(ocean_grid_type),   intent(in)  :: g    !< The ocean's grid structure
  type(verticalgrid_type), intent(in)  :: gv   !< The ocean's vertical grid structure
  type(unit_scale_type),   intent(in)  :: us   !< A dimensional unit scaling type
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                           intent(in)  :: h    !< Layer thicknesses [H ~> m or kg m-2]
  type(thermo_var_ptrs),   intent(in)  :: tv   !< A structure pointing to various thermodynamic variables
  real, dimension(SZIB_(G),SZJ_(G),SZK_(G)), &
                           intent(out) :: pfu  !< Zonal pressure force acceleration [L T-2 ~> m s-2]
  real, dimension(SZI_(G),SZJB_(G),SZK_(G)), &
                           intent(out) :: pfv  !< Meridional pressure force acceleration [L T-2 ~> m s-2]
  type(pressureforce_cs),  pointer     :: cs   !< Pressure force control structure
  type(ale_cs),            pointer     :: ale_csp !< ALE control structure
  real, dimension(:,:), &
                 optional, pointer     :: p_atm !< The pressure at the ice-ocean or
                                               !! atmosphere-ocean interface [R L2 T-2 ~> Pa].
  real, dimension(SZI_(G),SZJ_(G),SZK_(G)), &
                 optional, intent(out) :: pbce !< The baroclinic pressure anomaly in each layer
                                               !! due to eta anomalies [L2 T-2 H-1 ~> m s-2 or m4 s-2 kg-1].
  real, dimension(SZI_(G),SZJ_(G)), &
                 optional, intent(out) :: eta  !< The bottom mass used to calculate PFu and PFv,
                                               !! [H ~> m or kg m-2], with any tidal contributions.

  if (cs%Analytic_FV_PGF) then
    if (gv%Boussinesq) then
      call pressureforce_fv_bouss(h, tv, pfu, pfv, g, gv, us, cs%PressureForce_FV_CSp, &
                                   ale_csp, p_atm, pbce, eta)
    else
      call pressureforce_fv_nonbouss(h, tv, pfu, pfv, g, gv, us, cs%PressureForce_FV_CSp, &
                                      ale_csp, p_atm, pbce, eta)
    endif
  else
    if (gv%Boussinesq) then
      call pressureforce_mont_bouss(h, tv, pfu, pfv, g, gv, us, cs%PressureForce_Mont_CSp, &
                                    p_atm, pbce, eta)
    else
      call pressureforce_mont_nonbouss(h, tv, pfu, pfv, g, gv, us, cs%PressureForce_Mont_CSp, &
                                       p_atm, pbce, eta)
    endif
  endif

end subroutine pressureforce

!> Initialize the pressure force control structure
subroutine pressureforce_init(Time, G, GV, US, param_file, diag, CS, tides_CSp)
  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(unit_scale_type),   intent(in)    :: us   !< A dimensional unit scaling type
  type(param_file_type),   intent(in)    :: param_file !< Parameter file handles
  type(diag_ctrl), target, intent(inout) :: diag !< Diagnostics control structure
  type(pressureforce_cs),  pointer       :: cs   !< Pressure force control structure
  type(tidal_forcing_cs), optional, pointer :: tides_csp !< Tide control structure
#include "version_variable.h"
  character(len=40)  :: mdl = "MOM_PressureForce" ! This module's name.

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

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "ANALYTIC_FV_PGF", cs%Analytic_FV_PGF, &
                 "If true the pressure gradient forces are calculated "//&
                 "with a finite volume form that analytically integrates "//&
                 "the equations of state in pressure to avoid any "//&
                 "possibility of numerical thermobaric instability, as "//&
                 "described in Adcroft et al., O. Mod. (2008).", default=.true.)

  if (cs%Analytic_FV_PGF) then
    call pressureforce_fv_init(time, g, gv, us, param_file, diag, &
             cs%PressureForce_FV_CSp, tides_csp)
  else
    call pressureforce_mont_init(time, g, gv, us, param_file, diag, &
             cs%PressureForce_Mont_CSp, tides_csp)
  endif

end subroutine pressureforce_init

!> Deallocate the pressure force control structure
subroutine pressureforce_end(CS)
  type(pressureforce_cs), pointer :: cs !< Pressure force control structure

  if (cs%Analytic_FV_PGF) then
    call pressureforce_fv_end(cs%PressureForce_FV_CSp)
  else
    call pressureforce_mont_end(cs%PressureForce_Mont_CSp)
  endif

  if (associated(cs)) deallocate(cs)
end subroutine pressureforce_end

!> \namespace mom_pressureforce
!!
!! This thin module provides a branch to two forms of the horizontal accelerations
!! due to pressure gradients. The two options currently available are a
!! Montgomery potential form (used in traditional isopycnal layer models), and the
!! analytic finite volume form.

end module mom_pressureforce