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.
bdytides.F90 in trunk/NEMO/OPA_SRC/BDY – NEMO

source: trunk/NEMO/OPA_SRC/BDY/bdytides.F90 @ 1125

Last change on this file since 1125 was 1125, checked in by ctlod, 16 years ago

trunk: BDY package code review (coding rules), see ticket: #214

  • Property svn:executable set to *
File size: 28.9 KB
Line 
1MODULE bdytides
2   !!======================================================================
3   !!                       ***  MODULE  bdytides  ***
4   !! Ocean dynamics:   Tidal forcing at open boundaries
5   !!======================================================================
6   !! History :  2.0  !  2007-01  (D.Storkey)  Original code
7   !!            2.3  !  2008-01  (J.Holt)  Add date correction. Origins POLCOMS v6.3 2007
8   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
9   !!----------------------------------------------------------------------
10#if defined key_bdy
11   !!----------------------------------------------------------------------
12   !!   'key_bdy'     Unstructured Open Boundary Condition
13   !!----------------------------------------------------------------------
14   !!   PUBLIC
15   !!      tide_init     : read of namelist
16   !!      tide_data     : read in and initialisation of tidal constituents at boundary
17   !!      tide_update   : calculation of tidal forcing at each timestep
18   !!   PRIVATE
19   !!      uvset         :\
20   !!      vday          : | Routines to correct tidal harmonics forcing for
21   !!      shpen         : | start time of integration
22   !!      ufset         : |
23   !!      vset          :/
24   !!----------------------------------------------------------------------
25   USE oce             ! ocean dynamics and tracers
26   USE dom_oce         ! ocean space and time domain
27   USE iom
28   USE in_out_manager  ! I/O units
29   USE phycst          ! physical constants
30   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
31   USE bdy_par         ! Unstructured boundary parameters
32   USE bdy_oce         ! ocean open boundary conditions
33   USE daymod          ! calendar
34
35   IMPLICIT NONE
36   PRIVATE
37
38   PUBLIC   tide_init     ! routine called in bdyini
39   PUBLIC   tide_data     ! routine called in bdyini
40   PUBLIC   tide_update   ! routine called in bdydyn
41
42   LOGICAL, PUBLIC            ::   ln_tide_date            !: =T correct tide phases and amplitude for model start date
43
44   INTEGER, PARAMETER ::   jptides_max = 15      !: Max number of tidal contituents
45   INTEGER            ::   ntide                 !: Actual number of tidal constituents
46
47   CHARACTER(len=80), PUBLIC                         ::   filtide    !: Filename root for tidal input files
48   CHARACTER(len= 4), PUBLIC, DIMENSION(jptides_max) ::   tide_cpt   !: Names of tidal components used.
49
50   INTEGER , DIMENSION(jptides_max) ::   nindx        !: ???
51   REAL(wp), DIMENSION(jptides_max) ::   tide_speed   !: Phase speed of tidal constituent (deg/hr)
52   
53   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   ssh1, ssh2   !: Tidal constituents : SSH
54   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   u1  , u2     !: Tidal constituents : U
55   REAL(wp), DIMENSION(jpbdim,jptides_max)  ::   v1  , v2     !: Tidal constituents : V
56   
57   !!----------------------------------------------------------------------
58   !! NEMO/OPA 3.0 , LOCEAN-IPSL (2008)
59   !! $Id: $
60   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
61   !!----------------------------------------------------------------------
62
63CONTAINS
64
65   SUBROUTINE tide_init
66      !!----------------------------------------------------------------------
67      !!                    ***  SUBROUTINE tide_init  ***
68      !!                     
69      !! ** Purpose : - Read in namelist for tides
70      !!
71      !!----------------------------------------------------------------------
72      INTEGER ::   itide                  ! dummy loop index
73      !!
74      NAMELIST/namtide/ln_tide_date, filtide, tide_cpt, tide_speed
75      !!----------------------------------------------------------------------
76
77      IF(lwp) WRITE(numout,*)
78      IF(lwp) WRITE(numout,*) 'tide_init : initialization of tidal harmonic forcing at open boundaries'
79      IF(lwp) WRITE(numout,*) '~~~~~~~~~'
80
81      ! Namelist namtide : tidal harmonic forcing at open boundaries
82      ln_tide_date = .false.
83      filtide(:) = ''
84      tide_speed(:) = 0.0
85      tide_cpt(:) = ''
86      REWIND( numnam )                                ! Read namelist parameters
87      READ  ( numnam, namtide )
88      !                                               ! Count number of components specified
89      ntide = jptides_max
90      itide = 1
91      DO WHILE( tide_cpt(itide) /= '' )
92         ntide = itide
93         itide = itide + 1
94      END DO
95      !                                               ! find constituents in standard list
96      DO itide = 1, ntide
97         nindx(itide) = 0
98         IF( TRIM( tide_cpt(itide) ) == 'Q1'  )   nindx(itide) =  1
99         IF( TRIM( tide_cpt(itide) ) == 'O1'  )   nindx(itide) =  2
100         IF( TRIM( tide_cpt(itide) ) == 'P1'  )   nindx(itide) =  3
101         IF( TRIM( tide_cpt(itide) ) == 'S1'  )   nindx(itide) =  4
102         IF( TRIM( tide_cpt(itide) ) == 'K1'  )   nindx(itide) =  5
103         IF( TRIM( tide_cpt(itide) ) == '2N2' )   nindx(itide) =  6
104         IF( TRIM( tide_cpt(itide) ) == 'MU2' )   nindx(itide) =  7
105         IF( TRIM( tide_cpt(itide) ) == 'N2'  )   nindx(itide) =  8
106         IF( TRIM( tide_cpt(itide) ) == 'NU2' )   nindx(itide) =  9
107         IF( TRIM( tide_cpt(itide) ) == 'M2'  )   nindx(itide) = 10
108         IF( TRIM( tide_cpt(itide) ) == 'L2'  )   nindx(itide) = 11
109         IF( TRIM( tide_cpt(itide) ) == 'T2'  )   nindx(itide) = 12
110         IF( TRIM( tide_cpt(itide) ) == 'S2'  )   nindx(itide) = 13
111         IF( TRIM( tide_cpt(itide) ) == 'K2'  )   nindx(itide) = 14
112         IF( TRIM( tide_cpt(itide) ) == 'M4'  )   nindx(itide) = 15
113         IF( nindx(itide) == 0  .AND. lwp ) THEN
114            WRITE(ctmp1,*) 'constitunent', itide,':', tide_cpt(itide), 'not in standard list'
115            CALL ctl_warn( ctmp1 )
116         ENDIF
117      END DO
118      !                                               ! Parameter control and print
119      IF( ntide < 1 ) THEN
120         CALL ctl_stop( '          Did not find any tidal components in namelist namtide' )
121      ELSE
122         IF(lwp) WRITE(numout,*) '          Namelist namtide : tidal harmonic forcing at open boundaries'
123         IF(lwp) WRITE(numout,*) '             tidal components specified ', ntide
124         IF(lwp) WRITE(numout,*) '                ', tide_cpt(1:ntide)
125         IF(lwp) WRITE(numout,*) '             associated phase speeds (deg/hr) : '
126         IF(lwp) WRITE(numout,*) '                ', tide_speed(1:ntide)
127      ENDIF
128
129      ! Initialisation of tidal harmonics arrays
130      sshtide(:) = 0.e0
131      utide  (:) = 0.e0
132      vtide  (:) = 0.e0
133      !
134   END SUBROUTINE tide_init
135
136
137   SUBROUTINE tide_data
138      !!----------------------------------------------------------------------
139      !!                    ***  SUBROUTINE tide_data  ***
140      !!                   
141      !! ** Purpose : - Read in tidal harmonics data and adjust for the start
142      !!                time of the model run.
143      !!
144      !!----------------------------------------------------------------------
145      INTEGER ::   itide, igrd, ib        ! dummy loop indices
146      CHARACTER(len=80) :: clfile         ! full file name for tidal input file
147      INTEGER ::   ipi, ipj, inum, idvar  ! temporary integers (netcdf read)
148      INTEGER, DIMENSION(3) :: lendta     ! length of data in the file (note may be different from nblendta!)
149      REAL(wp) ::  z_arg, z_atde, z_btde, z1t, z2t           
150      REAL(wp), DIMENSION(jpbdta,1) ::   zdta   ! temporary array for data fields
151      REAL(wp), DIMENSION(jptides_max) :: z_vplu, z_ftc
152      !!------------------------------------------------------------------------------
153
154      ! Open files and read in tidal forcing data
155      ! -----------------------------------------
156
157      ipj = 1
158
159      DO itide = 1, ntide
160         !                                                              ! SSH fields
161         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_T.nc'
162         IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
163         CALL iom_open( clfile, inum )
164         igrd = 1
165         IF( nblendta(igrd) <= 0 ) THEN
166            idvar = iom_varid( inum,'z1' )
167            IF(lwp) WRITE(numout,*) 'iom_file(1)%ndims(idvar) : ',iom_file%ndims(idvar)
168            nblendta(igrd) = iom_file(inum)%dimsz(1,idvar)
169            WRITE(numout,*) 'Dim size for z1 is ', nblendta(igrd)
170         ENDIF
171         ipi = nblendta(igrd)
172         CALL iom_get( inum, jpdom_unknown, 'z1', zdta(1:ipi,1:ipj) )
173         DO ib = 1, nblenrim(igrd)
174            ssh1(ib,itide) = zdta(nbmap(ib,igrd),1)
175         END DO
176         CALL iom_get( inum, jpdom_unknown, 'z2', zdta(1:ipi,1:ipj) )
177         DO ib = 1, nblenrim(igrd)
178            ssh2(ib,itide) = zdta(nbmap(ib,igrd),1)
179         END DO
180         CALL iom_close( inum )
181         !
182         !                                                              ! U fields
183         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_U.nc'
184         IF(lwp) WRITE(numout,*) 'Reading data from file ', clfile
185         CALL iom_open( clfile, inum )
186         igrd = 2
187         IF( lendta(igrd) <= 0 ) THEN
188            idvar = iom_varid( inum,'u1' )
189            lendta(igrd) = iom_file(inum)%dimsz(1,idvar)
190            WRITE(numout,*) 'Dim size for u1 is ',lendta(igrd)
191         ENDIF
192         ipi = lendta(igrd)
193         CALL iom_get( inum, jpdom_unknown, 'u1', zdta(1:ipi,1:ipj) )
194         DO ib = 1, nblenrim(igrd)
195            u1(ib,itide) = zdta(nbmap(ib,igrd),1)
196         END DO
197         CALL iom_get( inum, jpdom_unknown, 'u2', zdta(1:ipi,1:ipj) )
198         DO ib = 1, nblenrim(igrd)
199            u2(ib,itide) = zdta(nbmap(ib,igrd),1)
200         END DO
201         CALL iom_close( inum )
202         !
203         !                                                              ! V fields
204         clfile = TRIM(filtide)//TRIM(tide_cpt(itide))//'_grid_V.nc'
205         if(lwp) write(numout,*) 'Reading data from file ', clfile
206         CALL iom_open( clfile, inum )
207         igrd = 3
208         IF( lendta(igrd) <= 0 ) THEN
209            idvar = iom_varid( inum,'v1' )
210            lendta(igrd) = iom_file(inum)%dimsz(1,idvar)
211            WRITE(numout,*) 'Dim size for v1 is ', lendta(igrd)
212         ENDIF
213         ipi = lendta(igrd)
214         CALL iom_get( inum, jpdom_unknown, 'v1', zdta(1:ipi,1:ipj) )
215         DO ib = 1, nblenrim(igrd)
216            v1(ib,itide) = zdta(nbmap(ib,igrd),1)
217         END DO
218         CALL iom_get( inum, jpdom_unknown, 'v2', zdta(1:ipi,1:ipj) )
219         DO ib=1, nblenrim(igrd)
220            v2(ib,itide) = zdta(nbmap(ib,igrd),1)
221         END DO
222         CALL iom_close( inum )
223         !
224      END DO ! end loop on tidal components
225
226      IF( ln_tide_date ) THEN      ! correct for date factors
227
228!! used nmonth, nyear and nday from daymod....
229         ! Calculate date corrects for 15 standard consituents
230         ! This is the initialisation step, so nday, nmonth, nyear are the
231         ! initial date/time of the integration.
232
233         CALL uvset( 0, nday, nmonth, nyear, z_ftc, z_vplu )
234
235         IF(lwp) WRITE(numout,*) 'Correcting tide for date:', nday, nmonth, nyear
236
237         DO itide = 1, ntide       ! loop on tidal components
238            !
239            IF( nindx(itide) /= 0 ) THEN
240!!gm use rpi  and rad global variable 
241               z_arg = 3.14159265d0 * z_vplu(nindx(itide)) / 180.0d0
242               z_atde=z_ftc(nindx(itide))*cos(z_arg)
243               z_btde=z_ftc(nindx(itide))*sin(z_arg)
244               IF(lwp) WRITE(numout,'(2i5,8f10.6)') itide, nindx(itide), tide_speed(itide),   &
245                  &                                 z_ftc(nindx(itide)), z_vplu(nindx(itide))
246            ELSE
247               z_atde = 1.0_wp
248               z_btde = 0.0_wp
249            ENDIF
250            !                                         !  elevation         
251            igrd = 1
252            DO ib = 1, nblenrim(igrd)               
253               z1t = z_atde * ssh1(ib,itide) + z_btde * ssh2(ib,itide)
254               z2t = z_atde * ssh2(ib,itide) - z_btde * ssh1(ib,itide)
255               ssh1(ib,itide) = z1t
256               ssh2(ib,itide) = z2t
257            END DO
258            !                                         !  u       
259            igrd = 2
260            DO ib = 1, nblenrim(igrd)               
261               z1t = z_atde * u1(ib,itide) + z_btde * u2(ib,itide)
262               z2t = z_atde * u2(ib,itide) - z_btde * u1(ib,itide)
263               u1(ib,itide) = z1t
264               u2(ib,itide) = z2t
265            END DO
266            !                                         !  v       
267            igrd = 3
268            DO ib = 1, nblenrim(igrd)               
269               z1t = z_atde * v1(ib,itide) + z_btde * v2(ib,itide)
270               z2t = z_atde * v2(ib,itide) - z_btde * v1(ib,itide)
271               v1(ib,itide) = z1t
272               v2(ib,itide) = z2t
273            END DO
274            !
275         END DO   ! end loop on tidal components
276         !
277      ENDIF ! date correction
278      !
279   END SUBROUTINE tide_data
280
281
282   SUBROUTINE tide_update ( kt, jit )
283      !!----------------------------------------------------------------------
284      !!                 ***  SUBROUTINE tide_update  ***
285      !!               
286      !! ** Purpose : - Add tidal forcing to sshbdy, ubtbdy and vbtbdy arrays.
287      !!               
288      !!----------------------------------------------------------------------
289      INTEGER, INTENT( in ) ::   kt      ! Main timestep counter
290!!gm doctor jit ==> kit
291      INTEGER, INTENT( in ) ::   jit     ! Barotropic timestep counter (for timesplitting option)
292      !!
293      INTEGER  ::   itide, igrd, ib      ! dummy loop indices
294      REAL(wp) ::   z_arg, z_sarg            !           
295      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost
296      !!----------------------------------------------------------------------
297
298      ! Note tide phase speeds are in deg/hour, so we need to convert the
299      ! elapsed time in seconds to hours by dividing by 3600.0
300      IF( jit == 0 ) THEN 
301         z_arg = kt * rdt * rad / 3600.0
302      ELSE                              ! we are in a barotropic subcycle (for timesplitting option)
303         z_arg = ( (kt-1) * rdt + jit * rdtbt ) * rad / 3600.0
304      ENDIF
305
306      DO itide = 1, ntide
307         z_sarg = z_arg * tide_speed(itide)
308         z_cost(itide) = COS( z_sarg )
309         z_sist(itide) = SIN( z_sarg )
310      END DO
311
312      ! summing of tidal constituents into BDY arrays
313      sshtide(:) = 0.0
314      utide (:) = 0.0
315      vtide (:) = 0.0
316      !
317      DO itide = 1, ntide
318         igrd=1                              ! SSH on tracer grid.
319         DO ib = 1, nblenrim(igrd)
320            sshtide(ib) =sshtide(ib)+ ssh1(ib,itide)*z_cost(itide) + ssh2(ib,itide)*z_sist(itide)
321            !    if(lwp) write(numout,*) 'z',ib,itide,sshtide(ib), ssh1(ib,itide),ssh2(ib,itide)
322         END DO
323         igrd=2                              ! U grid
324         DO ib=1, nblenrim(igrd)
325            utide(ib) = utide(ib)+ u1(ib,itide)*z_cost(itide) + u2(ib,itide)*z_sist(itide)
326            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), u1(ib,itide),u2(ib,itide)
327         END DO
328         igrd=3                              ! V grid
329         DO ib=1, nblenrim(igrd)
330            vtide(ib) = vtide(ib)+ v1(ib,itide)*z_cost(itide) + v2(ib,itide)*z_sist(itide)
331            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), v1(ib,itide),v2(ib,itide)
332         END DO
333      END DO
334      !
335   END SUBROUTINE tide_update
336
337!!gm  doctor naming of dummy argument variables!!!   and all local variables
338
339   SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu )
340      !!----------------------------------------------------------------------
341      !!                   ***  SUBROUTINE uvset  ***
342      !!                     
343      !! ** Purpose : - adjust tidal forcing for date factors
344      !!               
345      !!----------------------------------------------------------------------
346      implicit none
347      INTEGER, INTENT( in ) ::   ihs     ! Start time hours
348      INTEGER, INTENT( in ) ::   iday    ! start time days
349      INTEGER, INTENT( in ) ::   imnth   ! start time month
350      INTEGER, INTENT( in ) ::   iyr     ! start time year
351      !!
352!!gm  nc is jptides_max ????
353      INTEGER , PARAMETER     ::  nc =15    ! maximum number of constituents
354      CHARACTER(len=8), DIMENSION(nc) :: cname
355      INTEGER                 ::   year, vd, ivdy, ndc, i, k
356      REAL(wp)                ::   ss, h, p, en, p1, rtd
357      REAL(wp), DIMENSION(nc) ::   f                          ! nodal correction
358      REAL(wp), DIMENSION(nc) ::   z_vplu                     ! phase correction
359      REAL(wp), DIMENSION(nc) ::   u, v, zig
360      !!
361      DATA cname/  'q1'    ,    'o1'    ,     'p1'   ,    's1'    ,     'k1'    ,   &
362         &         '2n2'   ,    'mu2'   ,     'n2'   ,    'nu2'   ,     'm2'    ,   &
363         &         'l2'    ,    't2'    ,     's2'   ,    'k2'    ,     'm4'      /
364      DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878,  .2625161701,   &
365         &      .4868657873, .4881373225, .4963669182, .4976384533,  .5058680490,   &
366         &      .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098   /
367      !!----------------------------------------------------------------------
368!
369!     ihs             -  start time gmt on ...
370!     iday/imnth/iyr  -  date   e.g.   12/10/87
371!
372      CALL vday(iday,imnth,iyr,ivdy)
373      vd = ivdy
374!
375!rp   note change of year number for d. blackman shpen
376!rp   if(iyr.ge.1000) year=iyr-1900
377!rp   if(iyr.lt.1000) year=iyr
378      year = iyr
379!
380!.....year  =  year of required data
381!.....vd    =  day of required data..set up for 0000gmt day year
382      ndc = nc
383!.....ndc   =  number of constituents allowed
384!!gm use rpi ?
385      rtd = 360.0 / 6.2831852
386      DO i = 1, ndc
387         zig(i) = zig(i)*rtd
388         ! sigo(i)= zig(i)
389      END DO
390
391!!gm try to avoid RETURN  in F90
392      IF( year == 0 )   RETURN
393     
394         CALL shpen( year, vd, ss, h , p , en, p1 )
395         CALL ufset( p   , en, u , f )
396         CALL vset ( ss  , h , p , en, p1, v )
397         !
398         DO k = 1, nc
399            z_vplu(k) = v(k) + u(k)
400            z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k)
401            DO WHILE( z_vplu(k) < 0    )
402               z_vplu(k) = z_vplu(k) + 360.0
403            END DO
404            DO WHILE( z_vplu(k) > 360. )
405               z_vplu(k) = z_vplu(k) - 360.0
406            END DO
407         END DO
408         !
409      END SUBROUTINE uvset
410
411
412      SUBROUTINE vday( iday, imnth, iy, ivdy )
413      !!----------------------------------------------------------------------
414      !!                   *** SUBROUTINE vday  ***
415      !!                 
416      !! ** Purpose : - adjust tidal forcing for date factors
417      !!               
418      !!----------------------------------------------------------------------
419      INTEGER, INTENT(in   ) ::   iday, imnth, iy   ! ????
420      INTEGER, INTENT(  out) ::   ivdy              ! ???
421      !!
422      INTEGER ::   iyr
423      !!------------------------------------------------------------------------------
424
425!!gm   nday_year in day mode is the variable compiuted here, no?
426!!gm      nday_year ,   &  !: curent day counted from jan 1st of the current year
427
428      !calculate day number in year from day/month/year
429       if(imnth.eq.1) ivdy=iday
430       if(imnth.eq.2) ivdy=iday+31
431       if(imnth.eq.3) ivdy=iday+59
432       if(imnth.eq.4) ivdy=iday+90
433       if(imnth.eq.5) ivdy=iday+120
434       if(imnth.eq.6) ivdy=iday+151
435       if(imnth.eq.7) ivdy=iday+181
436       if(imnth.eq.8) ivdy=iday+212
437       if(imnth.eq.9) ivdy=iday+243
438       if(imnth.eq.10) ivdy=iday+273
439       if(imnth.eq.11) ivdy=iday+304
440       if(imnth.eq.12) ivdy=iday+334
441       iyr=iy
442       if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1
443       if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1
444       if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1
445      !
446   END SUBROUTINE vday
447
448!!doctor norme for dummy arguments
449
450   SUBROUTINE shpen( year, vd, s, h, p, en, p1 )
451      !!----------------------------------------------------------------------
452      !!                    ***  SUBROUTINE shpen  ***
453      !!                   
454      !! ** Purpose : - calculate astronomical arguments for tides
455      !!                this version from d. blackman 30 nove 1990
456      !!
457      !!----------------------------------------------------------------------
458!!gm add INTENT in, out or inout....    DOCTOR name....
459!!gm please do not use variable name with a single letter:  impossible to search in a code
460      INTEGER  ::   year,vd
461      REAL(wp) ::   s,h,p,en,p1
462      !!
463      INTEGER  ::   yr,ilc,icent,it,iday,ild,ipos,nn,iyd
464      REAL(wp) ::   cycle,t,td,delt(84),delta,deltat
465      !!
466      DATA delt /-5.04, -3.90, -2.87, -0.58,  0.71,  1.80,           &
467         &        3.08,  4.63,  5.86,  7.21,  8.58, 10.50, 12.10,    &
468         &       12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39,    &
469         &       19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68,    &
470         &       22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63,    &
471         &       23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80,    &
472         &       24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89,    &
473         &       27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96,    &
474         &       31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39,    &
475         &       33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87,    &
476         &       38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00,    &
477         &       45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81,    &
478         &       52.57 /
479      !!----------------------------------------------------------------------
480
481      cycle = 360.0
482      ilc = 0
483      icent = year / 100
484      yr = year - icent * 100
485      t = icent - 20
486!
487!     for the following equations
488!     time origin is fixed at 00 hr of jan 1st,2000.
489!     see notes by cartwright
490!
491!!gm  old coding style, use CASE instead  and avoid GOTO (obsolescence in fortran 90)
492!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
493      it = icent - 20
494      if (it) 1,2,2
495    1 iday = it/4 -it
496      go to 3
497    2 iday = (it+3)/4 - it
498!
499!     t is in julian century
500!     correction in gegorian calander where only century year divisible
501!     by 4 is leap year.
502!
503    3 continue
504!
505      td = 0.0
506!
507!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
508      if (yr) 4,5,4
509!
510    4 iyd = 365*yr
511      ild = (yr-1)/4
512      if((icent - (icent/4)*4) .eq. 0) ilc = 1
513      td = iyd + ild + ilc
514!
515    5 td = td + iday + vd -1.0 - 0.5
516      t = t + (td/36525.0)
517!
518      ipos=year-1899
519      if (ipos .lt. 0) go to 7
520      if (ipos .gt. 83) go to 6
521!
522      delta = (delt(ipos+1)+delt(ipos))/2.0
523      go to 7
524!
525    6 delta= (65.0-50.5)/20.0*(year-1980)+50.5
526!
527    7 deltat = delta * 1.0e-6
528!
529!!gm   precision of the computation   :  example for s it should be replace by:
530!!gm  s   = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat   ==>  more precise  modify the last digits results
531      s   = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat
532      h   = 280.4661 + 36000.7698 *t + 0.0003*t*t +  11.0*deltat
533      p   =  83.3535 + 4069.0139  *t - 0.0103*t*t +       deltat
534      en  = 234.9555 + 1934.1363  *t - 0.0021*t*t +       deltat
535      p1  = 282.9384 + 1.7195     *t + 0.0005*t*t
536      !
537      nn = s / cycle
538      s = s - nn * cycle
539      IF( s < 0.e0 )   s = s + cycle
540      !
541      nn = h / cycle
542      h = h - cycle * nn
543      IF( h < 0.e0 )   h = h + cycle
544      !
545      nn = p / cycle
546      p = p - cycle * nn
547      IF( p < 0.e0)   p = p + cycle
548      !
549      nn = en / cycle
550      en = en - cycle * nn
551      IF( en < 0.e0 )   en = en + cycle
552      en = cycle - en
553      !
554      nn = p1 / cycle
555      p1 = p1 - nn * cycle
556      !
557   END SUBROUTINE shpen
558
559
560   SUBROUTINE ufset( p, cn, b, a )
561      !!----------------------------------------------------------------------
562      !!                    ***  SUBROUTINE ufset  ***
563      !!                   
564      !! ** Purpose : - calculate nodal parameters for the tides
565      !!               
566      !!----------------------------------------------------------------------
567!!gm  doctor naming of dummy argument variables!!!   and all local variables
568!!gm  nc is jptides_max ????
569      integer nc
570      parameter (nc=15)
571      REAL(wp) p,cn
572      !!
573     
574!!gm  rad is already a public variable defined in phycst.F90 ....   ==> doctor norme  local real start with "z"
575      REAL(wp) ::   w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad
576      REAL(wp) ::   a(nc), b(nc)
577      INTEGER  ::   ndc, k
578      !!----------------------------------------------------------------------
579
580      ndc = nc
581
582!    a=f       ,  b =u
583!    t is  zero as compared to tifa.
584!! use rad defined in phycst   (i.e.  add a USE phycst at the begining of the module
585      rad = 6.2831852d0/360.0
586      pw = p  * rad
587      nw = cn * rad
588      w1 = cos(   nw )
589      w2 = cos( 2*nw )
590      w3 = cos( 3*nw )
591      w4 = sin(   nw )
592      w5 = sin( 2*nw )
593      w6 = sin( 3*nw )
594      w7 = 1. - 0.2505 * COS( 2*pw      ) - 0.1102 * COS(2*pw-nw )   &
595         &    - 0.156  * COS( 2*pw-2*nw ) - 0.037  * COS(     nw )
596      w8 =    - 0.2505 * SIN( 2*pw      ) - 0.1102 * SIN(2*pw-nw )   &
597         &    - 0.0156 * SIN( 2*pw-2*nw ) - 0.037  * SIN(     nw )
598      !
599      a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3
600      b(1) =          0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6
601      !   q1
602      a(2) = a(1)
603      b(2) = b(1)
604      !   o1
605      a(3) = 1.0
606      b(3) = 0.0
607      !   p1
608      a(4) = 1.0
609      b(4) = 0.0
610      !   s1
611      a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3
612      b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6
613      !   k1
614      a(6) =1.0004 -0.0373*w1+ 0.0002*w2
615      b(6) = -0.0374*w4
616      !  2n2
617      a(7) = a(6)
618      b(7) = b(6)
619      !  mu2
620      a(8) = a(6)
621      b(8) = b(6)
622      !   n2
623      a(9) = a(6)
624      b(9) = b(6)
625      !  nu2
626      a(10) = a(6)
627      b(10) = b(6)
628      !   m2
629      a(11) = SQRT( w7 * w7 + w8 * w8 )
630      b(11) = ATAN( w8 / w7 )
631!!gmuse rpi  instead of 3.141992 ???   true pi is rpi=3.141592653589793_wp  .....   ????
632      IF( w7 < 0.e0 )   b(11) = b(11) + 3.141992
633      !   l2
634      a(12) = 1.0
635      b(12) = 0.0
636      !   t2
637      a(13)= a(12)
638      b(13)= b(12)
639      !   s2
640      a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3
641      b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6
642      !   k2
643      a(15) = a(6)*a(6)
644      b(15) = 2*b(6)
645      !   m4
646!!gm  old coding,  remove GOTO and label of lines
647!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
648      DO 40 k = 1,ndc
649         b(k) = b(k)/rad
65032       if (b(k)) 34,35,35
65134       b(k) = b(k) + 360.0
652         go to 32
65335       if (b(k)-360.0) 40,37,37
65437       b(k) = b(k)-360.0
655         go to 35
65640    continue
657      !
658   END SUBROUTINE ufset
659   
660
661   SUBROUTINE vset( s,h,p,en,p1,v)
662      !!----------------------------------------------------------------------
663      !!                    ***  SUBROUTINE vset  ***
664      !!                   
665      !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run
666      !!               
667      !!----------------------------------------------------------------------
668!!gm  doctor naming of dummy argument variables!!!   and all local variables
669!!gm  nc is jptides_max ????
670!!gm  en argument is not used: suppress it ?
671      integer nc
672      parameter (nc=15)
673      real(wp) s,h,p,en,p1
674      real(wp)   v(nc)
675      !!
676      integer ndc, k
677      !!----------------------------------------------------------------------
678
679      ndc = nc
680      !   v s  are computed here.
681      v(1) =-3*s +h +p +270      ! Q1
682      v(2) =-2*s +h +270.0       ! O1
683      v(3) =-h +270              ! P1
684      v(4) =180                  ! S1
685      v(5) =h +90.0              ! K1
686      v(6) =-4*s +2*h +2*p       ! 2N2
687      v(7) =-4*(s-h)             ! MU2
688      v(8) =-3*s +2*h +p         ! N2
689      v(9) =-3*s +4*h -p         ! MU2
690      v(10) =-2*s +2*h           ! M2
691      v(11) =-s +2*h -p +180     ! L2
692      v(12) =-h +p1              ! T2
693      v(13) =0                   ! S2
694      v(14) =h+h                 ! K2
695      v(15) =2*v(10)             ! M4
696!
697!!gm  old coding,  remove GOTO and label of lines
698!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
699      do 72 k = 1, ndc
70069       if( v(k) )   70,71,71
70170       v(k) = v(k)+360.0
702         go to 69
70371       if( v(k) - 360.0 )   72,73,73
70473       v(k) = v(k)-360.0
705         go to 71
70672    continue
707      !
708   END SUBROUTINE vset
709
710#else
711   !!----------------------------------------------------------------------
712   !!   Dummy module         NO Unstruct Open Boundary Conditions for tides
713   !!----------------------------------------------------------------------
714!!gm  are you sure we need to define filtide and tide_cpt ?
715   CHARACTER(len=80), PUBLIC               ::   filtide                !: Filename root for tidal input files
716   CHARACTER(len=4 ), PUBLIC, DIMENSION(1) ::   tide_cpt               !: Names of tidal components used.
717
718CONTAINS
719   SUBROUTINE tide_init                ! Empty routine
720   END SUBROUTINE tide_init
721   SUBROUTINE tide_data                ! Empty routine
722   END SUBROUTINE tide_data
723   SUBROUTINE tide_update( kt, kit )   ! Empty routine
724      WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, kit
725   END SUBROUTINE tide_update
726#endif
727
728   !!======================================================================
729END MODULE bdytides
Note: See TracBrowser for help on using the repository browser.