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 NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/BDY – NEMO

source: NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019/src/OCE/BDY/bdytides.F90 @ 12184

Last change on this file since 12184 was 12184, checked in by smueller, 4 years ago

Merging of the developments in /NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis@12122 with respect to /NEMO/trunk@12072 into /NEMO/branches/2019/dev_r12072_MERGE_OPTION2_2019 (tickets #2175 and #2194)
  • Property svn:keywords set to Id
File size: 22.7 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  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods
11   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
12   !!----------------------------------------------------------------------
13   !!   bdytide_init  : read of namelist and initialisation of tidal harmonics data
14   !!   tide_update   : calculation of tidal forcing at each timestep
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and tracers
17   USE dom_oce        ! ocean space and time domain
18   USE phycst         ! physical constants
19   USE bdy_oce        ! ocean open boundary conditions
20   USE tide_mod       !
21   USE daymod         ! calendar
22   !
23   USE in_out_manager ! I/O units
24   USE iom            ! xIO server
25   USE fldread        !
26   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   bdytide_init     ! routine called in bdy_init
32   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts
33
34   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
35      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0     !: Tidal constituents : SSH0   (read in file)
36      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0, v0   !: Tidal constituents : U0, V0 (read in file)
37      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh      !: Tidal constituents : SSH    (after nodal cor.)
38      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u , v    !: Tidal constituents : U , V  (after nodal cor.)
39   END TYPE TIDES_DATA
40
41!$AGRIF_DO_NOT_TREAT
42   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
43!$AGRIF_END_DO_NOT_TREAT
44   TYPE(OBC_DATA)  , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component)
45
46   INTEGER ::   kt_tide
47
48   !!----------------------------------------------------------------------
49   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
50   !! $Id$
51   !! Software governed by the CeCILL license (see ./LICENSE)
52   !!----------------------------------------------------------------------
53CONTAINS
54
55   SUBROUTINE bdytide_init
56      !!----------------------------------------------------------------------
57      !!                    ***  SUBROUTINE bdytide_init  ***
58      !!                     
59      !! ** Purpose : - Read in namelist for tides and initialise external
60      !!                tidal harmonics data
61      !!
62      !!----------------------------------------------------------------------
63      !! namelist variables
64      !!-------------------
65      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
66      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data
67      !!
68      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices
69      INTEGER                                   ::   ii, ij              !: dummy loop indices
70      INTEGER                                   ::   inum, igrd
71      INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays)
72      INTEGER                                   ::   ios                 ! Local integer output status for namelist read
73      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file
74      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data
75      REAL(wp),ALLOCATABLE, DIMENSION(:,:)      ::   ztr, zti            !:  "     "    "   "   "   "        "      "
76      !!
77      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
78      !!
79      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta
80      !!----------------------------------------------------------------------
81      !
82      IF(lwp) WRITE(numout,*)
83      IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
84      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
85
86      REWIND(numnam_cfg)
87
88      DO ib_bdy = 1, nb_bdy
89         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
90            !
91            td => tides(ib_bdy)
92
93            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
94            filtide(:) = ''
95
96            REWIND( numnam_ref )
97            READ  ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
98901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_tide in reference namelist' )
99            ! Don't REWIND here - may need to read more than one of these namelists.
100            READ  ( numnam_cfg, nambdy_tide, IOSTAT = ios, ERR = 902 )
101902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist' )
102            IF(lwm) WRITE ( numond, nambdy_tide )
103            !                                               ! Parameter control and print
104            IF(lwp) WRITE(numout,*) '  '
105            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
106            IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta
107            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
108            IF(lwp) THEN
109                    WRITE(numout,*) '             Tidal components: ' 
110               DO itide = 1, nb_harmo
111                  WRITE(numout,*)  '                 ', tide_harmonics(itide)%cname_tide 
112               END DO
113            ENDIF
114            IF(lwp) WRITE(numout,*) ' '
115
116            ! Allocate space for tidal harmonics data - get size from OBC data arrays
117            ! -----------------------------------------------------------------------
118
119            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
120            ! relaxation area     
121            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = idx_bdy(ib_bdy)%nblen   (:)
122            ELSE                                   ;   ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:)
123            ENDIF
124
125            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
126            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
127
128            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
129            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
130
131            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
132            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
133
134            td%ssh0(:,:,:) = 0._wp
135            td%ssh (:,:,:) = 0._wp
136            td%u0  (:,:,:) = 0._wp
137            td%u   (:,:,:) = 0._wp
138            td%v0  (:,:,:) = 0._wp
139            td%v   (:,:,:) = 0._wp
140
141            IF( ln_bdytide_2ddta ) THEN
142               ! It is assumed that each data file contains all complex harmonic amplitudes
143               ! given on the global domain (ie global, jpiglo x jpjglo)
144               !
145               ALLOCATE( zti(jpi,jpj), ztr(jpi,jpj) )
146               !
147               ! SSH fields
148               clfile = TRIM(filtide)//'_grid_T.nc'
149               CALL iom_open( clfile , inum ) 
150               igrd = 1                       ! Everything is at T-points here
151               DO itide = 1, nb_harmo
152                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) )
153                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) ) 
154                  DO ib = 1, ilen0(igrd)
155                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
156                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
157                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
158                     td%ssh0(ib,itide,1) = ztr(ii,ij)
159                     td%ssh0(ib,itide,2) = zti(ii,ij)
160                  END DO
161               END DO
162               CALL iom_close( inum )
163               !
164               ! U fields
165               clfile = TRIM(filtide)//'_grid_U.nc'
166               CALL iom_open( clfile , inum ) 
167               igrd = 2                       ! Everything is at U-points here
168               DO itide = 1, nb_harmo
169                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) )
170                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u2', zti(:,:) )
171                  DO ib = 1, ilen0(igrd)
172                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
173                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
174                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
175                     td%u0(ib,itide,1) = ztr(ii,ij)
176                     td%u0(ib,itide,2) = zti(ii,ij)
177                  END DO
178               END DO
179               CALL iom_close( inum )
180               !
181               ! V fields
182               clfile = TRIM(filtide)//'_grid_V.nc'
183               CALL iom_open( clfile , inum ) 
184               igrd = 3                       ! Everything is at V-points here
185               DO itide = 1, nb_harmo
186                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) )
187                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) )
188                  DO ib = 1, ilen0(igrd)
189                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
190                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
191                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
192                     td%v0(ib,itide,1) = ztr(ii,ij)
193                     td%v0(ib,itide,2) = zti(ii,ij)
194                  END DO
195               END DO 
196               CALL iom_close( inum )
197               !
198               DEALLOCATE( ztr, zti ) 
199               !
200            ELSE           
201               !
202               ! Read tidal data only on bdy segments
203               !
204               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
205               !
206               ! Open files and read in tidal forcing data
207               ! -----------------------------------------
208
209               DO itide = 1, nb_harmo
210                  !                                                              ! SSH fields
211                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc'
212                  CALL iom_open( clfile, inum )
213                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
214                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
215                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
216                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
217                  CALL iom_close( inum )
218                  !                                                              ! U fields
219                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc'
220                  CALL iom_open( clfile, inum )
221                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
222                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
223                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
224                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
225                  CALL iom_close( inum )
226                  !                                                              ! V fields
227                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc'
228                  CALL iom_open( clfile, inum )
229                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
230                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
231                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
232                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
233                  CALL iom_close( inum )
234                  !
235               END DO ! end loop on tidal components
236               !
237               DEALLOCATE( dta_read )
238               !
239            ENDIF ! ln_bdytide_2ddta=.true.
240            !
241            ! Allocate slow varying data in the case of time splitting:
242            ! Do it anyway because at this stage knowledge of free surface scheme is unknown
243            ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
244            ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
245            ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
246            dta_bdy_s(ib_bdy)%ssh(:) = 0._wp
247            dta_bdy_s(ib_bdy)%u2d(:) = 0._wp
248            dta_bdy_s(ib_bdy)%v2d(:) = 0._wp
249            !
250         ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2
251         !
252      END DO ! loop on ib_bdy
253      !
254   END SUBROUTINE bdytide_init
255
256
257   SUBROUTINE bdy_dta_tides( kt, kit, kt_offset )
258      !!----------------------------------------------------------------------
259      !!                 ***  SUBROUTINE bdy_dta_tides  ***
260      !!               
261      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
262      !!               
263      !!----------------------------------------------------------------------
264      INTEGER,           INTENT(in) ::   kt          ! Main timestep counter
265      INTEGER, OPTIONAL, INTENT(in) ::   kit         ! Barotropic timestep counter (for timesplitting option)
266      INTEGER, OPTIONAL, INTENT(in) ::   kt_offset   ! time offset in units of timesteps. NB. if kit
267      !                                              ! is present then units = subcycle timesteps.
268      !                                              ! kt_offset = 0  => get data at "now"    time level
269      !                                              ! kt_offset = -1 => get data at "before" time level
270      !                                              ! kt_offset = +1 => get data at "after"  time level
271      !                                              ! etc.
272      !
273      LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step
274      INTEGER  ::   itide, ib_bdy, ib, igrd   ! loop indices
275      INTEGER  ::   time_add                  ! time offset in units of timesteps
276      INTEGER, DIMENSION(jpbgrd)   ::   ilen0 
277      INTEGER, DIMENSION(1:jpbgrd) ::   nblen, nblenrim  ! short cuts
278      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist     
279      !!----------------------------------------------------------------------
280      !
281      lk_first_btstp=.TRUE.
282      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
283
284      time_add = 0
285      IF( PRESENT(kt_offset) ) THEN
286         time_add = kt_offset
287      ENDIF
288     
289      ! Absolute time from model initialization:   
290      IF( PRESENT(kit) ) THEN 
291         z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt
292      ELSE                             
293         z_arg = ( kt + time_add ) * rdt
294      ENDIF
295
296      ! Linear ramp on tidal component at open boundaries
297      zramp = 1.
298      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rn_tide_ramp_dt*rday),0.),1.)
299
300      DO ib_bdy = 1,nb_bdy
301         !
302         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
303            !
304            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
305            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
306            !
307            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:)
308            ELSE                                   ;   ilen0(:) = nblenrim(:)
309            ENDIF     
310            !
311            ! We refresh nodal factors every day below
312            ! This should be done somewhere else
313            IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
314               !
315               kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
316               !
317               IF(lwp) THEN
318               WRITE(numout,*)
319               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
320               WRITE(numout,*) '~~~~~~~~~~~~~~ '
321               ENDIF
322               !
323               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
324               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
325               !
326            ENDIF
327            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
328            !
329            ! If time splitting, initialize arrays from slow varying open boundary data:
330            IF ( PRESENT(kit) ) THEN           
331               IF ( dta_bdy(ib_bdy)%lneed_ssh   ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
332               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
333               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
334            ENDIF
335            !
336            ! Update open boundary data arrays:
337            DO itide = 1, nb_harmo
338               !
339               z_sarg = (z_arg + zoff) * tide_harmonics(itide)%omega
340               z_cost = zramp * COS( z_sarg )
341               z_sist = zramp * SIN( z_sarg )
342               !
343               IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN
344                  igrd=1                              ! SSH on tracer grid
345                  DO ib = 1, ilen0(igrd)
346                     dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
347                        &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
348                        &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
349                  END DO
350               ENDIF
351               !
352               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
353                  igrd=2                              ! U grid
354                  DO ib = 1, ilen0(igrd)
355                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
356                        &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
357                        &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
358                  END DO
359                  igrd=3                              ! V grid
360                  DO ib = 1, ilen0(igrd) 
361                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
362                        &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
363                        &                        tides(ib_bdy)%v(ib,itide,2)*z_sist )
364                  END DO
365               ENDIF
366            END DO             
367         END IF
368      END DO
369      !
370   END SUBROUTINE bdy_dta_tides
371
372
373   SUBROUTINE tide_init_elevation( idx, td )
374      !!----------------------------------------------------------------------
375      !!                 ***  ROUTINE tide_init_elevation  ***
376      !!----------------------------------------------------------------------
377      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
378      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
379      !
380      INTEGER ::   itide, igrd, ib       ! dummy loop indices
381      INTEGER, DIMENSION(1) ::   ilen0   ! length of boundary data (from OBC arrays)
382      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
383      !!----------------------------------------------------------------------
384      !
385      igrd=1   
386                              ! SSH on tracer grid.
387      ilen0(1) =  SIZE(td%ssh0(:,1,1))
388      !
389      ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) )
390      !
391      DO itide = 1, nb_harmo
392         DO ib = 1, ilen0(igrd)
393            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
394            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
395         END DO
396         DO ib = 1 , ilen0(igrd)
397            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
398            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u
399         ENDDO
400         DO ib = 1 , ilen0(igrd)
401            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
402            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
403         ENDDO
404      END DO
405      !
406      DEALLOCATE( mod_tide, phi_tide )
407      !
408   END SUBROUTINE tide_init_elevation
409
410
411   SUBROUTINE tide_init_velocities( idx, td )
412      !!----------------------------------------------------------------------
413      !!                 ***  ROUTINE tide_init_elevation  ***
414      !!----------------------------------------------------------------------
415      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
416      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
417      !
418      INTEGER ::   itide, igrd, ib       ! dummy loop indices
419      INTEGER, DIMENSION(3) ::   ilen0   ! length of boundary data (from OBC arrays)
420      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
421      !!----------------------------------------------------------------------
422      !
423      ilen0(2) =  SIZE(td%u0(:,1,1))
424      ilen0(3) =  SIZE(td%v0(:,1,1))
425      !
426      igrd=2                                 ! U grid.
427      !
428      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
429      !
430      DO itide = 1, nb_harmo
431         DO ib = 1, ilen0(igrd)
432            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
433            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
434         END DO
435         DO ib = 1, ilen0(igrd)
436            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
437            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
438         ENDDO
439         DO ib = 1, ilen0(igrd)
440            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
441            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
442         ENDDO
443      END DO
444      !
445      DEALLOCATE( mod_tide , phi_tide )
446      !
447      igrd=3                                 ! V grid.
448      !
449      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
450
451      DO itide = 1, nb_harmo
452         DO ib = 1, ilen0(igrd)
453            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
454            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
455         END DO
456         DO ib = 1, ilen0(igrd)
457            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
458            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
459         ENDDO
460         DO ib = 1, ilen0(igrd)
461            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
462            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
463         ENDDO
464      END DO
465      !
466      DEALLOCATE( mod_tide, phi_tide )
467      !
468  END SUBROUTINE tide_init_velocities
469
470   !!======================================================================
471END MODULE bdytides
472
Note: See TracBrowser for help on using the repository browser.