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

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

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