src/equilibrium_tide.h

Equilibrium tide

Equilibrium tide

This is directly adapted from equilibrium_tide.F in the ROMS source code.

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

[roms2021]	Arango et al. Wiki ROMS. September 2021. [ http ]
[arbic2018]	Brian K Arbic, Matthew H Alford, Joseph K Ansong, Maarten C Buijsman, Robert B Ciotti, J Thomas Farrar, Robert W Hallberg, Christopher E Henze, Christopher N Hill, Conrad A Luecke, et al. A primer on global internal tide and internal gravity wave continuum modeling in HYCOM and MITgcm. New frontiers in operational oceanography, 2018. [ DOI ]
[egbert2017]	Gary D Egbert and Richard D Ray. Tidal prediction. Journal of Marine Research, 75:189–237, 2017.
[arbic2004]	Brian K Arbic, Stephen T Garner, Robert W Hallberg, and Harper L Simmons. 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):3069–3101, 2004.
[doodson1941]	Arthur Thomas Doodson and Harold Dreyer Warburg. Admiralty manual of tides. His Majesty’s Stationery Office, 1941.

Define equilibrium tide constituents structure

typedef struct {
  double Afl;   // product: amp*f*love
  double amp;   // amplitude (m)
  double chi;   // phase at Greenwich meridian (degrees)
  double f;     // f nodal factor (nondimensional)
  double love;  // tidal Love number factor
  double nu;    // nu nodal factor (degrees)
  double omega; // frequency (1/s)
} Etide_t;

struct {
  Etide_t Q1;   // Q1 component, diurnal
  Etide_t O1;   // O1 component, diurnal
  Etide_t K1;   // K1 component, diurnal
  Etide_t N2;   // N2 component, semi-diurnal
  Etide_t M2;   // M2 component, semi-diurnal
  Etide_t S2;   // S2 component, semi-diurnal
  Etide_t K2;   // N2 component, semi-diurnal  
} Etide = {0};

const double deg2rad = pi/180.;

#include <time.h>

Computation of the tidal constituents

The tidal constituents are computed relative to the optional starting date given as argument.

If Lnodal is true then lunar nodal correction is applied.

void equilibrium_tide_constituents (const char * date = "2000-01-01 12:00:00",
				    bool Lnodal = false)
{

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

  struct tm astro_tm;
  assert (strptime ("2000-01-01 12:00:00", "%Y-%m-%d %H:%M:%S", &astro_tm));
  time_t astro_t = mktime (&astro_tm);
  assert (astro_t != (time_t) -1);

Terrestial time (in centuries) since tide reference date number. Recall that the length of a year is 365.2425 days for the now called Gregorian Calendar (corrected after 15 October 1582).

We also assume that a day is 86400 seconds (need to check this). Note that the original code in ROMS seems to round T to the nearest day??, which we do not do here.

  struct tm date_tm;
  if (!strptime (date, "%Y-%m-%d %H:%M:%S", &date_tm)) {
    fprintf (stderr,
	     "%s:%d: error: date string '%s' is not formatted as '%%Y-%%m-%%d %%H:%%M:%%S'\n",
	     __FILE__, LINENO, date);
    exit (1);
  }
  time_t date_t = mktime (&date_tm);
  assert (date_t != (time_t) -1);
  double T  = difftime (date_t, astro_t)/(86400.*36524.25);

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

Mean longitude of lunar node (Period = 18.6 years).

    double N = -234.955 - 1934.1363*T; // degrees
    N *= deg2rad; // radians

    Etide.O1.f = 1.009 + 0.187*cos(N) - 0.015*cos(2.0*N);
    Etide.K1.f = 1.006 + 0.115*cos(N) - 0.009*cos(2.0*N);
    Etide.M2.f = 1.0 - 0.037*cos(N);
    Etide.S2.f = 1.0;
    Etide.K2.f = 1.024 + 0.286*cos(N) + 0.008*cos(2.0*N);
    
    Etide.O1.nu = 10.8*sin(N) - 1.3*sin(2.0*N);
    Etide.K1.nu = -8.9*sin(N) + 0.7*sin(2.0*N);
    Etide.M2.nu = -2.1*sin(N);
    Etide.S2.nu = 0.0;
    Etide.K2.nu = -17.7*sin(N) + 0.7*sin(2.0*N);    
  }
  else {
    Etide.O1.f = 1.0;
    Etide.K1.f = 1.0;
    Etide.M2.f = 1.0;
    Etide.S2.f = 1.0;
    Etide.K2.f = 1.0;
    
    Etide.O1.nu = 0.0;
    Etide.K1.nu = 0.0;
    Etide.M2.nu = 0.0;
    Etide.S2.nu = 0.0;
    Etide.K2.nu = 0.0;
  }
  Etide.Q1.f = Etide.O1.f;
  Etide.N2.f = Etide.M2.f;
  Etide.Q1.nu = Etide.O1.nu;
  Etide.N2.nu = Etide.M2.nu;

