New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
bdytides.F90 in trunk/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

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

Last change on this file since 5084 was 5084, checked in by timgraham, 10 years ago

Fix for ticket 1366
tide phase speeds were being printed before being initialised in bdytide_init. As they are printed anyway in sbc_tide I have simply removed the print from bdytide_init

  • Property svn:keywords set to Id
File size: 29.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#if defined key_bdy
14   !!----------------------------------------------------------------------
15   !!   'key_bdy'     Open Boundary Condition
16   !!----------------------------------------------------------------------
17   !!   PUBLIC
18   !!      bdytide_init     : read of namelist and initialisation of tidal harmonics data
19   !!      tide_update   : calculation of tidal forcing at each timestep
20   !!----------------------------------------------------------------------
21   USE timing          ! Timing
22   USE oce             ! ocean dynamics and tracers
23   USE dom_oce         ! ocean space and time domain
24   USE iom
25   USE in_out_manager  ! I/O units
26   USE phycst          ! physical constants
27   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
28   USE bdy_par         ! Unstructured boundary parameters
29   USE bdy_oce         ! ocean open boundary conditions
30   USE daymod          ! calendar
31   USE wrk_nemo        ! Memory allocation
32   USE tideini
33!   USE tide_mod       ! Useless ??
34   USE fldread, ONLY: fld_map
35   USE dynspg_oce, ONLY: lk_dynspg_ts
36
37   IMPLICIT NONE
38   PRIVATE
39
40   PUBLIC   bdytide_init     ! routine called in bdy_init
41   PUBLIC   bdytide_update   ! routine called in bdy_dta
42   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts
43
44   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
45      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0       !: Tidal constituents : SSH0 (read in file)
46      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0         !: Tidal constituents : U0   (read in file)
47      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v0         !: Tidal constituents : V0   (read in file)
48      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh        !: Tidal constituents : SSH  (after nodal cor.)
49      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u          !: Tidal constituents : U    (after nodal cor.)
50      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   v          !: Tidal constituents : V    (after nodal cor.)
51   END TYPE TIDES_DATA
52
53!$AGRIF_DO_NOT_TREAT
54   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
55!$AGRIF_END_DO_NOT_TREAT
56   TYPE(OBC_DATA)  , PRIVATE, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component)
57
58   !!----------------------------------------------------------------------
59   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
60   !! $Id$
61   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
62   !!----------------------------------------------------------------------
63CONTAINS
64
65   SUBROUTINE bdytide_init
66      !!----------------------------------------------------------------------
67      !!                    ***  SUBROUTINE bdytide_init  ***
68      !!                     
69      !! ** Purpose : - Read in namelist for tides and initialise external
70      !!                tidal harmonics data
71      !!
72      !!----------------------------------------------------------------------
73      !! namelist variables
74      !!-------------------
75      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
76      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data
77      LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data
78      !!
79      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices
80      INTEGER                                   ::   ii, ij              !: dummy loop indices
81      INTEGER                                   ::   inum, igrd
82      INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays)
83      INTEGER, POINTER, DIMENSION(:)            ::   nblen, nblenrim     ! short cuts
84      INTEGER                                   ::   ios                 ! Local integer output status for namelist read
85      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file
86      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data
87      REAL(wp), POINTER, DIMENSION(:,:)         ::   ztr, zti            !:  "     "    "   "   "   "        "      "
88      !!
89      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
90      !!
91      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
92      !!----------------------------------------------------------------------
93
94      IF( nn_timing == 1 ) CALL timing_start('bdytide_init')
95
96      IF (nb_bdy>0) THEN
97         IF(lwp) WRITE(numout,*)
98         IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
99         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
100      ENDIF
101
102      REWIND(numnam_cfg)
103
104      DO ib_bdy = 1, nb_bdy
105         IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN
106
107            td => tides(ib_bdy)
108            nblen => idx_bdy(ib_bdy)%nblen
109            nblenrim => idx_bdy(ib_bdy)%nblenrim
110
111            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
112            filtide(:) = ''
113
114            ! Don't REWIND here - may need to read more than one of these namelists.
115            READ  ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
116901         IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_tide in reference namelist', lwp )
117            READ  ( numnam_cfg, nambdy_tide, IOSTAT = ios, ERR = 902 )
118902         IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist', lwp )
119            IF(lwm) WRITE ( numond, nambdy_tide )
120            !                                               ! Parameter control and print
121            IF(lwp) WRITE(numout,*) '  '
122            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
123            IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta
124            IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj
125            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
126            IF(lwp) THEN
127                    WRITE(numout,*) '             Tidal components: ' 
128               DO itide = 1, nb_harmo
129                  WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide 
130               END DO
131            ENDIF
132            IF(lwp) WRITE(numout,*) ' '
133
134            ! Allocate space for tidal harmonics data - get size from OBC data arrays
135            ! -----------------------------------------------------------------------
136
137            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
138            ! relaxation area     
139            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN
140               ilen0(:)=nblen(:)
141            ELSE
142               ilen0(:)=nblenrim(:)
143            ENDIF
144
145            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
146            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
147
148            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
149            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
150
151            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
152            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
153
154            td%ssh0(:,:,:) = 0._wp
155            td%ssh (:,:,:) = 0._wp
156            td%u0  (:,:,:) = 0._wp
157            td%u   (:,:,:) = 0._wp
158            td%v0  (:,:,:) = 0._wp
159            td%v   (:,:,:) = 0._wp
160
161            IF (ln_bdytide_2ddta) THEN
162               ! It is assumed that each data file contains all complex harmonic amplitudes
163               ! given on the data domain (ie global, jpidta x jpjdta)
164               !
165               CALL wrk_alloc( jpi, jpj, zti, ztr )
166               !
167               ! SSH fields
168               clfile = TRIM(filtide)//'_grid_T.nc'
169               CALL iom_open (clfile , inum ) 
170               igrd = 1                       ! Everything is at T-points here
171               DO itide = 1, nb_harmo
172                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
173                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 
174                  DO ib = 1, ilen0(igrd)
175                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
176                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
177                     td%ssh0(ib,itide,1) = ztr(ii,ij)
178                     td%ssh0(ib,itide,2) = zti(ii,ij)
179                  END DO
180               END DO
181               CALL iom_close( inum )
182               !
183               ! U fields
184               clfile = TRIM(filtide)//'_grid_U.nc'
185               CALL iom_open (clfile , inum ) 
186               igrd = 2                       ! Everything is at U-points here
187               DO itide = 1, nb_harmo
188                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
189                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
190                  DO ib = 1, ilen0(igrd)
191                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
192                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
193                     td%u0(ib,itide,1) = ztr(ii,ij)
194                     td%u0(ib,itide,2) = zti(ii,ij)
195                  END DO
196               END DO
197               CALL iom_close( inum )
198               !
199               ! V fields
200               clfile = TRIM(filtide)//'_grid_V.nc'
201               CALL iom_open (clfile , inum ) 
202               igrd = 3                       ! Everything is at V-points here
203               DO itide = 1, nb_harmo
204                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
205                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
206                  DO ib = 1, ilen0(igrd)
207                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
208                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
209                     td%v0(ib,itide,1) = ztr(ii,ij)
210                     td%v0(ib,itide,2) = zti(ii,ij)
211                  END DO
212               END DO 
213               CALL iom_close( inum )
214               !
215               CALL wrk_dealloc( jpi, jpj, ztr, zti ) 
216               !
217            ELSE           
218               !
219               ! Read tidal data only on bdy segments
220               !
221               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
222
223               ! Open files and read in tidal forcing data
224               ! -----------------------------------------
225
226               DO itide = 1, nb_harmo
227                  !                                                              ! SSH fields
228                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
229                  CALL iom_open( clfile, inum )
230                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
231                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
232                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
233                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
234                  CALL iom_close( inum )
235                  !                                                              ! U fields
236                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
237                  CALL iom_open( clfile, inum )
238                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
239                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
240                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
241                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
242                  CALL iom_close( inum )
243                  !                                                              ! V fields
244                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
245                  CALL iom_open( clfile, inum )
246                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
247                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
248                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
249                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
250                  CALL iom_close( inum )
251                  !
252               END DO ! end loop on tidal components
253               !
254               DEALLOCATE( dta_read )
255            ENDIF ! ln_bdytide_2ddta=.true.
256            !
257            IF ( ln_bdytide_conj ) THEN ! assume complex conjugate in data files
258               td%ssh0(:,:,2) = - td%ssh0(:,:,2)
259               td%u0  (:,:,2) = - td%u0  (:,:,2)
260               td%v0  (:,:,2) = - td%v0  (:,:,2)
261            ENDIF
262            !
263            IF ( lk_dynspg_ts ) THEN ! Allocate arrays to save slowly varying boundary data during
264                                     ! time splitting integration
265               ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
266               ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
267               ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
268               dta_bdy_s(ib_bdy)%ssh(:) = 0.e0
269               dta_bdy_s(ib_bdy)%u2d(:) = 0.e0
270               dta_bdy_s(ib_bdy)%v2d(:) = 0.e0
271            ENDIF
272            !
273         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2
274         !
275      END DO ! loop on ib_bdy
276
277      IF( nn_timing == 1 ) CALL timing_stop('bdytide_init')
278
279   END SUBROUTINE bdytide_init
280
281   SUBROUTINE bdytide_update ( kt, idx, dta, td, jit, time_offset )
282      !!----------------------------------------------------------------------
283      !!                 ***  SUBROUTINE bdytide_update  ***
284      !!               
285      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
286      !!               
287      !!----------------------------------------------------------------------
288      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter
289      TYPE(OBC_INDEX), INTENT( in )    ::   idx         ! OBC indices
290      TYPE(OBC_DATA),  INTENT(inout)   ::   dta         ! OBC external data
291      TYPE(TIDES_DATA),INTENT( inout ) ::   td          ! tidal harmonics data
292      INTEGER,INTENT(in),OPTIONAL      ::   jit         ! Barotropic timestep counter (for timesplitting option)
293      INTEGER,INTENT( in ), OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if jit
294                                                        ! is present then units = subcycle timesteps.
295                                                        ! time_offset = 0  => get data at "now"    time level
296                                                        ! time_offset = -1 => get data at "before" time level
297                                                        ! time_offset = +1 => get data at "after"  time level
298                                                        ! etc.
299      !!
300      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays)
301      INTEGER                          :: itide, igrd, ib   ! dummy loop indices
302      INTEGER                          :: time_add          ! time offset in units of timesteps
303      REAL(wp)                         :: z_arg, z_sarg, zflag, zramp     
304      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
305      !!----------------------------------------------------------------------
306
307      IF( nn_timing == 1 ) CALL timing_start('bdytide_update')
308
309      ilen0(1) =  SIZE(td%ssh(:,1,1))
310      ilen0(2) =  SIZE(td%u(:,1,1))
311      ilen0(3) =  SIZE(td%v(:,1,1))
312
313      zflag=1
314      IF ( PRESENT(jit) ) THEN
315        IF ( jit /= 1 ) zflag=0
316      ENDIF
317
318      IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. zflag==1 ) THEN
319        !
320        kt_tide = kt
321        !
322        IF(lwp) THEN
323           WRITE(numout,*)
324           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
325           WRITE(numout,*) '~~~~~~~~~~~~~~ '
326        ENDIF
327        !
328        CALL tide_init_elevation ( idx, td )
329        CALL tide_init_velocities( idx, td )
330        !
331      ENDIF
332
333      time_add = 0
334      IF( PRESENT(time_offset) ) THEN
335         time_add = time_offset
336      ENDIF
337         
338      IF( PRESENT(jit) ) THEN 
339         z_arg = ((kt-kt_tide) * rdt + (jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) )
340      ELSE                             
341         z_arg = ((kt-kt_tide)+time_add) * rdt
342      ENDIF
343
344      ! Linear ramp on tidal component at open boundaries
345      zramp = 1._wp
346      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp)
347
348      DO itide = 1, nb_harmo
349         z_sarg = z_arg * omega_tide(itide)
350         z_cost(itide) = COS( z_sarg )
351         z_sist(itide) = SIN( z_sarg )
352      END DO
353
354      DO itide = 1, nb_harmo
355         igrd=1                              ! SSH on tracer grid
356         DO ib = 1, ilen0(igrd)
357            dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide))
358         END DO
359         igrd=2                              ! U grid
360         DO ib = 1, ilen0(igrd)
361            dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide))
362         END DO
363         igrd=3                              ! V grid
364         DO ib = 1, ilen0(igrd) 
365            dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide))
366         END DO
367      END DO
368      !
369      IF( nn_timing == 1 ) CALL timing_stop('bdytide_update')
370      !
371   END SUBROUTINE bdytide_update
372
373   SUBROUTINE bdy_dta_tides( kt, kit, time_offset )
374      !!----------------------------------------------------------------------
375      !!                 ***  SUBROUTINE bdy_dta_tides  ***
376      !!               
377      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
378      !!               
379      !!----------------------------------------------------------------------
380      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter
381      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Barotropic timestep counter (for timesplitting option)
382      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if kit
383                                                        ! is present then units = subcycle timesteps.
384                                                        ! time_offset = 0  => get data at "now"    time level
385                                                        ! time_offset = -1 => get data at "before" time level
386                                                        ! time_offset = +1 => get data at "after"  time level
387                                                        ! etc.
388      !!
389      LOGICAL  :: lk_first_btstp  ! =.TRUE. if time splitting and first barotropic step
390      INTEGER,          DIMENSION(jpbgrd) :: ilen0 
391      INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim  ! short cuts
392      INTEGER  :: itide, ib_bdy, ib, igrd                     ! loop indices
393      INTEGER  :: time_add                                    ! time offset in units of timesteps
394      REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist     
395      !!----------------------------------------------------------------------
396
397      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_tides')
398
399      lk_first_btstp=.TRUE.
400      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
401
402      time_add = 0
403      IF( PRESENT(time_offset) ) THEN
404         time_add = time_offset
405      ENDIF
406     
407      ! Absolute time from model initialization:   
408      IF( PRESENT(kit) ) THEN 
409         z_arg = ( kt + (kit+0.5_wp*(time_add-1)) / REAL(nn_baro,wp) ) * rdt
410      ELSE                             
411         z_arg = ( kt + time_add ) * rdt
412      ENDIF
413
414      ! Linear ramp on tidal component at open boundaries
415      zramp = 1.
416      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.)
417
418      DO ib_bdy = 1,nb_bdy
419
420         IF ( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN
421
422            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
423            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
424
425            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN
426               ilen0(:)=nblen(:)
427            ELSE
428               ilen0(:)=nblenrim(:)
429            ENDIF     
430
431            ! We refresh nodal factors every day below
432            ! This should be done somewhere else
433            IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. lk_first_btstp ) THEN
434               !
435               kt_tide = kt               
436               !
437               IF(lwp) THEN
438               WRITE(numout,*)
439               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
440               WRITE(numout,*) '~~~~~~~~~~~~~~ '
441               ENDIF
442               !
443               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
444               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
445               !
446            ENDIF
447            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
448            !
449            ! If time splitting, save data at first barotropic iteration
450            IF ( PRESENT(kit) ) THEN
451               IF ( lk_first_btstp ) THEN ! Save slow varying open boundary data:
452                  IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy(ib_bdy)%ssh(1:ilen0(1))
453                  IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy(ib_bdy)%u2d(1:ilen0(2))
454                  IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy(ib_bdy)%v2d(1:ilen0(3))
455
456               ELSE ! Initialize arrays from slow varying open boundary data:           
457                  IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
458                  IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
459                  IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
460               ENDIF
461            ENDIF
462            !
463            ! Update open boundary data arrays:
464            DO itide = 1, nb_harmo
465               !
466               z_sarg = (z_arg + zoff) * omega_tide(itide)
467               z_cost = zramp * COS( z_sarg )
468               z_sist = zramp * SIN( z_sarg )
469               !
470               IF ( dta_bdy(ib_bdy)%ll_ssh ) THEN
471                  igrd=1                              ! SSH on tracer grid
472                  DO ib = 1, ilen0(igrd)
473                     dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
474                        &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
475                        &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
476                  END DO
477               ENDIF
478               !
479               IF ( dta_bdy(ib_bdy)%ll_u2d ) THEN
480                  igrd=2                              ! U grid
481                  DO ib = 1, ilen0(igrd)
482                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
483                        &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
484                        &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
485                  END DO
486               ENDIF
487               !
488               IF ( dta_bdy(ib_bdy)%ll_v2d ) THEN
489                  igrd=3                              ! V grid
490                  DO ib = 1, ilen0(igrd) 
491                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
492                        &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
493                        &                        tides(ib_bdy)%v(ib,itide,2)*z_sist )
494                  END DO
495               ENDIF
496            END DO             
497         END IF
498      END DO
499      !
500      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_tides')
501      !
502   END SUBROUTINE bdy_dta_tides
503
504   SUBROUTINE tide_init_elevation( 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      !! * Local declarations
511      INTEGER, DIMENSION(1)            ::   ilen0       !: length of boundary data (from OBC arrays)
512      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
513      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
514
515      igrd=1   
516                              ! SSH on tracer grid.
517   
518      ilen0(1) =  SIZE(td%ssh0(:,1,1))
519
520      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
521
522      DO itide = 1, nb_harmo
523         DO ib = 1, ilen0(igrd)
524            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
525            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
526         END DO
527         DO ib = 1 , ilen0(igrd)
528            mod_tide(ib)=mod_tide(ib)*ftide(itide)
529            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
530         ENDDO
531         DO ib = 1 , ilen0(igrd)
532            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
533            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
534         ENDDO
535      END DO
536
537      DEALLOCATE(mod_tide,phi_tide)
538
539   END SUBROUTINE tide_init_elevation
540
541   SUBROUTINE tide_init_velocities( idx, td )
542      !!----------------------------------------------------------------------
543      !!                 ***  ROUTINE tide_init_elevation  ***
544      !!----------------------------------------------------------------------
545      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
546      TYPE(TIDES_DATA),INTENT( inout )      ::   td      ! tidal harmonics data
547      !! * Local declarations
548      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays)
549      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
550      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
551
552      ilen0(2) =  SIZE(td%u0(:,1,1))
553      ilen0(3) =  SIZE(td%v0(:,1,1))
554
555      igrd=2                                 ! U grid.
556
557      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
558
559      DO itide = 1, nb_harmo
560         DO ib = 1, ilen0(igrd)
561            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
562            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
563         END DO
564         DO ib = 1, ilen0(igrd)
565            mod_tide(ib)=mod_tide(ib)*ftide(itide)
566            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
567         ENDDO
568         DO ib = 1, ilen0(igrd)
569            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
570            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
571         ENDDO
572      END DO
573
574      DEALLOCATE(mod_tide,phi_tide)
575
576      igrd=3                                 ! V grid.
577
578      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
579
580      DO itide = 1, nb_harmo
581         DO ib = 1, ilen0(igrd)
582            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
583            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
584         END DO
585         DO ib = 1, ilen0(igrd)
586            mod_tide(ib)=mod_tide(ib)*ftide(itide)
587            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
588         ENDDO
589         DO ib = 1, ilen0(igrd)
590            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
591            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
592         ENDDO
593      END DO
594
595      DEALLOCATE(mod_tide,phi_tide)
596
597  END SUBROUTINE tide_init_velocities
598#else
599   !!----------------------------------------------------------------------
600   !!   Dummy module         NO Unstruct Open Boundary Conditions for tides
601   !!----------------------------------------------------------------------
602CONTAINS
603   SUBROUTINE bdytide_init             ! Empty routine
604      WRITE(*,*) 'bdytide_init: You should not have seen this print! error?'
605   END SUBROUTINE bdytide_init
606   SUBROUTINE bdytide_update( kt, jit )   ! Empty routine
607      WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit
608   END SUBROUTINE bdytide_update
609   SUBROUTINE bdy_dta_tides( kt, kit, time_offset )     ! Empty routine
610      INTEGER, INTENT( in )            ::   kt          ! Dummy argument empty routine     
611      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Dummy argument empty routine
612      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! Dummy argument empty routine
613      WRITE(*,*) 'bdy_dta_tides: You should not have seen this print! error?', kt, jit
614   END SUBROUTINE bdy_dta_tides
615#endif
616
617   !!======================================================================
618END MODULE bdytides
619
Note: See TracBrowser for help on using the repository browser.