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
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(-)

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.
  real    :: sal_scalar !< The constant of proportionality between sea surface
                      !! height (really it should be bottom pressure) anomalies
                      !! and bottom geopotential anomalies.
  integer :: nc       !< The number of tidal constituents in use.
  real, dimension(MAX_CONSTITUENTS) :: &
    freq, &           !< The frequency of a tidal constituent [s-1].
    phase0, &         !< The phase of a tidal constituent at time 0, in radians.
    amp, &            !< The amplitude of a tidal constituent at time 0 [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, pointer, dimension(:,:,:) :: &
    sin_struct => null(), &    !< The sine and cosine based structures that can
    cos_struct => null(), &    !< be associated with the astronomical forcing.
    cosphasesal => null(), &   !< The cosine and sine of the phase of the
    sinphasesal => null(), &   !< self-attraction and loading amphidromes.
    ampsal => null(), &        !< The amplitude of the SAL [m].
    cosphase_prev => null(), & !< The cosine and sine of the phase of the
    sinphase_prev => null(), & !< amphidromes in the previous tidal solutions.
    amp_prev => null()         !< The amplitude of the previous tidal solution [m].
end type tidal_forcing_cs

integer :: id_clock_tides !< CPU clock for tides

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, T                                       !> Date offsets
  real, parameter :: PI = 4.0 * atan(1.0)            !> 3.14159...
  ! 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, N, fn, un)
  character (len=2), intent(in) :: constit              !> Tidal constituent to find modulation for.
  real, intent(in) :: N                                 !> Longitude of ascending node [rad].
                                                        !! Calculate using astro_longitudes_init.
  real, parameter :: RADIANS = 4.0 * atan(1.0) / 180.0  !> Converts degrees to radians.
  real, intent(out) :: &
    fn, & !> Amplitude modulation [nondim]
    un    !> Phase modulation [rad]
  select case (constit)
    case ("M2")
      fn = 1.0 - 0.037 * cos(n)
      un = -2.1 * radians * sin(n)
    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(n)
      un = -2.1 * radians * sin(n)
    case ("K2")
      fn = 1.024 + 0.286 * cos(n)
      un = -17.7 * radians * sin(n)
    case ("K1")
      fn = 1.006 + 0.115 * cos(n)
      un = -8.9 * radians * sin(n)
    case ("O1")
      fn = 1.009 + 0.187 * cos(n)
      un = 10.8 * radians * sin(n)
    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(n)
      un = 10.8 * radians * sin(n)
    case ("MF")
      fn = 1.043 + 0.414 * cos(n)
      un = -23.7 * radians * sin(n)
    case ("MM")
      fn = 1.0 - 0.130 * cos(n)
      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, 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(param_file_type), intent(in)    :: param_file !< A structure to parse for run-time parameters.
  type(tidal_forcing_cs), pointer      :: CS   !< A pointer that is set to point to the control
                                               !! structure for this module.
  ! 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, phase0_def, amp_def, love_def
  integer, dimension(3) :: tide_ref_date !< Reference date (t = 0) for tidal forcing.
  logical :: use_const  ! True if a constituent is being used.
  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.
  logical :: FAIL_IF_MISSING = .true.
