MOM_tidal_forcing.F90

Go to the documentation of this file.

module mom_tidal_forcing
!***********************************************************************
!*                   GNU General Public License                        *
!* This file is a part of MOM.                                         *
!*                                                                     *
!* MOM is free software; you can redistribute it and/or modify it and  *
!* are expected to follow the terms of the GNU General Public License  *
!* as published by the Free Software Foundation; either version 2 of   *
!* the License, or (at your option) any later version.                 *
!*                                                                     *
!* MOM is distributed in the hope that it will be useful, but WITHOUT  *
!* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY  *
!* or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public    *
!* License for more details.                                           *
!*                                                                     *
!* For the full text of the GNU General Public License,                *
!* write to: Free Software Foundation, Inc.,                           *
!*           675 Mass Ave, Cambridge, MA 02139, USA.                   *
!* or see:   http://www.gnu.org/licenses/gpl.html                      *
!***********************************************************************

!********+*********+*********+*********+*********+*********+*********+**
!*                                                                     *
!*  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 a line '#define TIDES' in that file.*
!*  Tidal constituents must be individually enabled with lines like    *
!*  '#define TIDE_M2'.  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.   *
!*                                                                     *
!***********************************************************************

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, 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, read_data
use mom_time_manager,  only : time_type, time_type_to_real

implicit none ; private

public calc_tidal_forcing, tidal_forcing_init, tidal_forcing_end
public tidal_forcing_sensitivity

#include <MOM_memory.h>

integer, parameter :: max_constituents = 10 ! The maximum number of tidal
                                            ! constituents that could be used.

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.
  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, in s-1.
    phase0, &         ! The phase of a tidal constituent at time 0, in radians.
    amp, &            ! The amplitude of a tidal constituent at time 0, in m.
    love_no           ! The Love number of a tidal constituent at time 0, ND.
  integer :: struct(max_constituents)
  character (len=16) :: const_name(max_constituents)

  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, in 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, in m.
end type tidal_forcing_cs

integer :: id_clock_tides

contains

!> 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.

! 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.
!
! Arguments: Time - The current model time.
!  (in)      G - The ocean's grid structure.
!  (in)      param_file - A structure indicating the open file to parse for
!                         model parameter values.
!  (in/out)  CS - A pointer that is set to point to the control structure
!                 for this module.
  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
  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 \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "from input files, specified by TIDAL_INPUT_FILE. \n"//&
                 "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 \n"//&
                 "tides to facilitate convergent iteration. \n"//&
                 "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 \n"//&
                 "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 \n"//&
                 "height (really it should be bottom pressure) anomalies \n"//&
                 "and bottom geopotential anomalies. This is only used if \n"//&
                 "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 accomodate 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

  ! 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%freq(c) = 1.4051890e-4 ; cs%struct(c) = 2
    cs%love_no(c) = 0.693 ; cs%amp(c) = 0.242334 ; cs%phase0(c) = 0.0
    freq_def(c) = cs%freq(c) ; love_def(c) = cs%love_no(c)
    amp_def(c) = cs%amp(c) ; phase0_def(c) = cs%phase0(c)
  endif

  if (use_s2) then
    c=c+1 ; cs%const_name(c) = "S2" ; cs%freq(c) = 1.4544410e-4 ; cs%struct(c) = 2
    cs%love_no(c) = 0.693 ; cs%amp(c) = 0.112743 ; cs%phase0(c) = 0.0
    freq_def(c) = cs%freq(c) ; love_def(c) = cs%love_no(c)
    amp_def(c) = cs%amp(c) ; phase0_def(c) = cs%phase0(c)
  endif

  if (use_n2) then
    c=c+1 ; cs%const_name(c) = "N2" ; cs%freq(c) = 1.3787970e-4 ; cs%struct(c) = 2
    cs%love_no(c) = 0.693 ; cs%amp(c) = 0.046397 ; cs%phase0(c) = 0.0
    freq_def(c) = cs%freq(c) ; love_def(c) = cs%love_no(c)
    amp_def(c) = cs%amp(c) ; phase0_def(c) = cs%phase0(c)
  endif

  if (use_k2) then
    c=c+1 ; cs%const_name(c) = "K2" ; cs%freq(c) = 1.4584234e-4 ; cs%struct(c) = 2
    cs%love_no(c) = 0.693 ; cs%amp(c) = 0.030684 ; cs%phase0(c) = 0.0
    freq_def(c) = cs%freq(c) ; love_def(c) = cs%love_no(c)
    amp_def(c) = cs%amp(c) ; phase0_def(c) = cs%phase0(c)
  endif

  if (use_k1) then
    c=c+1 ; cs%const_name(c) = "K1" ; cs%freq(c) = 0.7292117e-4 ; cs%struct(c) = 1
    cs%love_no(c) = 0.736 ; cs%amp(c) = 0.141565 ; cs%phase0(c) = 0.0
    freq_def(c) = cs%freq(c) ; love_def(c) = cs%love_no(c)
    amp_def(c) = cs%amp(c) ; phase0_def(c) = cs%phase0(c)
  endif

  if (use_o1) then
    c=c+1 ; cs%const_name(c) = "O1" ; cs%freq(c) = 0.6759774e-4 ; cs%struct(c) = 1
    cs%love_no(c) = 0.695 ; cs%amp(c) = 0.100661 ; cs%phase0(c) = 0.0
    freq_def(c) = cs%freq(c) ; love_def(c) = cs%love_no(c)
    amp_def(c) = cs%amp(c) ; phase0_def(c) = cs%phase0(c)
  endif

  if (use_p1) then
    c=c+1 ; cs%const_name(c) = "P1" ; cs%freq(c) = 0.7252295e-4 ; cs%struct(c) = 1
    cs%love_no(c) = 0.706 ; cs%amp(c) = 0.046848 ; cs%phase0(c) = 0.0
    freq_def(c) = cs%freq(c) ; love_def(c) = cs%love_no(c)
    amp_def(c) = cs%amp(c) ; phase0_def(c) = cs%phase0(c)
  endif

  if (use_q1) then
    c=c+1 ; cs%const_name(c) = "Q1" ; cs%freq(c) = 0.6495854e-4 ; cs%struct(c) = 1
    cs%love_no(c) = 0.695 ; cs%amp(c) = 0.019273 ; cs%phase0(c) = 0.0
    freq_def(c) = cs%freq(c) ; love_def(c) = cs%love_no(c)
    amp_def(c) = cs%amp(c) ; phase0_def(c) = cs%phase0(c)
  endif

  if (use_mf) then
    c=c+1 ; cs%const_name(c) = "MF" ; cs%freq(c) = 0.053234e-4 ; cs%struct(c) = 3
    cs%love_no(c) = 0.693 ; cs%amp(c) = 0.042041 ; cs%phase0(c) = 0.0
    freq_def(c) = cs%freq(c) ; love_def(c) = cs%love_no(c)
    amp_def(c) = cs%amp(c) ; phase0_def(c) = cs%phase0(c)
  endif

  if (use_mm) then
    c=c+1 ; cs%const_name(c) = "MM" ; cs%freq(c) = 0.026392e-4 ; cs%struct(c) = 3
    cs%love_no(c) = 0.693 ; cs%amp(c) = 0.022191 ; cs%phase0(c) = 0.0
    freq_def(c) = cs%freq(c) ; love_def(c) = cs%love_no(c)
    amp_def(c) = cs%amp(c) ; phase0_def(c) = cs%phase0(c)
  endif

  !   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. \n"//&
                   "This is only used if TIDES and TIDE_"//trim(cs%const_name(c))// &
                   " are true.", 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. \n"//&
                   "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. \n"//&
                   "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 needs a doxygen description.
