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

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

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/BDY/bdytides.F90 @ 11960

Last change on this file since 11960 was 11960, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. Merge in changes from 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. (svn merge -r 11614:11954). Resolved tree conflicts and one actual conflict. Sette tested(these changes alter the ext/AGRIF reference; remember to update). See ticket #2341
  • Property svn:keywords set to Id
File size: 27.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   !!   bdytide_init  : read of namelist and initialisation of tidal harmonics data
14   !!   tide_update   : calculation of tidal forcing at each timestep
15   !!----------------------------------------------------------------------
16   USE oce            ! ocean dynamics and tracers
17   USE dom_oce        ! ocean space and time domain
18   USE phycst         ! physical constants
19   USE bdy_oce        ! ocean open boundary conditions
20   USE tideini        !
21   USE daymod         ! calendar
22   !
23   USE in_out_manager ! I/O units
24   USE iom            ! xIO server
25   USE fldread        !
26   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
27
28   IMPLICIT NONE
29   PRIVATE
30
31   PUBLIC   bdytide_init     ! routine called in bdy_init
32   PUBLIC   bdytide_update   ! routine called in bdy_dta
33   PUBLIC   bdy_dta_tides    ! routine called in dyn_spg_ts
34
35   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
36      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0     !: Tidal constituents : SSH0   (read in file)
37      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0, v0   !: Tidal constituents : U0, V0 (read in file)
38      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh      !: Tidal constituents : SSH    (after nodal cor.)
39      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u , v    !: Tidal constituents : U , V  (after nodal cor.)
40   END TYPE TIDES_DATA
41
42!$AGRIF_DO_NOT_TREAT
43   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
44!$AGRIF_END_DO_NOT_TREAT
45   TYPE(OBC_DATA)  , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component)
46
47   !!----------------------------------------------------------------------
48   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
49   !! $Id$
50   !! Software governed by the CeCILL license (see ./LICENSE)
51   !!----------------------------------------------------------------------
52CONTAINS
53
54   SUBROUTINE bdytide_init
55      !!----------------------------------------------------------------------
56      !!                    ***  SUBROUTINE bdytide_init  ***
57      !!                     
58      !! ** Purpose : - Read in namelist for tides and initialise external
59      !!                tidal harmonics data
60      !!
61      !!----------------------------------------------------------------------
62      !! namelist variables
63      !!-------------------
64      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
65      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic data
66      LOGICAL                                   ::   ln_bdytide_conj     !: If true, assume complex conjugate tidal data
67      !!
68      INTEGER                                   ::   ib_bdy, itide, ib   !: dummy loop indices
69      INTEGER                                   ::   ii, ij              !: dummy loop indices
70      INTEGER                                   ::   inum, igrd
71      INTEGER, DIMENSION(3)                     ::   ilen0       !: length of boundary data (from OBC arrays)
72      INTEGER                                   ::   ios                 ! Local integer output status for namelist read
73      INTEGER                                   ::   nbdy_rdstart, nbdy_loc
74      CHARACTER(LEN=50)                         ::   cerrmsg             !: error string
75      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file
76      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data
77      REAL(wp),ALLOCATABLE, DIMENSION(:,:)      ::   ztr, zti            !:  "     "    "   "   "   "        "      "
78      !!
79      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
80      !!
81      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta, ln_bdytide_conj
82      !!----------------------------------------------------------------------
83      !
84      IF(lwp) WRITE(numout,*)
85      IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
86      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
87
88
89      nbdy_rdstart = 1
90      DO ib_bdy = 1, nb_bdy
91         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
92            !
93            td => tides(ib_bdy)
94
95            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
96            filtide(:) = ''
97
98            READ  ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
99901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_tide in reference namelist' )
100            !
101            ! Need to support possibility of reading more than one
102            ! nambdy_tide from the namelist_cfg internal file.
103            ! Do this by finding the ib_bdy'th occurence of nambdy_tide in the
104            ! character buffer as the starting point.
105            !
106            nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_tide' )
107            IF( nbdy_loc .GT. 0 ) THEN
108               nbdy_rdstart = nbdy_rdstart + nbdy_loc
109            ELSE
110               WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',ib_bdy,' of nambdy_tide not found'
111               ios = -1
112               CALL ctl_nam ( ios , cerrmsg )
113            ENDIF
114            READ  ( numnam_cfg( MAX( 1, nbdy_rdstart - 2 ): ), nambdy_tide, IOSTAT = ios, ERR = 902)
115902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist' )
116            IF(lwm) WRITE ( numond, nambdy_tide )
117            !                                               ! Parameter control and print
118            IF(lwp) WRITE(numout,*) '  '
119            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
120            IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta
121            IF(lwp) WRITE(numout,*) '             assume complex conjugate   : ', ln_bdytide_conj
122            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
123            IF(lwp) THEN
124                    WRITE(numout,*) '             Tidal components: ' 
125               DO itide = 1, nb_harmo
126                  WRITE(numout,*)  '                 ', Wave(ntide(itide))%cname_tide 
127               END DO
128            ENDIF
129            IF(lwp) WRITE(numout,*) ' '
130
131            ! Allocate space for tidal harmonics data - get size from OBC data arrays
132            ! -----------------------------------------------------------------------
133
134            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
135            ! relaxation area     
136            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = idx_bdy(ib_bdy)%nblen   (:)
137            ELSE                                   ;   ilen0(:) = idx_bdy(ib_bdy)%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 global domain (ie global, jpiglo x jpjglo)
159               !
160               ALLOCATE( zti(jpi,jpj), ztr(jpi,jpj) )
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_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_z1', ztr(:,:) )
168                  CALL iom_get( inum, jpdom_autoglo, 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                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
173                     td%ssh0(ib,itide,1) = ztr(ii,ij)
174                     td%ssh0(ib,itide,2) = zti(ii,ij)
175                  END DO
176               END DO
177               CALL iom_close( inum )
178               !
179               ! U fields
180               clfile = TRIM(filtide)//'_grid_U.nc'
181               CALL iom_open( clfile , inum ) 
182               igrd = 2                       ! Everything is at U-points here
183               DO itide = 1, nb_harmo
184                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u1', ztr(:,:) )
185                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_u2', zti(:,:) )
186                  DO ib = 1, ilen0(igrd)
187                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
188                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
189                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
190                     td%u0(ib,itide,1) = ztr(ii,ij)
191                     td%u0(ib,itide,2) = zti(ii,ij)
192                  END DO
193               END DO
194               CALL iom_close( inum )
195               !
196               ! V fields
197               clfile = TRIM(filtide)//'_grid_V.nc'
198               CALL iom_open( clfile , inum ) 
199               igrd = 3                       ! Everything is at V-points here
200               DO itide = 1, nb_harmo
201                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v1', ztr(:,:) )
202                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(Wave(ntide(itide))%cname_tide)//'_v2', zti(:,:) )
203                  DO ib = 1, ilen0(igrd)
204                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
205                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
206                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
207                     td%v0(ib,itide,1) = ztr(ii,ij)
208                     td%v0(ib,itide,2) = zti(ii,ij)
209                  END DO
210               END DO 
211               CALL iom_close( inum )
212               !
213               DEALLOCATE( ztr, zti ) 
214               !
215            ELSE           
216               !
217               ! Read tidal data only on bdy segments
218               !
219               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
220               !
221               ! Open files and read in tidal forcing data
222               ! -----------------------------------------
223
224               DO itide = 1, nb_harmo
225                  !                                                              ! SSH fields
226                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
227                  CALL iom_open( clfile, inum )
228                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
229                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
230                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
231                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
232                  CALL iom_close( inum )
233                  !                                                              ! U fields
234                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
235                  CALL iom_open( clfile, inum )
236                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
237                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
238                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
239                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
240                  CALL iom_close( inum )
241                  !                                                              ! V fields
242                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
243                  CALL iom_open( clfile, inum )
244                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
245                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
246                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
247                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
248                  CALL iom_close( inum )
249                  !
250               END DO ! end loop on tidal components
251               !
252               DEALLOCATE( dta_read )
253               !
254            ENDIF ! ln_bdytide_2ddta=.true.
255            !
256            IF( ln_bdytide_conj ) THEN    ! assume complex conjugate in data files
257               td%ssh0(:,:,2) = - td%ssh0(:,:,2)
258               td%u0  (:,:,2) = - td%u0  (:,:,2)
259               td%v0  (:,:,2) = - td%v0  (:,:,2)
260            ENDIF
261            !
262            ! Allocate slow varying data in the case of time splitting:
263            ! Do it anyway because at this stage knowledge of free surface scheme is unknown
264            ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
265            ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
266            ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
267            dta_bdy_s(ib_bdy)%ssh(:) = 0._wp
268            dta_bdy_s(ib_bdy)%u2d(:) = 0._wp
269            dta_bdy_s(ib_bdy)%v2d(:) = 0._wp
270            !
271         ENDIF ! nn_dyn2d_dta(ib_bdy) >= 2
272         !
273      END DO ! loop on ib_bdy
274      !
275   END SUBROUTINE bdytide_init
276
277
278   SUBROUTINE bdytide_update( kt, idx, dta, td, kit, kt_offset )
279      !!----------------------------------------------------------------------
280      !!                 ***  SUBROUTINE bdytide_update  ***
281      !!               
282      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
283      !!               
284      !!----------------------------------------------------------------------
285      INTEGER          , INTENT(in   ) ::   kt          ! Main timestep counter
286      TYPE(OBC_INDEX)  , INTENT(in   ) ::   idx         ! OBC indices
287      TYPE(OBC_DATA)   , INTENT(inout) ::   dta         ! OBC external data
288      TYPE(TIDES_DATA) , INTENT(inout) ::   td          ! tidal harmonics data
289      INTEGER, OPTIONAL, INTENT(in   ) ::   kit         ! Barotropic timestep counter (for timesplitting option)
290      INTEGER, OPTIONAL, INTENT(in   ) ::   kt_offset   ! time offset in units of timesteps. NB. if kit
291      !                                                 ! is present then units = subcycle timesteps.
292      !                                                 ! kt_offset = 0  => get data at "now"    time level
293      !                                                 ! kt_offset = -1 => get data at "before" time level
294      !                                                 ! kt_offset = +1 => get data at "after"  time level
295      !                                                 ! etc.
296      !
297      INTEGER  ::   itide, igrd, ib       ! dummy loop indices
298      INTEGER  ::   time_add              ! time offset in units of timesteps
299      INTEGER, DIMENSION(3) ::   ilen0    ! length of boundary data (from OBC arrays)
300      REAL(wp) ::   z_arg, z_sarg, zflag, zramp   ! local scalars   
301      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
302      !!----------------------------------------------------------------------
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(kit) ) THEN
310        IF ( kit /= 1 ) zflag=0
311      ENDIF
312
313      IF ( (nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000) .AND. zflag==1 ) THEN
314        !
315        kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
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(kt_offset) ) THEN
330         time_add = kt_offset
331      ENDIF
332         
333      IF( PRESENT(kit) ) THEN 
334         z_arg = ((kt-kt_tide) * rdt + (kit+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   END SUBROUTINE bdytide_update
365
366
367   SUBROUTINE bdy_dta_tides( kt, kit, kt_offset )
368      !!----------------------------------------------------------------------
369      !!                 ***  SUBROUTINE bdy_dta_tides  ***
370      !!               
371      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
372      !!               
373      !!----------------------------------------------------------------------
374      INTEGER,           INTENT(in) ::   kt          ! Main timestep counter
375      INTEGER, OPTIONAL, INTENT(in) ::   kit         ! Barotropic timestep counter (for timesplitting option)
376      INTEGER, OPTIONAL, INTENT(in) ::   kt_offset   ! time offset in units of timesteps. NB. if kit
377      !                                              ! is present then units = subcycle timesteps.
378      !                                              ! kt_offset = 0  => get data at "now"    time level
379      !                                              ! kt_offset = -1 => get data at "before" time level
380      !                                              ! kt_offset = +1 => get data at "after"  time level
381      !                                              ! etc.
382      !
383      LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step
384      INTEGER  ::   itide, ib_bdy, ib, igrd   ! loop indices
385      INTEGER  ::   time_add                  ! time offset in units of timesteps
386      INTEGER, DIMENSION(jpbgrd)   ::   ilen0 
387      INTEGER, DIMENSION(1:jpbgrd) ::   nblen, nblenrim  ! short cuts
388      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist     
389      !!----------------------------------------------------------------------
390      !
391      lk_first_btstp=.TRUE.
392      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
393
394      time_add = 0
395      IF( PRESENT(kt_offset) ) THEN
396         time_add = kt_offset
397      ENDIF
398     
399      ! Absolute time from model initialization:   
400      IF( PRESENT(kit) ) THEN 
401         z_arg = ( kt + (kit+time_add-1) / REAL(nn_baro,wp) ) * rdt
402      ELSE                             
403         z_arg = ( kt + time_add ) * rdt
404      ENDIF
405
406      ! Linear ramp on tidal component at open boundaries
407      zramp = 1.
408      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.)
409
410      DO ib_bdy = 1,nb_bdy
411         !
412         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
413            !
414            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
415            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
416            !
417            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:)
418            ELSE                                   ;   ilen0(:) = nblenrim(:)
419            ENDIF     
420            !
421            ! We refresh nodal factors every day below
422            ! This should be done somewhere else
423            IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
424               !
425               kt_tide = kt - (nsec_day - 0.5_wp * rdt)/rdt
426               !
427               IF(lwp) THEN
428               WRITE(numout,*)
429               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
430               WRITE(numout,*) '~~~~~~~~~~~~~~ '
431               ENDIF
432               !
433               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
434               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
435               !
436            ENDIF
437            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
438            !
439            ! If time splitting, initialize arrays from slow varying open boundary data:
440            IF ( PRESENT(kit) ) THEN           
441               IF ( dta_bdy(ib_bdy)%lneed_ssh   ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
442               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
443               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
444            ENDIF
445            !
446            ! Update open boundary data arrays:
447            DO itide = 1, nb_harmo
448               !
449               z_sarg = (z_arg + zoff) * omega_tide(itide)
450               z_cost = zramp * COS( z_sarg )
451               z_sist = zramp * SIN( z_sarg )
452               !
453               IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN
454                  igrd=1                              ! SSH on tracer grid
455                  DO ib = 1, ilen0(igrd)
456                     dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
457                        &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
458                        &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
459                  END DO
460               ENDIF
461               !
462               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
463                  igrd=2                              ! U grid
464                  DO ib = 1, ilen0(igrd)
465                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
466                        &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
467                        &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
468                  END DO
469                  igrd=3                              ! V grid
470                  DO ib = 1, ilen0(igrd) 
471                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
472                        &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
473                        &                        tides(ib_bdy)%v(ib,itide,2)*z_sist )
474                  END DO
475               ENDIF
476            END DO             
477         END IF
478      END DO
479      !
480   END SUBROUTINE bdy_dta_tides
481
482
483   SUBROUTINE tide_init_elevation( idx, td )
484      !!----------------------------------------------------------------------
485      !!                 ***  ROUTINE tide_init_elevation  ***
486      !!----------------------------------------------------------------------
487      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
488      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
489      !
490      INTEGER ::   itide, igrd, ib       ! dummy loop indices
491      INTEGER, DIMENSION(1) ::   ilen0   ! length of boundary data (from OBC arrays)
492      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
493      !!----------------------------------------------------------------------
494      !
495      igrd=1   
496                              ! SSH on tracer grid.
497      ilen0(1) =  SIZE(td%ssh0(:,1,1))
498      !
499      ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) )
500      !
501      DO itide = 1, nb_harmo
502         DO ib = 1, ilen0(igrd)
503            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
504            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
505         END DO
506         DO ib = 1 , ilen0(igrd)
507            mod_tide(ib)=mod_tide(ib)*ftide(itide)
508            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
509         ENDDO
510         DO ib = 1 , ilen0(igrd)
511            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
512            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
513         ENDDO
514      END DO
515      !
516      DEALLOCATE( mod_tide, phi_tide )
517      !
518   END SUBROUTINE tide_init_elevation
519
520
521   SUBROUTINE tide_init_velocities( idx, td )
522      !!----------------------------------------------------------------------
523      !!                 ***  ROUTINE tide_init_elevation  ***
524      !!----------------------------------------------------------------------
525      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
526      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
527      !
528      INTEGER ::   itide, igrd, ib       ! dummy loop indices
529      INTEGER, DIMENSION(3) ::   ilen0   ! length of boundary data (from OBC arrays)
530      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
531      !!----------------------------------------------------------------------
532      !
533      ilen0(2) =  SIZE(td%u0(:,1,1))
534      ilen0(3) =  SIZE(td%v0(:,1,1))
535      !
536      igrd=2                                 ! U grid.
537      !
538      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
539      !
540      DO itide = 1, nb_harmo
541         DO ib = 1, ilen0(igrd)
542            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
543            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
544         END DO
545         DO ib = 1, ilen0(igrd)
546            mod_tide(ib)=mod_tide(ib)*ftide(itide)
547            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
548         ENDDO
549         DO ib = 1, ilen0(igrd)
550            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
551            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
552         ENDDO
553      END DO
554      !
555      DEALLOCATE( mod_tide , phi_tide )
556      !
557      igrd=3                                 ! V grid.
558      !
559      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
560
561      DO itide = 1, nb_harmo
562         DO ib = 1, ilen0(igrd)
563            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
564            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
565         END DO
566         DO ib = 1, ilen0(igrd)
567            mod_tide(ib)=mod_tide(ib)*ftide(itide)
568            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
569         ENDDO
570         DO ib = 1, ilen0(igrd)
571            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
572            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
573         ENDDO
574      END DO
575      !
576      DEALLOCATE( mod_tide, phi_tide )
577      !
578  END SUBROUTINE tide_init_velocities
579
580   !!======================================================================
581END MODULE bdytides
582
Note: See TracBrowser for help on using the repository browser.