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

source: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 3294

Last change on this file since 3294 was 3294, checked in by rblod, 11 years ago

Merge of 3.4beta into the trunk

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