source: NEMO/branches/UKMO/NEMO_4.0_surge/src/OCE/SBC/tide_mod.F90 @ 11180

Last change on this file since 11180 was 11180, checked in by clne, 22 months ago

Initial commit of code for 2d (surge) work in NEMO4.
This is aiming to replicate the 3.6 version in branches/UKMO/dev_r5518_Surge_Modelling

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