MOM_tidal_forcing.F90

!> Tidal contributions to geopotential
module MOM_tidal_forcing

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

use MOM_cpu_clock,     only : cpu_clock_id, cpu_clock_begin, cpu_clock_end, &
                              CLOCK_MODULE, CLOCK_ROUTINE
use MOM_domains,       only : pass_var
use MOM_error_handler, only : MOM_error, MOM_mesg, FATAL, WARNING
use MOM_file_parser,   only : get_param, log_version, param_file_type
use MOM_grid,          only : ocean_grid_type
use MOM_io,            only : field_exists, file_exists, MOM_read_data
use MOM_time_manager,  only : set_date, time_type, time_type_to_real, operator(-)
use MOM_unit_scaling,  only : unit_scale_type
use MOM_spherical_harmonics, only : spherical_harmonics_init, spherical_harmonics_end, order2index, calc_lmax
use MOM_spherical_harmonics, only : spherical_harmonics_forward, spherical_harmonics_inverse
use MOM_spherical_harmonics, only : sht_CS
use MOM_load_love_numbers, only : Love_Data

implicit none ; private

public calc_tidal_forcing, tidal_forcing_init, tidal_forcing_end
public tidal_forcing_sensitivity
! MOM_open_boundary uses the following to set tides on the boundary.
public astro_longitudes_init, eq_phase, nodal_fu, tidal_frequency

#include <MOM_memory.h>

integer, parameter :: MAX_CONSTITUENTS = 10 !< The maximum number of tidal
                                            !! constituents that could be used.
!> Simple type to store astronomical longitudes used to calculate tidal phases.
type, public :: astro_longitudes
  real :: &
    s, &  !< Mean longitude of moon [rad]
    h, &  !< Mean longitude of sun [rad]
    p, &  !< Mean longitude of lunar perigee [rad]
    N     !< Longitude of ascending node [rad]
end type astro_longitudes

!> The control structure for the MOM_tidal_forcing module
type, public :: tidal_forcing_CS ; private
  logical :: use_sal_scalar !< If true, use the scalar approximation when
                      !! calculating self-attraction and loading.
  logical :: tidal_sal_from_file !< If true, Read the tidal self-attraction
                      !! and loading from input files, specified
                      !! by TIDAL_INPUT_FILE.
  logical :: use_prev_tides !< If true, use the SAL from the previous
                      !! iteration of the tides to facilitate convergence.
  logical :: use_eq_phase !< If true, tidal forcing is phase-shifted to match
                      !! equilibrium tide. Set to false if providing tidal phases
                      !! that have already been shifted by the
                      !! astronomical/equilibrium argument.
  logical :: tidal_sal_sht !< If true, use online spherical harmonics to calculate SAL
  real    :: sal_scalar !< The constant of proportionality between sea surface
                      !! height (really it should be bottom pressure) anomalies
                      !! and bottom geopotential anomalies [nondim].
  integer :: nc       !< The number of tidal constituents in use.
  real, dimension(MAX_CONSTITUENTS) :: &
    freq, &           !< The frequency of a tidal constituent [T-1 ~> s-1].
    phase0, &         !< The phase of a tidal constituent at time 0 [rad].
    amp, &            !< The amplitude of a tidal constituent at time 0 [Z ~> m].
    love_no           !< The Love number of a tidal constituent at time 0 [nondim].
  integer :: struct(MAX_CONSTITUENTS) !< An encoded spatial structure for each constituent
  character (len=16) :: const_name(MAX_CONSTITUENTS) !< The name of each constituent

  type(time_type) :: time_ref !< Reference time (t = 0) used to calculate tidal forcing.
  type(astro_longitudes) :: tidal_longitudes !< Astronomical longitudes used to calculate
                                   !! tidal phases at t = 0.
  real, allocatable :: &
    sin_struct(:,:,:), &    !< The sine and cosine based structures that can
    cos_struct(:,:,:), &    !< be associated with the astronomical forcing [nondim].
    cosphasesal(:,:,:), &   !< The cosine and sine of the phase of the
    sinphasesal(:,:,:), &   !< self-attraction and loading amphidromes.
    ampsal(:,:,:), &        !< The amplitude of the SAL [Z ~> m].
    cosphase_prev(:,:,:), & !< The cosine and sine of the phase of the
    sinphase_prev(:,:,:), & !< amphidromes in the previous tidal solutions.
    amp_prev(:,:,:)         !< The amplitude of the previous tidal solution [Z ~> m].
  type(sht_CS) :: sht       !< Spherical harmonic transforms (SHT) for SAL
  integer :: sal_sht_Nd     !< Maximum degree for SHT [nodim]
  real, allocatable :: Love_Scaling(:)      !< Love number for each SHT mode [nodim]
  real, allocatable :: Snm_Re(:), & !< Real and imaginary SHT coefficient for SHT SAL
                       Snm_Im(:)    !< [Z ~> m]
end type tidal_forcing_CS

integer :: id_clock_tides !< CPU clock for tides
integer :: id_clock_SAL   !< CPU clock for self-attraction and loading

contains

