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/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/2015/dev_r5836_NOC3_vvl_by_default/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 5845

Last change on this file since 5845 was 5845, checked in by gm, 8 years ago

#1613: vvl by default: suppression of domzgr_substitute.h90

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