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

source: trunk/NEMOGCM/NEMO/OPA_SRC/BDY/bdytides.F90 @ 4673

Last change on this file since 4673 was 4624, checked in by acc, 10 years ago

#1305. Fix slow start-up problems on some systems by introducing and using lwm logical to restrict output of merged namelists to the first (or only) processor. lwm is true only on the first processor regardless of ln_ctl. Small changes to all flavours of nemogcm.F90 are also required to write namctl and namcfg after the call to mynode which now opens output.namelist.dyn and writes nammpp.
  • 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, ONLY: fld_map
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      !!
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 cpt name    -     Phase speed (deg/hr)'           
128               DO itide = 1, nb_harmo
129                  WRITE(numout,*)  '             ', Wave(ntide(itide))%cname_tide, omega_tide(itide)   
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               ! Open files and read in tidal forcing data
224               ! -----------------------------------------
225
226               DO itide = 1, nb_harmo
227                  !                                                              ! SSH fields
228                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_T.nc'
229                  CALL iom_open( clfile, inum )
230                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
231                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
232                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
233                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
234                  CALL iom_close( inum )
235                  !                                                              ! U fields
236                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_U.nc'
237                  CALL iom_open( clfile, inum )
238                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
239                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
240                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
241                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
242                  CALL iom_close( inum )
243                  !                                                              ! V fields
244                  clfile = TRIM(filtide)//TRIM(Wave(ntide(itide))%cname_tide)//'_grid_V.nc'
245                  CALL iom_open( clfile, inum )
246                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
247                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
248                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
249                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
250                  CALL iom_close( inum )
251                  !
252               END DO ! end loop on tidal components
253               !
254               DEALLOCATE( dta_read )
255            ENDIF ! ln_bdytide_2ddta=.true.
256            !
257            IF ( ln_bdytide_conj ) THEN ! assume complex conjugate in data files
258               td%ssh0(:,:,2) = - td%ssh0(:,:,2)
259               td%u0  (:,:,2) = - td%u0  (:,:,2)
260               td%v0  (:,:,2) = - td%v0  (:,:,2)
261            ENDIF
262            !
263            IF ( lk_dynspg_ts ) THEN ! Allocate arrays to save slowly varying boundary data during
264                                     ! time splitting integration
265               ALLOCATE( dta_bdy_s(ib_bdy)%ssh ( ilen0(1) ) )
266               ALLOCATE( dta_bdy_s(ib_bdy)%u2d ( ilen0(2) ) )
267               ALLOCATE( dta_bdy_s(ib_bdy)%v2d ( ilen0(3) ) )
268               dta_bdy_s(ib_bdy)%ssh(:) = 0.e0
269               dta_bdy_s(ib_bdy)%u2d(:) = 0.e0
270               dta_bdy_s(ib_bdy)%v2d(:) = 0.e0
271            ENDIF
272            !
273         ENDIF ! nn_dyn2d_dta(ib_bdy) .ge. 2
274         !
275      END DO ! loop on ib_bdy
276
277      IF( nn_timing == 1 ) CALL timing_stop('bdytide_init')
278
279   END SUBROUTINE bdytide_init
280
281   SUBROUTINE bdytide_update ( kt, idx, dta, td, jit, time_offset )
282      !!----------------------------------------------------------------------
283      !!                 ***  SUBROUTINE bdytide_update  ***
284      !!               
285      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
286      !!               
287      !!----------------------------------------------------------------------
288      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter
289      TYPE(OBC_INDEX), INTENT( in )    ::   idx         ! OBC indices
290      TYPE(OBC_DATA),  INTENT(inout)   ::   dta         ! OBC external data
291      TYPE(TIDES_DATA),INTENT( inout ) ::   td          ! tidal harmonics data
292      INTEGER,INTENT(in),OPTIONAL      ::   jit         ! Barotropic timestep counter (for timesplitting option)
293      INTEGER,INTENT( in ), OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if jit
294                                                        ! is present then units = subcycle timesteps.
295                                                        ! time_offset = 0  => get data at "now"    time level
296                                                        ! time_offset = -1 => get data at "before" time level
297                                                        ! time_offset = +1 => get data at "after"  time level
298                                                        ! etc.
299      !!
300      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays)
301      INTEGER                          :: itide, igrd, ib   ! dummy loop indices
302      INTEGER                          :: time_add          ! time offset in units of timesteps
303      REAL(wp)                         :: z_arg, z_sarg, zflag, zramp     
304      REAL(wp), DIMENSION(jpmax_harmo) :: z_sist, z_cost
305      !!----------------------------------------------------------------------
306
307      IF( nn_timing == 1 ) CALL timing_start('bdytide_update')
308
309      ilen0(1) =  SIZE(td%ssh(:,1,1))
310      ilen0(2) =  SIZE(td%u(:,1,1))
311      ilen0(3) =  SIZE(td%v(:,1,1))
312
313      zflag=1
314      IF ( PRESENT(jit) ) THEN
315        IF ( jit /= 1 ) zflag=0
316      ENDIF
317
318      IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. zflag==1 ) THEN
319        !
320        kt_tide = kt
321        !
322        IF(lwp) THEN
323           WRITE(numout,*)
324           WRITE(numout,*) 'bdytide_update : (re)Initialization of the tidal bdy forcing at kt=',kt
325           WRITE(numout,*) '~~~~~~~~~~~~~~ '
326        ENDIF
327        !
328        CALL tide_init_elevation ( idx, td )
329        CALL tide_init_velocities( idx, td )
330        !
331      ENDIF
332
333      time_add = 0
334      IF( PRESENT(time_offset) ) THEN
335         time_add = time_offset
336      ENDIF
337         
338      IF( PRESENT(jit) ) THEN 
339         z_arg = ((kt-kt_tide) * rdt + (jit+0.5_wp*(time_add-1)) * rdt / REAL(nn_baro,wp) )
340      ELSE                             
341         z_arg = ((kt-kt_tide)+time_add) * rdt
342      ENDIF
343
344      ! Linear ramp on tidal component at open boundaries
345      zramp = 1._wp
346      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg + (kt_tide-nit000)*rdt)/(rdttideramp*rday),0._wp),1._wp)
347
348      DO itide = 1, nb_harmo
349         z_sarg = z_arg * omega_tide(itide)
350         z_cost(itide) = COS( z_sarg )
351         z_sist(itide) = SIN( z_sarg )
352      END DO
353
354      DO itide = 1, nb_harmo
355         igrd=1                              ! SSH on tracer grid
356         DO ib = 1, ilen0(igrd)
357            dta%ssh(ib) = dta%ssh(ib) + zramp*(td%ssh(ib,itide,1)*z_cost(itide) + td%ssh(ib,itide,2)*z_sist(itide))
358         END DO
359         igrd=2                              ! U grid
360         DO ib = 1, ilen0(igrd)
361            dta%u2d(ib) = dta%u2d(ib) + zramp*(td%u  (ib,itide,1)*z_cost(itide) + td%u  (ib,itide,2)*z_sist(itide))
362         END DO
363         igrd=3                              ! V grid
364         DO ib = 1, ilen0(igrd) 
365            dta%v2d(ib) = dta%v2d(ib) + zramp*(td%v  (ib,itide,1)*z_cost(itide) + td%v  (ib,itide,2)*z_sist(itide))
366         END DO
367      END DO
368      !
369      IF( nn_timing == 1 ) CALL timing_stop('bdytide_update')
370      !
371   END SUBROUTINE bdytide_update
372
373   SUBROUTINE bdy_dta_tides( kt, kit, time_offset )
374      !!----------------------------------------------------------------------
375      !!                 ***  SUBROUTINE bdy_dta_tides  ***
376      !!               
377      !! ** Purpose : - Add tidal forcing to ssh, u2d and v2d OBC data arrays.
378      !!               
379      !!----------------------------------------------------------------------
380      INTEGER, INTENT( in )            ::   kt          ! Main timestep counter
381      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Barotropic timestep counter (for timesplitting option)
382      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! time offset in units of timesteps. NB. if kit
383                                                        ! is present then units = subcycle timesteps.
384                                                        ! time_offset = 0  => get data at "now"    time level
385                                                        ! time_offset = -1 => get data at "before" time level
386                                                        ! time_offset = +1 => get data at "after"  time level
387                                                        ! etc.
388      !!
389      LOGICAL  :: lk_first_btstp  ! =.TRUE. if time splitting and first barotropic step
390      INTEGER,          DIMENSION(jpbgrd) :: ilen0 
391      INTEGER, DIMENSION(1:jpbgrd) :: nblen, nblenrim  ! short cuts
392      INTEGER  :: itide, ib_bdy, ib, igrd                     ! loop indices
393      INTEGER  :: time_add                                    ! time offset in units of timesteps
394      REAL(wp) :: z_arg, z_sarg, zramp, zoff, z_cost, z_sist     
395      !!----------------------------------------------------------------------
396
397      IF( nn_timing == 1 ) CALL timing_start('bdy_dta_tides')
398
399      lk_first_btstp=.TRUE.
400      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
401
402      time_add = 0
403      IF( PRESENT(time_offset) ) THEN
404         time_add = time_offset
405      ENDIF
406     
407      ! Absolute time from model initialization:   
408      IF( PRESENT(kit) ) THEN 
409         z_arg = ( kt + (kit+0.5_wp*(time_add-1)) / REAL(nn_baro,wp) ) * rdt
410      ELSE                             
411         z_arg = ( kt + time_add ) * rdt
412      ENDIF
413
414      ! Linear ramp on tidal component at open boundaries
415      zramp = 1.
416      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - nit000*rdt)/(rdttideramp*rday),0.),1.)
417
418      DO ib_bdy = 1,nb_bdy
419
420         ! line below should be simplified (runoff case)
421!! CHANUT: TO BE SORTED OUT
422!!         IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(nn_tra(ib_bdy).NE.4)) THEN
423         IF ( nn_dyn2d_dta(ib_bdy) .ge. 2 ) THEN
424
425            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
426            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
427
428            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN
429               ilen0(:)=nblen(:)
430            ELSE
431               ilen0(:)=nblenrim(:)
432            ENDIF     
433
434            ! We refresh nodal factors every day below
435            ! This should be done somewhere else
436            IF ( nsec_day == NINT(0.5_wp * rdttra(1)) .AND. lk_first_btstp ) THEN
437               !
438               kt_tide = kt               
439               !
440               IF(lwp) THEN
441               WRITE(numout,*)
442               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
443               WRITE(numout,*) '~~~~~~~~~~~~~~ '
444               ENDIF
445               !
446               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
447               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
448               !
449            ENDIF
450            zoff = -kt_tide * rdt ! time offset relative to nodal factor computation time
451            !
452            ! If time splitting, save data at first barotropic iteration
453            IF ( PRESENT(kit) ) THEN
454               IF ( lk_first_btstp ) THEN ! Save slow varying open boundary data:
455                  dta_bdy_s(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy(ib_bdy)%ssh(1:ilen0(1))
456                  dta_bdy_s(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy(ib_bdy)%u2d(1:ilen0(2))
457                  dta_bdy_s(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy(ib_bdy)%v2d(1:ilen0(3))
458
459               ELSE ! Initialize arrays from slow varying open boundary data:           
460                  dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
461                  dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
462                  dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
463               ENDIF
464            ENDIF
465            !
466            ! Update open boundary data arrays:
467            DO itide = 1, nb_harmo
468               !
469               z_sarg = (z_arg + zoff) * omega_tide(itide)
470               z_cost = zramp * COS( z_sarg )
471               z_sist = zramp * SIN( z_sarg )
472               !
473               igrd=1                              ! SSH on tracer grid
474               DO ib = 1, ilen0(igrd)
475                  dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
476                     &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
477                     &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
478               END DO
479               !
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               !
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            END DO
494         END IF
495      END DO
496      !
497      IF( nn_timing == 1 ) CALL timing_stop('bdy_dta_tides')
498      !
499   END SUBROUTINE bdy_dta_tides
500
501   SUBROUTINE tide_init_elevation( idx, td )
502      !!----------------------------------------------------------------------
503      !!                 ***  ROUTINE tide_init_elevation  ***
504      !!----------------------------------------------------------------------
505      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
506      TYPE(TIDES_DATA),INTENT( inout )   ::   td      ! tidal harmonics data
507      !! * Local declarations
508      INTEGER, DIMENSION(1)            ::   ilen0       !: length of boundary data (from OBC arrays)
509      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
510      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
511
512      igrd=1   
513                              ! SSH on tracer grid.
514   
515      ilen0(1) =  SIZE(td%ssh0(:,1,1))
516
517      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
518
519      DO itide = 1, nb_harmo
520         DO ib = 1, ilen0(igrd)
521            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
522            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
523         END DO
524         DO ib = 1 , ilen0(igrd)
525            mod_tide(ib)=mod_tide(ib)*ftide(itide)
526            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
527         ENDDO
528         DO ib = 1 , ilen0(igrd)
529            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
530            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
531         ENDDO
532      END DO
533
534      DEALLOCATE(mod_tide,phi_tide)
535
536   END SUBROUTINE tide_init_elevation
537
538   SUBROUTINE tide_init_velocities( idx, td )
539      !!----------------------------------------------------------------------
540      !!                 ***  ROUTINE tide_init_elevation  ***
541      !!----------------------------------------------------------------------
542      TYPE(OBC_INDEX), INTENT( in )      ::   idx     ! OBC indices
543      TYPE(TIDES_DATA),INTENT( inout )      ::   td      ! tidal harmonics data
544      !! * Local declarations
545      INTEGER, DIMENSION(3)            ::   ilen0       !: length of boundary data (from OBC arrays)
546      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
547      INTEGER                            ::   itide, igrd, ib      ! dummy loop indices
548
549      ilen0(2) =  SIZE(td%u0(:,1,1))
550      ilen0(3) =  SIZE(td%v0(:,1,1))
551
552      igrd=2                                 ! U grid.
553
554      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
555
556      DO itide = 1, nb_harmo
557         DO ib = 1, ilen0(igrd)
558            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
559            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
560         END DO
561         DO ib = 1, ilen0(igrd)
562            mod_tide(ib)=mod_tide(ib)*ftide(itide)
563            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
564         ENDDO
565         DO ib = 1, ilen0(igrd)
566            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
567            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
568         ENDDO
569      END DO
570
571      DEALLOCATE(mod_tide,phi_tide)
572
573      igrd=3                                 ! V grid.
574
575      ALLOCATE(mod_tide(ilen0(igrd)),phi_tide(ilen0(igrd)))
576
577      DO itide = 1, nb_harmo
578         DO ib = 1, ilen0(igrd)
579            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
580            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
581         END DO
582         DO ib = 1, ilen0(igrd)
583            mod_tide(ib)=mod_tide(ib)*ftide(itide)
584            phi_tide(ib)=phi_tide(ib)+v0tide(itide)+utide(itide)
585         ENDDO
586         DO ib = 1, ilen0(igrd)
587            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
588            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
589         ENDDO
590      END DO
591
592      DEALLOCATE(mod_tide,phi_tide)
593
594  END SUBROUTINE tide_init_velocities
595#else
596   !!----------------------------------------------------------------------
597   !!   Dummy module         NO Unstruct Open Boundary Conditions for tides
598   !!----------------------------------------------------------------------
599CONTAINS
600   SUBROUTINE bdytide_init             ! Empty routine
601      WRITE(*,*) 'bdytide_init: You should not have seen this print! error?'
602   END SUBROUTINE bdytide_init
603   SUBROUTINE bdytide_update( kt, jit )   ! Empty routine
604      WRITE(*,*) 'bdytide_update: You should not have seen this print! error?', kt, jit
605   END SUBROUTINE bdytide_update
606   SUBROUTINE bdy_dta_tides( kt, kit, time_offset )     ! Empty routine
607      INTEGER, INTENT( in )            ::   kt          ! Dummy argument empty routine     
608      INTEGER, INTENT( in ),OPTIONAL   ::   kit         ! Dummy argument empty routine
609      INTEGER, INTENT( in ),OPTIONAL   ::   time_offset ! Dummy argument empty routine
610      WRITE(*,*) 'bdy_dta_tides: You should not have seen this print! error?', kt, jit
611   END SUBROUTINE bdy_dta_tides
612#endif
613
614   !!======================================================================
615END MODULE bdytides
616
Note: See TracBrowser for help on using the repository browser.