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/releases/r4.0/r4.0-HEAD/src/OCE/BDY – NEMO

source: NEMO/releases/r4.0/r4.0-HEAD/src/OCE/BDY/bdytides.F90 @ 12910

Last change on this file since 12910 was 12910, checked in by smasson, 4 years ago

r4.0-HEAD: bugfix potential out-of-bounds in bdydta/bdytides, see #2461
  • Property svn:keywords set to Id
File size: 27.4 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 tideini        !
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   bdytide_update   ! routine called in bdy_dta
33   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts
34
35   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
36      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0     !: Tidal constituents : SSH0   (read in file)
37      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0, v0   !: Tidal constituents : U0, V0 (read in file)
38      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh      !: Tidal constituents : SSH    (after nodal cor.)
39      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u , v    !: Tidal constituents : U , V  (after nodal cor.)
40   END TYPE TIDES_DATA
41
42!$AGRIF_DO_NOT_TREAT
43   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
44!$AGRIF_END_DO_NOT_TREAT
45   TYPE(OBC_DATA)  , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component)
46
47   !!----------------------------------------------------------------------
48   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
49   !! $Id$
50   !! Software governed by the CeCILL license (see ./LICENSE)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE bdytide_init
55      !!----------------------------------------------------------------------
56      !!                    ***  SUBROUTINE bdytide_init  ***
57      !!                     
58      !! ** Purpose : - Read in namelist for tides and initialise external
59      !!                tidal harmonics data
60      !!
61      !!----------------------------------------------------------------------
62      !! namelist variables
63      !!-------------------
64      CHARACTER(len=80)                         ::   filtide             ! Filename root for tidal input files
65      LOGICAL                                   ::   ln_bdytide_2ddta    ! If true, read 2d harmonic data
66      LOGICAL                                   ::   ln_bdytide_conj     ! If true, assume complex conjugate tidal data
67      !!
68      INTEGER                                   ::   ib_bdy, itide, ib   ! dummy loop indices
69      INTEGER                                   ::   ii, ij              ! dummy loop indices
70      INTEGER                                   ::   inum, igrd
71      INTEGER                                   ::   isz                 ! bdy data size
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      TYPE(  OBC_DATA), POINTER                 ::   dta                 ! local short cut
79      !!
80      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
81      !!----------------------------------------------------------------------
82      !
83      IF(lwp) WRITE(numout,*)
84      IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
85      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
86
87      REWIND(numnam_cfg)
88
89      DO ib_bdy = 1, nb_bdy
90         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
91            !
92            td  => tides(ib_bdy)
93            dta => dta_bdy(ib_bdy)
94         
95            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
96            filtide(:) = ''
97
98            REWIND( numnam_ref )
99            READ  ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
100901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_tide in reference namelist' )
101            ! Don't REWIND here - may need to read more than one of these namelists.
102            READ  ( numnam_cfg, nambdy_tide, IOSTAT = ios, ERR = 902 )
103902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist' )
104            IF(lwm) WRITE ( numond, nambdy_tide )
105            !                                               ! Parameter control and print
106            IF(lwp) WRITE(numout,*) '  '
107            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
108            IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta
109            IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj
110            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
111            IF(lwp) THEN
112                    WRITE(numout,*) '             Tidal components: ' 
113               DO itide = 1, nb_harmo
114                  WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide 
115               END DO
116            ENDIF
117            IF(lwp) WRITE(numout,*) ' '
118
119            ! Allocate space for tidal harmonics data - get size from BDY data arrays
120            ! Allocate also slow varying data in the case of time splitting:
121            ! Do it anyway because at this stage knowledge of free surface scheme is unknown
122            ! -----------------------------------------------------------------------
123            IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
124               isz = SIZE(dta%ssh)
125               ALLOCATE( td%ssh0( isz, nb_harmo, 2 ), td%ssh( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%ssh( isz ) )
126               dta_bdy_s(ib_bdy)%ssh(:) = 0._wp   ! needed?
127            ENDIF
128            IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
129               isz = SIZE(dta%u2d)
130               ALLOCATE( td%u0  ( isz, nb_harmo, 2 ), td%u  ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%u2d( isz ) )
131               dta_bdy_s(ib_bdy)%u2d(:) = 0._wp   ! needed?
132            ENDIF
133            IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
134               isz = SIZE(dta%v2d)
135               ALLOCATE( td%v0  ( isz, nb_harmo, 2 ), td%v  ( isz, nb_harmo, 2 ), dta_bdy_s(ib_bdy)%v2d( isz ) )
136               dta_bdy_s(ib_bdy)%v2d(:) = 0._wp   ! needed?
137            ENDIF
138
139            ! fill td%ssh0, td%u0, td%v0
140            ! -----------------------------------------------------------------------
141            IF( ln_bdytide_2ddta ) THEN
142               !
143               ! It is assumed that each data file contains all complex harmonic amplitudes
144               ! given on the global domain (ie global, jpiglo x jpjglo)
145               !
146               ALLOCATE( zti(jpi,jpj), ztr(jpi,jpj) )
147               !
148               ! SSH fields
149               IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
150                  clfile = TRIM(filtide)//'_grid_T.nc'
151                  CALL iom_open( clfile , inum ) 
152                  igrd = 1                       ! Everything is at T-points here
153                  DO itide = 1, nb_harmo
154                     CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
155                     CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 
156                     DO ib = 1, SIZE(dta%ssh)
157                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
158                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
159                        td%ssh0(ib,itide,1) = ztr(ii,ij)
160                        td%ssh0(ib,itide,2) = zti(ii,ij)
161                     END DO
162                  END DO
163                  CALL iom_close( inum )
164               END IF
165               !
166               ! U fields
167               IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
168                  clfile = TRIM(filtide)//'_grid_U.nc'
169                  CALL iom_open( clfile , inum ) 
170                  igrd = 2                       ! Everything is at U-points here
171                  DO itide = 1, nb_harmo
172                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
173                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
174                     DO ib = 1, SIZE(dta%u2d)
175                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
176                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
177                        td%u0(ib,itide,1) = ztr(ii,ij)
178                        td%u0(ib,itide,2) = zti(ii,ij)
179                     END DO
180                  END DO
181                  CALL iom_close( inum )
182               END IF
183               !
184               ! V fields
185               IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
186                  clfile = TRIM(filtide)//'_grid_V.nc'
187                  CALL iom_open( clfile , inum ) 
188                  igrd = 3                       ! Everything is at V-points here
189                  DO itide = 1, nb_harmo
190                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
191                     CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
192                     DO ib = 1, SIZE(dta%v2d)
193                        ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
194                        ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
195                        td%v0(ib,itide,1) = ztr(ii,ij)
196                        td%v0(ib,itide,2) = zti(ii,ij)
197                     END DO
198                  END DO
199                  CALL iom_close( inum )
200               END IF
201               !
202               DEALLOCATE( ztr, zti ) 
203               !
204            ELSE           
205               !
206               ! Read tidal data only on bdy segments
207               !
208               ALLOCATE( dta_read( MAXVAL( idx_bdy(ib_bdy)%nblen(:) ), 1, 1 ) )
209               !
210               ! Open files and read in tidal forcing data
211               ! -----------------------------------------
212
213               DO itide = 1, nb_harmo
214                  !                                                              ! SSH fields
215                  IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
216                     isz = SIZE(dta%ssh)
217                     clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
218                     CALL iom_open( clfile, inum )
219                     CALL fld_map( inum, 'z1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
220                     td%ssh0(:,itide,1) = dta_read(1:isz,1,1)
221                     CALL fld_map( inum, 'z2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
222                     td%ssh0(:,itide,2) = dta_read(1:isz,1,1)
223                     CALL iom_close( inum )
224                  ENDIF
225                  !                                                              ! U fields
226                  IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
227                     isz = SIZE(dta%u2d)
228                     clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
229                     CALL iom_open( clfile, inum )
230                     CALL fld_map( inum, 'u1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
231                     td%u0(:,itide,1) = dta_read(1:isz,1,1)
232                     CALL fld_map( inum, 'u2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
233                     td%u0(:,itide,2) = dta_read(1:isz,1,1)
234                     CALL iom_close( inum )
235                  ENDIF
236                  !                                                              ! V fields
237                  IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
238                     isz = SIZE(dta%v2d)
239                     clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
240                     CALL iom_open( clfile, inum )
241                     CALL fld_map( inum, 'v1', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
242                     td%v0(:,itide,1) = dta_read(1:isz,1,1)
243                     CALL fld_map( inum, 'v2', dta_read(1:isz,1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
244                     td%v0(:,itide,2) = dta_read(1:isz,1,1)
245                     CALL iom_close( inum )
246                  ENDIF
247                  !
248               END DO ! end loop on tidal components
249               !
250               DEALLOCATE( dta_read )
251               !
252            ENDIF ! ln_bdytide_2ddta=.true.
253            !
254            IF( ln_bdytide_conj ) THEN    ! assume complex conjugate in data files
255               IF( ASSOCIATED(dta%ssh) )   td%ssh0(:,:,2) = - td%ssh0(:,:,2)
256               IF( ASSOCIATED(dta%u2d) )   td%u0  (:,:,2) = - td%u0  (:,:,2)
257               IF( ASSOCIATED(dta%v2d) )   td%v0  (:,:,2) = - td%v0  (:,:,2)
258            ENDIF
259            !
260         ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2
261         !
262      END DO ! loop on ib_bdy
263      !
264   END SUBROUTINE bdytide_init
265
266
267   SUBROUTINE bdytide_update( kt, idx, dta, td, kit, kt_offset )
268      !!----------------------------------------------------------------------
269      !!                 ***  SUBROUTINE bdytide_update  ***
270      !!               
271      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
272      !!               
273      !!----------------------------------------------------------------------
274      INTEGER          , INTENT(in   ) ::   kt          ! Main timestep counter
275      TYPE(OBC_INDEX)  , INTENT(in   ) ::   idx         ! OBC indices
276      TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data
277      TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data
278      INTEGER, OPTIONAL, INTENT(in   ) ::   kit         ! Barotropic timestep counter (for timesplitting option)
279      INTEGER, OPTIONAL, INTENT(in   ) ::   kt_offset   ! time offset in units of timesteps. NB. if kit
280      !                                                 ! is present then units = subcycle timesteps.
281      !                                                 ! kt_offset = 0  => get data at "now"    time level
282      !                                                 ! kt_offset = -1 => get data at "before" time level
283      !                                                 ! kt_offset = +1 => get data at "after"  time level
284      !                                                 ! etc.
285      !
286      INTEGER  ::   itide, ib             ! dummy loop indices
287      INTEGER  ::   time_add              ! time offset in units of timesteps
288      INTEGER  ::   isz                   ! bdy data size
289      REAL(wp) ::   z_arg, z_sarg, zflag, zramp   ! local scalars   
290      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
291      !!----------------------------------------------------------------------
292      !
293      zflag=1
294      IF ( PRESENT(kit) ) THEN
295        IF ( kit /= 1 ) zflag=0
296      ENDIF
297      !
298      IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN
299        !
300        kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
301        !
302        IF(lwp) THEN
303           WRITE(numout,*)
304           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
305           WRITE(numout,*) '~~~~~~~~~~~~~~ '
306        ENDIF
307        !
308        CALL tide_init_elevation ( idx, td )
309        CALL tide_init_velocities( idx, td )
310        !
311      ENDIF
312
313      time_add = 0
314      IF( PRESENT(kt_offset) ) THEN
315         time_add = kt_offset
316      ENDIF
317         
318      IF( PRESENT(kit) ) THEN 
319         z_arg = ((kt-kt_tide) * rdt + (kit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) )
320      ELSE                             
321         z_arg = ((kt-kt_tide)+time_add) * rdt
322      ENDIF
323
324      ! Linear ramp on tidal component at open boundaries
325      zramp = 1._wp
326      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp)
327
328      DO itide = 1, nb_harmo
329         z_sarg = z_arg * omega_tide(itide)
330         z_cost(itide) = COS( z_sarg )
331         z_sist(itide) = SIN( z_sarg )
332      END DO
333
334      DO itide = 1, nb_harmo
335         ! SSH on tracer grid
336         IF( ASSOCIATED(dta%ssh) ) THEN   ! we use bdy ssh on this mpi subdomain
337           DO ib = 1, SIZE(dta%ssh)
338               dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide))
339            END DO
340         ENDIF
341         ! U grid
342         IF( ASSOCIATED(dta%u2d) ) THEN   ! we use bdy u2d on this mpi subdomain
343            DO ib = 1, SIZE(dta%u2d)
344               dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide))
345            END DO
346         ENDIF
347         ! V grid
348         IF( ASSOCIATED(dta%v2d) ) THEN   ! we use bdy v2d on this mpi subdomain
349            DO ib = 1, SIZE(dta%v2d) 
350               dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide))
351            END DO
352         ENDIF
353      END DO
354      !
355   END SUBROUTINE bdytide_update
356
357
358   SUBROUTINE bdy_dta_tides( kt, kit, kt_offset )
359      !!----------------------------------------------------------------------
360      !!                 ***  SUBROUTINE bdy_dta_tides  ***
361      !!               
362      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
363      !!               
364      !!----------------------------------------------------------------------
365      INTEGER,           INTENT(in) ::   kt          ! Main timestep counter
366      INTEGER, OPTIONAL, INTENT(in) ::   kit         ! Barotropic timestep counter (for timesplitting option)
367      INTEGER, OPTIONAL, INTENT(in) ::   kt_offset   ! time offset in units of timesteps. NB. if kit
368      !                                              ! is present then units = subcycle timesteps.
369      !                                              ! kt_offset = 0  => get data at "now"    time level
370      !                                              ! kt_offset = -1 => get data at "before" time level
371      !                                              ! kt_offset = +1 => get data at "after"  time level
372      !                                              ! etc.
373      !
374      LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step
375      INTEGER  ::   itide, ib_bdy, ib         ! loop indices
376      INTEGER  ::   time_add                  ! time offset in units of timesteps
377      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist     
378      !!----------------------------------------------------------------------
379      !
380      lk_first_btstp=.TRUE.
381      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
382
383      time_add = 0
384      IF( PRESENT(kt_offset) ) THEN
385         time_add = kt_offset
386      ENDIF
387     
388      ! Absolute time from model initialization:   
389      IF( PRESENT(kit) ) THEN 
390         z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt
391      ELSE                             
392         z_arg = ( kt + time_add ) * rdt
393      ENDIF
394
395      ! Linear ramp on tidal component at open boundaries
396      zramp = 1.
397      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.)
398
399      DO ib_bdy = 1,nb_bdy
400         !
401         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
402            !
403            ! We refresh nodal factors every day below
404            ! This should be done somewhere else
405            IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
406               !
407               kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
408               !
409               IF(lwp) THEN
410               WRITE(numout,*)
411               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
412               WRITE(numout,*) '~~~~~~~~~~~~~~ '
413               ENDIF
414               !
415               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
416               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
417               !
418            ENDIF
419            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
420            !
421            ! If time splitting, initialize arrays from slow varying open boundary data:
422            IF ( PRESENT(kit) ) THEN           
423               IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) dta_bdy(ib_bdy)%ssh(:) = dta_bdy_s(ib_bdy)%ssh(:)
424               IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) dta_bdy(ib_bdy)%u2d(:) = dta_bdy_s(ib_bdy)%u2d(:)
425               IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) dta_bdy(ib_bdy)%v2d(:) = dta_bdy_s(ib_bdy)%v2d(:)
426            ENDIF
427            !
428            ! Update open boundary data arrays:
429            DO itide = 1, nb_harmo
430               !
431               z_sarg = (z_arg + zoff) * omega_tide(itide)
432               z_cost = zramp * COS( z_sarg )
433               z_sist = zramp * SIN( z_sarg )
434               !
435               IF ( ASSOCIATED(dta_bdy(ib_bdy)%ssh) ) THEN   ! SSH on tracer grid
436                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%ssh)
437                     dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
438                        &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
439                        &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
440                  END DO
441               ENDIF
442               !
443               IF ( ASSOCIATED(dta_bdy(ib_bdy)%u2d) ) THEN  ! U grid
444                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%u2d)
445                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
446                        &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
447                        &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
448                  END DO
449               ENDIF
450               !
451               IF ( ASSOCIATED(dta_bdy(ib_bdy)%v2d) ) THEN   ! V grid
452                  DO ib = 1, SIZE(dta_bdy(ib_bdy)%v2d)
453                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
454                        &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
455                        &                        tides(ib_bdy)%v(ib,itide,2)*z_sist )
456                  END DO
457               ENDIF
458               !
459            END DO             
460         END IF
461      END DO
462      !
463   END SUBROUTINE bdy_dta_tides
464
465
466   SUBROUTINE tide_init_elevation( idx, td )
467      !!----------------------------------------------------------------------
468      !!                 ***  ROUTINE tide_init_elevation  ***
469      !!----------------------------------------------------------------------
470      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
471      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
472      !
473      INTEGER ::   itide, isz, ib       ! dummy loop indices
474      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
475      !!----------------------------------------------------------------------
476      !
477      IF( ASSOCIATED(td%ssh0) ) THEN   ! SSH on tracer grid.
478         !
479         isz = SIZE( td%ssh0, dim = 1 )
480         ALLOCATE( mod_tide(isz), phi_tide(isz) )
481         !
482         DO itide = 1, nb_harmo
483            DO ib = 1, isz
484               mod_tide(ib)=SQRT( td%ssh0(ib,itide,1)*td%ssh0(ib,itide,1) + td%ssh0(ib,itide,2)*td%ssh0(ib,itide,2) )
485               phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
486            END DO
487            DO ib = 1, isz
488               mod_tide(ib)=mod_tide(ib)*ftide(itide)
489               phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
490            END DO
491            DO ib = 1, isz
492               td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
493               td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
494            END DO
495         END DO
496         !
497         DEALLOCATE( mod_tide, phi_tide )
498         !
499      ENDIF
500      !
501   END SUBROUTINE tide_init_elevation
502
503
504   SUBROUTINE tide_init_velocities( idx, td )
505      !!----------------------------------------------------------------------
506      !!                 ***  ROUTINE tide_init_elevation  ***
507      !!----------------------------------------------------------------------
508      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
509      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
510      !
511      INTEGER ::   itide, isz, ib        ! dummy loop indices
512      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
513      !!----------------------------------------------------------------------
514      !
515      IF( ASSOCIATED(td%u0) ) THEN   ! U grid. we use bdy u2d on this mpi subdomain
516         !
517         isz = SIZE( td%u0, dim = 1 )
518         ALLOCATE( mod_tide(isz), phi_tide(isz) )
519         !
520         DO itide = 1, nb_harmo
521            DO ib = 1, isz
522               mod_tide(ib)=SQRT( td%u0(ib,itide,1)*td%u0(ib,itide,1) + td%u0(ib,itide,2)*td%u0(ib,itide,2) )
523               phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
524            END DO
525            DO ib = 1, isz
526               mod_tide(ib)=mod_tide(ib)*ftide(itide)
527               phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
528            END DO
529            DO ib = 1, isz
530               td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
531               td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
532            END DO
533         END DO
534         !
535         DEALLOCATE( mod_tide, phi_tide )
536         !
537      ENDIF
538      !
539      IF( ASSOCIATED(td%v0) ) THEN   ! V grid. we use bdy u2d on this mpi subdomain
540         !
541         isz = SIZE( td%v0, dim = 1 )
542         ALLOCATE( mod_tide(isz), phi_tide(isz) )
543         !
544         DO itide = 1, nb_harmo
545            DO ib = 1, isz
546               mod_tide(ib)=SQRT( td%v0(ib,itide,1)*td%v0(ib,itide,1) + td%v0(ib,itide,2)*td%v0(ib,itide,2) )
547               phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
548            END DO
549            DO ib = 1, isz
550               mod_tide(ib)=mod_tide(ib)*ftide(itide)
551               phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
552            END DO
553            DO ib = 1, isz
554               td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
555               td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
556            END DO
557         END DO
558         !
559         DEALLOCATE( mod_tide, phi_tide )
560         !
561      ENDIF
562      !
563   END SUBROUTINE tide_init_velocities
564
565   !!======================================================================
566END MODULE bdytides
567
Note: See TracBrowser for help on using the repository browser.