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

source: branches/2011/DEV_r2739_STFC_dCSE/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 4436

Last change on this file since 4436 was 4436, checked in by trackstand2, 10 years ago

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