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

source: branches/NERC/dev_r3874_FASTNEt/NEMOGCM/NEMO/OPA_SRC/SBC/tide_mod.F90 @ 6736

Last change on this file since 6736 was 6736, checked in by jamesharle, 8 years ago

FASTNEt code modifications

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