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/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC – NEMO

source: NEMO/branches/2018/dev_r10057_ENHANCE03_ZTILDE/tests/COMODO-IW/MY_SRC/tide_mod.F90 @ 10118

Last change on this file since 10118 was 10118, checked in by jchanut, 6 years ago

ztilde update (4): #2126
Add internal wave ALE benchmark (from COMODO group)

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