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/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE – NEMO

source: NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/TDE/tide_mod.F90 @ 10839

Last change on this file since 10839 was 10839, checked in by smueller, 5 years ago

Adjustments to enable the use of subroutine tide_init_harmonics for the initialisation of harmonic analysis of a number of tidal constituents that differs from the number of constituents selected for tidal forcing (ticket #2194)
  • Property svn:keywords set to Id
File size: 30.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 oce            ! ocean dynamics and tracers variables
9   USE dom_oce        ! ocean space and time domain
10   USE phycst         ! physical constant
11   USE daymod         ! calendar
12   !
13   USE in_out_manager ! I/O units
14   USE iom            ! xIOs server
15   USE ioipsl         ! NetCDF IPSL library
16   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
17
18   IMPLICIT NONE
19   PRIVATE
20
21   PUBLIC   tide_init
22   PUBLIC   tide_harmo       ! called by tideini and diaharm modules
23   PUBLIC   tide_init_components ! called internally and by module diaharm
24   PUBLIC   tide_init_harmonics  ! called internally and by module diaharm
25   PUBLIC   tide_init_load
26   PUBLIC   tide_init_potential
27   PUBLIC   upd_tide         ! called in dynspg_... modules
28
29   INTEGER, PUBLIC, PARAMETER ::   jpmax_harmo = 64   !: maximum number of harmonic components
30
31   TYPE, PUBLIC ::    tide
32      CHARACTER(LEN=4) ::   cname_tide = ''
33      REAL(wp)         ::   equitide
34      INTEGER          ::   nutide
35      INTEGER          ::   nt, ns, nh, np, np1, shift
36      INTEGER          ::   nksi, nnu0, nnu1, nnu2, R
37      INTEGER          ::   nformula
38   END TYPE tide
39
40   TYPE(tide), PUBLIC, DIMENSION(:), POINTER ::   tide_components !: Array of selected tidal component parameters
41
42   TYPE, PUBLIC ::   tide_harmonic       !:   Oscillation parameters of harmonic tidal components
43      CHARACTER(LEN=4) ::   cname_tide   !    Name of component
44      REAL(wp)         ::   equitide     !    Amplitude of equilibrium tide
45      REAL(wp)         ::   f            !    Node factor
46      REAL(wp)         ::   omega        !    Angular velocity
47      REAL(wp)         ::   v0           !    Initial phase at prime meridian
48      REAL(wp)         ::   u            !    Phase correction
49   END type tide_harmonic
50
51   TYPE(tide_harmonic), PUBLIC, DIMENSION(:), POINTER ::   tide_harmonics !: Oscillation parameters of selected tidal components
52
53   LOGICAL , PUBLIC ::   ln_tide         !:
54   LOGICAL , PUBLIC ::   ln_tide_pot     !:
55   LOGICAL , PUBLIC ::   ln_read_load    !:
56   LOGICAL , PUBLIC ::   ln_scal_load    !:
57   LOGICAL , PUBLIC ::   ln_tide_ramp    !:
58   INTEGER , PUBLIC ::   nb_harmo        !: Number of active tidal components
59   INTEGER , PUBLIC ::   kt_tide         !:
60   REAL(wp), PUBLIC ::   rn_tide_ramp_dt     !:
61   REAL(wp), PUBLIC ::   rn_scal_load    !:
62   CHARACTER(lc), PUBLIC ::   cn_tide_load   !:
63   REAL(wp)         ::   rn_tide_gamma   ! Tidal tilt factor
64
65   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:)   ::   pot_astro !: tidal potential
66   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   amp_pot, phi_pot
67   REAL(wp), PUBLIC, ALLOCATABLE, DIMENSION(:,:,:) ::   amp_load, phi_load
68
69
70   REAL(wp) ::   sh_T, sh_s, sh_h, sh_p, sh_p1             ! astronomic angles
71   REAL(wp) ::   sh_xi, sh_nu, sh_nuprim, sh_nusec, sh_R   !
72   REAL(wp) ::   sh_I, sh_x1ra, sh_N                       !
73
74   !!----------------------------------------------------------------------
75   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
76   !! $Id$
77   !! Software governed by the CeCILL license (see ./LICENSE)
78   !!----------------------------------------------------------------------
79CONTAINS
80
81   SUBROUTINE tide_init
82      !!----------------------------------------------------------------------
83      !!                 ***  ROUTINE tide_init  ***
84      !!----------------------------------------------------------------------     
85      INTEGER  :: ji, jk
86      CHARACTER(LEN=4), DIMENSION(jpmax_harmo) :: sn_tide_cnames ! Names of selected tidal components
87      INTEGER  ::   ios                 ! Local integer output status for namelist read
88      !
89      NAMELIST/nam_tide/ln_tide, ln_tide_pot, rn_tide_gamma, ln_scal_load, ln_read_load, cn_tide_load, &
90                  &     ln_tide_ramp, rn_scal_load, rn_tide_ramp_dt, sn_tide_cnames
91      !!----------------------------------------------------------------------
92      !
93      ! Initialise all array elements of sn_tide_cnames, as some of them
94      ! typically do not appear in namelist_ref or namelist_cfg
95      sn_tide_cnames(:) = ''
96      ! Read Namelist nam_tide
97      REWIND( numnam_ref )              ! Namelist nam_tide in reference namelist : Tides
98      READ  ( numnam_ref, nam_tide, IOSTAT = ios, ERR = 901)
99901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nam_tide in reference namelist', lwp )
100      !
101      REWIND( numnam_cfg )              ! Namelist nam_tide in configuration namelist : Tides
102      READ  ( numnam_cfg, nam_tide, IOSTAT = ios, ERR = 902 )
103902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nam_tide in configuration namelist', lwp )
104      IF(lwm) WRITE ( numond, nam_tide )
105      !
106      IF( ln_tide ) THEN
107         IF (lwp) THEN
108            WRITE(numout,*)
109            WRITE(numout,*) 'tide_init : Initialization of the tidal components'
110            WRITE(numout,*) '~~~~~~~~~ '
111            WRITE(numout,*) '   Namelist nam_tide'
112            WRITE(numout,*) '      Use tidal components                       ln_tide         = ', ln_tide
113            WRITE(numout,*) '         Apply astronomical potential            ln_tide_pot     = ', ln_tide_pot
114            WRITE(numout,*) '         Tidal tilt factor                       rn_tide_gamma   = ', rn_tide_gamma
115            WRITE(numout,*) '         Use scalar approx. for load potential   ln_scal_load    = ', ln_scal_load
116            WRITE(numout,*) '         Read load potential from file           ln_read_load    = ', ln_read_load
117            WRITE(numout,*) '         Apply ramp on tides at startup          ln_tide_ramp    = ', ln_tide_ramp
118            WRITE(numout,*) '         Fraction of SSH used in scal. approx.   rn_scal_load    = ', rn_scal_load
119            WRITE(numout,*) '         Duration (days) of ramp                 rn_tide_ramp_dt = ', rn_tide_ramp_dt
120         ENDIF
121      ELSE
122         rn_scal_load = 0._wp 
123
124         IF(lwp) WRITE(numout,*)
125         IF(lwp) WRITE(numout,*) 'tide_init : tidal components not used (ln_tide = F)'
126         IF(lwp) WRITE(numout,*) '~~~~~~~~~ '
127         RETURN
128      ENDIF
129      !
130      ! Initialise array of selected tidal components
131      CALL tide_init_components(sn_tide_cnames, tide_components)
132      ! Number of active tidal components
133      nb_harmo = size(tide_components)
134      !       
135      ! Ensure that tidal components have been set in namelist_cfg
136      IF( nb_harmo == 0 )   CALL ctl_stop( 'tide_init : No tidal components set in nam_tide' )
137      !
138      IF( ln_read_load.AND.(.NOT.ln_tide_pot) ) &
139          &   CALL ctl_stop('ln_read_load requires ln_tide_pot')
140      IF( ln_scal_load.AND.(.NOT.ln_tide_pot) ) &
141          &   CALL ctl_stop('ln_scal_load requires ln_tide_pot')
142      IF( ln_scal_load.AND.ln_read_load ) &
143          &   CALL ctl_stop('Choose between ln_scal_load and ln_read_load')
144      IF( ln_tide_ramp.AND.((nitend-nit000+1)*rdt/rday < rn_tide_ramp_dt) )   &
145         &   CALL ctl_stop('rn_tide_ramp_dt must be lower than run duration')
146      IF( ln_tide_ramp.AND.(rn_tide_ramp_dt<0.) ) &
147         &   CALL ctl_stop('rn_tide_ramp_dt must be positive')
148      !
149      ! Initialise array used to store tidal oscillation parameters (frequency,
150      ! amplitude, phase)
151      CALL tide_init_harmonics(tide_components, tide_harmonics)
152      !
153      ! Reference time step for time-dependent tidal parameters
154      kt_tide = nit000
155      !
156      IF (.NOT.ln_scal_load ) rn_scal_load = 0._wp
157      !
158      ALLOCATE( amp_pot(jpi,jpj,nb_harmo),                      &
159           &      phi_pot(jpi,jpj,nb_harmo), pot_astro(jpi,jpj)   )
160      IF( ln_read_load ) THEN
161         ALLOCATE( amp_load(jpi,jpj,nb_harmo), phi_load(jpi,jpj,nb_harmo) )
162      ENDIF
163      !
164   END SUBROUTINE tide_init
165
166
167   SUBROUTINE tide_init_components(pcnames, ptide_comp)
168      !!----------------------------------------------------------------------
169      !!                 ***  ROUTINE tide_init_components  ***
170      !!
171      !! Returns pointer to array of variables of type 'tide' that contain
172      !! information about the selected tidal components
173      !! ----------------------------------------------------------------------
174      CHARACTER(LEN=4),              DIMENSION(jpmax_harmo), INTENT(in)  ::   pcnames         ! Names of selected components
175      TYPE(tide),       POINTER,     DIMENSION(:),           INTENT(out) ::   ptide_comp      ! Selected components
176      INTEGER,          ALLOCATABLE, DIMENSION(:)                        ::   kcomppos        ! Indices of selected components
177      INTEGER                                                            ::   kcomp, jk, ji   ! Miscellaneous integers
178      TYPE(tide),       POINTER,     DIMENSION(:)                        ::   tide_components ! All available components
179     
180      ! Populate local array with information about all available tidal
181      ! components
182      !
183      ! Note, here 'tide_components' locally overrides the global module
184      ! variable of the same name to enable the use of the global name in the
185      ! include file that contains the initialisation of elements of array
186      ! 'tide_components'
187      ALLOCATE(tide_components(jpmax_harmo), kcomppos(jpmax_harmo))
188      ! Initialise array of indices of the selected componenents
189      kcomppos(:) = 0
190      ! Include tidal component parameters for all available components
191#include "tide.h90"
192     
193      ! Identify the selected components that are availble
194      kcomp = 0
195      DO jk = 1, jpmax_harmo
196         IF (TRIM(pcnames(jk)) /= '') THEN
197            DO ji = 1, jpmax_harmo
198               IF (TRIM(pcnames(jk)) == tide_components(ji)%cname_tide) THEN
199                  kcomp = kcomp + 1
200                  WRITE(numout, '(10X,"Tidal component #",I2.2,36X,"= ",A4)') kcomp, pcnames(jk)
201                  kcomppos(kcomp) = ji
202                  EXIT
203               END IF
204            END DO
205         END IF
206      END DO
207     
208      ! Allocate and populate reduced list of components
209      ALLOCATE(ptide_comp(kcomp))
210      DO jk = 1, kcomp
211         ptide_comp(jk) = tide_components(kcomppos(jk))
212      END DO
213     
214      ! Release local array of available components and list of selected
215      ! components
216      DEALLOCATE(tide_components, kcomppos)
217     
218   END SUBROUTINE tide_init_components
219
220
221   SUBROUTINE tide_init_harmonics(ptide_comp, ptide_harmo)
222      !!----------------------------------------------------------------------
223      !!                 ***  ROUTINE tide_init_harmonics  ***
224      !!
225      !! Returns pointer to array of variables of type 'tide_harmonics' that
226      !! contain oscillation parameters of the selected harmonic tidal
227      !! components
228      !! ----------------------------------------------------------------------
229      TYPE(tide),          POINTER,     DIMENSION(:) ::   ptide_comp      ! Selected components
230      TYPE(tide_harmonic), POINTER,     DIMENSION(:) ::   ptide_harmo     ! Oscillation parameters of selected components
231
232      ! Allocate and populate array of oscillation parameters
233      ALLOCATE(ptide_harmo(size(ptide_comp)))
234      ptide_harmo(:)%cname_tide = ptide_comp(:)%cname_tide
235      ptide_harmo(:)%equitide = ptide_comp(:)%equitide
236      CALL tide_harmo(ptide_comp, ptide_harmo)
237
238   END SUBROUTINE tide_init_harmonics
239
240
241   SUBROUTINE tide_init_potential
242      !!----------------------------------------------------------------------
243      !!                 ***  ROUTINE tide_init_potential  ***
244      !!----------------------------------------------------------------------
245      INTEGER  ::   ji, jj, jk   ! dummy loop indices
246      REAL(wp) ::   zcons, ztmp1, ztmp2, zlat, zlon, ztmp, zamp, zcs   ! local scalar
247      !!----------------------------------------------------------------------
248
249      DO jk = 1, nb_harmo
250         zcons = rn_tide_gamma * tide_components(jk)%equitide * tide_harmonics(jk)%f
251         DO ji = 1, jpi
252            DO jj = 1, jpj
253               ztmp1 =  tide_harmonics(jk)%f * amp_pot(ji,jj,jk) * COS( phi_pot(ji,jj,jk) &
254                  &                                                   + tide_harmonics(jk)%v0 + tide_harmonics(jk)%u )
255               ztmp2 = -tide_harmonics(jk)%f * amp_pot(ji,jj,jk) * SIN( phi_pot(ji,jj,jk) &
256                  &                                                   + tide_harmonics(jk)%v0 + tide_harmonics(jk)%u )
257               zlat = gphit(ji,jj)*rad !! latitude en radian
258               zlon = glamt(ji,jj)*rad !! longitude en radian
259               ztmp = tide_harmonics(jk)%v0 + tide_harmonics(jk)%u + tide_components(jk)%nutide * zlon
260               ! le potentiel est composé des effets des astres:
261               IF    ( tide_components(jk)%nutide == 1 )  THEN  ;  zcs = zcons * SIN( 2._wp*zlat )
262               ELSEIF( tide_components(jk)%nutide == 2 )  THEN  ;  zcs = zcons * COS( zlat )**2
263               ELSE                                         ;  zcs = 0._wp
264               ENDIF
265               ztmp1 = ztmp1 + zcs * COS( ztmp )
266               ztmp2 = ztmp2 - zcs * SIN( ztmp )
267               zamp = SQRT( ztmp1*ztmp1 + ztmp2*ztmp2 )
268               amp_pot(ji,jj,jk) = zamp
269               phi_pot(ji,jj,jk) = ATAN2( -ztmp2 / MAX( 1.e-10_wp , zamp ) ,   &
270                  &                        ztmp1 / MAX( 1.e-10_wp,  zamp )   )
271            END DO
272         END DO
273      END DO
274      !
275   END SUBROUTINE tide_init_potential
276
277
278   SUBROUTINE tide_init_load
279      !!----------------------------------------------------------------------
280      !!                 ***  ROUTINE tide_init_load  ***
281      !!----------------------------------------------------------------------
282      INTEGER :: inum                 ! Logical unit of input file
283      INTEGER :: ji, jj, itide        ! dummy loop indices
284      REAL(wp), DIMENSION(jpi,jpj) ::   ztr, zti   !: workspace to read in tidal harmonics data
285      !!----------------------------------------------------------------------
286      IF(lwp) THEN
287         WRITE(numout,*)
288         WRITE(numout,*) 'tide_init_load : Initialization of load potential from file'
289         WRITE(numout,*) '~~~~~~~~~~~~~~ '
290      ENDIF
291      !
292      CALL iom_open ( cn_tide_load , inum )
293      !
294      DO itide = 1, nb_harmo
295         CALL iom_get  ( inum, jpdom_data,TRIM(tide_components(itide)%cname_tide)//'_z1', ztr(:,:) )
296         CALL iom_get  ( inum, jpdom_data,TRIM(tide_components(itide)%cname_tide)//'_z2', zti(:,:) )
297         !
298         DO ji=1,jpi
299            DO jj=1,jpj
300               amp_load(ji,jj,itide) =  SQRT( ztr(ji,jj)**2. + zti(ji,jj)**2. )
301               phi_load(ji,jj,itide) = ATAN2(-zti(ji,jj), ztr(ji,jj) )
302            END DO
303         END DO
304         !
305      END DO
306      CALL iom_close( inum )
307      !
308   END SUBROUTINE tide_init_load
309
310
311   SUBROUTINE tide_harmo( ptide_comp, ptide_harmo )
312      !
313      TYPE(tide),          DIMENSION(:), POINTER ::   ptide_comp   ! Array of selected tidal component parameters
314      TYPE(tide_harmonic), DIMENSION(:), POINTER ::   ptide_harmo  ! Oscillation parameters of selected tidal components
315      !
316      CALL astronomic_angle
317      CALL tide_pulse( ptide_comp, ptide_harmo )
318      CALL tide_vuf(   ptide_comp, ptide_harmo )
319      !
320   END SUBROUTINE tide_harmo
321
322
323   SUBROUTINE astronomic_angle
324      !!----------------------------------------------------------------------
325      !!  tj is time elapsed since 1st January 1900, 0 hour, counted in julian
326      !!  century (e.g. time in days divide by 36525)
327      !!----------------------------------------------------------------------
328      REAL(wp) ::   cosI, p, q, t2, t4, sin2I, s2, tgI2, P1, sh_tgn2, at1, at2
329      REAL(wp) ::   zqy , zsy, zday, zdj, zhfrac
330      !!----------------------------------------------------------------------
331      !
332      zqy = AINT( (nyear-1901.)/4. )
333      zsy = nyear - 1900.
334      !
335      zdj  = dayjul( nyear, nmonth, nday )
336      zday = zdj + zqy - 1.
337      !
338      zhfrac = nsec_day / 3600.
339      !
340      !----------------------------------------------------------------------
341      !  Sh_n Longitude of ascending lunar node
342      !----------------------------------------------------------------------
343      sh_N=(259.1560564-19.328185764*zsy-.0529539336*zday-.0022064139*zhfrac)*rad
344      !----------------------------------------------------------------------
345      ! T mean solar angle (Greenwhich time)
346      !----------------------------------------------------------------------
347      sh_T=(180.+zhfrac*(360./24.))*rad
348      !----------------------------------------------------------------------
349      ! h mean solar Longitude
350      !----------------------------------------------------------------------
351      sh_h=(280.1895014-.238724988*zsy+.9856473288*zday+.0410686387*zhfrac)*rad
352      !----------------------------------------------------------------------
353      ! s mean lunar Longitude
354      !----------------------------------------------------------------------
355      sh_s=(277.0256206+129.38482032*zsy+13.176396768*zday+.549016532*zhfrac)*rad
356      !----------------------------------------------------------------------
357      ! p1 Longitude of solar perigee
358      !----------------------------------------------------------------------
359      sh_p1=(281.2208569+.01717836*zsy+.000047064*zday+.000001961*zhfrac)*rad
360      !----------------------------------------------------------------------
361      ! p Longitude of lunar perigee
362      !----------------------------------------------------------------------
363      sh_p=(334.3837214+40.66246584*zsy+.111404016*zday+.004641834*zhfrac)*rad
364
365      sh_N = MOD( sh_N ,2*rpi )
366      sh_s = MOD( sh_s ,2*rpi )
367      sh_h = MOD( sh_h, 2*rpi )
368      sh_p = MOD( sh_p, 2*rpi )
369      sh_p1= MOD( sh_p1,2*rpi )
370
371      cosI = 0.913694997 -0.035692561 *cos(sh_N)
372
373      sh_I = ACOS( cosI )
374
375      sin2I   = sin(sh_I)
376      sh_tgn2 = tan(sh_N/2.0)
377
378      at1=atan(1.01883*sh_tgn2)
379      at2=atan(0.64412*sh_tgn2)
380
381      sh_xi=-at1-at2+sh_N
382
383      IF( sh_N > rpi )   sh_xi=sh_xi-2.0*rpi
384
385      sh_nu = at1 - at2
386
387      !----------------------------------------------------------------------
388      ! For constituents l2 k1 k2
389      !----------------------------------------------------------------------
390
391      tgI2 = tan(sh_I/2.0)
392      P1   = sh_p-sh_xi
393
394      t2 = tgI2*tgI2
395      t4 = t2*t2
396      sh_x1ra = sqrt( 1.0-12.0*t2*cos(2.0*P1)+36.0*t4 )
397
398      p = sin(2.0*P1)
399      q = 1.0/(6.0*t2)-cos(2.0*P1)
400      sh_R = atan(p/q)
401
402      p = sin(2.0*sh_I)*sin(sh_nu)
403      q = sin(2.0*sh_I)*cos(sh_nu)+0.3347
404      sh_nuprim = atan(p/q)
405
406      s2 = sin(sh_I)*sin(sh_I)
407      p  = s2*sin(2.0*sh_nu)
408      q  = s2*cos(2.0*sh_nu)+0.0727
409      sh_nusec = 0.5*atan(p/q)
410      !
411   END SUBROUTINE astronomic_angle
412
413
414   SUBROUTINE tide_pulse( ptide_comp, ptide_harmo )
415      !!----------------------------------------------------------------------
416      !!                     ***  ROUTINE tide_pulse  ***
417      !!                     
418      !! ** Purpose : Compute tidal frequencies
419      !!----------------------------------------------------------------------
420      TYPE(tide),          DIMENSION(:), POINTER ::   ptide_comp   ! Array of selected tidal component parameters
421      TYPE(tide_harmonic), DIMENSION(:), POINTER ::   ptide_harmo  ! Oscillation parameters of selected tidal components
422      !
423      INTEGER  ::   jh
424      REAL(wp) ::   zscale
425      REAL(wp) ::   zomega_T =  13149000.0_wp
426      REAL(wp) ::   zomega_s =    481267.892_wp
427      REAL(wp) ::   zomega_h =     36000.76892_wp
428      REAL(wp) ::   zomega_p =      4069.0322056_wp
429      REAL(wp) ::   zomega_n =      1934.1423972_wp
430      REAL(wp) ::   zomega_p1=         1.719175_wp
431      !!----------------------------------------------------------------------
432      !
433      zscale =  rad / ( 36525._wp * 86400._wp ) 
434      !
435      DO jh = 1, size(ptide_harmo)
436         ptide_harmo(jh)%omega = (  zomega_T * ptide_comp( jh )%nT   &
437            &                         + zomega_s * ptide_comp( jh )%ns   &
438            &                         + zomega_h * ptide_comp( jh )%nh   &
439            &                         + zomega_p * ptide_comp( jh )%np   &
440            &                         + zomega_p1* ptide_comp( jh )%np1  ) * zscale
441      END DO
442      !
443   END SUBROUTINE tide_pulse
444
445
446   SUBROUTINE tide_vuf( ptide_comp, ptide_harmo )
447      !!----------------------------------------------------------------------
448      !!                     ***  ROUTINE tide_vuf  ***
449      !!                     
450      !! ** Purpose : Compute nodal modulation corrections
451      !!
452      !! ** Outputs : vt: Phase of tidal potential relative to Greenwich (radians)
453      !!              ut: Phase correction u due to nodal motion (radians)
454      !!              ft: Nodal correction factor
455      !!----------------------------------------------------------------------
456      TYPE(tide),          DIMENSION(:), POINTER ::   ptide_comp   ! Array of selected tidal component parameters
457      TYPE(tide_harmonic), DIMENSION(:), POINTER ::   ptide_harmo  ! Oscillation parameters of selected tidal components
458      !
459      INTEGER ::   jh   ! dummy loop index
460      !!----------------------------------------------------------------------
461      !
462      DO jh = 1, size(ptide_harmo)
463         !  Phase of the tidal potential relative to the Greenwhich
464         !  meridian (e.g. the position of the fictuous celestial body). Units are radian:
465         ptide_harmo(jh)%v0 = sh_T * ptide_comp( jh )%nT    &
466            &                   + sh_s * ptide_comp( jh )%ns    &
467            &                   + sh_h * ptide_comp( jh )%nh    &
468            &                   + sh_p * ptide_comp( jh )%np    &
469            &                   + sh_p1* ptide_comp( jh )%np1   &
470            &                   +        ptide_comp( jh )%shift * rad
471         !
472         !  Phase correction u due to nodal motion. Units are radian:
473         ptide_harmo(jh)%u = sh_xi     * ptide_comp( jh )%nksi   &
474            &                  + sh_nu     * ptide_comp( jh )%nnu0   &
475            &                  + sh_nuprim * ptide_comp( jh )%nnu1   &
476            &                  + sh_nusec  * ptide_comp( jh )%nnu2   &
477            &                  + sh_R      * ptide_comp( jh )%R
478
479         !  Nodal correction factor:
480         ptide_harmo(jh)%f = nodal_factort( ptide_comp( jh )%nformula )
481      END DO
482      !
483   END SUBROUTINE tide_vuf
484
485
486   RECURSIVE FUNCTION nodal_factort( kformula ) RESULT( zf )
487      !!----------------------------------------------------------------------
488      !!----------------------------------------------------------------------
489      INTEGER, INTENT(in) :: kformula
490      !
491      REAL(wp) :: zf
492      REAL(wp) :: zs, zf1, zf2
493      !!----------------------------------------------------------------------
494      !
495      SELECT CASE( kformula )
496      !
497      CASE( 0 )                  !==  formule 0, solar waves
498         zf = 1.0
499         !
500      CASE( 1 )                  !==  formule 1, compound waves (78 x 78)
501         zf=nodal_factort(78)
502         zf = zf * zf
503         !
504      CASE ( 2 )                 !==  formule 2, compound waves (78 x 0)  ===  (78)
505       zf1= nodal_factort(78)
506       zf = nodal_factort( 0)
507       zf = zf1 * zf
508       !
509      CASE ( 4 )                 !==  formule 4,  compound waves (78 x 235)
510         zf1 = nodal_factort( 78)
511         zf  = nodal_factort(235)
512         zf  = zf1 * zf
513         !
514      CASE ( 5 )                 !==  formule 5,  compound waves (78 *78 x 235)
515         zf1 = nodal_factort( 78)
516         zf  = nodal_factort(235)
517         zf  = zf * zf1 * zf1
518         !
519      CASE ( 6 )                 !==  formule 6,  compound waves (78 *78 x 0)
520         zf1 = nodal_factort(78)
521         zf  = nodal_factort( 0)
522         zf  = zf * zf1 * zf1 
523         !
524      CASE( 7 )                  !==  formule 7, compound waves (75 x 75)
525         zf = nodal_factort(75)
526         zf = zf * zf
527         !
528      CASE( 8 )                  !==  formule 8,  compound waves (78 x 0 x 235)
529         zf  = nodal_factort( 78)
530         zf1 = nodal_factort(  0)
531         zf2 = nodal_factort(235)
532         zf  = zf * zf1 * zf2
533         !
534      CASE( 9 )                  !==  formule 9,  compound waves (78 x 0 x 227)
535         zf  = nodal_factort( 78)
536         zf1 = nodal_factort(  0)
537         zf2 = nodal_factort(227)
538         zf  = zf * zf1 * zf2
539         !
540      CASE( 10 )                 !==  formule 10,  compound waves (78 x 227)
541         zf  = nodal_factort( 78)
542         zf1 = nodal_factort(227)
543         zf  = zf * zf1
544         !
545      CASE( 11 )                 !==  formule 11,  compound waves (75 x 0)
546!!gm bug???? zf 2 fois !
547         zf = nodal_factort(75)
548         zf1 = nodal_factort( 0)
549         zf = zf * zf1
550         !
551      CASE( 12 )                 !==  formule 12,  compound waves (78 x 78 x 78 x 0)
552         zf1 = nodal_factort(78)
553         zf  = nodal_factort( 0)
554         zf  = zf * zf1 * zf1 * zf1
555         !
556      CASE( 13 )                 !==  formule 13, compound waves (78 x 75)
557         zf1 = nodal_factort(78)
558         zf  = nodal_factort(75)
559         zf  = zf * zf1
560         !
561      CASE( 14 )                 !==  formule 14, compound waves (235 x 0)  ===  (235)
562         zf  = nodal_factort(235)
563         zf1 = nodal_factort(  0)
564         zf  = zf * zf1
565         !
566      CASE( 15 )                 !==  formule 15, compound waves (235 x 75)
567         zf  = nodal_factort(235)
568         zf1 = nodal_factort( 75)
569         zf  = zf * zf1
570         !
571      CASE( 16 )                 !==  formule 16, compound waves (78 x 0 x 0)  ===  (78)
572         zf  = nodal_factort(78)
573         zf1 = nodal_factort( 0)
574         zf  = zf * zf1 * zf1
575         !
576      CASE( 17 )                 !==  formule 17,  compound waves (227 x 0)
577         zf1 = nodal_factort(227)
578         zf  = nodal_factort(  0)
579         zf  = zf * zf1
580         !
581      CASE( 18 )                 !==  formule 18,  compound waves (78 x 78 x 78 )
582         zf1 = nodal_factort(78)
583         zf  = zf1 * zf1 * zf1
584         !
585      CASE( 19 )                 !==  formule 19, compound waves (78 x 0 x 0 x 0)  ===  (78)
586!!gm bug2 ==>>>   here identical to formule 16,  a third multiplication by zf1 is missing
587         zf  = nodal_factort(78)
588         zf1 = nodal_factort( 0)
589         zf = zf * zf1 * zf1
590         !
591      CASE( 73 )                 !==  formule 73
592         zs = sin(sh_I)
593         zf = (2./3.-zs*zs)/0.5021
594         !
595      CASE( 74 )                 !==  formule 74
596         zs = sin(sh_I)
597         zf = zs * zs / 0.1578
598         !
599      CASE( 75 )                 !==  formule 75
600         zs = cos(sh_I/2)
601         zf = sin(sh_I) * zs * zs / 0.3800
602         !
603      CASE( 76 )                 !==  formule 76
604         zf = sin(2*sh_I) / 0.7214
605         !
606      CASE( 77 )                 !==  formule 77
607         zs = sin(sh_I/2)
608         zf = sin(sh_I) * zs * zs / 0.0164
609         !
610      CASE( 78 )                 !==  formule 78
611         zs = cos(sh_I/2)
612         zf = zs * zs * zs * zs / 0.9154
613         !
614      CASE( 79 )                 !==  formule 79
615         zs = sin(sh_I)
616         zf = zs * zs / 0.1565
617         !
618      CASE( 144 )                !==  formule 144
619         zs = sin(sh_I/2)
620         zf = ( 1-10*zs*zs+15*zs*zs*zs*zs ) * cos(sh_I/2) / 0.5873
621         !
622      CASE( 149 )                !==  formule 149
623         zs = cos(sh_I/2)
624         zf = zs*zs*zs*zs*zs*zs / 0.8758
625         !
626      CASE( 215 )                !==  formule 215
627         zs = cos(sh_I/2)
628         zf = zs*zs*zs*zs / 0.9154 * sh_x1ra
629         !
630      CASE( 227 )                !==  formule 227
631         zs = sin(2*sh_I)
632         zf = sqrt( 0.8965*zs*zs+0.6001*zs*cos (sh_nu)+0.1006 )
633         !
634      CASE ( 235 )               !==  formule 235
635         zs = sin(sh_I)
636         zf = sqrt( 19.0444*zs*zs*zs*zs + 2.7702*zs*zs*cos(2*sh_nu) + .0981 )
637         !
638      END SELECT
639      !
640   END FUNCTION nodal_factort
641
642
643   FUNCTION dayjul( kyr, kmonth, kday )
644      !!----------------------------------------------------------------------
645      !!  *** THIS ROUTINE COMPUTES THE JULIAN DAY (AS A REAL VARIABLE)
646      !!----------------------------------------------------------------------
647      INTEGER,INTENT(in) ::   kyr, kmonth, kday
648      !
649      INTEGER,DIMENSION(12) ::  idayt, idays
650      INTEGER  ::   inc, ji
651      REAL(wp) ::   dayjul, zyq
652      !
653      DATA idayt/0.,31.,59.,90.,120.,151.,181.,212.,243.,273.,304.,334./
654      !!----------------------------------------------------------------------
655      !
656      idays(1) = 0.
657      idays(2) = 31.
658      inc = 0.
659      zyq = MOD( kyr-1900. , 4. )
660      IF( zyq == 0.)   inc = 1.
661      DO ji = 3, 12
662         idays(ji)=idayt(ji)+inc
663      END DO
664      dayjul = idays(kmonth) + kday
665      !
666   END FUNCTION dayjul
667
668
669   SUBROUTINE upd_tide( kt, kit, time_offset )
670      !!----------------------------------------------------------------------
671      !!                 ***  ROUTINE upd_tide  ***
672      !!
673      !! ** Purpose :   provide at each time step the astronomical potential
674      !!
675      !! ** Method  :   computed from pulsation and amplitude of all tide components
676      !!
677      !! ** Action  :   pot_astro   actronomical potential
678      !!----------------------------------------------------------------------     
679      INTEGER, INTENT(in)           ::   kt      ! ocean time-step index
680      INTEGER, INTENT(in), OPTIONAL ::   kit     ! external mode sub-time-step index (lk_dynspg_ts=T)
681      INTEGER, INTENT(in), OPTIONAL ::   time_offset ! time offset in number
682                                                     ! of internal steps             (lk_dynspg_ts=F)
683                                                     ! of external steps             (lk_dynspg_ts=T)
684      !
685      INTEGER  ::   joffset      ! local integer
686      INTEGER  ::   ji, jj, jk   ! dummy loop indices
687      REAL(wp) ::   zt, zramp    ! local scalar
688      REAL(wp), DIMENSION(nb_harmo) ::   zwt 
689      !!----------------------------------------------------------------------     
690      !
691      !                               ! tide pulsation at model time step (or sub-time-step)
692      zt = ( kt - kt_tide ) * rdt
693      !
694      joffset = 0
695      IF( PRESENT( time_offset ) )   joffset = time_offset
696      !
697      IF( PRESENT( kit ) )   THEN
698         zt = zt + ( kit +  joffset - 1 ) * rdt / REAL( nn_baro, wp )
699      ELSE
700         zt = zt + joffset * rdt
701      ENDIF
702      !
703      zwt(:) = tide_harmonics(:)%omega * zt
704
705      pot_astro(:,:) = 0._wp          ! update tidal potential (sum of all harmonics)
706      DO jk = 1, nb_harmo   
707         pot_astro(:,:) = pot_astro(:,:) + amp_pot(:,:,jk) * COS( zwt(jk) + phi_pot(:,:,jk) )     
708      END DO
709      !
710      IF( ln_tide_ramp ) THEN         ! linear increase if asked
711         zt = ( kt - nit000 ) * rdt
712         IF( PRESENT( kit ) )   zt = zt + ( kit + joffset -1) * rdt / REAL( nn_baro, wp )
713         zramp = MIN(  MAX( zt / (rn_tide_ramp_dt*rday) , 0._wp ) , 1._wp  )
714         pot_astro(:,:) = zramp * pot_astro(:,:)
715      ENDIF
716      !
717   END SUBROUTINE upd_tide
718
719   !!======================================================================
720END MODULE tide_mod
Note: See TracBrowser for help on using the repository browser.