!> Finds astronomical longitudes s, h, p, and N,
!! the mean longitude of the moon, sun, lunar perigee, and ascending node, respectively,
!! at the specified reference time time_ref.
!! These formulas were obtained from
!! Kowalik and Luick, "Modern Theory and Practice of Tide Analysis and Tidal Power", 2019
!! (their Equation I.71), which are based on Schureman, 1958.
!! For simplicity, the time associated with time_ref should
!! be at midnight. These formulas also only make sense if
!! the calendar is gregorian.
subroutine astro_longitudes_init(time_ref, longitudes)
  type(time_type), intent(in) :: time_ref            !> Time to calculate longitudes for.
  type(astro_longitudes), intent(out) :: longitudes  !> Lunar and solar longitudes at time_ref.
  real :: D                                          !> Time since the reference date [days]
  real :: T                                          !> Time in Julian centuries [centuries]
  real, parameter :: PI = 4.0 * atan(1.0)            !> 3.14159... [nondim]
  ! Find date at time_ref in days since 1900-01-01
  D = time_type_to_real(time_ref - set_date(1900, 1, 1)) / (24.0 * 3600.0)
  ! Time since 1900-01-01 in Julian centuries
  ! Kowalik and Luick use 36526, but Schureman uses 36525 which I think is correct.
  T = D / 36525.0
  ! Calculate longitudes, including converting to radians on [0, 2pi)
  ! s: Mean longitude of moon
  longitudes%s = mod((277.0248 + 481267.8906 * T) + 0.0011 * (T**2), 360.0) * PI / 180.0
  ! h: Mean longitude of sun
  longitudes%h = mod((280.1895 + 36000.7689 * T) + 3.0310e-4 * (T**2), 360.0) * PI / 180.0
  ! p: Mean longitude of lunar perigee
  longitudes%p = mod((334.3853 + 4069.0340 * T) - 0.0103 * (T**2), 360.0) * PI / 180.0
  ! n: Longitude of ascending node
  longitudes%N = mod((259.1568 - 1934.142 * T) + 0.0021 * (T**2), 360.0) * PI / 180.0
end subroutine astro_longitudes_init

!> Calculates the equilibrium phase argument for the given tidal
!! constituent constit and the astronomical longitudes and the reference time.
!! These formulas follow Table I.4 of Kowalik and Luick,
!! "Modern Theory and Practice of Tide Analysis and Tidal Power", 2019.
function eq_phase(constit, longitudes)
  character (len=2), intent(in) :: constit !> Name of constituent (e.g., M2).
  type(astro_longitudes), intent(in) :: longitudes   !> Mean longitudes calculated using astro_longitudes_init
  real, parameter :: PI = 4.0 * atan(1.0)  !> 3.14159...
  real :: eq_phase                         !> The equilibrium phase argument for the constituent [rad].

  select case (constit)
    case ("M2")
      eq_phase = 2 * (longitudes%h - longitudes%s)
    case ("S2")
      eq_phase = 0.0
    case ("N2")
      eq_phase = (- 3 * longitudes%s + 2 * longitudes%h) + longitudes%p
    case ("K2")
      eq_phase = 2 * longitudes%h
    case ("K1")
      eq_phase = longitudes%h + PI / 2.0
    case ("O1")
      eq_phase = (- 2 * longitudes%s + longitudes%h) - PI / 2.0
    case ("P1")
      eq_phase = - longitudes%h - PI / 2.0
    case ("Q1")
      eq_phase = ((- 3 * longitudes%s + longitudes%h) + longitudes%p) - PI / 2.0
    case ("MF")
      eq_phase = 2 * longitudes%s
    case ("MM")
      eq_phase = longitudes%s - longitudes%p
    case default
      call MOM_error(FATAL, "eq_phase: unrecognized constituent")
  end select
end function eq_phase

!> Looks up angular frequencies for the main tidal constituents.
!! Values used here are from previous versions of MOM.
function tidal_frequency(constit)
  character (len=2), intent(in) :: constit !> Constituent to look up
  real :: tidal_frequency                  !> Angular frequency [s-1]

  select case (constit)
    case ("M2")
      tidal_frequency = 1.4051890e-4
    case ("S2")
      tidal_frequency = 1.4544410e-4
    case ("N2")
      tidal_frequency = 1.3787970e-4
    case ("K2")
      tidal_frequency = 1.4584234e-4
    case ("K1")
      tidal_frequency = 0.7292117e-4
    case ("O1")
      tidal_frequency = 0.6759774e-4
    case ("P1")
      tidal_frequency = 0.7252295e-4
    case ("Q1")
      tidal_frequency = 0.6495854e-4
    case ("MF")
      tidal_frequency = 0.053234e-4
    case ("MM")
      tidal_frequency = 0.026392e-4
    case default
      call MOM_error(FATAL, "tidal_frequency: unrecognized constituent")
  end select
end function tidal_frequency

!> Find amplitude (f) and phase (u) modulation of tidal constituents by the 18.6
!! year nodal cycle. Values here follow Table I.6 in Kowalik and Luick,
!! "Modern Theory and Practice of Tide Analysis and Tidal Power", 2019.
subroutine nodal_fu(constit, nodelon, fn, un)
  character (len=2), intent(in)  :: constit !> Tidal constituent to find modulation for.
  real,              intent(in)  :: nodelon !> Longitude of ascending node [rad], which
                                            !! can be calculated using astro_longitudes_init.
  real,              intent(out) :: fn      !> Amplitude modulation [nondim]
  real,              intent(out) :: un      !> Phase modulation [rad]

  real, parameter :: RADIANS = 4.0 * atan(1.0) / 180.0  !> Converts degrees to radians [nondim]

  select case (constit)
    case ("M2")
      fn = 1.0 - 0.037 * cos(nodelon)
      un = -2.1 * RADIANS * sin(nodelon)
    case ("S2")
      fn = 1.0  ! Solar S2 has no amplitude modulation.
      un = 0.0  ! S2 has no phase modulation.
    case ("N2")
      fn = 1.0 - 0.037 * cos(nodelon)
      un = -2.1 * RADIANS * sin(nodelon)
    case ("K2")
      fn = 1.024 + 0.286 * cos(nodelon)
      un = -17.7 * RADIANS * sin(nodelon)
    case ("K1")
      fn = 1.006 + 0.115 * cos(nodelon)
      un = -8.9 * RADIANS * sin(nodelon)
    case ("O1")
      fn = 1.009 + 0.187 * cos(nodelon)
      un = 10.8 * RADIANS * sin(nodelon)
    case ("P1")
      fn = 1.0  ! P1 has no amplitude modulation.
      un = 0.0  ! P1 has no phase modulation.
    case ("Q1")
      fn = 1.009 + 0.187 * cos(nodelon)
      un = 10.8 * RADIANS * sin(nodelon)
    case ("MF")
      fn = 1.043 + 0.414 * cos(nodelon)
      un = -23.7 * RADIANS * sin(nodelon)
    case ("MM")
      fn = 1.0 - 0.130 * cos(nodelon)
      un = 0.0  ! MM has no phase modulation.
    case default
      call MOM_error(FATAL, "nodal_fu: unrecognized constituent")
  end select

