Commit d53cb277 authored by thomas.forbriger's avatar thomas.forbriger

libseife [TASK]: add M3 frequency to C-version of tide removal

parent 06c5a73b
......@@ -30,6 +30,7 @@
* REVISIONS and CHANGES
* - 28/06/2005 V1.0 Thomas Forbriger
* - 22/06/2020 V1.1 handle case of nstep=1 explicitly
* V1.2 make array size flexible and add M3 frequency
*
* ============================================================================
*/
......@@ -51,7 +52,10 @@ void seife_tid(double* x, int n, double dt, int nstep)
{
/* Initialized data */
static double omega[6] = { 1.93227,.92954,2.,1.00274,1.89567,.89293 };
#define C_MFREQ 7
#define C_MDIM (C_MFREQ*2+1)
static double omega[C_MFREQ] = { 2.898450424,
1.93227,.92954,2.,1.00274,1.89567,.89293 };
static double zero = 0.;
static double one = 1.;
static double two = 2.;
......@@ -63,18 +67,19 @@ void seife_tid(double* x, int n, double dt, int nstep)
// double atan(double), cos(double), sin(double);
/* Local variables */
static double tdif, omeg[6];
static double tdif, omeg[C_MFREQ];
static int ndim;
static double step, tint, omeg0, xgez1, xgez2, xgez3, a[169] /*
was [13][13] */, c__[13], d__[13], e[13], f[13];
static double step, tint, omeg0, xgez1, xgez2, xgez3,
a[C_MDIM*C_MDIM], c__[C_MDIM], d__[C_MDIM],
e[C_MDIM], f[C_MDIM];
static int i__, j, k;
static double t;
static int nfreq, k2, k3;
static double tstep2;
static int jj;
static double rs[13], sx, dth;
static double rs[C_MDIM], sx, dth;
static int nco;
static double cor[6], dur;
static double cor[C_MFREQ], dur;
// extern /* Subroutine */ int seife_gauss__(doublereal *, integer *,
// integer *, doublereal *, doublereal *);
......@@ -84,7 +89,7 @@ void seife_tid(double* x, int n, double dt, int nstep)
--x;
/* Function Body */
ndim = 13;
ndim = C_MDIM;
dth = dt / two;
if (nstep == 0) {
......@@ -101,18 +106,18 @@ void seife_tid(double* x, int n, double dt, int nstep)
}
/* determine the number of frequencies required for a good fit */
dur = n * dt / 3600.;
nfreq = 6;
nfreq = C_MFREQ;
if (dur < 35.) {
nfreq = 5;
nfreq = 6;
}
if (dur < 18.) {
nfreq = 4;
nfreq = 5;
}
if (dur < 14.) {
nfreq = 3;
nfreq = 4;
}
if (dur < 5.) {
nfreq = 2;
nfreq = 3;
}
nco = (nfreq << 1) + 1;
omeg0 = atan(one) * 8. / 86400.;
......@@ -128,7 +133,7 @@ void seife_tid(double* x, int n, double dt, int nstep)
i__2 = nco;
for (k = 1; k <= i__2; ++k) {
/* L1: */
a[i__ + k * 13 - 14] = 0.;
a[i__ + k * C_MDIM - (C_MDIM+1)] = 0.;
}
}
if (nstep == 1) {
......@@ -160,7 +165,7 @@ void seife_tid(double* x, int n, double dt, int nstep)
rs[i__ - 1] += sx * c__[i__ - 1];
i__3 = nco;
for (k = 1; k <= i__3; ++k) {
a[i__ + k * 13 - 14] += c__[i__ - 1] * c__[k - 1];
a[i__ + k * C_MDIM - (C_MDIM+1)] += c__[i__ - 1] * c__[k - 1];
}
}
}
......@@ -218,7 +223,7 @@ void seife_tid(double* x, int n, double dt, int nstep)
i__4 = nco;
for (k = 1; k <= i__4; ++k) {
/* L2: */
a[i__ + k * 13 - 14] += c__[i__ - 1] * c__[k - 1];
a[i__ + k * C_MDIM - (C_MDIM+1)] += c__[i__ - 1] * c__[k - 1];
}
}
}
......
Markdown is supported
0% or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment