equilibrium_tide.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE equilibrium_tide_mod
#if defined TIDE_GENERATING_FORCES
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                John Wilkin   !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This module computes the equilibrium tide defined as the shape      !
!  the sea surface (m) would assume if it were motionless and in       !
!  equilibrium with the tide generating forces on a fluid planet.      !
!  It is used to compute the tide generation force (TGF) terms for     !
!  the pressure gradient.                                              !
!                                                                      !
!  References:                                                         !
!                                                                      !
!  Arbic, B.K., Garner, S.T., Hallberg, R.W. and Simmons, H.L., 2004:  !
!    The accuracy of surface elevations in forward global barotropic   !
!    and baroclinic tide models. Deep Sea Research Part II: Topical    !
!    Studies in Oceanography, 51(25-26), pp. 3069-3101.                !
!                                                                      !
!  Arbic, B.K., Alford, M.H., Ansong, J.K., Buijsman, M.C., Ciotti,    !
!    R.B., Farrar, J.T., Hallberg, R.W., Henze, C.E., Hill, C.N.,      !
!    Luecke, C.A. and Menemenlis, D., 2018: Primer on Global Internal  !
!    Tide and Internal Gravity Wave Continuum Modeling in HYCOM and    !
!    MITgcm. In: New Frontiers in Operational Oceanography, E.         !
!    Chassignet, A. Pascual, J. Tintore and J. Verron (Eds.), GODAE    !
!    OceanView, 307-392, doi: 10.17125/gov2018.ch13.                   !
!                                                                      !
!  Doodson, A.T. and Warburg, H.D., 1941: Admiralty Manual of Tides.   !
!    His Majesty's Stationery Office, London, UK, 270 pp.              !
!                                                                      !
!  Egbert, G.D. and Ray, R.D., 2017. Tidal prediction, J. Mar. Res.,   !
!    75(3), pp.189-237.                                                !
!                                                                      !
!=======================================================================
!
      USE mod_kinds
!
      USE exchange_2d_mod
      USE dateclock_mod,   ONLY : datenum
# ifdef DISTRIBUTE
      USE mp_exchange_mod, ONLY : mp_exchange2d
# endif
      USE mod_param,       ONLY : bounds, iadm, inlm, irpm, itlm,       &
     &                            nghostpoints
      USE mod_scalars,     ONLY : ewperiodic, nsperiodic, lnodal,       &
     &                            deg2rad, tide_start, time, rclock
!
      implicit none
!
!  Define equilibrium tide constituents structure.
!
      TYPE t_etide
        real(r8) :: Afl        ! product: amp*f*love
        real(r8) :: amp        ! amplitude (m)
        real(r8) :: chi        ! phase at Greenwich meridian (degrees)
        real(r8) :: f          ! f nodal factor (nondimensional)
        real(r8) :: love       ! tidal Love number factor
        real(r8) :: nu         ! nu nodal factor (degrees)
        real(r8) :: omega      ! frequency (1/s)
      END TYPE t_etide
!
      TYPE (T_ETIDE) :: Q1     ! Q1 component, diurnal
      TYPE (T_ETIDE) :: O1     ! O1 component, diurnal
      TYPE (T_ETIDE) :: K1     ! K1 component, diurnal
      TYPE (T_ETIDE) :: N2     ! N2 component, semi-diurnal
      TYPE (T_ETIDE) :: M2     ! M2 component, semi-diurnal
      TYPE (T_ETIDE) :: S2     ! S2 component, semi-diurnal
      TYPE (T_ETIDE) :: K2     ! N2 component, semi-diurnal
!
      PUBLIC :: equilibrium_tide
      PUBLIC :: harmonic_constituents
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE equilibrium_tide (ng, tile, model)
!***********************************************************************
!
      USE mod_grid
      USE mod_ocean
!
!  Imported variable declarations.
!
      integer, intent(in) :: ng, tile, model
!
!  Local variable declarations.
!
      character (len=*), parameter :: MyFile =                          &
     &  __FILE__
!
# include "tile.h"
!
# ifdef PROFILE
      CALL wclock_on (ng, inlm, 11, __line__, myfile)
