rdb_ocean_tide_astro.F90 Source File

!! Equilibrium-tide astronomy generator (mean longitudes, equilibrium
!! arguments, nodal corrections, constituent catalog).
module rdb_ocean_tide_astro
   !! Clean-room astronomy for the equilibrium (astronomical) body-force
   !! tide.  Computes the four mean longitudes (moon `s`, sun `h`, lunar
   !! perigee `p`, ascending node `N`) from Schureman's polynomials, the
   !! per-constituent equilibrium argument `V_c` (with the load-bearing
   !! ±pi/2 diurnal signs), and the slowly-varying nodal amplitude/phase
   !! corrections `f_c(N)`, `u_c(N)`.  Host-side scalar generator — called
   !! once at init (to bake the ref-date offset into `phase0`) and once per
   !! outer step (to advance the running argument); NOT a device kernel.
   !!
   !! References (the recipe is built from these, not from any model source):
   !!   * Doodson (1921), Proc. R. Soc. Lond. A 100.
   !!   * Schureman (1958), "Manual of harmonic analysis and prediction of
   !!     tides", US C&GS Spec. Pub. 98 (mean-longitude polynomials; the
   !!     36525-day Julian century).
   !!   * Cartwright & Tayler (1971) / Cartwright & Edden (1973) (amplitudes).
   !!   * Kowalik & Luick (2019), "Modern Theory and Practice of Tide
   !!     Analysis and Prediction" (Tables I.4 argument, I.6 nodal).
   use, intrinsic :: iso_fortran_env, only: int64
   use rdb_constants, only: wp
   implicit none
   private

   public :: TIDE_NAME, TIDE_SPECIES, TIDE_AMP, TIDE_LOVE, TIDE_OMEGA
   public :: TIDES_CATALOG_SIZE
   public :: gregorian_day_number, days_since_1900
   public :: mean_longitudes, equilibrium_arguments, nodal_fu
   public :: tide_name_index, parse_date_string

   real(wp), parameter, public :: TIDE_PI = 4.0_wp*atan(1.0_wp)
      !! Pi to working precision.
   real(wp), parameter, public :: TIDE_DEG2RAD = TIDE_PI/180.0_wp
      !! Degrees -> radians.

   integer, parameter :: TIDES_CATALOG_SIZE = 10
      !! Full catalog: 8 default (M2 S2 N2 K2 K1 O1 P1 Q1) + MF MM branch.

   character(len=2), parameter :: TIDE_NAME(TIDES_CATALOG_SIZE) = &
                                  ["M2", "S2", "N2", "K2", "K1", "O1", "P1", "Q1", "MF", "MM"]
      !! Constituent names (catalog order).
   integer, parameter :: TIDE_SPECIES(TIDES_CATALOG_SIZE) = &
                         [2, 2, 2, 2, 1, 1, 1, 1, 3, 3]
      !! Species / structure-slice index: 1 diurnal, 2 semidiurnal,
      !! 3 long-period.  Directly indexes the (nx,ny,3) struct slices.
   real(wp), parameter :: TIDE_AMP(TIDES_CATALOG_SIZE) = &
                          [0.242334_wp, 0.112743_wp, 0.046397_wp, 0.030684_wp, &
                           0.141565_wp, 0.100661_wp, 0.046848_wp, 0.019273_wp, &
                           0.042041_wp, 0.022191_wp]
      !! Equilibrium amplitudes A (m).
   real(wp), parameter :: TIDE_LOVE(TIDES_CATALOG_SIZE) = &
                          [0.693_wp, 0.693_wp, 0.693_wp, 0.693_wp, &
                           0.736_wp, 0.695_wp, 0.706_wp, 0.695_wp, 0.693_wp, 0.693_wp]
      !! Effective Love-number factors (1 + k - h).
   real(wp), parameter :: TIDE_OMEGA(TIDES_CATALOG_SIZE) = &
                          [1.4051890e-4_wp, 1.4544410e-4_wp, 1.3787970e-4_wp, 1.4584234e-4_wp, &
                           0.7292117e-4_wp, 0.6759774e-4_wp, 0.7252295e-4_wp, 0.6495854e-4_wp, &
                           0.053234e-4_wp, 0.026392e-4_wp]
      !! Angular frequencies omega_c (rad/s).

