source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 6079

Last change on this file since 6079 was 6079, checked in by jamesharle, 6 years ago

merge to trunk@5936

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