New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

tide_mod.F90 in NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/SBC – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/SBC/tide_mod.F90 @ 15547

Last change on this file since 15547 was 15547, checked in by hadjt, 3 years ago
Moving ftide, utide, v0tide output to sbctides, adding tide_t output to diaharm_fast
File size: 75.8 KB

Line
1	MODULE tide_mod
2	!!======================================================================
3	!! * MODULE tide_mod *
4	!! Compute nodal modulations corrections and pulsations
5	!!======================================================================
6	!! History : 1.0 ! 2007 (O. Le Galloudec) Original code
7	!!----------------------------------------------------------------------
8	!JT USE oce ! ocean dynamics and tracers
9	USE lib_mpp ! distributed memory computing library
10
11	USE dom_oce ! ocean space and time domain
12	USE phycst ! physical constant
13	USE daymod ! calendar
14	USE in_out_manager ! I/O units
15	USE ioipsl , ONLY : ymds2ju ! for calendar
16
17	USE iom
18
19	IMPLICIT NONE
20	PRIVATE
21
22	PUBLIC tide_harmo ! called by tideini and diaharm modules
23	PUBLIC tide_init_Wave ! called by tideini and diaharm modules
24	PUBLIC tide_init_calendar_options ! called by tideini and diaharm modules
25
26	! davbyr: increase maximum number of harmonics from 19 to 34
27	INTEGER, PUBLIC, PARAMETER :: jpmax_harmo = 34 !: maximum number of harmonic
28
29
30	REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) :: v0linearintercept
31
32
33
34	TYPE, PUBLIC :: tide
35	CHARACTER(LEN=4) :: cname_tide
36	REAL(wp) :: equitide
37	INTEGER :: nutide
38	INTEGER :: nt, ns, nh, np, np1, shift
39	INTEGER :: nksi, nnu0, nnu1, nnu2, R
40	INTEGER :: nformula
41	END TYPE tide
42
43	TYPE(tide), PUBLIC, DIMENSION(jpmax_harmo) :: Wave !:
44
45	REAL(wp) :: sh_T, sh_s, sh_h, sh_p, sh_p1 ! astronomic angles
46	REAL(wp) :: sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R !
47	REAL(wp) :: sh_I, sh_x1ra, sh_N !
48
49
50	!JT origin angles
51	REAL(wp) :: sh_T_o, sh_s_o, sh_h_o, sh_p_o, sh_p1_o, sh_N_o ! astronomic angles
52	!REAL(wp) :: sh_xi_o, sh_nu_o, sh_nuprim_o, sh_nusec_o, sh_R_o !
53	!REAL(wp) :: sh_I_o, sh_x1ra_o !
54	!JT origin angles
55
56	!!JT
57	INTEGER(KIND=8) :: days_since_origin
58	LOGICAL :: ln_astro_verbose
59	!LOGICAL :: ln_tide_360_cal
60	!LOGICAL :: ln_tide_drift_time_cont_manual
61	LOGICAL :: ln_tide_drift ! Do we want to run with "drifting" tides? (Namelist)
62	LOGICAL :: ln_tide_compress ! Do we want to run with "compressed" tides? (Namelist)
63	INTEGER :: nn_tide_orig_yr,nn_tide_orig_mn,nn_tide_orig_dy !JT
64
65	!!JT
66
67	!!----------------------------------------------------------------------
68	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
69	!! $Id$
70	!! Software governed by the CeCILL license (see ./LICENSE)
71	!!----------------------------------------------------------------------
72	CONTAINS
73
74	SUBROUTINE tide_init_Wave
75	# include "tide.h90"
76	END SUBROUTINE tide_init_Wave
77
78
79	SUBROUTINE tide_init_calendar_options
80
81	INTEGER :: ios
82
83
84	ln_tide_drift = .FALSE.
85	ln_tide_compress = .FALSE.
86
87	NAMELIST/nam_tides360/ ln_tide_drift,ln_tide_compress,ln_astro_verbose,&
88	& nn_tide_orig_yr,nn_tide_orig_mn,nn_tide_orig_dy
89
90	! read in Namelist.
91	!!----------------------------------------------------------------------
92	!
93	REWIND ( numnam_ref ) ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics
94	READ ( numnam_ref, nam_tides360, IOSTAT=ios, ERR= 901 )
95	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_tides360 in reference namelist' )
96
97	REWIND( numnam_cfg ) ! Namelist nam_diatmb in configuration namelist TMB diagnostics
98	READ ( numnam_cfg, nam_tides360, IOSTAT = ios, ERR = 902 )
99	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_tides360 in configuration namelist' )
100	IF(lwm) WRITE ( numond, nam_tides360 )
101
102
103	IF( lwp ) THEN
104	WRITE(numout,*) " "
105	WRITE(numout,*) "tide_harmo: nam_tides360 - 360 day tides "
106	WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
107	WRITE(numout,*) " tides360: allow tides to drift through year: ln_tide_drift = ",ln_tide_drift
108	WRITE(numout,*) " tides360: Compress tides, so around a 360 day year: ln_tide_compress = ",ln_tide_compress
109	WRITE(numout,*) " tides360: USE ln_tide_compress WITH CARE. INCOMPLETE."
110	WRITE(numout,*) " tides360: Increase output verbosity: ln_astro_verbose = ",ln_astro_verbose
111	!WRITE(numout,*) " tides360: Calculate time between origin and gregorian and 360 manually: ln_tide_drift_time_cont_manual = ",ln_tide_drift_time_cont_manual
112	WRITE(numout,*) " tides360: 360 day origin date year: nn_tide_orig_yr = ",nn_tide_orig_yr
113	WRITE(numout,*) " tides360: 360 day origin date month: nn_tide_orig_mn = ",nn_tide_orig_mn
114	WRITE(numout,*) " tides360: 360 day origin date day: nn_tide_orig_dy = ",nn_tide_orig_dy
115	WRITE(numout,*) " "
116	ENDIF
117
118
119	IF( nleapy == 30 ) THEN
120	IF ( ln_tide_drift .AND. ln_tide_compress ) THEN
121	CALL ctl_stop( 'tide_harmo: nam_tides360: if 360 day calendar ln_tide_drift and ln_tide_compress cannot be true' )
122	ENDIF
123
124
125	IF ( ln_tide_drift ) THEN
126	WRITE(numout,*) " tides360: Tides continuous so equinoctal tides drift through the year,"
127	WRITE(numout,*) " as the S2-K2 beating occurs 5 days later every year."
128	ENDIF
129
130	IF ( ln_tide_compress ) THEN
131	WRITE(numout,*) " tides360: The Tropical Year (and so some tidal periods) are compressed,"
132	WRITE(numout,*) " so the tides repeat with an annual cycle, so the "
133	WRITE(numout,*) " the S2-K2 beating is fixed relative to the calendar, but the "
134	WRITE(numout,*) " M2 period varies slightly."
135	WRITE(numout,*) " Use with care, as this requires more work."
136	ENDIF
137
138	IF ( ( .NOT. ln_tide_drift ) .AND. ( .NOT. ln_tide_compress ) ) THEN
139	WRITE(numout,*) " tides360: Use the default NEMO tide code, where the tides are reset "
140	WRITE(numout,*) " at the beginning of each month, leading to a slight discontinuity"
141	WRITE(numout,*) " in the tides, and making tidal analysis difficult."
142	ENDIF
143
144	ELSE
145	WRITE(numout,*) " tides360: Gregorian calendar so using standard tides"
146	ENDIF
147
148	!IF ( ln_tide_compress )
149	CALL astronomic_angle_origin
150
151	END SUBROUTINE tide_init_calendar_options
152
153
154	SUBROUTINE tide_harmo( pomega, pvt, put , pcor, ktide ,kc)
155
156	!! Externally called by sbctide.F90/sbc_tide
157	!! Externally named: omega_tide, v0tide, utide, ftide, ntide, nb_harmo
158	!!----------------------------------------------------------------------
159	!!----------------------------------------------------------------------
160	INTEGER , DIMENSION(kc), INTENT(in ) :: ktide ! Indice of tidal constituents
161	INTEGER , INTENT(in ) :: kc ! Total number of tidal constituents
162	REAL(wp), DIMENSION(kc), INTENT(out) :: pomega ! pulsation in radians/s
163	REAL(wp), DIMENSION(kc), INTENT(out) :: pvt, put, pcor !
164	!!----------------------------------------------------------------------
165	!
166
167	! INTEGER :: ios
168
169
170	! ln_tide_drift = .FALSE.
171	! ln_tide_compress = .FALSE.
172
173	! NAMELIST/nam_tides360/ ln_tide_drift,ln_tide_compress,ln_astro_verbose,&
174	! & nn_tide_orig_yr,nn_tide_orig_mn,nn_tide_orig_dy
175
176	! ! read in Namelist.
177	! !!----------------------------------------------------------------------
178	! !
179	! REWIND ( numnam_ref ) ! Read Namelist nam_diatmb in referdiatmbence namelist : TMB diagnostics
180	! READ ( numnam_ref, nam_tides360, IOSTAT=ios, ERR= 901 )
181	!901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'nam_tides360 in reference namelist' )
182
183	! REWIND( numnam_cfg ) ! Namelist nam_diatmb in configuration namelist TMB diagnostics
184	! READ ( numnam_cfg, nam_tides360, IOSTAT = ios, ERR = 902 )
185	!902 IF( ios > 0 ) CALL ctl_nam ( ios , 'nam_tides360 in configuration namelist' )
186	! IF(lwm) WRITE ( numond, nam_tides360 )
187
188
189	! IF( lwp ) THEN
190	! WRITE(numout,*) " "
191	! WRITE(numout,*) "tide_harmo: nam_tides360 - 360 day tides "
192	! WRITE(numout,*) "~~~~~~~~~~~~~~~~~~~~~"
193	! WRITE(numout,*) " tides360: allow tides to drift through year: ln_tide_drift = ",ln_tide_drift
194	! WRITE(numout,*) " tides360: Compress tides, so around a 360 day year: ln_tide_compress = ",ln_tide_compress
195	! WRITE(numout,*) " tides360: USE ln_tide_compress WITH CARE. INCOMPLETE."
196	! WRITE(numout,*) " tides360: Increase output verbosity: ln_astro_verbose = ",ln_astro_verbose
197	! !WRITE(numout,*) " tides360: Calculate time between origin and gregorian and 360 manually: ln_tide_drift_time_cont_manual = ",ln_tide_drift_time_cont_manual
198	! WRITE(numout,*) " tides360: 360 day origin date year: nn_tide_orig_yr = ",nn_tide_orig_yr
199	! WRITE(numout,*) " tides360: 360 day origin date month: nn_tide_orig_mn = ",nn_tide_orig_mn
200	! WRITE(numout,*) " tides360: 360 day origin date day: nn_tide_orig_dy = ",nn_tide_orig_dy
201	! WRITE(numout,*) " "
202	! ENDIF
203
204	!
205	! IF( nleapy == 30 ) THEN
206	! IF ( ln_tide_drift .AND. ln_tide_compress ) THEN
207	! CALL ctl_stop( 'tide_harmo: nam_tides360: if 360 day calendar ln_tide_drift and ln_tide_compress cannot be true' )
208	! ENDIF
209	!
210
211	! IF ( ln_tide_drift ) THEN
212	! WRITE(numout,*) " tides360: Tides continuous so equinoctal tides drift through the year,"
213	! WRITE(numout,*) " as the S2-K2 beating occurs 5 days later every year."
214	! ENDIF
215
216	! IF ( ln_tide_compress ) THEN
217	! WRITE(numout,*) " tides360: The Tropical Year (and so some tidal periods) are compressed,"
218	! WRITE(numout,*) " so the tides repeat with an annual cycle, so the "
219	! WRITE(numout,*) " the S2-K2 beating is fixed relative to the calendar, but the "
220	! WRITE(numout,*) " M2 period varies slightly."
221	! WRITE(numout,*) " Use with care, as this requires more work."
222	! ENDIF
223
224	! IF ( ( .NOT. ln_tide_drift ) .AND. ( .NOT. ln_tide_compress ) ) THEN
225	! WRITE(numout,*) " tides360: Use the default NEMO tide code, where the tides are reset "
226	! WRITE(numout,*) " at the beginning of each month, leading to a slight discontinuity"
227	! WRITE(numout,*) " in the tides, and making tidal analysis difficult."
228	! ENDIF
229
230	! ELSE
231	! WRITE(numout,*) " tides360: Gregorian calendar so using standard tides"
232	! ENDIF
233
234	CALL astronomic_angle
235	CALL tide_pulse( pomega, ktide ,kc )
236	CALL tide_vuf ( pvt, put, pcor, ktide ,kc )
237
238	!
239	END SUBROUTINE tide_harmo
240
241
242	SUBROUTINE astronomic_angle
243	!!----------------------------------------------------------------------
244	!! tj is time elapsed since 1st January 1900, 0 hour, counted in julian
245	!! century (e.g. time in days divide by 36525)
246	!!----------------------------------------------------------------------
247	REAL(wp) :: cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2
248	REAL(wp) :: zqy , zsy, zday, zdj, zhfrac
249
250
251	! JT
252	! Tides are added as boundary conditions, and as tidal potential.
253	!
254	! For the Boundaries, the complex tide amplitudes are give for each point.
255	! This gives the amplitude and the phase for each consititent.
256	! The tidal frequency is calculated from Wave in tide.h90 via tide_pulse below.
257	! A the start(ish) of everyday, astronomic_angle is called via tide_harmo
258	! from SBC/sbctide.F90.
259	! The astronomic_angle specifies the location of the moon and sun etc at the given
260	! model time. these are used to update the tidal phase.
261	!
262	! In the bdytides.F90 the function bdy_dta_tides calls
263	! tide_init_elevation and tide_init_velocities (also in bdytides.F90)
264	! every day. This uses the values from astro angles to update the phase of the
265	! tidal constiuents as read in from the boundary files.
266	!
267	! The tidal potential in (re)initialised every day in sbstide.F90 in the function
268	! tide_init_potential. This uses the values from astro angles:
269	! (v0tide + utide) and produces amp_pot and phi_pot.
270	! These are then used in SBC/updtide.F90 (every timestep?) to set pot_astro.
271	!
272	! Both SBC/sbctide.F90 and bdy_dta_tides calculate zoff+z_arg which is the number of seconds since the beginning of the day.
273	! the tidal phases are then corrected for this reset with the v0tide parameter, calucate by tide_vuf.
274	! nodal correction is much smaller, with ftide (which affects the amplitude), and utide (which affects the phase).
275	!
276	!
277	! As the phase of the tidal constituents for both the boundaries and the tidal potential
278	! are adjusted by the astronomic_angle, we can adapt this one module to adapt the tides
279	! for 360 day calendars.
280	!
281	! There are different approaches to tides in a 360 day calendar.
282	! 1) (current), the tides are effectively reset to the first of the month.
283	! therefore skip 31st's and repeat 29th and 30th of Feb
284	! its happy with extra days of the month (doesn't crash for 30th Feb)
285	! Tide is anchored to correct part of the year, but extra/missing days
286	! are unrealistic, add noise to the system, and make least square tidal analysis fail
287	!
288	! 2) Start the tides at the begining of the run and then let run continuously.
289	! The tides drift throughout the year, so the equinox's are not at the correct part of the year.
290	!
291	! This is the approach set up below
292	!
293	! 2b) Adapt the equations to use decimal years (they sort of do, as they use day of year)
294	! This would make the counting forward and backward from the origin easier
295	! (the final step (going from DOY to mon and yr) would be removed)
296	!
297	! 4) Adapt the equations that affect the location of the moon and tides.
298	! This very likely to be a hugely complex job, that would affect the amphidromic systems,
299	! As you're likely to need to change many/all of the tidal constants. this is then likely
300	! to change the tidal frequencies, and so the the tidal wave speed, and hence the amphidromes,
301	! co-tide and co-phase lines.
302	!
303	! This approach is not followed
304	!
305	!
306	! To make the tide continueous for 360 and 365.25 day calendars,
307	! firstly, I make temporary working yr/mn/day integers.
308	! for the gregorian calendar these are simply set to nyear, nmonth and nday.
309	!
310	! For a 360 day calendar, I count the days from 1900/1/1 to the current day
311	! according to the the 360 day calendar.
312	! I then count forward that many days according to the gregorian calendar.
313	! therefore every 30 day month of the 360d model run, the tides move forward 30 days.
314
315
316
317	! Questions:
318	! Are the better ways of adding offets to dates in Nemo/Fortran? i.e. python timedelta from datetime
319	! is there a leap year function in NEMO/Fortran?
320	! Not sure if the algorithm is very stable for different origins.
321	! nleap is corrected for 1900 not being a leap year. Needs updated for a different origin year
322	! Does it work if its not starting on a leap year?
323	! Does it work if its after the 28th Feb?
324	! Does it work if the origin is after the start date of the run?
325	! When adjusting the DOY and Y for the number of leap years, what happens if its more than 365 leap years?
326	!
327	! h mean solar Longitude and s mean lunar Longitude and are functions of zhfrac, zday and zsy,
328	! but the coeffiencents are not 1/86400:1:365
329	! zday is the DOY corrected for the number of leap years since 1900.
330	! So can run from 20 to 385. this wouldn't matter if the coefficients were 1/86400:1:365
331	! should zsy and zday be updated so zday is between e.g. 1 and 365??
332	!
333	! What are the impacts on the NWS if the tide drifts?
334	! What is the impact on the NWS if the tide repeats/skips days?
335	! can this make the model go unstable?
336
337
338
339	! New variables defined for new code
340	INTEGER :: yr_org,mn_org,dy_org !JT
341	REAL(wp) :: sec_grg !JT
342	INTEGER :: yr_grg,mn_grg,dy_grg !JT
343	INTEGER :: yr_360,mn_360,dy_360 !JT
344	INTEGER :: yr_wrk,mn_wrk,dy_wrk !JT
345
346
347
348	LOGICAL :: ln_tide_drift_time_cont ! Do we correct for a continueous time
349
350
351	! INTEGER(KIND=8) :: days_since_origin ! added to module
352	INTEGER :: init_yr, day_in_init_yr,nleap,init_doy
353	INTEGER :: init_doy_inc_l,yg_is_leap_mod,doy_grg
354	INTEGER,DIMENSION(12) :: idayt, idays
355	!INTEGER :: ji
356
357	!REAL(wp) :: fjulday_org !: current julian day
358	! REAL(wp) :: days_since_origin_ymds2ju
359	INTEGER(KIND=8) :: days_since_origin_ymds2ju_int
360
361
362
363	REAL(wp) :: jul_org_greg,jul_org_360,jul_pres_360
364
365	DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
366
367
368	! Currently hardcode the verbosity and the of the code
369	! how to I read the calendar type
370	!ln_tide_360_cal = .TRUE.
371	!IF ( nleapy == 30) ln_tide_360_cal = .TRUE.
372
373
374
375	!! Nameslist values
376
377
378	IF (ln_astro_verbose .AND. lwp) THEN
379	WRITE(numout,*) 'astro '
380	WRITE(numout,*) 'astro -------------------------------------------------'
381	ENDIF
382
383
384	ln_tide_drift_time_cont = .FALSE. ! the same the original code
385
386	IF( nleapy == 30 ) THEN
387	IF ( ln_tide_drift ) THEN
388	ln_tide_drift_time_cont = .TRUE.
389	ENDIF
390	IF ( ln_tide_compress ) THEN
391	! ##################################################################
392	! ##################################################################
393	! ##################################################################
394	! #####
395	! ##### For the 360 day tide constituents,
396	! ##### We only use days_since_origin for v0tide in tide_vuf.
397	! #####
398	! ##### To use the correct tide nodal correction (utide)
399	! ##### (which is a small ajustment)
400	! ##### use keep that linked to the gregorian dates.
401	! #####
402	! #####
403	! #####
404	! ##### Therefore, we keep yr_wrk, mn_wrk, dy_wrk to equal
405	! ##### nyear, nmonth, nday
406	! #####
407	! ##################################################################
408	! ##################################################################
409	! ##################################################################
410
411	ln_tide_drift_time_cont = .FALSE.
412
413	! ##################################################################
414	! ##################################################################
415	! ##################################################################
416	! #####
417	! ##### NEMO4.0.4
418	! ##### BUT!!! need to calc days_since_origin, so need to
419	! ##### set ln_tide_drift_time_cont too true, then reset
420	! ##### yr_wrk, mn_wrk, dy_wrk to equal nyear, nmonth, nday
421	! #####
422	! #####
423	! ##################################################################
424	! ##################################################################
425	! ##################################################################
426	ln_tide_drift_time_cont = .TRUE.
427	ENDIF
428
429	ELSE
430	ln_tide_drift_time_cont = .FALSE.
431	ENDIF
432
433
434	IF (ln_astro_verbose .AND. lwp) THEN
435	WRITE(numout,*) 'astro ln_tide_drift_time_cont = ',ln_tide_drift_time_cont
436	ENDIF
437
438	!IF( ln_tide_360_cal ) THEN
439	!IF( nleapy == 30 ) THEN
440	IF( ln_tide_drift_time_cont ) THEN
441
442
443	! clear and set current dates.
444	yr_360 = nyear
445	mn_360 = nmonth
446	dy_360 = nday
447
448
449	yr_grg = 0
450	mn_grg = 0
451	dy_grg = 0
452
453	yr_wrk = 0
454	mn_wrk = 0
455	dy_wrk = 0
456
457	! Set the origin in the name list
458
459	yr_org = nn_tide_orig_yr
460	mn_org = nn_tide_orig_mn
461	dy_org = nn_tide_orig_dy
462
463
464	!IF (ln_tide_drift_time_cont_manual) THEN
465
466
467
468	! IF (ln_astro_verbose .AND. lwp) THEN
469	! WRITE(numout,) 'astro: yr_360,yr_org,((yr_360-yr_org)360)', yr_360,yr_org,((yr_360-yr_org)*360)
470	! WRITE(numout,) 'astro: mn_360,mn_org,((mn_360-mn_org)30)', mn_360,mn_org,((mn_360-mn_org)*30)
471	! WRITE(numout,*) 'astro: dy_360,dy_org,(dy_360-dy_org)', dy_360,dy_org,(dy_360-dy_org)
472	! ENDIF
473	!
474	! ! how many days from 1900 in the 360 day calendar
475	! days_since_origin = ((yr_360-yr_org)360) + ((mn_360-mn_org)30) + (dy_360-dy_org)
476	!
477	! ! first guess of what year this would be for the same numbers of days from 1/1/1900 in a gregorian calendar
478	! init_yr = yr_org + days_since_origin/365
479	!
480	! ! was the initial estimated year a leap year? how many days in this year?
481	! day_in_init_yr = 365
482	! if (MOD(init_yr,4) == 0) day_in_init_yr = 366
483	!
484	!
485	!
486	! !CALL ymds2ju_JT (yr_org, mn_org, dy_org, 0.0, fjulday_org,360.)
487	!
488	! !IF (ln_astro_verbose) THEN
489	! ! IF(lwp) THEN
490	! ! WRITE(numout,*) 'astro: ymds2ju_JT yr_org, mn_org, dy_org,fjulday_org', yr_org, mn_org, dy_org,fjulday_org
491	! ! ENDIF
492	! !ENDIF
493	!
494	!
495	! !CALL ymds2ju( yr_org, mn_org, dy_org, 0.0, fjulday_org ) ! we assume that we start run at 00:00
496	! !IF( ABS(fjulday_org - REAL(NINT(fjulday_org),wp)) < 0.1 / rday ) fjulday_org = REAL(NINT(fjulday_org),wp) ! avoid truncation error
497	! !fjulday_org = fjulday_org + 1. ! move back to the day at nit000 (and not at nit000 - 1)
498	!
499	! !days_since_origin_ymds2ju_int = AINT(fjulday - fjulday_org)
500	!
501	! IF (ln_astro_verbose .AND. lwp) THEN
502	! WRITE(numout,*) 'astro: days_since_origin,init_yr,day_in_init_yr', days_since_origin,init_yr,day_in_init_yr
503	! !WRITE(numout,*) 'astro: fjulday_org', fjulday_org
504	! !WRITE(numout,*) 'astro: fjulday', fjulday
505	! !WRITE(numout,*) 'astro: fjulday - fjulday_org', fjulday - fjulday_org
506	! !WRITE(numout,*) 'astro: days_since_origin_ymds2ju_int', days_since_origin_ymds2ju_int
507	! ENDIF
508	!
509	!
510	! ! how many leap years since the origin.
511	! nleap = (yr_360-yr_org)/4 - 1 !1900 is not a leap year
512	!
513	! ! initial estimate of the day of year
514	! init_doy = MOD(days_since_origin,365)
515	!
516	! ! correct the initial estimate for the DOY for the number of leap days since the origin
517	! init_doy_inc_l = init_doy - nleap
518	!
519	!
520	! IF (ln_astro_verbose .AND. lwp) THEN
521	! WRITE(numout,*) 'astro: nleap,init_doy,init_doy_inc_l',nleap,init_doy,init_doy_inc_l
522	! ENDIF
523	!
524	!
525	! ! The number of leap days could pull the DOY before 0.
526	! ! in which case decrement the year, and reset the DOY.
527	! ! of the origin is 365 leap years ago, and initial DOY could be adjusted by more than one year..
528	! ! Unlikely to be a prob, but need to remember if planning very long control runs. Need to think about this.
529	!
530	! IF (init_doy_inc_l .LT. 0) THEN
531	! init_doy_inc_l = init_doy_inc_l+365
532	! init_yr = init_yr - 1
533	! IF (MOD(init_yr, 4) == 0 ) THEN
534	! init_doy_inc_l = init_doy_inc_l + 1
535	! ENDIF
536	! ENDIF
537	!
538	!
539	! ! This gives the year and the day of year in the gregorian calendar
540	! yr_grg = init_yr
541	! doy_grg = init_doy_inc_l
542	! yg_is_leap_mod = MOD(yr_grg, 4)
543	!
544	! IF (ln_astro_verbose .AND. lwp) THEN
545	! WRITE(numout,*) 'astro: yr_grg,doy_grg,yg_is_leap_mod',yr_grg,doy_grg,yg_is_leap_mod
546	! ENDIF
547	!
548	!
549	! ! Convert from day of year to month and day in the gregorian calendar.
550	! ! dayjul code adapted
551	! ! this perhaps should be a function, but not sure how to write one
552	! ! there may be this code functionality elsewhere in NEMO
553	! !!----------------------------------------------------------------------
554	!
555	!
556	! ! what is the DOY of the first day of the month for each month.
557	! ! correct for leap years.
558	!
559	! idays(1) = 0.
560	! idays(2) = 31.
561	! inc = 0.
562	! IF( yg_is_leap_mod == 0.) inc = 1.
563	!
564	! DO ji = 3, 12
565	! idays(ji)=idayt(ji)+inc
566	! END DO
567	!
568	! ! cycle through the months.
569	! ! if the DOY is greater than the DOY of the first Day of Month
570	! ! Note the month. Calculate day of month by subtraction.
571	! ! Once beyond the correct month, the if statement won't be true, so wont calculate.
572	!
573	! DO ji = 1, 12
574	! IF ( doy_grg .GE. idays(ji) ) THEN
575	! mn_grg = ji
576	! dy_grg = doy_grg-idays(ji) +1
577	! ENDIF
578	! END DO
579	!
580	!
581	!
582	!
583	!
584	! IF(ln_astro_verbose .AND. lwp) THEN
585	! WRITE(numout,*) 'astro: mn_grg,dy_grg',mn_grg,dy_grg
586	! WRITE(numout,*) ' '
587	! WRITE(numout,*) 'tide_mod_astro_ang 360_corr : yr_360,mn_360,dy_360,yr_grg,mn_grg,dy_grg,doy_grg =',yr_360,mn_360,dy_360,yr_grg,mn_grg,dy_grg,doy_grg
588	!
589	! WRITE(numout,*) ' '
590	! ENDIF
591	!
592	!
593	!
594	! IF (ln_astro_verbose .AND. lwp) WRITE(numout,*) 'tide_mod_astro_ang_meth_1,',yr_grg, mn_grg, dy_grg
595
596
597	!ELSE ! ln_tide_drift_time_cont_manual
598
599
600	! number of days since 15th October 1582, for namelist origin, in both calendars, and for current model day.
601
602	CALL ymds2ju_JT( yr_org,mn_org,dy_org, 0. ,jul_org_greg,365.24 )
603	CALL ymds2ju_JT( yr_org,mn_org,dy_org, 0. ,jul_org_360,360. )
604	CALL ymds2ju_JT( yr_360,mn_360,dy_360, 0. ,jul_pres_360,360. )
605
606	! Calculate the days since the origin: days_since_origin_ymds2ju_int
607	! How many days between the current day, and the origin, in the 360 day calendar.
608	days_since_origin_ymds2ju_int = jul_pres_360 - jul_org_360
609
610	IF (ln_astro_verbose .AND. lwp) THEN
611	WRITE(numout,*) 'tide_mod_astro_ang 360_corr : jul_org_360,jul_pres_360,jul_pres_360 - jul_org_360 =',jul_org_360,jul_pres_360,jul_pres_360 - jul_org_360
612	WRITE(numout,*) 'tide_mod_astro_ang 360_corr : days_since_origin_ymds2ju_int, days_since_origin_ymds2ju_int mod 360 =',days_since_origin_ymds2ju_int,MOD( days_since_origin_ymds2ju_int ,360 )
613	WRITE(numout,*) 'tide_mod_astro_ang 360_corr : yr_org,mn_org,dy_org, jul_org_greg =',yr_org,mn_org,dy_org, jul_org_greg
614	ENDIF
615
616	!add days_since_origin_ymds2ju_int days to the origin in the gregorian calendar.
617	CALL ju2ymds_JT( days_since_origin_ymds2ju_int + jul_org_greg, yr_grg, mn_grg, dy_grg, sec_grg,365.24 )
618
619	IF (ln_astro_verbose .AND. lwp) THEN
620	WRITE(numout,*) 'tide_mod_astro_ang 360_corr : yr_grg, mn_grg, dy_grg =',yr_grg, mn_grg, dy_grg
621	WRITE(numout,*) 'tide_mod_astro_ang 360_corr : yr_360, mn_360, dy_360 =',yr_360, mn_360, dy_360
622	WRITE(numout,*) 'tide_mod_astro_ang 360_corr : yr_org, mn_org, dy_org =',yr_org, mn_org, dy_org
623	ENDIF
624
625
626
627
628	IF (ln_astro_verbose .AND. lwp) WRITE(numout,*) 'tide_mod_astro_ang_meth_2,',yr_grg, mn_grg, dy_grg
629
630	!ENDIF !ln_tide_drift_time_cont_manual
631
632	! for 360 calendars, work with the pseudo gregorian dates
633	yr_wrk = yr_grg
634	mn_wrk = mn_grg
635	dy_wrk = dy_grg
636
637	days_since_origin = days_since_origin_ymds2ju_int
638
639
640	!IF (ln_tide_compress) THEN
641	! yr_wrk = nyear
642	! mn_wrk = nmonth
643	! dy_wrk = nday
644	!ENDIF
645
646	ELSE
647
648	! for gregorian calendars, work with the model gregorian dates
649	yr_wrk = nyear
650	mn_wrk = nmonth
651	dy_wrk = nday
652
653	ENDIF
654
655	! continue with original code, using working year, month and day.
656
657	!
658	zqy = AINT( (yr_wrk-1901.)/4. ) ! leap years since 1901
659	zsy = yr_wrk - 1900. ! years since 1900
660	!
661	zdj = dayjul( yr_wrk, mn_wrk, dy_wrk ) ! day number of year
662	zday = zdj + zqy - 1. ! day number of year + No of leap yrs
663	! i.e. what would doy if every year = 365 day??
664	!
665	zhfrac = nsec_day / 3600. ! The seconds of the day/3600
666
667
668	!
669	!----------------------------------------------------------------------
670	! Sh_n Longitude of ascending lunar node
671	!----------------------------------------------------------------------
672	sh_N=(259.1560564-19.328185764zsy-.0529539336zday-.0022064139zhfrac)rad
673	!----------------------------------------------------------------------
674	! T mean solar angle (Greenwhich time)
675	!----------------------------------------------------------------------
676	sh_T=(180.+zhfrac(360./24.))rad
677	!----------------------------------------------------------------------
678	! h mean solar Longitude
679	!----------------------------------------------------------------------
680	sh_h=(280.1895014-.238724988zsy+.9856473288zday+.0410686387zhfrac)rad
681	!----------------------------------------------------------------------
682	! s mean lunar Longitude
683	!----------------------------------------------------------------------
684	sh_s=(277.0256206+129.38482032zsy+13.176396768zday+.549016532zhfrac)rad
685	!----------------------------------------------------------------------
686	! p1 Longitude of solar perigee
687	!----------------------------------------------------------------------
688	sh_p1=(281.2208569+.01717836zsy+.000047064zday+.000001961zhfrac)rad
689	!----------------------------------------------------------------------
690	! p Longitude of lunar perigee
691	!----------------------------------------------------------------------
692	sh_p=(334.3837214+40.66246584zsy+.111404016zday+.004641834zhfrac)rad
693
694
695
696	IF(ln_astro_verbose .AND. lwp) THEN
697	WRITE(numout,*)
698	WRITE(numout,*) 'tide_mod_astro_ang : yr_wrk,mn_wrk,dy_wrk=',yr_wrk,mn_wrk,dy_wrk
699	WRITE(numout,*) 'tide_mod_astro_ang : nyear, nmonth, nday,nsec_day=',nyear, nmonth, nday,nsec_day
700	WRITE(numout,*) 'tide_mod_astro_ang : sh_N,sh_T,sh_h,sh_s,sh_p1,sh_p=', sh_N,sh_T,sh_h,sh_s,sh_p1,sh_p
701	WRITE(numout,*) 'tide_mod_astro_ang : zsy,zday,zhfrac,rad=', zsy,zday,zhfrac,rad
702	WRITE(numout,*) 'tide_mod_astro_ang : zqy ,zdj,yr_wrk, mn_wrk, dy_wrk =', zqy ,zdj,yr_wrk, mn_wrk, dy_wrk
703	WRITE(numout,*) '~~~~~~~~~~~~~~ '
704	ENDIF
705
706
707
708	sh_N = MOD( sh_N ,2*rpi )
709	sh_s = MOD( sh_s ,2*rpi )
710	sh_h = MOD( sh_h, 2*rpi )
711	sh_p = MOD( sh_p, 2*rpi )
712	sh_p1= MOD( sh_p1,2*rpi )
713
714	cosI = 0.913694997 -0.035692561 *cos(sh_N)
715
716	sh_I = ACOS( cosI )
717
718	sin2I = sin(sh_I)
719	sh_tgn2 = tan(sh_N/2.0)
720
721	at1=atan(1.01883*sh_tgn2)
722	at2=atan(0.64412*sh_tgn2)
723
724	sh_xi=-at1-at2+sh_N
725
726	IF( sh_N > rpi ) sh_xi=sh_xi-2.0*rpi
727
728	sh_nu = at1 - at2
729
730	!----------------------------------------------------------------------
731	! For constituents l2 k1 k2
732	!----------------------------------------------------------------------
733
734	tgI2 = tan(sh_I/2.0)
735	P1 = sh_p-sh_xi
736
737	t2 = tgI2*tgI2
738	t4 = t2*t2
739	sh_x1ra = sqrt( 1.0-12.0t2cos(2.0P1)+36.0t4 )
740
741	p = sin(2.0*P1)
742	q = 1.0/(6.0t2)-cos(2.0P1)
743	sh_R = atan(p/q)
744
745	p = sin(2.0sh_I)sin(sh_nu)
746	q = sin(2.0sh_I)cos(sh_nu)+0.3347
747	sh_nuprim = atan(p/q)
748
749	s2 = sin(sh_I)*sin(sh_I)
750	p = s2sin(2.0sh_nu)
751	q = s2cos(2.0sh_nu)+0.0727
752	sh_nusec = 0.5*atan(p/q)
753	!
754	END SUBROUTINE astronomic_angle
755
756
757
758
759
760
761	SUBROUTINE astronomic_angle_origin
762	!!----------------------------------------------------------------------
763	!! tj is time elapsed since 1st January 1900, 0 hour, counted in julian
764	!! century (e.g. time in days divide by 36525)
765	!!----------------------------------------------------------------------
766	REAL(wp) :: cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2
767	REAL(wp) :: zqy , zsy, zday, zdj, zhfrac
768
769
770
771
772
773	! New variables defined for new code
774	INTEGER :: yr_wrk,mn_wrk,dy_wrk !JT
775
776
777
778	! for gregorian calendars, work with the model gregorian dates
779	yr_wrk = nn_tide_orig_yr
780	mn_wrk = nn_tide_orig_mn
781	dy_wrk = nn_tide_orig_dy
782
783
784	!
785	zqy = AINT( (yr_wrk-1901.)/4. ) ! leap years since 1901
786	zsy = yr_wrk - 1900. ! years since 1900
787	!
788	zdj = dayjul( yr_wrk, mn_wrk, dy_wrk ) ! day number of year
789	zday = zdj + zqy - 1. ! day number of year + No of leap yrs
790	! i.e. what would doy if every year = 365 day??
791	!
792	zhfrac = nsec_day / 3600. ! The seconds of the day/3600
793
794
795	!
796	!----------------------------------------------------------------------
797	! Sh_n Longitude of ascending lunar node
798	!----------------------------------------------------------------------
799	sh_N_o=(259.1560564-19.328185764zsy-.0529539336zday-.0022064139zhfrac)rad
800	!----------------------------------------------------------------------
801	! T mean solar angle (Greenwhich time)
802	!----------------------------------------------------------------------
803	sh_T_o=(180.+zhfrac(360./24.))rad
804	!----------------------------------------------------------------------
805	! h mean solar Longitude
806	!----------------------------------------------------------------------
807	sh_h_o=(280.1895014-.238724988zsy+.9856473288zday+.0410686387zhfrac)rad
808	!----------------------------------------------------------------------
809	! s mean lunar Longitude
810	!----------------------------------------------------------------------
811	sh_s_o=(277.0256206+129.38482032zsy+13.176396768zday+.549016532zhfrac)rad
812	!----------------------------------------------------------------------
813	! p1 Longitude of solar perigee
814	!----------------------------------------------------------------------
815	sh_p1_o=(281.2208569+.01717836zsy+.000047064zday+.000001961zhfrac)rad
816	!----------------------------------------------------------------------
817	! p Longitude of lunar perigee
818	!----------------------------------------------------------------------
819	sh_p_o=(334.3837214+40.66246584zsy+.111404016zday+.004641834zhfrac)rad
820
821
822
823	IF(ln_astro_verbose .AND. lwp) THEN
824	WRITE(numout,*)
825	WRITE(numout,*) 'tide_mod_astro_ang_orig,yr_wrk,mn_wrk,dy_wrk,nsec_day,=',yr_wrk,mn_wrk,dy_wrk,nsec_day
826	WRITE(numout,*) 'tide_mod_astro_ang_orig,sh_N_o,sh_T_o,sh_h_o,sh_s_o,sh_p1_o,sh_p_o,', sh_N_o,sh_T_o,sh_h_o,sh_s_o,sh_p1_o,sh_p_o
827	WRITE(numout,*) 'tide_mod_astro_ang_orig,zqy ,zdj,zsy,zday,zhfrac,rad,', zqy ,zdj,zsy,zday,zhfrac,rad
828
829	WRITE(numout,*) '~~~~~~~~~~~~~~ '
830	ENDIF
831
832
833
834	sh_N_o = MOD( sh_N_o ,2*rpi )
835	sh_s_o = MOD( sh_s_o ,2*rpi )
836	sh_h_o = MOD( sh_h_o, 2*rpi )
837	sh_p_o = MOD( sh_p_o, 2*rpi )
838	sh_p1_o= MOD( sh_p1_o,2*rpi )
839
840
841	! cosI = 0.913694997 -0.035692561 *cos(sh_N_o)
842	!
843	!
844	!
845	!! REAL(wp) :: cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2
846	!! REAL(wp) :: zqy , zsy, zday, zdj, zhfrac
847	!
848	!
849	! sh_I_o = ACOS( cosI )
850	!
851	! sin2I = sin(sh_I_o)
852	! sh_tgn2 = tan(sh_N_o/2.0)
853	!
854	! at1=atan(1.01883*sh_tgn2)
855	! at2=atan(0.64412*sh_tgn2)
856	!
857	! sh_xi_o=-at1-at2+sh_N
858	!
859	! IF( sh_N_o > rpi ) sh_xi_o=sh_xi_o-2.0*rpi
860	!
861	! sh_nu_o = at1 - at2
862	!
863	! !----------------------------------------------------------------------
864	! ! For constituents l2 k1 k2
865	! !----------------------------------------------------------------------
866	!
867	! tgI2 = tan(sh_I_o/2.0)
868	! P1 = sh_p_o-sh_xi_o
869	!
870	! t2 = tgI2*tgI2
871	! t4 = t2*t2
872	! sh_x1ra_o = sqrt( 1.0-12.0t2cos(2.0P1)+36.0t4 )
873	!
874	! p = sin(2.0*P1)
875	! q = 1.0/(6.0t2)-cos(2.0P1)
876	! sh_R = atan(p/q)
877	!
878	! p = sin(2.0sh_I)sin(sh_nu)
879	! q = sin(2.0sh_I)cos(sh_nu)+0.3347
880	! sh_nuprim_o = atan(p/q)
881	!
882	! s2 = sin(sh_I_o)*sin(sh_I_o)
883	! p = s2sin(2.0sh_nu_o)
884	! q = s2cos(2.0sh_nu_o)+0.0727
885	! sh_nusec_o = 0.5*atan(p/q)
886
887
888
889	!
890	END SUBROUTINE astronomic_angle_origin
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906	SUBROUTINE tide_pulse( pomega, ktide ,kc )
907	!!----------------------------------------------------------------------
908	!! * ROUTINE tide_pulse *
909	!!
910	!! ** Purpose : Compute tidal frequencies
911	!!----------------------------------------------------------------------
912	INTEGER , INTENT(in ) :: kc ! Total number of tidal constituents
913	INTEGER , DIMENSION(kc), INTENT(in ) :: ktide ! Indice of tidal constituents
914	REAL(wp), DIMENSION(kc), INTENT(out) :: pomega ! pulsation in radians/s
915	!
916	INTEGER :: jh
917	REAL(wp) :: zscale
918	REAL(wp) :: zomega_T = 13149000.0_wp ! Mean Solar Day ! degrees/century
919	REAL(wp) :: zomega_s = 481267.892_wp ! Sidereal Month ! degrees/century
920	REAL(wp) :: zomega_h != 36000.76892_wp ! Tropical Year ! degrees/century
921	REAL(wp) :: zomega_p = 4069.0322056_wp ! Moons Perigee ! degrees/century
922	REAL(wp) :: zomega_n = 1934.1423972_wp ! Regression of Lunar Nodes ! degrees/century
923	REAL(wp) :: zomega_p1= 1.719175_wp ! Perihelion ! degrees/century
924	!!----------------------------------------------------------------------
925
926	zomega_h = 36000.76892_wp ! Tropical Year ! degrees/century
927	IF (( nleapy == 30 ) .AND. ln_tide_compress) zomega_h = 36525.0_wp ! 360 day Tropical Year ! degrees/century (360 deg/360 days= 1deg/day. cent = 36525
928
929	!
930	zscale = rad / ( 36525._wp * 86400._wp ) ! Convert to radians per second.
931	!
932	DO jh = 1, kc
933	pomega(jh) = ( zomega_T * Wave( ktide(jh) )%nT &
934	& + zomega_s * Wave( ktide(jh) )%ns &
935	& + zomega_h * Wave( ktide(jh) )%nh &
936	& + zomega_p * Wave( ktide(jh) )%np &
937	& + zomega_p1* Wave( ktide(jh) )%np1 ) * zscale
938	END DO
939
940	IF (ln_astro_verbose .AND. lwp) THEN
941
942	WRITE(numout,*) 'astro tide_pulse nleapy:',nleapy
943	WRITE(numout,*) 'astro tide_pulse zomega_h:',zomega_h
944	WRITE(numout,*) 'astro tide_pulse if zomega_h = 36000.76892 for 365.24 day year'
945	WRITE(numout,*) 'astro tide_pulse if zomega_h = 36525.00000 for 360.00 day year'
946
947
948	DO jh = 1, kc
949	WRITE(numout,) 'astro tide_pulse const, pomega, period(hr):',Wave(ktide(jh))%cname_tide, pomega(jh),2rpi/(3600.0_wp*pomega(jh))
950	END DO
951
952	ENDIF
953	!
954	END SUBROUTINE tide_pulse
955
956
957	SUBROUTINE tide_vuf( pvt, put, pcor, ktide ,kc )
958	!!----------------------------------------------------------------------
959	!! * ROUTINE tide_vuf *
960	!!
961	!! ** Purpose : Compute nodal modulation corrections
962	!!
963	!! ** Outputs : vt: Phase of tidal potential relative to Greenwich (radians)
964	!! ut: Phase correction u due to nodal motion (radians)
965	!! ft: Nodal correction factor
966	!!----------------------------------------------------------------------
967	INTEGER , INTENT(in ) :: kc ! Total number of tidal constituents
968	INTEGER , DIMENSION(kc), INTENT(in ) :: ktide ! Indice of tidal constituents
969	REAL(wp), DIMENSION(kc), INTENT(out) :: pvt, put, pcor !
970
971	!
972	INTEGER :: jh ! dummy loop index
973	!!----------------------------------------------------------------------
974
975	!! JT for compress
976	REAL(wp) :: hours_since_origin
977	REAL(wp), DIMENSION(kc) :: pomega ! pulsation in radians/s
978	REAL(wp), DIMENSION(kc) :: freq_per_day, v0linearslope ! ,v0linearintercept ! pulsation in radians/s !offset,cycle_reset,freq,per_hr
979	CHARACTER (len=40) :: tmp_name
980	!! JT for compress
981
982
983	IF( .NOT. ALLOCATED(v0linearintercept) ) ALLOCATE( v0linearintercept(kc) )
984
985
986	IF ( ln_tide_compress) THEN
987
988	CALL tide_pulse( pomega, ktide ,kc )
989
990	DO jh = 1, kc
991	!per_hr(jh) = (2*rpi/pomega(jh))/3600.
992	!freq(jh) = (2*rpi/per_hr(jh))
993	!freq_per_day(jh) = freq(jh)*24
994	freq_per_day(jh) = pomega(jh) * 86400.0_wp
995	!cycle_reset(jh) = mod(hours_since_originfreq(jh),2.rpi)
996	v0linearslope(jh) = - mod ( (-freq_per_day(jh)), (2*rpi) )
997	IF(ln_astro_verbose .AND. lwp) WRITE(numout,) 'astro tide_vuf 1:',jh,kc,ktide(jh),v0linearslope(jh),freq_per_day(jh), pomega(jh),(2rpi/pomega(jh))/3600.! * 86400.0_wp,freq(jh)*24,per_hr(jh),freq(jh)
998	ENDDO
999
1000
1001	! !offset(1) = 0.10789890_wp
1002	! !offset(2) = 1.10897897_wp
1003	! !offset(3) = 2.11005903_wp
1004	! !offset(4) = 0.00000000_wp
1005	! !offset(5) = 3.47632710_wp
1006	! !offset(6) = 0.16751976_wp
1007	! !offset(7) = -0.05503165_wp
1008	! !offset(8) = 0.94604842_wp
1009	! !offset(9) = 6.10534877_wp
1010	! !offset(10) = 0.21579780_wp
1011	! !offset(11) = 0.00000000_wp
1012	! !offset(12) = 0.00000000_wp
1013	! !offset(13) = 0.00000000_wp
1014	! !offset(14) = 0.00000000_wp
1015	! !offset(15) = 3.14159265_wp
1016	! !offset(16) = 0.21833313_wp
1017	! !offset(17) = 5.50043837_wp
1018	! !offset(18) = 2.24841149_wp
1019	! !offset(19) = 0.01800173_wp
1020
1021	! !v0linearintercept(1) = 0.11044027_wp
1022	! !v0linearintercept(2) = 1.11152799_wp
1023	! !v0linearintercept(3) = 2.11261570_wp
1024	! !v0linearintercept(4) = 0.00000000_wp
1025	! !v0linearintercept(5) = 3.49727335_wp
1026	! !v0linearintercept(6) = 0.17784035_wp
1027	! !v0linearintercept(7) = 6.21578523_wp
1028	! !v0linearintercept(8) = 0.93368764_wp
1029	! !v0linearintercept(9) = 6.10534496_wp
1030	! !v0linearintercept(10) = 0.22088055_wp
1031	! !v0linearintercept(11) = 0.00000000_wp
1032	! !v0linearintercept(12) = 0.00000000_wp
1033	! !v0linearintercept(13) = 0.00000000_wp
1034	! !v0linearintercept(14) = 0.00000000_wp
1035	! !v0linearintercept(15) = 3.14159265_wp
1036
1037	! !v0linearintercept(1) = v0linearintercept(1) - 0.000000_wp
1038	! !v0linearintercept(2) = v0linearintercept(2) - 0.000000_wp
1039	! !v0linearintercept(3) = v0linearintercept(3) - 0_wp
1040	! !v0linearintercept(4) = v0linearintercept(4) - 0.165795_wp
1041	! !v0linearintercept(5) = v0linearintercept(5) + 2.821252_wp
1042	! !v0linearintercept(6) = v0linearintercept(6) + 0.479504_wp
1043	! !v0linearintercept(7) = v0linearintercept(7) - 2.175621_wp
1044	! !v0linearintercept(8) = v0linearintercept(8) + 1.900267_wp
1045	! !v0linearintercept(9) = v0linearintercept(9) + 0.107633_wp
1046	! !v0linearintercept(10) = v0linearintercept(10) - 0.000000_wp
1047	! !v0linearintercept(11) = v0linearintercept(11) - 0.000000_wp
1048	! !v0linearintercept(12) = v0linearintercept(12) - 0.225730_wp
1049	! !v0linearintercept(13) = v0linearintercept(13) - 0.238641_wp
1050	! !v0linearintercept(14) = v0linearintercept(14) - 3.005851_wp
1051	! !v0linearintercept(15) = v0linearintercept(15) - 0.000000_wp
1052
1053	! !v0linearintercept(1) = 0.11044026999999999_wp
1054	! !v0linearintercept(2) = 1.11152798999999990_wp
1055	! !v0linearintercept(3) = 2.11261570000000010_wp
1056	! !v0linearintercept(4) = -0.16579500000000000_wp
1057	! !v0linearintercept(5) = 6.31852534999999980_wp
1058	! !v0linearintercept(6) = 0.65734435000000002_wp
1059	! !v0linearintercept(7) = 4.04016423000000020_wp
1060	! !v0linearintercept(8) = 2.83395464000000000_wp
1061	! !v0linearintercept(9) = 6.21297795999999990_wp
1062	! !v0linearintercept(10) = 0.22088055000000001_wp
1063	! !v0linearintercept(11) = 0.00000000000000000_wp
1064	! !v0linearintercept(12) = -0.22572999999999999_wp
1065	! !v0linearintercept(13) = -0.23864099999999999_wp
1066	! !v0linearintercept(14) = -3.00585099999999980_wp
1067	! !v0linearintercept(15) = 3.14159265000000020_wp
1068
1069	! v0linearintercept( 1) = 0.2208805500_wp - (rpi* 68.0_wp/180.0_wp) ! M2 1
1070	! v0linearintercept( 2) = 3.1186126191_wp ! N2 2
1071	! v0linearintercept( 3) = 0.9305155436_wp ! 2N2 3
1072	! v0linearintercept( 4) = 0.0194858941_wp ! S2 4
1073	! v0linearintercept( 5) = -2.5213114949_wp ! K2 5
1074	! v0linearintercept( 6) = 6.5970532125_wp ! K1 6
1075	! v0linearintercept( 7) = 1.1115279900_wp ! O1 7
1076	! v0linearintercept( 8) = 0.1104402700_wp ! Q1 8
1077	! ! v0linearintercept( 9) = 4.2269096542_wp ! P1 9
1078	! !v0linearintercept( 9) = -2.0351042402_wp ! P1 9 compress3
1079	! !v0linearintercept( 9) = -2.0351042402_wp - 2.6179938779914944 ! P1 9 compress4
1080
1081	! v0linearintercept( 9) = rpi* 345.0_wp/180.0_wp - (rpi* 140.0_wp/180.0_wp) ! P1 9 compress4
1082
1083	! v0linearintercept(10) = 3.1415926500_wp ! M4 10
1084	! v0linearintercept(11) = 0.0000000000_wp ! Mf 11
1085	! v0linearintercept(12) = 0.0000000000_wp ! Mm 12
1086	! v0linearintercept(13) = 0.0000000000_wp ! Msqm 13
1087	! v0linearintercept(14) = 0.0000000000_wp ! Mtm 14
1088	! v0linearintercept(15) = -0.0230244122_wp ! S1 15
1089	! v0linearintercept(16) = 4.2565208698_wp ! MU2 16
1090	! v0linearintercept(17) = 6.5001767059_wp ! NU2 17
1091	! v0linearintercept(18) = 0.0000000000_wp - (rpi* 113.0_wp/180.0_wp) ! L2 18
1092	! v0linearintercept(19) = 0.0092971808_wp ! T2 19 + rpi/2.
1093
1094	! !v0linearintercept(1) = v0linearintercept(1) - 0.034975_wp ! M2
1095	! !v0linearintercept(2) = v0linearintercept(2) - 0.030244_wp ! N2
1096	! !v0linearintercept(3) = v0linearintercept(3) - 0.036046_wp ! 2N2
1097	! !v0linearintercept(4) = v0linearintercept(4) + 0.002092_wp ! S2
1098	! !v0linearintercept(5) = v0linearintercept(5) - 0.273826_wp ! K2
1099	! !v0linearintercept(6) = v0linearintercept(6) - 0.144677_wp ! K1
1100	! !v0linearintercept(7) = v0linearintercept(7) + 0.031938_wp ! O1
1101	! !v0linearintercept(8) = v0linearintercept(8) - 0.812030_wp ! Q1
1102	! !v0linearintercept(9) = v0linearintercept(9) + 2.109118_wp ! P1
1103	! !v0linearintercept(10) = v0linearintercept(10) + 0.070021_wp ! M4
1104	! !v0linearintercept(11) = v0linearintercept(11) - 0.000000_wp ! Mf
1105	! !v0linearintercept(12) = v0linearintercept(12) - 0.000000_wp ! Mm
1106	! !v0linearintercept(13) = v0linearintercept(13) - 0.000000_wp ! Msqm
1107	! !v0linearintercept(14) = v0linearintercept(14) - 0.000000_wp ! Mtm
1108	! !v0linearintercept(15) = v0linearintercept(15) - 0.035676_wp ! S1
1109	! !v0linearintercept(16) = v0linearintercept(16) + 0.007598_wp ! MU2
1110	! !v0linearintercept(17) = v0linearintercept(17) - 0.043060_wp ! NU2
1111	! !v0linearintercept(18) = v0linearintercept(18) + 0.023561_wp ! L2
1112	! !v0linearintercept(19) = v0linearintercept(19) + 0.025624_wp ! T2
1113
1114	! v0linearintercept(1) = v0linearintercept(1) - (rpi*2.003909_wp/180.0_wp) ! M2
1115	! v0linearintercept(2) = v0linearintercept(2) - (rpi*1.732874_wp/180.0_wp) ! N2
1116	! v0linearintercept(3) = v0linearintercept(3) - (rpi*2.065265_wp/180.0_wp) ! 2N2
1117	! v0linearintercept(4) = v0linearintercept(4) + (rpi*0.119842_wp/180.0_wp) ! S2
1118	! v0linearintercept(5) = v0linearintercept(5) - (rpi*15.689068_wp/180.0_wp) ! K2
1119	! v0linearintercept(6) = v0linearintercept(6) - (rpi*8.289390_wp/180.0_wp) ! K1
1120	! v0linearintercept(7) = v0linearintercept(7) + (rpi*1.829931_wp/180.0_wp) ! O1
1121	! v0linearintercept(8) = v0linearintercept(8) - (rpi*46.525902_wp/180.0_wp) ! Q1
1122	! v0linearintercept(9) = v0linearintercept(9) + (rpi*120.843575_wp/180.0_wp) ! P1
1123	! v0linearintercept(10) = v0linearintercept(10) + (rpi*4.011896_wp/180.0_wp) ! M4
1124	! v0linearintercept(11) = v0linearintercept(11) - (rpi*0.000000_wp/180.0_wp) ! Mf
1125	! v0linearintercept(12) = v0linearintercept(12) - (rpi*0.000000_wp/180.0_wp) ! Mm
1126	! v0linearintercept(13) = v0linearintercept(13) - (rpi*0.000000_wp/180.0_wp) ! Msqm
1127	! v0linearintercept(14) = v0linearintercept(14) - (rpi*0.000000_wp/180.0_wp) ! Mtm
1128	! v0linearintercept(15) = v0linearintercept(15) - (rpi*2.044069_wp/180.0_wp) ! S1
1129	! v0linearintercept(16) = v0linearintercept(16) + (rpi*0.435315_wp/180.0_wp) ! MU2
1130	! v0linearintercept(17) = v0linearintercept(17) - (rpi*2.467160_wp/180.0_wp) ! NU2
1131	! v0linearintercept(18) = v0linearintercept(18) + (rpi*1.349939_wp/180.0_wp) ! L2
1132	! v0linearintercept(19) = v0linearintercept(19) + (rpi*1.468170_wp/180.0_wp) ! T2
1133
1134
1135	! ! wave data.
1136
1137	! !Wave( 1) = tide( 'M2' , 0.242297 , 2 , 2 , -2 , 2 , 0 , 0 , 0 , 2 , -2 , 0 , 0 , 0 , 78 )
1138	! !Wave( 2) = tide( 'N2' , 0.046313 , 2 , 2 , -3 , 2 , 1 , 0 , 0 , 2 , -2 , 0 , 0 , 0 , 78 )
1139	! !Wave( 3) = tide( '2N2' , 0.006184 , 2 , 2 , -4 , 2 , 2 , 0 , 0 , 2 , -2 , 0 , 0 , 0 , 78 )
1140	! !Wave( 4) = tide( 'S2' , 0.113572 , 2 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 )
1141	! !Wave( 5) = tide( 'K2' , 0.030875 , 2 , 2 , 0 , 2 , 0 , 0 , 0 , 0 , 0 , 0 , -2 , 0 , 235 )
1142	! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
1143	! !Wave( 6) = tide( 'K1' , 0.142408 , 1 , 1 , 0 , 1 , 0 , 0 , -90 , 0 , 0 , -1 , 0 , 0 , 227 )
1144	! !Wave( 7) = tide( 'O1' , 0.101266 , 1 , 1 , -2 , 1 , 0 , 0 , +90 , 2 , -1 , 0 , 0 , 0 , 75 )
1145	! !Wave( 8) = tide( 'Q1' , 0.019387 , 1 , 1 , -3 , 1 , 1 , 0 , +90 , 2 , -1 , 0 , 0 , 0 , 75 )
1146	! !Wave( 9) = tide( 'P1' , 0.047129 , 1 , 1 , 0 , -1 , 0 , 0 , +90 , 0 , 0 , 0 , 0 , 0 , 0 )
1147	! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
1148	! !Wave(10) = tide( 'M4' , 0.000000 , 4 , 4 , -4 , 4 , 0 , 0 , 0 , 4 , -4 , 0 , 0 , 0 , 1 )
1149	! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
1150	! !Wave(11) = tide( 'Mf' , 0.042017 , 0 , 0 , 2 , 0 , 0 , 0 , 0 , -2 , 0 , 0 , 0 , 0 , 74 )
1151	! !Wave(12) = tide( 'Mm' , 0.022191 , 0 , 0 , 1 , 0 , -1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 73 )
1152	! !Wave(13) = tide( 'Msqm' , 0.000667 , 0 , 0 , 4 , -2 , 0 , 0 , 0 , -2 , 0 , 0 , 0 , 0 , 74 )
1153	! !Wave(14) = tide( 'Mtm' , 0.008049 , 0 , 0 , 3 , 0 , -1 , 0 , 0 , -2 , 0 , 0 , 0 , 0 , 74 )
1154	! !! ! ! ! ! ! ! ! ! ! ! ! ! ! ! ! !
1155	! !Wave(15) = tide( 'S1' , 0.000000 , 1 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 , 0 )
1156	! !Wave(16) = tide( 'MU2' , 0.005841 , 2 , 2 , -4 , 4 , 0 , 0 , 0 , 2 , -2 , 0 , 0 , 0 , 78 )
1157	! !Wave(17) = tide( 'NU2' , 0.009094 , 2 , 2 , -3 , 4 , -1 , 0 , 0 , 2 , -2 , 0 , 0 , 0 , 78 )
1158	! !Wave(18) = tide( 'L2' , 0.006694 , 2 , 2 , -1 , 2 , -1 , 0 , +180 , 2 , -2 , 0 , 0 , 0 , 215 )
1159	! !Wave(19) = tide( 'T2' , 0.006614 , 2 , 2 , 0 , -1 , 0 , 1 , 0 , 0 , 0 , 0 , 0 , 0 , 0 )
1160
1161	! !name list
1162	! ! clname(1)='Q1'
1163	! ! clname(2)='O1'
1164	! ! clname(3)='P1'
1165	! ! clname(4)='S1'
1166	! ! clname(5)='K1'
1167	! ! clname(6)='2N2'
1168	! ! clname(7)='MU2'
1169	! ! clname(8)='N2'
1170	! ! clname(9)='NU2'
1171	! ! clname(10)='M2'
1172	! ! clname(11)='L2'
1173	! ! clname(12)='T2'
1174	! ! clname(13)='S2'
1175	! ! clname(14)='K2'
1176	! ! clname(15)='M4'
1177
1178	! ! ktide 8,7,9,15
1179
1180	! !ktide =
1181	! !8
1182	! !7
1183	! !9
1184	! !15
1185	! !6
1186	! !3
1187	! !16
1188	! !2
1189	! !17
1190	! !1
1191	! !18
1192	! !19
1193	! !4
1194	! !5
1195	! !10
1196
1197
1198
1199
1200
1201
1202
1203
1204
1205
1206
1207
1208	! !NEMO4
1209
1210	!! clname(1)='Q1'
1211	!! clname(2)='O1'
1212	!! clname(3)='P1'
1213	!! clname(4)='S1'
1214	!! clname(5)='K1'
1215	!! clname(6)='2N2'
1216	!! clname(7)='MU2'
1217	!! clname(8)='N2'
1218	!! clname(9)='NU2'
1219	!! clname(10)='M2'
1220	!! clname(11)='L2'
1221	!! clname(12)='T2'
1222	!! clname(13)='S2'
1223	!! clname(14)='K2'
1224	!! clname(15)='M4'
1225	!! ktide = [10,9,11,12,8,23,21,15,22,14,18,19,16,17,28]
1226
1227
1228	! v0linearintercept( 1) = 0.1104402700_wp ! Q1 8
1229	! v0linearintercept( 2) = 1.1115279900_wp ! O1 7
1230	! v0linearintercept( 3) = rpi* 345.0_wp/180.0_wp - (rpi* 140.0_wp/180.0_wp) ! P1 9 compress4
1231	! v0linearintercept( 4) = -0.0230244122_wp ! S1 15
1232	! v0linearintercept( 5) = 6.5970532125_wp ! K1 6
1233	! v0linearintercept( 6) = 0.9305155436_wp ! 2N2 3
1234	! v0linearintercept( 7) = 4.2565208698_wp ! MU2 16
1235	! v0linearintercept( 8) = 3.1186126191_wp ! N2 2
1236	! v0linearintercept( 9) = 6.5001767059_wp ! NU2 17
1237	! v0linearintercept(10) = 0.2208805500_wp - (rpi* 68.0_wp/180.0_wp) ! M2 1
1238	! v0linearintercept(11) = 0.0000000000_wp - (rpi* 113.0_wp/180.0_wp) ! L2 18
1239	! v0linearintercept(12) = 0.0092971808_wp ! T2 19 + rpi/2.
1240	! v0linearintercept(13) = 0.0194858941_wp ! S2 4
1241	! v0linearintercept(14) = -2.5213114949_wp ! K2 5
1242	! v0linearintercept(15) = 3.1415926500_wp ! M4 10
1243
1244
1245
1246	! v0linearintercept( 1) = v0linearintercept( 1) - (rpi*46.525902_wp/180.0_wp) ! Q1
1247	! v0linearintercept( 2) = v0linearintercept( 2) + (rpi*1.829931_wp/180.0_wp) ! O1
1248	! v0linearintercept( 3) = v0linearintercept( 3) + (rpi*120.843575_wp/180.0_wp) ! P1
1249	! v0linearintercept( 4) = v0linearintercept( 4) - (rpi*2.044069_wp/180.0_wp) ! S1
1250	! v0linearintercept( 5) = v0linearintercept( 5) - (rpi*8.289390_wp/180.0_wp) ! K1
1251	! v0linearintercept( 6) = v0linearintercept( 6) - (rpi*2.065265_wp/180.0_wp) ! 2N2
1252	! v0linearintercept( 7) = v0linearintercept( 7) + (rpi*0.435315_wp/180.0_wp) ! MU2
1253	! v0linearintercept( 8) = v0linearintercept( 8) - (rpi*1.732874_wp/180.0_wp) ! N2
1254	! v0linearintercept( 9) = v0linearintercept( 9) - (rpi*2.467160_wp/180.0_wp) ! NU2
1255	! v0linearintercept(10) = v0linearintercept(10) - (rpi*2.003909_wp/180.0_wp) ! M2
1256	! v0linearintercept(11) = v0linearintercept(11) + (rpi*1.349939_wp/180.0_wp) ! L2
1257	! v0linearintercept(12) = v0linearintercept(12) + (rpi*1.468170_wp/180.0_wp) ! T2
1258	! v0linearintercept(13) = v0linearintercept(13) + (rpi*0.119842_wp/180.0_wp) ! S2
1259	! v0linearintercept(14) = v0linearintercept(14) - (rpi*15.689068_wp/180.0_wp) ! K2
1260	! v0linearintercept(14) = v0linearintercept(15) + (rpi*4.011896_wp/180.0_wp) ! M4
1261
1262
1263
1264
1265
1266
1267	DO jh = 1, kc
1268	IF(ln_astro_verbose .AND. lwp) WRITE(numout,*) 'astro tide_vuf 2:',jh,days_since_origin,v0linearslope(jh),v0linearintercept(ktide(jh))!,cycle_reset(jh)! ,offset(jh)
1269	ENDDO
1270	ENDIF
1271
1272
1273
1274	!!----------------------------------------------------------------------
1275	!!----------------------------------------------------------------------
1276	!! JT
1277	!!! phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
1278	!!----------------------------------------------------------------------
1279	!!----------------------------------------------------------------------
1280
1281	DO jh = 1, kc
1282	! Phase of the tidal potential relative to the Greenwhich
1283	! meridian (e.g. the position of the fictuous celestial body). Units are radian:
1284	! Linear with time
1285
1286	v0linearintercept(jh) = sh_T_o * Wave( ktide(jh) )%nT &
1287	& + sh_s_o * Wave( ktide(jh) )%ns &
1288	& + sh_h_o * Wave( ktide(jh) )%nh &
1289	& + sh_p_o * Wave( ktide(jh) )%np &
1290	& + sh_p1_o* Wave( ktide(jh) )%np1 &
1291	& + Wave( ktide(jh) )%shift * rad
1292
1293
1294	IF ( ln_tide_compress ) THEN
1295
1296	!JT pvt(jh) = mod( ( (v0linearslope(jh)days_since_origin) + v0linearintercept( ktide(jh) ) ), 2rpi)-(2*rpi)
1297	pvt(jh) = mod( ( (v0linearslope(jh)days_since_origin) + v0linearintercept( jh ) ), 2rpi)-(2*rpi)
1298	ELSE
1299	pvt(jh) = sh_T * Wave( ktide(jh) )%nT &
1300	& + sh_s * Wave( ktide(jh) )%ns &
1301	& + sh_h * Wave( ktide(jh) )%nh &
1302	& + sh_p * Wave( ktide(jh) )%np &
1303	& + sh_p1* Wave( ktide(jh) )%np1 &
1304	& + Wave( ktide(jh) )%shift * rad
1305	ENDIF
1306
1307
1308
1309	!
1310	! Phase correction u due to nodal motion. Units are radian:
1311	! Cyclical with time. Much smaller terms than pvt.
1312	put(jh) = sh_xi * Wave( ktide(jh) )%nksi &
1313	& + sh_nu * Wave( ktide(jh) )%nnu0 &
1314	& + sh_nuprim * Wave( ktide(jh) )%nnu1 &
1315	& + sh_nusec * Wave( ktide(jh) )%nnu2 &
1316	& + sh_R * Wave( ktide(jh) )%R
1317
1318
1319
1320
1321
1322	!JT
1323	pvt(jh) = mod( pvt(jh), 2*rpi)
1324	put(jh) = mod( put(jh), 2*rpi)
1325	!JT
1326
1327
1328	! Nodal correction factor:
1329	pcor(jh) = nodal_factort( Wave( ktide(jh) )%nformula )
1330
1331
1332	END DO
1333
1334
1335	IF(ln_astro_verbose .AND. lwp) THEN
1336	DO jh = 1, kc
1337	WRITE(numout,) 'astro tide_vuf 3,',jh,Wave(ktide(jh))%cname_tide, pomega(jh),2rpi/(3600.0_wppomega(jh)),pvt(jh), put(jh), pcor(jh), mod(pvt(jh),2rpi), mod(put(jh),2rpi), 2rpi
1338	END DO
1339	ENDIF
1340
1341
1342
1343
1344	DO jh = 1, kc
1345
1346	tmp_name=TRIM(Wave(ktide(jh))%cname_tide)//'_utide'
1347	IF( iom_use(TRIM(tmp_name)) ) THEN
1348	! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anau(jh) )
1349	CALL iom_put( TRIM(tmp_name), put(jh) )
1350	!ELSE
1351	! IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
1352	ENDIF
1353
1354	tmp_name=TRIM(Wave(ktide(jh))%cname_tide)//'_v0tide'
1355	IF( iom_use(TRIM(tmp_name)) ) THEN
1356	! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anav(jh) )
1357	CALL iom_put( TRIM(tmp_name), pvt(jh) )
1358	!ELSE
1359	! IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
1360	ENDIF
1361
1362	tmp_name=TRIM(Wave(ktide(jh))%cname_tide)//'_v0tide_origin'
1363	IF( iom_use(TRIM(tmp_name)) ) THEN
1364	! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anav(jh) )
1365	CALL iom_put( TRIM(tmp_name), v0linearintercept(jh) )
1366	!ELSE
1367	! IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
1368	ENDIF
1369
1370
1371
1372
1373	tmp_name=TRIM(Wave(ktide(jh))%cname_tide)//'_ftide'
1374	IF( iom_use(TRIM(tmp_name)) ) THEN
1375	! IF(lwp) WRITE(numout,*) "harm_ana_out: iom_put: ",TRIM(tmp_name),'; shape = ', SHAPE(anaf(jh) )
1376	CALL iom_put( TRIM(tmp_name), pcor(jh) )
1377	!ELSE
1378	! IF(lwp) WRITE(numout,*) "harm_ana_out: not requested: ",TRIM(tmp_name)
1379	ENDIF
1380
1381	END DO
1382
1383
1384
1385
1386
1387
1388
1389
1390
1391
1392
1393
1394	!
1395	END SUBROUTINE tide_vuf
1396
1397
1398	RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf )
1399	!!----------------------------------------------------------------------
1400	!!----------------------------------------------------------------------
1401	INTEGER, INTENT(in) :: kformula
1402	!
1403	REAL(wp) :: zf
1404	REAL(wp) :: zs, zf1, zf2
1405	!!----------------------------------------------------------------------
1406	!
1407	SELECT CASE( kformula )
1408	!
1409	CASE( 0 ) !== formule 0, solar waves
1410	zf = 1.0
1411	!
1412	CASE( 1 ) !== formule 1, compound waves (78 x 78)
1413	zf=nodal_factort(78)
1414	zf = zf * zf
1415	!
1416	CASE ( 2 ) !== formule 2, compound waves (78 x 0) === (78)
1417	zf1= nodal_factort(78)
1418	zf = nodal_factort( 0)
1419	zf = zf1 * zf
1420	!
1421	CASE ( 4 ) !== formule 4, compound waves (78 x 235)
1422	zf1 = nodal_factort( 78)
1423	zf = nodal_factort(235)
1424	zf = zf1 * zf
1425	!
1426	CASE ( 5 ) !== formule 5, compound waves (78 *78 x 235)
1427	zf1 = nodal_factort( 78)
1428	zf = nodal_factort(235)
1429	zf = zf * zf1 * zf1
1430	!
1431	CASE ( 6 ) !== formule 6, compound waves (78 *78 x 0)
1432	zf1 = nodal_factort(78)
1433	zf = nodal_factort( 0)
1434	zf = zf * zf1 * zf1
1435	!
1436	CASE( 7 ) !== formule 7, compound waves (75 x 75)
1437	zf = nodal_factort(75)
1438	zf = zf * zf
1439	!
1440	CASE( 8 ) !== formule 8, compound waves (78 x 0 x 235)
1441	zf = nodal_factort( 78)
1442	zf1 = nodal_factort( 0)
1443	zf2 = nodal_factort(235)
1444	zf = zf * zf1 * zf2
1445	!
1446	CASE( 9 ) !== formule 9, compound waves (78 x 0 x 227)
1447	zf = nodal_factort( 78)
1448	zf1 = nodal_factort( 0)
1449	zf2 = nodal_factort(227)
1450	zf = zf * zf1 * zf2
1451	!
1452	CASE( 10 ) !== formule 10, compound waves (78 x 227)
1453	zf = nodal_factort( 78)
1454	zf1 = nodal_factort(227)
1455	zf = zf * zf1
1456	!
1457	CASE( 11 ) !== formule 11, compound waves (75 x 0)
1458	!!gm bug???? zf 2 fois !
1459	zf = nodal_factort(75)
1460	zf1 = nodal_factort( 0)
1461	zf = zf * zf1
1462	!
1463	CASE( 12 ) !== formule 12, compound waves (78 x 78 x 78 x 0)
1464	zf1 = nodal_factort(78)
1465	zf = nodal_factort( 0)
1466	zf = zf * zf1 * zf1 * zf1
1467	!
1468	CASE( 13 ) !== formule 13, compound waves (78 x 75)
1469	zf1 = nodal_factort(78)
1470	zf = nodal_factort(75)
1471	zf = zf * zf1
1472	!
1473	CASE( 14 ) !== formule 14, compound waves (235 x 0) === (235)
1474	zf = nodal_factort(235)
1475	zf1 = nodal_factort( 0)
1476	zf = zf * zf1
1477	!
1478	CASE( 15 ) !== formule 15, compound waves (235 x 75)
1479	zf = nodal_factort(235)
1480	zf1 = nodal_factort( 75)
1481	zf = zf * zf1
1482	!
1483	CASE( 16 ) !== formule 16, compound waves (78 x 0 x 0) === (78)
1484	zf = nodal_factort(78)
1485	zf1 = nodal_factort( 0)
1486	zf = zf * zf1 * zf1
1487	!
1488	CASE( 17 ) !== formule 17, compound waves (227 x 0)
1489	zf1 = nodal_factort(227)
1490	zf = nodal_factort( 0)
1491	zf = zf * zf1
1492	!
1493	CASE( 18 ) !== formule 18, compound waves (78 x 78 x 78 )
1494	zf1 = nodal_factort(78)
1495	zf = zf1 * zf1 * zf1
1496	!
1497	CASE( 19 ) !== formule 19, compound waves (78 x 0 x 0 x 0) === (78)
1498	!!gm bug2 ==>>> here identical to formule 16, a third multiplication by zf1 is missing
1499	zf = nodal_factort(78)
1500	zf1 = nodal_factort( 0)
1501	zf = zf * zf1 * zf1
1502	!
1503
1504	!--- davbyr 11/2017
1505	CASE( 20 ) !== formule 20, compound waves ( 78 x 78 x 78 x 78 )
1506	zf1 = nodal_factort(78)
1507	zf = zf1 * zf1 * zf1 * zf1
1508	!--- END davbyr
1509	CASE( 73 ) !== formule 73
1510	zs = sin(sh_I)
1511	zf = (2./3.-zs*zs)/0.5021
1512	!
1513	CASE( 74 ) !== formule 74
1514	zs = sin(sh_I)
1515	zf = zs * zs / 0.1578
1516	!
1517	CASE( 75 ) !== formule 75
1518	zs = cos(sh_I/2)
1519	zf = sin(sh_I) * zs * zs / 0.3800
1520	!
1521	CASE( 76 ) !== formule 76
1522	zf = sin(2*sh_I) / 0.7214
1523	!
1524	CASE( 77 ) !== formule 77
1525	zs = sin(sh_I/2)
1526	zf = sin(sh_I) * zs * zs / 0.0164
1527	!
1528	CASE( 78 ) !== formule 78
1529	zs = cos(sh_I/2)
1530	zf = zs * zs * zs * zs / 0.9154
1531	!
1532	CASE( 79 ) !== formule 79
1533	zs = sin(sh_I)
1534	zf = zs * zs / 0.1565
1535	!
1536	CASE( 144 ) !== formule 144
1537	zs = sin(sh_I/2)
1538	zf = ( 1-10zszs+15zszszszs ) * cos(sh_I/2) / 0.5873
1539	!
1540	CASE( 149 ) !== formule 149
1541	zs = cos(sh_I/2)
1542	zf = zszszszszs*zs / 0.8758
1543	!
1544	CASE( 215 ) !== formule 215
1545	zs = cos(sh_I/2)
1546	zf = zszszszs / 0.9154 sh_x1ra
1547	!
1548	CASE( 227 ) !== formule 227
1549	zs = sin(2*sh_I)
1550	zf = sqrt( 0.8965zszs+0.6001zscos (sh_nu)+0.1006 )
1551	!
1552	CASE ( 235 ) !== formule 235
1553	zs = sin(sh_I)
1554	zf = sqrt( 19.0444zszszszs + 2.7702zszscos(2sh_nu) + .0981 )
1555	!
1556	END SELECT
1557	!
1558	END FUNCTION nodal_factort
1559
1560
1561	FUNCTION dayjul( kyr, kmonth, kday )
1562	!!----------------------------------------------------------------------
1563	!! *** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE)
1564	!!----------------------------------------------------------------------
1565	INTEGER,INTENT(in) :: kyr, kmonth, kday
1566	!
1567	INTEGER,DIMENSION(12) :: idayt, idays
1568	INTEGER :: inc, ji
1569	REAL(wp) :: dayjul, zyq
1570	!
1571	DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
1572	!!----------------------------------------------------------------------
1573	!
1574	idays(1) = 0.
1575	idays(2) = 31.
1576	inc = 0.
1577	zyq = MOD( kyr-1900. , 4. )
1578	IF( zyq == 0.) inc = 1.
1579	DO ji = 3, 12
1580	idays(ji)=idayt(ji)+inc
1581	END DO
1582	dayjul = idays(kmonth) + kday
1583	!
1584	END FUNCTION dayjul
1585
1586
1587
1588
1589
1590
1591	SUBROUTINE ju2ymds_JT (julian,year,month,day,sec,one_year)
1592	!---------------------------------------------------------------------
1593	! IMPLICIT NONE
1594	!-
1595	! REAL,INTENT(IN) :: julian
1596	!-
1597	! INTEGER,INTENT(OUT) :: year,month,day
1598	! REAL,INTENT(OUT) :: sec
1599	!-
1600	! INTEGER :: julian_day
1601	! REAL :: julian_sec
1602	!---------------------------------------------------------------------
1603	! julian_day = INT(julian)
1604	! julian_sec = (julian-julian_day)*one_day
1605	!-
1606	! CALL ju2ymds_internal(julian_day,julian_sec,year,month,day,sec)
1607	!---------------------
1608	!END SUBROUTINE ju2ymds
1609	!-
1610	!===
1611	!-
1612	!SUBROUTINE ju2ymds_internal (julian_day,julian_sec,year,month,day,sec)
1613	!---------------------------------------------------------------------
1614	!- This subroutine computes from the julian day the year,
1615	!- month, day and seconds
1616	!-
1617	!- In 1968 in a letter to the editor of Communications of the ACM
1618	!- (CACM, volume 11, number 10, October 1968, p.657) Henry F. Fliegel
1619	!- and Thomas C. Van Flandern presented such an algorithm.
1620	!-
1621	!- See also : http://www.magnet.ch/serendipity/hermetic/cal_stud/jdn.htm
1622	!-
1623	!- In the case of the Gregorian calendar we have chosen to use
1624	!- the Lilian day numbers. This is the day counter which starts
1625	!- on the 15th October 1582. This is the day at which Pope
1626	!- Gregory XIII introduced the Gregorian calendar.
1627	!- Compared to the true Julian calendar, which starts some 7980
1628	!- years ago, the Lilian days are smaler and are dealt with easily
1629	!- on 32 bit machines. With the true Julian days you can only the
1630	!- fraction of the day in the real part to a precision of a 1/4 of
1631	!- a day with 32 bits.
1632	!---------------------------------------------------------------------
1633	IMPLICIT NONE
1634
1635
1636	!-
1637
1638	REAL,INTENT(IN) :: julian,one_year
1639	!-
1640	INTEGER,INTENT(OUT) :: year,month,day
1641	REAL,INTENT(OUT) :: sec
1642	!-
1643	INTEGER :: l,n,i,jd,j,d,m,y,ml
1644	INTEGER :: add_day
1645	REAL :: eps_day
1646
1647	REAL,PARAMETER :: one_day = 86400.0
1648	!---------------------------------------------------------------------
1649
1650	INTEGER :: julian_day
1651	REAL :: julian_sec
1652	INTEGER :: mon_len(12)
1653	!---------------------------------------------------------------------
1654
1655	IF ( (one_year > 365.0).AND.(one_year < 366.0) ) THEN
1656	mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
1657	ELSE IF ( ABS(one_year-365.0) <= EPSILON(one_year) ) THEN
1658	mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
1659	ELSE IF ( ABS(one_year-366.0) <= EPSILON(one_year) ) THEN
1660	mon_len(:)=(/31,29,31,30,31,30,31,31,30,31,30,31/)
1661	ELSE IF ( ABS(one_year-360.0) <= EPSILON(one_year) ) THEN
1662	mon_len(:)=(/30,30,30,30,30,30,30,30,30,30,30,30/)
1663	ENDIF
1664
1665
1666
1667	julian_day = INT(julian)
1668	julian_sec = (julian-julian_day)*one_day
1669
1670	eps_day = SPACING(one_day)
1671	! lock_one_year = .TRUE.
1672	!-
1673	jd = julian_day
1674	sec = julian_sec
1675	IF (sec > (one_day-eps_day)) THEN
1676	add_day = INT(sec/one_day)
1677	sec = sec-add_day*one_day
1678	jd = jd+add_day
1679	ENDIF
1680	IF (sec < -eps_day) THEN
1681	sec = sec+one_day
1682	jd = jd-1
1683	ENDIF
1684	!-
1685	IF ( (one_year > 365.0).AND.(one_year < 366.0) ) THEN
1686	!-- Gregorian
1687	jd = jd+2299160
1688	!-
1689	l = jd+68569
1690	n = (4*l)/146097
1691	l = l-(146097*n+3)/4
1692	i = (4000*(l+1))/1461001
1693	l = l-(1461*i)/4+31
1694	j = (80*l)/2447
1695	d = l-(2447*j)/80
1696	l = j/11
1697	m = j+2-(12*l)
1698	y = 100*(n-49)+i+l
1699	ELSE IF ( (ABS(one_year-365.0) <= EPSILON(one_year)) &
1700	& .OR.(ABS(one_year-366.0) <= EPSILON(one_year)) ) THEN
1701	!-- No leap or All leap
1702	!if ( ABS(one_year-365.0) <= EPSILON(one_year) ) mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
1703	!if ( ABS(one_year-366.0) <= EPSILON(one_year) ) mon_len(:)=(/31,29,31,30,31,30,31,31,30,31,30,31/)
1704	y = jd/NINT(one_year)
1705	l = jd-y*NINT(one_year)
1706	m = 1
1707	ml = 0
1708	DO WHILE (ml+mon_len(m) <= l)
1709	ml = ml+mon_len(m)
1710	m = m+1
1711	ENDDO
1712	d = l-ml+1
1713	ELSE
1714	!-- others
1715	ml = NINT(one_year/12.)
1716	y = jd/NINT(one_year)
1717	l = jd-y*NINT(one_year)
1718	m = (l/ml)+1
1719	d = l-(m-1)*ml+1
1720	ENDIF
1721	!-
1722	day = d
1723	month = m
1724	year = y
1725	!------------------------------
1726
1727
1728	END SUBROUTINE ju2ymds_JT
1729
1730
1731
1732
1733
1734
1735
1736
1737	!-
1738	SUBROUTINE ymds2ju_JT (year,month,day,sec,julian,one_year)
1739	!---------------------------------------------------------------------
1740	! IMPLICIT NONE
1741	!-
1742	! INTEGER,INTENT(IN) :: year,month,day
1743	! REAL,INTENT(IN) :: sec
1744	!-
1745	! REAL,INTENT(OUT) :: julian
1746	!-
1747	! INTEGER :: julian_day
1748	! REAL :: julian_sec
1749	!---------------------------------------------------------------------
1750	! CALL ymds2ju_internal (year,month,day,sec,julian_day,julian_sec)
1751	!-
1752	!---------------------
1753	!END SUBROUTINE ymds2ju
1754	!-
1755	!===
1756	!-
1757	!SUBROUTINE ymds2ju_internal (year,month,day,sec,julian_day,julian_sec)
1758	!---------------------------------------------------------------------
1759	!- Converts year, month, day and seconds into a julian day
1760	!-
1761	!- In 1968 in a letter to the editor of Communications of the ACM
1762	!- (CACM, volume 11, number 10, October 1968, p.657) Henry F. Fliegel
1763	!- and Thomas C. Van Flandern presented such an algorithm.
1764	!-
1765	!- See also : http://www.magnet.ch/serendipity/hermetic/cal_stud/jdn.htm
1766	!-
1767	!- In the case of the Gregorian calendar we have chosen to use
1768	!- the Lilian day numbers. This is the day counter which starts
1769	!- on the 15th October 1582.
1770	!- This is the day at which Pope Gregory XIII introduced the
1771	!- Gregorian calendar.
1772	!- Compared to the true Julian calendar, which starts some
1773	!- 7980 years ago, the Lilian days are smaler and are dealt with
1774	!- easily on 32 bit machines. With the true Julian days you can only
1775	!- the fraction of the day in the real part to a precision of
1776	!- a 1/4 of a day with 32 bits.
1777	!---------------------------------------------------------------------
1778	IMPLICIT NONE
1779	!-
1780	INTEGER,INTENT(IN) :: year,month,day
1781	REAL,INTENT(IN) :: sec
1782	REAL,INTENT(IN) :: one_year
1783	!-
1784	! INTEGER,INTENT(OUT) :: julian_day
1785	! REAL,INTENT(OUT) :: julian_sec
1786	REAL,INTENT(OUT) :: julian
1787	INTEGER :: julian_day
1788	REAL :: julian_sec
1789	!-
1790	INTEGER :: jd,m,y,d,ml
1791	REAL,PARAMETER :: one_day = 86400.0
1792
1793
1794	INTEGER :: mon_len(12)
1795
1796	IF ( (one_year > 365.0).AND.(one_year < 366.0) ) THEN
1797	mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
1798	ELSE IF ( ABS(one_year-365.0) <= EPSILON(one_year) ) THEN
1799	mon_len(:)=(/31,28,31,30,31,30,31,31,30,31,30,31/)
1800	ELSE IF ( ABS(one_year-366.0) <= EPSILON(one_year) ) THEN
1801	mon_len(:)=(/31,29,31,30,31,30,31,31,30,31,30,31/)
1802	ELSE IF ( ABS(one_year-360.0) <= EPSILON(one_year) ) THEN
1803	mon_len(:)=(/30,30,30,30,30,30,30,30,30,30,30,30/)
1804	ENDIF
1805
1806
1807	!---------------------------------------------------------------------
1808	! lock_one_year = .TRUE.
1809	!-
1810
1811
1812	m = month
1813	y = year
1814	d = day
1815
1816
1817	!---------------------------------------------------------------------
1818
1819
1820	!-
1821	!- We deduce the calendar from the length of the year as it
1822	!- is faster than an INDEX on the calendar variable.
1823	!-
1824	IF ( (one_year > 365.0).AND.(one_year < 366.0) ) THEN
1825	!-- "Gregorian"
1826	jd = (1461*(y+4800+INT((m-14)/12)))/4 &
1827	& +(367(m-2-12(INT((m-14)/12))))/12 &
1828	& -(3*((y+4900+INT((m-14)/12))/100))/4 &
1829	& +d-32075
1830	jd = jd-2299160
1831	ELSE IF ( (ABS(one_year-365.0) <= EPSILON(one_year)) &
1832	& .OR.(ABS(one_year-366.0) <= EPSILON(one_year)) ) THEN
1833	!-- "No leap" or "All leap"
1834	ml = SUM(mon_len(1:m-1))
1835	jd = y*NINT(one_year)+ml+(d-1)
1836	ELSE
1837	!-- Calendar with regular month
1838	ml = NINT(one_year/12.)
1839	jd = yNINT(one_year)+(m-1)ml+(d-1)
1840	ENDIF
1841	!-
1842	julian_day = jd
1843	julian_sec = sec
1844
1845	julian = julian_day+julian_sec/one_day
1846
1847	!------------------------------
1848	END SUBROUTINE ymds2ju_JT
1849
1850
1851
1852
1853
1854
1855
1856
1857	!!======================================================================
1858	END MODULE tide_mod

Note: See TracBrowser for help on using the repository browser.

Download in other formats: