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_r12114_ticket_2263/src/OCE/BDY – NEMO

source: NEMO/branches/2019/dev_r12114_ticket_2263/src/OCE/BDY/bdytides.F90 @ 12127

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

dev_r12114_ticket_2263: replace integer kt_offset by real pt_offset, see #2263
  • Property svn:keywords set to Id
File size: 26.3 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, 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, ln_bdytide_conj
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,*) '             assume complex conjugate   : ', ln_bdytide_conj
108            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
109            IF(lwp) THEN
110                    WRITE(numout,*) '             Tidal components: ' 
111               DO itide = 1, nb_harmo
112                  WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide 
113               END DO
114            ENDIF
115            IF(lwp) WRITE(numout,*) ' '
116
117            ! Allocate space for tidal harmonics data - get size from OBC data arrays
118            ! -----------------------------------------------------------------------
119
120            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
121            ! relaxation area     
122            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = idx_bdy(ib_bdy)%nblen   (:)
123            ELSE                                   ;   ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:)
124            ENDIF
125
126            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
127            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
128
129            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
130            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
131
132            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
133            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
134
135            td%ssh0(:,:,:) = 0._wp
136            td%ssh (:,:,:) = 0._wp
137            td%u0  (:,:,:) = 0._wp
138            td%u   (:,:,:) = 0._wp
139            td%v0  (:,:,:) = 0._wp
140            td%v   (:,:,:) = 0._wp
141
142            IF( ln_bdytide_2ddta ) THEN
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               clfile = TRIM(filtide)//'_grid_T.nc'
150               CALL iom_open( clfile , inum ) 
151               igrd = 1                       ! Everything is at T-points here
152               DO itide = 1, nb_harmo
153                  CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
154                  CALL iom_get( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 
155                  DO ib = 1, ilen0(igrd)
156                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
157                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
158                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
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               !
165               ! U fields
166               clfile = TRIM(filtide)//'_grid_U.nc'
167               CALL iom_open( clfile , inum ) 
168               igrd = 2                       ! Everything is at U-points here
169               DO itide = 1, nb_harmo
170                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
171                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
172                  DO ib = 1, ilen0(igrd)
173                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
174                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
175                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
176                     td%u0(ib,itide,1) = ztr(ii,ij)
177                     td%u0(ib,itide,2) = zti(ii,ij)
178                  END DO
179               END DO
180               CALL iom_close( inum )
181               !
182               ! V fields
183               clfile = TRIM(filtide)//'_grid_V.nc'
184               CALL iom_open( clfile , inum ) 
185               igrd = 3                       ! Everything is at V-points here
186               DO itide = 1, nb_harmo
187                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
188                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
189                  DO ib = 1, ilen0(igrd)
190                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
191                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
192                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
193                     td%v0(ib,itide,1) = ztr(ii,ij)
194                     td%v0(ib,itide,2) = zti(ii,ij)
195                  END DO
196               END DO 
197               CALL iom_close( inum )
198               !
199               DEALLOCATE( ztr, zti ) 
200               !
201            ELSE           
202               !
203               ! Read tidal data only on bdy segments
204               !
205               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
206               !
207               ! Open files and read in tidal forcing data
208               ! -----------------------------------------
209
210               DO itide = 1, nb_harmo
211                  !                                                              ! SSH fields
212                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
213                  CALL iom_open( clfile, inum )
214                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
215                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
216                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
217                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
218                  CALL iom_close( inum )
219                  !                                                              ! U fields
220                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
221                  CALL iom_open( clfile, inum )
222                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
223                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
224                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
225                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
226                  CALL iom_close( inum )
227                  !                                                              ! V fields
228                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
229                  CALL iom_open( clfile, inum )
230                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
231                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
232                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
233                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
234                  CALL iom_close( inum )
235                  !
236               END DO ! end loop on tidal components
237               !
238               DEALLOCATE( dta_read )
239               !
240            ENDIF ! ln_bdytide_2ddta=.true.
241            !
242            IF( ln_bdytide_conj ) THEN    ! assume complex conjugate in data files
243               td%ssh0(:,:,2) = - td%ssh0(:,:,2)
244               td%u0  (:,:,2) = - td%u0  (:,:,2)
245               td%v0  (:,:,2) = - td%v0  (:,:,2)
246            ENDIF
247            !
248            ! Allocate slow varying data in the case of time splitting:
249            ! Do it anyway because at this stage knowledge of free surface scheme is unknown
250            ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
251            ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
252            ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
253            dta_bdy_s(ib_bdy)%ssh(:) = 0._wp
254            dta_bdy_s(ib_bdy)%u2d(:) = 0._wp
255            dta_bdy_s(ib_bdy)%v2d(:) = 0._wp
256            !
257         ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2
258         !
259      END DO ! loop on ib_bdy
260      !
261   END SUBROUTINE bdytide_init
262
263
264   SUBROUTINE bdytide_update( kt, idx, dta, td, kit, pt_offset )
265      !!----------------------------------------------------------------------
266      !!                 ***  SUBROUTINE bdytide_update  ***
267      !!               
268      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
269      !!               
270      !!----------------------------------------------------------------------
271      INTEGER          , INTENT(in   ) ::   kt          ! Main timestep counter
272      TYPE(OBC_INDEX)  , INTENT(in   ) ::   idx         ! OBC indices
273      TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data
274      TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data
275      INTEGER, OPTIONAL, INTENT(in   ) ::   kit         ! Barotropic timestep counter (for timesplitting option)
276      REAL(wp), OPTIONAL,INTENT(in   ) ::   pt_offset   ! time offset in units of timesteps. NB. if kit
277      !                                                 ! is present then units = subcycle timesteps.
278      !
279      INTEGER  ::   itide, igrd, ib       ! dummy loop indices
280      INTEGER, DIMENSION(3) ::   ilen0    ! length of boundary data (from OBC arrays)
281      REAL(wp) ::   z_arg, z_sarg, zflag, zramp, zt_offset  ! local scalars   
282      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
283      !!----------------------------------------------------------------------
284      !
285      ilen0(1) =  SIZE(td%ssh(:,1,1))
286      ilen0(2) =  SIZE(td%u(:,1,1))
287      ilen0(3) =  SIZE(td%v(:,1,1))
288
289      zflag=1
290      IF ( PRESENT(kit) ) THEN
291        IF ( kit /= 1 ) zflag=0
292      ENDIF
293
294      IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN
295        !
296        kt_tide = kt - NINT((REAL(nsec_day,wp) - 0.5_wp * rdt)/rdt)
297        !
298        IF(lwp) THEN
299           WRITE(numout,*)
300           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
301           WRITE(numout,*) '~~~~~~~~~~~~~~ '
302        ENDIF
303        !
304        CALL tide_init_elevation ( idx, td )
305        CALL tide_init_velocities( idx, td )
306        !
307      ENDIF
308
309      zt_offset = 0._wp
310      IF( PRESENT(pt_offset) )   zt_offset = pt_offset
311         
312      IF( PRESENT(kit) ) THEN 
313         z_arg =   REAL( kt - kt_tide, wp) * rdt + ( REAL(kit, wp) + 0.5_wp * ( zt_offset - 1._wp ) ) * rdt / REAL(nn_baro, wp)
314      ELSE                             
315         z_arg = ( REAL( kt - kt_tide, wp) + zt_offset ) * rdt
316      ENDIF
317
318      ! Linear ramp on tidal component at open boundaries
319      zramp = 1._wp
320      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + REAL(kt_tide-nit000,wp)*rdt)/(rdttideramp*rday),0._wp),1._wp)
321
322      DO itide = 1, nb_harmo
323         z_sarg = z_arg * omega_tide(itide)
324         z_cost(itide) = COS( z_sarg )
325         z_sist(itide) = SIN( z_sarg )
326      END DO
327
328      DO itide = 1, nb_harmo
329         igrd=1                              ! SSH on tracer grid
330         DO ib = 1, ilen0(igrd)
331            dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide))
332         END DO
333         igrd=2                              ! U grid
334         DO ib = 1, ilen0(igrd)
335            dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide))
336         END DO
337         igrd=3                              ! V grid
338         DO ib = 1, ilen0(igrd) 
339            dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide))
340         END DO
341      END DO
342      !
343   END SUBROUTINE bdytide_update
344
345
346   SUBROUTINE bdy_dta_tides( kt, kit, pt_offset )
347      !!----------------------------------------------------------------------
348      !!                 ***  SUBROUTINE bdy_dta_tides  ***
349      !!               
350      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
351      !!               
352      !!----------------------------------------------------------------------
353      INTEGER,           INTENT(in) ::   kt          ! Main timestep counter
354      INTEGER, OPTIONAL, INTENT(in) ::   kit         ! Barotropic timestep counter (for timesplitting option)
355      REAL(wp),OPTIONAL, INTENT(in) ::   pt_offset   ! time offset in units of timesteps. NB. if kit
356      !                                              ! is present then units = subcycle timesteps.
357      !
358      LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step
359      INTEGER  ::   itide, ib_bdy, ib, igrd   ! loop indices
360      INTEGER, DIMENSION(jpbgrd)   ::   ilen0 
361      INTEGER, DIMENSION(1:jpbgrd) ::   nblen, nblenrim  ! short cuts
362      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist, zt_offset   
363      !!----------------------------------------------------------------------
364      !
365      lk_first_btstp=.TRUE.
366      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
367
368      zt_offset = 0._wp
369      IF( PRESENT(pt_offset) )   zt_offset = pt_offset
370     
371      ! Absolute time from model initialization:   
372      IF( PRESENT(kit) ) THEN 
373         z_arg = ( REAL(kt, wp) + ( REAL(kit, wp) + zt_offset - 1. ) / REAL(nn_baro, wp) ) * rdt
374      ELSE                             
375         z_arg = ( REAL(kt, wp) + zt_offset ) * rdt
376      ENDIF
377
378      ! Linear ramp on tidal component at open boundaries
379      zramp = 1.
380      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - REAL(nit000, wp)*rdt)/(rdttideramp*rday),0.),1.)
381
382      DO ib_bdy = 1,nb_bdy
383         !
384         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
385            !
386            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
387            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
388            !
389            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:)
390            ELSE                                   ;   ilen0(:) = nblenrim(:)
391            ENDIF     
392            !
393            ! We refresh nodal factors every day below
394            ! This should be done somewhere else
395            IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
396               !
397               kt_tide = kt - NINT((REAL(nsec_day,wp) - 0.5_wp * rdt)/rdt)
398               !
399               IF(lwp) THEN
400               WRITE(numout,*)
401               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
402               WRITE(numout,*) '~~~~~~~~~~~~~~ '
403               ENDIF
404               !
405               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
406               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
407               !
408            ENDIF
409            zoff = REAL(-kt_tide,wp) * rdt ! time offset relative to nodal factor computation time
410            !
411            ! If time splitting, initialize arrays from slow varying open boundary data:
412            IF ( PRESENT(kit) ) THEN           
413               IF ( dta_bdy(ib_bdy)%lneed_ssh   ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
414               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
415               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
416            ENDIF
417            !
418            ! Update open boundary data arrays:
419            DO itide = 1, nb_harmo
420               !
421               z_sarg = (z_arg + zoff) * omega_tide(itide)
422               z_cost = zramp * COS( z_sarg )
423               z_sist = zramp * SIN( z_sarg )
424               !
425               IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN
426                  igrd=1                              ! SSH on tracer grid
427                  DO ib = 1, ilen0(igrd)
428                     dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
429                        &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
430                        &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
431                  END DO
432               ENDIF
433               !
434               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
435                  igrd=2                              ! U grid
436                  DO ib = 1, ilen0(igrd)
437                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
438                        &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
439                        &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
440                  END DO
441                  igrd=3                              ! V grid
442                  DO ib = 1, ilen0(igrd) 
443                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
444                        &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
445                        &                        tides(ib_bdy)%v(ib,itide,2)*z_sist )
446                  END DO
447               ENDIF
448            END DO             
449         END IF
450      END DO
451      !
452   END SUBROUTINE bdy_dta_tides
453
454
455   SUBROUTINE tide_init_elevation( idx, td )
456      !!----------------------------------------------------------------------
457      !!                 ***  ROUTINE tide_init_elevation  ***
458      !!----------------------------------------------------------------------
459      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
460      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
461      !
462      INTEGER ::   itide, igrd, ib       ! dummy loop indices
463      INTEGER, DIMENSION(1) ::   ilen0   ! length of boundary data (from OBC arrays)
464      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
465      !!----------------------------------------------------------------------
466      !
467      igrd=1   
468                              ! SSH on tracer grid.
469      ilen0(1) =  SIZE(td%ssh0(:,1,1))
470      !
471      ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) )
472      !
473      DO itide = 1, nb_harmo
474         DO ib = 1, ilen0(igrd)
475            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
476            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
477         END DO
478         DO ib = 1 , ilen0(igrd)
479            mod_tide(ib)=mod_tide(ib)*ftide(itide)
480            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
481         ENDDO
482         DO ib = 1 , ilen0(igrd)
483            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
484            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
485         ENDDO
486      END DO
487      !
488      DEALLOCATE( mod_tide, phi_tide )
489      !
490   END SUBROUTINE tide_init_elevation
491
492
493   SUBROUTINE tide_init_velocities( idx, td )
494      !!----------------------------------------------------------------------
495      !!                 ***  ROUTINE tide_init_elevation  ***
496      !!----------------------------------------------------------------------
497      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
498      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
499      !
500      INTEGER ::   itide, igrd, ib       ! dummy loop indices
501      INTEGER, DIMENSION(3) ::   ilen0   ! length of boundary data (from OBC arrays)
502      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
503      !!----------------------------------------------------------------------
504      !
505      ilen0(2) =  SIZE(td%u0(:,1,1))
506      ilen0(3) =  SIZE(td%v0(:,1,1))
507      !
508      igrd=2                                 ! U grid.
509      !
510      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
511      !
512      DO itide = 1, nb_harmo
513         DO ib = 1, ilen0(igrd)
514            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
515            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
516         END DO
517         DO ib = 1, ilen0(igrd)
518            mod_tide(ib)=mod_tide(ib)*ftide(itide)
519            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
520         ENDDO
521         DO ib = 1, ilen0(igrd)
522            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
523            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
524         ENDDO
525      END DO
526      !
527      DEALLOCATE( mod_tide , phi_tide )
528      !
529      igrd=3                                 ! V grid.
530      !
531      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
532
533      DO itide = 1, nb_harmo
534         DO ib = 1, ilen0(igrd)
535            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
536            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
537         END DO
538         DO ib = 1, ilen0(igrd)
539            mod_tide(ib)=mod_tide(ib)*ftide(itide)
540            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
541         ENDDO
542         DO ib = 1, ilen0(igrd)
543            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
544            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
545         ENDDO
546      END DO
547      !
548      DEALLOCATE( mod_tide, phi_tide )
549      !
550  END SUBROUTINE tide_init_velocities
551
552   !!======================================================================
553END MODULE bdytides
554
Note: See TracBrowser for help on using the repository browser.