Context Navigation

tide_mod.F90 @ 15490

Last change on this file since 15490 was 15490, checked in by hadjt, 8 months ago

SBC/tide_mod.F90

Adding updated tide module to handle 360 day tides. This is working, but needs some more work.

Care need for compress option, as unpublished, work in progress. Worked needed to set phase offsets (in v0tide).

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

Download in other formats: