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 @ 10773

Last change on this file since 10773 was 10773, checked in by smueller, 5 years ago

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