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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
tide_mod.F90 in NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/SBC – NEMO

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

Last change on this file since 15500 was 15500, checked in by hadjt, 8 months 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 
1MODULE 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   !!----------------------------------------------------------------------
66CONTAINS
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 )
89901   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 )
93902   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.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*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-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad
674      !----------------------------------------------------------------------
675      ! s mean lunar Longitude
676      !----------------------------------------------------------------------
677      sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad
678      !----------------------------------------------------------------------
679      ! p1 Longitude of solar perigee
680      !----------------------------------------------------------------------
681      sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad
682      !----------------------------------------------------------------------
683      ! p Longitude of lunar perigee
684      !----------------------------------------------------------------------
685      sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*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.0*t2*cos(2.0*P1)+36.0*t4 )
733
734      p = sin(2.0*P1)
735      q = 1.0/(6.0*t2)-cos(2.0*P1)
736      sh_R = atan(p/q)
737
738      p = sin(2.0*sh_I)*sin(sh_nu)
739      q = sin(2.0*sh_I)*cos(sh_nu)+0.3347
740      sh_nuprim = atan(p/q)
741
742      s2 = sin(sh_I)*sin(sh_I)
743      p  = s2*sin(2.0*sh_nu)
744      q  = s2*cos(2.0*sh_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.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*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-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad
801      !----------------------------------------------------------------------
802      ! s mean lunar Longitude
803      !----------------------------------------------------------------------
804      sh_s_o=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad
805      !----------------------------------------------------------------------
806      ! p1 Longitude of solar perigee
807      !----------------------------------------------------------------------
808      sh_p1_o=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad
809      !----------------------------------------------------------------------
810      ! p Longitude of lunar perigee
811      !----------------------------------------------------------------------
812      sh_p_o=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*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.0*t2*cos(2.0*P1)+36.0*t4 )
866!
867!      p = sin(2.0*P1)
868!      q = 1.0/(6.0*t2)-cos(2.0*P1)
869!      sh_R = atan(p/q)
870!
871!      p = sin(2.0*sh_I)*sin(sh_nu)
872!      q = sin(2.0*sh_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  = s2*sin(2.0*sh_nu_o)
877!      q  = s2*cos(2.0*sh_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),2*rpi/(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_origin*freq(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),(2*rpi/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)  ) ),  2*rpi)-(2*rpi)
1288             pvt(jh) = mod(  ( (v0linearslope(jh)*days_since_origin) + v0linearintercept(  jh  ) ),  2*rpi)-(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),2*rpi/(3600.0_wp*pomega(jh)),pvt(jh), put(jh), pcor(jh), mod(pvt(jh),2*rpi), mod(put(jh),2*rpi), 2*rpi
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-10*zs*zs+15*zs*zs*zs*zs ) * cos(sh_I/2) / 0.5873
1479         !
1480      CASE( 149 )                !==  formule 149
1481         zs = cos(sh_I/2)
1482         zf = zs*zs*zs*zs*zs*zs / 0.8758
1483         !
1484      CASE( 215 )                !==  formule 215
1485         zs = cos(sh_I/2)
1486         zf = zs*zs*zs*zs / 0.9154 * sh_x1ra
1487         !
1488      CASE( 227 )                !==  formule 227
1489         zs = sin(2*sh_I)
1490         zf = sqrt( 0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006 )
1491         !
1492      CASE ( 235 )               !==  formule 235
1493         zs = sin(sh_I)
1494         zf = sqrt( 19.0444*zs*zs*zs*zs + 2.7702*zs*zs*cos(2*sh_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 = y*NINT(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   !!======================================================================
1798END MODULE tide_mod
Note: See TracBrowser for help on using the repository browser.