# endif
      CALL equilibrium_tide_tile (ng, tile, model,                      &
     &                            lbi, ubi, lbj, ubj,                   &
     &                            grid(ng) % lonr,                      &
     &                            grid(ng) % Cos2Lat,                   &
     &                            grid(ng) % SinLat2,                   &
# ifdef ADJOINT
     &                            ocean(ng) % ad_eq_tide,               &
# endif
# if defined TANGENT || defined TL_IOMS
     &                            ocean(ng) % tl_eq_tide,               &
# endif
     &                            ocean(ng) % eq_tide)
# ifdef PROFILE
      CALL wclock_off (ng, inlm, 11, __line__, myfile)
# endif
!
      RETURN
      END SUBROUTINE equilibrium_tide
!
!***********************************************************************
      SUBROUTINE equilibrium_tide_tile (ng, tile, model,                &
     &                                  LBi, UBi, LBj, UBj,             &
     &                                  lonr, Cos2Lat, SinLat2,         &
# ifdef ADJOINT
     &                                  ad_eq_tide,                     &
# endif
# if defined TANGENT || defined TL_IOMS
     &                                  tl_eq_tide,                     &
# endif
     &                                  eq_tide)
!***********************************************************************
!
!  Imported variables declarations.
!
      integer, intent(in) :: ng, tile, model
      integer, intent(in) :: LBi, UBi, LBj, UBj
!
# ifdef ASSUMED_SHAPE
      real(r8), intent(in) :: lonr(LBi:,LBj:)
      real(r8), intent(in) :: SinLat2(LBi:,LBj:)          ! SIN(2*latr)
      real(r8), intent(in) :: Cos2Lat(LBi:,LBj:)          ! COS2(latr)
#  ifdef ADJOINT
      real(r8), intent(out) :: ad_eq_tide(LBi:,LBj:)
#  endif
#  if defined TANGENT || defined TL_IOMS
      real(r8), intent(out) :: tl_eq_tide(LBi:,LBj:)
#  endif
      real(r8), intent(out) :: eq_tide(LBi:,LBj:)
# else
      real(r8), intent(in) :: lonr(LBi:UBi,LBj:UBj)
      real(r8), intent(in) :: SinLat2(LBi:UBi,LBj:UBj)    ! SIN(2*latr)
      real(r8), intent(in) :: Cos2Lat(LBi:UBi,LBj:UBj)    ! COS2(latr)
#  ifdef ADJOINT
      real(r8), intent(out) :: ad_eq_tide(LBi:UBi,LBj:UBj)
#  endif
#  if defined TANGENT || defined TL_IOMS
      real(r8), intent(out) :: tl_eq_tide(LBi:UBi,LBj:UBj)
#  endif
      real(r8), intent(out) :: eq_tide(LBi:UBi,LBj:UBj)
# endif
!
!  Local variables declarations.
!
      integer :: i, j
!
      real(dp) :: t_time_day, t_time_sec
 
# include "set_bounds.h"
!
!-----------------------------------------------------------------------
!  Computes sea surface associated with the equilibrium tide
!-----------------------------------------------------------------------
!
!  Set equilibrium tide time in seconds. It can be negative if ROMS time
!  (seconds since reference date, Rclock%DateNumber) is earlier than the
!  date of zero phase (Rclock%tide_DateNumber) forcing tidal boundary
!  data. This zero phase date is set up when preparing the tidal NetCDF
!  file.
!
      t_time_sec=(rclock%DateNumber(2)+time(ng))-                       &
     &           rclock%tide_DateNumber(2)
!
!  Compute astronomial equilibrium tide (m) with diurnal and semidiurnal
!  constituents. The long-period constituents and high-order harmonics
!  are neglected.
!
      DO j=jstrt,jendt
        DO i=istrt,iendt
          eq_tide(i,j)=q1%Afl*sinlat2(i,j)*                             &
     &                 cos(q1%omega*t_time_sec+                         &
     &                     deg2rad*(lonr(i,j)+q1%chi+q1%nu))+           &
     &                 o1%Afl*sinlat2(i,j)*                             &
     &                 cos(o1%omega*t_time_sec+                         &
     &                     deg2rad*(lonr(i,j)+o1%chi+o1%nu))+           &
     &                 k1%Afl*sinlat2(i,j)*                             &
     &                 cos(k1%omega*t_time_sec+                         &
     &                     deg2rad*(lonr(i,j)+k1%chi+k1%nu))+           &
     &                 n2%Afl*cos2lat(i,j)*                             &
     &                 cos(n2%omega*t_time_sec+                         &
     &                     deg2rad*(2.0_r8*lonr(i,j)+n2%chi+n2%nu))+    &
     &                 m2%Afl*cos2lat(i,j)*                             &
     &                 cos(m2%omega*t_time_sec+                         &
     &                     deg2rad*(2.0_r8*lonr(i,j)+m2%chi+m2%nu))+    &
     &                 s2%Afl*cos2lat(i,j)*                             &
     &                 cos(s2%omega*t_time_sec+                         &
     &                     deg2rad*(2.0_r8*lonr(i,j)+s2%chi+s2%nu))+    &
     &                 k2%Afl*cos2lat(i,j)*                             &
     &                 cos(k2%omega*t_time_sec+                         &
     &                     deg2rad*(2.0_r8*lonr(i,j)+k2%chi+k2%nu))
        END DO
      END DO
!
      IF (ewperiodic(ng).or.nsperiodic(ng)) THEN
        CALL exchange_r2d_tile (ng, tile,                               &
     &                          lbi, ubi, lbj, ubj,                     &
     &                          eq_tide)
      END IF
# ifdef DISTRIBUTE
      CALL mp_exchange2d (ng, tile, inlm, 1,                            &
     &                    lbi, ubi, lbj, ubj,                           &
     &                    nghostpoints,                                 &
     &                    ewperiodic(ng), nsperiodic(ng),               &
     &                    eq_tide)
# endif
# ifdef ADJOINT
!
!  Set adjoint of equilibrium tide.
!
      IF (model.eq.iadm) THEN
        ad_eq_tide=0.0_r8
      END IF
# endif
# if defined TANGENT || defined TL_IOMS
!
!  Set tangent linear of equilibrium tide.
!
      IF ((model.eq.itlm).or.(model.eq.irpm)) THEN
        tl_eq_tide=0.0_r8
      END IF
# endif
!
      RETURN
      END SUBROUTINE equilibrium_tide_tile
!
!***********************************************************************
      SUBROUTINE harmonic_constituents (Lnodal)
!***********************************************************************
!
!  Imported variables declarations.
!
      logical,  intent(in) :: Lnodal      ! Apply lunar nodal correction
!
!  Local variables declarations.
!
      real(r8) :: N, T, h, p, s
 
      real(dp) :: Astro_DateNumber(2)     ! Jan 1, 2000 12:00:00
!
!-----------------------------------------------------------------------
!  Compute fundamental astronomical parameters. Adapted from Egbert and
!  Ray, 2017, Table 1.
!-----------------------------------------------------------------------
!
!  Set Egbert and Ray time reference for their astronomical parameters
!  (Table 1): days since 2000-01-01:12:00:00.
!
      CALL datenum (astro_datenumber, 2000, 1, 1, 12, 0, 0.0_dp)
!
!  Terrestial time (in centuries) since tide reference date number. It
!  is set-up according to the chosen "tide_start" in standard input file
!  or the tidal forcing NetCDF file if found (default, recommended).
!
!  Recall that the length of a year is 365.2425 days for the now called
!  Gregorian Calendar (corrected after 15 October 1582).
!
      t=(rclock%tide_DateNumber(1)-astro_datenumber(1))/36524.25_r8
!
!  Mean longitude of the moon (Period = tropical month).
!
      s=218.316_r8+481267.8812_r8*t
!
!  Mean longitude of the sun (Period = tropical year).
!
      h=280.466_r8+36000.7698_r8*t
!
!  Mean longitude of lunar perigee (Period = 8.85 years).
!
      p=83.353_r8+4069.0137_r8*t
!
!  Mean longitude of lunar node (Period = 18.6 years).
!
      n=-234.955_r8-1934.1363_r8*t                       ! degrees
      n=n*deg2rad                                        ! radians
