source: NEMO/branches/UKMO/r12083_India_uncoupled/src/OCE/SBC/tide_mod.F90 @ 12453

Last change on this file since 12453 was 12453, checked in by jcastill, 12 months ago

First implementation of the branch - compiling after merge

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