subroutine find_in_files(tidal_input_files,varname,array,G)
  character(len=*),                 intent(in)  :: tidal_input_files(:)
  character(len=*),                 intent(in)  :: varname
  type(ocean_grid_type),            intent(in)  :: G    !< The ocean's grid structure
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: array

  integer :: nf

  do nf=1,size(tidal_input_files)
    if (len_trim(tidal_input_files(nf)) == 0) cycle
    if (field_exists(tidal_input_files(nf), varname, g%Domain%mpp_domain)) then
      call read_data(tidal_input_files(nf), varname, array, g%Domain%mpp_domain)
      return
    endif
  enddo

  do nf=size(tidal_input_files),1,-1
    if (file_exists(tidal_input_files(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(tidal_input_files(nf)))
    endif
  enddo

  call mom_error(fatal, "MOM_tidal_forcing.F90: Unable to find any of the "// &
                  "tidal input files, including "//trim(tidal_input_files(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.
!   This subroutine calculates returns the partial derivative of the local
! geopotential height with the input sea surface height due to self-attraction
! and loading.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 tidal_forcing_init.
!  (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 m, 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)
  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 in m.
  real, dimension(SZI_(G),SZJ_(G)), intent(out) :: eta_tidal !< The tidal forcing geopotential
                                                             !! anomalies, in 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.

!   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 m, but
! probably the input for eta should really be replaced with the column mass
! anomalies.
!
! Arguments: Time - The time for the caluculation.
!  (in)      eta - The sea surface height anomaly from a time-mean geoid in m.
!  (out)     eta_tidal - The tidal forcing geopotential anomalies, in m.
!  (in)      G - The ocean's grid structure.
!  (in)      CS - The control structure returned by a previous call to
!                 tidal_forcing_init.
!  (out)     deta_tidal_deta - the partial derivative of eta_tidal with the
!                              local value of eta, nondim.

  real :: eta_astro(szi_(g),szj_(g))
  real :: eta_SAL(szi_(g),szj_(g))
  real :: now       ! The relative time in seconds.
  real :: amp_cosomegat, amp_sinomegat
  real :: cosomegat, sinomegat
  real :: eta_prop
  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)

  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

  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

  call cpu_clock_end(id_clock_tides)

end subroutine calc_tidal_forcing

subroutine tidal_forcing_end(CS)
  type(tidal_forcing_cs), pointer   :: CS

  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

end module mom_tidal_forcing