Context Navigation

tide_mod.F90 @ 15500

Last change on this file since 15500 was 15500, checked in by hadjt, 3 years ago

diaharm_fast.F90

Found an important bug. it wasn't set up for monthly cycles, so was only reseting the "seconds of the day" count once a month. Therefore, the v0tide was introducing a jump every day, so was fitting poorly.

It is still not working properly, as the amplitudes output are too high, and the phase maps are not correct.

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

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

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

Context Navigation

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

Download in other formats: