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 branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/2013/dev_r3858_NOC_ZTC/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 4263

Last change on this file since 4263 was 4263, checked in by acc, 10 years ago

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