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 branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: branches/UKMO/dev_r5518_GO6_under_ice_relax_dr_hook/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90 @ 11738

Last change on this file since 11738 was 11738, checked in by marc, 5 years ago

The Dr Hook changes from my perl code.
File size: 18.1 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 dom_oce        ! ocean space and time domain
9   USE phycst         ! physical constant
10   USE daymod         ! calendar
11
12   USE yomhook, ONLY: lhook, dr_hook
13   USE parkind1, ONLY: jprb, jpim
14
15   IMPLICIT NONE
16   PRIVATE
17
18   PUBLIC   tide_harmo       ! called by tideini and diaharm modules
19   PUBLIC   tide_init_Wave   ! called by tideini and diaharm modules
20
21   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 19   !: maximum number of harmonic
22
23   TYPE, PUBLIC ::    tide
24      CHARACTER(LEN=4) ::   cname_tide
25      REAL(wp)         ::   equitide
26      INTEGER          ::   nutide
27      INTEGER          ::   nt, ns, nh, np, np1, shift
28      INTEGER          ::   nksi, nnu0, nnu1, nnu2, R
29      INTEGER          ::   nformula
30   END TYPE tide
31
32   TYPE(tide), PUBLIC, DIMENSION(jpmax_harmo) ::   Wave   !:
33
34   REAL(wp) ::   sh_T, sh_s, sh_h, sh_p, sh_p1             ! astronomic angles
35   REAL(wp) ::   sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R   !
36   REAL(wp) ::   sh_I, sh_x1ra, sh_N                       !
37
38   !!----------------------------------------------------------------------
39   !! NEMO/OPA 3.3 , LOCEAN-IPSL (2010)
40   !! $Id$
41   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
42   !!----------------------------------------------------------------------
43CONTAINS
44
45   SUBROUTINE tide_init_Wave
46   INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
47   INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
48   REAL(KIND=jprb)               :: zhook_handle
49
50   CHARACTER(LEN=*), PARAMETER :: RoutineName='TIDE_INIT_WAVE'
51
52   IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
53
54#     include "tide.h90"
55   IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
56   END SUBROUTINE tide_init_Wave
57
58
59   SUBROUTINE tide_harmo( pomega, pvt, put , pcor, ktide ,kc)
60      !!----------------------------------------------------------------------
61      !!----------------------------------------------------------------------
62      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide            ! Indice of tidal constituents
63      INTEGER                , INTENT(in ) ::   kc               ! Total number of tidal constituents
64      REAL(wp), DIMENSION(kc), INTENT(out) ::   pomega           ! pulsation in radians/s
65      REAL(wp), DIMENSION(kc), INTENT(out) ::   pvt, put, pcor   !
66      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
67      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
68      REAL(KIND=jprb)               :: zhook_handle
69
70      CHARACTER(LEN=*), PARAMETER :: RoutineName='TIDE_HARMO'
71
72      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
73
74      !!----------------------------------------------------------------------
75      !
76      CALL astronomic_angle
77      CALL tide_pulse( pomega, ktide ,kc )
78      CALL tide_vuf  ( pvt, put, pcor, ktide ,kc )
79      !
80      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
81   END SUBROUTINE tide_harmo
82
83
84   SUBROUTINE astronomic_angle
85      !!----------------------------------------------------------------------
86      !!  tj is time elapsed since 1st January 1900, 0 hour, counted in julian
87      !!  century (e.g. time in days divide by 36525)
88      !!----------------------------------------------------------------------
89      REAL(wp) ::   cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2
90      REAL(wp) ::   zqy , zsy, zday, zdj, zhfrac
91      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
92      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
93      REAL(KIND=jprb)               :: zhook_handle
94
95      CHARACTER(LEN=*), PARAMETER :: RoutineName='ASTRONOMIC_ANGLE'
96
97      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
98
99      !!----------------------------------------------------------------------
100      !
101      zqy = AINT( (nyear-1901.)/4. )
102      zsy = nyear - 1900.
103      !
104      zdj  = dayjul( nyear, nmonth, nday )
105      zday = zdj + zqy - 1.
106      !
107      zhfrac = nsec_day / 3600.
108      !
109      !----------------------------------------------------------------------
110      !  Sh_n Longitude of ascending lunar node
111      !----------------------------------------------------------------------
112      sh_N=(259.1560564-19.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*rad
113      !----------------------------------------------------------------------
114      ! T mean solar angle (Greenwhich time)
115      !----------------------------------------------------------------------
116      sh_T=(180.+zhfrac*(360./24.))*rad
117      !----------------------------------------------------------------------
118      ! h mean solar Longitude
119      !----------------------------------------------------------------------
120      sh_h=(280.1895014-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad
121      !----------------------------------------------------------------------
122      ! s mean lunar Longitude
123      !----------------------------------------------------------------------
124      sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad
125      !----------------------------------------------------------------------
126      ! p1 Longitude of solar perigee
127      !----------------------------------------------------------------------
128      sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad
129      !----------------------------------------------------------------------
130      ! p Longitude of lunar perigee
131      !----------------------------------------------------------------------
132      sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad
133
134      sh_N = MOD( sh_N ,2*rpi )
135      sh_s = MOD( sh_s ,2*rpi )
136      sh_h = MOD( sh_h, 2*rpi )
137      sh_p = MOD( sh_p, 2*rpi )
138      sh_p1= MOD( sh_p1,2*rpi )
139
140      cosI = 0.913694997 -0.035692561 *cos(sh_N)
141
142      sh_I = ACOS( cosI )
143
144      sin2I   = sin(sh_I)
145      sh_tgn2 = tan(sh_N/2.0)
146
147      at1=atan(1.01883*sh_tgn2)
148      at2=atan(0.64412*sh_tgn2)
149
150      sh_xi=-at1-at2+sh_N
151
152      IF( sh_N > rpi )   sh_xi=sh_xi-2.0*rpi
153
154      sh_nu = at1 - at2
155
156      !----------------------------------------------------------------------
157      ! For constituents l2 k1 k2
158      !----------------------------------------------------------------------
159
160      tgI2 = tan(sh_I/2.0)
161      P1   = sh_p-sh_xi
162
163      t2 = tgI2*tgI2
164      t4 = t2*t2
165      sh_x1ra = sqrt( 1.0-12.0*t2*cos(2.0*P1)+36.0*t4 )
166
167      p = sin(2.0*P1)
168      q = 1.0/(6.0*t2)-cos(2.0*P1)
169      sh_R = atan(p/q)
170
171      p = sin(2.0*sh_I)*sin(sh_nu)
172      q = sin(2.0*sh_I)*cos(sh_nu)+0.3347
173      sh_nuprim = atan(p/q)
174
175      s2 = sin(sh_I)*sin(sh_I)
176      p  = s2*sin(2.0*sh_nu)
177      q  = s2*cos(2.0*sh_nu)+0.0727
178      sh_nusec = 0.5*atan(p/q)
179      !
180      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
181   END SUBROUTINE astronomic_angle
182
183
184   SUBROUTINE tide_pulse( pomega, ktide ,kc )
185      !!----------------------------------------------------------------------
186      !!                     ***  ROUTINE tide_pulse  ***
187      !!                     
188      !! ** Purpose : Compute tidal frequencies
189      !!----------------------------------------------------------------------
190      INTEGER                , INTENT(in ) ::   kc       ! Total number of tidal constituents
191      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide    ! Indice of tidal constituents
192      REAL(wp), DIMENSION(kc), INTENT(out) ::   pomega   ! pulsation in radians/s
193      !
194      INTEGER  ::   jh
195      REAL(wp) ::   zscale
196      REAL(wp) ::   zomega_T =  13149000.0_wp
197      REAL(wp) ::   zomega_s =    481267.892_wp
198      REAL(wp) ::   zomega_h =     36000.76892_wp
199      REAL(wp) ::   zomega_p =      4069.0322056_wp
200      REAL(wp) ::   zomega_n =      1934.1423972_wp
201      REAL(wp) ::   zomega_p1=         1.719175_wp
202      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
203      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
204      REAL(KIND=jprb)               :: zhook_handle
205
206      CHARACTER(LEN=*), PARAMETER :: RoutineName='TIDE_PULSE'
207
208      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
209
210      !!----------------------------------------------------------------------
211      !
212      zscale =  rad / ( 36525._wp * 86400._wp ) 
213      !
214      DO jh = 1, kc
215         pomega(jh) = (  zomega_T * Wave( ktide(jh) )%nT   &
216            &          + zomega_s * Wave( ktide(jh) )%ns   &
217            &          + zomega_h * Wave( ktide(jh) )%nh   &
218            &          + zomega_p * Wave( ktide(jh) )%np   &
219            &          + zomega_p1* Wave( ktide(jh) )%np1  ) * zscale
220      END DO
221      !
222      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
223   END SUBROUTINE tide_pulse
224
225
226   SUBROUTINE tide_vuf( pvt, put, pcor, ktide ,kc )
227      !!----------------------------------------------------------------------
228      !!                     ***  ROUTINE tide_vuf  ***
229      !!                     
230      !! ** Purpose : Compute nodal modulation corrections
231      !!
232      !! ** Outputs : vt: Phase of tidal potential relative to Greenwich (radians)
233      !!              ut: Phase correction u due to nodal motion (radians)
234      !!              ft: Nodal correction factor
235      !!----------------------------------------------------------------------
236      INTEGER                , INTENT(in ) ::   kc               ! Total number of tidal constituents
237      INTEGER , DIMENSION(kc), INTENT(in ) ::   ktide            ! Indice of tidal constituents
238      REAL(wp), DIMENSION(kc), INTENT(out) ::   pvt, put, pcor   !
239      !
240      INTEGER ::   jh   ! dummy loop index
241      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
242      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
243      REAL(KIND=jprb)               :: zhook_handle
244
245      CHARACTER(LEN=*), PARAMETER :: RoutineName='TIDE_VUF'
246
247      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
248
249      !!----------------------------------------------------------------------
250      !
251      DO jh = 1, kc
252         !  Phase of the tidal potential relative to the Greenwhich
253         !  meridian (e.g. the position of the fictuous celestial body). Units are radian:
254         pvt(jh) = sh_T * Wave( ktide(jh) )%nT    &
255            &    + sh_s * Wave( ktide(jh) )%ns    &
256            &    + sh_h * Wave( ktide(jh) )%nh    &
257            &    + sh_p * Wave( ktide(jh) )%np    &
258            &    + sh_p1* Wave( ktide(jh) )%np1   &
259            &    +        Wave( ktide(jh) )%shift * rad
260         !
261         !  Phase correction u due to nodal motion. Units are radian:
262         put(jh) = sh_xi     * Wave( ktide(jh) )%nksi   &
263            &    + sh_nu     * Wave( ktide(jh) )%nnu0   &
264            &    + sh_nuprim * Wave( ktide(jh) )%nnu1   &
265            &    + sh_nusec  * Wave( ktide(jh) )%nnu2   &
266            &    + sh_R      * Wave( ktide(jh) )%R
267
268         !  Nodal correction factor:
269         pcor(jh) = nodal_factort( Wave( ktide(jh) )%nformula )
270      END DO
271      !
272      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
273   END SUBROUTINE tide_vuf
274
275
276   RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf )
277      !!----------------------------------------------------------------------
278      !!----------------------------------------------------------------------
279      INTEGER, INTENT(in) :: kformula
280      !
281      REAL(wp) :: zf
282      REAL(wp) :: zs, zf1, zf2
283      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
284      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
285      REAL(KIND=jprb)               :: zhook_handle
286
287      CHARACTER(LEN=*), PARAMETER :: RoutineName='NODAL_FACTORT'
288
289      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
290
291      !!----------------------------------------------------------------------
292      !
293      SELECT CASE( kformula )
294      !
295      CASE( 0 )                  !==  formule 0, solar waves
296         zf = 1.0
297         !
298      CASE( 1 )                  !==  formule 1, compound waves (78 x 78)
299         zf=nodal_factort(78)
300         zf = zf * zf
301         !
302      CASE ( 2 )                 !==  formule 2, compound waves (78 x 0)  ===  (78)
303       zf1= nodal_factort(78)
304       zf = nodal_factort( 0)
305       zf = zf1 * zf
306       !
307      CASE ( 4 )                 !==  formule 4,  compound waves (78 x 235)
308         zf1 = nodal_factort( 78)
309         zf  = nodal_factort(235)
310         zf  = zf1 * zf
311         !
312      CASE ( 5 )                 !==  formule 5,  compound waves (78 *78 x 235)
313         zf1 = nodal_factort( 78)
314         zf  = nodal_factort(235)
315         zf  = zf * zf1 * zf1
316         !
317      CASE ( 6 )                 !==  formule 6,  compound waves (78 *78 x 0)
318         zf1 = nodal_factort(78)
319         zf  = nodal_factort( 0)
320         zf  = zf * zf1 * zf1 
321         !
322      CASE( 7 )                  !==  formule 7, compound waves (75 x 75)
323         zf = nodal_factort(75)
324         zf = zf * zf
325         !
326      CASE( 8 )                  !==  formule 8,  compound waves (78 x 0 x 235)
327         zf  = nodal_factort( 78)
328         zf1 = nodal_factort(  0)
329         zf2 = nodal_factort(235)
330         zf  = zf * zf1 * zf2
331         !
332      CASE( 9 )                  !==  formule 9,  compound waves (78 x 0 x 227)
333         zf  = nodal_factort( 78)
334         zf1 = nodal_factort(  0)
335         zf2 = nodal_factort(227)
336         zf  = zf * zf1 * zf2
337         !
338      CASE( 10 )                 !==  formule 10,  compound waves (78 x 227)
339         zf  = nodal_factort( 78)
340         zf1 = nodal_factort(227)
341         zf  = zf * zf1
342         !
343      CASE( 11 )                 !==  formule 11,  compound waves (75 x 0)
344!!gm bug???? zf 2 fois !
345         zf = nodal_factort(75)
346         zf = nodal_factort( 0)
347         zf = zf * zf1
348         !
349      CASE( 12 )                 !==  formule 12,  compound waves (78 x 78 x 78 x 0)
350         zf1 = nodal_factort(78)
351         zf  = nodal_factort( 0)
352         zf  = zf * zf1 * zf1 * zf1
353         !
354      CASE( 13 )                 !==  formule 13, compound waves (78 x 75)
355         zf1 = nodal_factort(78)
356         zf  = nodal_factort(75)
357         zf  = zf * zf1
358         !
359      CASE( 14 )                 !==  formule 14, compound waves (235 x 0)  ===  (235)
360         zf  = nodal_factort(235)
361         zf1 = nodal_factort(  0)
362         zf  = zf * zf1
363         !
364      CASE( 15 )                 !==  formule 15, compound waves (235 x 75)
365         zf  = nodal_factort(235)
366         zf1 = nodal_factort( 75)
367         zf  = zf * zf1
368         !
369      CASE( 16 )                 !==  formule 16, compound waves (78 x 0 x 0)  ===  (78)
370         zf  = nodal_factort(78)
371         zf1 = nodal_factort( 0)
372         zf  = zf * zf1 * zf1
373         !
374      CASE( 17 )                 !==  formule 17,  compound waves (227 x 0)
375         zf1 = nodal_factort(227)
376         zf  = nodal_factort(  0)
377         zf  = zf * zf1
378         !
379      CASE( 18 )                 !==  formule 18,  compound waves (78 x 78 x 78 )
380         zf1 = nodal_factort(78)
381         zf  = zf1 * zf1 * zf1
382         !
383      CASE( 19 )                 !==  formule 19, compound waves (78 x 0 x 0 x 0)  ===  (78)
384!!gm bug2 ==>>>   here identical to formule 16,  a third multiplication by zf1 is missing
385         zf  = nodal_factort(78)
386         zf1 = nodal_factort( 0)
387         zf = zf * zf1 * zf1
388         !
389      CASE( 73 )                 !==  formule 73
390         zs = sin(sh_I)
391         zf = (2./3.-zs*zs)/0.5021
392         !
393      CASE( 74 )                 !==  formule 74
394         zs = sin(sh_I)
395         zf = zs * zs / 0.1578
396         !
397      CASE( 75 )                 !==  formule 75
398         zs = cos(sh_I/2)
399         zf = sin(sh_I) * zs * zs / 0.3800
400         !
401      CASE( 76 )                 !==  formule 76
402         zf = sin(2*sh_I) / 0.7214
403         !
404      CASE( 77 )                 !==  formule 77
405         zs = sin(sh_I/2)
406         zf = sin(sh_I) * zs * zs / 0.0164
407         !
408      CASE( 78 )                 !==  formule 78
409         zs = cos(sh_I/2)
410         zf = zs * zs * zs * zs / 0.9154
411         !
412      CASE( 79 )                 !==  formule 79
413         zs = sin(sh_I)
414         zf = zs * zs / 0.1565
415         !
416      CASE( 144 )                !==  formule 144
417         zs = sin(sh_I/2)
418         zf = ( 1-10*zs*zs+15*zs*zs*zs*zs ) * cos(sh_I/2) / 0.5873
419         !
420      CASE( 149 )                !==  formule 149
421         zs = cos(sh_I/2)
422         zf = zs*zs*zs*zs*zs*zs / 0.8758
423         !
424      CASE( 215 )                !==  formule 215
425         zs = cos(sh_I/2)
426         zf = zs*zs*zs*zs / 0.9154 * sh_x1ra
427         !
428      CASE( 227 )                !==  formule 227
429         zs = sin(2*sh_I)
430         zf = sqrt( 0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006 )
431         !
432      CASE ( 235 )               !==  formule 235
433         zs = sin(sh_I)
434         zf = sqrt( 19.0444*zs*zs*zs*zs + 2.7702*zs*zs*cos(2*sh_nu) + .0981 )
435         !
436      END SELECT
437      !
438      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
439   END FUNCTION nodal_factort
440
441
442   FUNCTION dayjul( kyr, kmonth, kday )
443      !!----------------------------------------------------------------------
444      !!  *** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE)
445      !!----------------------------------------------------------------------
446      INTEGER,INTENT(in) ::   kyr, kmonth, kday
447      !
448      INTEGER,DIMENSION(12) ::  idayt, idays
449      INTEGER  ::   inc, ji
450      REAL(wp) ::   dayjul, zyq
451      !
452      DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
453      INTEGER(KIND=jpim), PARAMETER :: zhook_in = 0
454      INTEGER(KIND=jpim), PARAMETER :: zhook_out = 1
455      REAL(KIND=jprb)               :: zhook_handle
456
457      CHARACTER(LEN=*), PARAMETER :: RoutineName='DAYJUL'
458
459      IF (lhook) CALL dr_hook(RoutineName,zhook_in,zhook_handle)
460
461      !!----------------------------------------------------------------------
462      !
463      idays(1) = 0.
464      idays(2) = 31.
465      inc = 0.
466      zyq = MOD( kyr-1900. , 4. )
467      IF( zyq == 0.)   inc = 1.
468      DO ji = 3, 12
469         idays(ji)=idayt(ji)+inc
470      END DO
471      dayjul = idays(kmonth) + kday
472      !
473      IF (lhook) CALL dr_hook(RoutineName,zhook_out,zhook_handle)
474   END FUNCTION dayjul
475
476   !!======================================================================
477END MODULE tide_mod
Note: See TracBrowser for help on using the repository browser.