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

source: NEMO/branches/UKMO/r14075_read_tides_fix/src/OCE/BDY/bdytides.F90 @ 15654

Last change on this file since 15654 was 15654, checked in by jcastill, 2 years ago

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