Mean longitude of the moon (Period = tropical month).

  double s = 218.316 + 481267.8812*T;

Mean longitude of the sun (Period = tropical year).

  double h = 280.466 + 36000.7698*T;

Mean longitude of lunar perigee (Period = 8.85 years).

  double p = 83.353 + 4069.0137*T;

Compute tidal constituent phase “chi” (degrees) at Greenwich meridian.

  Etide.Q1.chi = h - 3.0*s + p - 90.0;
  Etide.O1.chi = h - 2.0*s - 90.0;
  Etide.K1.chi = h + 90.0;
  Etide.N2.chi = 2.0*h - 3.0*s + p;
  Etide.M2.chi = 2.0*h - 2.0*s;
  Etide.S2.chi = 0.0;
  Etide.K2.chi = 2.0*h;

Set tidal constituent frequency (1/s).

  Etide.Q1.omega = 0.6495854E-4;
  Etide.O1.omega = 0.6759774E-4;
  Etide.K1.omega = 0.7292117E-4;
  Etide.N2.omega = 1.378797E-4;
  Etide.M2.omega = 1.405189E-4;
  Etide.S2.omega = 1.454441E-4;
  Etide.K2.omega = 1.458423E-4;

Set tidal constituent amplitude (m).

  Etide.Q1.amp =  1.9273E-2;
  Etide.O1.amp = 10.0661E-2;
  Etide.K1.amp = 14.1565E-2;
  Etide.N2.amp =  4.6397E-2;
  Etide.M2.amp = 24.2334E-2;
  Etide.S2.amp = 11.2743E-2;
  Etide.K2.amp =  3.0684E-2;

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

  Etide.Q1.love = 0.695;
  Etide.O1.love = 0.695;
  Etide.K1.love = 0.736;
  Etide.N2.love = 0.693;
  Etide.M2.love = 0.693;
  Etide.S2.love = 0.693;
  Etide.K2.love = 0.693;

Compute product of ampflove.

  Etide.Q1.Afl = Etide.Q1.amp*Etide.Q1.f*Etide.Q1.love;
  Etide.O1.Afl = Etide.O1.amp*Etide.O1.f*Etide.O1.love;
  Etide.K1.Afl = Etide.K1.amp*Etide.K1.f*Etide.K1.love;
  Etide.N2.Afl = Etide.N2.amp*Etide.N2.f*Etide.N2.love;
  Etide.M2.Afl = Etide.M2.amp*Etide.M2.f*Etide.M2.love;
  Etide.S2.Afl = Etide.S2.amp*Etide.S2.f*Etide.S2.love;
  Etide.K2.Afl = Etide.K2.amp*Etide.K2.f*Etide.K2.love;
}

Computes surface elevation associated with the equilibrium tide

Assumes that the coordinate system is longitude in degrees (x) and latitude in degrees (y).

Note that the tidal constituents must first be defined by calling equilibrium_tide_constituents().

The time t given (in seconds) is relative to the date/time specified when calling equilibrium_tide_constituents().

void equilibrium_tide (scalar tide, double t)
{
  if (!Etide.Q1.omega) {
    fprintf (stderr,
	     "%s:%d: error: the tidal constituents must first be defined by calling "
	     "equilibrium_tide_constituents()\n", __FILE__, LINENO);
    exit (1);
  }

Compute astronomical equilibrium tide (m) with diurnal and semidiurnal constituents. The long-period constituents and high-order harmonics are neglected.

  foreach()
    tide[] =
    sin(2.*deg2rad*y)*(Etide.Q1.Afl*cos(Etide.Q1.omega*t +
					deg2rad*(x + Etide.Q1.chi + Etide.Q1.nu)) + 
		       Etide.O1.Afl*cos(Etide.O1.omega*t +
					deg2rad*(x + Etide.O1.chi + Etide.O1.nu)) + 
		       Etide.K1.Afl*cos(Etide.K1.omega*t +
					deg2rad*(x + Etide.K1.chi + Etide.K1.nu))) + 
    sq(cos(deg2rad*y))*(Etide.N2.Afl*cos(Etide.N2.omega*t +
					 deg2rad*(2.0*x + Etide.N2.chi + Etide.N2.nu)) +  
			Etide.M2.Afl*cos(Etide.M2.omega*t +
					 deg2rad*(2.0*x + Etide.M2.chi + Etide.M2.nu)) +
			Etide.S2.Afl*cos(Etide.S2.omega*t +
					 deg2rad*(2.0*x + Etide.S2.chi + Etide.S2.nu)) +
			Etide.K2.Afl*cos(Etide.K2.omega*t +
					 deg2rad*(2.0*x + Etide.K2.chi + Etide.K2.nu)));
}

Usage

Examples

Global Tides