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_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/DEV_r1879_FCM/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 2013

Last change on this file since 2013 was 2013, checked in by smasson, 14 years ago

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