! 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
  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 (associated(cs)) then
    call mom_error(warning, "tidal_forcing_init called with an associated "// &
                            "control structure.")
    return
  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, "TIDES", tides, &
                 "If true, apply tidal momentum forcing.", default=.false.)

  if (.not.tides) return

  allocate(cs)

  ! Set up the spatial structure functions for the diurnal, semidiurnal, and
  ! low-frequency tidal components.
  allocate(cs%sin_struct(isd:ied,jsd:jed,3)) ; cs%sin_struct(:,:,:) = 0.0
  allocate(cs%cos_struct(isd:ied,jsd:jed,3)) ; cs%cos_struct(:,:,:) = 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.)

  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
  ! Set the parameters for all components that are in use.
  ! 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)
  c=0
  if (use_m2) then
    c=c+1 ; cs%const_name(c) = "M2" ; cs%struct(c) = 2
    cs%love_no(c) = 0.693 ; cs%amp(c) = 0.242334
  endif

  if (use_s2) then
    c=c+1 ; cs%const_name(c) = "S2" ; cs%struct(c) = 2
    cs%love_no(c) = 0.693 ; cs%amp(c) = 0.112743
  endif

  if (use_n2) then
    c=c+1 ; cs%const_name(c) = "N2" ; cs%struct(c) = 2
    cs%love_no(c) = 0.693 ; cs%amp(c) = 0.046397
  endif

  if (use_k2) then
    c=c+1 ; cs%const_name(c) = "K2" ; cs%struct(c) = 2
    cs%love_no(c) = 0.693 ; cs%amp(c) = 0.030684
  endif

  if (use_k1) then
    c=c+1 ; cs%const_name(c) = "K1" ; cs%struct(c) = 1
    cs%love_no(c) = 0.736 ; cs%amp(c) = 0.141565
  endif

  if (use_o1) then
    c=c+1 ; cs%const_name(c) = "O1" ; cs%struct(c) = 1
    cs%love_no(c) = 0.695 ; cs%amp(c) = 0.100661
  endif

  if (use_p1) then
    c=c+1 ; cs%const_name(c) = "P1" ; cs%struct(c) = 1
    cs%love_no(c) = 0.706 ; cs%amp(c) = 0.046848
  endif

  if (use_q1) then
    c=c+1 ; cs%const_name(c) = "Q1" ; cs%struct(c) = 1
    cs%love_no(c) = 0.695 ; cs%amp(c) = 0.019273
  endif

  if (use_mf) then
    c=c+1 ; cs%const_name(c) = "MF" ; cs%struct(c) = 3
    cs%love_no(c) = 0.693 ; cs%amp(c) = 0.042041
  endif

  if (use_mm) then
    c=c+1 ; cs%const_name(c) = "MM" ; cs%struct(c) = 3
    cs%love_no(c) = 0.693 ; cs%amp(c) = 0.022191
  endif

  ! Set defaults for all included constituents
  ! and things that can be set by functions
  do c=1,nc
    cs%freq(c) = tidal_frequency(cs%const_name(c))
    freq_def(c) = cs%freq(c)
    love_def(c) = cs%love_no(c)
    amp_def(c) = cs%amp(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))
    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))
    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)
      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)
      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

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

end subroutine tidal_forcing_init

!> 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)
  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
  ! Local variables
  integer :: nf

  do nf=1,size(filenames)
    if (len_trim(filenames(nf)) == 0) cycle
    if (field_exists(filenames(nf), varname, g%Domain%mpp_domain)) then
      call mom_read_data(filenames(nf), varname, array, g%Domain)
      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), pointer     :: 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, CS, deta_tidal_deta, m_to_Z)
  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(tidal_forcing_cs),           pointer     :: CS        !< The control structure returned by a
                                                             !! previous call to tidal_forcing_init.
  real, optional,                   intent(out) :: deta_tidal_deta !< The partial derivative of
                                                             !! eta_tidal with the local value of
                                                             !! eta [nondim].
  real, optional,                   intent(in)  :: m_to_Z    !< A scaling factor from m to the units of eta.

  ! Local variables
  real :: now       ! The relative time in seconds.
  real :: amp_cosomegat, amp_sinomegat
  real :: cosomegat, sinomegat
  real :: m_Z       ! A scaling factor from m to depth units.
  real :: eta_prop  ! The nondimenional constant of proportionality beteen eta and eta_tidal.
  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

  if (.not.associated(cs)) return

  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 = 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

  if (present(deta_tidal_deta)) then
    deta_tidal_deta = eta_prop
    do j=jsq,jeq+1 ; do i=isq,ieq+1 ; eta_tidal(i,j) = 0.0 ; enddo ; enddo
  else
    do j=jsq,jeq+1 ; do i=isq,ieq+1
      eta_tidal(i,j) = eta_prop*eta(i,j)
    enddo ; enddo
  endif

  m_z = 1.0 ; if (present(m_to_z)) m_z = m_to_z

  do c=1,cs%nc
    m = cs%struct(c)
    amp_cosomegat = m_z*cs%amp(c)*cs%love_no(c) * cos(cs%freq(c)*now + cs%phase0(c))
    amp_sinomegat = m_z*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) + m_z*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) - m_z*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

  call cpu_clock_end(id_clock_tides)

end subroutine calc_tidal_forcing

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

  if (associated(cs%sin_struct)) deallocate(cs%sin_struct)
  if (associated(cs%cos_struct)) deallocate(cs%cos_struct)

  if (associated(cs%cosphasesal)) deallocate(cs%cosphasesal)
  if (associated(cs%sinphasesal)) deallocate(cs%sinphasesal)
  if (associated(cs%ampsal))      deallocate(cs%ampsal)

  if (associated(cs%cosphase_prev)) deallocate(cs%cosphase_prev)
  if (associated(cs%sinphase_prev)) deallocate(cs%sinphase_prev)
  if (associated(cs%amp_prev))      deallocate(cs%amp_prev)

  if (associated(cs)) deallocate(cs)

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.

end module mom_tidal_forcing