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/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/2014/dev_r4704_NOC5_MPP_BDY_UPDATE/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 6225

Last change on this file since 6225 was 6225, checked in by jamesharle, 8 years ago

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