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

Last change on this file since 12045 was 12045, checked in by smueller, 10 months ago

Modifications related to coding style and conventions, and addition of a reference in module tide_mod (ticket #2194)

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