end subroutine nodal_fu

!> This subroutine allocates space for the static variables used
!! by this module.  The metrics may be effectively 0, 1, or 2-D arrays,
!! while fields like the background viscosities are 2-D arrays.
!! ALLOC is a macro defined in MOM_memory.h for allocate or nothing with
!! static memory.
subroutine tidal_forcing_init(Time, G, US, param_file, CS)
  type(time_type),        intent(in)    :: Time !< The current model time.
  type(ocean_grid_type),  intent(inout) :: G    !< The ocean's grid structure.
  type(unit_scale_type),  intent(in)    :: US   !< A dimensional unit scaling type
  type(param_file_type),  intent(in)    :: param_file !< A structure to parse for run-time parameters.
  type(tidal_forcing_CS), intent(inout) :: CS   !< Tidal forcing control struct

  ! Local variables
  real, dimension(SZI_(G), SZJ_(G)) :: &
    phase, &          ! The phase of some tidal constituent.
    lat_rad, lon_rad  ! Latitudes and longitudes of h-points in radians.
  real :: deg_to_rad
  real, dimension(MAX_CONSTITUENTS) :: freq_def ! Default frequency for each tidal constituent [s-1]
  real, dimension(MAX_CONSTITUENTS) :: phase0_def ! Default reference phase for each tidal constituent [rad]
  real, dimension(MAX_CONSTITUENTS) :: amp_def  ! Default amplitude for each tidal constituent [m]
  real, dimension(MAX_CONSTITUENTS) :: love_def ! Default love number for each constituent [nondim]
  integer, dimension(3) :: tide_ref_date !< Reference date (t = 0) for tidal forcing.
  logical :: use_M2, use_S2, use_N2, use_K2, use_K1, use_O1, use_P1, use_Q1
  logical :: use_MF, use_MM
  logical :: tides      ! True if a tidal forcing is to be used.
  ! This include declares and sets the variable "version".
