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 branches/DEV_r1986_BDY_updates/NEMO/OPA_SRC/BDY – NEMO

source: branches/DEV_r1986_BDY_updates/NEMO/OPA_SRC/BDY/bdytides.F90 @ 2168

Last change on this file since 2168 was 2168, checked in by rblod, 14 years ago

Cosmetic changes on BDY branch

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