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_r11943_MERGE_2019/src/OCE/BDY – NEMO

source: NEMO/branches/2019/dev_r11943_MERGE_2019/src/OCE/BDY/bdytides.F90 @ 12340

Last change on this file since 12340 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.
  • Property svn:keywords set to Id
File size: 23.0 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   bdy_dta_tides    ! routine called in dyn_spg_ts
33
34   TYPE, PUBLIC ::   TIDES_DATA     !: Storage for external tidal harmonics data
35      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh0     !: Tidal constituents : SSH0   (read in file)
36      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u0, v0   !: Tidal constituents : U0, V0 (read in file)
37      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   ssh      !: Tidal constituents : SSH    (after nodal cor.)
38      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   u , v    !: Tidal constituents : U , V  (after nodal cor.)
39   END TYPE TIDES_DATA
40
41!$AGRIF_DO_NOT_TREAT
42   TYPE(TIDES_DATA), PUBLIC, DIMENSION(jp_bdy), TARGET :: tides  !: External tidal harmonics data
43!$AGRIF_END_DO_NOT_TREAT
44   TYPE(OBC_DATA)  , PUBLIC, DIMENSION(jp_bdy) :: dta_bdy_s  !: bdy external data (slow component)
45
46   INTEGER ::   kt_tide
47
48   !! * Substitutions
49#  include "do_loop_substitute.h90"
50   !!----------------------------------------------------------------------
51   !! NEMO/OCE 4.0 , NEMO Consortium (2018)
52   !! $Id$
53   !! Software governed by the CeCILL license (see ./LICENSE)
54   !!----------------------------------------------------------------------
55CONTAINS
56
57   SUBROUTINE bdytide_init
58      !!----------------------------------------------------------------------
59      !!                    ***  SUBROUTINE bdytide_init  ***
60      !!                     
61      !! ** Purpose : - Read in namelist for tides and initialise external
62      !!                tidal harmonics data
63      !!
64      !!----------------------------------------------------------------------
65      !! namelist variables
66      !!-------------------
67      CHARACTER(len=80)                         ::   filtide             !: Filename root for tidal input files
68      LOGICAL                                   ::   ln_bdytide_2ddta    !: If true, read 2d harmonic 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                                   ::   ios                 ! Local integer output status for namelist read
75      INTEGER                                   ::   nbdy_rdstart, nbdy_loc
76      CHARACTER(LEN=50)                         ::   cerrmsg             !: error string
77      CHARACTER(len=80)                         ::   clfile              !: full file name for tidal input file
78      REAL(wp),ALLOCATABLE, DIMENSION(:,:,:)    ::   dta_read            !: work space to read in tidal harmonics data
79      REAL(wp),ALLOCATABLE, DIMENSION(:,:)      ::   ztr, zti            !:  "     "    "   "   "   "        "      "
80      !!
81      TYPE(TIDES_DATA),  POINTER                ::   td                  !: local short cut   
82      !!
83      NAMELIST/nambdy_tide/filtide, ln_bdytide_2ddta
84      !!----------------------------------------------------------------------
85      !
86      IF(lwp) WRITE(numout,*)
87      IF(lwp) WRITE(numout,*) 'bdytide_init : initialization of tidal harmonic forcing at open boundaries'
88      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~'
89
90
91      nbdy_rdstart = 1
92      DO ib_bdy = 1, nb_bdy
93         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
94            !
95            td => tides(ib_bdy)
96
97            ! Namelist nambdy_tide : tidal harmonic forcing at open boundaries
98            filtide(:) = ''
99
100            READ  ( numnam_ref, nambdy_tide, IOSTAT = ios, ERR = 901)
101901         IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_tide in reference namelist' )
102            !
103            ! Need to support possibility of reading more than one
104            ! nambdy_tide from the namelist_cfg internal file.
105            ! Do this by finding the ib_bdy'th occurence of nambdy_tide in the
106            ! character buffer as the starting point.
107            !
108            nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_tide' )
109            IF( nbdy_loc .GT. 0 ) THEN
110               nbdy_rdstart = nbdy_rdstart + nbdy_loc
111            ELSE
112               WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',ib_bdy,' of nambdy_tide not found'
113               ios = -1
114               CALL ctl_nam ( ios , cerrmsg )
115            ENDIF
116            READ  ( numnam_cfg( MAX( 1, nbdy_rdstart - 2 ): ), nambdy_tide, IOSTAT = ios, ERR = 902)
117902         IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy_tide in configuration namelist' )
118            IF(lwm) WRITE ( numond, nambdy_tide )
119            !                                               ! Parameter control and print
120            IF(lwp) WRITE(numout,*) '  '
121            IF(lwp) WRITE(numout,*) '          Namelist nambdy_tide : tidal harmonic forcing at open boundaries'
122            IF(lwp) WRITE(numout,*) '             read tidal data in 2d files: ', ln_bdytide_2ddta
123            IF(lwp) WRITE(numout,*) '             Number of tidal components to read: ', nb_harmo
124            IF(lwp) THEN
125                    WRITE(numout,*) '             Tidal components: ' 
126               DO itide = 1, nb_harmo
127                  WRITE(numout,*)  '                 ', tide_harmonics(itide)%cname_tide 
128               END DO
129            ENDIF
130            IF(lwp) WRITE(numout,*) ' '
131
132            ! Allocate space for tidal harmonics data - get size from OBC data arrays
133            ! -----------------------------------------------------------------------
134
135            ! JC: If FRS scheme is used, we assume that tidal is needed over the whole
136            ! relaxation area     
137            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = idx_bdy(ib_bdy)%nblen   (:)
138            ELSE                                   ;   ilen0(:) = idx_bdy(ib_bdy)%nblenrim(:)
139            ENDIF
140
141            ALLOCATE( td%ssh0( ilen0(1), nb_harmo, 2 ) )
142            ALLOCATE( td%ssh ( ilen0(1), nb_harmo, 2 ) )
143
144            ALLOCATE( td%u0( ilen0(2), nb_harmo, 2 ) )
145            ALLOCATE( td%u ( ilen0(2), nb_harmo, 2 ) )
146
147            ALLOCATE( td%v0( ilen0(3), nb_harmo, 2 ) )
148            ALLOCATE( td%v ( ilen0(3), nb_harmo, 2 ) )
149
150            td%ssh0(:,:,:) = 0._wp
151            td%ssh (:,:,:) = 0._wp
152            td%u0  (:,:,:) = 0._wp
153            td%u   (:,:,:) = 0._wp
154            td%v0  (:,:,:) = 0._wp
155            td%v   (:,:,:) = 0._wp
156
157            IF( ln_bdytide_2ddta ) THEN
158               ! It is assumed that each data file contains all complex harmonic amplitudes
159               ! given on the global domain (ie global, jpiglo x jpjglo)
160               !
161               ALLOCATE( zti(jpi,jpj), ztr(jpi,jpj) )
162               !
163               ! SSH fields
164               clfile = TRIM(filtide)//'_grid_T.nc'
165               CALL iom_open( clfile , inum ) 
166               igrd = 1                       ! Everything is at T-points here
167               DO itide = 1, nb_harmo
168                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z1', ztr(:,:) )
169                  CALL iom_get( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_z2', zti(:,:) ) 
170                  DO ib = 1, ilen0(igrd)
171                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
172                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
173                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
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_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_u1', ztr(:,:) )
186                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(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                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
191                     td%u0(ib,itide,1) = ztr(ii,ij)
192                     td%u0(ib,itide,2) = zti(ii,ij)
193                  END DO
194               END DO
195               CALL iom_close( inum )
196               !
197               ! V fields
198               clfile = TRIM(filtide)//'_grid_V.nc'
199               CALL iom_open( clfile , inum ) 
200               igrd = 3                       ! Everything is at V-points here
201               DO itide = 1, nb_harmo
202                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v1', ztr(:,:) )
203                  CALL iom_get  ( inum, jpdom_autoglo, TRIM(tide_harmonics(itide)%cname_tide)//'_v2', zti(:,:) )
204                  DO ib = 1, ilen0(igrd)
205                     ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
206                     ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
207                     IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE   ! to remove?
208                     td%v0(ib,itide,1) = ztr(ii,ij)
209                     td%v0(ib,itide,2) = zti(ii,ij)
210                  END DO
211               END DO 
212               CALL iom_close( inum )
213               !
214               DEALLOCATE( ztr, zti ) 
215               !
216            ELSE           
217               !
218               ! Read tidal data only on bdy segments
219               !
220               ALLOCATE( dta_read( MAXVAL(ilen0(1:3)), 1, 1 ) )
221               !
222               ! Open files and read in tidal forcing data
223               ! -----------------------------------------
224
225               DO itide = 1, nb_harmo
226                  !                                                              ! SSH fields
227                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_T.nc'
228                  CALL iom_open( clfile, inum )
229                  CALL fld_map( inum, 'z1' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
230                  td%ssh0(:,itide,1) = dta_read(1:ilen0(1),1,1)
231                  CALL fld_map( inum, 'z2' , dta_read(1:ilen0(1),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,1) )
232                  td%ssh0(:,itide,2) = dta_read(1:ilen0(1),1,1)
233                  CALL iom_close( inum )
234                  !                                                              ! U fields
235                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_U.nc'
236                  CALL iom_open( clfile, inum )
237                  CALL fld_map( inum, 'u1' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
238                  td%u0(:,itide,1) = dta_read(1:ilen0(2),1,1)
239                  CALL fld_map( inum, 'u2' , dta_read(1:ilen0(2),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,2) )
240                  td%u0(:,itide,2) = dta_read(1:ilen0(2),1,1)
241                  CALL iom_close( inum )
242                  !                                                              ! V fields
243                  clfile = TRIM(filtide)//TRIM(tide_harmonics(itide)%cname_tide)//'_grid_V.nc'
244                  CALL iom_open( clfile, inum )
245                  CALL fld_map( inum, 'v1' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
246                  td%v0(:,itide,1) = dta_read(1:ilen0(3),1,1)
247                  CALL fld_map( inum, 'v2' , dta_read(1:ilen0(3),1:1,1:1) , 1, idx_bdy(ib_bdy)%nbmap(:,3) )
248                  td%v0(:,itide,2) = dta_read(1:ilen0(3),1,1)
249                  CALL iom_close( inum )
250                  !
251               END DO ! end loop on tidal components
252               !
253               DEALLOCATE( dta_read )
254               !
255            ENDIF ! ln_bdytide_2ddta=.true.
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 bdy_dta_tides( kt, kit, pt_offset )
274      !!----------------------------------------------------------------------
275      !!                 ***  SUBROUTINE bdy_dta_tides  ***
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      INTEGER, OPTIONAL, INTENT(in) ::   kit         ! Barotropic timestep counter (for timesplitting option)
282      REAL(wp),OPTIONAL, INTENT(in) ::   pt_offset   ! time offset in units of timesteps
283      !
284      LOGICAL  ::   lk_first_btstp            ! =.TRUE. if time splitting and first barotropic step
285      INTEGER  ::   itide, ib_bdy, ib, igrd   ! loop indices
286      INTEGER, DIMENSION(jpbgrd)   ::   ilen0 
287      INTEGER, DIMENSION(1:jpbgrd) ::   nblen, nblenrim  ! short cuts
288      REAL(wp) ::   z_arg, z_sarg, zramp, zoff, z_cost, z_sist, zt_offset   
289      !!----------------------------------------------------------------------
290      !
291      lk_first_btstp=.TRUE.
292      IF ( PRESENT(kit).AND.( kit /= 1 ) ) THEN ; lk_first_btstp=.FALSE. ; ENDIF
293
294      zt_offset = 0._wp
295      IF( PRESENT(pt_offset) )   zt_offset = pt_offset
296     
297      ! Absolute time from model initialization:   
298      IF( PRESENT(kit) ) THEN 
299         z_arg = ( REAL(kt, wp) + ( REAL(kit, wp) + zt_offset - 1. ) / REAL(nn_baro, wp) ) * rdt
300      ELSE                             
301         z_arg = ( REAL(kt, wp) + zt_offset ) * rdt
302      ENDIF
303
304      ! Linear ramp on tidal component at open boundaries
305      zramp = 1.
306      IF (ln_tide_ramp) zramp = MIN(MAX( (z_arg - REAL(nit000,wp)*rdt)/(rn_tide_ramp_dt*rday),0.),1.)
307
308      DO ib_bdy = 1,nb_bdy
309         !
310         IF( nn_dyn2d_dta(ib_bdy) >= 2 ) THEN
311            !
312            nblen(1:jpbgrd) = idx_bdy(ib_bdy)%nblen(1:jpbgrd)
313            nblenrim(1:jpbgrd) = idx_bdy(ib_bdy)%nblenrim(1:jpbgrd)
314            !
315            IF( cn_dyn2d(ib_bdy) == 'frs' ) THEN   ;   ilen0(:) = nblen   (:)
316            ELSE                                   ;   ilen0(:) = nblenrim(:)
317            ENDIF     
318            !
319            ! We refresh nodal factors every day below
320            ! This should be done somewhere else
321            IF ( ( nsec_day == NINT(0.5_wp * rdt) .OR. kt==nit000 ) .AND. lk_first_btstp ) THEN
322               !
323               kt_tide = kt - NINT((REAL(nsec_day,wp) - 0.5_wp * rdt)/rdt)
324               !
325               IF(lwp) THEN
326               WRITE(numout,*)
327               WRITE(numout,*) 'bdy_tide_dta : Refresh nodal factors for tidal open bdy data at kt=',kt
328               WRITE(numout,*) '~~~~~~~~~~~~~~ '
329               ENDIF
330               !
331               CALL tide_init_elevation ( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
332               CALL tide_init_velocities( idx=idx_bdy(ib_bdy), td=tides(ib_bdy) )
333               !
334            ENDIF
335            zoff = REAL(-kt_tide,wp) * rdt ! time offset relative to nodal factor computation time
336            !
337            ! If time splitting, initialize arrays from slow varying open boundary data:
338            IF ( PRESENT(kit) ) THEN           
339               IF ( dta_bdy(ib_bdy)%lneed_ssh   ) dta_bdy(ib_bdy)%ssh(1:ilen0(1)) = dta_bdy_s(ib_bdy)%ssh(1:ilen0(1))
340               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%u2d(1:ilen0(2)) = dta_bdy_s(ib_bdy)%u2d(1:ilen0(2))
341               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) dta_bdy(ib_bdy)%v2d(1:ilen0(3)) = dta_bdy_s(ib_bdy)%v2d(1:ilen0(3))
342            ENDIF
343            !
344            ! Update open boundary data arrays:
345            DO itide = 1, nb_harmo
346               !
347               z_sarg = (z_arg + zoff) * tide_harmonics(itide)%omega
348               z_cost = zramp * COS( z_sarg )
349               z_sist = zramp * SIN( z_sarg )
350               !
351               IF ( dta_bdy(ib_bdy)%lneed_ssh ) THEN
352                  igrd=1                              ! SSH on tracer grid
353                  DO ib = 1, ilen0(igrd)
354                     dta_bdy(ib_bdy)%ssh(ib) = dta_bdy(ib_bdy)%ssh(ib) + &
355                        &                      ( tides(ib_bdy)%ssh(ib,itide,1)*z_cost + &
356                        &                        tides(ib_bdy)%ssh(ib,itide,2)*z_sist )
357                  END DO
358               ENDIF
359               !
360               IF ( dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
361                  igrd=2                              ! U grid
362                  DO ib = 1, ilen0(igrd)
363                     dta_bdy(ib_bdy)%u2d(ib) = dta_bdy(ib_bdy)%u2d(ib) + &
364                        &                      ( tides(ib_bdy)%u(ib,itide,1)*z_cost + &
365                        &                        tides(ib_bdy)%u(ib,itide,2)*z_sist )
366                  END DO
367                  igrd=3                              ! V grid
368                  DO ib = 1, ilen0(igrd) 
369                     dta_bdy(ib_bdy)%v2d(ib) = dta_bdy(ib_bdy)%v2d(ib) + &
370                        &                      ( tides(ib_bdy)%v(ib,itide,1)*z_cost + &
371                        &                        tides(ib_bdy)%v(ib,itide,2)*z_sist )
372                  END DO
373               ENDIF
374            END DO             
375         END IF
376      END DO
377      !
378   END SUBROUTINE bdy_dta_tides
379
380
381   SUBROUTINE tide_init_elevation( idx, td )
382      !!----------------------------------------------------------------------
383      !!                 ***  ROUTINE tide_init_elevation  ***
384      !!----------------------------------------------------------------------
385      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
386      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
387      !
388      INTEGER ::   itide, igrd, ib       ! dummy loop indices
389      INTEGER, DIMENSION(1) ::   ilen0   ! length of boundary data (from OBC arrays)
390      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
391      !!----------------------------------------------------------------------
392      !
393      igrd=1   
394                              ! SSH on tracer grid.
395      ilen0(1) =  SIZE(td%ssh0(:,1,1))
396      !
397      ALLOCATE( mod_tide(ilen0(igrd)), phi_tide(ilen0(igrd)) )
398      !
399      DO itide = 1, nb_harmo
400         DO ib = 1, ilen0(igrd)
401            mod_tide(ib)=SQRT(td%ssh0(ib,itide,1)**2.+td%ssh0(ib,itide,2)**2.)
402            phi_tide(ib)=ATAN2(-td%ssh0(ib,itide,2),td%ssh0(ib,itide,1))
403         END DO
404         DO ib = 1 , ilen0(igrd)
405            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
406            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0+tide_harmonics(itide)%u
407         ENDDO
408         DO ib = 1 , ilen0(igrd)
409            td%ssh(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
410            td%ssh(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
411         ENDDO
412      END DO
413      !
414      DEALLOCATE( mod_tide, phi_tide )
415      !
416   END SUBROUTINE tide_init_elevation
417
418
419   SUBROUTINE tide_init_velocities( idx, td )
420      !!----------------------------------------------------------------------
421      !!                 ***  ROUTINE tide_init_elevation  ***
422      !!----------------------------------------------------------------------
423      TYPE(OBC_INDEX) , INTENT(in   ) ::   idx   ! OBC indices
424      TYPE(TIDES_DATA), INTENT(inout) ::   td    ! tidal harmonics data
425      !
426      INTEGER ::   itide, igrd, ib       ! dummy loop indices
427      INTEGER, DIMENSION(3) ::   ilen0   ! length of boundary data (from OBC arrays)
428      REAL(wp),ALLOCATABLE, DIMENSION(:) ::   mod_tide, phi_tide
429      !!----------------------------------------------------------------------
430      !
431      ilen0(2) =  SIZE(td%u0(:,1,1))
432      ilen0(3) =  SIZE(td%v0(:,1,1))
433      !
434      igrd=2                                 ! U grid.
435      !
436      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
437      !
438      DO itide = 1, nb_harmo
439         DO ib = 1, ilen0(igrd)
440            mod_tide(ib)=SQRT(td%u0(ib,itide,1)**2.+td%u0(ib,itide,2)**2.)
441            phi_tide(ib)=ATAN2(-td%u0(ib,itide,2),td%u0(ib,itide,1))
442         END DO
443         DO ib = 1, ilen0(igrd)
444            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
445            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
446         ENDDO
447         DO ib = 1, ilen0(igrd)
448            td%u(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
449            td%u(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
450         ENDDO
451      END DO
452      !
453      DEALLOCATE( mod_tide , phi_tide )
454      !
455      igrd=3                                 ! V grid.
456      !
457      ALLOCATE( mod_tide(ilen0(igrd)) , phi_tide(ilen0(igrd)) )
458
459      DO itide = 1, nb_harmo
460         DO ib = 1, ilen0(igrd)
461            mod_tide(ib)=SQRT(td%v0(ib,itide,1)**2.+td%v0(ib,itide,2)**2.)
462            phi_tide(ib)=ATAN2(-td%v0(ib,itide,2),td%v0(ib,itide,1))
463         END DO
464         DO ib = 1, ilen0(igrd)
465            mod_tide(ib)=mod_tide(ib)*tide_harmonics(itide)%f
466            phi_tide(ib)=phi_tide(ib)+tide_harmonics(itide)%v0 + tide_harmonics(itide)%u
467         ENDDO
468         DO ib = 1, ilen0(igrd)
469            td%v(ib,itide,1)= mod_tide(ib)*COS(phi_tide(ib))
470            td%v(ib,itide,2)=-mod_tide(ib)*SIN(phi_tide(ib))
471         ENDDO
472      END DO
473      !
474      DEALLOCATE( mod_tide, phi_tide )
475      !
476  END SUBROUTINE tide_init_velocities
477
478   !!======================================================================
479END MODULE bdytides
480
Note: See TracBrowser for help on using the repository browser.