source: branches/UKMO/r6232_collate_bgc_diagnostics/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 11132

Last change on this file since 11132 was 11132, checked in by jcastill, 17 months ago

Remove svn keywords

File size: 29.7 KB
Line 
1MODULE bdytides
2   !!======================================================================
3   !!                       ***  MODULE  bdytides  ***
4   !! Ocean dynamics:   Tidal forcing at open boundaries
5   !!======================================================================
6   !! History :  2.0  !  2007-01  (D.Storkey)  Original code
7   !!            2.3  !  2008-01  (J.Holt)  Add date correction. Origins POLCOMS v6.3 2007
8   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
9   !!            3.3  !  2010-09  (D.Storkey and E.O'Dea)  bug fixes
10   !!            3.4  !  2012-09  (G. Reffray and J. Chanut) New inputs + mods
11   !!            3.5  !  2013-07  (J. Chanut) Compliant with time splitting changes
12   !!----------------------------------------------------------------------
13#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
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      TYPE(MAP_POINTER), DIMENSION(jpbgrd)      ::   ibmap_ptr           !: array of pointers to nbmap
91      !!
92      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
93      !!----------------------------------------------------------------------
94
95      IF( nn_timing == 1 ) CALL timing_start('bdytide_init')
96
97      IF (nb_bdy>0) THEN
98         IF(lwp) WRITE(numout,*)
99         IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
100         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
101      ENDIF
102
103      REWIND(numnam_cfg)
104
105      DO ib_bdy = 1, nb_bdy
106         IF( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN
107
108            td => tides(ib_bdy)
109            nblen => idx_bdy(ib_bdy)%nblen
110            nblenrim => idx_bdy(ib_bdy)%nblenrim
111
112            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
113            filtide(:) = ''
114
115            ! Don't REWIND here - may need to read more than one of these namelists.
116            READ  ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
117901         IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_tide in reference namelist', lwp )
118            READ  ( numnam_cfg, nambdy_tide, IOSTAT = ios, ERR = 902 )
119902         IF( ios /= 0 ) CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist', lwp )
120            IF(lwm) WRITE ( numond, nambdy_tide )
121            !                                               ! Parameter control and print
122            IF(lwp) WRITE(numout,*) '  '
123            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
124            IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta
125            IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj
126            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
127            IF(lwp) THEN
128                    WRITE(numout,*) '             Tidal components: ' 
129               DO itide = 1, nb_harmo
130                  WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide 
131               END DO
132            ENDIF
133            IF(lwp) WRITE(numout,*) ' '
134
135            ! Allocate space for tidal harmonics data - get size from OBC data arrays
136            ! -----------------------------------------------------------------------
137
138            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
139            ! relaxation area     
140            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN
141               ilen0(:)=nblen(:)
142            ELSE
143               ilen0(:)=nblenrim(:)
144            ENDIF
145
146            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
147            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
148
149            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
150            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
151
152            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
153            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
154
155            td%ssh0(:,:,:) = 0._wp
156            td%ssh (:,:,:) = 0._wp
157            td%u0  (:,:,:) = 0._wp
158            td%u   (:,:,:) = 0._wp
159            td%v0  (:,:,:) = 0._wp
160            td%v   (:,:,:) = 0._wp
161
162            IF (ln_bdytide_2ddta) THEN
163               ! It is assumed that each data file contains all complex harmonic amplitudes
164               ! given on the data domain (ie global, jpidta x jpjdta)
165               !
166               CALL wrk_alloc( jpi, jpj, zti, ztr )
167               !
168               ! SSH fields
169               clfile = TRIM(filtide)//'_grid_T.nc'
170               CALL iom_open (clfile , inum ) 
171               igrd = 1                       ! Everything is at T-points here
172               DO itide = 1, nb_harmo
173                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
174                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_z2', zti(:,:) ) 
175                  DO ib = 1, ilen0(igrd)
176                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
177                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
178                     td%ssh0(ib,itide,1) = ztr(ii,ij)
179                     td%ssh0(ib,itide,2) = zti(ii,ij)
180                  END DO
181               END DO
182               CALL iom_close( inum )
183               !
184               ! U fields
185               clfile = TRIM(filtide)//'_grid_U.nc'
186               CALL iom_open (clfile , inum ) 
187               igrd = 2                       ! Everything is at U-points here
188               DO itide = 1, nb_harmo
189                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
190                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
191                  DO ib = 1, ilen0(igrd)
192                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
193                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
194                     td%u0(ib,itide,1) = ztr(ii,ij)
195                     td%u0(ib,itide,2) = zti(ii,ij)
196                  END DO
197               END DO
198               CALL iom_close( inum )
199               !
200               ! V fields
201               clfile = TRIM(filtide)//'_grid_V.nc'
202               CALL iom_open (clfile , inum ) 
203               igrd = 3                       ! Everything is at V-points here
204               DO itide = 1, nb_harmo
205                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
206                  CALL iom_get  ( inum, jpdom_data, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
207                  DO ib = 1, ilen0(igrd)
208                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
209                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
210                     td%v0(ib,itide,1) = ztr(ii,ij)
211                     td%v0(ib,itide,2) = zti(ii,ij)
212                  END DO
213               END DO 
214               CALL iom_close( inum )
215               !
216               CALL wrk_dealloc( jpi, jpj, ztr, zti ) 
217               !
218            ELSE           
219               !
220               ! Read tidal data only on bdy segments
221               !
222               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
223               !
224               ! Set map structure
225               ibmap_ptr(1)%ptr => idx_bdy(ib_bdy)%nbmap(:,1)
226               ibmap_ptr(1)%ll_unstruc = ln_coords_file(ib_bdy)
227               ibmap_ptr(2)%ptr => idx_bdy(ib_bdy)%nbmap(:,2)
228               ibmap_ptr(2)%ll_unstruc = ln_coords_file(ib_bdy)
229               ibmap_ptr(3)%ptr => idx_bdy(ib_bdy)%nbmap(:,3)
230               ibmap_ptr(3)%ll_unstruc = ln_coords_file(ib_bdy)
231
232               ! Open files and read in tidal forcing data
233               ! -----------------------------------------
234
235               DO itide = 1, nb_harmo
236                  !                                                              ! SSH fields
237                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
238                  CALL iom_open( clfile, inum )
239                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1,  ibmap_ptr(1) )
240                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
241                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1,  ibmap_ptr(1) )
242                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
243                  CALL iom_close( inum )
244                  !                                                              ! U fields
245                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
246                  CALL iom_open( clfile, inum )
247                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, ibmap_ptr(2) )
248                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
249                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, ibmap_ptr(2) )
250                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
251                  CALL iom_close( inum )
252                  !                                                              ! V fields
253                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
254                  CALL iom_open( clfile, inum )
255                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, ibmap_ptr(3) )
256                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
257                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, ibmap_ptr(3) )
258                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
259                  CALL iom_close( inum )
260                  !
261               END DO ! end loop on tidal components
262               !
263               DEALLOCATE( dta_read )
264            ENDIF ! ln_bdytide_2ddta=.true.
265            !
266            IF ( ln_bdytide_conj ) THEN ! assume complex conjugate in data files
267               td%ssh0(:,:,2) = - td%ssh0(:,:,2)
268               td%u0  (:,:,2) = - td%u0  (:,:,2)
269               td%v0  (:,:,2) = - td%v0  (:,:,2)
270            ENDIF
271            !
272            IF ( lk_dynspg_ts ) THEN ! Allocate arrays to save slowly varying boundary data during
273                                     ! time splitting integration
274               ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
275               ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
276               ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
277               dta_bdy_s(ib_bdy)%ssh(:) = 0.e0
278               dta_bdy_s(ib_bdy)%u2d(:) = 0.e0
279               dta_bdy_s(ib_bdy)%v2d(:) = 0.e0
280            ENDIF
281            !
282         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2
283         !
284      END DO ! loop on ib_bdy
285
286      IF( nn_timing == 1 ) CALL timing_stop('bdytide_init')
287
288   END SUBROUTINE bdytide_init
289
290   SUBROUTINE bdytide_update ( kt, idx, dta, td, jit, time_offset )
291      !!----------------------------------------------------------------------
292      !!                 ***  SUBROUTINE bdytide_update  ***
293      !!               
294      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
295      !!               
296      !!----------------------------------------------------------------------
297      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter
298      TYPE(OBC_INDEX), INTENT( in )    ::   idx         ! OBC indices
299      TYPE(OBC_DATA),  INTENT(inout)   ::   dta         ! OBC external data
300      TYPE(TIDES_DATA),INTENT( inout ) ::   td          ! tidal harmonics data
301      INTEGER,INTENT(in),OPTIONAL      ::   jit         ! Barotropic timestep counter (for timesplitting option)
302      INTEGER,INTENT( in ), OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if jit
303                                                        ! is present then units = subcycle timesteps.
304                                                        ! time_offset = 0  => get data at "now"    time level
305                                                        ! time_offset = -1 => get data at "before" time level
306                                                        ! time_offset = +1 => get data at "after"  time level
307                                                        ! etc.
308      !!
309      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays)
310      INTEGER                          :: itide, igrd, ib   ! dummy loop indices
311      INTEGER                          :: time_add          ! time offset in units of timesteps
312      REAL(wp)                         :: z_arg, z_sarg, zflag, zramp     
313      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
314      !!----------------------------------------------------------------------
315
316      IF( nn_timing == 1 ) CALL timing_start('bdytide_update')
317
318      ilen0(1) =  SIZE(td%ssh(:,1,1))
319      ilen0(2) =  SIZE(td%u(:,1,1))
320      ilen0(3) =  SIZE(td%v(:,1,1))
321
322      zflag=1
323      IF ( PRESENT(jit) ) THEN
324        IF ( jit /= 1 ) zflag=0
325      ENDIF
326
327      IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. zflag==1 ) THEN
328        !
329        kt_tide = kt
330        !
331        IF(lwp) THEN
332           WRITE(numout,*)
333           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
334           WRITE(numout,*) '~~~~~~~~~~~~~~ '
335        ENDIF
336        !
337        CALL tide_init_elevation ( idx, td )
338        CALL tide_init_velocities( idx, td )
339        !
340      ENDIF
341
342      time_add = 0
343      IF( PRESENT(time_offset) ) THEN
344         time_add = time_offset
345      ENDIF
346         
347      IF( PRESENT(jit) ) THEN 
348         z_arg = ((kt-kt_tide) * rdt + (jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) )
349      ELSE                             
350         z_arg = ((kt-kt_tide)+time_add) * rdt
351      ENDIF
352
353      ! Linear ramp on tidal component at open boundaries
354      zramp = 1._wp
355      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp)
356
357      DO itide = 1, nb_harmo
358         z_sarg = z_arg * omega_tide(itide)
359         z_cost(itide) = COS( z_sarg )
360         z_sist(itide) = SIN( z_sarg )
361      END DO
362
363      DO itide = 1, nb_harmo
364         igrd=1                              ! SSH on tracer grid
365         DO ib = 1, ilen0(igrd)
366            dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide))
367         END DO
368         igrd=2                              ! U grid
369         DO ib = 1, ilen0(igrd)
370            dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide))
371         END DO
372         igrd=3                              ! V grid
373         DO ib = 1, ilen0(igrd) 
374            dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide))
375         END DO
376      END DO
377      !
378      IF( nn_timing == 1 ) CALL timing_stop('bdytide_update')
379      !
380   END SUBROUTINE bdytide_update
381
382   SUBROUTINE bdy_dta_tides( kt, kit, time_offset )
383      !!----------------------------------------------------------------------
384      !!                 ***  SUBROUTINE bdy_dta_tides  ***
385      !!               
386      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
387      !!               
388      !!----------------------------------------------------------------------
389      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter
390      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Barotropic timestep counter (for timesplitting option)
391      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if kit
392                                                        ! is present then units = subcycle timesteps.
393                                                        ! time_offset = 0  => get data at "now"    time level
394                                                        ! time_offset = -1 => get data at "before" time level
395                                                        ! time_offset = +1 => get data at "after"  time level
396                                                        ! etc.
397      !!
398      LOGICAL  :: lk_first_btstp  ! =.TRUE. if time splitting and first barotropic step
399      INTEGER,          DIMENSION(jpbgrd) :: ilen0 
400      INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim  ! short cuts
401      INTEGER  :: itide, ib_bdy, ib, igrd                     ! loop indices
402      INTEGER  :: time_add                                    ! time offset in units of timesteps
403      REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist     
404      !!----------------------------------------------------------------------
405
406      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_tides')
407
408      lk_first_btstp=.TRUE.
409      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
410
411      time_add = 0
412      IF( PRESENT(time_offset) ) THEN
413         time_add = time_offset
414      ENDIF
415     
416      ! Absolute time from model initialization:   
417      IF( PRESENT(kit) ) THEN 
418         z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt
419      ELSE                             
420         z_arg = ( kt + time_add ) * rdt
421      ENDIF
422
423      ! Linear ramp on tidal component at open boundaries
424      zramp = 1.
425      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.)
426
427      DO ib_bdy = 1,nb_bdy
428
429         IF ( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN
430
431            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
432            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
433
434            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN
435               ilen0(:)=nblen(:)
436            ELSE
437               ilen0(:)=nblenrim(:)
438            ENDIF     
439
440            ! We refresh nodal factors every day below
441            ! This should be done somewhere else
442            IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. lk_first_btstp ) THEN
443               !
444               kt_tide = kt               
445               !
446               IF(lwp) THEN
447               WRITE(numout,*)
448               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
449               WRITE(numout,*) '~~~~~~~~~~~~~~ '
450               ENDIF
451               !
452               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
453               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
454               !
455            ENDIF
456            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
457            !
458            ! If time splitting, save data at first barotropic iteration
459            IF ( PRESENT(kit) ) THEN
460               IF ( lk_first_btstp ) THEN ! Save slow varying open boundary data:
461                  IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy(ib_bdy)%ssh(1:ilen0(1))
462                  IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy(ib_bdy)%u2d(1:ilen0(2))
463                  IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy(ib_bdy)%v2d(1:ilen0(3))
464
465               ELSE ! Initialize arrays from slow varying open boundary data:           
466                  IF ( dta_bdy(ib_bdy)%ll_ssh ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
467                  IF ( dta_bdy(ib_bdy)%ll_u2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
468                  IF ( dta_bdy(ib_bdy)%ll_v2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
469               ENDIF
470            ENDIF
471            !
472            ! Update open boundary data arrays:
473            DO itide = 1, nb_harmo
474               !
475               z_sarg = (z_arg + zoff) * omega_tide(itide)
476               z_cost = zramp * COS( z_sarg )
477               z_sist = zramp * SIN( z_sarg )
478               !
479               IF ( dta_bdy(ib_bdy)%ll_ssh ) THEN
480                  igrd=1                              ! SSH on tracer grid
481                  DO ib = 1, ilen0(igrd)
482                     dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
483                        &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
484                        &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
485                  END DO
486               ENDIF
487               !
488               IF ( dta_bdy(ib_bdy)%ll_u2d ) THEN
489                  igrd=2                              ! U grid
490                  DO ib = 1, ilen0(igrd)
491                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
492                        &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
493                        &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
494                  END DO
495               ENDIF
496               !
497               IF ( dta_bdy(ib_bdy)%ll_v2d ) THEN
498                  igrd=3                              ! V grid
499                  DO ib = 1, ilen0(igrd) 
500                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
501                        &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
502                        &                        tides(ib_bdy)%v(ib,itide,2)*z_sist )
503                  END DO
504               ENDIF
505            END DO             
506         END IF
507      END DO
508      !
509      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_tides')
510      !
511   END SUBROUTINE bdy_dta_tides
512
513   SUBROUTINE tide_init_elevation( idx, td )
514      !!----------------------------------------------------------------------
515      !!                 ***  ROUTINE tide_init_elevation  ***
516      !!----------------------------------------------------------------------
517      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
518      TYPE(TIDES_DATA),INTENT( inout )   ::   td      ! tidal harmonics data
519      !! * Local declarations
520      INTEGER, DIMENSION(1)            ::   ilen0       !: length of boundary data (from OBC arrays)
521      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
522      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
523
524      igrd=1   
525                              ! SSH on tracer grid.
526   
527      ilen0(1) =  SIZE(td%ssh0(:,1,1))
528
529      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
530
531      DO itide = 1, nb_harmo
532         DO ib = 1, ilen0(igrd)
533            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
534            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
535         END DO
536         DO ib = 1 , ilen0(igrd)
537            mod_tide(ib)=mod_tide(ib)*ftide(itide)
538            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
539         ENDDO
540         DO ib = 1 , ilen0(igrd)
541            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
542            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
543         ENDDO
544      END DO
545
546      DEALLOCATE(mod_tide,phi_tide)
547
548   END SUBROUTINE tide_init_elevation
549
550   SUBROUTINE tide_init_velocities( idx, td )
551      !!----------------------------------------------------------------------
552      !!                 ***  ROUTINE tide_init_elevation  ***
553      !!----------------------------------------------------------------------
554      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
555      TYPE(TIDES_DATA),INTENT( inout )      ::   td      ! tidal harmonics data
556      !! * Local declarations
557      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays)
558      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
559      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
560
561      ilen0(2) =  SIZE(td%u0(:,1,1))
562      ilen0(3) =  SIZE(td%v0(:,1,1))
563
564      igrd=2                                 ! U grid.
565
566      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
567
568      DO itide = 1, nb_harmo
569         DO ib = 1, ilen0(igrd)
570            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
571            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
572         END DO
573         DO ib = 1, ilen0(igrd)
574            mod_tide(ib)=mod_tide(ib)*ftide(itide)
575            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
576         ENDDO
577         DO ib = 1, ilen0(igrd)
578            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
579            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
580         ENDDO
581      END DO
582
583      DEALLOCATE(mod_tide,phi_tide)
584
585      igrd=3                                 ! V grid.
586
587      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
588
589      DO itide = 1, nb_harmo
590         DO ib = 1, ilen0(igrd)
591            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
592            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
593         END DO
594         DO ib = 1, ilen0(igrd)
595            mod_tide(ib)=mod_tide(ib)*ftide(itide)
596            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
597         ENDDO
598         DO ib = 1, ilen0(igrd)
599            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
600            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
601         ENDDO
602      END DO
603
604      DEALLOCATE(mod_tide,phi_tide)
605
606  END SUBROUTINE tide_init_velocities
607#else
608   !!----------------------------------------------------------------------
609   !!   Dummy module         NO Unstruct Open Boundary Conditions for tides
610   !!----------------------------------------------------------------------
611CONTAINS
612   SUBROUTINE bdytide_init             ! Empty routine
613      WRITE(*,*) 'bdytide_init: You should not have seen this print! error?'
614   END SUBROUTINE bdytide_init
615   SUBROUTINE bdytide_update( kt, jit )   ! Empty routine
616      WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit
617   END SUBROUTINE bdytide_update
618   SUBROUTINE bdy_dta_tides( kt, kit, time_offset )     ! Empty routine
619      INTEGER, INTENT( in )            ::   kt          ! Dummy argument empty routine     
620      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Dummy argument empty routine
621      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! Dummy argument empty routine
622      WRITE(*,*) 'bdy_dta_tides: You should not have seen this print! error?', kt, jit
623   END SUBROUTINE bdy_dta_tides
624#endif
625
626   !!======================================================================
627END MODULE bdytides
628
Note: See TracBrowser for help on using the repository browser.