!
!-----------------------------------------------------------------------
!  Compute the harmonic constituents of the equilibrium tide at
!  Greenwich. Adapted from Doodson and Warburg, 1941, Table 1.
!-----------------------------------------------------------------------
!
!  Nodal factors "f" and "nu" account for the slow modulation of the
!  tidal constituents due (principally) to the 18.6-year lunar nodal
!  cycle.
!
      IF (lnodal) THEN                                 ! f nodal factor
        o1%f=1.009_r8+0.187_r8*cos(n)-0.015_r8*cos(2.0_r8*n)
        k1%f=1.006_r8+0.115_r8*cos(n)-0.009_r8*cos(2.0_r8*n)
        m2%f=1.0_r8-0.037_r8*cos(n)
        s2%f=1.0_r8
        k2%f=1.024_r8+0.286_r8*cos(n)+0.008_r8*cos(2.0_r8*n)
      ELSE
        o1%f=1.0_r8
        k1%f=1.0_r8
        m2%f=1.0_r8
        s2%f=1.0_r8
        k2%f=1.0_r8
      END IF
      q1%f=o1%f
      n2%f=m2%f
!
      IF (lnodal) THEN                                 ! nu nodal factor
        o1%nu=10.8_r8*sin(n)-1.3_r8*sin(2.0_r8*n)
        k1%nu=-8.9_r8*sin(n)+0.7_r8*sin(2.0_r8*n)
        m2%nu=-2.1_r8*sin(n)
        s2%nu=0.0_r8
        k2%nu=-17.7_r8*sin(n)+0.7_r8*sin(2.0_r8*n)
      ELSE
        o1%nu=0.0_r8
        k1%nu=0.0_r8
        m2%nu=0.0_r8
        s2%nu=0.0_r8
        k2%nu=0.0_r8
      END IF
      q1%nu=o1%nu
      n2%nu=m2%nu
!
!  Compute tidal constituent phase "chi" (degrees) at Greenwich meridian,
!  lambda = 0.
!
      q1%chi=h-3.0_r8*s+p-90.0_r8
      o1%chi=h-2.0_r8*s-90.0_r8
      k1%chi=h+90.0_r8
      n2%chi=2.0_r8*h-3.0_r8*s+p
      m2%chi=2.0_r8*h-2.0_r8*s
      s2%chi=0.0_r8
      k2%chi=2.0_r8*h
!
!  Compute tidal constituent frequency (1/s).
!
      q1%omega=0.6495854e-4_r8
      o1%omega=0.6759774e-4_r8
      k1%omega=0.7292117e-4_r8
      n2%omega=1.378797e-4_r8
      m2%omega=1.405189e-4_r8
      s2%omega=1.454441e-4_r8
      k2%omega=1.458423e-4_r8
!
!  Compute tidal constituent amplitude (m).
!
      q1%amp= 1.9273e-2_r8
      o1%amp=10.0661e-2_r8
      k1%amp=14.1565e-2_r8
      n2%amp= 4.6397e-2_r8
      m2%amp=24.2334e-2_r8
      s2%amp=11.2743e-2_r8
      k2%amp= 3.0684e-2_r8
!
!  Compute tidal constituent Love number factors (1+k2-h2).  The Love
!  number is defined as the ratio of the body tide to the height of
!  the static equilibrium tide.
!
      q1%love=0.695_r8
      o1%love=0.695_r8
      k1%love=0.736_r8
      n2%love=0.693_r8
      m2%love=0.693_r8
      s2%love=0.693_r8
      k2%love=0.693_r8
!
!  Compute product of amp*f*love.
!
      q1%Afl=q1%amp*q1%f*q1%love
      o1%Afl=o1%amp*o1%f*o1%love
      k1%Afl=k1%amp*k1%f*k1%love
      n2%Afl=n2%amp*n2%f*n2%love
      m2%Afl=m2%amp*m2%f*m2%love
      s2%Afl=s2%amp*s2%f*s2%love
      k2%Afl=k2%amp*k2%f*k2%love
!
      RETURN
      END SUBROUTINE harmonic_constituents
#endif
      END MODULE equilibrium_tide_mod