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 @ 11764

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

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