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

Last change on this file since 3211 was 3211, checked in by spickles2, 12 years ago

Stephen Pickles, 11 Dec 2011

Commit to bring the rest of the DCSE NEMO development branch
in line with the latest development version. This includes
array index re-ordering of all OPA_SRC/.

  • Property svn:keywords set to Id
File size: 29.5 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      !!----------------------------------------------------------------------
293      !!                 ***  SUBROUTINE tide_update  ***
294      !!               
295      !! ** Purpose : - Add tidal forcing to sshbdy, ubtbdy and vbtbdy arrays.
296      !!               
297      !!----------------------------------------------------------------------
298      INTEGER, INTENT( in ) ::   kt      ! Main timestep counter
299!!gm doctor jit ==> kit
300      INTEGER, INTENT( in ) ::   jit     ! Barotropic timestep counter (for timesplitting option)
301      !!
302      INTEGER  ::   itide, igrd, ib      ! dummy loop indices
303      REAL(wp) ::   z_arg, z_sarg            !           
304      REAL(wp), DIMENSION(jptides_max) :: z_sist, z_cost
305      !!----------------------------------------------------------------------
306
307      ! Note tide phase speeds are in deg/hour, so we need to convert the
308      ! elapsed time in seconds to hours by dividing by 3600.0
309      IF( jit == 0 ) THEN 
310         z_arg = kt * rdt * rad / 3600.0
311      ELSE                              ! we are in a barotropic subcycle (for timesplitting option)
312!         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,lwp) ) * rad / 3600.0
313         z_arg = ( (kt-1) * rdt + jit * rdt / REAL(nn_baro,wp) ) * rad / 3600.0
314      ENDIF
315
316      DO itide = 1, ntide
317         z_sarg = z_arg * tide_speed(itide)
318         z_cost(itide) = COS( z_sarg )
319         z_sist(itide) = SIN( z_sarg )
320      END DO
321
322      ! summing of tidal constituents into BDY arrays
323      sshtide(:) = 0.0
324      utide (:) = 0.0
325      vtide (:) = 0.0
326      !
327      DO itide = 1, ntide
328         igrd=4                              ! SSH on tracer grid.
329         DO ib = 1, nblenrim(igrd)
330            sshtide(ib) =sshtide(ib)+ ssh1(ib,itide)*z_cost(itide) + ssh2(ib,itide)*z_sist(itide)
331            !    if(lwp) write(numout,*) 'z',ib,itide,sshtide(ib), ssh1(ib,itide),ssh2(ib,itide)
332         END DO
333         igrd=5                              ! U grid
334         DO ib=1, nblenrim(igrd)
335            utide(ib) = utide(ib)+ u1(ib,itide)*z_cost(itide) + u2(ib,itide)*z_sist(itide)
336            !    if(lwp) write(numout,*) 'u',ib,itide,utide(ib), u1(ib,itide),u2(ib,itide)
337         END DO
338         igrd=6                              ! V grid
339         DO ib=1, nblenrim(igrd)
340            vtide(ib) = vtide(ib)+ v1(ib,itide)*z_cost(itide) + v2(ib,itide)*z_sist(itide)
341            !    if(lwp) write(numout,*) 'v',ib,itide,vtide(ib), v1(ib,itide),v2(ib,itide)
342         END DO
343      END DO
344      !
345   END SUBROUTINE tide_update
346
347!!gm  doctor naming of dummy argument variables!!!   and all local variables
348
349   SUBROUTINE uvset( ihs, iday, imnth, iyr, f, z_vplu )
350      !!----------------------------------------------------------------------
351      !!                   ***  SUBROUTINE uvset  ***
352      !!                     
353      !! ** Purpose : - adjust tidal forcing for date factors
354      !!               
355      !!----------------------------------------------------------------------
356      implicit none
357      INTEGER, INTENT( in ) ::   ihs     ! Start time hours
358      INTEGER, INTENT( in ) ::   iday    ! start time days
359      INTEGER, INTENT( in ) ::   imnth   ! start time month
360      INTEGER, INTENT( in ) ::   iyr     ! start time year
361      !!
362!!gm  nc is jptides_max ????
363      INTEGER , PARAMETER     ::  nc =15    ! maximum number of constituents
364      CHARACTER(len=8), DIMENSION(nc) :: cname
365      INTEGER                 ::   year, vd, ivdy, ndc, i, k
366      REAL(wp)                ::   ss, h, p, en, p1, rtd
367      REAL(wp), DIMENSION(nc) ::   f                          ! nodal correction
368      REAL(wp), DIMENSION(nc) ::   z_vplu                     ! phase correction
369      REAL(wp), DIMENSION(nc) ::   u, v, zig
370      !!
371      DATA cname/  'q1'    ,    'o1'    ,     'p1'   ,    's1'    ,     'k1'    ,   &
372         &         '2n2'   ,    'mu2'   ,     'n2'   ,    'nu2'   ,     'm2'    ,   &
373         &         'l2'    ,    't2'    ,     's2'   ,    'k2'    ,     'm4'      /
374      DATA zig/ .2338507481, .2433518789, .2610826055, .2617993878,  .2625161701,   &
375         &      .4868657873, .4881373225, .4963669182, .4976384533,  .5058680490,   &
376         &      .5153691799, .5228820265, .5235987756, .5250323419, 1.011736098   /
377      !!----------------------------------------------------------------------
378!
379!     ihs             -  start time gmt on ...
380!     iday/imnth/iyr  -  date   e.g.   12/10/87
381!
382      CALL vday(iday,imnth,iyr,ivdy)
383      vd = ivdy
384!
385!rp   note change of year number for d. blackman shpen
386!rp   if(iyr.ge.1000) year=iyr-1900
387!rp   if(iyr.lt.1000) year=iyr
388      year = iyr
389!
390!.....year  =  year of required data
391!.....vd    =  day of required data..set up for 0000gmt day year
392      ndc = nc
393!.....ndc   =  number of constituents allowed
394!!gm use rpi ?
395      rtd = 360.0 / 6.2831852
396      DO i = 1, ndc
397         zig(i) = zig(i)*rtd
398         ! sigo(i)= zig(i)
399      END DO
400
401!!gm try to avoid RETURN  in F90
402      IF( year == 0 )   RETURN
403     
404         CALL shpen( year, vd, ss, h , p , en, p1 )
405         CALL ufset( p   , en, u , f )
406         CALL vset ( ss  , h , p , en, p1, v )
407         !
408         DO k = 1, nc
409            z_vplu(k) = v(k) + u(k)
410            z_vplu(k) = z_vplu(k) + dble(ihs) * zig(k)
411            DO WHILE( z_vplu(k) < 0    )
412               z_vplu(k) = z_vplu(k) + 360.0
413            END DO
414            DO WHILE( z_vplu(k) > 360. )
415               z_vplu(k) = z_vplu(k) - 360.0
416            END DO
417         END DO
418         !
419      END SUBROUTINE uvset
420
421
422      SUBROUTINE vday( iday, imnth, iy, ivdy )
423      !!----------------------------------------------------------------------
424      !!                   *** SUBROUTINE vday  ***
425      !!                 
426      !! ** Purpose : - adjust tidal forcing for date factors
427      !!               
428      !!----------------------------------------------------------------------
429      INTEGER, INTENT(in   ) ::   iday, imnth, iy   ! ????
430      INTEGER, INTENT(  out) ::   ivdy              ! ???
431      !!
432      INTEGER ::   iyr
433      !!------------------------------------------------------------------------------
434
435!!gm   nday_year in day mode is the variable compiuted here, no?
436!!gm      nday_year ,   &  !: curent day counted from jan 1st of the current year
437
438      !calculate day number in year from day/month/year
439       if(imnth.eq.1) ivdy=iday
440       if(imnth.eq.2) ivdy=iday+31
441       if(imnth.eq.3) ivdy=iday+59
442       if(imnth.eq.4) ivdy=iday+90
443       if(imnth.eq.5) ivdy=iday+120
444       if(imnth.eq.6) ivdy=iday+151
445       if(imnth.eq.7) ivdy=iday+181
446       if(imnth.eq.8) ivdy=iday+212
447       if(imnth.eq.9) ivdy=iday+243
448       if(imnth.eq.10) ivdy=iday+273
449       if(imnth.eq.11) ivdy=iday+304
450       if(imnth.eq.12) ivdy=iday+334
451       iyr=iy
452       if(mod(iyr,4).eq.0.and.imnth.gt.2) ivdy=ivdy+1
453       if(mod(iyr,100).eq.0.and.imnth.gt.2) ivdy=ivdy-1
454       if(mod(iyr,400).eq.0.and.imnth.gt.2) ivdy=ivdy+1
455      !
456   END SUBROUTINE vday
457
458!!doctor norme for dummy arguments
459
460   SUBROUTINE shpen( year, vd, s, h, p, en, p1 )
461      !!----------------------------------------------------------------------
462      !!                    ***  SUBROUTINE shpen  ***
463      !!                   
464      !! ** Purpose : - calculate astronomical arguments for tides
465      !!                this version from d. blackman 30 nove 1990
466      !!
467      !!----------------------------------------------------------------------
468!!gm add INTENT in, out or inout....    DOCTOR name....
469!!gm please do not use variable name with a single letter:  impossible to search in a code
470      INTEGER  ::   year,vd
471      REAL(wp) ::   s,h,p,en,p1
472      !!
473      INTEGER  ::   yr,ilc,icent,it,iday,ild,ipos,nn,iyd
474      REAL(wp) ::   cycle,t,td,delt(84),delta,deltat
475      !!
476      DATA delt /-5.04, -3.90, -2.87, -0.58,  0.71,  1.80,           &
477         &        3.08,  4.63,  5.86,  7.21,  8.58, 10.50, 12.10,    &
478         &       12.49, 14.41, 15.59, 15.81, 17.52, 19.01, 18.39,    &
479         &       19.55, 20.36, 21.01, 21.81, 21.76, 22.35, 22.68,    &
480         &       22.94, 22.93, 22.69, 22.94, 23.20, 23.31, 23.63,    &
481         &       23.47, 23.68, 23.62, 23.53, 23.59, 23.99, 23.80,    &
482         &       24.20, 24.99, 24.97, 25.72, 26.21, 26.37, 26.89,    &
483         &       27.68, 28.13, 28.94, 29.42, 29.66, 30.29, 30.96,    &
484         &       31.09, 31.59, 31.52, 31.92, 32.45, 32.91, 33.39,    &
485         &       33.80, 34.23, 34.73, 35.40, 36.14, 36.99, 37.87,    &
486         &       38.75, 39.70, 40.70, 41.68, 42.82, 43.96, 45.00,    &
487         &       45.98, 47.00, 48.03, 49.10, 50.10, 50.97, 51.81,    &
488         &       52.57 /
489      !!----------------------------------------------------------------------
490
491      cycle = 360.0
492      ilc = 0
493      icent = year / 100
494      yr = year - icent * 100
495      t = icent - 20
496!
497!     for the following equations
498!     time origin is fixed at 00 hr of jan 1st,2000.
499!     see notes by cartwright
500!
501!!gm  old coding style, use CASE instead  and avoid GOTO (obsolescence in fortran 90)
502!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
503      it = icent - 20
504      if (it) 1,2,2
505    1 iday = it/4 -it
506      go to 3
507    2 iday = (it+3)/4 - it
508!
509!     t is in julian century
510!     correction in gegorian calander where only century year divisible
511!     by 4 is leap year.
512!
513    3 continue
514!
515      td = 0.0
516!
517!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
518      if (yr) 4,5,4
519!
520    4 iyd = 365*yr
521      ild = (yr-1)/4
522      if((icent - (icent/4)*4) .eq. 0) ilc = 1
523      td = iyd + ild + ilc
524!
525    5 td = td + iday + vd -1.0 - 0.5
526      t = t + (td/36525.0)
527!
528      ipos=year-1899
529      if (ipos .lt. 0) go to 7
530      if (ipos .gt. 83) go to 6
531!
532      delta = (delt(ipos+1)+delt(ipos))/2.0
533      go to 7
534!
535    6 delta= (65.0-50.5)/20.0*(year-1980)+50.5
536!
537    7 deltat = delta * 1.0e-6
538!
539!!gm   precision of the computation   :  example for s it should be replace by:
540!!gm  s   = 218.3165 + (481267.8813 - 0.0016*t)*t + 152.0*deltat   ==>  more precise  modify the last digits results
541      s   = 218.3165 + 481267.8813*t - 0.0016*t*t + 152.0*deltat
542      h   = 280.4661 + 36000.7698 *t + 0.0003*t*t +  11.0*deltat
543      p   =  83.3535 + 4069.0139  *t - 0.0103*t*t +       deltat
544      en  = 234.9555 + 1934.1363  *t - 0.0021*t*t +       deltat
545      p1  = 282.9384 + 1.7195     *t + 0.0005*t*t
546      !
547      nn = s / cycle
548      s = s - nn * cycle
549      IF( s < 0.e0 )   s = s + cycle
550      !
551      nn = h / cycle
552      h = h - cycle * nn
553      IF( h < 0.e0 )   h = h + cycle
554      !
555      nn = p / cycle
556      p = p - cycle * nn
557      IF( p < 0.e0)   p = p + cycle
558      !
559      nn = en / cycle
560      en = en - cycle * nn
561      IF( en < 0.e0 )   en = en + cycle
562      en = cycle - en
563      !
564      nn = p1 / cycle
565      p1 = p1 - nn * cycle
566      !
567   END SUBROUTINE shpen
568
569
570   SUBROUTINE ufset( p, cn, b, a )
571      !!----------------------------------------------------------------------
572      !!                    ***  SUBROUTINE ufset  ***
573      !!                   
574      !! ** Purpose : - calculate nodal parameters for the tides
575      !!               
576      !!----------------------------------------------------------------------
577!!gm  doctor naming of dummy argument variables!!!   and all local variables
578!!gm  nc is jptides_max ????
579      integer nc
580      parameter (nc=15)
581      REAL(wp) p,cn
582      !!
583     
584!!gm  rad is already a public variable defined in phycst.F90 ....   ==> doctor norme  local real start with "z"
585      REAL(wp) ::   w1, w2, w3, w4, w5, w6, w7, w8, nw, pw, rad
586      REAL(wp) ::   a(nc), b(nc)
587      INTEGER  ::   ndc, k
588      !!----------------------------------------------------------------------
589
590      ndc = nc
591
592!    a=f       ,  b =u
593!    t is  zero as compared to tifa.
594!! use rad defined in phycst   (i.e.  add a USE phycst at the begining of the module
595      rad = 6.2831852d0/360.0
596      pw = p  * rad
597      nw = cn * rad
598      w1 = cos(   nw )
599      w2 = cos( 2*nw )
600      w3 = cos( 3*nw )
601      w4 = sin(   nw )
602      w5 = sin( 2*nw )
603      w6 = sin( 3*nw )
604      w7 = 1. - 0.2505 * COS( 2*pw      ) - 0.1102 * COS(2*pw-nw )   &
605         &    - 0.156  * COS( 2*pw-2*nw ) - 0.037  * COS(     nw )
606      w8 =    - 0.2505 * SIN( 2*pw      ) - 0.1102 * SIN(2*pw-nw )   &
607         &    - 0.0156 * SIN( 2*pw-2*nw ) - 0.037  * SIN(     nw )
608      !
609      a(1) = 1.0089 + 0.1871 * w1 - 0.0147 * w2 + 0.0014 * w3
610      b(1) =          0.1885 * w4 - 0.0234 * w5 + 0.0033 * w6
611      !   q1
612      a(2) = a(1)
613      b(2) = b(1)
614      !   o1
615      a(3) = 1.0
616      b(3) = 0.0
617      !   p1
618      a(4) = 1.0
619      b(4) = 0.0
620      !   s1
621      a(5) = 1.0060+0.1150*w1- 0.0088*w2 +0.0006*w3
622      b(5) = -0.1546*w4 + 0.0119*w5 -0.0012*w6
623      !   k1
624      a(6) =1.0004 -0.0373*w1+ 0.0002*w2
625      b(6) = -0.0374*w4
626      !  2n2
627      a(7) = a(6)
628      b(7) = b(6)
629      !  mu2
630      a(8) = a(6)
631      b(8) = b(6)
632      !   n2
633      a(9) = a(6)
634      b(9) = b(6)
635      !  nu2
636      a(10) = a(6)
637      b(10) = b(6)
638      !   m2
639      a(11) = SQRT( w7 * w7 + w8 * w8 )
640      b(11) = ATAN( w8 / w7 )
641!!gmuse rpi  instead of 3.141992 ???   true pi is rpi=3.141592653589793_wp  .....   ????
642      IF( w7 < 0.e0 )   b(11) = b(11) + 3.141992
643      !   l2
644      a(12) = 1.0
645      b(12) = 0.0
646      !   t2
647      a(13)= a(12)
648      b(13)= b(12)
649      !   s2
650      a(14) = 1.0241+0.2863*w1+0.0083*w2 -0.0015*w3
651      b(14) = -0.3096*w4 + 0.0119*w5 - 0.0007*w6
652      !   k2
653      a(15) = a(6)*a(6)
654      b(15) = 2*b(6)
655      !   m4
656!!gm  old coding,  remove GOTO and label of lines
657!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
658      DO 40 k = 1,ndc
659         b(k) = b(k)/rad
66032       if (b(k)) 34,35,35
66134       b(k) = b(k) + 360.0
662         go to 32
66335       if (b(k)-360.0) 40,37,37
66437       b(k) = b(k)-360.0
665         go to 35
66640    continue
667      !
668   END SUBROUTINE ufset
669   
670
671   SUBROUTINE vset( s,h,p,en,p1,v)
672      !!----------------------------------------------------------------------
673      !!                    ***  SUBROUTINE vset  ***
674      !!                   
675      !! ** Purpose : - calculate tidal phases for 0000gmt on start day of run
676      !!               
677      !!----------------------------------------------------------------------
678!!gm  doctor naming of dummy argument variables!!!   and all local variables
679!!gm  nc is jptides_max ????
680!!gm  en argument is not used: suppress it ?
681      integer nc
682      parameter (nc=15)
683      real(wp) s,h,p,en,p1
684      real(wp)   v(nc)
685      !!
686      integer ndc, k
687      !!----------------------------------------------------------------------
688
689      ndc = nc
690      !   v s  are computed here.
691      v(1) =-3*s +h +p +270      ! Q1
692      v(2) =-2*s +h +270.0       ! O1
693      v(3) =-h +270              ! P1
694      v(4) =180                  ! S1
695      v(5) =h +90.0              ! K1
696      v(6) =-4*s +2*h +2*p       ! 2N2
697      v(7) =-4*(s-h)             ! MU2
698      v(8) =-3*s +2*h +p         ! N2
699      v(9) =-3*s +4*h -p         ! MU2
700      v(10) =-2*s +2*h           ! M2
701      v(11) =-s +2*h -p +180     ! L2
702      v(12) =-h +p1              ! T2
703      v(13) =0                   ! S2
704      v(14) =h+h                 ! K2
705      v(15) =2*v(10)             ! M4
706!
707!!gm  old coding,  remove GOTO and label of lines
708!!gm obsol(   1): Arithmetic IF statement is used   ===>  remove this in Fortran 90
709      do 72 k = 1, ndc
71069       if( v(k) )   70,71,71
71170       v(k) = v(k)+360.0
712         go to 69
71371       if( v(k) - 360.0 )   72,73,73
71473       v(k) = v(k)-360.0
715         go to 71
71672    continue
717      !
718   END SUBROUTINE vset
719
720#else
721   !!----------------------------------------------------------------------
722   !!   Dummy module         NO Unstruct Open Boundary Conditions for tides
723   !!----------------------------------------------------------------------
724!!gm  are you sure we need to define filtide and tide_cpt ?
725   CHARACTER(len=80), PUBLIC               ::   filtide                !: Filename root for tidal input files
726   CHARACTER(len=4 ), PUBLIC, DIMENSION(1) ::   tide_cpt               !: Names of tidal components used.
727
728CONTAINS
729   SUBROUTINE tide_init                ! Empty routine
730   END SUBROUTINE tide_init
731   SUBROUTINE tide_data                ! Empty routine
732   END SUBROUTINE tide_data
733   SUBROUTINE tide_update( kt, kit )   ! Empty routine
734      WRITE(*,*) 'tide_update: You should not have seen this print! error?', kt, kit
735   END SUBROUTINE tide_update
736#endif
737
738   !!======================================================================
739END MODULE bdytides
Note: See TracBrowser for help on using the repository browser.