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

Last change on this file since 2093 was 2093, checked in by davestorkey, 14 years ago

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