source: NEMO/branches/2020/KERNEL-03_Storkey_Coward_RK3_stage2/src/OCE/TDE/tide_mod.F90 @ 12443

Last change on this file since 12443 was 12443, checked in by davestorkey, 9 months ago

2020/KERNEL-03_Storkey_Coward_RK3_stage2: More variable renaming:
atfp → rn_atfp (use namelist parameter everywhere)
rdtbt → rDt_e
nn_baro → nn_e
rn_scal_load → rn_load
rau0 → rho0

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