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/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/2016/dev_merge_2016/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 7412

Last change on this file since 7412 was 7412, checked in by lovato, 8 years ago

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