WAVEWATCH III  beta 0.0.1
w3tidemd.F90
Go to the documentation of this file.
1 #include "w3macros.h"
2 !/ ------------------------------------------------------------------- /
3 MODULE w3tidemd
4  !/
5  !/ +-----------------------------------+
6  !/ | WAVEWATCH III |
7  !/ ! M. Foreman, IOS !
8  !/ | FORTRAN 90 |
9  !/ | Last update : 21-Apr-2020 |
10  !/ +-----------------------------------+
11  !/
12  !/ 01-Sep-2012 : Origination. ( version 4.07 )
13  !/ 04-Mar-2013 : Correction of FAST and new VFAST ( version 4.08 )
14  !/ 21-Apr-2020 : Correction of time and implicit none( version 7.13 )
15  !/
16  ! 1. Purpose :
17  !
18  ! Tidal analysis of time series for storage of tidal constituents
19  ! only. This module is built around the versatile tidal analysis
20  ! package : http://www.pac.dfo-mpo.gc.ca/science/oceans/tidal-marees/index-eng.htm
21  ! by Mike Foreman et al. (see publication in J. Ocean Atmos. Tech.: vol.
22  ! http://journals.ametsoc.org/doi/pdf/10.1175/2008JTECHO615.1 )
23  ! Adaptation to WAVEWATCH III was performed by F. Ardhuin
24  !
25  ! STILL TO BE DONE:
26  ! - adding a namelist in ww3_grid to allow adjustment of TIDE_DT
27  ! - check on constituents (M2, S2, N2 ...) when running ww3_shel,
28  ! in order to allow use of different sets of constituents
29  ! - add residual currents (or geostrophic ...) ...
30  ! - make this work with multigrids
31  !
32  ! 2. Variables and types :
33  !
34  ! Name Type Scope Description
35  ! ----------------------------------------------------------------
36  ! ----------------------------------------------------------------
37  !
38  ! 3. Subroutines and functions :
39  !
40  ! Name Type Scope Description
41  ! ----------------------------------------------------------------
42  ! TIDE_WRITE_RESULTS Subr. Public Writes tidal results
43  ! (with M. Forman's format)
44  ! TIDE_PREDICT Subr. Public Predicts tide from amp. & phases
45  ! TIDE_SETTINGS_FULL Subr. Public Choice of constituents
46  ! TIDE_SETTINGS_FAST Subr. Public Choice of constituents
47  ! TIDE_SETTINGS_VFAST Subr. Public Choice of constituents
48  ! TIDE_SET_INDICES Subr. Public
49  ! SETVUF_FAST Subr. Public Calculates the V,u,f values
50  ! TIDE_READ_SETTINGS Subr. Public Reads data from file (IOS format)
51  ! TIDE_READ_ANAPAR Subr. Public Reads data from file (IOS format)
52  ! TIDE_READ_TIMESERIES Subr. Public Reads data from file (IOS format)
53  ! ASTR Subr. Public Calculates the ephermides
54  ! JULDAYT Func. Public Julian day
55  ! CALDATT Subr. Public
56  ! dsvbksb Subr. Public
57  ! dsvdcmp Subr. Public
58  ! svd Subr. Public Matrix singular value decomposition
59  ! VUF_SET_PARAMETERS Subr. Public
60  ! VUF Subr. Public
61  ! OPNVUF Subr. Public
62  ! SETVUF Subr. Public
63  ! flex_tidana_webpage Subr. Public
64  ! TIDE_PREDICT Subr. Public Tide prediction and error estimate
65  ! TIDE_PREDICT_ONLY Subr. Public Tide prediction only
66  ! ----------------------------------------------------------------
67  !
68  ! 4. Subroutines and functions used :
69  !
70  ! Name Type Module Description
71  ! ----------------------------------------------------------------
72  ! ----------------------------------------------------------------
73  !
74  ! 5. Remarks :
75  !
76  ! 6. Switches :
77  !
78  ! 7. Source code :
79  !/
80  !/ ------------------------------------------------------------------- /
81  !/
82  ! PUBLIC
83  !/
84  !/ Private variables
85  !/
86  DOUBLE PRECISION, PARAMETER :: twpi=3.1415926535898*2.
87  DOUBLE PRECISION, PARAMETER :: fac=twpi/360.
88 
89  !
90  ! Array sizes
91  !
92  INTEGER, PARAMETER :: mc=70,nr=106000,nmaxp1=mc*2,nmaxpm=nr*2+nmaxp1
93  INTEGER, PARAMETER :: mc2=mc*2
94  CHARACTER*5, PARAMETER :: kblank=' '
95  !
96  INTEGER :: ntidal_con, ntotal_con, nkonco
97  CHARACTER*5, ALLOCATABLE :: tidecon_allnames(:) ! array of names of tidal constituents
98  CHARACTER*5, ALLOCATABLE :: konco_con(:)
99  INTEGER, ALLOCATABLE, PRIVATE :: ii(:),jj(:),kk(:),ll(:),mm(:),nn(:),nj(:)
100  REAL, ALLOCATABLE :: semi(:),coef_con(:)
101  REAL , ALLOCATABLE :: v_arg(:,:),f_arg(:,:),u_arg(:,:)
102  REAL :: ee(180),ph(180)
103  INTEGER :: ldel(180),mdel(180),ndel(180),ir(180)
104  ! these two index table are used in VUF ...
105  INTEGER , ALLOCATABLE :: tide_indexj(:),tide_indexjk(:)
106  !
107  ! Parameters for tidal analysis
108  !
109  INTEGER :: tide_mf, tide_nx, tide_ny
110  REAL, ALLOCATABLE :: tide_freqc(:) ! array of freq. of tidal constituents
111  CHARACTER(LEN=5), ALLOCATABLE:: tidecon_namei(:) ! array of names of tidal constituents
112  CHARACTER(LEN=5), ALLOCATABLE:: tidecon_name(:) ! array of names of tidal constituents
113  CHARACTER(LEN=5) :: tide_konan(10), tide_konin(10,10)
114  REAL :: tide_r(10,10), tide_zeta(10,10)
115  REAL :: tide_sigan(10),tide_sigin(10,10) ! these two are only read from files and written out
116  INTEGER :: tide_nin,tide_ninf(10)
117  REAL, ALLOCATABLE :: tidal_const(:,:,:,:,:) ! array of freq. of tidal constituents
118  !
119  ! Data to be analyzed
120  !
121  INTEGER(KIND=4) :: tide_nti
122  REAL, ALLOCATABLE :: tide_data(:,:)
123  INTEGER(KIND=4), ALLOCATABLE :: tide_days(:), tide_secs(:)
124  REAL(kind=8), allocatable :: tide_hours(:)
125  REAL, PARAMETER :: tide_dt = 1800. ! time step used for forcing
126  !
127  ! Analysis result
128  !
129  REAL :: tide_ampc(mc,2), tide_phg(mc,2), &
130  tide_sig1(mc,2), tide_sig2(mc,2), &
131  tide_sig3(mc,2), tide_ttest(mc,2)
132  REAL :: tide_ampci(10,10,2), tide_phgi(10,10,2)
133  INTEGER :: tide_index(mc),tide_index2(mc)
134 
135 
136  INTEGER :: ndset, tide_verbose = 0
137 
138  !PUBLIC :: TIDE_MF, TIDECON_NAME
139 
140  !/
141 CONTAINS
142 
143  !/ ------------------------------------------------------------------- /
144  SUBROUTINE tide_write_results(LP,filename,ndef,KD1, KD2, ITZ, xlat,xlon, &
145  RES, SSQ, RMSR0, SDEV0,SDEV,RMSR, RESMAX, IMAX, RMSRP)
146 
147  IMPLICIT NONE
148  !
149  CHARACTER*256, INTENT(IN) :: filename
150  INTEGER, INTENT(IN) :: LP, NDEF, IMAX(NDEF)
151  INTEGER(KIND=4),INTENT(IN) :: KD1,KD2
152  CHARACTER*4 , INTENT(IN) :: ITZ
153  REAL(KIND=8), intent(in) :: rmsr0(ndef), &
154  sdev0(ndef), sdev(ndef), rmsr(ndef), resmax(ndef), rmsrp(ndef)
155  REAL , INTENT(IN) :: RES(NDEF), SSQ(NDEF), XLAT,XLON
156  !
157  INTEGER :: IDEF, I, K, K2, L, I1, INFTOT
158  INTEGER :: ID1,IM1,IY1,ID2,IM2,IY2
159 
160  open(unit=lp,file=filename,status='unknown',form='formatted')
161 
162  CALL caldatt(kd1,id1,im1,iy1)
163  CALL caldatt(kd2,id2,im2,iy2)
164 
165  WRITE(lp,15) id1,im1,iy1,id2,im2,iy2,itz
166 15 FORMAT(/'THE ANALYSIS PERIOD IS FROM',i3,'/',i2,'/',i4, &
167  ' TO ',i2,'/',i2,'/',i4,' IN THE TIME ZONE ',a4)
168  WRITE(lp,*)'USING SVD TO SOLVE THE OVERDETERMINED SYSTEM'
169  ! write(lp,150)ID1,IM1,IC1,IY1,ID2,IM2,IC2,IY2
170 150 format(2i3,2i2,5x,2i3,2i2)
171 
172  WRITE(lp,255)tide_nti
173 255 FORMAT('NUMBER OF POINTS IN THE ANALYSIS =',i6)
174  write(lp,*) ' nin=',tide_nin
175 
176  DO idef=1,ndef
177 
178  WRITE(lp,52) res(idef),ssq(idef)
179 52 FORMAT('LARGEST RESIDUAL MAGNITUDE & RESIDUAL SUM OF SQUARES:' &
180  ,2e12.5)
181  WRITE(lp,66) sdev0(idef),rmsr0(idef)
182 66 FORMAT( &
183  'ST. DEV. OF RIGHT HAND SIDES OF ORIGINAL OVERDETERMINED SYSTEM:' &
184  ,e12.5/ &
185  ' AND THE ROOT MEAN SQUARE RESIDUAL ERROR:' &
186  ,e12.5)
187 
188 
189  write(lp,*) ' rms residual: brute force =',rmsr(idef)
190  write(lp,*) ' max residual: ',resmax(idef),imax(idef)
191 
192  WRITE(lp,41)
193 41 FORMAT('HARMONIC ANALYSIS RESULTS: AMPLITUDES, PHASE LAGS, C, S, &
194  & amp SD estimates, t-test value')
195  ! write out results for constant term & linear trend
196 
197  DO i=1,tide_mf
198  WRITE(lp,43) tidecon_name(i),tide_freqc(i),tide_ampc(i,idef),tide_phg(i,idef), &
199  tide_sig1(i,idef),tide_sig2(i,idef),tide_sig3(i,idef),tide_ttest(i,idef)
200  END DO
201 
202 43 FORMAT(5x,a5,4x,f12.9,2x,f10.5,2x,f10.3,5x,4f8.3)
203 
204  !
205  !* INFERENCE results are given now
206  !
207  IF (tide_nin.GE.0) THEN
208  write(lp,*) ' INFERENCE RESULTS'
209  l=0
210  do k=1,tide_nin
211  do i=2,tide_mf
212  IF (tidecon_name(i).eq.tide_konan(k)) EXIT
213  END DO
214  i1=i
215  do k2=1,tide_ninf(k)
216  l=l+1
217  write(lp,79) tide_konin(k,k2),tide_sigin(k,k2),tide_ampci(k,k2,idef), &
218  tide_phgi(k,k2,idef)
219 79 format(5x,a5,4x,f12.9,15x,f10.4,5x,f10.4)
220  END DO
221  END DO
222  inftot=l
223  END IF
224 
225 
226  WRITE(lp,70)tide_nti,tide_mf*2,' ',xlat,xlon,sngl(sdev0),sngl(sdev)
227 70 format('N,m,LAT,LON,SDEV0,SDEV: ',2i10,a4,f9.4,f10.4,2f10.2)
228  !
229  IF (tide_nin.GT.0) THEN
230  WRITE(lp,71) rmsrp(idef)
231 71 FORMAT('ROOT MEAN SQUARE RESIDUAL ERROR AFTER INFERENCE IS', &
232  e15.6, //)
233  ELSE
234  WRITE(lp,72) rmsrp(idef)
235 72 FORMAT('RECALCULATED ROOT MEAN SQUARE RESIDUAL ERROR IS ', &
236  e15.6, //)
237  ENDIF
238  END DO
239 
240  END SUBROUTINE tide_write_results
241 
242  !/ ------------------------------------------------------------------- /
243  SUBROUTINE tide_find_indices_prediction(LIST,INDS,TIDE_PRMF)
244  !/ +-----------------------------------+
245  !/ | F. Ardhuin |
246  !/ | FORTRAN 90 |
247  !/ | Last update : 28-Feb-2013 |
248  !/ +-----------------------------------+
249  !/
250  !/ 29-Jun-2013 : Creation ( version 4.11 )
251  !/
252  ! 1. Purpose :
253  !
254  ! Finds indices of tidal constituents to be used for prediction
255  !
256  ! 3. Parameters :
257  !
258  ! Parameter list
259  ! ----------------------------------------------------------------
260  ! LIST Char I Array of tidal constituents names to be used
261  ! INDS I.A. O Array of indices
262  ! TIDE_PRMF I.A. O number of constituents to be used
263  ! ----------------------------------------------------------------
264  !
265  !
266  ! 4. Subroutines used :
267  !
268  ! None
269  !
270  ! 5. Called by :
271  !
272  ! ww3_prtide
273  !
274  ! 6. Error messages :
275  !
276  ! None.
277  !
278  ! 7. Remarks :
279  !
280  ! 8. Structure :
281  !
282  ! See source code.
283  !
284  ! 9. Switches :
285  !
286  ! !/S Enable subroutine tracing.
287  ! !/T Enable test output.
288  !
289  ! 10. Source code :
290  !
291  USE w3odatmd, ONLY: iaproc, naproc, naperr, napout
292  USE w3odatmd, ONLY: ndse, ndso
293 #ifdef W3_S
294  USE w3servmd, ONLY: strace
295 #endif
296  !/ ------------------------------------------------------------------- /
297  IMPLICIT NONE
298  !/
299  !/ ------------------------------------------------------------------- /
300  !/ Parameter list
301  !/
302  CHARACTER(LEN=100), INTENT(IN) :: LIST(70)
303  INTEGER, INTENT(OUT) :: INDS(70), TIDE_PRMF
304 
305  INTEGER J, FOUND
306 #ifdef W3_S
307  INTEGER, SAVE :: IENT = 0
308 #endif
309  !
310 #ifdef W3_S
311  CALL strace (ient, 'TIDE_FIND_INDICES_PREDICTION')
312 #endif
313  !
314  tide_prmf=0
315  IF (trim(list(1)).EQ.'VFAST' .OR. trim(list(1)).EQ.'FAST') THEN
316  DO j=1,tide_mf
317  inds(j)=j
318  END DO
319  tide_prmf = tide_mf
320  RETURN
321  END IF
322  !
323  DO WHILE (len_trim(list(tide_prmf+1)).NE.0)
324  tide_prmf=tide_prmf+1
325  found = 0
326  DO j=1,tide_mf
327  IF (trim(tidecon_name(j)).EQ.trim(list(tide_prmf))) THEN
328  inds(tide_prmf)=j
329  found=1
330  IF (iaproc.EQ.napout) WRITE(ndso,'(A,A,E12.2)') 'Tidal constituent to be used in pre:', &
331  trim(list(tide_prmf)),tide_freqc(j)
332  END IF
333  END DO
334  IF (found.EQ.0 .AND. iaproc.EQ.napout) WRITE(ndso,'(3A)') 'Tidal constituent ',trim(list(tide_prmf)), &
335  ' not available.'
336  END DO
337  !
338  END SUBROUTINE tide_find_indices_prediction
339 
340 
341  !/ ------------------------------------------------------------------- /
342  SUBROUTINE tide_find_indices_analysis(LIST)
343  !/ +-----------------------------------+
344  !/ | F. Ardhuin |
345  !/ | FORTRAN 90 |
346  !/ | Last update : 21-Apr-2020 |
347  !/ +-----------------------------------+
348  !/
349  !/ 29-Jun-2013 : Creation ( version 4.11 )
350  !/ 21-Apr-2020 : Add 5 additional tidal const. ( version 7.13 )
351  !/
352  ! 1. Purpose :
353  !
354  ! Finds indices of tidal constituents to be used for analysis
355  !
356  ! 3. Parameters :
357  !
358  ! Parameter list
359  ! ----------------------------------------------------------------
360  ! LIST Char I Array of tidal constituents names to be used
361  ! ----------------------------------------------------------------
362  !
363  !
364  ! 4. Subroutines used :
365  !
366  ! None
367  !
368  ! 5. Called by :
369  !
370  ! ww3_prtide
371  !
372  ! 6. Error messages :
373  !
374  ! None.
375  !
376  ! 7. Remarks :
377  !
378  ! 8. Structure :
379  !
380  ! See source code.
381  !
382  ! 9. Switches :
383  !
384  ! !/S Enable subroutine tracing.
385  ! !/T Enable test output.
386  !
387  ! 10. Source code :
388  !
389  USE w3odatmd, ONLY: iaproc, naproc, naperr, napout
390  USE w3odatmd, ONLY: ndse, ndso
391 #ifdef W3_S
392  USE w3servmd, ONLY: strace
393 #endif
394  !
395  !/ ------------------------------------------------------------------- /
396  IMPLICIT NONE
397  !/
398  !/ ------------------------------------------------------------------- /
399  !/ Parameter list
400  !/
401  CHARACTER(LEN=100), INTENT(IN) :: LIST(70)
402  !
403  INTEGER TIDE_MF_ALL
404  CHARACTER(LEN=5) :: TIDECON_NAME_ALL(65) ! array of names of tidal constituents
405  REAL :: TIDE_FREQC_ALL(65) ! array of freq. of tidal constituents
406  INTEGER :: INDS(65), J, FOUND, NTIDES
407 #ifdef W3_S
408  INTEGER, SAVE :: IENT = 0
409 #endif
410  !
411 #ifdef W3_S
412  CALL strace (ient, 'TIDE_FIND_INDICES_PREDICTION')
413 #endif
414  !
415  tidecon_name_all(:)=(/ &
416  'Z0 ', 'SSA ', 'MSM ', 'MM ', 'MSF ', 'MF ', 'ALP1 ', '2Q1 ', 'SIG1 ', 'Q1 ', &
417  'RHO1 ', 'O1 ', 'TAU1 ', 'BET1 ', 'NO1 ', 'CHI1 ', 'P1 ', 'K1 ', 'PHI1 ', 'THE1 ', &
418  'J1 ', 'SO1 ', 'OO1 ', 'UPS1 ', 'OQ2 ', 'EPS2 ', '2N2 ', 'MU2 ', 'N2 ', 'NU2 ', &
419  'M2 ', 'MKS2 ', 'LDA2 ', 'L2 ', 'S2 ', 'K2 ', 'MSN2 ', 'ETA2 ', 'MO3 ', 'M3 ', &
420  'SO3 ', 'MK3 ', 'SK3 ', 'MN4 ', 'M4 ', 'SN4 ', 'MS4 ', 'MK4 ', 'S4 ', 'SK4 ', &
421  '2MK5 ', '2SK5 ', '2MN6 ', 'M6 ', '2MS6 ', '2MK6 ', '2SM6 ', 'MSK6 ', '3MK7 ', 'M8 ', &
422  'N4 ', 'R2 ', 'S1 ', 'SA ', 'T2 ' /)
423  !
424  tide_freqc_all(:)=(/0.0000000000, 0.0002281591, 0.0013097807, 0.0015121520, 0.0028219327, 0.0030500918, &
425  0.0343965698, 0.0357063505, 0.0359087218, 0.0372185025, 0.0374208738, &
426  0.0387306544, 0.0389588136, 0.0400404351, 0.0402685943, 0.0404709655, &
427  0.0415525871, 0.0417807462, 0.0420089053, 0.0430905269, 0.0432928982, &
428  0.0446026789, 0.0448308380, 0.0463429900, 0.0759749448, 0.0761773160, &
429  0.0774870967, 0.0776894680, 0.0789992487, 0.0792016200, 0.0805114007, &
430  0.0807395598, 0.0818211814, 0.0820235526, 0.0833333333, 0.0835614924, &
431  0.0848454853, 0.0850736444, 0.1192420551, 0.1207671010, 0.1220639878, &
432  0.1222921469, 0.1251140796, 0.1595106494, 0.1610228013, 0.1623325820, &
433  0.1638447340, 0.1640728931, 0.1666666667, 0.1668948258, 0.2028035476, &
434  0.2084474129, 0.2400220500, 0.2415342020, 0.2443561347, 0.2445842938, &
435  0.2471780673, 0.2474062264, 0.2833149482, 0.3220456027, &
436  0.157998497 , 0.083447407 , 0.041666667 , 0.000114080 , 0.083219259 /)
437 
438  inds(:) = 0
439  !
440  IF (trim(list(1)).EQ.'FAST') THEN
441  tide_mf = 44
442  inds(1:44)= (/ 1, 2, 3, 4, 5, 6, 12, 17, 18, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, &
443  37, 38, 39, 40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 50, 51, 52, 53, 54,&
444  55, 56, 57, 58, 59, 60 /)
445  !
446  ELSE IF (trim(list(1)).EQ.'VFAST') THEN
447  tide_mf = 20
448  inds(1:20)= (/ 1, 2, 3, 5, 6, 27, 28, 29, 30, 31, 35, 36, 37, 44, 45, 47, 49, 54, 55, 60 /)
449  !
450  ELSE
451  tide_mf=0
452  ntides=0
453  !
454  DO WHILE (len_trim(list(tide_mf+1)).NE.0)
455  !
456  tide_mf=tide_mf+1
457  found = 0
458  DO j=1,65
459  IF (trim(tidecon_name_all(j)).EQ.trim(list(tide_mf))) THEN
460  ntides=ntides+1
461  inds(ntides)=j
462  found = 1
463  IF (iaproc.EQ.napout) WRITE(ndso,'(A,I4,2A,E12.2)') &
464  'Tidal constituent in analysis:', j, ' ', &
465  trim(tidecon_name_all(j)),tide_freqc_all(j)
466  END IF
467  END DO
468  IF (found.EQ.0 .AND. iaproc.EQ.napout) WRITE(ndso,'(A,I4,A,A)') &
469  'Tidal constituent ',tide_mf,trim(list(tide_mf)),' not available.'
470  !
471  END DO
472  !
473  tide_mf=ntides
474  END IF
475  !
476  ! Defines names and frequencies
477  !
478  IF (ALLOCATED(tide_freqc)) DEALLOCATE(tide_freqc)
479  ALLOCATE(tide_freqc(tide_mf),tidecon_name(tide_mf))
480 
481  DO j=1,tide_mf
482  tidecon_name(j) = tidecon_name_all(inds(j))
483  tide_freqc(j) = tide_freqc_all(inds(j))
484  END DO
485  CALL tide_set_indices
486 
487  !
488  END SUBROUTINE tide_find_indices_analysis
489  !/ ------------------------------------------------------------------- /
490 
491 
492 
493  !/ ------------------------------------------------------------------- /
494  SUBROUTINE tide_set_indices
495  !
496  IMPLICIT NONE
497  !
498  INTEGER J, K, K1, L, J1, JL, L2, KM1, JBASE
499  !
500  DO l=1,tide_mf
501  DO k=1,ntotal_con
502  IF (tidecon_allnames(k).EQ.tidecon_name(l)) tide_index2(l)=k
503  END DO
504  END DO
505  !
506  tide_indexj(:)=0
507  tide_indexjk(:)=0
508  jbase=0
509  k1=ntidal_con+1
510  !
511  DO k=k1,ntotal_con
512  j1=jbase+1
513  tide_indexj(k)=j1
514  jl=jbase+nj(k)
515  DO j=j1,jl
516  km1=k-1
517  l2=0
518  DO l2=1,km1
519  IF (tidecon_allnames(l2).EQ.konco_con(j)) THEN
520  tide_indexjk(j)=l2
521  END IF
522  END DO ! L2
523  END DO ! J
524  jbase=jl
525  END DO ! K
526  !
527  END SUBROUTINE tide_set_indices
528 
529 
530  !/ ------------------------------------------------------------------- /
531  SUBROUTINE setvuf_fast(h,pp,s,p,enp,dh,dpp,ds,dp,dnp,tau,XLAT,F,U,V)
532  ! setvuf calculates the V,u,f values at time hr for all constituents
533  IMPLICIT NONE
534  REAL, INTENT(IN) :: XLAT
535  REAL(KIND=8), intent(in) :: h,pp,s,p,enp,dh,dpp,ds,dp,dnp,tau
536  INTEGER, PARAMETER :: NTIL=44
537  REAL , INTENT(OUT) :: F(NTIL),U(NTIL),V(NTIL)
538  REAL :: FA(170),UA(170),VA(170)
539  !
540  ! Local variables
541  !
542  INTEGER :: K, L, L2, JBASE, J1, J, JL, K1
543  REAL(KIND=4), parameter :: pi=3.1415926536
544  REAL(KIND=4), parameter :: twopi=2.*3.1415926536
545 
546  REAL :: SLAT, VDBL, RR, SUMC, SUMS, UUDBL, UU, CXLAT
547  INTEGER :: IUU, IV
548  ! This comment was taken from t_tide, a matlab tidal prediction suite
549  !
550  ! Apparently the second-order terms in the tidal potential go to zero
551  ! at the equator, but the third-order terms do not. Hence when trying
552  ! to infer the third-order terms from the second-order terms, the
553  ! nodal correction factors blow up. In order to prevent this, it is
554  ! assumed that the equatorial forcing is due to second-order forcing
555  ! OFF the equator, from about the 5 degree location. Latitudes are
556  ! hence (somewhat arbitrarily) forced to be no closer than 5 deg to
557  ! the equator, as per note in Foreman.
558  cxlat = max(abs(xlat), 5.)
559  slat=sin(pi*cxlat/180.)
560 
561  jbase=0
562  ! All
563  ! 1'Z0 *','SA ','SSA *','MSM *','MM *,'MSF *','MF *','ALP1*','2Q1 ','SIG1 ', &
564  !11'Q1 ','RHO1 ','O1 *','TAU1 ','BET1 ','NO1 ','CHI1 ','PI1 ','P1 *','S1 ', &
565  !21'K1 *','PSI1 ','PHI1 ','THE1 ','J1 ','OO1 ','UPS1 ','OQ2 ','EPS2*','2N2 *', &
566  !31'MU2 *','N2 *','NU2 *','GAM2 ','H1 ','M2 *','H2 ','LDA2*','L2 *','T2 ', &
567  !41'S2 *','R2 ','K2 *','ETA2 ','M3 ','2PO1 ','SO1 ','ST36 ','2NS2 ','ST37 ', &
568  !51'ST1 ','ST2 ','ST3 ','O2 ','ST4 ','SNK2 ','OP2 ','MKS2*','ST5 ','ST6 ', &
569  !61'2SK2 ','MSN2 ','ST7 ','2SM2 ','ST38 ','SKM2 ','2SN2 ','NO3 ','MO3 ','NK3 ', &
570  !71'SO3 ','MK3 ','SP3 ','SK3 ','ST8 ','N4 ','3MS4 ','ST39 ','MN4 ','ST40 ', &
571  !81'ST9 ','M4 ','ST10 ','SN4 ','KN4 ','MS4 ','MK4 ','SL4 ','S4 ','SK4 ', &
572  !91'MNO5 ','2MO5 ','3MP5 ','MNK5 ','2MP5 ','2MK5 ','MSK5 ','3KM5 ','2SK5 ','ST11 ', &
573  !01'2NM6 ','ST12 ','ST41 ','2MN6 ','ST13 ','M6 ','MSN6 ','MKN6 ','2MS6 ','2MK6 ', &
574  !11'NSK6 ','2SM6 ','MSK6 ','ST42 ','S6 ','ST14 ','ST15 ','M7 ','ST16 ','3MK7 ', &
575  !21'ST17 ','ST18 ','3MN8 ','ST19 ','M8 ','ST20 ','ST21 ','3MS8 ','3MK8 ','ST22 ', &
576  !31'ST23 ','ST24 ','ST25 ','ST26 ','4MK9 ','ST27 ','ST28 ','M10 ','ST29 ','ST30 ', &
577  !41'ST31 ','ST32 ','ST33 ','M12 ','ST34 ','ST35 '/)
578 
579  ! Possible
580  ! 'Z0 ', 'SSA ', 'MSM ', 'MM ', 'MSF ', 'MF ', 'ALP1 ', '2Q1 ', 'SIG1 ', 'Q1 ', &
581  ! 'RHO1 ', 'O1 ', 'TAU1 ', 'BET1 ', 'NO1 ', 'CHI1 ', 'P1 ', 'K1 ', 'PHI1 ', 'THE1 ', &
582  ! 'J1 ', 'SO1 ', 'OO1 ', 'UPS1 ', 'OQ2 ', 'EPS2 ', '2N2 ', 'MU2 ', 'N2 ', 'NU2 ', &
583  ! 'M2 ', 'MKS2 ', 'LDA2 ', 'L2 ', 'S2 ', 'K2 ', 'MSN2 ', 'ETA2 ', 'MO3 ', 'M3 ', &
584  ! 'SO3 ', 'MK3 ', 'SK3 ', 'MN4 ', 'M4 ', 'SN4 ', 'MS4 ', 'MK4 ', 'S4 ', 'SK4 ', &
585  ! '2MK5 ', '2SK5 ', '2MN6 ', 'M6 ', '2MS6 ', '2MK6 ', '2SM6 ', 'MSK6 ', '3MK7 ', 'M8 ' /)
586 
587  ! Subset
588  ! 'Z0 ', 'SSA ', 'MSM ', 'MM ', 'MSF ', 'MF ', &
589  ! 'O1 ', 'P1 ', 'K1 ' , &
590  ! 'EPS2 ', '2N2 ', 'MU2 ', 'N2 ', 'NU2 ', &
591  ! 'M2 ', 'MKS2 ', 'LDA2 ', 'L2 ', 'S2 ', 'K2 ', 'MSN2 ', 'ETA2 ', 'MO3 ', 'M3 ', &
592  ! 'SO3 ', 'MK3 ', 'SK3 ', 'MN4 ', 'M4 ', 'SN4 ', 'MS4 ', 'MK4 ', 'S4 ', 'SK4 ', &
593  ! '2MK5 ', '2SK5 ', '2MN6 ', 'M6 ', '2MS6 ', '2MK6 ', '2SM6 ', 'MSK6 ', '3MK7 ', 'M8 ' /)
594 
595 
596  !
597  jbase=0
598 
599  ! initialize arrays to avoid NaN values
600  fa(:)=0
601  ua(:)=0
602  va(:)=0
603 
604  DO k=1,ntidal_con
605  j1=jbase+1
606  jl=jbase+nj(k)
607  DO l=1,tide_mf
608  IF (tide_index2(l).EQ.k) THEN
609  vdbl=ii(k)*tau+jj(k)*s+kk(k)*h+ll(k)*p+mm(k)*enp+nn(k)*pp+semi(k)
610  iv=vdbl
611  iv=(iv/2)*2
612  sumc=1.
613  sums=0.
614  DO j=j1,jl
615  ! ITIME ???
616  !***********************************************************************
617  !* HERE THE SATELLITE AMPLITUDE RATIO ADJUSTMENT FOR LATITUDE IS MADE
618  !
619  rr=ee(j)
620  l2=ir(j)+1
621  IF (l2.EQ.2) THEN
622  rr=ee(j)*0.36309*(1.-5.*slat*slat)/slat
623  ELSE IF (l2.EQ.3) THEN
624  rr=ee(j)*2.59808*slat
625  END IF
626  uudbl=ldel(j)*p+mdel(j)*enp+ndel(j)*pp+ph(j)
627  iuu=uudbl
628  uu=uudbl-iuu
629  sumc=sumc+rr*cos(uu*twopi)
630  sums=sums+rr*sin(uu*twopi)
631  END DO
632  !
633  fa(k)=sqrt(sumc*sumc+sums*sums)
634  va(k)=vdbl-iv
635  ua(k)=atan2(sums,sumc)/twopi
636  END IF
637  jbase=jl ! (indx(L).EQ.K)
638  END DO ! L
639  END DO ! K
640  !
641  ! HERE F AND V+U OF THE SHALLOW WATER CONSTITUENTS ARE COMPUTED FROM
642  ! THE VALUES OF THE MAIN CONSTITUENT FROM WHICH THEY ARE DERIVED.
643  !
644  k1=ntidal_con+1
645  IF (k1.GT.ntotal_con) RETURN
646  !
647  DO k=k1,ntotal_con
648  fa(k)=1.0
649  va(k)=0.0
650  ua(k)=0.
651  DO j=tide_indexj(k),tide_indexj(k)+nj(k)-1
652  l2=tide_indexjk(j)
653  fa(k)=fa(k)*fa(l2)**abs(coef_con(j))
654  va(k)=va(k)+coef_con(j)*va(l2)
655  ua(k)=ua(k)+coef_con(j)*ua(l2)
656  END DO ! J
657  END DO ! K
658  !
659  DO l=1,tide_mf
660  f(l)=fa(tide_index2(l))
661  u(l)=ua(tide_index2(l))
662  v(l)=va(tide_index2(l))
663  END DO ! L
664 
665  RETURN
666  !
667  END SUBROUTINE setvuf_fast
668 
669 
670  !/ ------------------------------------------------------------------- /
671  SUBROUTINE tide_read_settings(filename,fnam6,fnam7,fnam8,fnam9,fnam11)
672  IMPLICIT NONE
673  CHARACTER*256, INTENT(IN) :: filename
674  CHARACTER*256, INTENT(OUT) :: fnam6,fnam7,fnam8,fnam9,fnam11
675 
676  INTEGER KIN
677 
678  ! Parameters for reading KR1
679  ! FILE I/O
680  ! KIN is the master input file.
681  ! fnam6 is the file to which the output is sent. It is assigned the number
682  ! lp.
683  ! fnam7 is file containing the constituents to be included in the analysis,
684  ! the analysis period, inference parameters, the flag controlling
685  ! height or current analyses, and site information. It is assigned the
686  ! number kr1.
687  ! fnam8 is the file containing all the astronomical argument information
688  ! (it should not have to be changed)
689  ! fnam11 is a file to which information on the SVD matrix fit is output when
690  !
691  ! Original, fitted, and residual time series are output to file 25 while
692  ! the same are also output to file 26 in a format that could be input to
693  ! Excel for plotting.
694  !
695  kin=40 ! Input file assigned to unit 4
696 
697  ! OPEN(UNIT=KIN,FILE='tuk75_tidana.inp',STATUS='OLD')
698  ! OPEN(UNIT=KIN,FILE='victoria_2008_test.inp',STATUS='OLD')
699  OPEN(unit=kin,file=filename,status='OLD')
700  ! OPEN(UNIT=KIN,FILE='kiw05_mar2008.inp',STATUS='OLD')
701  ! OPEN(UNIT=KIN,FILE='tcs05_sep07-mar08.inp',STATUS='OLD')
702  read(kin,'(a)') fnam6
703  read(kin,'(a)') fnam7
704  read(kin,'(a)') fnam8
705  read(kin,'(a)') fnam9
706  read(kin,'(a)') fnam11
707  ! open(unit=11,file=fnam11,status='unknown',form='formatted')
708  ! unit 25 stores the residual time series
709  !read(KIN,'(a)') fname
710  !open(unit=25,file=fname,status='unknown',form='formatted')
711  !read(KIN,'(a)') fname
712  !open(unit=26,file=fname,status='unknown',form='formatted')
713  END SUBROUTINE tide_read_settings
714 
715  !/ ------------------------------------------------------------------- /
716  SUBROUTINE tide_read_anapar(KR1,LP,filename,KD1,KD2,XLON,XLAT,NDEF,ITREND,ITZ)
717  ! Parameters for reading KR1
718  IMPLICIT NONE
719 
720  INTEGER, INTENT(IN) :: KR1, LP
721  CHARACTER*256, INTENT(IN) :: filename
722  INTEGER(KIND=4), INTENT(OUT) :: KD1, KD2
723  INTEGER , INTENT(OUT) :: NDEF,ITREND
724  REAL, INTENT(OUT) :: XLON, XLAT
725  CHARACTER*4 , INTENT(OUT) :: ITZ
726  !
727  INTEGER :: I, IY
728  INTEGER ID1,IM1,IY1,ID2,IM2,IY2,IC1,IC2
729  INTEGER JSTN, LATD,LATM,LOND,LONM, K, K2
730  CHARACTER*4 NSTN(5)
731 
732  open(unit=kr1,file=filename,status='old',form='formatted')
733 
734  !
735  !*
736  !***********************************************************************
737  !* READ FROM DEVICE KR1 THE ANALYSIS TYPE AND TIDAL STATION DETAILS.
738  !*
739  !* 1)ONE RECORD FOR THE VARIABLES TIDE_MF
740  !* TIDE_MF = THE NUMBER OF CONSTITUENTS, INCLUDING THE CONSTANT TERM Z0,
741  !* TO BE IN THE LEAST SQUARES FIT.
742  !*
743  !* 2)ONE RECORD FOR EACH OF THE TIDE_MF CONSTITUENTS TO BE INCLUDED IN THE
744  !* FIT. EACH RECORD CONTAINS THE VARIABLES NAME AND TIDE_FREQC IN THE
745  !* FORMAT (A5,2X,F13.10). NAME IS THE CONSTITUENT NAME, WHICH SHOULD
746  !* BE LEFT JUSTIFIED IN THE ALPHANUMERIC FIELD, WHILE TIDE_FREQC IS ITS
747  !* FREQUENCY MEASURED IN CYCLES PER HOUR.
748  !*
749  !* 3) ONE RECORD IN THE FORMAT (8I5) CONTAINING THE FOLLOWING
750  !* INFORMATION ON THE TIME PERIOD OF THE ANALYSIS.
751  !* ID1,IM1,IY1 - DAY,MONTH,YEAR OF THE BEGINNING OF THE ANALYSIS
752  !* PERIOD,
753  !* ID2,IM2,IY2 - DAY,MONTH,YEAR OF THE END OF THE ANALYSIS PERIOD.
754  !* IC1,IC2 - CENTURY OFR THE BEGINNING AND END OF THE ANALYSIS
755  !* PERIOD (ZERO VALUES ARE RESET TO 19)
756  !*
757  !*
758  !* 4)ONE RECORD IN THE FORMAT (I5,5A4,1X,A4,4I5) CONTAINING THE
759  !* FOLLOWING TIDAL STATION INFORMATION.
760  !* JSTN = TIDAL STATION NUMBER,
761  !* (NSTN(I),I=1,5) = TIDAL STATION NAME,
762  ! * ITZ = TIME ZONE IN WHICH THE OBSERVATIONS WERE RECORDED,
763  ! * LATD,LATM = STATION LATITUDE IN DEGREES AND MINUTES,
764  ! * LOND,LONM = STATION LONGITUDE IN DEGREES AND MINUTES.
765  ! *
766  ! * 5)ONE SET RECORDS FOR EACH POSSIBLE INFERENCE. THE FIRST RECORD HAS THE
767  ! * CONSITUENT NAME, ITS FREQUENCY, AND THE NUMBER OF CONSTITUENTS TO BE
768  ! * INFERRED (4X,A5,E16.10,i5), WHILE THERE IS ONE RECORD FOR EACH OF THE
769  ! * CONSTITUENTS TO BE INFERRED WITH THE NAME, FREQUENCY, AMPLITUDE RATIO
770  ! * (INFERRED TO REFERENCE) AND PHASE DIFFERENCE (GREENWICH PHASE LAG OF
771  ! * THE INFERRED CONSTITUENT SUBTRACTED FROM THE GREENWICH PHASE LAG OF THE
772  ! * (ANALYSED CONSTITUENT IN THE FORMAT(4X,A5,E16.10,2F10.3)
773  ! *
774  ! * FOR KR1 INPUT, ALL CONSTITUENT NAMES SHOULD BE LEFT JUSTIFIED IN
775  ! * THE ALPHANUMERIC FIELD, FREQUENCIES ARE MEASURED IN CYCLES/HOUR, AND
776  ! * ALL CONSTITUENTS MUST BE INCLUDED IN THE LIST IN READ FROM FNAM8.
777  !
778  ! write(6,*) ' reading from unit kr1'
779  READ(kr1,*) tide_mf,ndef,itrend
780  ALLOCATE(tide_freqc(tide_mf),tidecon_name(tide_mf))
781  ! ndef=1 if only 1D field to be analysed (eg., elevations)
782  ! ndef=2 if 2D field: velocity components, EW followed by NS : this is now de-activated
783 #ifdef W3_T
784  WRITE(6,*) ' number of constituents & degrees of freedom=',tide_mf,ndef
785 #endif
786  IF (itrend.eq.1) then
787 #ifdef W3_T
788  WRITE(6,*) ' a linear trend is included in the analysis'
789 #endif
790  else
791 #ifdef W3_T
792  WRITE(6,*) ' no linear trend is included'
793 #endif
794  END IF
795  ! TIDE_MF= number of consituents, excluding linear trend. The constant
796  ! term, Z0 should be first in the list.
797  ! itrend= 1 if include linear trend
798  ! itrend= otherwise, no trend
799  ! number of unknowns, M, depends on whether we have a linear trend
800 
801 10 FORMAT(2i5,f5.2)
802  READ(kr1,11) (tidecon_name(i),tide_freqc(i),i=1,tide_mf)
803 #ifdef W3_T
804  WRITE(6,*) (tidecon_name(i),tide_freqc(i),i=1,tide_mf)
805 #endif
806 11 FORMAT(4x,a5,f16.10)
807  READ(kr1,7) id1,im1,iy1,id2,im2,iy2,ic1,ic2
808 #ifdef W3_T
809  WRITE(6,*) id1,im1,iy1,id2,im2,iy2,ic1,ic2
810 #endif
811  IF (ic1.EQ.0) ic1=19
812  IF (ic2.EQ.0) ic2=19
813 7 FORMAT(16i5)
814  READ(kr1,9) jstn,nstn(1:5),itz,latd,latm,lond,lonm
815 9 FORMAT(i5,5a4,1x,a4,4i5)
816 
817  iy=ic1*100+iy1
818  kd1=juldayt(id1,im1,iy)
819 
820  iy=ic2*100+iy2
821  kd2=juldayt(id2,im2,iy)
822  !
823  ! read in inference information now as it will be used in the lsq matrix
824  !
825  DO k=1,10
826  READ(kr1,'(4X,A5,E17.10,i5)')tide_konan(k),tide_sigan(k),tide_ninf(k)
827  ! write(6,1010)TIDE_KONAN(K),TIDE_SIGAN(K),TIDE_NINF(k)
828  IF (tide_konan(k).EQ.kblank) EXIT
829  do k2=1,tide_ninf(k)
830  read(kr1,'(4X,A5,E17.10,2F10.3)') tide_konin(k,k2),tide_sigin(k,k2),tide_r(k,k2),tide_zeta(k,k2)
831  END DO
832  END DO
833  tide_nin=k-1
834  CLOSE(kr1)
835 
836  xlat=latd+latm/60.
837  xlon=lond+lonm/60.
838 
839  RETURN
840  END SUBROUTINE tide_read_anapar
841 
842  !/ ------------------------------------------------------------------- /
843  SUBROUTINE tide_read_timeseries(KR2,filename,KD1,KD2,TIDE_NTI,NDEF)
844  !
845  IMPLICIT NONE
846  !
847  INTEGER, INTENT(IN) :: KR2, NDEF
848  CHARACTER*256, INTENT(IN) :: filename
849  INTEGER(KIND=4), INTENT(IN) :: KD1,KD2
850  INTEGER(KIND=4), INTENT(OUT) :: TIDE_NTI
851  !
852  INTEGER :: I, idd,imm,icc,iyy,ihh,imin,isec,iy
853  REAL, ALLOCATABLE :: TIDE_DATATMP(:,:)
854  INTEGER(KIND=4), ALLOCATABLE :: TIDE_DAYSTMP(:), TIDE_SECSTMP(:)
855  INTEGER(KIND=4) :: KDD
856  INTEGER :: ICODE
857  REAL :: htt(NDEF)
858  !
859  ! Initialize Variables
860  !
861  ALLOCATE( tide_datatmp(nr,ndef), tide_daystmp(nr), tide_secstmp(nr) )
862 
863  ! Reads in data between dates kd1 and kd2 in file kr2
864  !
865  OPEN(unit=kr2,file=filename,status='old',form='formatted')
866 
867  icode = 0
868  kdd = kd1
869  i=0
870  DO WHILE(icode.EQ.0.AND.kdd.LE.kd2)
871  !
872  ! reads with the original Foreman's format
873  !
874  READ(kr2,145,iostat=icode) idd,imm,icc,iyy,ihh,imin,htt(1:ndef)
875 145 format(6i2,4f10.4)
876  isec=0
877  iy=icc*100+iyy
878 
879  kdd=juldayt(idd,imm,iy)
880 
881  IF (kdd.lt.kd1) then
882  WRITE(*,*) icc,iyy,imm,idd,ihh,imin
883  WRITE(*,*)'kd, kd1, kd2 =',kdd,kd1,kd2
884  write(*,*) ' observation before analysis period'
885  ELSE
886  !
887  ! Fills in data array
888  !
889  IF (icode.EQ.0.AND.kdd.LE.kd2) THEN
890  i=i+1
891  tide_datatmp(i,:)=htt(:)
892  tide_daystmp(i)=kdd
893  tide_secstmp(i)=ihh*3600+imin*60+isec
894  END IF
895  END IF
896  END DO
897 
898  tide_nti=i
899  ALLOCATE( tide_data(tide_nti,ndef) )
900  ALLOCATE( tide_days(tide_nti), tide_secs(tide_nti), tide_hours(tide_nti) )
901  tide_data(1:tide_nti,1:ndef)=tide_datatmp(1:tide_nti,1:ndef)
902  tide_days(1:tide_nti)=tide_daystmp(1:tide_nti)
903  tide_secs(1:tide_nti)=tide_secstmp(1:tide_nti)
904  tide_hours(:)=24.d0*dfloat(tide_days(:))+dfloat(tide_secs(:))/3600.d0
905  CLOSE(kr2)
906  RETURN
907  END SUBROUTINE tide_read_timeseries
908 
909  !/ ------------------------------------------------------------------- /
910  SUBROUTINE astr(d1,h,pp,s,p,np,dh,dpp,ds,dp,dnp)
911  ! this subroutine calculates the following five ephermides
912  ! of the sun and moon
913  ! h = mean longitude of the sum
914  ! pp = mean longitude of the solar perigee
915  ! s = mean longitude of the moon
916  ! p = mean longitude of the lunar perigee
917  ! np = negative of the longitude of the mean ascending node
918  ! and their rates of change.
919  ! Units for the ephermides are cycles and for their derivatives
920  ! are cycles/365 days
921  ! The formulae for calculating this ephermides were taken from
922  ! pages 98 and 107 of the Explanatory Supplement to the
923  ! Astronomical Ephermeris and the American Ephermis and
924  ! Nautical Almanac (1961)
925  !
926  implicit none
927  REAL(KIND=8), intent(in ):: d1
928  REAL(KIND=8), intent(out):: h,pp,s,p,np,dh,dpp,ds,dp,dnp
929  !
930  ! Local variables
931  !
932  REAL(KIND=8) :: d2,f,f2
933 
934  d2=d1*1.d-4
935  f=360.d0
936  f2=f/365.d0
937  h=279.696678d0+.9856473354d0*d1+.00002267d0*d2*d2
938  pp=281.220833d0+.0000470684d0*d1+.0000339d0*d2*d2+&
939  .00000007d0*d2**3
940  s=270.434164d0+13.1763965268d0*d1-.000085d0*d2*d2+&
941  .000000039d0*d2**3
942  p=334.329556d0+.1114040803d0*d1-.0007739d0*d2*d2-&
943  .00000026d0*d2**3
944  np=-259.183275d0+.0529539222d0*d1-.0001557d0*d2*d2-&
945  .00000005d0*d2**3
946  h=h/f
947  pp=pp/f
948  s=s/f
949  p=p/f
950  np=np/f
951  h=h-dint(h)
952  pp=pp-dint(pp)
953  s=s-dint(s)
954  p=p-dint(p)
955  np=np-dint(np)
956  dh=.9856473354d0+2.d-8*.00002267d0*d1
957  dpp=.0000470684d0+2.d-8*.0000339d0*d1&
958  +3.d-12*.00000007d0*d1**2
959  ds=13.1763965268d0-2.d-8*.000085d0*d1+&
960  3.d-12*.000000039d0*d1**2
961  dp=.1114040803d0-2.d-8*.0007739d0*d1-&
962  3.d-12*.00000026d0*d1**2
963  dnp=+.0529539222d0-2.d-8*.0001557d0*d1-&
964  3.d-12*.00000005d0*d1**2
965  dh=dh/f2
966  dpp=dpp/f2
967  ds=ds/f2
968  dp=dp/f2
969  dnp=dnp/f2
970  return
971  end SUBROUTINE astr
972 
973 
974  !/ ------------------------------------------------------------------- /
975  ! Note by FA: should try to replace with standard distance d= sqrt(a^2+b^2)
976  FUNCTION dpythag(a,b)
977  DOUBLE PRECISION a,b,dpythag
978  DOUBLE PRECISION absa,absb
979  absa=abs(a)
980  absb=abs(b)
981  IF (absa.gt.absb)then
982  dpythag=absa*sqrt(1.0d0+(absb/absa)**2)
983  else
984  IF (absb.eq.0.0d0)then
985  dpythag=0.0d0
986  else
987  dpythag=absb*sqrt(1.0d0+(absa/absb)**2)
988  endif
989  endif
990  return
991  END FUNCTION dpythag
992 
993  !/ ------------------------------------------------------------------- /
994  ! (C) Copr. 1986-92 Numerical Recipes Software '%1&&Yw^2.
995  SUBROUTINE dsvbksb(u,w,v,m,n,mp,np,b,x)
996  INTEGER m,mp,n,np,NMAX
997  DOUBLE PRECISION b(mp),u(mp,np),v(np,np),w(np),x(np)
998  parameter(nmax=500)
999  INTEGER i,j,jj
1000  DOUBLE PRECISION s,tmp(NMAX)
1001  do j=1,n
1002  s=0.0d0
1003  IF (w(j).ne.0.0d0)then
1004  do i=1,m
1005  s=s+u(i,j)*b(i)
1006  end do
1007  s=s/w(j)
1008  endif
1009  tmp(j)=s
1010  end do
1011  do j=1,n
1012  s=0.0d0
1013  do jj=1,n
1014  s=s+v(j,jj)*tmp(jj)
1015  end do
1016  x(j)=s
1017  end do
1018  return
1019  END SUBROUTINE dsvbksb
1020 
1021 
1022  !/ ------------------------------------------------------------------- /
1023  SUBROUTINE dsvdcmp(a,m,n,mp,np,w,v)
1024  INTEGER m,mp,n,np,NMAX
1025  DOUBLE PRECISION a(mp,np),v(np,np),w(np)
1026  parameter(nmax=500)
1027 
1028  INTEGER i,its,j,jj,k,l,nm
1029  DOUBLE PRECISION anorm,c,f,g,h,s,scale,x,y,z,rv1(NMAX)
1030 
1031  g=0.0d0
1032  scale=0.0d0
1033  anorm=0.0d0
1034  do i=1,n
1035  l=i+1
1036  rv1(i)=scale*g
1037  g=0.0d0
1038  s=0.0d0
1039  scale=0.0d0
1040  IF (i.le.m)then
1041  DO k=i,m
1042  scale=scale+abs(a(k,i))
1043  END DO
1044  IF (scale.ne.0.0d0) THEN
1045  DO k=i,m
1046  a(k,i)=a(k,i)/scale
1047  s=s+a(k,i)*a(k,i)
1048  END DO
1049  f=a(i,i)
1050  g=-sign(sqrt(s),f)
1051  h=f*g-s
1052  a(i,i)=f-g
1053  DO j=l,n
1054  s=0.0d0
1055  DO k=i,m
1056  s=s+a(k,i)*a(k,j)
1057  END DO
1058  f=s/h
1059  DO k=i,m
1060  a(k,j)=a(k,j)+f*a(k,i)
1061  END DO
1062  END DO
1063  DO k=i,m
1064  a(k,i)=scale*a(k,i)
1065  END DO
1066  !
1067  END IF
1068  !
1069  END IF
1070  !
1071  w(i)=scale *g
1072  g=0.0d0
1073  s=0.0d0
1074  scale=0.0d0
1075  IF ((i.le.m).and.(i.ne.n))then
1076  do k=l,n
1077  scale=scale+abs(a(i,k))
1078  end do
1079  IF (scale.ne.0.0d0)then
1080  do k=l,n
1081  a(i,k)=a(i,k)/scale
1082  s=s+a(i,k)*a(i,k)
1083  end do
1084  f=a(i,l)
1085  g=-sign(sqrt(s),f)
1086  h=f*g-s
1087  a(i,l)=f-g
1088  do k=l,n
1089  rv1(k)=a(i,k)/h
1090  end do
1091  do j=l,m
1092  s=0.0d0
1093  do k=l,n
1094  s=s+a(j,k)*a(i,k)
1095  end do
1096  do k=l,n
1097  a(j,k)=a(j,k)+s*rv1(k)
1098  end do
1099  end do
1100  do k=l,n
1101  a(i,k)=scale*a(i,k)
1102  end do
1103  endif
1104  endif
1105  anorm=max(anorm,(abs(w(i))+abs(rv1(i))))
1106  end do
1107  do i=n,1,-1
1108  IF (i.lt.n)then
1109  IF (g.ne.0.0d0)then
1110  do j=l,n
1111  v(j,i)=(a(i,j)/a(i,l))/g
1112  end do
1113  do j=l,n
1114  s=0.0d0
1115  do k=l,n
1116  s=s+a(i,k)*v(k,j)
1117  end do
1118  do k=l,n
1119  v(k,j)=v(k,j)+s*v(k,i)
1120  end do
1121  end do
1122  endif
1123  do j=l,n
1124  v(i,j)=0.0d0
1125  v(j,i)=0.0d0
1126  end do
1127  endif
1128  v(i,i)=1.0d0
1129  g=rv1(i)
1130  l=i
1131  end do
1132  do i=min(m,n),1,-1
1133  l=i+1
1134  g=w(i)
1135  do j=l,n
1136  a(i,j)=0.0d0
1137  end do
1138  IF (g.ne.0.0d0)then
1139  g=1.0d0/g
1140  do j=l,n
1141  s=0.0d0
1142  do k=l,m
1143  s=s+a(k,i)*a(k,j)
1144  end do
1145  f=(s/a(i,i))*g
1146  do k=i,m
1147  a(k,j)=a(k,j)+f*a(k,i)
1148  end do
1149  end do
1150  do j=i,m
1151  a(j,i)=a(j,i)*g
1152  end do
1153  else
1154  do j= i,m
1155  a(j,i)=0.0d0
1156  end do
1157  endif
1158  a(i,i)=a(i,i)+1.0d0
1159  end do
1160  do k=n,1,-1
1161  do its=1,30
1162  do l=k,1,-1
1163  nm=l-1
1164  IF ((abs(rv1(l))+anorm).eq.anorm) goto 2
1165  IF ((abs(w(nm))+anorm).eq.anorm) goto 1
1166  end do
1167 1 c=0.0d0
1168  s=1.0d0
1169  do i=l,k
1170  f=s*rv1(i)
1171  rv1(i)=c*rv1(i)
1172  IF ((abs(f)+anorm).eq.anorm) goto 2
1173  g=w(i)
1174  h=dpythag(f,g)
1175  w(i)=h
1176  h=1.0d0/h
1177  c= (g*h)
1178  s=-(f*h)
1179  do j=1,m
1180  y=a(j,nm)
1181  z=a(j,i)
1182  a(j,nm)=(y*c)+(z*s)
1183  a(j,i)=-(y*s)+(z*c)
1184  end do
1185  end do
1186 2 z=w(k)
1187  IF (l.eq.k)then
1188  IF (z.lt.0.0d0)then
1189  w(k)=-z
1190  do j=1,n
1191  v(j,k)=-v(j,k)
1192  end do
1193  endif
1194  goto 3
1195  endif
1196  IF (its.eq.30) THEN
1197  WRITE(6,*) 'no convergence in svdcmp'
1198  stop
1199  END IF
1200  x=w(l)
1201  nm=k-1
1202  y=w(nm)
1203  g=rv1(nm)
1204  h=rv1(k)
1205  f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0d0*h*y)
1206  g=dpythag(f,1.0d0)
1207  f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x
1208  c=1.0d0
1209  s=1.0d0
1210  do j=l,nm
1211  i=j+1
1212  g=rv1(i)
1213  y=w(i)
1214  h=s*g
1215  g=c*g
1216  z=dpythag(f,h)
1217  rv1(j)=z
1218  c=f/z
1219  s=h/z
1220  f= (x*c)+(g*s)
1221  g=-(x*s)+(g*c)
1222  h=y*s
1223  y=y*c
1224  do jj=1,n
1225  x=v(jj,j)
1226  z=v(jj,i)
1227  v(jj,j)= (x*c)+(z*s)
1228  v(jj,i)=-(x*s)+(z*c)
1229  end do
1230  z=dpythag(f,h)
1231  w(j)=z
1232  IF (z.ne.0.0d0)then
1233  z=1.0d0/z
1234  c=f*z
1235  s=h*z
1236  endif
1237  f= (c*g)+(s*y)
1238  x=-(s*g)+(c*y)
1239  do jj=1,m
1240  y=a(jj,j)
1241  z=a(jj,i)
1242  a(jj,j)= (y*c)+(z*s)
1243  a(jj,i)=-(y*s)+(z*c)
1244  end do
1245  end do
1246  rv1(l)=0.0d0
1247  rv1(k)=f
1248  w(k)=x
1249  end do
1250 3 continue
1251  end do
1252  return
1253  END SUBROUTINE dsvdcmp
1254 
1255  !/ ------------------------------------------------------------------- /
1256  FUNCTION juldayt(id,mm,iyyy)
1257  ! See numerical recipes 2nd ed. The order of month and day have been swapped!
1258  !*********************************************************************
1259  IMPLICIT NONE
1260  INTEGER id,mm,iyyy
1261  INTEGER igreg
1262  INTEGER*4 ja,jm,jy
1263  INTEGER*4 juldayt
1264  igreg=15+31*(10+12*1582)
1265  jy=iyyy
1266  IF (jy.EQ.0) WRITE(6,*) 'There is no zero year !!'
1267  IF (jy.LT.0) jy=jy+1
1268  IF (mm.GT.2) THEN
1269  jm=mm+1
1270  ELSE
1271  jy=jy-1
1272  jm=mm+13
1273  ENDIF
1274  juldayt=int(365.25*jy)+int(30.6001*jm)+id+1720995
1275  IF (id+31*(mm+12*iyyy).GE.igreg) THEN
1276  ja=int(0.01*jy)
1277  juldayt=juldayt+2-ja+int(0.25*ja)
1278  ENDIF
1279  RETURN
1280  END FUNCTION juldayt
1281 
1282  !/ ------------------------------------------------------------------- /
1283  SUBROUTINE caldatt(julian,id,mm,iyyy)
1284  ! See numerical recipes 2nd ed. The order of month and day have been swapped!
1285  ! Should be removed : same is now in W3TIMEMD
1286  !*********************************************************************
1287  IMPLICIT NONE
1288  INTEGER(KIND=4), INTENT(in) :: julian
1289  INTEGER(KIND=4), INTENT(out) :: id,mm,iyyy
1290  INTEGER(KIND=4), PARAMETER :: IGREG=2299161
1291  INTEGER(KIND=4) ja,jalpha,jb,jc,jd,je
1292  if (julian.GE.igreg) THEN
1293  jalpha=int(((julian-1867216)-0.25)/36524.25)
1294  ja=julian+1+jalpha-int(0.25*jalpha)
1295  ELSE
1296  ja=julian
1297  ENDIF
1298  jb=ja+1524
1299  jc=int(6680.+((jb-2439870)-122.1)/365.25)
1300  jd=365*jc+int(0.25*jc)
1301  je=int((jb-jd)/30.6001)
1302  id=jb-jd-int(30.6001*je)
1303  mm=je-1
1304  IF (mm.GT.12) mm=mm-12
1305  iyyy=jc-4715
1306  IF (mm.GT.2) iyyy=iyyy-1
1307  IF (iyyy.LE.0) iyyy=iyyy-1
1308  RETURN
1309  END SUBROUTINE caldatt
1310 
1311  !/ ------------------------------------------------------------------- /
1312  subroutine svd(q,u,v,cov,w,p,b,sig,ic,m,n,mm,N2,toler,jc &
1313  ,ssq,res)
1314  !-----------------------------------------------------------------------
1315  ! svd uses singular-value-decomposition to calculate the least-squares
1316  ! solution p to an overdetermined system of linear equations with
1317  ! coefficient matrix q, which includes right hand side vector b.
1318  !
1319  ! there are two ways to use svd:
1320  ! 1 given an overdetermined system, svd will orthogonalize
1321  ! a and b and produce the least-squares solution.
1322  ! 2 given an orthogonalized a (i.e. output from 1),
1323  ! svd will orthogonalize b with respect to a and produce
1324  ! the least-squares solution. this allows the use of
1325  ! multiple r.h.s. without reorthogonalizing a.
1326  !-----------------------------------------------------------------------
1327  ! description of parameters:
1328  ! ic an input code which must be set to 1 or 2
1329  ! m the number of equations (rows of q) to solve.
1330  ! n the total number of columns of q to be used (<N2)
1331  ! N2 the number of columns of q, at least n+1.
1332  ! mm the number of rows of q, at least N2+m.
1333  ! q an mm-by-N2 array which on entry must contain the
1334  ! matrix a in its first m rows and n columns, and the vector
1335  ! b in its first m rows of column mm. on exit the residual
1336  ! of equation i is stored in q(i,N2), i=1 to m.
1337  ! sig measurement error (standard deviation) for the rhs
1338  ! from the calling program (can be set to 1.)
1339  ! u an mm-by-N2 matrix used in svd, replaced by "left" matrix u
1340  ! v an N2-by-N2 matrix of n "right" eigenvectors of q (=u)
1341  ! w an N2-diagonal vector(matrix) of n eigenvalues of q (u)
1342  ! cov an output covariance matrix between eigenvectors of q (u)
1343  ! toler an input tolerance, for acceptable matrix condition number
1344  ! p an output array of dimension at least N2 containing
1345  ! the least-squares solution in positions 1 to n.
1346  ! jc an output code which is set to the index of the
1347  ! 1st dependent column, if such a column is detected.
1348  ! ssq sum of squares of the residuals
1349  ! res the largest residual in magnitude
1350  !-----------------------------------------------------------------------
1351  ! history:
1352  ! written by j. cherniawsky, august 1997
1353  ! last modified august 1998
1354  !-----------------------------------------------------------------------
1355  IMPLICIT NONE
1356 
1357  ! PARAMETERS
1358  integer, parameter :: nwt=302 ! make nwt > expected value of N2
1359 
1360  ! I/O VARIABLES
1361  integer, intent(IN) :: ic, m, n, N2, mm
1362  real, intent(IN) :: toler
1363  real, intent(INOUT) :: q(mm,N2), ssq, res
1364  integer, intent(INOUT) :: jc
1365  double precision, intent(INOUT) :: sig(mm)
1366  double precision, intent(OUT) :: u(mm,N2),v(N2,N2),cov(N2,N2), &
1367  w(N2),b(mm),p(N2)
1368 
1369  ! LOCAL VARIABLES
1370  double precision :: wti(nwt)
1371  integer :: i, j, k
1372  real :: wmax, thresh
1373  double precision :: eps, sum, resi
1374 
1375  jc=0
1376  do i=1,mm
1377  b(i)=q(i,n2)
1378  enddo
1379  ! no need to solve if only rhs has changed
1380  IF (ic.eq.2) go to 10
1381  ! define a "design matrix" u(=a) and set-up working arrays
1382  do j=1,n2
1383  do i=1,mm
1384  u(i,j)=q(i,j)
1385  enddo
1386  enddo
1387  ! compute svd decomposition of u(=a), with a being replaced by its upper
1388  ! matrix u, viz a=u*w*transpose(v), and vector w is output of a diagonal
1389  ! matrix of singular values w(i), i=1,n.
1390  call dsvdcmp(u,m,n,mm,n2,w,v)
1391  ! check for small singular values
1392  wmax=0.
1393  do j=1,n
1394  IF (w(j).gt.wmax) wmax=w(j)
1395  enddo
1396  thresh=toler*wmax
1397  do j=1,n
1398  IF (w(j).lt.thresh) then
1399  w(j)=0.d0
1400  IF (jc.lt.1) jc=j
1401  endif
1402  enddo
1403 10 eps=1.d-10
1404  ! compute summation weights (wti, used below)
1405  do j=1,n
1406  wti(j)=0.d0
1407  IF (w(j).gt.eps) then
1408  ! wti(j)=sig(j)*sig(j)/(w(j)*w(j))
1409  wti(j)=1.d0/(w(j)*w(j))
1410  endif
1411  enddo
1412  ! use back-substitution to compute the solution p(i), i=1,n
1413  call dsvbksb(u,w,v,m,n,mm,n2,b,p)
1414  ! compute chisq (=ssq) and the largest residual (res)
1415  ssq=0.
1416  res=0.
1417  do i=1,m
1418  sum=0.d0
1419  do j=1,n
1420  sum=sum+p(j)*q(i,j)
1421  enddo
1422  resi=abs(b(i)-sum)
1423  ! TIDE_MF addition
1424  q(i,n2)=b(i)-sum
1425  res=max(res,resi)
1426  ssq=ssq+resi**2
1427  enddo
1428  ! compute variances, covariances, these may need to be given dimension
1429  ! of b(i), e.g., using sig(i), but this is better done after return to main
1430  do i=1,n
1431  do j=1,i
1432  sum=0.d0
1433  do k=1,n
1434  sum=sum+v(i,k)*v(j,k)*wti(k)
1435  enddo
1436  cov(i,j)=sum
1437  cov(j,i)=sum
1438  enddo
1439  enddo
1440  return
1441  end subroutine svd
1442 
1443  !/ ------------------------------------------------------------------- /
1444  SUBROUTINE vuf_set_parameters
1445  !
1446  !***********************************************************************
1447  !* HERE THE MAIN CONSTITUENTS AND THEIR DOODSON NUMBERS ARE SET
1448  !* FORMAT (6X,A5,1X,6I3,F5.2,I4). THE VALUES ARE RESPECTIVELY
1449  !* TIDECON_ALLNAMES = CONSTITUENT NAME
1450  !* II,JJ,KK,LL,MM,NN = THE SIX DOODSON NUMBERS
1451  !* SEMI = PHASE CORRECTION
1452  !* NJ = THE NUMBER OF SATELLITES FOR THIS CONSTITUENT.
1453  !* THE END OF ALL MAIN CONSTITUENTS IS DENOTED BY A BLANK CARD.
1454  !
1455 
1456  IMPLICIT NONE
1457 
1458  INTEGER :: JLM
1459 
1460  ntidal_con = 45
1461  ntotal_con = 45+101
1462  nkonco = 251
1463  jlm = 170
1464 
1465  ALLOCATE(tide_indexj(ntotal_con),tide_indexjk(nkonco))
1466  !
1467  ALLOCATE(tidecon_allnames(ntotal_con))
1468  ALLOCATE(ii(ntidal_con),jj(ntidal_con),kk(ntidal_con),ll(ntidal_con),mm(ntidal_con), &
1469  nn(ntidal_con),semi(ntidal_con), nj(ntotal_con))
1470 
1471  ALLOCATE(konco_con(nkonco),coef_con(nkonco))
1472 
1473  tidecon_allnames(:)=(/ &
1474  'Z0 ','SA ','SSA ','MSM ','MM ','MSF ','MF ','ALP1 ','2Q1 ','SIG1 ', &
1475  'Q1 ','RHO1 ','O1 ','TAU1 ','BET1 ','NO1 ','CHI1 ','PI1 ','P1 ','S1 ', &
1476  'K1 ','PSI1 ','PHI1 ','THE1 ','J1 ','OO1 ','UPS1 ','OQ2 ','EPS2 ','2N2 ', &
1477  'MU2 ','N2 ','NU2 ','GAM2 ','H1 ','M2 ','H2 ','LDA2 ','L2 ','T2 ', &
1478  'S2 ','R2 ','K2 ','ETA2 ','M3 ','2PO1 ','SO1 ','ST36 ','2NS2 ','ST37 ', &
1479  'ST1 ','ST2 ','ST3 ','O2 ','ST4 ','SNK2 ','OP2 ','MKS2 ','ST5 ','ST6 ', &
1480  '2SK2 ','MSN2 ','ST7 ','2SM2 ','ST38 ','SKM2 ','2SN2 ','NO3 ','MO3 ','NK3 ', &
1481  'SO3 ','MK3 ','SP3 ','SK3 ','ST8 ','N4 ','3MS4 ','ST39 ','MN4 ','ST40 ', &
1482  'ST9 ','M4 ','ST10 ','SN4 ','KN4 ','MS4 ','MK4 ','SL4 ','S4 ','SK4 ', &
1483  'MNO5 ','2MO5 ','3MP5 ','MNK5 ','2MP5 ','2MK5 ','MSK5 ','3KM5 ','2SK5 ','ST11 ', &
1484  '2NM6 ','ST12 ','ST41 ','2MN6 ','ST13 ','M6 ','MSN6 ','MKN6 ','2MS6 ','2MK6 ', &
1485  'NSK6 ','2SM6 ','MSK6 ','ST42 ','S6 ','ST14 ','ST15 ','M7 ','ST16 ','3MK7 ', &
1486  'ST17 ','ST18 ','3MN8 ','ST19 ','M8 ','ST20 ','ST21 ','3MS8 ','3MK8 ','ST22 ', &
1487  'ST23 ','ST24 ','ST25 ','ST26 ','4MK9 ','ST27 ','ST28 ','M10 ','ST29 ','ST30 ', &
1488  'ST31 ','ST32 ','ST33 ','M12 ','ST34 ','ST35 '/)
1489 
1490  ii(:)=(/ 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, &
1491  1, 1, 1, 1, 1, 1, 1, 1, 1, 1, &
1492  1, 1, 1, 1, 1, 1, 1, 2, 2, 2, &
1493  2, 2, 2, 2, 2, 2, 2, 2, 2, 2, &
1494  2, 2, 2, 2, 3 /)
1495  jj(:)=(/ 0, 0, 0, 1, 1, 2, 2, -4, -3, -3, &
1496  -2, -2, -1, -1, 0, 0, 0, 1, 1, 1, &
1497  1, 1, 1, 2, 2, 3, 4, -3, -3, -2, &
1498  -2, -1, -1, 0, 0, 0, 0, 1, 1, 2, &
1499  2, 2, 2, 3, 0 /)
1500  kk(:)=(/ 0, 1, 2, -2, 0, -2, 0, 2, 0, 2, &
1501  0, 2, 0, 2, -2, 0, 2, -3, -2, -1, &
1502  0, 1, 2, -2, 0, 0, 0, 0, 2, 0, &
1503  2, 0, 2, -2, -1, 0, 1, -2, 0, -3, &
1504  -2, -1, 0, 0, 0 /)
1505  ll(:)=(/ 0, 0, 0, 1, -1, 0, 0, 1, 2, 0, &
1506  1, -1, 0, 0, 1, 1, -1, 0, 0, 0, &
1507  0, 0, 0, 1, -1, 0, -1, 3, 1, 2, &
1508  0, 1, -1, 2, 0, 0, 0, 1, -1, 0, &
1509  0, 0, 0, -1, 0 /)
1510  mm(:)=(/ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1511  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1512  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1513  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1514  0, 0, 0, 0, 0 /)
1515  nn(:)=(/ 0, -1, 0, 0, 0, 0, 0, 0, 0, 0, &
1516  0, 0, 0, 0, 0, 0, 0, 1, 0, 1, &
1517  0, -1, 0, 0, 0, 0, 0, 0, 0, 0, &
1518  0, 0, 0, 0, 1, 0, -1, 0, 0, 1, &
1519  0, -1, 0, 0, 0 /)
1520  semi(:)=(/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00,-0.25,-0.25,-0.25, &
1521  -0.25,-0.25,-0.25,-0.75,-0.75,-0.75,-0.75,-0.25,-0.25,-0.75, &
1522  -0.75,-0.75,-0.75,-0.75,-0.75,-0.75,-0.75, 0.00, 0.00, 0.00, &
1523  0.00, 0.00, 0.00,-0.50,-0.50, 0.00, 0.00,-0.50,-0.50, 0.00, &
1524  0.00,-0.50, 0.00, 0.00,-0.50 /)
1525  nj(:)=(/ 1, 1, 1, 1, 1, 1, 1, 2, 5, 4, &
1526  10, 5, 8, 5, 1, 9, 2, 1, 6, 2, &
1527  10, 1, 5, 4, 10, 8, 5, 2, 3, 4, &
1528  3, 4, 4, 3, 2, 9, 1, 1, 5, 1, &
1529  3, 2, 5, 7, 1, 2, 2, 3, 2, 2, &
1530  3, 4, 3, 1, 3, 3, 2, 3, 3, 4, &
1531  2, 3, 4, 2, 3, 3, 2, 2, 2, 2, &
1532  2, 2, 2, 2, 3, 1, 2, 4, 2, 3, &
1533  4, 1, 3, 2, 2, 2, 2, 2, 1, 2, &
1534  3, 2, 2, 3, 2, 2, 3, 3, 2, 3, &
1535  2, 4, 3, 2, 4, 1, 3, 3, 2, 2, &
1536  3, 2, 3, 3, 1, 3, 3, 1, 3, 2, &
1537  4, 2, 2, 4, 1, 3, 3, 2, 2, 4, &
1538  2, 3, 3, 3, 2, 3, 2, 1, 3, 2, &
1539  4, 2, 3, 1, 2, 4/)
1540 
1541  ldel(1:jlm)=(/ 0, 0, 0, 0, 0, 0, 0, -1, 0, -2, &
1542  -1, -1, 0, 0, -1, 0, 0, 2, -2, -2, &
1543  -1, -1, -1, 0, -1, 0, 1, 2, 0, 0, &
1544  1, 2, 2, -1, 0, 0, 1, 1, 1, 2, &
1545  2, -2, -1, 0, 0, 0, 0, -2, -2, -2, &
1546  -1, -1, -1, 0, 0, 0, 0, 0, 0, 0, &
1547  0, 0, 1, 2, 2, 0, 0, -2, -1, -1, &
1548  -1, 0, 0, 0, 0, 1, 1, 0, -2, -2, &
1549  0, 0, 0, -2, -1, 0, 0, 0, 0, 0, &
1550  1, 1, 1, 1, 2, 2, 2, -2, -2, -2, &
1551  -1, -1, 0, 0, 0, -2, 0, 0, 1, 1, &
1552  -1, 0, -1, -1, 0, -2, -1, -1, 0, -1, &
1553  -1, 0, -2, -1, 0, 0, 0, 1, 2, 2, &
1554  -2, -1, 0, 0, 1, -1, -1, 0, 0, 1, &
1555  1, 1, 2, 2, 0, 0, 0, 2, 2, 2, &
1556  2, 0, 0, 1, 2, 0, 0, -1, -1, 0, &
1557  0, 0, 0, 0, 0, 1, 1, 1, 2, 0/)
1558  mdel(1:jlm)=(/ 0, 0, 0, 0, 0, 0, 0, 0, -1, -2, &
1559  -1, 0, -2, -1, 0, -2, -1, 0, -3, -2, &
1560  -2, -1, 0, -2, 0, -1, 0, 0, -2, -1, &
1561  0, 0, 1, 0, -2, -1, -1, 0, 1, 0, &
1562  1, 0, 0, -1, 1, 2, -1, -2, -1, 0, &
1563  -1, 0, 1, -1, 1, 2, -1, 1, -1, -2, &
1564  -1, 0, 0, 0, 1, 0, 1, -1, -1, 0, &
1565  1, -2, -1, 1, 2, 0, 1, 1, 0, 1, &
1566  0, 1, 2, -1, 0, -1, 1, -1, 1, 2, &
1567  -1, 0, 1, 2, 0, 1, 2, -1, 0, 1, &
1568  0, 1, 1, 2, 3, 0, 1, 2, 0, 1, &
1569  0, -1, -1, 0, -1, -2, -1, 0, -1, -1, &
1570  0, -1, -2, 0, -2, -1, -1, 0, 0, 1, &
1571  -2, 0, -1, -1, 0, -1, 0, -2, -1, -1, &
1572  0, 1, 0, 1, -1, -1, -1, -1, 0, 1, &
1573  2, 0, -1, 0, 0, 0, 1, 0, 1, -1, &
1574  1, 2, -1, 1, 2, 0, 1, 2, 0, -1 /)
1575  ndel(1:jlm)=(/ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1576  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1577  0, 0, 0, 0, 1, 0, 0, 0, 0, 0, &
1578  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1579  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1580  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1581  0, 2, 0, 0, 0, -2, 0, 0, 0, 0, &
1582  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1583  -2, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1584  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1585  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1586  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1587  0, 0, 0, 1, 0, 0, 0, 0, 0, 0, &
1588  0, 0, 0, 0, -1, 0, 0, 0, 0, 0, &
1589  0, 0, 0, 0, 0, 0, 0, 0, 0, 0, &
1590  0, 0, 0, 0, 0, 2, 2, 0, 0, 0, &
1591  0, 0, 0, 0, 0, 0, 0, 0, 0, 0 /)
1592  ph(1:jlm)=(/ 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.75, 0.00, 0.50, &
1593  0.75, 0.75, 0.50, 0.00, 0.75, 0.50, 0.00, 0.50, 0.50, 0.50, &
1594  0.75, 0.75, 0.75, 0.50, 0.00, 0.00, 0.75, 0.50, 0.50, 0.00, &
1595  0.75, 0.50, 0.00, 0.25, 0.50, 0.00, 0.25, 0.75, 0.25, 0.50, &
1596  0.50, 0.00, 0.25, 0.50, 0.50, 0.50, 0.00, 0.50, 0.00, 0.00, &
1597  0.75, 0.25, 0.75, 0.50, 0.00, 0.50, 0.50, 0.00, 0.50, 0.00, &
1598  0.50, 0.50, 0.75, 0.50, 0.50, 0.00, 0.50, 0.00, 0.75, 0.25, &
1599  0.75, 0.00, 0.50, 0.00, 0.50, 0.25, 0.25, 0.00, 0.00, 0.00, &
1600  0.00, 0.50, 0.50, 0.00, 0.25, 0.50, 0.00, 0.50, 0.00, 0.50, &
1601  0.75, 0.25, 0.25, 0.25, 0.50, 0.50, 0.50, 0.50, 0.00, 0.00, &
1602  0.25, 0.25, 0.00, 0.00, 0.00, 0.00, 0.00, 0.00, 0.25, 0.25, &
1603  0.25, 0.50, 0.25, 0.25, 0.50, 0.50, 0.25, 0.25, 0.50, 0.25, &
1604  0.25, 0.50, 0.50, 0.00, 0.00, 0.50, 0.50, 0.75, 0.00, 0.50, &
1605  0.00, 0.25, 0.50, 0.50, 0.50, 0.75, 0.75, 0.00, 0.50, 0.25, &
1606  0.75, 0.75, 0.00, 0.00, 0.50, 0.50, 0.50, 0.00, 0.50, 0.50, &
1607  0.50, 0.00, 0.00, 0.75, 0.00, 0.50, 0.00, 0.75, 0.75, 0.50, &
1608  0.00, 0.00, 0.50, 0.00, 0.00, 0.75, 0.75, 0.75, 0.50, 0.50 /)
1609 
1610  ee(1:jlm)=(/ 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0000, 0.0360, 0.1906, 0.0063, &
1611  0.0241, 0.0607, 0.0063, 0.1885, 0.0095, 0.0061, 0.1884, 0.0087, 0.0007, 0.0039, &
1612  0.0010, 0.0115, 0.0292, 0.0057, 0.0008, 0.1884, 0.0018, 0.0028, 0.0058, 0.1882, &
1613  0.0131, 0.0576, 0.0175, 0.0003, 0.0058, 0.1885, 0.0004, 0.0029, 0.0004, 0.0064, &
1614  0.0010, 0.0446, 0.0426, 0.0284, 0.2170, 0.0142, 0.2266, 0.0057, 0.0665, 0.3596, &
1615  0.0331, 0.2227, 0.0290, 0.0290, 0.2004, 0.0054, 0.0282, 0.2187, 0.0078, 0.0008, &
1616  0.0112, 0.0004, 0.0004, 0.0015, 0.0003, 0.3534, 0.0264, 0.0002, 0.0001, 0.0007, &
1617  0.0001, 0.0001, 0.0198, 0.1356, 0.0029, 0.0002, 0.0001, 0.0190, 0.0344, 0.0106, &
1618  0.0132, 0.0384, 0.0185, 0.0300, 0.0141, 0.0317, 0.1993, 0.0294, 0.1980, 0.0047, &
1619  0.0027, 0.0816, 0.0331, 0.0027, 0.0152, 0.0098, 0.0057, 0.0037, 0.1496, 0.0296, &
1620  0.0240, 0.0099, 0.6398, 0.1342, 0.0086, 0.0611, 0.6399, 0.1318, 0.0289, 0.0257, &
1621  0.1042, 0.0386, 0.0075, 0.0402, 0.0373, 0.0061, 0.0117, 0.0678, 0.0374, 0.0018, &
1622  0.0104, 0.0375, 0.0039, 0.0008, 0.0005, 0.0373, 0.0373, 0.0042, 0.0042, 0.0036, &
1623  0.1429, 0.0293, 0.0330, 0.0224, 0.0447, 0.0001, 0.0004, 0.0005, 0.0373, 0.0001, &
1624  0.0009, 0.0002, 0.0006, 0.0002, 0.0217, 0.0448, 0.0366, 0.0047, 0.2505, 0.1102, &
1625  0.0156, 0.0000, 0.0022, 0.0001, 0.0001, 0.2535, 0.0141, 0.0024, 0.0004, 0.0128, &
1626  0.2980, 0.0324, 0.0187, 0.4355, 0.0467, 0.0747, 0.0482, 0.0093, 0.0078, 0.0564 /)
1627 
1628  ir(1:jlm)=(/ 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, &
1629  1, 1, 0, 0, 1, 0, 0, 0, 0, 0, &
1630  1, 1, 1, 0, 0, 0, 1, 0, 0, 0, &
1631  1, 0, 0, 1, 0, 0, 1, 1, 1, 0, &
1632  0, 0, 1, 0, 0, 0, 0, 0, 0, 0, &
1633  1, 1, 1, 0, 0, 0, 0, 0, 0, 0, &
1634  0, 0, 1, 0, 0, 0, 0, 0, 1, 1, &
1635  1, 0, 0, 0, 0, 1, 1, 0, 0, 0, &
1636  0, 0, 0, 0, 1, 0, 0, 0, 0, 0, &
1637  1, 1, 1, 1, 0, 0, 0, 0, 0, 0, &
1638  1, 1, 0, 0, 0, 0, 0, 0, 1, 1, &
1639  2, 0, 2, 2, 0, 0, 2, 2, 0, 2, &
1640  2, 0, 0, 0, 0, 0, 0, 2, 0, 0, &
1641  0, 2, 0, 0, 0, 2, 2, 0, 0, 2, &
1642  2, 2, 0, 0, 0, 0, 0, 0, 0, 0, &
1643  0, 0, 0, 2, 0, 0, 0, 2, 2, 0, &
1644  0, 0, 0, 0, 0, 2, 2, 2, 0, 0 /)
1645 
1646  coef_con(:)=(/ 2.00,-1.00, 1.00,-1.00, 2.00, 1.00,-2.00, 2.00,-1.00, 3.00, &
1647  -2.00, 2.00, 1.00,-2.00, 1.00, 1.00, 1.00,-2.00, 2.00, 1.00, &
1648  -2.00, 2.00, 2.00, 1.00,-2.00, 1.00, 1.00,-1.00, 1.00, 1.00, &
1649  1.00, 1.00,-1.00, 1.00, 2.00,-2.00, 2.00, 1.00,-1.00,-1.00, &
1650  2.00,-1.00, 1.00, 1.00,-1.00, 2.00, 1.00,-1.00,-1.00, 2.00, &
1651  -1.00, 2.00, 1.00,-2.00, 1.00, 1.00,-1.00, 2.00,-1.00, 1.00, &
1652  1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, &
1653  1.00, 1.00, 1.00, 2.00, 1.00,-1.00, 2.00, 3.00,-1.00, 1.00, &
1654  1.00, 1.00,-1.00, 1.00, 1.00, 2.00, 1.00,-1.00, 1.00, 1.00, &
1655  1.00,-1.00, 2.00, 2.00, 1.00,-1.00, 1.00, 1.00, 1.00, 1.00, &
1656  1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 2.00, 1.00, 1.00, 1.00, &
1657  1.00, 1.00, 2.00, 1.00, 3.00,-1.00, 1.00, 1.00, 1.00, 2.00, &
1658  1.00, 2.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 2.00, &
1659  1.00, 3.00, 1.00,-1.00, 2.00, 1.00, 2.00, 1.00, 1.00,-1.00, &
1660  3.00, 1.00,-1.00, 2.00, 1.00, 2.00, 1.00, 1.00,-1.00, 3.00, &
1661  1.00, 1.00, 1.00, 1.00, 1.00, 1.00, 2.00, 1.00, 2.00, 1.00, &
1662  1.00, 1.00, 1.00, 2.00, 1.00, 1.00, 1.00, 1.00, 2.00, 2.00, &
1663  -1.00, 3.00, 2.00, 1.00, 1.00, 2.00, 1.00, 1.00, 3.50, 2.00, &
1664  1.00, 1.00, 3.00, 1.00, 1.00, 1.00, 1.00, 1.00, 2.00, 2.00, &
1665  3.00, 1.00, 3.00, 1.00, 1.00,-1.00, 4.00, 2.00, 1.00, 1.00, &
1666  2.00, 1.00, 1.00, 3.00, 1.00, 3.00, 1.00, 1.00, 1.00, 1.00, &
1667  1.00, 2.00, 2.00, 2.00, 1.00, 1.00, 2.00, 2.00, 1.00, 3.00, &
1668  1.00, 1.00, 4.00, 1.00, 3.00, 1.00, 1.00, 4.00, 1.00, 5.00, &
1669  3.00, 1.00, 1.00, 4.00, 1.00, 2.00, 1.00, 1.00, 1.00, 3.00, &
1670  2.00, 4.00, 1.00, 1.00, 6.00, 5.00, 1.00, 3.00, 1.00, 1.00, &
1671  1.00 /)
1672  konco_con(:)=(/ &
1673  'P1 ','O1 ','S2 ','O1 ','M2 ','N2 ','S2 ','N2 ','S2 ','M2 ', &
1674  'S2 ','N2 ','K2 ','S2 ','M2 ','N2 ','K2 ','S2 ','M2 ','S2 ', &
1675  'K2 ','O1 ','K2 ','N2 ','S2 ','S2 ','N2 ','K2 ','O1 ','P1 ', &
1676  'M2 ','K2 ','S2 ','M2 ','K2 ','S2 ','S2 ','N2 ','M2 ','K2 ', &
1677  'S2 ','K2 ','M2 ','S2 ','N2 ','K2 ','M2 ','S2 ','N2 ','S2 ', &
1678  'M2 ','M2 ','S2 ','N2 ','S2 ','K2 ','M2 ','S2 ','N2 ','N2 ', &
1679  'O1 ','M2 ','O1 ','N2 ','K1 ','S2 ','O1 ','M2 ','K1 ','S2 ', &
1680  'P1 ','S2 ','K1 ','M2 ','N2 ','S2 ','N2 ','M2 ','S2 ','M2 ', &
1681  'S2 ','N2 ','K2 ','M2 ','N2 ','M2 ','S2 ','K2 ','M2 ','N2 ', &
1682  'K2 ','S2 ','M2 ','M2 ','K2 ','S2 ','S2 ','N2 ','K2 ','N2 ', &
1683  'M2 ','S2 ','M2 ','K2 ','S2 ','L2 ','S2 ','S2 ','K2 ','M2 ', &
1684  'N2 ','O1 ','M2 ','O1 ','M2 ','P1 ','M2 ','N2 ','K1 ','M2 ', &
1685  'P1 ','M2 ','K1 ','M2 ','S2 ','K1 ','K2 ','K1 ','M2 ','S2 ', &
1686  'K1 ','N2 ','K2 ','S2 ','N2 ','M2 ','N2 ','M2 ','K2 ','S2 ', &
1687  'M2 ','S2 ','K2 ','M2 ','N2 ','M2 ','N2 ','K2 ','S2 ','M2 ', &
1688  'M2 ','S2 ','N2 ','M2 ','K2 ','N2 ','M2 ','S2 ','M2 ','K2 ', &
1689  'N2 ','S2 ','K2 ','S2 ','M2 ','M2 ','S2 ','K2 ','M2 ','S2 ', &
1690  'K2 ','S2 ','M2 ','N2 ','O1 ','N2 ','M2 ','K1 ','M2 ','M2 ', &
1691  'S2 ','O1 ','M2 ','K1 ','M2 ','S2 ','K2 ','O1 ','M2 ','N2 ', &
1692  'M2 ','N2 ','M2 ','N2 ','K2 ','S2 ','M2 ','M2 ','S2 ','N2 ', &
1693  'M2 ','N2 ','K2 ','M2 ','S2 ','M2 ','K2 ','M2 ','S2 ','N2 ', &
1694  'K2 ','M2 ','S2 ','M2 ','S2 ','K2 ','M2 ','N2 ','K1 ','M2 ', &
1695  'N2 ','K1 ','M2 ','K1 ','M2 ','S2 ','K1 ','M2 ','N2 ','M2 ', &
1696  'M2 ','N2 ','S2 ','M2 ','S2 ','M2 ','N2 ','S2 ','K2 ','M2 ', &
1697  'S2 ','M2 ','S2 ','K1 ','M2 ','M2 ','S2 ','M2 ','N2 ','K2 ', &
1698  'S2 '/)
1699 
1700 
1701  END SUBROUTINE vuf_set_parameters
1702 
1703  !/ ------------------------------------------------------------------- /
1704  SUBROUTINE vuf (KONX,Vx,ux,FX,ITIME)
1705  !
1706  !
1707  !* GIVEN CONSTITUENT KONX , THE NODAL CORRECTIONS V+U AND F ARE RETURNED
1708  ! That subroutine can now be replaced by look-up table to get K from J
1709  !
1710  !
1711  IMPLICIT NONE
1712  !
1713  character*5, INTENT(IN) :: KONX
1714  REAL(KIND=8), intent(out) :: vx,ux,fx
1715  INTEGER, INTENT(IN) :: ITIME
1716 
1717  INTEGER :: K
1718 
1719  DO k=1,ntotal_con
1720  IF (tidecon_allnames(k).eq.konx) go to 40
1721  END DO
1722  WRITE(ndset,30) konx
1723 30 FORMAT('ERROR IN VUF: STOP.',a5)
1724  stop
1725 40 vx=v_arg(k,itime)
1726  ux=u_arg(k,itime)
1727  fx=f_arg(k,itime)
1728  RETURN
1729  !
1730  !***********************************************************************
1731  !* THE ASTRONOMICAL ARGUMENTS AND THEIR RATES OF CHANGE,
1732  !* S0,H0,P0,ENP0,PP0,DS,DH,DP,DNP,DPP, ARE READ FROM TWO RECORDS IN
1733  !* THE FORMAT(5F13.10):
1734  !* S0 = MEAN LONGITUDE OF THE MOON (CYCLES) AT 000 ET 1/1/1976.
1735  !* H0 = MEAN LONGITUDE OF THE SUN.
1736  !* P0 = MEAN LONGITUDE OF THE LUNAR PERIGEE.
1737  !* ENP0= NEGATIVE OF THE MEAN LONGITUDE OF THE ASCENDING NODE.
1738  !* PP0 = MEAN LONGITUDE OF THE SOLAR PERIGEE (PERIHELION).
1739  !* DS,DH,DP,DNP,DPP ARE THEIR RESPECTIVE RATES OF CHANGE OVER A 365
1740  !* DAY PERIOD AS OF 000 ET 1/1/1976.
1741  !
1742 
1743  END SUBROUTINE vuf
1744 
1745  !/ ------------------------------------------------------------------- /
1746  SUBROUTINE opnvuf(filename)
1747 
1748  IMPLICIT NONE
1749 
1750  CHARACTER*256, INTENT(IN) :: filename
1751 
1752  INTEGER :: J, J1, JBASE, J4, K, K1, JL, JLM, KR
1753 
1754  DOUBLE PRECISION, PARAMETER :: TWOPI=3.1415926535898*2.
1755  DOUBLE PRECISION, PARAMETER :: FAC=twpi/360.
1756 
1757  kr=8
1758  !
1759  !***********************************************************************
1760  !* HERE THE MAIN CONSTITUENTS AND THEIR DOODSON NUMBERS ARE READ IN
1761  !* FORMAT (6X,A5,1X,6I3,F5.2,I4). THE VALUES ARE RESPECTIVELY
1762  !* TIDECON_ALLNAMES = CONSTITUENT NAME
1763  !* II,JJ,KK,LL,MM,NN = THE SIX DOODSON NUMBERS
1764  !* SEMI = PHASE CORRECTION
1765  !* NJ = THE NUMBER OF SATELLITES FOR THIS CONSTITUENT.
1766  !* THE END OF ALL MAIN CONSTITUENTS IS DENOTED BY A BLANK CARD.
1767  !
1768  ALLOCATE(tidecon_allnames(170))
1769  ALLOCATE(ii(mc2),jj(mc2),kk(mc2),ll(mc2),mm(mc2),nn(mc2), &
1770  semi(mc2), nj(170))
1771  ALLOCATE(konco_con(320),coef_con(320))
1772 
1773  kr=8
1774  open(unit=kr,file=filename,status='old',form='formatted')
1775 
1776  jbase=0
1777  DO k=1,1000
1778  READ(kr,60)tidecon_allnames(k),ii(k),jj(k),kk(k),ll(k),mm(k),nn(k),semi(k), &
1779  nj(k)
1780 60 FORMAT(6x,a5,1x,6i3,f5.2,i4)
1781  !WRITE(995,'(I4,A5,1X,6I3,F5.2,I4)') K,TIDECON_ALLNAMES(K),II(K),JJ(K),KK(K),LL(K),MM(K),NN(K),SEMI(K), &
1782  ! NJ(K)
1783  IF (tidecon_allnames(k).eq.kblank) go to 100
1784 70 j1=jbase+1
1785  IF (nj(k).LT.1) THEN
1786  nj(k)=1
1787  jl=j1
1788  ph(j1)=0.
1789  ee(j1)=0.
1790  ldel(j1)=0
1791  mdel(j1)=0
1792  ndel(j1)=0
1793  ir(j1)=0
1794  ELSE
1795  jl=jbase+nj(k)
1796  !
1797  !***********************************************************************
1798  !* IF NJ>0, INFORMATION ON THE SATELLITE CONSTITUENTS IS READ , THREE
1799  !* SATELLITES PER CARD, IN THE FORMAT (11X,3(3I3,F4.2,F7.4,1X,I1,1X)).
1800  !* FOR EACH SATELLITE THE VALUES READ ARE
1801  !* LDEL,MDEL,NDEL = THE CHANGES IN THE LAST THREE DOODSON NUMBERS
1802  !* FROM THOSE OF THE MAIN CONSTITUENT.
1803  !* PH = THE PHASE CORRECTION
1804  !* EE = THE AMPLITUDE RATIO OF THE SATELLITE TIDAL POTENTIAL TO
1805  !* THAT OF THE MAIN CONSTITUENT.
1806  !* IR = 1 IF THE AMPLITUDE RATIO HAS TO BE MULTIPLIED BY THE
1807  !* LATITUDE CORRECTION FACTOR FOR DIURNAL CONSTITUENTS
1808  !* 2 IF THE AMPLITUDE RATIO HAS TO BE MULTIPLIED BY THE
1809  !* LATITUDE CORRECTION FACTOR FOR SEMI-DIURNAL CONSTI-
1810  !* TUENTS.
1811  !* OTHERWISE IF NO CORRECTION IS REQUIRED TO THE AMPLITUDE
1812  !* RATIO.
1813  !
1814  READ(kr,80)(ldel(j),mdel(j),ndel(j),ph(j),ee(j),ir(j),j=j1,jl)
1815 80 FORMAT((11x,3(3i3,f4.2,f7.4,1x,i1,1x)))
1816  END IF
1817  jbase=jl
1818  end do
1819 100 ntidal_con=k-1
1820  jlm=jl
1821 
1822  !
1823  !***********************************************************************
1824  !* THE SHALLOW WATER CONSTITUENTS AND THE MAIN CONSTITUENTS FROM WHICH
1825  !* THEY ARE DERIVED ARE READ IN HERE WITH THE FORMAT
1826  !* (6X,A5,I1,2X,4(F5.2,A5,5X)). THE VALUES ARE RESPECTIVELY
1827  !* TIDECON_ALLNAMES = NAME OF THE SHALLOW WATER CONSTITUENT
1828  !* NJ = NUMBER OF MAIN CONSTITUENTS FROM WHICH IT IS DERIVED.
1829  !* COEF_CON,KONCO_CON = COMBINATION NUMBER AND NAME OF THESE MAIN
1830  !* CONSTITUENTS.
1831  !* THE END OF THESE CONSTITUENTS IS DENOTED BY A BLANK CARD.
1832  !
1833  jbase=0
1834  k1=ntidal_con+1
1835  DO k=k1,1000
1836  j1=jbase+1
1837  j4=j1+3
1838  READ(kr,130)tidecon_allnames(k),nj(k),(coef_con(j),konco_con(j),j=j1,j4)
1839 130 FORMAT(6x,a5,i1,2x,4(f5.2,a5,5x))
1840  !WRITE(995,130)TIDECON_ALLNAMES(K),NJ(K),(COEF_CON(J),KONCO_CON(J),J=J1,J4)
1841  IF (tidecon_allnames(k).eq.kblank) go to 170
1842  jbase=jbase+nj(k)
1843  end do
1844 
1845 170 ntotal_con=k-1
1846 
1847  ! Write out for cut and paste ...
1848  ! WRITE(6,*) 'Numbers:',NTIDAL_CON, NTOTAL_CON, JLM, J1, J4
1849  ! WRITE(996,*) NTIDAL, NTOTAL_CON, JLM, J1, J4
1850  !999 FORMAT(10("'",A5,"',"),' &')
1851  ! WRITE(996,999) TIDECON_ALLNAMES(1:NTOTAL_CON)
1852  !998 FORMAT(10(I3,","),' &')
1853  !997 FORMAT(10(I4,","),' &')
1854  !991 FORMAT(10(F5.2,","),' &')
1855  !990 FORMAT(10(F7.4,","),' &')
1856  ! WRITE(996,998) II(1:NTIDAL_CON)
1857  ! WRITE(996,998) JJ(1:NTIDAL_CON)
1858  ! WRITE(996,998) KK(1:NTIDAL_CON)
1859  ! WRITE(996,998) LL(1:NTIDAL_CON)
1860  ! WRITE(996,998) MM(1:NTIDAL_CON)
1861  ! WRITE(996,998) NN(1:NTIDAL_CON)
1862  ! WRITE(996,991) SEMI(1:NTIDAL_CON)
1863  ! WRITE(996,997) NJ(1:NTOTAL_CON)
1864 
1865  ! WRITE(996,998) LDEL(1:JLM)
1866  ! WRITE(996,998) MDEL(1:JLM)
1867  ! WRITE(996,998) NDEL(1:JLM)
1868  ! WRITE(996,991) PH(1:JLM)
1869  ! WRITE(996,990) EE(1:JLM)
1870  ! WRITE(996,998) IR(1:JLM)
1871 
1872  ! WRITE(996,991) COEF_CON(1:J4)
1873  ! WRITE(996,999) KONCO_CON(1:J4)
1874  RETURN
1875  !
1876  !***********************************************************************
1877  !* NTIDAL_CON IS THE NUMBER OF MAIN CONSTITUENTS
1878  !* NTOTAL_CON IS THE NUMBER OF CONSTITUENTS (MAIN + SHALLOW WATER)
1879  !* FOR THE GIVEN TIME hr, THE TABLE OF F AND V+U VALUES IS
1880  !* CALCULATED FOR ALL THE CONSTITUENTS.
1881  !* F IS THE NODAL MODULATION ADJUSTMENT FACTOR FOR AMPLITUDE
1882  !* U IS THE NODAL MODULATION ADJUSTMENT FACTOR FOR PHASE
1883  !* V IS THE ASTRONOMICAL ARGUMENT ADJUSTMENT FOR PHASE.
1884  !
1885  ! setvuf calculates the V,u,f values at time hr for all constituents
1886  !
1887  END SUBROUTINE opnvuf
1888 
1889  !/ ------------------------------------------------------------------- /
1890  SUBROUTINE setvuf(hr,XLAT,ITIME)
1891  ! setvuf calculates the V,u,f values at time hr for all constituents
1892 
1893  IMPLICIT NONE
1894 
1895  REAL(KIND=8), intent(in) :: hr
1896  REAL, INTENT(IN) :: XLAT
1897  INTEGER, INTENT(IN) :: ITIME
1898 
1899  !
1900  ! Local variables
1901  !
1902  INTEGER :: KD0, INT24, INTDYS, &
1903  JBASE, J, K, L, J1, K1, JL, LK, iflag
1904  INTEGER :: IV, IUU
1905  !
1906  REAL(KIND=4), parameter :: pi=3.1415926536
1907  REAL(KIND=4), parameter :: twopi=2.*3.1415926536
1908  !
1909  REAL :: SLAT, VDBL, VV, SUMC, SUMS, RR, &
1910  UUDBL, UU, CXLAT
1911  REAL(KIND=8) :: d1,h,pp,s,p,enp,dh,dpp,ds,dp,dnp,hh,tau
1912 
1913  INTEGER :: indx(170)
1914 
1915  cxlat = max(abs(xlat), 5.)
1916  slat=sin(pi*cxlat/180.)
1917  !
1918  !***********************************************************************
1919  !* THE ASTRONOMICAL ARGUMENTS ARE CALCULATED BY LINEAR APPROXIMATION
1920  !* AT THE MID POINT OF THE ANALYSIS PERIOD.
1921  !
1922  ! day number measured from January 0.5 1900 (i.e.,
1923  ! 1200 UT December 31, 1899
1924  d1=hr/24.d0
1925  ! This was with "gregorian days from KDAY"
1926  !call gday(31,12,99,18,kd0) ! CALL GDAY(IDd,IMm,IYy,ICc,KDd)
1927  ! Now uses "julian days from JULDAYT"
1928  ! KD0= 693961
1929  !KD0=JULDAYT(31,12,1899) ! JULDAYT(id,mm,iyyy)
1930  kd0= 2415020
1931  ! substracting 0.5day is not necessary anymore with new time functions
1932  d1=d1-dfloat(kd0)
1933  call astr(d1,h,pp,s,p,enp,dh,dpp,ds,dp,dnp)
1934  int24=24
1935  intdys=int((hr+0.00001)/int24)
1936  hh=hr-dfloat(intdys*int24)
1937  tau=hh/24.d0+h-s
1938  !
1939  !***********************************************************************
1940  !* ONLY THE FRACTIONAL PART OF A SOLAR DAY NEED BE RETAINED FOR COMPU-
1941  !* TING THE LUNAR TIME TAU.
1942  !
1943  jbase=0
1944  DO k=1,ntidal_con
1945  do l=1,tide_mf
1946  IF (tidecon_allnames(k).eq.tidecon_name(l)) then
1947  indx(k)=l
1948  END IF
1949  end do
1950  vdbl=ii(k)*tau+jj(k)*s+kk(k)*h+ll(k)*p+mm(k)*enp+nn(k)*pp+semi(k)
1951  iv=vdbl
1952  iv=(iv/2)*2
1953  vv=vdbl-iv
1954  j1=jbase+1
1955  jl=jbase+nj(k)
1956  sumc=1.
1957  sums=0.
1958  DO j=j1,jl
1959  !
1960  !***********************************************************************
1961  !* HERE THE SATELLITE AMPLITUDE RATIO ADJUSTMENT FOR LATITUDE IS MADE
1962  !
1963  rr=ee(j)
1964  l=ir(j)+1
1965  IF (l.EQ.2) THEN
1966  rr=ee(j)*0.36309*(1.-5.*slat*slat)/slat
1967  ELSE IF (l.EQ.3) THEN
1968  rr=ee(j)*2.59808*slat
1969  END IF
1970  uudbl=ldel(j)*p+mdel(j)*enp+ndel(j)*pp+ph(j)
1971  iuu=uudbl
1972  uu=uudbl-iuu
1973  sumc=sumc+rr*cos(uu*twopi)
1974  sums=sums+rr*sin(uu*twopi)
1975  end do
1976  f_arg(k,itime)=sqrt(sumc*sumc+sums*sums)
1977  v_arg(k,itime)=vv
1978  u_arg(k,itime)=atan2(sums,sumc)/twopi
1979  jbase=jl
1980  end do
1981  !
1982  !***********************************************************************
1983  !* HERE F AND V+U OF THE SHALLOW WATER CONSTITUENTS ARE COMPUTED FROM
1984  !* THE VALUES OF THE MAIN CONSTITUENT FROM WHICH THEY ARE DERIVED.
1985  !
1986  jbase=0
1987  k1=ntidal_con+1
1988  IF (k1.GT.ntotal_con) RETURN
1989  !
1990  DO k=k1,ntotal_con
1991  f_arg(k,itime)=1.0
1992  v_arg(k,itime)=0.0
1993  u_arg(k,itime)=0.
1994  iflag=0
1995  DO lk=1,tide_mf
1996  IF (tidecon_allnames(k).eq.tidecon_name(lk)) then
1997  iflag=1
1998  EXIT
1999  END IF
2000  END DO ! lk
2001 
2002  DO j=tide_indexj(k),tide_indexj(k)+nj(k)-1
2003  l=tide_indexjk(j)
2004  f_arg(k,itime)=f_arg(k,itime)*f_arg(l,itime)**abs(coef_con(j))
2005  v_arg(k,itime)=v_arg(k,itime)+coef_con(j)*v_arg(l,itime)
2006  u_arg(k,itime)=u_arg(k,itime)+coef_con(j)*u_arg(l,itime)
2007  END DO ! J
2008  END DO ! K
2009 
2010  ! Test output for verification purposes
2011  IF (itime.EQ.-1) THEN
2012  WRITE(992,'(A,F20.2,13F8.3)') 'TEST ISEA 0:', &
2013  d1,h,s,tau,pp,s,p,enp,dh,dpp,ds,dp,dnp,xlat
2014  do l=1,tide_mf
2015  do k=1,ntotal_con
2016  IF (tidecon_allnames(k).eq.tidecon_name(l)) then
2017  tide_index(l)=k
2018  WRITE(992,'(A,4I9,F12.0,3F8.3,I4,X,A)') 'TEST ISEA 1:',1,l,20071201,0,hr, &
2019  f_arg(k,itime),u_arg(k,itime),v_arg(k,itime),tide_index(l),tidecon_name(l)
2020  END IF
2021  END DO
2022  ENDDO
2023  ENDIF
2024 
2025  RETURN
2026  !
2027  END SUBROUTINE setvuf
2028 
2029 
2030  !/ ------------------------------------------------------------------- /
2031  subroutine flex_tidana_webpage(IX,IY,XLON,XLAT,KD1,KD2,ndef, itrend, RES, SSQ, RMSR0, SDEV0, &
2032  RMSR, RESMAX, IMAX, ITEST)
2033  !
2034  !***********************************************************************
2035  !*
2036  !* THIS PROGRAM DOES A TIDAL HEIGHTS 'HARMONIC' ANALYSIS OF IRREGULARLY
2037  !* SAMPLED OBSERVATIONS. THE ANALYSIS METHOD IS A LEAST SQUARES FIT
2038  !* USING SVD COUPLED WITH NODAL MODULATION AND INFERENCE(IF SO REQUESTED).
2039  !*
2040  !* The code is based on TOPEX analysis code originally developed by Josef
2041  !* Cherniawsky (JAOT, 2001, 18(4): 649-664) and modified by Rob Bell and
2042  !* Mike Foreman. Enhancements to that version include
2043  !*
2044  !* 1. Provision for multi-constituent inferences computed directly within
2045  !* the least squares matrix rather than as post fit corrections. This
2046  !* means that the inferred constituents will affect all constituents, not
2047  !* just the reference constituent.
2048  !*
2049  !* 2. An extension to permit the analysis of current observations.
2050  !*
2051  !* 3. Removal of a central time as the basis for the calculation of the
2052  !* astronomical arguments V. Now the V value for each observation, as well
2053  !* as those for the nodal corrections f and u (done for the JAOT analysis),
2054  !* are incorporated directly into the overdetermined matrix. These changes
2055  !* mean that analyses no longer need be restricted to periods of a year or
2056  !* less. (Though as the period approaches 18.6 years, using another "long
2057  !* period" analysis program that solves for the "nodal satellites" directly
2058  !* is advisable.)
2059  !*
2060  !***********************************************************************
2061  !*
2062  !* FILE REFERENCE NUMBERS OF DEVICES REQUIRED BY THIS PROGRAM.
2063  !* KR - INPUT FILE - CONTAINS THE TIDAL CONSTITUENT INFORMATION.
2064  !* KR1 - INPUT FILE - GIVES ANALYSIS TYPE AND TIDAL STATION
2065  !* DETAILS.
2066  !* KR2 - INPUT FILE - CONTAINS THE OBSERVED TIMES AND HEIGHTS.
2067  !* PRESENTLY KR,KR1,KR2, AND LP ARE ASSIGNED THE RESPECTIVE VALUES
2068  !* 8,5,9, AND 6. SEE THE MANUAL OR COMMENT STATEMENTS WITHIN THIS
2069  !* PROGRAM FOR FURTHER DETAILS ON THEIR USE.
2070  !*
2071  !***********************************************************************
2072  !*
2073  !* ARRAY DEFINITIONS AND DIMENSION GUIDELINES.
2074  !*
2075  !* LET MC BE THE TOTAL NUMBER OF CONSTITUENTS, INCLUDING Z0
2076  !* AND ANY INFERRED CONSTITUENTS, TO BE INCLUDED IN THE
2077  !* ANALYSIS; (For T/P, MC=30 > NUMBER OF CONSTITUENTS)
2078  !* TIDE_NTI BE THE NUMBER OF TIDAL HEIGHT OBSERVATIONS;
2079  !* NR BE THE NUMBER OF INPUT RECORDS OF OBSERVED TIDAL
2080  !* HEIGHTS (For T/P data, same as TIDE_NTI, set to 200)
2081  !* MPAR BE 2*MC-1;
2082  !* NEQ BE TIDE_NTI*2 IF ALL THE OBSERVATIONS ARE EXTREMES AND
2083  !* THE DERIVATIVE CONDITION IS TO BE INCLUDED FOR EACH,
2084  !* AND TIDE_NTI OTHERWISE (= NR for T/P data).
2085  !* THEN PARAMETERS NMAXP1, AND NMAXPM SHOULD BE AT LEAST MPAR+1, AND
2086  !* NEQ+MPAR RESPECTIVELY. THEY ARE CURRENTLY SET TO 40 AND 240.
2087  !*
2088  !* TIDECON_NAME(I) IS THE ARRAY CONTAINING ALL THE CONSTITUENT NAMES,
2089  !* INCLUDING Z0 AND ANY INFERRED CONSTITUENTS, TO BE IN
2090  !* THE ANALYSIS. IT SHOULD BE DIMENSIONED AT LEAST MC.
2091  !* TIDE_FREQC(I), ARE THE ARRAYS OF FREQUENCIES IN CYCLES/HR AND
2092  !* FREQ(I) RADIANS/HR RESPECTIVELY CORRESPONDING TO THE
2093  !* CONSTITUENT NAME(I). THEY SHOULD BE DIMENSIONED AT
2094  !* LEAST MC.
2095  !* AMP(I),PH(I) ARE ARRAYS CONTAINING THE RAW AMPLITUDE AND PHASE FOR
2096  !* CONSTITUENT NAME(I) AS FOUND VIA THE LEAST SQUARES
2097  !* ANALYSIS. THEY SHOULD BE DIMENSIONED AT LEAST MC.
2098  !* AMPC(I),PHG(I) ARE ARRAYS CONTAINING THE AMPLITUDE AND PHASE FOR
2099  !* CONSTITUENT NAME(I) AFTER CORRECTIONS FOR NODAL
2100  !* MODULATION, ASTRONOMICAL ARGUMENT AND INFERRED
2101  !* CONSTITUENTS. THEIR MINIMUM DIMENSION SHOULD BE MC.
2102  !* TIDE_DATA(I) AND HEIGHTS, OF THE OBSERVED DATA AS IT IS INPUT BY
2103  !* RECORD. THEY SHOULD BE DIMENSIONED ACCORDINGLY( AT
2104  !* PRESENT ONLY 6 OBSERVATIONS ARE EXPECTED PER RECORD).
2105  !* X(I) ARRAY CONTAINING ALL THE TIMES(IN HOURS AS
2106  !* MEASURED FROM THE CENTRE OF THE ANALYSIS PERIOD) ITS MINIMUM
2107  !* DIMENSION SHOULD BE TIDE_NTI.
2108  !* NSTN(I) IS THE ARRAY CONTAINING THE TIDAL STATION TIDECON_NAME. IT
2109  !* SHOULD HAVE MINIMUM DIMENSION 5.
2110  !* Q(I) IS THE OVERDETERMINED ARRAY OF EQUATIONS THAT IS
2111  !* SOLVED IN THE LEAST SQUARES SENSE BY THE MODIFIED
2112  !* GRAM-SCHMIDT ALGORITHM. IT SHOULD HAVE THE EXACT
2113  !* DIMENSION OF NMAXPM BY NMAXP1.
2114  !* P(I) IS THE ARRAY CONTAINING THE TIDAL CONSTITUENT SINE
2115  !* AND COSINE COEFICIENTS AS FOUND WITH THE LEAST
2116  !* SQUARES FIT. IT SHOULD HAVE MINIMUM DIMENSION MPAR.
2117  !
2118  !***********************************************************************
2119  !
2120  IMPLICIT NONE
2121 
2122  INTEGER, INTENT(IN) :: NDEF, ITREND, ITEST
2123  INTEGER(KIND=4), INTENT(IN) :: KD1, KD2, IX, IY
2124  REAL , INTENT(IN) :: XLON, XLAT
2125  REAL(KIND=8) , intent(out):: sdev0(ndef), rmsr0(ndef), rmsr(ndef), resmax(ndef)
2126  REAL , INTENT(OUT):: SSQ(NDEF), RES(NDEF)
2127  INTEGER , INTENT(OUT):: IMAX(NDEF)
2128  !
2129  INTEGER :: I, I1, I2, I21, II1, IDEF, ICODE, INFLAG, &
2130  J, INFTOT, IREP, J2, JCODE, JJ1, K, K2, KH1, &
2131  KH2, KHM, KINF, L, M, MEQ, N, NCOL, NEW, NMAX
2132  REAL(KIND=8) :: aamp, arg, arg1, arg2, arg3, c, c2, c3, &
2133  fx, fxi, s, s2, s3, ux, vx, uxi, vxi, &
2134  wmin, wmax, xmid
2135  REAL :: TOLER
2136  REAL(KIND=8) :: av, sdev, sum2, hrm
2137  DOUBLE PRECISION :: X(NR),Y(NR), TIME(NR)
2138  REAL :: Q(NMAXPM,NMAXP1),FREQ(MC),AMP(MC),PH(MC)
2139  DOUBLE PRECISION :: P(NMAXP1),CENHR,CUMHR
2140  DOUBLE PRECISION :: yy
2141  !
2142  ! Additional arrays, for use in the SVD routine (J.Ch., Aug. 1997)
2143  DOUBLE PRECISION :: U(NMAXPM,NMAXP1),V(NMAXP1,NMAXP1), &
2144  COV(NMAXP1,NMAXP1),B(NMAXPM),W(NMAXP1),SIG(NMAXPM)
2145  !
2146  !***********************************************************************
2147  !
2148 
2149  kh1=24*kd1
2150  kh2=24*(kd2+1)
2151  khm=(kh1+kh2)/2
2152  hrm=khm
2153  cenhr=dfloat((kh2-kh1)/2)
2154 
2155  IF (itrend.eq.1) then
2156  m=2*tide_mf
2157  else
2158  m=2*tide_mf-1
2159  END IF
2160 
2161 
2162  i=0
2163  DO i=1,tide_mf
2164  freq(i)=tide_freqc(i)*twpi
2165  END DO
2166  !
2167  !***********************************************************************
2168  !* DETERMINE THE CENTRAL HOUR OF THE ANALYSIS PERIOD AND SET UP THE
2169  !* DEPENDENT AND INDEPENDENT VARIABLES, Y AND X.
2170  ! actually, CUMHR=24.d0*(KD-KHM) (check), but keep same notation as before
2171  !
2172  k=1
2173  DO i=1,tide_nti
2174  cumhr=-cenhr+24.d0*(tide_days(i)-kd1)
2175  time(i)=cumhr+dfloat(tide_secs(i))/3600.d0
2176  x(k)=time(i)-time(1)
2177  n=k
2178  k=k+1
2179  END DO
2180 
2181  !
2182  !***********************************************************************
2183  !* SETTING UP THE OVERDETERMINED MATRIX AND SOLVING WITH MODIFIED SVD
2184  !
2185  irep=0
2186 
2187  DO idef=1,ndef ! loop thru once or twice
2188 
2189  ! Modifies the time reference xmid and puts it at 0 : modification by FA, 2012/09/26
2190  ! the impact of that change has not been verified ...
2191  xmid=0. !0.5*(TIDE_HOURS(1)+TIDE_HOURS(TIDE_NTI))
2192  q(1:nmaxpm,1:nmaxp1)=0.0
2193  DO i=1,n
2194  ! if itrend=1, then
2195  ! first 2 parameters are constant and linear trend (per 365 days)
2196  ! fitted as const+trend(t-tmid) where tmid (=xmid) is the middle time
2197  ! of the analysis period (This makes the constant consistent with z0
2198  ! in the old analysis program)
2199  ! If itrend=0 then the second parameter is is associated with the next
2200  ! constituent
2201  q(i,1)=1.
2202  IF (itrend.eq.1) then
2203  ! Q(I,2)=(x(i)-xmid)/(24.*365.)
2204  q(i,2)=x(i)/(24.*365.)
2205  END IF
2206  q(i,nmaxp1)=tide_data(i,idef)
2207  icode=1
2208  ! should only have to assemble lhs of matrix when idef=1
2209  ! but something is not right if don't do it 2nd time too
2210  ! CALL SETVUF(TIDE_HOURS(I),xlat,I)
2211  DO j=2,tide_mf
2212  CALL vuf(tidecon_name(j),vx,ux,fx,i)
2213  ! check to see if this constituent is to be used for inference
2214  inflag=0
2215  kinf=0
2216  IF (tide_nin.GE.0) THEN
2217  do k=1,tide_nin
2218  IF (tidecon_name(j).eq.tide_konan(k)) then
2219  inflag=1
2220  kinf=k
2221  EXIT
2222  END IF
2223  END DO
2224  END IF
2225  !
2226  IF (inflag.eq.0) then
2227  arg=(vx+ux)*twpi
2228  IF (itrend.eq.1) then
2229  j2=2*(j-1)+1
2230  else
2231  j2=2*(j-1)
2232  END IF
2233  jj1=j2+1
2234  q(i,j2)=cos(arg)*fx
2235  q(i,jj1)=sin(arg)*fx
2236  else
2237  IF (itrend.eq.1) then
2238  j2=2*(j-1)+1
2239  else
2240  j2=2*(j-1)
2241  END IF
2242  jj1=j2+1
2243  arg1=(vx+ux)*twpi
2244  q(i,j2)=cos(arg1)*fx
2245  q(i,jj1)=sin(arg1)*fx
2246  do k2=1,tide_ninf(kinf)
2247  CALL vuf(tide_konin(kinf,k2),vxi,uxi,fxi,i)
2248  ! freq is radians/hr but sigin is cycles/hr
2249  arg2=(vxi+uxi)*twpi
2250  c2=cos(arg2)
2251  s2=sin(arg2)
2252  arg3=tide_zeta(kinf,k2)*fac
2253  c3=cos(arg3)
2254  s3=sin(arg3)
2255  q(i,j2)=q(i,j2)+fxi*tide_r(kinf,k2)*(c2*c3-s2*s3)
2256  q(i,jj1)=q(i,jj1)+fxi*tide_r(kinf,k2)*(c2*s3+s2*c3)
2257  END DO
2258  END IF
2259  END DO ! j
2260  END DO !i
2261 
2262 #ifdef W3_T
2263  WRITE(6,*) 'assembled overdetermined matrix and/or rhs'
2264 #endif
2265  nmax=m
2266  meq=n
2267  ssq(idef)=1.0
2268  res(idef)=1.0
2269  ncol=nmax
2270  new=nmax
2271  !
2272  !***********************************************************************
2273  !* CALCULATION OF THE STANDARD DEVIATION OF THE RIGHT HAND SIDES OF
2274  !* THE OVERDETERMINED SYSTEM
2275  !
2276  av=0.d0
2277  DO i=1,meq
2278  av=av+q(i,nmaxp1)
2279  END DO
2280  av=av/meq
2281  sdev=0.d0
2282  DO i=1,meq
2283  sdev=sdev+(q(i,nmaxp1)-av)**2
2284  END DO
2285  sdev=sdev/(meq-1)
2286  sdev=sqrt(sdev)
2287  sdev0(idef)=sdev
2288 109 CONTINUE
2289  !
2290  ! USE SINGULAR-VALUE-DECOMPOSITION TO SOLVE THE OVERDETERMINED SYSTEM
2291  !
2292  toler=1.e-5
2293  DO i=1,nmaxpm
2294  sig(i)=1.d0
2295  END DO
2296  !
2297  ! no solution if meq lt m. ie underdetermined system
2298  ! go to next time series
2299  IF (meq.le.m) then
2300  write(ndset,*) ' underdetermined system: no svd solution',ix,iy,meq,m
2301  stop
2302  END IF
2303 #ifdef W3_T
2304  WRITE(6,*) ' applying svd'
2305 #endif
2306  CALL svd(q,u,v,cov,w,p,b,sig,icode,meq,nmax,nmaxpm,nmaxp1,toler &
2307  ,jcode,ssq(idef),res(idef))
2308  ! IF (JCODE.GT.0) WRITE(LP,55)JCODE
2309  ! 55 FORMAT('COLUMN',I5,' IS THE 1ST DEPENDENT COLUMNS IN SVD')
2310  ! write out eigenvalues
2311  wmax=-1000.
2312  wmin=1000.
2313  do i=1,nmax
2314  IF (w(i).gt.wmax) wmax=w(i)
2315  IF (w(i).lt.wmin) wmin=w(i)
2316  end do
2317  ! write(6,*) ' max, min eigenvalues =',wmax,wmin
2318  ! write(6,*) ' all eigenvalues'
2319  ! write(6,56) (w(i),i=1,nmax)
2320 56 format(10e12.5)
2321  !***********************************************************************
2322  IF (ssq(idef).gt.1.e-10) then
2323  rmsr0(idef)=sqrt(ssq(idef)/(meq-m))
2324  else
2325  rmsr0(idef)=0.
2326  END IF
2327 
2328 
2329  rmsr(idef)=0.d0
2330  resmax(idef)=0.
2331  do i=1,n
2332  yy=q(i,nmaxp1)
2333  rmsr(idef)=rmsr(idef)+yy*yy
2334  IF (abs(yy).gt.resmax(idef)) then
2335  resmax(idef)=abs(yy)
2336  imax(idef)=i
2337  END IF
2338  end do
2339 160 format(' ',7i2,f15.5,f10.5,i6)
2340  IF (rmsr(idef).gt.1.e-10) then
2341  rmsr(idef)=dsqrt(rmsr(idef)/(n-m))
2342  else
2343  rmsr(idef)=0.
2344  END IF
2345  ! close(unit=25)
2346  !
2347  !***********************************************************************
2348  !* CALCULATE AMPLITUDES AND PHASES
2349  !
2350  ! if itrend=1 then the linear trend is shown as the phase of the constant
2351  ! Z0 term (& the true phase of Z0 is zero)
2352  ! otherwise, the phase of Z0 is shown as zero
2353  amp(1)=p(1)
2354  IF (itrend.eq.1) then
2355  ph(1)=p(2)
2356  else
2357  ph(1)=0.
2358  END IF
2359  DO i=2,tide_mf
2360  !
2361  IF (itrend.eq.1) then
2362  i2=2*(i-1)+1
2363  else
2364  i2=2*(i-1)
2365  END IF
2366  i21=i2+1
2367  c=p(i2)
2368  s=p(i21)
2369  aamp=sqrt(c*c+s*s)
2370  IF (aamp.LT.1.e-5) THEN
2371  ph(i)=0.
2372  ELSE
2373  ph(i)=atan2(s,c)/fac
2374  IF (ph(i).LT.0.) ph(i)=ph(i)+360.
2375  END IF
2376  amp(i)=aamp
2377  END DO ! end of loop on TIDE_MF
2378  !***********************************************************************
2379  ! Note that with f & u included in the lsq fit, we only need V from routine VUF
2380  ! but we don't want to correct with V for a central hour. Better to include
2381  ! the right V in the lsq fit. This has been done.
2382 
2383  tide_ampc(1,idef)=amp(1)
2384  tide_phg(1,idef)=ph(1)
2385  DO i=2,tide_mf
2386  tide_ampc(i,idef)=amp(i)
2387  tide_phg(i,idef)=ph(i)
2388  END DO
2389  !
2390  tide_sig3(:,idef)=0.
2391  tide_ttest(:,idef)=0.
2392  !
2393  IF (itest.GE.1) THEN
2394  !---------------------------------------------------
2395  !
2396  i=1
2397  IF (cov(1,1).gt.1.e-8) then
2398  tide_sig1(i,idef)=sqrt(cov(1,1))*rmsr0(idef)
2399  else
2400  tide_sig1(i,idef)=0.
2401  END IF
2402  IF (itrend.eq.1.and.cov(2,2).gt.1.e-8) then
2403  tide_sig2(i,idef)=sqrt(cov(2,2))*rmsr0(idef)
2404  else
2405  tide_sig2(i,idef)=0.
2406  END IF
2407  tide_sig3(i,idef)=0.
2408  tide_ttest(i,idef)=0.
2409  !
2410  ! results for the other constituents
2411  !
2412  DO i=2,tide_mf
2413  IF (itrend.eq.1) then
2414  i2=2*(i-1)+1
2415  else
2416  i2=2*(i-1)
2417  END IF
2418  ii1=i2+1
2419  !
2420  ! multiply cov values with residual standard deviation, as described in equation
2421  ! (6) of Cherniasky et al. (2001)
2422  !
2423  IF (cov(i2,i2).gt.1.e-8) then
2424  tide_sig1(i,idef)=sqrt(cov(i2,i2))*rmsr0(idef)
2425  else
2426  tide_sig1(i,idef)=0.
2427  END IF
2428  IF (cov(ii1,ii1).gt.1.e-8) then
2429  tide_sig2(i,idef)=sqrt(cov(ii1,ii1))*rmsr0(idef)
2430  else
2431  tide_sig2(i,idef)=0.
2432  END IF
2433  ! from equation 11 in Pawlowicz et al (2002)
2434  c=tide_ampc(i,idef)*cos(tide_phg(i,idef)*fac)
2435  s=tide_ampc(i,idef)*sin(tide_phg(i,idef)*fac)
2436  tide_sig3(i,idef)=sqrt(((c*tide_sig1(i,idef))**2+(s*tide_sig2(i,idef))**2)/(c**2+s**2))
2437  tide_ttest(i,idef)=tide_ampc(i,idef)/tide_sig3(i,idef)
2438  END DO
2439  !---------------------------------------------------
2440  END IF ! (ITEST.GE.1)
2441  !
2442  ! now inferred constituents
2443  !
2444  IF (tide_nin.GE.0) THEN
2445  l=0
2446  do k=1,tide_nin
2447  do i=2,tide_mf
2448  IF (tidecon_name(i).eq.tide_konan(k)) EXIT
2449  END DO
2450  i1=i
2451  do k2=1,tide_ninf(k)
2452  l=l+1
2453  tide_ampci(k,k2,idef)=tide_ampc(i1,idef)*tide_r(k,k2)
2454  tide_phgi(k,k2,idef)=tide_phg(i1,idef)-tide_zeta(k,k2)
2455  END DO
2456  END DO
2457  inftot=l
2458  END IF
2459 
2460 
2461  !
2462  !***********************************************************************
2463  ! compute (Cherniawsky et al (2001), page 653) and rank correlation coefficients
2464  ! largest niter value are computed and shown
2465  ! if itrend=1, then the second part of Z0 is the linear trend coefficient
2466  !
2467  ! do i=1,m
2468  ! do j=1,i
2469  ! cor(i,j)=cov(i,j)/sqrt(cov(i,i)*cov(j,j))
2470  ! END DO
2471  ! END DO
2472  !
2473  ! niter=20
2474  ! do 81 iter=1,niter
2475  ! cormax=0.
2476  ! do i=2,m
2477  ! im1=i-1
2478  ! do j=1,im1
2479  ! ac=abs(cor(i,j))
2480  ! if (ac.gt.cormax) then
2481  ! cormax=ac
2482  ! imax=i
2483  ! jmax=j
2484  ! END IF
2485  ! END DO
2486  ! END DO
2487  ! if (itrend.eq.1) then
2488  ! iconst=(imax+1)/2
2489  ! jconst=(jmax+1)/2
2490  ! else
2491  ! iconst=(imax+2)/2
2492  ! jconst=(jmax+2)/2
2493  ! END IF
2494  ! write(lp,83) iter,cormax,imax,jmax,TIDECON_NAME(iconst),TIDECON_NAME(jconst)
2495  !83 format(i5,' largest correlation coefficient is ',f8.3,' at (i,j)=' &
2496  ! ,2i5,' for constituents ',a5,' and ',a5)
2497  ! cor(imax,jmax)=0.
2498  ! END DO
2499  !
2500  END DO ! end of loop on IDEF
2501  RETURN
2502 
2503  END SUBROUTINE flex_tidana_webpage
2504 
2505 
2506 
2507 
2508 
2509 
2510  !/ ------------------------------------------------------------------- /
2511  SUBROUTINE tide_predict(itrend,ndef,N,HOURS, DATAIN, PREDICTED, RESID, XLAT,SDEV,RMSR)
2512  !
2513  IMPLICIT NONE
2514  !
2515  INTEGER, INTENT(IN) :: itrend, NDEF, N
2516  REAL(KIND=8) , intent(in) :: hours(n)
2517  REAL, INTENT(IN) :: XLAT, DATAIN(N,NDEF)
2518  REAL(KIND=8), intent(out) :: sdev(ndef),rmsr(ndef)
2519  REAL, INTENT(OUT) :: PREDICTED(N,NDEF), RESID(N,NDEF)
2520  !
2521  INTEGER :: IDEF, I, K, K2
2522  INTEGER :: J, M
2523  REAL :: ARG, ADD, SUM1, SSQ
2524  REAL(KIND=8) :: vx, ux, fx
2525 
2526  m = 2*tide_mf
2527 
2528  DO idef=1,ndef
2529  sdev=0.d0
2530  DO i=1,n
2531  IF (itrend.eq.1) THEN
2532  sum1=tide_ampc(1,idef)+tide_phg(1,idef)*hours(i)/(365.*24.)
2533  ELSE
2534  sum1=tide_ampc(1,idef)
2535  END IF
2536  !
2537  DO j=2,tide_mf
2538  CALL vuf(tidecon_name(j),vx,ux,fx,i)
2539  arg=(vx+ux)*twpi-tide_phg(j,idef)*fac
2540  add=fx*tide_ampc(j,idef)*cos(arg)
2541  sum1=sum1+add
2542  END DO
2543  !
2544  IF (tide_nin.NE.0) THEN
2545  DO k=1,tide_nin
2546  DO k2=1,tide_ninf(k)
2547  CALL vuf(tide_konin(k,k2),vx,ux,fx,i)
2548  arg=(vx+ux)*twpi-tide_phgi(k,k2,idef)*fac
2549  add=fx*tide_ampci(k,k2,idef)*cos(arg)
2550  sum1=sum1+add
2551  END DO
2552  END DO
2553  END IF
2554  !
2555  predicted(i,idef)=sum1
2556  !
2557  resid(i,idef)=datain(i,idef)-sum1
2558  sdev(idef)=sdev(idef)+resid(i,idef)**2
2559  END DO
2560  !
2561  ssq=sdev(idef)
2562  rmsr(idef)=sqrt(ssq/(n-m))
2563  sdev(idef)=sqrt(ssq/n)
2564  ENDDO
2565  !
2566  END SUBROUTINE tide_predict
2567 
2568  !/ ------------------------------------------------------------------- /
2569  SUBROUTINE tide_predict_only(itrend,ndef,N,TIDE_HOURS, PREDICTED, XLAT)
2570  !
2571  IMPLICIT NONE
2572  !
2573  INTEGER, INTENT(IN) :: itrend, NDEF, N
2574  REAL(KIND=8) , intent(in) :: tide_hours(n)
2575  REAL, INTENT(IN) :: XLAT
2576  REAL, INTENT(OUT) :: PREDICTED(N,NDEF)
2577  !
2578  INTEGER :: IDEF, I, K, K2
2579  INTEGER :: J
2580  REAL :: ARG, ADD, SUM1
2581  REAL(KIND=8) :: vx, ux, fx
2582 
2583  DO idef=1,ndef
2584  DO i=1,n
2585  IF (itrend.eq.1) THEN
2586  sum1=tide_ampc(1,idef)+tide_phg(1,idef)*tide_hours(i)/(365.*24.)
2587  ELSE
2588  sum1=tide_ampc(1,idef)
2589  END IF
2590  !
2591  DO j=2,tide_mf
2592  CALL vuf(tidecon_name(j),vx,ux,fx,i)
2593  arg=(vx+ux)*twpi-tide_phg(j,idef)*fac
2594  add=fx*tide_ampc(j,idef)*cos(arg)
2595  sum1=sum1+add
2596  END DO
2597  !
2598  IF (tide_nin.NE.0) THEN
2599  DO k=1,tide_nin
2600  DO k2=1,tide_ninf(k)
2601  CALL vuf(tide_konin(k,k2),vx,ux,fx,i)
2602  arg=(vx+ux)*twpi-tide_phgi(k,k2,idef)*fac
2603  add=fx*tide_ampci(k,k2,idef)*cos(arg)
2604  sum1=sum1+add
2605  END DO
2606  END DO
2607  END IF
2608  !
2609  predicted(i,idef)=sum1
2610  !
2611  END DO
2612  !
2613  ENDDO
2614  !
2615  END SUBROUTINE tide_predict_only
2616 
2617  !/
2618  !/ End of module WMTIDEMD -------------------------------------------- /
2619  !/
2620 END MODULE w3tidemd
w3tidemd::vuf_set_parameters
subroutine vuf_set_parameters
Definition: w3tidemd.F90:1445
w3tidemd::ee
real, dimension(180) ee
Definition: w3tidemd.F90:102
w3tidemd::tide_ttest
real, dimension(mc, 2) tide_ttest
Definition: w3tidemd.F90:129
constants::pi
real, parameter pi
PI Value of Pi.
Definition: constants.F90:71
w3tidemd::tidecon_namei
character(len=5), dimension(:), allocatable tidecon_namei
Definition: w3tidemd.F90:111
w3tidemd::setvuf
subroutine setvuf(hr, XLAT, ITIME)
Definition: w3tidemd.F90:1891
w3tidemd::tide_indexj
integer, dimension(:), allocatable tide_indexj
Definition: w3tidemd.F90:105
w3tidemd::vuf
subroutine vuf(KONX, Vx, ux, FX, ITIME)
Definition: w3tidemd.F90:1705
w3tidemd::nr
integer, parameter nr
Definition: w3tidemd.F90:92
w3tidemd::ldel
integer, dimension(180) ldel
Definition: w3tidemd.F90:103
w3tidemd::tide_verbose
integer tide_verbose
Definition: w3tidemd.F90:136
w3tidemd::ndset
integer ndset
Definition: w3tidemd.F90:136
w3tidemd::tide_nx
integer tide_nx
Definition: w3tidemd.F90:109
w3tidemd::tide_sigan
real, dimension(10) tide_sigan
Definition: w3tidemd.F90:115
w3tidemd
Definition: w3tidemd.F90:3
w3tidemd::tide_index
integer, dimension(mc) tide_index
Definition: w3tidemd.F90:133
w3tidemd::v_arg
real, dimension(:,:), allocatable v_arg
Definition: w3tidemd.F90:101
w3tidemd::flex_tidana_webpage
subroutine flex_tidana_webpage(IX, IY, XLON, XLAT, KD1, KD2, ndef, itrend, RES, SSQ, RMSR0, SDEV0, RMSR, RESMAX, IMAX, ITEST)
Definition: w3tidemd.F90:2033
w3tidemd::nkonco
integer nkonco
Definition: w3tidemd.F90:96
w3tidemd::tide_find_indices_prediction
subroutine tide_find_indices_prediction(LIST, INDS, TIDE_PRMF)
Definition: w3tidemd.F90:244
w3tidemd::tidal_const
real, dimension(:,:,:,:,:), allocatable tidal_const
Definition: w3tidemd.F90:117
w3tidemd::tide_sigin
real, dimension(10, 10) tide_sigin
Definition: w3tidemd.F90:115
w3tidemd::tide_zeta
real, dimension(10, 10) tide_zeta
Definition: w3tidemd.F90:114
w3tidemd::tide_ampci
real, dimension(10, 10, 2) tide_ampci
Definition: w3tidemd.F90:132
w3odatmd::iaproc
integer, pointer iaproc
Definition: w3odatmd.F90:457
w3tidemd::tide_write_results
subroutine tide_write_results(LP, filename, ndef, KD1, KD2, ITZ, xlat, xlon, RES, SSQ, RMSR0, SDEV0, SDEV, RMSR, RESMAX, IMAX, RMSRP)
Definition: w3tidemd.F90:146
w3tidemd::nmaxpm
integer, parameter nmaxpm
Definition: w3tidemd.F90:92
w3tidemd::ntidal_con
integer ntidal_con
Definition: w3tidemd.F90:96
w3tidemd::konco_con
character *5, dimension(:), allocatable konco_con
Definition: w3tidemd.F90:98
w3tidemd::tide_hours
real(kind=8), dimension(:), allocatable tide_hours
Definition: w3tidemd.F90:124
w3tidemd::ph
real, dimension(180) ph
Definition: w3tidemd.F90:102
w3tidemd::mc
integer, parameter mc
Definition: w3tidemd.F90:92
w3tidemd::tide_set_indices
subroutine tide_set_indices
Definition: w3tidemd.F90:495
w3tidemd::tide_r
real, dimension(10, 10) tide_r
Definition: w3tidemd.F90:114
w3tidemd::tide_read_settings
subroutine tide_read_settings(filename, fnam6, fnam7, fnam8, fnam9, fnam11)
Definition: w3tidemd.F90:672
w3tidemd::tide_sig2
real, dimension(mc, 2) tide_sig2
Definition: w3tidemd.F90:129
w3tidemd::ndel
integer, dimension(180) ndel
Definition: w3tidemd.F90:103
w3odatmd::ndse
integer, pointer ndse
Definition: w3odatmd.F90:456
w3tidemd::tide_index2
integer, dimension(mc) tide_index2
Definition: w3tidemd.F90:133
w3tidemd::tide_ny
integer tide_ny
Definition: w3tidemd.F90:109
w3tidemd::coef_con
real, dimension(:), allocatable coef_con
Definition: w3tidemd.F90:100
w3tidemd::tide_sig3
real, dimension(mc, 2) tide_sig3
Definition: w3tidemd.F90:129
w3odatmd::naperr
integer, pointer naperr
Definition: w3odatmd.F90:457
w3tidemd::tide_nin
integer tide_nin
Definition: w3tidemd.F90:116
w3tidemd::juldayt
integer *4 function juldayt(id, mm, iyyy)
Definition: w3tidemd.F90:1257
w3servmd
Definition: w3servmd.F90:3
w3tidemd::tide_indexjk
integer, dimension(:), allocatable tide_indexjk
Definition: w3tidemd.F90:105
w3tidemd::tidecon_allnames
character *5, dimension(:), allocatable tidecon_allnames
Definition: w3tidemd.F90:97
w3tidemd::dpythag
double precision function dpythag(a, b)
Definition: w3tidemd.F90:977
w3tidemd::astr
subroutine astr(d1, h, pp, s, p, np, dh, dpp, ds, dp, dnp)
Definition: w3tidemd.F90:911
w3tidemd::u_arg
real, dimension(:,:), allocatable u_arg
Definition: w3tidemd.F90:101
w3tidemd::tide_nti
integer(kind=4) tide_nti
Definition: w3tidemd.F90:121
w3odatmd
Definition: w3odatmd.F90:3
w3tidemd::dsvbksb
subroutine dsvbksb(u, w, v, m, n, mp, np, b, x)
Definition: w3tidemd.F90:996
w3tidemd::tide_predict
subroutine tide_predict(itrend, ndef, N, HOURS, DATAIN, PREDICTED, RESID, XLAT, SDEV, RMSR)
Definition: w3tidemd.F90:2512
w3tidemd::kblank
character *5, parameter kblank
Definition: w3tidemd.F90:94
w3tidemd::tide_konan
character(len=5), dimension(10) tide_konan
Definition: w3tidemd.F90:113
w3odatmd::naproc
integer, pointer naproc
Definition: w3odatmd.F90:457
w3tidemd::tide_read_anapar
subroutine tide_read_anapar(KR1, LP, filename, KD1, KD2, XLON, XLAT, NDEF, ITREND, ITZ)
Definition: w3tidemd.F90:717
w3tidemd::tide_find_indices_analysis
subroutine tide_find_indices_analysis(LIST)
Definition: w3tidemd.F90:343
file
file(STRINGS ${CMAKE_BINARY_DIR}/switch switch_strings) separate_arguments(switches UNIX_COMMAND $
Definition: CMakeLists.txt:3
w3tidemd::caldatt
subroutine caldatt(julian, id, mm, iyyy)
Definition: w3tidemd.F90:1284
w3tidemd::tide_predict_only
subroutine tide_predict_only(itrend, ndef, N, TIDE_HOURS, PREDICTED, XLAT)
Definition: w3tidemd.F90:2570
w3servmd::strace
subroutine strace(IENT, SNAME)
Definition: w3servmd.F90:148
w3tidemd::setvuf_fast
subroutine setvuf_fast(h, pp, s, p, enp, dh, dpp, ds, dp, dnp, tau, XLAT, F, U, V)
Definition: w3tidemd.F90:532
w3odatmd::ndso
integer, pointer ndso
Definition: w3odatmd.F90:456
w3tidemd::tide_sig1
real, dimension(mc, 2) tide_sig1
Definition: w3tidemd.F90:129
w3odatmd::napout
integer, pointer napout
Definition: w3odatmd.F90:457
w3tidemd::tidecon_name
character(len=5), dimension(:), allocatable tidecon_name
Definition: w3tidemd.F90:112
w3tidemd::tide_ninf
integer, dimension(10) tide_ninf
Definition: w3tidemd.F90:116
w3tidemd::tide_mf
integer tide_mf
Definition: w3tidemd.F90:109
w3tidemd::svd
subroutine svd(q, u, v, cov, w, p, b, sig, ic, m, n, mm, N2, toler, jc, ssq, res)
Definition: w3tidemd.F90:1314
w3tidemd::f_arg
real, dimension(:,:), allocatable f_arg
Definition: w3tidemd.F90:101
w3tidemd::semi
real, dimension(:), allocatable semi
Definition: w3tidemd.F90:100
w3tidemd::mdel
integer, dimension(180) mdel
Definition: w3tidemd.F90:103
w3tidemd::tide_read_timeseries
subroutine tide_read_timeseries(KR2, filename, KD1, KD2, TIDE_NTI, NDEF)
Definition: w3tidemd.F90:844
w3tidemd::tide_days
integer(kind=4), dimension(:), allocatable tide_days
Definition: w3tidemd.F90:123
w3tidemd::tide_freqc
real, dimension(:), allocatable tide_freqc
Definition: w3tidemd.F90:110
w3tidemd::ir
integer, dimension(180) ir
Definition: w3tidemd.F90:103
w3tidemd::dsvdcmp
subroutine dsvdcmp(a, m, n, mp, np, w, v)
Definition: w3tidemd.F90:1024
w3tidemd::nmaxp1
integer, parameter nmaxp1
Definition: w3tidemd.F90:92
w3tidemd::twpi
double precision, parameter twpi
Definition: w3tidemd.F90:86
w3tidemd::tide_phgi
real, dimension(10, 10, 2) tide_phgi
Definition: w3tidemd.F90:132
w3tidemd::tide_ampc
real, dimension(mc, 2) tide_ampc
Definition: w3tidemd.F90:129
w3tidemd::ntotal_con
integer ntotal_con
Definition: w3tidemd.F90:96
w3tidemd::fac
double precision, parameter fac
Definition: w3tidemd.F90:87
w3tidemd::tide_data
real, dimension(:,:), allocatable tide_data
Definition: w3tidemd.F90:122
w3tidemd::tide_phg
real, dimension(mc, 2) tide_phg
Definition: w3tidemd.F90:129
w3tidemd::tide_dt
real, parameter tide_dt
Definition: w3tidemd.F90:125
w3tidemd::mc2
integer, parameter mc2
Definition: w3tidemd.F90:93
w3tidemd::opnvuf
subroutine opnvuf(filename)
Definition: w3tidemd.F90:1747
w3tidemd::tide_konin
character(len=5), dimension(10, 10) tide_konin
Definition: w3tidemd.F90:113
w3tidemd::tide_secs
integer(kind=4), dimension(:), allocatable tide_secs
Definition: w3tidemd.F90:123