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/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE – NEMO

source: NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE/tide_mod.F90 @ 10772

Last change on this file since 10772 was 10772, checked in by smueller, 3 years ago

Removal of module tideini and inclusion of its contents into module tide_mod (ticket #2194)

  • Property svn:keywords set to Id
File size: 20.6 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   USE oce            ! ocean dynamics and tracers variables
9   USE dom_oce        ! ocean space and time domain
10   USE phycst         ! physical constant
11   USE daymod         ! calendar
12   !
13   USE in_out_manager ! I/O units
14   USE iom            ! xIOs server
15   USE ioipsl         ! NetCDF IPSL library
16   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
17
18   IMPLICIT NONE
19   PRIVATE
20
21   PUBLIC   tide_init
22   PUBLIC   tide_harmo       ! called by tideini and diaharm modules
23   PUBLIC   tide_init_Wave   ! called by tideini and diaharm modules
24
25   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 19   !: 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), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   omega_tide   !:
39   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   v0tide       !:
40   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   utide        !:
41   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:) ::   ftide        !:
42
43   LOGICAL , PUBLIC ::   ln_tide         !:
44   LOGICAL , PUBLIC ::   ln_tide_pot     !:
45   LOGICAL , PUBLIC ::   ln_read_load    !:
46   LOGICAL , PUBLIC ::   ln_scal_load    !:
47   LOGICAL , PUBLIC ::   ln_tide_ramp    !:
48   INTEGER , PUBLIC ::   nb_harmo        !:
49   INTEGER , PUBLIC ::   kt_tide         !:
50   REAL(wp), PUBLIC ::   rdttideramp     !:
51   REAL(wp), PUBLIC ::   rn_scal_load    !:
52   CHARACTER(lc), PUBLIC ::   cn_tide_load   !:
53
54   INTEGER , PUBLIC, ALLOCATABLE, DIMENSION(:) ::   ntide   !:
55
56   REAL(wp) ::   sh_T, sh_s, sh_h, sh_p, sh_p1             ! astronomic angles
57   REAL(wp) ::   sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R   !
58   REAL(wp) ::   sh_I, sh_x1ra, sh_N                       !
59
60   !!----------------------------------------------------------------------
61   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
62   !! $Id$
63   !! Software governed by the CeCILL license (see ./LICENSE)
64   !!----------------------------------------------------------------------
65CONTAINS
66
67   SUBROUTINE tide_init
68      !!----------------------------------------------------------------------
69      !!                 ***  ROUTINE tide_init  ***
70      !!----------------------------------------------------------------------     
71      INTEGER  :: ji, jk
72      CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: clname
73      INTEGER  ::   ios                 ! Local integer output status for namelist read
74      !
75      NAMELIST/nam_tide/ln_tide, ln_tide_pot, ln_scal_load, ln_read_load, cn_tide_load, &
76                  &     ln_tide_ramp, rn_scal_load, rdttideramp, clname
77      !!----------------------------------------------------------------------
78      !
79      ! Read Namelist nam_tide
80      REWIND( numnam_ref )              ! Namelist nam_tide in reference namelist : Tides
81      READ  ( numnam_ref, nam_tide, IOSTAT = ios, ERR = 901)
82901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_tide in reference namelist', lwp )
83      !
84      REWIND( numnam_cfg )              ! Namelist nam_tide in configuration namelist : Tides
85      READ  ( numnam_cfg, nam_tide, IOSTAT = ios, ERR = 902 )
86902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_tide in configuration namelist', lwp )
87      IF(lwm) WRITE ( numond, nam_tide )
88      !
89      IF( ln_tide ) THEN
90         IF (lwp) THEN
91            WRITE(numout,*)
92            WRITE(numout,*) 'tide_init : Initialization of the tidal components'
93            WRITE(numout,*) '~~~~~~~~~ '
94            WRITE(numout,*) '   Namelist nam_tide'
95            WRITE(numout,*) '      Use tidal components                       ln_tide      = ', ln_tide
96            WRITE(numout,*) '         Apply astronomical potential            ln_tide_pot  = ', ln_tide_pot
97            WRITE(numout,*) '         Use scalar approx. for load potential   ln_scal_load = ', ln_scal_load
98            WRITE(numout,*) '         Read load potential from file           ln_read_load = ', ln_read_load
99            WRITE(numout,*) '         Apply ramp on tides at startup          ln_tide_ramp = ', ln_tide_ramp
100            WRITE(numout,*) '         Fraction of SSH used in scal. approx.   rn_scal_load = ', rn_scal_load
101            WRITE(numout,*) '         Duration (days) of ramp                 rdttideramp  = ', rdttideramp
102         ENDIF
103      ELSE
104         rn_scal_load = 0._wp 
105
106         IF(lwp) WRITE(numout,*)
107         IF(lwp) WRITE(numout,*) 'tide_init : tidal components not used (ln_tide = F)'
108         IF(lwp) WRITE(numout,*) '~~~~~~~~~ '
109         RETURN
110      ENDIF
111      !
112      CALL tide_init_Wave
113      !
114      nb_harmo=0
115      DO jk = 1, jpmax_harmo
116         DO ji = 1,jpmax_harmo
117            IF( TRIM(clname(jk)) == Wave(ji)%cname_tide )   nb_harmo = nb_harmo + 1
118         END DO
119      END DO
120      !       
121      ! Ensure that tidal components have been set in namelist_cfg
122      IF( nb_harmo == 0 )   CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' )
123      !
124      IF( ln_read_load.AND.(.NOT.ln_tide_pot) ) &
125          &   CALL ctl_stop('ln_read_load requires ln_tide_pot')
126      IF( ln_scal_load.AND.(.NOT.ln_tide_pot) ) &
127          &   CALL ctl_stop('ln_scal_load requires ln_tide_pot')
128      IF( ln_scal_load.AND.ln_read_load ) &
129          &   CALL ctl_stop('Choose between ln_scal_load and ln_read_load')
130      IF( ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rdttideramp) )   &
131         &   CALL ctl_stop('rdttideramp must be lower than run duration')
132      IF( ln_tide_ramp.AND.(rdttideramp<0.) ) &
133         &   CALL ctl_stop('rdttideramp must be positive')
134      !
135      ALLOCATE( ntide(nb_harmo) )
136      DO jk = 1, nb_harmo
137         DO ji = 1, jpmax_harmo
138            IF( TRIM(clname(jk)) == Wave(ji)%cname_tide ) THEN
139               ntide(jk) = ji
140               EXIT
141            ENDIF
142         END DO
143      END DO
144      !
145      ALLOCATE( omega_tide(nb_harmo), v0tide    (nb_harmo),   &
146         &      utide     (nb_harmo), ftide     (nb_harmo)  )
147      kt_tide = nit000
148      !
149      IF (.NOT.ln_scal_load ) rn_scal_load = 0._wp
150      !
151   END SUBROUTINE tide_init
152
153
154   SUBROUTINE tide_init_Wave
155#     include "tide.h90"
156   END SUBROUTINE tide_init_Wave
157
158
159   SUBROUTINE tide_harmo( pomega, pvt, put , pcor, ktide ,kc)
160      !!----------------------------------------------------------------------
161      !!----------------------------------------------------------------------
162      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide            ! Indice of tidal constituents
163      INTEGER                , INTENT(in ) ::   kc               ! Total number of tidal constituents
164      REAL(wp), DIMENSION(kc), INTENT(out) ::   pomega           ! pulsation in radians/s
165      REAL(wp), DIMENSION(kc), INTENT(out) ::   pvt, put, pcor   !
166      !!----------------------------------------------------------------------
167      !
168      CALL astronomic_angle
169      CALL tide_pulse( pomega, ktide ,kc )
170      CALL tide_vuf  ( pvt, put, pcor, ktide ,kc )
171      !
172   END SUBROUTINE tide_harmo
173
174
175   SUBROUTINE astronomic_angle
176      !!----------------------------------------------------------------------
177      !!  tj is time elapsed since 1st January 1900, 0 hour, counted in julian
178      !!  century (e.g. time in days divide by 36525)
179      !!----------------------------------------------------------------------
180      REAL(wp) ::   cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2
181      REAL(wp) ::   zqy , zsy, zday, zdj, zhfrac
182      !!----------------------------------------------------------------------
183      !
184      zqy = AINT( (nyear-1901.)/4. )
185      zsy = nyear - 1900.
186      !
187      zdj  = dayjul( nyear, nmonth, nday )
188      zday = zdj + zqy - 1.
189      !
190      zhfrac = nsec_day / 3600.
191      !
192      !----------------------------------------------------------------------
193      !  Sh_n Longitude of ascending lunar node
194      !----------------------------------------------------------------------
195      sh_N=(259.1560564-19.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*rad
196      !----------------------------------------------------------------------
197      ! T mean solar angle (Greenwhich time)
198      !----------------------------------------------------------------------
199      sh_T=(180.+zhfrac*(360./24.))*rad
200      !----------------------------------------------------------------------
201      ! h mean solar Longitude
202      !----------------------------------------------------------------------
203      sh_h=(280.1895014-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad
204      !----------------------------------------------------------------------
205      ! s mean lunar Longitude
206      !----------------------------------------------------------------------
207      sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad
208      !----------------------------------------------------------------------
209      ! p1 Longitude of solar perigee
210      !----------------------------------------------------------------------
211      sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad
212      !----------------------------------------------------------------------
213      ! p Longitude of lunar perigee
214      !----------------------------------------------------------------------
215      sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad
216
217      sh_N = MOD( sh_N ,2*rpi )
218      sh_s = MOD( sh_s ,2*rpi )
219      sh_h = MOD( sh_h, 2*rpi )
220      sh_p = MOD( sh_p, 2*rpi )
221      sh_p1= MOD( sh_p1,2*rpi )
222
223      cosI = 0.913694997 -0.035692561 *cos(sh_N)
224
225      sh_I = ACOS( cosI )
226
227      sin2I   = sin(sh_I)
228      sh_tgn2 = tan(sh_N/2.0)
229
230      at1=atan(1.01883*sh_tgn2)
231      at2=atan(0.64412*sh_tgn2)
232
233      sh_xi=-at1-at2+sh_N
234
235      IF( sh_N > rpi )   sh_xi=sh_xi-2.0*rpi
236
237      sh_nu = at1 - at2
238
239      !----------------------------------------------------------------------
240      ! For constituents l2 k1 k2
241      !----------------------------------------------------------------------
242
243      tgI2 = tan(sh_I/2.0)
244      P1   = sh_p-sh_xi
245
246      t2 = tgI2*tgI2
247      t4 = t2*t2
248      sh_x1ra = sqrt( 1.0-12.0*t2*cos(2.0*P1)+36.0*t4 )
249
250      p = sin(2.0*P1)
251      q = 1.0/(6.0*t2)-cos(2.0*P1)
252      sh_R = atan(p/q)
253
254      p = sin(2.0*sh_I)*sin(sh_nu)
255      q = sin(2.0*sh_I)*cos(sh_nu)+0.3347
256      sh_nuprim = atan(p/q)
257
258      s2 = sin(sh_I)*sin(sh_I)
259      p  = s2*sin(2.0*sh_nu)
260      q  = s2*cos(2.0*sh_nu)+0.0727
261      sh_nusec = 0.5*atan(p/q)
262      !
263   END SUBROUTINE astronomic_angle
264
265
266   SUBROUTINE tide_pulse( pomega, ktide ,kc )
267      !!----------------------------------------------------------------------
268      !!                     ***  ROUTINE tide_pulse  ***
269      !!                     
270      !! ** Purpose : Compute tidal frequencies
271      !!----------------------------------------------------------------------
272      INTEGER                , INTENT(in ) ::   kc       ! Total number of tidal constituents
273      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide    ! Indice of tidal constituents
274      REAL(wp), DIMENSION(kc), INTENT(out) ::   pomega   ! pulsation in radians/s
275      !
276      INTEGER  ::   jh
277      REAL(wp) ::   zscale
278      REAL(wp) ::   zomega_T =  13149000.0_wp
279      REAL(wp) ::   zomega_s =    481267.892_wp
280      REAL(wp) ::   zomega_h =     36000.76892_wp
281      REAL(wp) ::   zomega_p =      4069.0322056_wp
282      REAL(wp) ::   zomega_n =      1934.1423972_wp
283      REAL(wp) ::   zomega_p1=         1.719175_wp
284      !!----------------------------------------------------------------------
285      !
286      zscale =  rad / ( 36525._wp * 86400._wp ) 
287      !
288      DO jh = 1, kc
289         pomega(jh) = (  zomega_T * Wave( ktide(jh) )%nT   &
290            &          + zomega_s * Wave( ktide(jh) )%ns   &
291            &          + zomega_h * Wave( ktide(jh) )%nh   &
292            &          + zomega_p * Wave( ktide(jh) )%np   &
293            &          + zomega_p1* Wave( ktide(jh) )%np1  ) * zscale
294      END DO
295      !
296   END SUBROUTINE tide_pulse
297
298
299   SUBROUTINE tide_vuf( pvt, put, pcor, ktide ,kc )
300      !!----------------------------------------------------------------------
301      !!                     ***  ROUTINE tide_vuf  ***
302      !!                     
303      !! ** Purpose : Compute nodal modulation corrections
304      !!
305      !! ** Outputs : vt: Phase of tidal potential relative to Greenwich (radians)
306      !!              ut: Phase correction u due to nodal motion (radians)
307      !!              ft: Nodal correction factor
308      !!----------------------------------------------------------------------
309      INTEGER                , INTENT(in ) ::   kc               ! Total number of tidal constituents
310      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide            ! Indice of tidal constituents
311      REAL(wp), DIMENSION(kc), INTENT(out) ::   pvt, put, pcor   !
312      !
313      INTEGER ::   jh   ! dummy loop index
314      !!----------------------------------------------------------------------
315      !
316      DO jh = 1, kc
317         !  Phase of the tidal potential relative to the Greenwhich
318         !  meridian (e.g. the position of the fictuous celestial body). Units are radian:
319         pvt(jh) = sh_T * Wave( ktide(jh) )%nT    &
320            &    + sh_s * Wave( ktide(jh) )%ns    &
321            &    + sh_h * Wave( ktide(jh) )%nh    &
322            &    + sh_p * Wave( ktide(jh) )%np    &
323            &    + sh_p1* Wave( ktide(jh) )%np1   &
324            &    +        Wave( ktide(jh) )%shift * rad
325         !
326         !  Phase correction u due to nodal motion. Units are radian:
327         put(jh) = sh_xi     * Wave( ktide(jh) )%nksi   &
328            &    + sh_nu     * Wave( ktide(jh) )%nnu0   &
329            &    + sh_nuprim * Wave( ktide(jh) )%nnu1   &
330            &    + sh_nusec  * Wave( ktide(jh) )%nnu2   &
331            &    + sh_R      * Wave( ktide(jh) )%R
332
333         !  Nodal correction factor:
334         pcor(jh) = nodal_factort( Wave( ktide(jh) )%nformula )
335      END DO
336      !
337   END SUBROUTINE tide_vuf
338
339
340   RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf )
341      !!----------------------------------------------------------------------
342      !!----------------------------------------------------------------------
343      INTEGER, INTENT(in) :: kformula
344      !
345      REAL(wp) :: zf
346      REAL(wp) :: zs, zf1, zf2
347      !!----------------------------------------------------------------------
348      !
349      SELECT CASE( kformula )
350      !
351      CASE( 0 )                  !==  formule 0, solar waves
352         zf = 1.0
353         !
354      CASE( 1 )                  !==  formule 1, compound waves (78 x 78)
355         zf=nodal_factort(78)
356         zf = zf * zf
357         !
358      CASE ( 2 )                 !==  formule 2, compound waves (78 x 0)  ===  (78)
359       zf1= nodal_factort(78)
360       zf = nodal_factort( 0)
361       zf = zf1 * zf
362       !
363      CASE ( 4 )                 !==  formule 4,  compound waves (78 x 235)
364         zf1 = nodal_factort( 78)
365         zf  = nodal_factort(235)
366         zf  = zf1 * zf
367         !
368      CASE ( 5 )                 !==  formule 5,  compound waves (78 *78 x 235)
369         zf1 = nodal_factort( 78)
370         zf  = nodal_factort(235)
371         zf  = zf * zf1 * zf1
372         !
373      CASE ( 6 )                 !==  formule 6,  compound waves (78 *78 x 0)
374         zf1 = nodal_factort(78)
375         zf  = nodal_factort( 0)
376         zf  = zf * zf1 * zf1 
377         !
378      CASE( 7 )                  !==  formule 7, compound waves (75 x 75)
379         zf = nodal_factort(75)
380         zf = zf * zf
381         !
382      CASE( 8 )                  !==  formule 8,  compound waves (78 x 0 x 235)
383         zf  = nodal_factort( 78)
384         zf1 = nodal_factort(  0)
385         zf2 = nodal_factort(235)
386         zf  = zf * zf1 * zf2
387         !
388      CASE( 9 )                  !==  formule 9,  compound waves (78 x 0 x 227)
389         zf  = nodal_factort( 78)
390         zf1 = nodal_factort(  0)
391         zf2 = nodal_factort(227)
392         zf  = zf * zf1 * zf2
393         !
394      CASE( 10 )                 !==  formule 10,  compound waves (78 x 227)
395         zf  = nodal_factort( 78)
396         zf1 = nodal_factort(227)
397         zf  = zf * zf1
398         !
399      CASE( 11 )                 !==  formule 11,  compound waves (75 x 0)
400!!gm bug???? zf 2 fois !
401         zf = nodal_factort(75)
402         zf1 = nodal_factort( 0)
403         zf = zf * zf1
404         !
405      CASE( 12 )                 !==  formule 12,  compound waves (78 x 78 x 78 x 0)
406         zf1 = nodal_factort(78)
407         zf  = nodal_factort( 0)
408         zf  = zf * zf1 * zf1 * zf1
409         !
410      CASE( 13 )                 !==  formule 13, compound waves (78 x 75)
411         zf1 = nodal_factort(78)
412         zf  = nodal_factort(75)
413         zf  = zf * zf1
414         !
415      CASE( 14 )                 !==  formule 14, compound waves (235 x 0)  ===  (235)
416         zf  = nodal_factort(235)
417         zf1 = nodal_factort(  0)
418         zf  = zf * zf1
419         !
420      CASE( 15 )                 !==  formule 15, compound waves (235 x 75)
421         zf  = nodal_factort(235)
422         zf1 = nodal_factort( 75)
423         zf  = zf * zf1
424         !
425      CASE( 16 )                 !==  formule 16, compound waves (78 x 0 x 0)  ===  (78)
426         zf  = nodal_factort(78)
427         zf1 = nodal_factort( 0)
428         zf  = zf * zf1 * zf1
429         !
430      CASE( 17 )                 !==  formule 17,  compound waves (227 x 0)
431         zf1 = nodal_factort(227)
432         zf  = nodal_factort(  0)
433         zf  = zf * zf1
434         !
435      CASE( 18 )                 !==  formule 18,  compound waves (78 x 78 x 78 )
436         zf1 = nodal_factort(78)
437         zf  = zf1 * zf1 * zf1
438         !
439      CASE( 19 )                 !==  formule 19, compound waves (78 x 0 x 0 x 0)  ===  (78)
440!!gm bug2 ==>>>   here identical to formule 16,  a third multiplication by zf1 is missing
441         zf  = nodal_factort(78)
442         zf1 = nodal_factort( 0)
443         zf = zf * zf1 * zf1
444         !
445      CASE( 73 )                 !==  formule 73
446         zs = sin(sh_I)
447         zf = (2./3.-zs*zs)/0.5021
448         !
449      CASE( 74 )                 !==  formule 74
450         zs = sin(sh_I)
451         zf = zs * zs / 0.1578
452         !
453      CASE( 75 )                 !==  formule 75
454         zs = cos(sh_I/2)
455         zf = sin(sh_I) * zs * zs / 0.3800
456         !
457      CASE( 76 )                 !==  formule 76
458         zf = sin(2*sh_I) / 0.7214
459         !
460      CASE( 77 )                 !==  formule 77
461         zs = sin(sh_I/2)
462         zf = sin(sh_I) * zs * zs / 0.0164
463         !
464      CASE( 78 )                 !==  formule 78
465         zs = cos(sh_I/2)
466         zf = zs * zs * zs * zs / 0.9154
467         !
468      CASE( 79 )                 !==  formule 79
469         zs = sin(sh_I)
470         zf = zs * zs / 0.1565
471         !
472      CASE( 144 )                !==  formule 144
473         zs = sin(sh_I/2)
474         zf = ( 1-10*zs*zs+15*zs*zs*zs*zs ) * cos(sh_I/2) / 0.5873
475         !
476      CASE( 149 )                !==  formule 149
477         zs = cos(sh_I/2)
478         zf = zs*zs*zs*zs*zs*zs / 0.8758
479         !
480      CASE( 215 )                !==  formule 215
481         zs = cos(sh_I/2)
482         zf = zs*zs*zs*zs / 0.9154 * sh_x1ra
483         !
484      CASE( 227 )                !==  formule 227
485         zs = sin(2*sh_I)
486         zf = sqrt( 0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006 )
487         !
488      CASE ( 235 )               !==  formule 235
489         zs = sin(sh_I)
490         zf = sqrt( 19.0444*zs*zs*zs*zs + 2.7702*zs*zs*cos(2*sh_nu) + .0981 )
491         !
492      END SELECT
493      !
494   END FUNCTION nodal_factort
495
496
497   FUNCTION dayjul( kyr, kmonth, kday )
498      !!----------------------------------------------------------------------
499      !!  *** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE)
500      !!----------------------------------------------------------------------
501      INTEGER,INTENT(in) ::   kyr, kmonth, kday
502      !
503      INTEGER,DIMENSION(12) ::  idayt, idays
504      INTEGER  ::   inc, ji
505      REAL(wp) ::   dayjul, zyq
506      !
507      DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
508      !!----------------------------------------------------------------------
509      !
510      idays(1) = 0.
511      idays(2) = 31.
512      inc = 0.
513      zyq = MOD( kyr-1900. , 4. )
514      IF( zyq == 0.)   inc = 1.
515      DO ji = 3, 12
516         idays(ji)=idayt(ji)+inc
517      END DO
518      dayjul = idays(kmonth) + kday
519      !
520   END FUNCTION dayjul
521
522   !!======================================================================
523END MODULE tide_mod
Note: See TracBrowser for help on using the repository browser.