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

Last change on this file since 10852 was 10852, checked in by smueller, 2 years ago

Renaming of subroutine sbc_tide of module sbctide to tide_update, transfer of this subroutine to module tide_mod, removal of the emptied module sbctide, and related adjustments (ticket #2194)

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