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 trunk/NEMOGCM/NEMO/OPA_SRC/SBC – NEMO

source: trunk/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 12 years ago

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