contains

   pure function gregorian_day_number(y, m, d) result(jdn)
      !! Proleptic-Gregorian Julian Day Number (integer, at 00:00 UT).
      !! Fliegel & Van Flandern algorithm.
      integer, intent(in) :: y, m, d
      integer(int64) :: jdn
      integer(int64) :: a, yy, mm
      a = int((14 - m)/12, int64)
      yy = int(y, int64) + 4800_int64 - a
      mm = int(m, int64) + 12_int64*a - 3_int64
      jdn = int(d, int64) + (153_int64*mm + 2_int64)/5_int64 + 365_int64*yy &
            + yy/4_int64 - yy/100_int64 + yy/400_int64 - 32045_int64
   end function gregorian_day_number

   pure function days_since_1900(y, m, d) result(dnum)
      !! Days since the astronomical origin 1900-01-01 00:00 UT.
      integer, intent(in) :: y, m, d
      real(wp) :: dnum
      dnum = real(gregorian_day_number(y, m, d) &
                  - gregorian_day_number(1900, 1, 1), wp)
   end function days_since_1900

   pure subroutine mean_longitudes(dnum, s_deg, h_deg, p_deg, n_deg)
      !! Mean longitudes at day number `dnum` (Schureman polynomials).
      !! Returned in DEGREES, folded to [0,360).  T = dnum/36525 (Julian
      !! centuries, Schureman's 36525-day century).  Multiply by
      !! TIDE_DEG2RAD to get radians for the equilibrium arguments.
      real(wp), intent(in) :: dnum
      real(wp), intent(out) :: s_deg, h_deg, p_deg, n_deg
      real(wp) :: t
      t = dnum/36525.0_wp
      s_deg = wrap360(277.0248_wp + 481267.8906_wp*t + 0.0011_wp*t**2)
      h_deg = wrap360(280.1895_wp + 36000.7689_wp*t + 3.0310e-4_wp*t**2)
      p_deg = wrap360(334.3853_wp + 4069.0340_wp*t - 0.0103_wp*t**2)
      n_deg = wrap360(259.1568_wp - 1934.142_wp*t + 0.0021_wp*t**2)
   end subroutine mean_longitudes

   pure function wrap360(x) result(y)
      !! Fold an angle in degrees to the canonical residue [0,360).
      !! (Fortran `mod` returns the sign of the argument, so a negative
      !! polynomial value needs the +360 canonicalisation.)
      real(wp), intent(in) :: x
      real(wp) :: y
      y = mod(x, 360.0_wp)
      if (y < 0.0_wp) y = y + 360.0_wp
   end function wrap360

   pure subroutine equilibrium_arguments(dnum, v_arg)
      !! Equilibrium argument `V_c` (radians) for every catalog
      !! constituent at day number `dnum`.  s,h,p,N are taken in radians
      !! (deg-mod-360 -> rad); the result is left un-modded (cos/sin are
      !! periodic).  The ±pi/2 diurnal signs are load-bearing.
      real(wp), intent(in) :: dnum
      real(wp), intent(out) :: v_arg(TIDES_CATALOG_SIZE)
      real(wp) :: s_deg, h_deg, p_deg, n_deg
      real(wp) :: s, h, p
      real(wp) :: half_pi
      call mean_longitudes(dnum, s_deg, h_deg, p_deg, n_deg)
      s = s_deg*TIDE_DEG2RAD
      h = h_deg*TIDE_DEG2RAD
      p = p_deg*TIDE_DEG2RAD
      half_pi = 0.5_wp*TIDE_PI
      v_arg(1) = 2.0_wp*(h - s)                    ! M2
      v_arg(2) = 0.0_wp                            ! S2
      v_arg(3) = -3.0_wp*s + 2.0_wp*h + p          ! N2
      v_arg(4) = 2.0_wp*h                          ! K2
      v_arg(5) = h + half_pi                       ! K1  (+pi/2)
      v_arg(6) = -2.0_wp*s + h - half_pi           ! O1  (-pi/2)
      v_arg(7) = -h - half_pi                      ! P1  (-pi/2)
      v_arg(8) = -3.0_wp*s + h + p - half_pi       ! Q1  (-pi/2)
      v_arg(9) = 2.0_wp*s                          ! MF
      v_arg(10) = s - p                            ! MM
   end subroutine equilibrium_arguments

   pure subroutine nodal_fu(dnum, add_nodal, f_nodal, u_nodal)
      !! Nodal amplitude factor `f_c(N)` (nondim) and phase `u_c(N)`
      !! (radians), fixed at the nodal reference date's `N`.
      !! `add_nodal = .false.` returns f=1, u=0 for every constituent.
      real(wp), intent(in) :: dnum
      logical, intent(in) :: add_nodal
      real(wp), intent(out) :: f_nodal(TIDES_CATALOG_SIZE)
      real(wp), intent(out) :: u_nodal(TIDES_CATALOG_SIZE)
      real(wp) :: s_deg, h_deg, p_deg, n_deg
      real(wp) :: cos_n, sin_n
      integer :: c
      if (.not. add_nodal) then
         do c = 1, TIDES_CATALOG_SIZE
            f_nodal(c) = 1.0_wp
            u_nodal(c) = 0.0_wp
         end do
         return
      end if
      call mean_longitudes(dnum, s_deg, h_deg, p_deg, n_deg)
      cos_n = cos(n_deg*TIDE_DEG2RAD)
      sin_n = sin(n_deg*TIDE_DEG2RAD)
      ! M2, N2
      f_nodal(1) = 1.0_wp - 0.037_wp*cos_n
      u_nodal(1) = (-2.1_wp)*TIDE_DEG2RAD*sin_n
      f_nodal(3) = 1.0_wp - 0.037_wp*cos_n
      u_nodal(3) = (-2.1_wp)*TIDE_DEG2RAD*sin_n
      ! K2
      f_nodal(4) = 1.024_wp + 0.286_wp*cos_n
      u_nodal(4) = (-17.7_wp)*TIDE_DEG2RAD*sin_n
      ! K1
      f_nodal(5) = 1.006_wp + 0.115_wp*cos_n
      u_nodal(5) = (-8.9_wp)*TIDE_DEG2RAD*sin_n
      ! O1, Q1  (NOTE +sign on u)
      f_nodal(6) = 1.009_wp + 0.187_wp*cos_n
      u_nodal(6) = (10.8_wp)*TIDE_DEG2RAD*sin_n
      f_nodal(8) = 1.009_wp + 0.187_wp*cos_n
      u_nodal(8) = (10.8_wp)*TIDE_DEG2RAD*sin_n
      ! S2, P1
      f_nodal(2) = 1.0_wp
      u_nodal(2) = 0.0_wp
      f_nodal(7) = 1.0_wp
      u_nodal(7) = 0.0_wp
      ! MF
      f_nodal(9) = 1.043_wp + 0.414_wp*cos_n
      u_nodal(9) = (-23.7_wp)*TIDE_DEG2RAD*sin_n
      ! MM
      f_nodal(10) = 1.0_wp - 0.130_wp*cos_n
      u_nodal(10) = 0.0_wp
   end subroutine nodal_fu

   pure function tide_name_index(name) result(idx)
      !! Catalog index (1..TIDES_CATALOG_SIZE) of a constituent name,
      !! matched case-insensitively; -1 if unknown.
      character(len=*), intent(in) :: name
      integer :: idx
      integer :: c
      character(len=2) :: up
      up = upcase2(adjustl(name))
      idx = -1
      do c = 1, TIDES_CATALOG_SIZE
         if (up == TIDE_NAME(c)) then
            idx = c
            return
         end if
      end do
   end function tide_name_index

   pure function upcase2(s) result(u)
      !! Uppercase the first two characters of a token.
      character(len=*), intent(in) :: s
      character(len=2) :: u
      integer :: k, ic
      u = "  "
      do k = 1, min(2, len_trim(s))
         ic = iachar(s(k:k))
         if (ic >= iachar("a") .and. ic <= iachar("z")) ic = ic - 32
         u(k:k) = achar(ic)
      end do
   end function upcase2

   pure subroutine parse_date_string(str, y, m, d, ok)
      !! Parse a "YYYY-MM-DD" calendar date.  `ok = .false.` on any
      !! malformed field (caller decides fail-loud policy).
      character(len=*), intent(in) :: str
      integer, intent(out) :: y, m, d
      logical, intent(out) :: ok
      character(len=:), allocatable :: t
      character(len=128) :: emsg
      integer :: p1, p2, ios
      y = 0
      m = 0
      d = 0
      ok = .false.
      t = trim(adjustl(str))
      p1 = index(t, "-")
      if (p1 <= 1) return
      p2 = index(t(p1 + 1:), "-")
      if (p2 <= 1) return
      p2 = p1 + p2
      read (t(1:p1 - 1), *, iostat=ios, iomsg=emsg) y
      if (ios /= 0) return
      read (t(p1 + 1:p2 - 1), *, iostat=ios, iomsg=emsg) m
      if (ios /= 0) return
      read (t(p2 + 1:), *, iostat=ios, iomsg=emsg) d
      if (ios /= 0) return
      if (m < 1 .or. m > 12 .or. d < 1 .or. d > 31) return
      ok = .true.
   end subroutine parse_date_string

end module rdb_ocean_tide_astro