# include "version_variable.h"
  character(len=40)  :: mdl = "MOM_tidal_forcing" ! This module's name.
  character(len=128) :: mesg
  character(len=200) :: tidal_input_files(4*MAX_CONSTITUENTS)
  integer :: i, j, c, is, ie, js, je, isd, ied, jsd, jed, nc
  integer :: lmax ! Total modes of the real spherical harmonics [nondim]
  real :: rhoW    ! The average density of sea water [R ~> kg m-3].
  real :: rhoE    ! The average density of Earth [R ~> kg m-3].

  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
  isd = G%isd ; ied = G%ied ; jsd = G%jsd; jed = G%jed

  ! Read all relevant parameters and write them to the model log.
  call log_version(param_file, mdl, version, "")
  call get_param(param_file, mdl, "TIDES", tides, &
                 "If true, apply tidal momentum forcing.", default=.false.)

  if (.not.tides) return

  ! Set up the spatial structure functions for the diurnal, semidiurnal, and
  ! low-frequency tidal components.
  allocate(CS%sin_struct(isd:ied,jsd:jed,3), source=0.0)
  allocate(CS%cos_struct(isd:ied,jsd:jed,3), source=0.0)
  deg_to_rad = 4.0*ATAN(1.0)/180.0
  do j=js-1,je+1 ; do i=is-1,ie+1
    lat_rad(i,j) = G%geoLatT(i,j)*deg_to_rad
    lon_rad(i,j) = G%geoLonT(i,j)*deg_to_rad
  enddo ; enddo
  do j=js-1,je+1 ; do i=is-1,ie+1
    CS%sin_struct(i,j,1) = -sin(2.0*lat_rad(i,j)) * sin(lon_rad(i,j))
    CS%cos_struct(i,j,1) =  sin(2.0*lat_rad(i,j)) * cos(lon_rad(i,j))
    CS%sin_struct(i,j,2) = -cos(lat_rad(i,j))**2 * sin(2.0*lon_rad(i,j))
    CS%cos_struct(i,j,2) =  cos(lat_rad(i,j))**2 * cos(2.0*lon_rad(i,j))
    CS%sin_struct(i,j,3) =  0.0
    CS%cos_struct(i,j,3) = (0.5-1.5*sin(lat_rad(i,j))**2)
  enddo ; enddo

  call get_param(param_file, mdl, "TIDE_M2", use_M2, &
                 "If true, apply tidal momentum forcing at the M2 "//&
                 "frequency. This is only used if TIDES is true.", &
                 default=.false.)
  call get_param(param_file, mdl, "TIDE_S2", use_S2, &
                 "If true, apply tidal momentum forcing at the S2 "//&
                 "frequency. This is only used if TIDES is true.", &
                 default=.false.)
  call get_param(param_file, mdl, "TIDE_N2", use_N2, &
                 "If true, apply tidal momentum forcing at the N2 "//&
                 "frequency. This is only used if TIDES is true.", &
                 default=.false.)
  call get_param(param_file, mdl, "TIDE_K2", use_K2, &
                 "If true, apply tidal momentum forcing at the K2 "//&
                 "frequency. This is only used if TIDES is true.", &
                 default=.false.)
  call get_param(param_file, mdl, "TIDE_K1", use_K1, &
                 "If true, apply tidal momentum forcing at the K1 "//&
                 "frequency. This is only used if TIDES is true.", &
                 default=.false.)
  call get_param(param_file, mdl, "TIDE_O1", use_O1, &
                 "If true, apply tidal momentum forcing at the O1 "//&
                 "frequency. This is only used if TIDES is true.", &
                 default=.false.)
  call get_param(param_file, mdl, "TIDE_P1", use_P1, &
                 "If true, apply tidal momentum forcing at the P1 "//&
                 "frequency. This is only used if TIDES is true.", &
                 default=.false.)
  call get_param(param_file, mdl, "TIDE_Q1", use_Q1, &
                 "If true, apply tidal momentum forcing at the Q1 "//&
                 "frequency. This is only used if TIDES is true.", &
                 default=.false.)
  call get_param(param_file, mdl, "TIDE_MF", use_MF, &
                 "If true, apply tidal momentum forcing at the MF "//&
                 "frequency. This is only used if TIDES is true.", &
                 default=.false.)
  call get_param(param_file, mdl, "TIDE_MM", use_MM, &
                 "If true, apply tidal momentum forcing at the MM "//&
                 "frequency. This is only used if TIDES is true.", &
                 default=.false.)

  ! Determine how many tidal components are to be used.
  nc = 0
  if (use_M2) nc=nc+1 ; if (use_S2) nc=nc+1
  if (use_N2) nc=nc+1 ; if (use_K2) nc=nc+1
  if (use_K1) nc=nc+1 ; if (use_O1) nc=nc+1
  if (use_P1) nc=nc+1 ; if (use_Q1) nc=nc+1
  if (use_MF) nc=nc+1 ; if (use_MM) nc=nc+1
  CS%nc = nc

  if (nc == 0) then
    call MOM_error(FATAL, "tidal_forcing_init: "// &
        "TIDES are defined, but no tidal constituents are used.")
    return
  endif

  call get_param(param_file, mdl, "TIDAL_SAL_FROM_FILE", CS%tidal_sal_from_file, &
                 "If true, read the tidal self-attraction and loading "//&
                 "from input files, specified by TIDAL_INPUT_FILE. "//&
                 "This is only used if TIDES is true.", default=.false.)
  call get_param(param_file, mdl, "USE_PREVIOUS_TIDES", CS%use_prev_tides, &
                 "If true, use the SAL from the previous iteration of the "//&
                 "tides to facilitate convergent iteration. "//&
                 "This is only used if TIDES is true.", default=.false.)
  call get_param(param_file, mdl, "TIDE_USE_SAL_SCALAR", CS%use_sal_scalar, &
                 "If true and TIDES is true, use the scalar approximation "//&
                 "when calculating self-attraction and loading.", &
                 default=.not.CS%tidal_sal_from_file)
  ! If it is being used, sal_scalar MUST be specified in param_file.
  if (CS%use_sal_scalar .or. CS%use_prev_tides) &
    call get_param(param_file, mdl, "TIDE_SAL_SCALAR_VALUE", CS%sal_scalar, &
                 "The constant of proportionality between sea surface "//&
                 "height (really it should be bottom pressure) anomalies "//&
                 "and bottom geopotential anomalies. This is only used if "//&
                 "TIDES and TIDE_USE_SAL_SCALAR are true.", units="m m-1", &
                 fail_if_missing=.true.)

  call get_param(param_file, mdl, "TIDAL_SAL_SHT", CS%tidal_sal_sht, &
                 "If true, use the online spherical harmonics method to calculate "//&
                 "self-attraction and loading term in tides.", default=.false.)

  if (nc > MAX_CONSTITUENTS) then
    write(mesg,'("Increase MAX_CONSTITUENTS in MOM_tidal_forcing.F90 to at least",I3, &
                &"to accommodate all the registered tidal constituents.")') nc
    call MOM_error(FATAL, "MOM_tidal_forcing"//mesg)
  endif

  do c=1,4*MAX_CONSTITUENTS ; tidal_input_files(c) = "" ; enddo

  if (CS%tidal_sal_from_file .or. CS%use_prev_tides) then
    call get_param(param_file, mdl, "TIDAL_INPUT_FILE", tidal_input_files, &
                   "A list of input files for tidal information.",         &
                   default = "", fail_if_missing=.true.)
  endif

  call get_param(param_file, mdl, "TIDE_REF_DATE", tide_ref_date, &
                 "Year,month,day to use as reference date for tidal forcing. "//&
                 "If not specified, defaults to 0.", &
                 default=0)

  call get_param(param_file, mdl, "TIDE_USE_EQ_PHASE", CS%use_eq_phase, &
                 "Correct phases by calculating equilibrium phase arguments for TIDE_REF_DATE. ", &
                 default=.false., fail_if_missing=.false.)

  if (sum(tide_ref_date) == 0) then  ! tide_ref_date defaults to 0.
    CS%time_ref = set_date(1, 1, 1)
  else
    if (.not. CS%use_eq_phase) then
      ! Using a reference date but not using phase relative to equilibrium.
      ! This makes sense as long as either phases are overridden, or
      ! correctly simulating tidal phases is not desired.
      call MOM_mesg('Tidal phases will *not* be corrected with equilibrium arguments.')
    endif
    CS%time_ref = set_date(tide_ref_date(1), tide_ref_date(2), tide_ref_date(3))
  endif

  ! Initialize reference time for tides and find relevant lunar and solar
  ! longitudes at the reference time.
  if (CS%use_eq_phase) call astro_longitudes_init(CS%time_ref, CS%tidal_longitudes)

  ! Set the parameters for all components that are in use.
  c=0
  if (use_M2) then
    c=c+1 ; CS%const_name(c) = "M2" ; CS%struct(c) = 2
    CS%love_no(c) = 0.693 ; amp_def(c) = 0.242334 ! Default amplitude in m.
  endif

  if (use_S2) then
    c=c+1 ; CS%const_name(c) = "S2" ; CS%struct(c) = 2
    CS%love_no(c) = 0.693 ; amp_def(c) = 0.112743 ! Default amplitude in m.
  endif

  if (use_N2) then
    c=c+1 ; CS%const_name(c) = "N2" ; CS%struct(c) = 2
    CS%love_no(c) = 0.693 ; amp_def(c) = 0.046397 ! Default amplitude in m.
  endif

  if (use_K2) then
    c=c+1 ; CS%const_name(c) = "K2" ; CS%struct(c) = 2
    CS%love_no(c) = 0.693 ; amp_def(c) = 0.030684 ! Default amplitude in m.
  endif

  if (use_K1) then
    c=c+1 ; CS%const_name(c) = "K1" ; CS%struct(c) = 1
    CS%love_no(c) = 0.736 ; amp_def(c) = 0.141565 ! Default amplitude in m.
  endif

  if (use_O1) then
    c=c+1 ; CS%const_name(c) = "O1" ; CS%struct(c) = 1
    CS%love_no(c) = 0.695 ; amp_def(c) = 0.100661 ! Default amplitude in m.
  endif

  if (use_P1) then
    c=c+1 ; CS%const_name(c) = "P1" ; CS%struct(c) = 1
    CS%love_no(c) = 0.706 ; amp_def(c) = 0.046848 ! Default amplitude in m.
  endif

  if (use_Q1) then
    c=c+1 ; CS%const_name(c) = "Q1" ; CS%struct(c) = 1
    CS%love_no(c) = 0.695 ; amp_def(c) = 0.019273 ! Default amplitude in m.
  endif

  if (use_MF) then
    c=c+1 ; CS%const_name(c) = "MF" ; CS%struct(c) = 3
    CS%love_no(c) = 0.693 ; amp_def(c) = 0.042041 ! Default amplitude in m.
  endif

  if (use_MM) then
    c=c+1 ; CS%const_name(c) = "MM" ; CS%struct(c) = 3
    CS%love_no(c) = 0.693 ; amp_def(c) = 0.022191 ! Default amplitude in m.
  endif

  ! Set defaults for all included constituents
  ! and things that can be set by functions
  do c=1,nc
    freq_def(c) = tidal_frequency(CS%const_name(c))
    love_def(c) = CS%love_no(c)
    CS%phase0(c) = 0.0
    if (CS%use_eq_phase) then
      phase0_def(c) = eq_phase(CS%const_name(c), CS%tidal_longitudes)
    else
      phase0_def(c) = 0.0
    endif
  enddo

  !  Parse the input file to potentially override the default values for the
  ! frequency, amplitude and initial phase of each constituent, and log the
  ! values that are actually used.
  do c=1,nc
    call get_param(param_file, mdl, "TIDE_"//trim(CS%const_name(c))//"_FREQ", CS%freq(c), &
                   "Frequency of the "//trim(CS%const_name(c))//" tidal constituent. "//&
                   "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// &
                   " are true, or if OBC_TIDE_N_CONSTITUENTS > 0 and "//trim(CS%const_name(c))// &
                   " is in OBC_TIDE_CONSTITUENTS.", units="s-1", default=freq_def(c), scale=US%T_to_s)
    call get_param(param_file, mdl, "TIDE_"//trim(CS%const_name(c))//"_AMP", CS%amp(c), &
                   "Amplitude of the "//trim(CS%const_name(c))//" tidal constituent. "//&
                   "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// &
                   " are true.", units="m", default=amp_def(c), scale=US%m_to_Z)
    call get_param(param_file, mdl, "TIDE_"//trim(CS%const_name(c))//"_PHASE_T0", CS%phase0(c), &
                   "Phase of the "//trim(CS%const_name(c))//" tidal constituent at time 0. "//&
                   "This is only used if TIDES and TIDE_"//trim(CS%const_name(c))// &
                   " are true.", units="radians", default=phase0_def(c))
  enddo

  if (CS%tidal_sal_from_file) then
    allocate(CS%cosphasesal(isd:ied,jsd:jed,nc))
    allocate(CS%sinphasesal(isd:ied,jsd:jed,nc))
    allocate(CS%ampsal(isd:ied,jsd:jed,nc))
    do c=1,nc
      ! Read variables with names like PHASE_SAL_M2 and AMP_SAL_M2.
      call find_in_files(tidal_input_files, "PHASE_SAL_"//trim(CS%const_name(c)), phase, G)
      call find_in_files(tidal_input_files, "AMP_SAL_"//trim(CS%const_name(c)), CS%ampsal(:,:,c), &
                         G, scale=US%m_to_Z)
      call pass_var(phase,           G%domain,complete=.false.)
      call pass_var(CS%ampsal(:,:,c),G%domain,complete=.true.)
      do j=js-1,je+1 ; do i=is-1,ie+1
        CS%cosphasesal(i,j,c) = cos(phase(i,j)*deg_to_rad)
        CS%sinphasesal(i,j,c) = sin(phase(i,j)*deg_to_rad)
      enddo ; enddo
    enddo
  endif

  if (CS%USE_PREV_TIDES) then
    allocate(CS%cosphase_prev(isd:ied,jsd:jed,nc))
    allocate(CS%sinphase_prev(isd:ied,jsd:jed,nc))
    allocate(CS%amp_prev(isd:ied,jsd:jed,nc))
    do c=1,nc
      ! Read variables with names like PHASE_PREV_M2 and AMP_PREV_M2.
      call find_in_files(tidal_input_files, "PHASE_PREV_"//trim(CS%const_name(c)), phase, G)
      call find_in_files(tidal_input_files, "AMP_PREV_"//trim(CS%const_name(c)), CS%amp_prev(:,:,c), &
                         G, scale=US%m_to_Z)
      call pass_var(phase,             G%domain,complete=.false.)
      call pass_var(CS%amp_prev(:,:,c),G%domain,complete=.true.)
      do j=js-1,je+1 ; do i=is-1,ie+1
        CS%cosphase_prev(i,j,c) = cos(phase(i,j)*deg_to_rad)
        CS%sinphase_prev(i,j,c) = sin(phase(i,j)*deg_to_rad)
      enddo ; enddo
    enddo
  endif

  if (CS%tidal_sal_sht) then
    call get_param(param_file, mdl, "TIDAL_SAL_SHT_DEGREE", CS%sal_sht_Nd, &
                   "The maximum degree of the spherical harmonics transformation used for "// &
                   "calculating the self-attraction and loading term for tides.", &
                   default=0, do_not_log=.not. CS%tidal_sal_sht)
    call get_param(param_file, mdl, "RHO_0", rhoW, &
                   "The mean ocean density used with BOUSSINESQ true to "//&
                   "calculate accelerations and the mass for conservation "//&
                   "properties, or with BOUSSINSEQ false to convert some "//&
                   "parameters from vertical units of m to kg m-2.", &
                   units="kg m-3", default=1035.0, scale=US%kg_m3_to_R, do_not_log=.True.)
    call get_param(param_file, mdl, "RHO_E", rhoE, &
                   "The mean solid earth density.  This is used for calculating the "// &
                   "self-attraction and loading term.", units="kg m-3", &
                   default=5517.0, scale=US%kg_m3_to_R, do_not_log=.not. CS%tidal_sal_sht)
    lmax = calc_lmax(CS%sal_sht_Nd)
    allocate(CS%Snm_Re(lmax)); CS%Snm_Re(:) = 0.0
    allocate(CS%Snm_Im(lmax)); CS%Snm_Im(:) = 0.0

    allocate(CS%Love_Scaling(lmax)); CS%Love_Scaling(:) = 0.0
    call calc_love_scaling(CS%sal_sht_Nd, rhoW, rhoE, CS%Love_Scaling)
    call spherical_harmonics_init(G, param_file, CS%sht)
    id_clock_SAL = cpu_clock_id('(Ocean SAL)', grain=CLOCK_ROUTINE)
  endif

  id_clock_tides = cpu_clock_id('(Ocean tides)', grain=CLOCK_MODULE)

end subroutine tidal_forcing_init

!> This subroutine calculates coefficients of the spherical harmonic modes for self-attraction and loading.
!! The algorithm is based on the SAL implementation in MPAS-ocean, which was modified by Kristin Barton from
!! routine written by K. Quinn (March 2010) and modified by M. Schindelegger (May 2017).
subroutine calc_love_scaling(nlm, rhoW, rhoE, Love_Scaling)
  integer, intent(in) :: nlm  !< Maximum spherical harmonics degree [nondim]
  real,    intent(in) :: rhoW !< The average density of sea water [R ~> kg m-3]
  real,    intent(in) :: rhoE !< The average density of Earth [R ~> kg m-3]
  real, dimension(:), intent(out) :: Love_Scaling !< Scaling factors for inverse SHT [nondim]

  ! Local variables
  real, dimension(:), allocatable :: HDat, LDat, KDat ! Love numbers converted in CF reference frames
  real :: H1, L1, K1 ! Temporary variables to store degree 1 Love numbers
  integer :: n_tot ! Size of the stored Love numbers
  integer :: n, m, l

  n_tot = size(Love_Data, dim=2)

  if (nlm+1 > n_tot) call MOM_error(FATAL, "MOM_tidal_forcing " // &
    "calc_love_scaling: maximum spherical harmonics degree is larger than " // &
    "the size of the stored Love numbers in MOM_load_love_number.")

  allocate(HDat(nlm+1), LDat(nlm+1), KDat(nlm+1))
  HDat(:) = Love_Data(2,1:nlm+1) ; LDat(:) = Love_Data(3,1:nlm+1) ; KDat(:) = Love_Data(4,1:nlm+1)

  ! Convert reference frames from CM to CF
  if (nlm > 0) then
    H1 = HDat(2) ; L1 = LDat(2) ;  K1 = KDat(2)
    HDat(2) = ( 2.0 / 3.0) * (H1 - L1)
    LDat(2) = (-1.0 / 3.0) * (H1 - L1)
    KDat(2) = (-1.0 / 3.0) * H1 - (2.0 / 3.0) * L1 - 1.0
  endif

  do m=0,nlm ; do n=m,nlm
    l = order2index(m,nlm)
    Love_Scaling(l+n-m) = (3.0 / real(2*n+1)) * (rhoW / rhoE) * (1.0 + KDat(n+1) - HDat(n+1))
  enddo ; enddo
end subroutine calc_love_scaling

!> This subroutine finds a named variable in a list of files and reads its
!! values into a domain-decomposed 2-d array
subroutine find_in_files(filenames, varname, array, G, scale)
  character(len=*), dimension(:),   intent(in)  :: filenames !< The names of the files to search for the named variable
  character(len=*),                 intent(in)  :: varname   !< The name of the variable to read
  type(ocean_grid_type),            intent(in)  :: G         !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: array     !< The array to fill with the data
  real,                   optional, intent(in)  :: scale     !< A factor by which to rescale the array.
  ! Local variables
  integer :: nf

  do nf=1,size(filenames)
    if (LEN_TRIM(filenames(nf)) == 0) cycle
    if (field_exists(filenames(nf), varname, MOM_domain=G%Domain)) then
      call MOM_read_data(filenames(nf), varname, array, G%Domain, scale=scale)
      return
    endif
  enddo

  do nf=size(filenames),1,-1
    if (file_exists(filenames(nf), G%Domain)) then
      call MOM_error(FATAL, "MOM_tidal_forcing.F90: Unable to find "// &
         trim(varname)//" in any of the tidal input files, last tried "// &
         trim(filenames(nf)))
    endif
  enddo

  call MOM_error(FATAL, "MOM_tidal_forcing.F90: Unable to find any of the "// &
                  "tidal input files, including "//trim(filenames(1)))

end subroutine find_in_files

!>   This subroutine calculates returns the partial derivative of the local
!! geopotential height with the input sea surface height due to self-attraction
!! and loading.
subroutine tidal_forcing_sensitivity(G, CS, deta_tidal_deta)
  type(ocean_grid_type),  intent(in)  :: G  !< The ocean's grid structure.
  type(tidal_forcing_CS), intent(in)  :: CS !< The control structure returned by a previous call to tidal_forcing_init.
  real,                   intent(out) :: deta_tidal_deta !< The partial derivative of eta_tidal with
                                            !! the local value of eta [nondim].

  if (CS%USE_SAL_SCALAR .and. CS%USE_PREV_TIDES) then
    deta_tidal_deta = 2.0*CS%SAL_SCALAR
  elseif (CS%USE_SAL_SCALAR .or. CS%USE_PREV_TIDES) then
    deta_tidal_deta = CS%SAL_SCALAR
  else
    deta_tidal_deta = 0.0
  endif
end subroutine tidal_forcing_sensitivity

!>   This subroutine calculates the geopotential anomalies that drive the tides,
!! including self-attraction and loading.  Optionally, it also returns the
!! partial derivative of the local geopotential height with the input sea surface
!! height.  For now, eta and eta_tidal are both geopotential heights in depth
!! units, but probably the input for eta should really be replaced with the
!! column mass anomalies.
subroutine calc_tidal_forcing(Time, eta, eta_tidal, G, US, CS)
  type(ocean_grid_type),            intent(in)  :: G         !< The ocean's grid structure.
  type(time_type),                  intent(in)  :: Time      !< The time for the caluculation.
  real, dimension(SZI_(G),SZJ_(G)), intent(in)  :: eta       !< The sea surface height anomaly from
                                                             !! a time-mean geoid [Z ~> m].
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_tidal !< The tidal forcing geopotential height
                                                             !! anomalies [Z ~> m].
  type(unit_scale_type),            intent(in)  :: US        !< A dimensional unit scaling type
  type(tidal_forcing_CS),           intent(inout)  :: CS        !< The control structure returned by a
                                                             !! previous call to tidal_forcing_init.

  ! Local variables
  real, dimension(SZI_(G),SZJ_(G))  :: eta_sal  !< SAL calculated by spherical harmonics
  real :: now       ! The relative time compared with the tidal reference [T ~> s]
  real :: amp_cosomegat, amp_sinomegat ! The tidal amplitudes times the components of phase [Z ~> m]
  real :: cosomegat, sinomegat ! The components of the phase [nondim]
  real :: eta_prop  ! The nondimenional constant of proportionality beteen eta and eta_tidal [nondim]
  integer :: i, j, c, m, is, ie, js, je, Isq, Ieq, Jsq, Jeq
  is = G%isc ; ie = G%iec ; js = G%jsc ; je = G%jec
  Isq = G%IscB ; Ieq = G%IecB ; Jsq = G%JscB ; Jeq = G%JecB

  call cpu_clock_begin(id_clock_tides)

  if (CS%nc == 0) then
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1 ; eta_tidal(i,j) = 0.0 ; enddo ; enddo
    return
  endif

  now = US%s_to_T * time_type_to_real(Time - cs%time_ref)

  if (CS%USE_SAL_SCALAR .and. CS%USE_PREV_TIDES) then
    eta_prop = 2.0*CS%SAL_SCALAR
  elseif (CS%USE_SAL_SCALAR .or. CS%USE_PREV_TIDES) then
    eta_prop = CS%SAL_SCALAR
  else
    eta_prop = 0.0
  endif

  do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
    eta_tidal(i,j) = eta_prop*eta(i,j)
  enddo ; enddo

  do c=1,CS%nc
    m = CS%struct(c)
    amp_cosomegat = CS%amp(c)*CS%love_no(c) * cos(CS%freq(c)*now + CS%phase0(c))
    amp_sinomegat = CS%amp(c)*CS%love_no(c) * sin(CS%freq(c)*now + CS%phase0(c))
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      eta_tidal(i,j) = eta_tidal(i,j) + (amp_cosomegat*CS%cos_struct(i,j,m) + &
                                         amp_sinomegat*CS%sin_struct(i,j,m))
    enddo ; enddo
  enddo

  if (CS%tidal_sal_from_file) then ; do c=1,CS%nc
    cosomegat = cos(CS%freq(c)*now)
    sinomegat = sin(CS%freq(c)*now)
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      eta_tidal(i,j) = eta_tidal(i,j) + CS%ampsal(i,j,c) * &
           (cosomegat*CS%cosphasesal(i,j,c) + sinomegat*CS%sinphasesal(i,j,c))
    enddo ; enddo
  enddo ; endif

  if (CS%USE_PREV_TIDES) then ; do c=1,CS%nc
    cosomegat = cos(CS%freq(c)*now)
    sinomegat = sin(CS%freq(c)*now)
    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      eta_tidal(i,j) = eta_tidal(i,j) - CS%SAL_SCALAR*CS%amp_prev(i,j,c) * &
          (cosomegat*CS%cosphase_prev(i,j,c) + sinomegat*CS%sinphase_prev(i,j,c))
    enddo ; enddo
  enddo ; endif

  if (CS%tidal_sal_sht) then
    eta_sal(:,:) = 0.0
    call calc_SAL_sht(eta, eta_sal, G, CS)

    do j=Jsq,Jeq+1 ; do i=Isq,Ieq+1
      eta_tidal(i,j) = eta_tidal(i,j) + eta_sal(i,j)
    enddo ; enddo
  endif
  call cpu_clock_end(id_clock_tides)

end subroutine calc_tidal_forcing

!> This subroutine calculates self-attraction and loading using the spherical harmonics method.
subroutine calc_SAL_sht(eta, eta_sal, G, CS)
  type(ocean_grid_type),   intent(in) :: G !< The ocean's grid structure.
  real, dimension(SZI_(G),SZJ_(G)), intent(in)  :: eta  !< The sea surface height anomaly from
                                                        !! a time-mean geoid [Z ~> m].
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_sal !< The sea surface height anomaly from
                                                           !! self-attraction and loading [Z ~> m].
  type(tidal_forcing_CS), intent(inout) :: CS !< Tidal forcing control struct

  ! Local variables
  integer :: n, m, l

  call cpu_clock_begin(id_clock_SAL)

  call spherical_harmonics_forward(G, CS%sht, eta, CS%Snm_Re, CS%Snm_Im, CS%sal_sht_Nd)

  ! Multiply scaling factors to each mode
  do m = 0,CS%sal_sht_Nd
    l = order2index(m, CS%sal_sht_Nd)
    do n = m,CS%sal_sht_Nd
      CS%Snm_Re(l+n-m) = CS%Snm_Re(l+n-m) * CS%Love_Scaling(l+n-m)
      CS%Snm_Im(l+n-m) = CS%Snm_Im(l+n-m) * CS%Love_Scaling(l+n-m)
    enddo
  enddo

  call spherical_harmonics_inverse(G, CS%sht, CS%Snm_Re, CS%Snm_Im, eta_sal, CS%sal_sht_Nd)

  call pass_var(eta_sal, G%domain)

  call cpu_clock_end(id_clock_SAL)
end subroutine calc_SAL_sht

!> This subroutine deallocates memory associated with the tidal forcing module.
subroutine tidal_forcing_end(CS)
  type(tidal_forcing_CS), intent(inout) :: CS !< The control structure returned by a previous call
                                              !! to tidal_forcing_init; it is deallocated here.

  if (allocated(CS%sin_struct)) deallocate(CS%sin_struct)
  if (allocated(CS%cos_struct)) deallocate(CS%cos_struct)

  if (allocated(CS%cosphasesal)) deallocate(CS%cosphasesal)
  if (allocated(CS%sinphasesal)) deallocate(CS%sinphasesal)
  if (allocated(CS%ampsal))      deallocate(CS%ampsal)

  if (allocated(CS%cosphase_prev)) deallocate(CS%cosphase_prev)
  if (allocated(CS%sinphase_prev)) deallocate(CS%sinphase_prev)
  if (allocated(CS%amp_prev))      deallocate(CS%amp_prev)

  if (CS%tidal_sal_sht) then
    if (allocated(CS%Love_Scaling)) deallocate(CS%Love_Scaling)
    if (allocated(CS%Snm_Re)) deallocate(CS%Snm_Re)
    if (allocated(CS%Snm_Im)) deallocate(CS%Snm_Im)
    call spherical_harmonics_end(CS%sht)
  endif
end subroutine tidal_forcing_end

!> \namespace tidal_forcing
!!
!! Code by Robert Hallberg, August 2005, based on C-code by Harper
!! Simmons, February, 2003, in turn based on code by Brian Arbic.
!!
!!   The main subroutine in this file calculates the total tidal
!! contribution to the geopotential, including self-attraction and
!! loading terms and the astronomical contributions.  All options
!! are selected with entries in a file that is parsed at run-time.
!! Overall tides are enabled with the run-time parameter 'TIDES=True'.
!! Tidal constituents must be individually enabled with lines like
!! 'TIDE_M2=True'.  This file has default values of amplitude,
!! frequency, Love number, and phase at time 0 for the Earth's M2,
!! S2, N2, K2, K1, O1, P1, Q1,  MF, and MM tidal constituents, but
!! the frequency, amplitude and phase ant time 0 for each constituent
!! can be changed at run time by setting variables like TIDE_M2_FREQ,
!! TIDE_M2_AMP and TIDE_M2_PHASE_T0 (for M2).
!!
!!   In addition, the approach to calculating self-attraction and
!! loading is set at run time.  The default is to use the scalar
!! approximation, with a coefficient TIDE_SAL_SCALAR_VALUE that must
!! be set in the run-time file (for global runs, 0.094 is typical).
!! Alternately, TIDAL_SAL_FROM_FILE can be set to read the SAL from
!! a file containing the results of a previous simulation. To iterate
!! the SAL to convergence, USE_PREVIOUS_TIDES may be useful (for
!! details, see Arbic et al., 2004, DSR II). With TIDAL_SAL_FROM_FILE
!! or USE_PREVIOUS_TIDES,a list of input files must be provided to
!! describe each constituent's properties from a previous solution.
!!
!!   This module also contains a method to calculate self-attraction
!! and loading using spherical harmonic transforms. The algorithm is
!! based on SAL calculation in Model for Prediction Across Scales
!! (MPAS)-Ocean developed by Los Alamos National Laboratory and
!! University of Michigan (Barton et al. (2022) and Brus et al. (2022)).
!!
!! Barton, K.N., Nairita, P., Brus, S.R., Petersen, M.R., Arbic, B.K., Engwirda, D., Roberts, A.F., Westerink, J.,
!! Wirasaet, D., and Schindelegger, M., 2022: Performance of Model for Prediction Across Scales (MPAS) Ocean as a
!! Global Barotropic Tide Model. Journal of Advances in Modeling Earth Systems, in review.
!!
!! Brus, S.R., Barton, K.N., Nairita, P., Roberts, A.F., Engwirda, D., Petersen, M.R., Arbic, B.K., Wirasaet, D.,
!! Westerink, J., and Schindelegger, M., 2022: Scalable self attraction and loading calculations for unstructured ocean
!! models. Ocean Modelling, in review.
end module MOM_tidal_forcing