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 @ 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 
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
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   !!----------------------------------------------------------------------
58CONTAINS
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 )
92901   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 )
96902   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.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*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-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad
597      !----------------------------------------------------------------------
598      ! s mean lunar Longitude
599      !----------------------------------------------------------------------
600      sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad
601      !----------------------------------------------------------------------
602      ! p1 Longitude of solar perigee
603      !----------------------------------------------------------------------
604      sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad
605      !----------------------------------------------------------------------
606      ! p Longitude of lunar perigee
607      !----------------------------------------------------------------------
608      sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*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.0*t2*cos(2.0*P1)+36.0*t4 )
656
657      p = sin(2.0*P1)
658      q = 1.0/(6.0*t2)-cos(2.0*P1)
659      sh_R = atan(p/q)
660
661      p = sin(2.0*sh_I)*sin(sh_nu)
662      q = sin(2.0*sh_I)*cos(sh_nu)+0.3347
663      sh_nuprim = atan(p/q)
664
665      s2 = sin(sh_I)*sin(sh_I)
666      p  = s2*sin(2.0*sh_nu)
667      q  = s2*cos(2.0*sh_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),2*rpi/(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_origin*freq(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),(2*rpi/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)  ) ),  2*rpi)-(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-10*zs*zs+15*zs*zs*zs*zs ) * cos(sh_I/2) / 0.5873
1230         !
1231      CASE( 149 )                !==  formule 149
1232         zs = cos(sh_I/2)
1233         zf = zs*zs*zs*zs*zs*zs / 0.8758
1234         !
1235      CASE( 215 )                !==  formule 215
1236         zs = cos(sh_I/2)
1237         zf = zs*zs*zs*zs / 0.9154 * sh_x1ra
1238         !
1239      CASE( 227 )                !==  formule 227
1240         zs = sin(2*sh_I)
1241         zf = sqrt( 0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006 )
1242         !
1243      CASE ( 235 )               !==  formule 235
1244         zs = sin(sh_I)
1245         zf = sqrt( 19.0444*zs*zs*zs*zs + 2.7702*zs*zs*cos(2*sh_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 = y*NINT(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   !!======================================================================
1549END MODULE tide_mod
Note: See TracBrowser for help on using the repository browser.