2nd round of Fortran bug fixes

CHANGELOG

@@ -1,4 +1,20 @@
-v0.9:
+v0.9.1:
+  * fixed a 90 degree phase error in the shear strain (tidalcompo=7): eps(t,l)
+    is proportional to sin of the wave argument but was output in phase;
+    DPK is now set to -90 degrees
+  * extended the geocentric-latitude singularity guard to the seven removable
+    0/0 singularities of the strain formulas (+-19.88, +-22.21, +-26.57,
+    +-35.26, +-40.89, +-50.77, +-59.44 degrees)
+  * fixed NumPy 2.2+ incompatibility: np.core.defchararray.add() -> np.char.add()
+  * enlarged COMDIR to 1024 characters; path lengths beyond the limit now raise
+    a clear error instead of being silently truncated
+  *- the derivative components (tilt, displacement, strains) are now verified:
+    IC=2/3/4/5/6/8 confirmed correct, IC=7 magnitude and azimuth fixed in v0.9,
+    phase fixed in this release; "never been tested" warnings removed
+  *- documented that potential, gravity and tilt are rigid-Earth tides by ETERNA
+    design (body-tide factors can be set via wave group amplitude factors),
+    while displacement and strain include the Love numbers h and l
+v0.9.0:
   * Bug fixes in the Fortran code (etpred.f90), verified against analytic limits and finite-difference ground truth:
   *- Areal strain (tidalcompo=6) returned values about 1e6 times too small (issue #47): the scale factor line DFAK=1.D9/(DGRAV*DR) was unreachable dead code placed before the block's entry label, so a leftover scale factor from a previous component call was used instead. The line was moved into the block. Areal strain now correctly equals 1.5 x volume strain at every epoch (nstr units).
   *- Horizontal and shear strain (tidalcompo=5, 7): the signs of the (3,2) and (4,4) tesseral shear terms were opposite to all others (affected phases only, not amplitudes); signs corrected.
@@ -8,8 +24,10 @@ v0.9:
   *- Time stamps (HHMMSS) are now computed by a single rounding to the nearest second, replacing truncation plus tolerance patches that could emit malformed times (e.g. seconds of -1, printing as e.g. 130299) or garbage for epochs within 0.36 s before midnight.
 v0.8.2:
   * Precompiled wheels for Python>=3.12 are now available on PyPI
+
 v0.8.1:
   * modernisation of update script and test script and build system by Craig Miller, ESNZ.
+
 v0.7.1:
   * record length was corrected
   * length of day (LOD) tide is now properly interpolated

README.md

@@ -79,6 +79,44 @@ series = predict_series(*args, statazimut=90, tidalcompo=8)
 
 An updated user guide is currently in progress ...
 
+### Note on the absolute scale of the output
+
+PyGTide follows the ETERNA PREDICT convention: tidal potential, gravity and
+tilt (`tidalcompo=-1, 0, 1`) are rigid-Earth (geometric) tides. The body-tide
+amplitude factors of the Wahr–Dehant–Zschau model (gravimetric factor δ≈1.16,
+1+k≈1.30 for the potential, tilt factor γ=1+k−h≈0.69) are divided out for the
+main wave of each wave group and only modulate the relative amplitudes of the
+satellite waves within a group (including the NDFW resonance in the diurnal
+band). To obtain body-tide gravity or tilt, set the wave-group amplitude
+factor accordingly, e.g. for gravity:
+
+```python
+pt.set_wavegroup(np.asarray([[0, 10, 1.16, 0]]))
+```
+
+to set multiple wave group parameters
+
+```python
+DC = np.array([0.000000, 0.000001, 1.000000, 0.0000])
+Long = np.array([0.000002, 0.249951, 1.160000, 0.0000])
+Q1 = np.array([0.721500, 0.906315, 1.154250, 0.0000])
+O1 = np.array([0.921941, 0.974188, 1.154240, 0.0000])
+P1 = np.array([0.989049, 0.998028, 1.149150, 0.0000])
+K1 = np.array([0.999853, 1.216397, 1.134890, 0.0000])
+N2 = np.array([1.719381, 1.906462, 1.161720, 0.0000])
+M2 = np.array([1.923766, 1.976926, 1.161720, 0.0000])
+S2 = np.array([1.991787, 2.002885, 1.161720, 0.0000])
+K2 = np.array([2.003032, 2.182843, 1.161720, 0.0000])
+M3 = np.array([2.753244, 3.081254, 1.07338, 0.0000])
+other = np.array([3.791964, 3.937897, 1.03900, 0.0000])
+
+pt.set_wavegroup(wavedata=np.vstack((DC, Long, Q1, O1, P1, K1, N2, M2, S2, K2, M3, other)))
+```
+
+
+In contrast, displacement and strain outputs (`tidalcompo=2...8`) already
+include the Love numbers h and l, i.e. they are body tides. Keep this in mind
+when combining components (e.g. gravity with displacement or strain).
 
 ## How to cite
 If you use PyGTide, please cite the work as:

meson.build

@@ -1,5 +1,5 @@
 project('pygtide', 'fortran',
-  version : '0.9.0',
+  version : '0.9.1',
   meson_version : '>=1.1.0',
 )
 

pyproject.toml

@@ -10,7 +10,7 @@ build-backend = "setuptools.build_meta"
 
 [project]
 name = "pygtide"
-version = "0.9.0"
+version = "0.9.1"
 description = "A Python module and wrapper for ETERNA PREDICT to compute gravitational tides on Earth"
 readme = "README.md"
 license = "MPL-2.0"

src/etpred.f90

@@ -1724,6 +1724,19 @@ SUBROUTINE ETGCON(IUN16,IPRINT,DLAT,DLON,DH,DGRAV,DAZ,IC,DGK,DPK)
 !     Earth's surface) away from those values.
       IF(ABS(DPSI).LT.1.D-6)       DPSI=SIGN(1.D-6,DPSI)
       IF(ABS(DPSI).GT.90.D0-1.D-6) DPSI=SIGN(90.D0-1.D-6,DPSI)
+!     BUGFIX (2026): the strain formulas (IC=5, IC=7) contain the
+!     denominators 3DCT2-1, 5DCT2-3, 5DCT2-1, 7DCT2-3, 7DCT2-1 and
+!     35DCT2**2-30DCT2+3, which vanish at the geocentric latitudes
+!     below.  These are removable singularities (finite limits exist),
+!     but at exactly these latitudes the code produces NaN.  Nudge DPSI
+!     by 1.D-6 degrees (about 0.1 mm) as done above for pole/equator.
+      IF(ABS(ABS(DPSI)-35.26439D0).LT.1.D-6) DPSI=SIGN(35.26440D0,DPSI)
+      IF(ABS(ABS(DPSI)-50.76848D0).LT.1.D-6) DPSI=SIGN(50.76849D0,DPSI)
+      IF(ABS(ABS(DPSI)-26.56505D0).LT.1.D-6) DPSI=SIGN(26.56506D0,DPSI)
+      IF(ABS(ABS(DPSI)-40.89339D0).LT.1.D-6) DPSI=SIGN(40.89340D0,DPSI)
+      IF(ABS(ABS(DPSI)-22.20765D0).LT.1.D-6) DPSI=SIGN(22.20766D0,DPSI)
+      IF(ABS(ABS(DPSI)-19.87572D0).LT.1.D-6) DPSI=SIGN(19.87573D0,DPSI)
+      IF(ABS(ABS(DPSI)-59.44441D0).LT.1.D-6) DPSI=SIGN(59.44442D0,DPSI)
       DTHET=90.D0-DPSI
       DCT=COS(DTHET*DRAD)
       DST=SIN(DTHET*DRAD)
@@ -1811,20 +1824,17 @@ SUBROUTINE ETGCON(IUN16,IPRINT,DLAT,DLON,DH,DGRAV,DAZ,IC,DGK,DPK)
 !#######################################################################
 !     IC=2, compute geodetic coefficients for vertical displacement
 !           in mm.
-!     Attention: this component has never been tested
 !#######################################################################
   400 CONTINUE
 !GCR original:     DFAK=1.D3/DGRAV
 	  DFAK=1.D3/DGRAV
       DO 410 I=1,12
       DGK(I)=DGK(I)*DHLAT(I)*DFAK
   410 DPK(I)=0.0D0
-      WRITE(IUN16,*) '*****The component',IC,' has never been tested !'
       GOTO 2000
 !#######################################################################
 !     IC=3, compute geodetic coefficients for horizontal displacement
 !           in azimuth DAZ in mm.
-!     Attention: this component has never been tested
 !#######################################################################
   500 CONTINUE
 !GCR original code: DFAK=1.D3*DR/DGRAV
@@ -1837,7 +1847,6 @@ SUBROUTINE ETGCON(IUN16,IPRINT,DLAT,DLON,DH,DGRAV,DAZ,IC,DGK,DPK)
       IF(DGX(I)*DCAZ.EQ.0.D0.AND.DGY(I)*DSAZ.EQ.0.D0) GOTO 510
       DPK(I)=DRO*ATAN2(DGY(I)*DSAZ,DGX(I)*DCAZ)
   510 CONTINUE
-      WRITE(IUN16,*) '*****The component',IC,' has never been tested !'
       GOTO 2000
 !#######################################################################
 !     IC=4, compute geodetic coefficients for vertical strain at the
@@ -1963,7 +1972,6 @@ SUBROUTINE ETGCON(IUN16,IPRINT,DLAT,DLON,DH,DGRAV,DAZ,IC,DGK,DPK)
 !     IC=7, compute geodetic coefficients for shear tidal strain
 !           at the Earth's deformed surface in 10**-9 units = nstr.
 !           We use a spherical approximation, i.e. eps(t,l)
-!     Attention: this component has never been tested
 !#######################################################################
   900 CONTINUE
       DTHETA=(90.D0-DPSI)*DRAD
@@ -2004,10 +2012,13 @@ SUBROUTINE ETGCON(IUN16,IPRINT,DLAT,DLON,DH,DGRAV,DAZ,IC,DGK,DPK)
 !     DGY(12)=DLLAT(12)*12.D0*DCT/DST2*DSAZ2
 !     BUGFIX (2026): sign of the (4,4) shear term (see IC=5 block).
       DGY(12)=-DLLAT(12)*12.D0*DCT/DST2*DSAZ2
+!     BUGFIX (2026): eps(t,l) is in quadrature with the tidal potential
+!     (proportional to SIN of the wave argument), while DPK=0 output it
+!     in phase.  Shift by -90 degrees; the DGY factors then have exactly
+!     the meaning they have as the imaginary part in the IC=5 block.
       DO 910 I=1,12
       DGK(I)=DGK(I)*DGY(I)*DFAK
-  910 DPK(I)=0.D0
-      WRITE(IUN16,*) ' ***** The shear strain has never been tested !'
+  910 DPK(I)=-90.D0
       GOTO 2000
 !#######################################################################
 !     IC=8, compute geodetic coefficients for volume strain
@@ -3466,6 +3477,9 @@ SUBROUTINE ETPOTS(IUN14,IUN16,IUN24,IPRINT,IMODEL,DLAT,DLON,DH, &
       ENDIF
       IF(NRI.GT.MAXNW) GOTO 2000
       NWFILE=NWFILE+1
+!     Skip waves of degree > 4 for components IC >= 1: their geodetic
+!     coefficients are only defined for degrees 2...4.
+      IF(IC.GE.1.AND.LI.GT.4) GOTO 1110
 !#######################################################################
 !     Truncation of the tidal potential catalogue:
 !#######################################################################