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.
bdyini.F90 in branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY – NEMO

source: branches/2012/dev_MERCATOR_2012_rev3555/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 3583

Last change on this file since 3583 was 3583, checked in by cbricaud, 11 years ago

add modification from dev_r3327_MERCATOR1_BDY branch in dev_MERCATOR_2012_rev3555 branch
  • Property svn:keywords set to Id
File size: 61.8 KB
Line 
1MODULE bdyini
2   !!======================================================================
3   !!                       ***  MODULE  bdyini  ***
4   !! Unstructured open boundaries : initialisation
5   !!======================================================================
6   !! History :  1.0  !  2005-01  (J. Chanut, A. Sellar)  Original code
7   !!             -   !  2007-01  (D. Storkey) Update to use IOM module
8   !!             -   !  2007-01  (D. Storkey) Tidal forcing
9   !!            3.0  !  2008-04  (NEMO team)  add in the reference version
10   !!            3.3  !  2010-09  (E.O'Dea) updates for Shelf configurations
11   !!            3.3  !  2010-09  (D.Storkey) add ice boundary conditions
12   !!            3.4  !  2011     (D. Storkey) rewrite in preparation for OBC-BDY merge
13   !!            3.4  !  2012     (J. Chanut) straight open boundary case update
14   !!----------------------------------------------------------------------
15#if defined key_bdy
16   !!----------------------------------------------------------------------
17   !!   'key_bdy'                     Unstructured Open Boundary Conditions
18   !!----------------------------------------------------------------------
19   !!   bdy_init       : Initialization of unstructured open boundaries
20   !!----------------------------------------------------------------------
21   USE timing          ! Timing
22   USE oce             ! ocean dynamics and tracers variables
23   USE dom_oce         ! ocean space and time domain
24   USE bdy_oce         ! unstructured open boundary conditions
25   USE in_out_manager  ! I/O units
26   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
27   USE lib_mpp         ! for mpp_sum 
28   USE iom             ! I/O
29   USE sbctide, ONLY: lk_tide ! Tidal forcing or not
30   USE phycst, ONLY: rday
31
32   IMPLICIT NONE
33   PRIVATE
34
35   PUBLIC   bdy_init   ! routine called in nemo_init
36
37   INTEGER, PARAMETER          :: jp_nseg = 100
38   INTEGER, PARAMETER          :: nrimmax = 20 ! maximum rimwidth in structured
39                                               ! open boundary data files
40   ! Straight open boundary segment parameters:
41   INTEGER  :: nbdysege, nbdysegw, nbdysegn, nbdysegs 
42   INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft, npckge
43   INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft, npckgw
44   INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft, npckgn
45   INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft, npckgs
46   !!----------------------------------------------------------------------
47   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
48   !! $Id$
49   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
50   !!----------------------------------------------------------------------
51CONTAINS
52   
53   SUBROUTINE bdy_init
54      !!----------------------------------------------------------------------
55      !!                 ***  ROUTINE bdy_init  ***
56      !!         
57      !! ** Purpose :   Initialization of the dynamics and tracer fields with
58      !!              unstructured open boundaries.
59      !!
60      !! ** Method  :   Read initialization arrays (mask, indices) to identify
61      !!              an unstructured open boundary
62      !!
63      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
64      !!----------------------------------------------------------------------     
65      ! namelist variables
66      !-------------------
67      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile
68      CHARACTER(LEN=1)   ::   ctypebdy
69      INTEGER :: nbdyind, nbdybeg, nbdyend
70
71      ! local variables
72      !-------------------
73      INTEGER  ::   ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices
74      INTEGER  ::   icount, icountr, ibr_max, ilen1, ibm1  ! local integers
75      INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       -
76      INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       -
77      INTEGER  ::   jpbdtau, jpbdtas                       !   -       -
78      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       -
79      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts
80      REAL   , POINTER  ::  flagu, flagv                   !    -   -
81      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
82      INTEGER, DIMENSION (2)                  ::   kdimsz
83      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays
84      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta
85      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
86      CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid
87      !!
88      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             &
89         &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, &
90         &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         & 
91         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp,              &
92#if defined key_lim2
93         &             nn_ice_lim2, nn_ice_lim2_dta,                       &
94#endif
95         &             ln_vol, nn_volctl, nn_rimwidth
96      !!
97      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
98
99      !!----------------------------------------------------------------------
100
101      IF( nn_timing == 1 ) CALL timing_start('bdy_init')
102
103      IF( bdy_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' )
104
105      IF(lwp) WRITE(numout,*)
106      IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
107      IF(lwp) WRITE(numout,*) '~~~~~~~~'
108      !
109
110      IF( jperio /= 0 )   CALL ctl_stop( 'Cyclic or symmetric,',   &
111         &                               ' and general open boundary condition are not compatible' )
112
113      cgrid= (/'t','u','v'/)
114     
115      ! -----------------------------------------
116      ! Initialise and read namelist parameters
117      ! -----------------------------------------
118
119      nb_bdy            = 0
120      ln_coords_file(:) = .false.
121      cn_coords_file(:) = ''
122      ln_mask_file      = .false.
123      cn_mask_file(:)   = ''
124      nn_dyn2d(:)       = 0
125      nn_dyn2d_dta(:)   = -1  ! uninitialised flag
126      nn_dyn3d(:)       = 0
127      nn_dyn3d_dta(:)   = -1  ! uninitialised flag
128      nn_tra(:)         = 0
129      nn_tra_dta(:)     = -1  ! uninitialised flag
130      ln_tra_dmp(:)     = .false.
131      ln_dyn3d_dmp(:)   = .false.
132      rn_time_dmp(:)    = 1.
133#if defined key_lim2
134      nn_ice_lim2(:)    = 0
135      nn_ice_lim2_dta(:)= -1  ! uninitialised flag
136#endif
137      ln_vol            = .false.
138      nn_volctl         = -1  ! uninitialised flag
139      nn_rimwidth(:)    = -1  ! uninitialised flag
140
141      REWIND( numnam )                   
142      READ  ( numnam, nambdy )
143
144      ! -----------------------------------------
145      ! Check and write out namelist parameters
146      ! -----------------------------------------
147      !                                   ! control prints
148      IF(lwp) WRITE(numout,*) '         nambdy'
149
150      IF( nb_bdy .eq. 0 ) THEN
151        IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
152      ELSE
153        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_bdy
154      ENDIF
155
156      DO ib_bdy = 1,nb_bdy
157        IF(lwp) WRITE(numout,*) ' ' 
158        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------' 
159
160        IF( ln_coords_file(ib_bdy) ) THEN
161           IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
162        ELSE
163           IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'
164        ENDIF
165        IF(lwp) WRITE(numout,*)
166
167        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
168        SELECT CASE( nn_dyn2d(ib_bdy) )                 
169          CASE(jp_none)         ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
170          CASE(jp_frs)          ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
171          CASE(jp_flather)      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition'
172          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
173        END SELECT
174        IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN
175           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !
176              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
177              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
178              CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file'
179              CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files'
180              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
181           END SELECT
182           IF (( nn_dyn2d_dta(ib_bdy) .ge. 2 ).AND.(.NOT.lk_tide)) THEN
183             CALL ctl_stop( 'You must activate key_tide to add tidal forcing at open boundaries' )
184           ENDIF
185        ENDIF
186        IF(lwp) WRITE(numout,*)
187
188        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
189        SELECT CASE( nn_dyn3d(ib_bdy) )                 
190          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
191          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
192          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
193          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
194          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' )
195        END SELECT
196        IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN
197           SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !
198              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
199              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
200              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
201           END SELECT
202        ENDIF
203
204        IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
205           IF ( nn_dyn3d(ib_bdy).EQ.0 ) THEN
206              IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'
207              ln_dyn3d_dmp(ib_bdy)=.false.
208           ELSEIF ( nn_dyn3d(ib_bdy).EQ.1 ) THEN
209              CALL ctl_stop( 'Use FRS OR relaxation' )
210           ELSE
211              IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone'
212              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
213              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
214           ENDIF
215        ELSE
216           IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities'
217        ENDIF
218        IF(lwp) WRITE(numout,*)
219
220        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
221        SELECT CASE( nn_tra(ib_bdy) )                 
222          CASE(jp_none)  ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
223          CASE(jp_frs)   ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
224          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Specified value'
225          CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      Neumann conditions'
226          CASE( 4 )      ;   IF(lwp) WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity'
227          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
228        END SELECT
229        IF( nn_tra(ib_bdy) .gt. 0 ) THEN
230           SELECT CASE( nn_tra_dta(ib_bdy) )                   !
231              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
232              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
233              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
234           END SELECT
235        ENDIF
236
237        IF ( ln_tra_dmp(ib_bdy) ) THEN
238           IF ( nn_tra(ib_bdy).EQ.0 ) THEN
239              IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
240              ln_tra_dmp(ib_bdy)=.false.
241           ELSEIF ( nn_tra(ib_bdy).EQ.1 ) THEN
242              CALL ctl_stop( 'Use FRS OR relaxation' )
243           ELSE
244              IF(lwp) WRITE(numout,*) '      + T/S relaxation zone'
245              IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
246              IF((lwp).AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
247           ENDIF
248        ELSE
249           IF(lwp) WRITE(numout,*) '      NO T/S relaxation'
250        ENDIF
251        IF(lwp) WRITE(numout,*)
252
253#if defined key_lim2
254        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
255        SELECT CASE( nn_ice_lim2(ib_bdy) )                 
256          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
257          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
258          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
259        END SELECT
260        IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN
261           SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   !
262              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
263              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
264              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' )
265           END SELECT
266        ENDIF
267        IF(lwp) WRITE(numout,*)
268#endif
269
270        IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
271        IF(lwp) WRITE(numout,*)
272
273      ENDDO
274
275     IF (nb_bdy .gt. 0) THEN
276        IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
277          IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
278          IF(lwp) WRITE(numout,*)
279          SELECT CASE ( nn_volctl )
280            CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant'
281            CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
282            CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
283          END SELECT
284          IF(lwp) WRITE(numout,*)
285        ELSE
286          IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
287          IF(lwp) WRITE(numout,*)
288        ENDIF
289     ENDIF
290
291      ! -------------------------------------------------
292      ! Initialise indices arrays for open boundaries
293      ! -------------------------------------------------
294
295      ! Work out global dimensions of boundary data
296      ! ---------------------------------------------
297      REWIND( numnam )                   
298               
299      nblendta(:,:) = 0
300      nbdysege = 0
301      nbdysegw = 0
302      nbdysegn = 0
303      nbdysegs = 0
304      icount   = 0 ! count user defined segments
305      ! Dimensions below are used to allocate arrays to read external data
306      jpbdtas = 1 ! Maximum size of boundary data (structured case)
307      jpbdtau = 1 ! Maximum size of boundary data (unstructured case)
308
309      DO ib_bdy = 1, nb_bdy
310
311         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
312 
313            icount = icount + 1
314            ! No REWIND here because may need to read more than one nambdy_index namelist.
315            READ  ( numnam, nambdy_index )
316
317            SELECT CASE ( TRIM(ctypebdy) )
318              CASE( 'N' )
319                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
320                    nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain.
321                    nbdybeg  = 2
322                    nbdyend  = jpiglo - 1
323                 ENDIF
324                 nbdysegn = nbdysegn + 1
325                 npckgn(nbdysegn) = ib_bdy ! Save bdy package number
326                 jpjnob(nbdysegn) = nbdyind
327                 jpindt(nbdysegn) = nbdybeg
328                 jpinft(nbdysegn) = nbdyend
329                 !
330              CASE( 'S' )
331                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
332                    nbdyind  = 2           ! set boundary to whole side of model domain.
333                    nbdybeg  = 2
334                    nbdyend  = jpiglo - 1
335                 ENDIF
336                 nbdysegs = nbdysegs + 1
337                 npckgs(nbdysegs) = ib_bdy ! Save bdy package number
338                 jpjsob(nbdysegs) = nbdyind
339                 jpisdt(nbdysegs) = nbdybeg
340                 jpisft(nbdysegs) = nbdyend
341                 !
342              CASE( 'E' )
343                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
344                    nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain.
345                    nbdybeg  = 2
346                    nbdyend  = jpjglo - 1
347                 ENDIF
348                 nbdysege = nbdysege + 1 
349                 npckge(nbdysege) = ib_bdy ! Save bdy package number
350                 jpieob(nbdysege) = nbdyind
351                 jpjedt(nbdysege) = nbdybeg
352                 jpjeft(nbdysege) = nbdyend
353                 !
354              CASE( 'W' )
355                 IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
356                    nbdyind  = 2           ! set boundary to whole side of model domain.
357                    nbdybeg  = 2
358                    nbdyend  = jpjglo - 1
359                 ENDIF
360                 nbdysegw = nbdysegw + 1
361                 npckgw(nbdysegw) = ib_bdy ! Save bdy package number
362                 jpiwob(nbdysegw) = nbdyind
363                 jpjwdt(nbdysegw) = nbdybeg
364                 jpjwft(nbdysegw) = nbdyend
365                 !
366              CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
367            END SELECT
368
369            ! For simplicity we assume that in case of straight bdy, arrays have the same length
370            ! (even if it is true that last tangential velocity points
371            ! are useless). This simplifies a little bit boundary data format (and agrees with format
372            ! used so far in obc package)
373
374            nblendta(1:jpbgrd,ib_bdy) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(ib_bdy)
375            jpbdtas = MAX(jpbdtas, (nbdyend - nbdybeg + 1))
376            IF (lwp.and.(nn_rimwidth(ib_bdy)>nrimmax)) &
377            & CALL ctl_stop( 'rimwidth must be lower than nrimmax' )
378
379         ELSE            ! Read size of arrays in boundary coordinates file.
380            CALL iom_open( cn_coords_file(ib_bdy), inum )
381            DO igrd = 1, jpbgrd
382               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 
383               nblendta(igrd,ib_bdy) = kdimsz(1)
384               jpbdtau = MAX(jpbdtau, kdimsz(1))
385            ENDDO
386            CALL iom_close( inum )
387
388         ENDIF
389
390      ENDDO ! ib_bdy
391
392      IF (nb_bdy>0) THEN
393         jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
394
395         ! Allocate arrays
396         !---------------
397         ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    &
398            &      nbrdta(jpbdta, jpbgrd, nb_bdy) )
399
400         ALLOCATE( dta_global(jpbdtau, 1, jpk) )
401         IF ( icount>0 ) ALLOCATE( dta_global2(jpbdtas, nrimmax, jpk) )
402         !
403      ENDIF
404
405      ! Now look for crossings in user (namelist) defined open boundary segments:
406      !--------------------------------------------------------------------------
407      IF ( icount>0 ) CALL bdy_ctl_seg
408
409      ! Calculate global boundary index arrays or read in from file
410      !------------------------------------------------------------               
411      ! 1. Read global index arrays from boundary coordinates file.
412      DO ib_bdy = 1, nb_bdy
413
414         IF( ln_coords_file(ib_bdy) ) THEN
415
416            CALL iom_open( cn_coords_file(ib_bdy), inum )
417            DO igrd = 1, jpbgrd
418               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
419               DO ii = 1,nblendta(igrd,ib_bdy)
420                  nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
421               END DO
422               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
423               DO ii = 1,nblendta(igrd,ib_bdy)
424                  nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
425               END DO
426               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
427               DO ii = 1,nblendta(igrd,ib_bdy)
428                  nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
429               END DO
430
431               ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
432               IF(lwp) WRITE(numout,*)
433               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
434               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
435               IF (ibr_max < nn_rimwidth(ib_bdy))   &
436                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
437            END DO
438            CALL iom_close( inum )
439
440         ENDIF
441
442      ENDDO     
443   
444      ! 2. Now fill indices corresponding to straight open boundary arrays:
445      ! East
446      !-----
447      DO iseg = 1, nbdysege
448         ib_bdy = npckge(iseg)
449         !
450         ! ------------ T points -------------
451         igrd=1
452         icount=0
453         DO ir = 1, nn_rimwidth(ib_bdy)
454            DO ij = jpjedt(iseg), jpjeft(iseg)
455               icount = icount + 1
456               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
457               nbjdta(icount, igrd, ib_bdy) = ij
458               nbrdta(icount, igrd, ib_bdy) = ir
459            ENDDO
460         ENDDO
461         !
462         ! ------------ U points -------------
463         igrd=2
464         icount=0
465         DO ir = 1, nn_rimwidth(ib_bdy)
466            DO ij = jpjedt(iseg), jpjeft(iseg)
467               icount = icount + 1
468               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir
469               nbjdta(icount, igrd, ib_bdy) = ij
470               nbrdta(icount, igrd, ib_bdy) = ir
471            ENDDO
472         ENDDO
473         !
474         ! ------------ V points -------------
475         igrd=3
476         icount=0
477         DO ir = 1, nn_rimwidth(ib_bdy)
478!            DO ij = jpjedt(iseg), jpjeft(iseg) - 1
479            DO ij = jpjedt(iseg), jpjeft(iseg)
480               icount = icount + 1
481               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
482               nbjdta(icount, igrd, ib_bdy) = ij
483               nbrdta(icount, igrd, ib_bdy) = ir
484            ENDDO
485            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
486            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
487         ENDDO
488      ENDDO
489      !
490      ! West
491      !-----
492      DO iseg = 1, nbdysegw
493         ib_bdy = npckgw(iseg)
494         !
495         ! ------------ T points -------------
496         igrd=1
497         icount=0
498         DO ir = 1, nn_rimwidth(ib_bdy)
499            DO ij = jpjwdt(iseg), jpjwft(iseg)
500               icount = icount + 1
501               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
502               nbjdta(icount, igrd, ib_bdy) = ij
503               nbrdta(icount, igrd, ib_bdy) = ir
504            ENDDO
505         ENDDO
506         !
507         ! ------------ U points -------------
508         igrd=2
509         icount=0
510         DO ir = 1, nn_rimwidth(ib_bdy)
511            DO ij = jpjwdt(iseg), jpjwft(iseg)
512               icount = icount + 1
513               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
514               nbjdta(icount, igrd, ib_bdy) = ij
515               nbrdta(icount, igrd, ib_bdy) = ir
516            ENDDO
517         ENDDO
518         !
519         ! ------------ V points -------------
520         igrd=3
521         icount=0
522         DO ir = 1, nn_rimwidth(ib_bdy)
523!            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
524            DO ij = jpjwdt(iseg), jpjwft(iseg)
525               icount = icount + 1
526               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
527               nbjdta(icount, igrd, ib_bdy) = ij
528               nbrdta(icount, igrd, ib_bdy) = ir
529            ENDDO
530            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
531            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
532         ENDDO
533      ENDDO
534      !
535      ! North
536      !-----
537      DO iseg = 1, nbdysegn
538         ib_bdy = npckgn(iseg)
539         !
540         ! ------------ T points -------------
541         igrd=1
542         icount=0
543         DO ir = 1, nn_rimwidth(ib_bdy)
544            DO ii = jpindt(iseg), jpinft(iseg)
545               icount = icount + 1
546               nbidta(icount, igrd, ib_bdy) = ii
547               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
548               nbrdta(icount, igrd, ib_bdy) = ir
549            ENDDO
550         ENDDO
551         !
552         ! ------------ U points -------------
553         igrd=2
554         icount=0
555         DO ir = 1, nn_rimwidth(ib_bdy)
556!            DO ii = jpindt(iseg), jpinft(iseg) - 1
557            DO ii = jpindt(iseg), jpinft(iseg)
558               icount = icount + 1
559               nbidta(icount, igrd, ib_bdy) = ii
560               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
561               nbrdta(icount, igrd, ib_bdy) = ir
562            ENDDO
563            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
564            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
565         ENDDO
566         !
567         ! ------------ V points -------------
568         igrd=3
569         icount=0
570         DO ir = 1, nn_rimwidth(ib_bdy)
571            DO ii = jpindt(iseg), jpinft(iseg)
572               icount = icount + 1
573               nbidta(icount, igrd, ib_bdy) = ii
574               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir
575               nbrdta(icount, igrd, ib_bdy) = ir
576            ENDDO
577         ENDDO
578      ENDDO
579      !
580      ! South
581      !-----
582      DO iseg = 1, nbdysegs
583         ib_bdy = npckgs(iseg)
584         !
585         ! ------------ T points -------------
586         igrd=1
587         icount=0
588         DO ir = 1, nn_rimwidth(ib_bdy)
589            DO ii = jpisdt(iseg), jpisft(iseg)
590               icount = icount + 1
591               nbidta(icount, igrd, ib_bdy) = ii
592               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
593               nbrdta(icount, igrd, ib_bdy) = ir
594            ENDDO
595         ENDDO
596         !
597         ! ------------ U points -------------
598         igrd=2
599         icount=0
600         DO ir = 1, nn_rimwidth(ib_bdy)
601!            DO ii = jpisdt(iseg), jpisft(iseg) - 1
602            DO ii = jpisdt(iseg), jpisft(iseg)
603               icount = icount + 1
604               nbidta(icount, igrd, ib_bdy) = ii
605               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
606               nbrdta(icount, igrd, ib_bdy) = ir
607            ENDDO
608            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
609            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
610         ENDDO
611         !
612         ! ------------ V points -------------
613         igrd=3
614         icount=0
615         DO ir = 1, nn_rimwidth(ib_bdy)
616            DO ii = jpisdt(iseg), jpisft(iseg)
617               icount = icount + 1
618               nbidta(icount, igrd, ib_bdy) = ii
619               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
620               nbrdta(icount, igrd, ib_bdy) = ir
621            ENDDO
622         ENDDO
623      ENDDO
624
625      !  Deal with duplicated points
626      !-----------------------------
627      ! We assign negative indices to duplicated points (to remove them from bdy points to be updated)
628      ! if their distance to the bdy is greater than the other
629      ! If their distance are the same, just keep only one to avoid updating a point twice
630      DO igrd = 1, jpbgrd
631         DO ib_bdy1 = 1, nb_bdy
632            DO ib_bdy2 = 1, nb_bdy
633               IF (ib_bdy1/=ib_bdy2) THEN
634                  DO ib1 = 1, nblendta(igrd,ib_bdy1)
635                     DO ib2 = 1, nblendta(igrd,ib_bdy2)
636                        IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. &
637                        &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
638!                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &
639!                                                       &              nbidta(ib1, igrd, ib_bdy1),      &
640!                                                       &              nbjdta(ib2, igrd, ib_bdy2)
641                           ! keep only points with the lowest distance to boundary:
642                           IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN
643                             nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2
644                             nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2
645                           ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN
646                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
647                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
648                           ! Arbitrary choice if distances are the same:
649                           ELSE
650                             nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
651                             nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
652                           ENDIF
653                        END IF
654                     END DO
655                  END DO
656               ENDIF
657            END DO
658         END DO
659      END DO
660
661      ! Work out dimensions of boundary data on each processor
662      ! ------------------------------------------------------
663
664      ! Rather assume that boundary data indices are given on global domain
665      ! TO BE DISCUSSED ?
666!      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2
667!      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1
668!      is = mjg(1) + 1            ! if monotasking and no zoom, is=2
669!      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1     
670      iw = mig(1) - jpizoom + 2         ! if monotasking and no zoom, iw=2
671      ie = mig(1) + nlci - jpizoom - 1  ! if monotasking and no zoom, ie=jpim1
672      is = mjg(1) - jpjzoom + 2         ! if monotasking and no zoom, is=2
673      in = mjg(1) + nlcj - jpjzoom - 1  ! if monotasking and no zoom, in=jpjm1
674
675      DO ib_bdy = 1, nb_bdy
676         DO igrd = 1, jpbgrd
677            icount  = 0
678            icountr = 0
679            idx_bdy(ib_bdy)%nblen(igrd)    = 0
680            idx_bdy(ib_bdy)%nblenrim(igrd) = 0
681            DO ib = 1, nblendta(igrd,ib_bdy)
682               ! check that data is in correct order in file
683               ibm1 = MAX(1,ib-1)
684               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc...
685                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
686                     CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', &
687                    'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory')
688                  ENDIF   
689               ENDIF
690               ! check if point is in local domain
691               IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
692                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN
693                  !
694                  icount = icount  + 1
695                  !
696                  IF( nbrdta(ib,igrd,ib_bdy) == 1 )   icountr = icountr+1
697               ENDIF
698            ENDDO
699            idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
700            idx_bdy(ib_bdy)%nblen   (igrd) = icount  !: length of boundary data on each proc       
701         ENDDO  ! igrd
702
703         ! Allocate index arrays for this boundary set
704         !--------------------------------------------
705         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
706         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
707         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
708         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) )
709         ALLOCATE( idx_bdy(ib_bdy)%nbd(ilen1,jpbgrd) )
710         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) )
711         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
712         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) )
713         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )
714
715         ! Dispatch mapping indices and discrete distances on each processor
716         ! -----------------------------------------------------------------
717
718         DO igrd = 1, jpbgrd
719            icount  = 0
720            ! Loop on rimwidth to ensure outermost points come first in the local arrays.
721            DO ir=1, nn_rimwidth(ib_bdy)
722               DO ib = 1, nblendta(igrd,ib_bdy)
723                  ! check if point is in local domain and equals ir
724                  IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
725                     & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   &
726                     & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
727                     !
728                     icount = icount  + 1
729
730                     ! Rather assume that boundary data indices are given on global domain
731                     ! TO BE DISCUSSED ?
732!                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1
733!                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
734                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+jpizoom
735                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+jpjzoom
736                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
737                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
738                  ENDIF
739               ENDDO
740            ENDDO
741         ENDDO 
742
743         ! Compute rim weights for FRS scheme
744         ! ----------------------------------
745         DO igrd = 1, jpbgrd
746            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
747               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
748               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
749!               idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.  ! quadratic
750!               idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy))       ! linear
751            END DO
752         END DO 
753
754         ! Compute damping coefficients
755         ! ----------------------------
756         DO igrd = 1, jpbgrd
757            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
758               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
759               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 
760               & *(FLOAT(nn_rimwidth(ib_bdy)+1-nbr)/FLOAT(nn_rimwidth(ib_bdy)))**2.   ! quadratic
761            END DO
762         END DO
763
764      ENDDO
765
766      ! ------------------------------------------------------
767      ! Initialise masks and find normal/tangential directions
768      ! ------------------------------------------------------
769
770      ! Read global 2D mask at T-points: bdytmask
771      ! -----------------------------------------
772      ! bdytmask = 1  on the computational domain AND on open boundaries
773      !          = 0  elsewhere   
774 
775      IF( ln_mask_file ) THEN
776         CALL iom_open( cn_mask_file, inum )
777         CALL iom_get ( inum, jpdom_data, 'bdy_msk', bdytmask(:,:) )
778         CALL iom_close( inum )
779
780         ! Derive mask on U and V grid from mask on T grid
781         bdyumask(:,:) = 0.e0
782         bdyvmask(:,:) = 0.e0
783         DO ij=1, jpjm1
784            DO ii=1, jpim1
785               bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
786               bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1) 
787            END DO
788         END DO
789         CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
790
791
792         ! Mask corrections
793         ! ----------------
794         DO ik = 1, jpkm1
795            DO ij = 1, jpj
796               DO ii = 1, jpi
797                  tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
798                  umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
799                  vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
800                  bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij)
801               END DO     
802            END DO
803         END DO
804
805         DO ik = 1, jpkm1
806            DO ij = 2, jpjm1
807               DO ii = 2, jpim1
808                  fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   &
809                     &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
810               END DO     
811            END DO
812         END DO
813
814         tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)
815
816      ENDIF ! ln_mask_file=.TRUE.
817     
818      bdytmask(:,:) = tmask(:,:,1)
819
820      ! bdy masks and bmask are now set to zero on boundary points:
821      igrd = 1       ! In the free surface case, bmask is at T-points
822      DO ib_bdy = 1, nb_bdy
823        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
824          bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
825        ENDDO
826      ENDDO
827      !
828      igrd = 1
829      DO ib_bdy = 1, nb_bdy
830        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
831          bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
832        ENDDO
833      ENDDO
834      !
835      igrd = 2
836      DO ib_bdy = 1, nb_bdy
837        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
838          bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
839        ENDDO
840      ENDDO
841      !
842      igrd = 3
843      DO ib_bdy = 1, nb_bdy
844        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
845          bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
846        ENDDO
847      ENDDO
848
849      ! Lateral boundary conditions
850      CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
851      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
852
853      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
854
855         idx_bdy(ib_bdy)%flagu(:) = 0.e0
856         idx_bdy(ib_bdy)%flagv(:) = 0.e0
857         icount = 0 
858
859         !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
860         !flagu =  0 : u is tangential
861         !flagu =  1 : u is normal to the boundary and is direction is inward
862 
863         igrd = 2      ! u-component
864         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
865            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
866            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
867            zefl = bdytmask(nbi  ,nbj)
868            zwfl = bdytmask(nbi+1,nbj)
869            IF( zefl + zwfl == 2 ) THEN
870               icount = icount + 1
871            ELSE
872               idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl
873            ENDIF
874         END DO
875
876         !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
877         !flagv =  0 : u is tangential
878         !flagv =  1 : u is normal to the boundary and is direction is inward
879
880         igrd = 3      ! v-component
881         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
882            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
883            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
884            znfl = bdytmask(nbi,nbj  )
885            zsfl = bdytmask(nbi,nbj+1)
886            IF( znfl + zsfl == 2 ) THEN
887               icount = icount + 1
888            ELSE
889               idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl
890            END IF
891         END DO
892
893         IF( icount /= 0 ) THEN
894            IF(lwp) WRITE(numout,*)
895            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   &
896               ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy
897            IF(lwp) WRITE(numout,*) ' ========== '
898            IF(lwp) WRITE(numout,*)
899            nstop = nstop + 1
900         ENDIF
901   
902      ENDDO
903
904      ! Compute total lateral surface for volume correction:
905      ! ----------------------------------------------------
906      ! JC: this must be done at each time step with key_vvl
907      bdysurftot = 0.e0 
908      IF( ln_vol ) THEN 
909         igrd = 2      ! Lateral surface at U-points
910         DO ib_bdy = 1, nb_bdy
911            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
912               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
913               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
914               flagu => idx_bdy(ib_bdy)%flagu(ib)
915               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           &
916                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) &
917                  &                    * tmask_i(nbi  , nbj)                           &
918                  &                    * tmask_i(nbi+1, nbj)                   
919            ENDDO
920         ENDDO
921
922         igrd=3 ! Add lateral surface at V-points
923         DO ib_bdy = 1, nb_bdy
924            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
925               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
926               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
927               flagv => idx_bdy(ib_bdy)%flagv(ib)
928               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           &
929                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) &
930                  &                    * tmask_i(nbi, nbj  )                           &
931                  &                    * tmask_i(nbi, nbj+1)
932            ENDDO
933         ENDDO
934         !
935         IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain
936      END IF   
937      !
938      ! Tidy up
939      !--------
940      IF (nb_bdy>0) THEN
941         DEALLOCATE(nbidta, nbjdta, nbrdta)
942      ENDIF
943
944      IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
945
946   END SUBROUTINE bdy_init
947
948   SUBROUTINE bdy_ctl_seg
949      !!----------------------------------------------------------------------
950      !!                 ***  ROUTINE bdy_ctl_seg  ***
951      !!
952      !! ** Purpose :   Check straight open boundary segments location
953      !!
954      !! ** Method  :   - Look for open boundary corners
955      !!                - Check that segments start or end on land
956      !!----------------------------------------------------------------------
957      INTEGER  ::   ib, ib1, ib2, ji ,jj, itest 
958      INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns 
959      REAL(wp), DIMENSION(2) ::   ztestmask
960      !!----------------------------------------------------------------------
961      !
962      IF (lwp) WRITE(numout,*) ' '
963      IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
964      IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
965      !
966      IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege
967      IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw
968      IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn
969      IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs
970      ! 1. Check bounds
971      !----------------
972      DO ib = 1, nbdysegn
973         IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib)
974         IF ((jpjnob(ib).ge.jpjglo-1).or.& 
975            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
976         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
977         IF (jpindt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
978         IF (jpinft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' )
979      END DO
980      !
981      DO ib = 1, nbdysegs
982         IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib)
983         IF ((jpjsob(ib).ge.jpjglo-1).or.& 
984            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
985         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
986         IF (jpisdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
987         IF (jpisft(ib).ge.jpiglo)     CALL ctl_stop( 'End index out of domain' )
988      END DO
989      !
990      DO ib = 1, nbdysege
991         IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib)
992         IF ((jpieob(ib).ge.jpiglo-1).or.& 
993            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
994         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
995         IF (jpjedt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
996         IF (jpjeft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' )
997      END DO
998      !
999      DO ib = 1, nbdysegw
1000         IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib)
1001         IF ((jpiwob(ib).ge.jpiglo-1).or.& 
1002            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1003         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1004         IF (jpjwdt(ib).le.1     )     CALL ctl_stop( 'Start index out of domain' )
1005         IF (jpjwft(ib).ge.jpjglo)     CALL ctl_stop( 'End index out of domain' )
1006      ENDDO
1007      !
1008      !     
1009      ! 2. Look for segment crossings
1010      !------------------------------
1011      IF (lwp) WRITE(numout,*) '**Look for segments corners  :'
1012      !
1013      itest = 0 ! corner number
1014      !
1015      ! flag to detect if start or end of open boundary belongs to a corner
1016      ! if not (=0), it must be on land.
1017      ! if a corner is detected, save bdy package number for further tests
1018      icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1019      ! South/West crossings
1020      IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1021         DO ib1 = 1, nbdysegw       
1022            DO ib2 = 1, nbdysegs
1023               IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1024                &  ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1025                &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1026                &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1027                  IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN 
1028                     ! We have a possible South-West corner                     
1029!                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1030!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1031                     icornw(ib1,1) = npckgs(ib2)
1032                     icorns(ib2,1) = npckgw(ib1)
1033                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1034                     IF(lwp) WRITE(numout,*)
1035                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1036                     &                                     jpisft(ib2), jpjwft(ib1)
1037                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1038                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), & 
1039                     &                                                    ' and South segment: ',npckgs(ib2)
1040                     IF(lwp) WRITE(numout,*)
1041                     nstop = nstop + 1
1042                  ELSE
1043                     IF(lwp) WRITE(numout,*)
1044                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and West Open boundary indices'
1045                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1) , &
1046                     &                                                    ' and South segment: ',npckgs(ib2)
1047                     IF(lwp) WRITE(numout,*)
1048                     nstop = nstop+1
1049                  END IF
1050               END IF
1051            END DO
1052         END DO
1053      END IF
1054      !
1055      ! South/East crossings
1056      IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1057         DO ib1 = 1, nbdysege
1058            DO ib2 = 1, nbdysegs
1059               IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1060                &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1061                &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. &
1062                &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN
1063                  IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1064                     ! We have a possible South-East corner
1065!                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1066!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1067                     icorne(ib1,1) = npckgs(ib2)
1068                     icorns(ib2,2) = npckge(ib1)
1069                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1070                     IF(lwp) WRITE(numout,*)
1071                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1072                     &                                     jpisdt(ib2), jpjeft(ib1)
1073                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1074                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), &
1075                     &                                                    ' and South segment: ',npckgs(ib2)
1076                     IF(lwp) WRITE(numout,*)
1077                     nstop = nstop + 1
1078                  ELSE
1079                     IF(lwp) WRITE(numout,*)
1080                     IF(lwp) WRITE(numout,*) ' E R R O R : Check South and East Open boundary indices'
1081                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), &
1082                     &                                                    ' and South segment: ',npckgs(ib2)
1083                     IF(lwp) WRITE(numout,*)
1084                     nstop = nstop + 1
1085                  END IF
1086               END IF
1087            END DO
1088         END DO
1089      END IF
1090      !
1091      ! North/West crossings
1092      IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1093         DO ib1 = 1, nbdysegw       
1094            DO ib2 = 1, nbdysegn
1095               IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. &
1096                &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. &
1097                &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1098                &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1099                  IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1100                     ! We have a possible North-West corner
1101!                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1102!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1103                     icornw(ib1,2) = npckgn(ib2)
1104                     icornn(ib2,1) = npckgw(ib1)
1105                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1106                     IF(lwp) WRITE(numout,*)
1107                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1108                     &                                     jpinft(ib2), jpjwdt(ib1)
1109                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1110                     IF(lwp) WRITE(numout,*) '             Crossing problem with West segment: ',npckgw(ib1), &
1111                     &                                                    ' and North segment: ',npckgn(ib2)
1112                     IF(lwp) WRITE(numout,*)
1113                     nstop = nstop + 1
1114                  ELSE
1115                     IF(lwp) WRITE(numout,*)
1116                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and West Open boundary indices'
1117                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1), &
1118                     &                                                    ' and North segment: ',npckgn(ib2)
1119                     IF(lwp) WRITE(numout,*)
1120                     nstop = nstop + 1
1121                  END IF
1122               END IF
1123            END DO
1124         END DO
1125      END IF
1126      !
1127      ! North/East crossings
1128      IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1129         DO ib1 = 1, nbdysege       
1130            DO ib2 = 1, nbdysegn
1131               IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1132                &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1133                &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1134                &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1135                  IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1136                     ! We have a possible North-East corner
1137!                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1138!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1139                     icorne(ib1,2) = npckgn(ib2)
1140                     icornn(ib2,2) = npckge(ib1)
1141                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1142                     IF(lwp) WRITE(numout,*)
1143                     IF(lwp) WRITE(numout,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
1144                     &                                     jpindt(ib2), jpjedt(ib1)
1145                     IF(lwp) WRITE(numout,*) ' ==========  Not allowed yet'
1146                     IF(lwp) WRITE(numout,*) '             Crossing problem with East segment: ',npckge(ib1), &
1147                     &                                                    ' and North segment: ',npckgn(ib2)
1148                     IF(lwp) WRITE(numout,*)
1149                     nstop = nstop + 1
1150                  ELSE
1151                     IF(lwp) WRITE(numout,*)
1152                     IF(lwp) WRITE(numout,*) ' E R R O R : Check North and East Open boundary indices'
1153                     IF(lwp) WRITE(numout,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), &
1154                     &                                                    ' and North segment: ',npckgn(ib2)
1155                     IF(lwp) WRITE(numout,*)
1156                     nstop = nstop + 1
1157                  END IF
1158               END IF
1159            END DO
1160         END DO
1161      END IF
1162      !
1163      ! 3. Check if segment extremities are on land
1164      !--------------------------------------------
1165      !
1166      ! West segments
1167      DO ib = 1, nbdysegw
1168         ! get mask at boundary extremities:
1169         ztestmask(1:2)=0.
1170         DO ji = 1, jpi
1171            DO jj = 1, jpj             
1172              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
1173               &  ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1174              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
1175               &  ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1176            END DO
1177         END DO
1178         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1179
1180         IF (ztestmask(1)==1) THEN
1181            IF (icornw(ib,1)==0) THEN
1182               IF(lwp) WRITE(numout,*)
1183               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
1184               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                 
1185               IF(lwp) WRITE(numout,*)
1186               nstop = nstop + 1
1187            ELSE
1188               ! This is a corner
1189               WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1190               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1191               itest=itest+1
1192            ENDIF
1193         ENDIF
1194         IF (ztestmask(2)==1) THEN
1195            IF (icornw(ib,2)==0) THEN
1196               IF(lwp) WRITE(numout,*)
1197               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgw(ib)
1198               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                 
1199               IF(lwp) WRITE(numout,*)
1200               nstop = nstop + 1
1201            ELSE
1202               ! This is a corner
1203               WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
1204               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
1205               itest=itest+1
1206            ENDIF
1207         ENDIF
1208      END DO
1209      !
1210      ! East segments
1211      DO ib = 1, nbdysege
1212         ! get mask at boundary extremities:
1213         ztestmask(1:2)=0.
1214         DO ji = 1, jpi
1215            DO jj = 1, jpj             
1216              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
1217               &  ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
1218              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
1219               &  ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1220            END DO
1221         END DO
1222         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1223
1224         IF (ztestmask(1)==1) THEN
1225            IF (icorne(ib,1)==0) THEN
1226               IF(lwp) WRITE(numout,*)
1227               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
1228               IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                 
1229               IF(lwp) WRITE(numout,*)
1230               nstop = nstop + 1 
1231            ELSE
1232               ! This is a corner
1233               WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1234               CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1235               itest=itest+1
1236            ENDIF
1237         ENDIF
1238         IF (ztestmask(2)==1) THEN
1239            IF (icorne(ib,2)==0) THEN
1240               IF(lwp) WRITE(numout,*)
1241               IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckge(ib)
1242               IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                 
1243               IF(lwp) WRITE(numout,*)
1244               nstop = nstop + 1
1245            ELSE
1246               ! This is a corner
1247               WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1248               CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1249               itest=itest+1
1250            ENDIF
1251         ENDIF
1252      END DO
1253      !
1254      ! South segments
1255      DO ib = 1, nbdysegs
1256         ! get mask at boundary extremities:
1257         ztestmask(1:2)=0.
1258         DO ji = 1, jpi
1259            DO jj = 1, jpj             
1260              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
1261               &  ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1262              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
1263               &  ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1264            END DO
1265         END DO
1266         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1267
1268         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1269            IF(lwp) WRITE(numout,*)
1270            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
1271            IF(lwp) WRITE(numout,*) ' ==========  does not start on land or on a corner'                                                 
1272            IF(lwp) WRITE(numout,*)
1273            nstop = nstop + 1
1274         ENDIF
1275         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1276            IF(lwp) WRITE(numout,*)
1277            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgs(ib)
1278            IF(lwp) WRITE(numout,*) ' ==========  does not end on land or on a corner'                                                 
1279            IF(lwp) WRITE(numout,*)
1280            nstop = nstop + 1
1281         ENDIF
1282      END DO
1283      !
1284      ! North segments
1285      DO ib = 1, nbdysegn
1286         ! get mask at boundary extremities:
1287         ztestmask(1:2)=0.
1288         DO ji = 1, jpi
1289            DO jj = 1, jpj             
1290              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
1291               &  ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
1292              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
1293               &  ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1294            END DO
1295         END DO
1296         IF( lk_mpp )   CALL mpp_sum( ztestmask, 2 )   ! sum over the global domain
1297
1298         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1299            IF(lwp) WRITE(numout,*)
1300            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
1301            IF(lwp) WRITE(numout,*) ' ==========  does not start on land'                                                 
1302            IF(lwp) WRITE(numout,*)
1303            nstop = nstop + 1
1304         ENDIF
1305         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1306            IF(lwp) WRITE(numout,*)
1307            IF(lwp) WRITE(numout,*) ' E R R O R : Open boundary segment ', npckgn(ib)
1308            IF(lwp) WRITE(numout,*) ' ==========  does not end on land'                                                 
1309            IF(lwp) WRITE(numout,*)
1310            nstop = nstop + 1
1311         ENDIF
1312      END DO
1313      !
1314      IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1315      !
1316      ! Other tests TBD:
1317      ! segments completly on land
1318      ! optimized open boundary array length according to landmask
1319      ! Nudging layers that overlap with interior domain
1320      !
1321   END SUBROUTINE bdy_ctl_seg
1322
1323   SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1324      !!----------------------------------------------------------------------
1325      !!                 ***  ROUTINE bdy_ctl_corn  ***
1326      !!
1327      !! ** Purpose :   Check numerical schemes consistency between
1328      !!                segments having a common corner
1329      !!
1330      !! ** Method  :   
1331      !!----------------------------------------------------------------------
1332      INTEGER, INTENT(in)  ::   ib1, ib2
1333      INTEGER :: itest
1334      !!----------------------------------------------------------------------
1335      itest = 0
1336
1337      IF (nn_dyn2d(ib1)/=nn_dyn2d(ib2)) itest = itest + 1
1338      IF (nn_dyn3d(ib1)/=nn_dyn3d(ib2)) itest = itest + 1
1339      IF (nn_tra(ib1)/=nn_tra(ib2)) itest = itest + 1
1340      !
1341      IF (nn_dyn2d_dta(ib1)/=nn_dyn2d_dta(ib2)) itest = itest + 1
1342      IF (nn_dyn3d_dta(ib1)/=nn_dyn3d_dta(ib2)) itest = itest + 1
1343      IF (nn_tra_dta(ib1)/=nn_tra_dta(ib2)) itest = itest + 1
1344      !
1345      IF (nn_rimwidth(ib1)/=nn_rimwidth(ib2)) itest = itest + 1   
1346      !
1347      IF ( itest>0 ) THEN
1348         IF(lwp) WRITE(numout,*) ' E R R O R : Segments ', ib1, 'and ', ib2
1349         IF(lwp) WRITE(numout,*) ' ==========  have different open bdy schemes'                                                 
1350         IF(lwp) WRITE(numout,*)
1351         nstop = nstop + 1
1352      ENDIF
1353      !
1354   END SUBROUTINE bdy_ctl_corn
1355
1356#else
1357   !!---------------------------------------------------------------------------------
1358   !!   Dummy module                                   NO open boundaries
1359   !!---------------------------------------------------------------------------------
1360CONTAINS
1361   SUBROUTINE bdy_init      ! Dummy routine
1362   END SUBROUTINE bdy_init
1363#endif
1364
1365   !!=================================================================================
1366END MODULE bdyini
Note: See TracBrowser for help on using the repository browser.