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

Last change on this file since 10793 was 10793, checked in by smueller, 19 months ago

Replacement of a hard-coded coefficient by namelist parameter rn_tide_gamma (ticket #2194)

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