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 NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY – NEMO

source: NEMO/branches/2019/dev_r10742_ENHANCE-12_SimonM-Tides/src/OCE/BDY/bdytides.F90 @ 10860

Last change on this file since 10860 was 10860, checked in by smueller, 5 years ago

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