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

source: branches/2012/dev_CMCC_INGV_2012/NEMOGCM/NEMO/OPA_SRC/BDY/bdyini.F90 @ 3646

Last change on this file since 3646 was 3646, checked in by vichi, 11 years ago

Add the resulting merged branch from CMCC and INGV 2012 developments
  • Property svn:keywords set to Id
File size: 46.9 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.5  !  2012     (S. Mocavero, I. Epicoco) Updates for the
14   !!                             optimization of BDY communications
15   !!----------------------------------------------------------------------
16#if defined key_bdy
17   !!----------------------------------------------------------------------
18   !!   'key_bdy'                     Unstructured Open Boundary Conditions
19   !!----------------------------------------------------------------------
20   !!   bdy_init       : Initialization of unstructured open boundaries
21   !!----------------------------------------------------------------------
22   USE timing          ! Timing
23   USE oce             ! ocean dynamics and tracers variables
24   USE dom_oce         ! ocean space and time domain
25   USE bdy_oce         ! unstructured open boundary conditions
26   USE in_out_manager  ! I/O units
27   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
28   USE lib_mpp         ! for mpp_sum 
29   USE iom             ! I/O
30
31   IMPLICIT NONE
32   PRIVATE
33
34   PUBLIC   bdy_init   ! routine called in nemo_init
35
36   !!----------------------------------------------------------------------
37   !! NEMO/OPA 4.0 , NEMO Consortium (2011)
38   !! $Id$
39   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
40   !!----------------------------------------------------------------------
41CONTAINS
42   
43   SUBROUTINE bdy_init
44      !!----------------------------------------------------------------------
45      !!                 ***  ROUTINE bdy_init  ***
46      !!         
47      !! ** Purpose :   Initialization of the dynamics and tracer fields with
48      !!              unstructured open boundaries.
49      !!
50      !! ** Method  :   Read initialization arrays (mask, indices) to identify
51      !!              an unstructured open boundary
52      !!
53      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
54      !!----------------------------------------------------------------------     
55      ! namelist variables
56      !-------------------
57      INTEGER, PARAMETER          :: jp_nseg = 100
58      INTEGER                     :: nbdysege, nbdysegw, nbdysegn, nbdysegs 
59      INTEGER, DIMENSION(jp_nseg) :: jpieob, jpjedt, jpjeft
60      INTEGER, DIMENSION(jp_nseg) :: jpiwob, jpjwdt, jpjwft
61      INTEGER, DIMENSION(jp_nseg) :: jpjnob, jpindt, jpinft
62      INTEGER, DIMENSION(jp_nseg) :: jpjsob, jpisdt, jpisft
63
64      ! local variables
65      !-------------------
66      INTEGER  ::   ib_bdy, ii, ij, ik, igrd, ib, ir, iseg ! dummy loop indices
67      INTEGER  ::   icount, icountr, ibr_max, ilen1, ibm1  ! local integers
68      INTEGER  ::   iw, ie, is, in, inum, id_dummy         !   -       -
69      INTEGER  ::   igrd_start, igrd_end, jpbdta           !   -       -
70      INTEGER, POINTER  ::  nbi, nbj, nbr                  ! short cuts
71      REAL   , POINTER  ::  flagu, flagv                   !    -   -
72      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
73      INTEGER, DIMENSION (2)                ::   kdimsz
74      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays
75      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta
76      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
77      REAL(wp), DIMENSION(jpidta,jpjdta)    ::   zmask            ! global domain mask
78      CHARACTER(LEN=80),DIMENSION(jpbgrd)  ::   clfile
79      CHARACTER(LEN=1),DIMENSION(jpbgrd)   ::   cgrid
80      INTEGER :: com_east, com_west, com_south, com_north          ! Flags for boundaries sending
81      INTEGER :: com_east_b, com_west_b, com_south_b, com_north_b  ! Flags for boundaries receiving
82      INTEGER :: iw_b(4), ie_b(4), is_b(4), in_b(4)                ! Arrays for neighbours coordinates
83
84      !!
85      NAMELIST/nambdy/ nb_bdy, ln_coords_file, cn_coords_file,             &
86         &             ln_mask_file, cn_mask_file, nn_dyn2d, nn_dyn2d_dta, &
87         &             nn_dyn3d, nn_dyn3d_dta, nn_tra, nn_tra_dta,         & 
88#if defined key_lim2
89         &             nn_ice_lim2, nn_ice_lim2_dta,                       &
90#endif
91         &             ln_vol, nn_volctl, nn_rimwidth
92      !!
93      NAMELIST/nambdy_index/ nbdysege, jpieob, jpjedt, jpjeft,             &
94                             nbdysegw, jpiwob, jpjwdt, jpjwft,             &
95                             nbdysegn, jpjnob, jpindt, jpinft,             &
96                             nbdysegs, jpjsob, jpisdt, jpisft
97
98      !!----------------------------------------------------------------------
99
100      IF( nn_timing == 1 ) CALL timing_start('bdy_init')
101
102      IF( bdy_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'bdy_init : unable to allocate oce arrays' )
103
104      IF(lwp) WRITE(numout,*)
105      IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
106      IF(lwp) WRITE(numout,*) '~~~~~~~~'
107      !
108
109      IF( jperio /= 0 )   CALL ctl_stop( 'Cyclic or symmetric,',   &
110         &                               ' and general open boundary condition are not compatible' )
111
112      cgrid= (/'t','u','v'/)
113
114      ! -----------------------------------------
115      ! Initialise and read namelist parameters
116      ! -----------------------------------------
117
118      nb_bdy            = 0
119      ln_coords_file(:) = .false.
120      cn_coords_file(:) = ''
121      ln_mask_file      = .false.
122      cn_mask_file(:)   = ''
123      nn_dyn2d(:)       = 0
124      nn_dyn2d_dta(:)   = -1  ! uninitialised flag
125      nn_dyn3d(:)       = 0
126      nn_dyn3d_dta(:)   = -1  ! uninitialised flag
127      nn_tra(:)         = 0
128      nn_tra_dta(:)     = -1  ! uninitialised flag
129#if defined key_lim2
130      nn_ice_lim2(:)    = 0
131      nn_ice_lim2_dta(:)= -1  ! uninitialised flag
132#endif
133      ln_vol            = .false.
134      nn_volctl         = -1  ! uninitialised flag
135      nn_rimwidth(:)    = -1  ! uninitialised flag
136
137      REWIND( numnam )                   
138      READ  ( numnam, nambdy )
139
140      ! -----------------------------------------
141      ! Check and write out namelist parameters
142      ! -----------------------------------------
143
144      !                                   ! control prints
145      IF(lwp) WRITE(numout,*) '         nambdy'
146
147      IF( nb_bdy .eq. 0 ) THEN
148        IF(lwp) WRITE(numout,*) 'nb_bdy = 0, NO OPEN BOUNDARIES APPLIED.'
149      ELSE
150        IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ',nb_bdy
151      ENDIF
152
153      DO ib_bdy = 1,nb_bdy
154        IF(lwp) WRITE(numout,*) ' ' 
155        IF(lwp) WRITE(numout,*) '------ Open boundary data set ',ib_bdy,'------' 
156
157        IF( ln_coords_file(ib_bdy) ) THEN
158           IF(lwp) WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
159        ELSE
160           IF(lwp) WRITE(numout,*) 'Boundary defined in namelist.'
161        ENDIF
162        IF(lwp) WRITE(numout,*)
163
164        IF(lwp) WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
165        SELECT CASE( nn_dyn2d(ib_bdy) )                 
166          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
167          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
168          CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      Flather radiation condition'
169          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn2d' )
170        END SELECT
171        IF( nn_dyn2d(ib_bdy) .gt. 0 ) THEN
172           SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !
173              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
174              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
175              CASE( 2 )      ;   IF(lwp) WRITE(numout,*) '      tidal harmonic forcing taken from file'
176              CASE( 3 )      ;   IF(lwp) WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files'
177              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
178           END SELECT
179        ENDIF
180        IF(lwp) WRITE(numout,*)
181
182        IF(lwp) WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
183        SELECT CASE( nn_dyn3d(ib_bdy) )                 
184          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
185          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
186          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_dyn3d' )
187        END SELECT
188        IF( nn_dyn3d(ib_bdy) .gt. 0 ) THEN
189           SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !
190              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
191              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
192              CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
193           END SELECT
194        ENDIF
195        IF(lwp) WRITE(numout,*)
196
197        IF(lwp) WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
198        SELECT CASE( nn_tra(ib_bdy) )                 
199          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
200          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
201          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
202        END SELECT
203        IF( nn_tra(ib_bdy) .gt. 0 ) THEN
204           SELECT CASE( nn_tra_dta(ib_bdy) )                   !
205              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
206              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
207              CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
208           END SELECT
209        ENDIF
210        IF(lwp) WRITE(numout,*)
211
212#if defined key_lim2
213        IF(lwp) WRITE(numout,*) 'Boundary conditions for sea ice:  '
214        SELECT CASE( nn_ice_lim2(ib_bdy) )                 
215          CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      no open boundary condition'       
216          CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      Flow Relaxation Scheme'
217          CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for nn_tra' )
218        END SELECT
219        IF( nn_ice_lim2(ib_bdy) .gt. 0 ) THEN
220           SELECT CASE( nn_ice_lim2_dta(ib_bdy) )                   !
221              CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      initial state used for bdy data'       
222              CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      boundary data taken from file'
223              CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_lim2_dta must be 0 or 1' )
224           END SELECT
225        ENDIF
226        IF(lwp) WRITE(numout,*)
227#endif
228
229        IF(lwp) WRITE(numout,*) 'Boundary rim width for the FRS scheme = ', nn_rimwidth(ib_bdy)
230        IF(lwp) WRITE(numout,*)
231
232      ENDDO
233
234     IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
235       IF(lwp) WRITE(numout,*) 'Volume correction applied at open boundaries'
236       IF(lwp) WRITE(numout,*)
237       SELECT CASE ( nn_volctl )
238         CASE( 1 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will be constant'
239         CASE( 0 )      ;   IF(lwp) WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
240         CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
241       END SELECT
242       IF(lwp) WRITE(numout,*)
243     ELSE
244       IF(lwp) WRITE(numout,*) 'No volume correction applied at open boundaries'
245       IF(lwp) WRITE(numout,*)
246     ENDIF
247
248      ! -------------------------------------------------
249      ! Initialise indices arrays for open boundaries
250      ! -------------------------------------------------
251
252      ! Work out global dimensions of boundary data
253      ! ---------------------------------------------
254      REWIND( numnam )                   
255      DO ib_bdy = 1, nb_bdy
256
257         jpbdta = 1
258         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Work out size of global arrays from namelist parameters
259 
260            ! No REWIND here because may need to read more than one nambdy_index namelist.
261            READ  ( numnam, nambdy_index )
262
263            ! Automatic boundary definition: if nbdysegX = -1
264            ! set boundary to whole side of model domain.
265            IF( nbdysege == -1 ) THEN
266               nbdysege = 1
267               jpieob(1) = jpiglo - 1
268               jpjedt(1) = 2
269               jpjeft(1) = jpjglo - 1
270            ENDIF
271            IF( nbdysegw == -1 ) THEN
272               nbdysegw = 1
273               jpiwob(1) = 2
274               jpjwdt(1) = 2
275               jpjwft(1) = jpjglo - 1
276            ENDIF
277            IF( nbdysegn == -1 ) THEN
278               nbdysegn = 1
279               jpjnob(1) = jpjglo - 1
280               jpindt(1) = 2
281               jpinft(1) = jpiglo - 1
282            ENDIF
283            IF( nbdysegs == -1 ) THEN
284               nbdysegs = 1
285               jpjsob(1) = 2
286               jpisdt(1) = 2
287               jpisft(1) = jpiglo - 1
288            ENDIF
289
290            nblendta(:,ib_bdy) = 0
291            DO iseg = 1, nbdysege
292               igrd = 1
293               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1               
294               igrd = 2
295               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg) + 1               
296               igrd = 3
297               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjeft(iseg) - jpjedt(iseg)               
298            ENDDO
299            DO iseg = 1, nbdysegw
300               igrd = 1
301               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1               
302               igrd = 2
303               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg) + 1               
304               igrd = 3
305               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpjwft(iseg) - jpjwdt(iseg)               
306            ENDDO
307            DO iseg = 1, nbdysegn
308               igrd = 1
309               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1               
310               igrd = 2
311               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg)               
312               igrd = 3
313               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpinft(iseg) - jpindt(iseg) + 1
314            ENDDO
315            DO iseg = 1, nbdysegs
316               igrd = 1
317               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1               
318               igrd = 2
319               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg)
320               igrd = 3
321               nblendta(igrd,ib_bdy) = nblendta(igrd,ib_bdy) + jpisft(iseg) - jpisdt(iseg) + 1               
322            ENDDO
323
324            nblendta(:,ib_bdy) = nblendta(:,ib_bdy) * nn_rimwidth(ib_bdy)
325            jpbdta = MAXVAL(nblendta(:,ib_bdy))               
326
327
328         ELSE            ! Read size of arrays in boundary coordinates file.
329
330
331            CALL iom_open( cn_coords_file(ib_bdy), inum )
332            jpbdta = 1
333            DO igrd = 1, jpbgrd
334               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 
335               nblendta(igrd,ib_bdy) = kdimsz(1)
336               jpbdta = MAX(jpbdta, kdimsz(1))
337            ENDDO
338
339         ENDIF
340
341      ENDDO ! ib_bdy
342
343      ! Allocate arrays
344      !---------------
345      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy),    &
346         &      nbrdta(jpbdta, jpbgrd, nb_bdy) )
347
348      ALLOCATE( dta_global(jpbdta, 1, jpk) )
349
350      ! Calculate global boundary index arrays or read in from file
351      !------------------------------------------------------------
352      REWIND( numnam )                   
353      DO ib_bdy = 1, nb_bdy
354
355         IF( .NOT. ln_coords_file(ib_bdy) ) THEN ! Calculate global index arrays from namelist parameters
356
357            ! No REWIND here because may need to read more than one nambdy_index namelist.
358            READ  ( numnam, nambdy_index )
359
360            ! Automatic boundary definition: if nbdysegX = -1
361            ! set boundary to whole side of model domain.
362            IF( nbdysege == -1 ) THEN
363               nbdysege = 1
364               jpieob(1) = jpiglo - 1
365               jpjedt(1) = 2
366               jpjeft(1) = jpjglo - 1
367            ENDIF
368            IF( nbdysegw == -1 ) THEN
369               nbdysegw = 1
370               jpiwob(1) = 2
371               jpjwdt(1) = 2
372               jpjwft(1) = jpjglo - 1
373            ENDIF
374            IF( nbdysegn == -1 ) THEN
375               nbdysegn = 1
376               jpjnob(1) = jpjglo - 1
377               jpindt(1) = 2
378               jpinft(1) = jpiglo - 1
379            ENDIF
380            IF( nbdysegs == -1 ) THEN
381               nbdysegs = 1
382               jpjsob(1) = 2
383               jpisdt(1) = 2
384               jpisft(1) = jpiglo - 1
385            ENDIF
386
387            ! ------------ T points -------------
388            igrd = 1 
389            icount = 0
390            DO ir = 1, nn_rimwidth(ib_bdy)
391               ! east
392               DO iseg = 1, nbdysege
393                  DO ij = jpjedt(iseg), jpjeft(iseg)
394                     icount = icount + 1
395                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
396                     nbjdta(icount, igrd, ib_bdy) = ij
397                     nbrdta(icount, igrd, ib_bdy) = ir
398                  ENDDO
399               ENDDO
400               ! west
401               DO iseg = 1, nbdysegw
402                  DO ij = jpjwdt(iseg), jpjwft(iseg)
403                     icount = icount + 1
404                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
405                     nbjdta(icount, igrd, ib_bdy) = ij
406                     nbrdta(icount, igrd, ib_bdy) = ir
407                  ENDDO
408               ENDDO
409               ! north
410               DO iseg = 1, nbdysegn
411                  DO ii = jpindt(iseg), jpinft(iseg)
412                     icount = icount + 1
413                     nbidta(icount, igrd, ib_bdy) = ii
414                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
415                     nbrdta(icount, igrd, ib_bdy) = ir
416                  ENDDO
417               ENDDO
418               ! south
419               DO iseg = 1, nbdysegs
420                  DO ii = jpisdt(iseg), jpisft(iseg)
421                     icount = icount + 1
422                     nbidta(icount, igrd, ib_bdy) = ii
423                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
424                     nbrdta(icount, igrd, ib_bdy) = ir
425                  ENDDO
426               ENDDO
427            ENDDO
428
429            ! ------------ U points -------------
430            igrd = 2 
431            icount = 0
432            DO ir = 1, nn_rimwidth(ib_bdy)
433               ! east
434               DO iseg = 1, nbdysege
435                  DO ij = jpjedt(iseg), jpjeft(iseg)
436                     icount = icount + 1
437                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir
438                     nbjdta(icount, igrd, ib_bdy) = ij
439                     nbrdta(icount, igrd, ib_bdy) = ir
440                  ENDDO
441               ENDDO
442               ! west
443               DO iseg = 1, nbdysegw
444                  DO ij = jpjwdt(iseg), jpjwft(iseg)
445                     icount = icount + 1
446                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
447                     nbjdta(icount, igrd, ib_bdy) = ij
448                     nbrdta(icount, igrd, ib_bdy) = ir
449                  ENDDO
450               ENDDO
451               ! north
452               DO iseg = 1, nbdysegn
453                  DO ii = jpindt(iseg), jpinft(iseg) - 1
454                     icount = icount + 1
455                     nbidta(icount, igrd, ib_bdy) = ii
456                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir + 1
457                     nbrdta(icount, igrd, ib_bdy) = ir
458                  ENDDO
459               ENDDO
460               ! south
461               DO iseg = 1, nbdysegs
462                  DO ii = jpisdt(iseg), jpisft(iseg) - 1
463                     icount = icount + 1
464                     nbidta(icount, igrd, ib_bdy) = ii
465                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
466                     nbrdta(icount, igrd, ib_bdy) = ir
467                  ENDDO
468               ENDDO
469            ENDDO
470
471            ! ------------ V points -------------
472            igrd = 3 
473            icount = 0
474            DO ir = 1, nn_rimwidth(ib_bdy)
475               ! east
476               DO iseg = 1, nbdysege
477                  DO ij = jpjedt(iseg), jpjeft(iseg) - 1
478                     icount = icount + 1
479                     nbidta(icount, igrd, ib_bdy) = jpieob(iseg) - ir + 1
480                     nbjdta(icount, igrd, ib_bdy) = ij
481                     nbrdta(icount, igrd, ib_bdy) = ir
482                  ENDDO
483               ENDDO
484               ! west
485               DO iseg = 1, nbdysegw
486                  DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
487                     icount = icount + 1
488                     nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
489                     nbjdta(icount, igrd, ib_bdy) = ij
490                     nbrdta(icount, igrd, ib_bdy) = ir
491                  ENDDO
492               ENDDO
493               ! north
494               DO iseg = 1, nbdysegn
495                  DO ii = jpindt(iseg), jpinft(iseg)
496                     icount = icount + 1
497                     nbidta(icount, igrd, ib_bdy) = ii
498                     nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) - ir
499                     nbrdta(icount, igrd, ib_bdy) = ir
500                  ENDDO
501               ENDDO
502               ! south
503               DO iseg = 1, nbdysegs
504                  DO ii = jpisdt(iseg), jpisft(iseg)
505                     icount = icount + 1
506                     nbidta(icount, igrd, ib_bdy) = ii
507                     nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
508                     nbrdta(icount, igrd, ib_bdy) = ir
509                  ENDDO
510               ENDDO
511            ENDDO
512
513         ELSE            ! Read global index arrays from boundary coordinates file.
514
515            DO igrd = 1, jpbgrd
516               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
517               DO ii = 1,nblendta(igrd,ib_bdy)
518                  nbidta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
519               END DO
520               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
521               DO ii = 1,nblendta(igrd,ib_bdy)
522                  nbjdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
523               END DO
524               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), dta_global(1:nblendta(igrd,ib_bdy),:,1) )
525               DO ii = 1,nblendta(igrd,ib_bdy)
526                  nbrdta(ii,igrd,ib_bdy) = INT( dta_global(ii,1,1) )
527               END DO
528
529               ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
530               IF(lwp) WRITE(numout,*)
531               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
532               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
533               IF (ibr_max < nn_rimwidth(ib_bdy))   &
534                     CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
535
536            END DO
537            CALL iom_close( inum )
538
539         ENDIF
540
541      ENDDO 
542
543      ! Work out dimensions of boundary data on each processor
544      ! ------------------------------------------------------
545     
546      iw = mig(1) + 1            ! if monotasking and no zoom, iw=2
547      ie = mig(1) + nlci-1 - 1   ! if monotasking and no zoom, ie=jpim1
548      is = mjg(1) + 1            ! if monotasking and no zoom, is=2
549      in = mjg(1) + nlcj-1 - 1   ! if monotasking and no zoom, in=jpjm1
550
551      ALLOCATE( nbondi_bdy(nb_bdy))
552      ALLOCATE( nbondj_bdy(nb_bdy))
553      nbondi_bdy(:)=2
554      nbondj_bdy(:)=2
555      ALLOCATE( nbondi_bdy_b(nb_bdy))
556      ALLOCATE( nbondj_bdy_b(nb_bdy))
557      nbondi_bdy_b(:)=2
558      nbondj_bdy_b(:)=2
559
560      ! Work out dimensions of boundary data on each neighbour process
561      IF(nbondi .eq. 0) THEN
562         iw_b(1) = jpizoom + nimppt(nowe+1)
563         ie_b(1) = jpizoom + nimppt(nowe+1)+nlcit(nowe+1)-3
564         is_b(1) = jpjzoom + njmppt(nowe+1)
565         in_b(1) = jpjzoom + njmppt(nowe+1)+nlcjt(nowe+1)-3
566
567         iw_b(2) = jpizoom + nimppt(noea+1)
568         ie_b(2) = jpizoom + nimppt(noea+1)+nlcit(noea+1)-3
569         is_b(2) = jpjzoom + njmppt(noea+1)
570         in_b(2) = jpjzoom + njmppt(noea+1)+nlcjt(noea+1)-3
571      ELSEIF(nbondi .eq. 1) THEN
572         iw_b(1) = jpizoom + nimppt(nowe+1)
573         ie_b(1) = jpizoom + nimppt(nowe+1)+nlcit(nowe+1)-3
574         is_b(1) = jpjzoom + njmppt(nowe+1)
575         in_b(1) = jpjzoom + njmppt(nowe+1)+nlcjt(nowe+1)-3
576      ELSEIF(nbondi .eq. -1) THEN
577         iw_b(2) = jpizoom + nimppt(noea+1)
578         ie_b(2) = jpizoom + nimppt(noea+1)+nlcit(noea+1)-3
579         is_b(2) = jpjzoom + njmppt(noea+1)
580         in_b(2) = jpjzoom + njmppt(noea+1)+nlcjt(noea+1)-3
581      ENDIF
582
583      IF(nbondj .eq. 0) THEN
584         iw_b(3) = jpizoom + nimppt(noso+1)
585         ie_b(3) = jpizoom + nimppt(noso+1)+nlcit(noso+1)-3
586         is_b(3) = jpjzoom + njmppt(noso+1)
587         in_b(3) = jpjzoom + njmppt(noso+1)+nlcjt(noso+1)-3
588
589         iw_b(4) = jpizoom + nimppt(nono+1)
590         ie_b(4) = jpizoom + nimppt(nono+1)+nlcit(nono+1)-3
591         is_b(4) = jpjzoom + njmppt(nono+1)
592         in_b(4) = jpjzoom + njmppt(nono+1)+nlcjt(nono+1)-3
593      ELSEIF(nbondj .eq. 1) THEN
594         iw_b(3) = jpizoom + nimppt(noso+1)
595         ie_b(3) = jpizoom + nimppt(noso+1)+nlcit(noso+1)-3
596         is_b(3) = jpjzoom + njmppt(noso+1)
597         in_b(3) = jpjzoom + njmppt(noso+1)+nlcjt(noso+1)-3
598      ELSEIF(nbondj .eq. -1) THEN
599         iw_b(4) = jpizoom + nimppt(nono+1)
600         ie_b(4) = jpizoom + nimppt(nono+1)+nlcit(nono+1)-3
601         is_b(4) = jpjzoom + njmppt(nono+1)
602         in_b(4) = jpjzoom + njmppt(nono+1)+nlcjt(nono+1)-3
603      ENDIF
604
605      DO ib_bdy = 1, nb_bdy
606         DO igrd = 1, jpbgrd
607            icount  = 0
608            icountr = 0
609            idx_bdy(ib_bdy)%nblen(igrd)    = 0
610            idx_bdy(ib_bdy)%nblenrim(igrd) = 0
611            DO ib = 1, nblendta(igrd,ib_bdy)
612               ! check that data is in correct order in file
613               ibm1 = MAX(1,ib-1)
614               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc...
615                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
616                     CALL ctl_stop('bdy_init : ERROR : boundary data in file must be defined in order of distance from edge nbr.', &
617                    'A utility for re-ordering boundary coordinates and data files exists in the TOOLS/OBC directory')
618                  ENDIF   
619               ENDIF
620               ! check if point is in local domain
621               IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
622                  & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in       ) THEN
623                  !
624                  icount = icount  + 1
625                  !
626                  IF( nbrdta(ib,igrd,ib_bdy) == 1 )   icountr = icountr+1
627               ENDIF
628            ENDDO
629            idx_bdy(ib_bdy)%nblenrim(igrd) = icountr !: length of rim boundary data on each proc
630            idx_bdy(ib_bdy)%nblen   (igrd) = icount  !: length of boundary data on each proc       
631         ENDDO  ! igrd
632
633         ! Allocate index arrays for this boundary set
634         !--------------------------------------------
635         ilen1 = MAXVAL(idx_bdy(ib_bdy)%nblen(:))
636         ALLOCATE( idx_bdy(ib_bdy)%nbi(ilen1,jpbgrd) )
637         ALLOCATE( idx_bdy(ib_bdy)%nbj(ilen1,jpbgrd) )
638         ALLOCATE( idx_bdy(ib_bdy)%nbr(ilen1,jpbgrd) )
639         ALLOCATE( idx_bdy(ib_bdy)%nbmap(ilen1,jpbgrd) )
640         ALLOCATE( idx_bdy(ib_bdy)%nbw(ilen1,jpbgrd) )
641         ALLOCATE( idx_bdy(ib_bdy)%flagu(ilen1) )
642         ALLOCATE( idx_bdy(ib_bdy)%flagv(ilen1) )
643
644         ! Dispatch mapping indices and discrete distances on each processor
645         ! -----------------------------------------------------------------
646
647         com_east = 0
648         com_west = 0
649         com_south = 0
650         com_north = 0
651
652         com_east_b = 0
653         com_west_b = 0
654         com_south_b = 0
655         com_north_b = 0
656         DO igrd = 1, jpbgrd
657            icount  = 0
658            ! Loop on rimwidth to ensure outermost points come first in the local arrays.
659            DO ir=1, nn_rimwidth(ib_bdy)
660               DO ib = 1, nblendta(igrd,ib_bdy)
661                  ! check if point is in local domain and equals ir
662                  IF(  nbidta(ib,igrd,ib_bdy) >= iw .AND. nbidta(ib,igrd,ib_bdy) <= ie .AND.   &
663                     & nbjdta(ib,igrd,ib_bdy) >= is .AND. nbjdta(ib,igrd,ib_bdy) <= in .AND.   &
664                     & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
665                     !
666                     icount = icount  + 1
667                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1
668                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1
669                     ! check if point has to be sent
670                     ii = idx_bdy(ib_bdy)%nbi(icount,igrd)
671                     ij = idx_bdy(ib_bdy)%nbj(icount,igrd)
672                     if((com_east .ne. 1) .and. (ii .eq. (nlci-1)) .and. (nbondi .le. 0)) then
673                        com_east = 1
674                     elseif((com_west .ne. 1) .and. (ii .eq. 2) .and. (nbondi .ge. 0) .and. (nbondi .ne. 2)) then
675                        com_west = 1
676                     endif
677                     if((com_south .ne. 1) .and. (ij .eq. 2) .and. (nbondj .ge. 0) .and. (nbondj .ne. 2)) then
678                        com_south = 1
679                     elseif((com_north .ne. 1) .and. (ij .eq. (nlcj-1)) .and. (nbondj .le. 0)) then
680                        com_north = 1
681                     endif
682                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
683                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
684                  ENDIF
685                  ! check if point has to be received from a neighbour
686                  IF(nbondi .eq. 0) THEN
687                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND.   &
688                       & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND.   &
689                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
690                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
691                       if((com_west_b .ne. 1) .and. (ii .eq. (nlcit(nowe+1)-1))) then
692                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
693                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
694                            com_south = 1
695                          elseif((ij .eq. nlcjt(nowe+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
696                            com_north = 1
697                          endif
698                          com_west_b = 1
699                       endif
700                     ENDIF
701                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND.   &
702                       & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND.   &
703                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
704                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
705                       if((com_east_b .ne. 1) .and. (ii .eq. 2)) then
706                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
707                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
708                            com_south = 1
709                          elseif((ij .eq. nlcjt(noea+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
710                            com_north = 1
711                          endif
712                          com_east_b = 1
713                       endif
714                     ENDIF
715                  ELSEIF(nbondi .eq. 1) THEN
716                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(1) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(1) .AND.   &
717                       & nbjdta(ib,igrd,ib_bdy) >= is_b(1) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(1) .AND.   &
718                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
719                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(1)+2
720                       if((com_west_b .ne. 1) .and. (ii .eq. (nlcit(nowe+1)-1))) then
721                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(1)+2
722                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
723                            com_south = 1
724                          elseif((ij .eq. nlcjt(nowe+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
725                            com_north = 1
726                          endif
727                          com_west_b = 1
728                       endif
729                     ENDIF
730                  ELSEIF(nbondi .eq. -1) THEN
731                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(2) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(2) .AND.   &
732                       & nbjdta(ib,igrd,ib_bdy) >= is_b(2) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(2) .AND.   &
733                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
734                       ii = nbidta(ib,igrd,ib_bdy)- iw_b(2)+2
735                       if((com_east_b .ne. 1) .and. (ii .eq. 2)) then
736                          ij = nbjdta(ib,igrd,ib_bdy) - is_b(2)+2
737                          if((ij .eq. 2) .and. (nbondj .eq. 0 .or. nbondj .eq. 1)) then
738                            com_south = 1
739                          elseif((ij .eq. nlcjt(noea+1)-1) .and. (nbondj .eq. 0 .or. nbondj .eq. -1)) then
740                            com_north = 1
741                          endif
742                          com_east_b = 1
743                       endif
744                     ENDIF
745                  ENDIF
746                  IF(nbondj .eq. 0) THEN
747                     IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1 .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
748                       & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
749                       com_north_b = 1 
750                     ENDIF
751                     IF(com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 .OR. nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
752                       & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
753                       com_south_b = 1 
754                     ENDIF
755                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND.   &
756                       & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND.   &
757                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
758                       ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
759                       if((com_south_b .ne. 1) .and. (ij .eq. (nlcjt(noso+1)-1))) then
760                          com_south_b = 1
761                       endif
762                     ENDIF
763                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND.   &
764                       & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND.   &
765                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
766                       ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
767                       if((com_north_b .ne. 1) .and. (ij .eq. 2)) then
768                          com_north_b = 1
769                       endif
770                     ENDIF
771                  ELSEIF(nbondj .eq. 1) THEN
772                     IF(com_south_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(3)-1 .OR. nbidta(ib,igrd,ib_bdy) == ie_b(3)+1) .AND. &
773                       & nbjdta(ib,igrd,ib_bdy) == in_b(3) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
774                       com_south_b = 1 
775                     ENDIF
776                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(3) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(3) .AND.   &
777                       & nbjdta(ib,igrd,ib_bdy) >= is_b(3) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(3) .AND.   &
778                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
779                       ij = nbjdta(ib,igrd,ib_bdy)- is_b(3)+2
780                       if((com_south_b .ne. 1) .and. (ij .eq. (nlcjt(noso+1)-1))) then
781                          com_south_b = 1
782                       endif
783                     ENDIF
784                  ELSEIF(nbondj .eq. -1) THEN
785                     IF(com_north_b .ne. 1 .AND. (nbidta(ib,igrd,ib_bdy) == iw_b(4)-1 .OR. nbidta(ib,igrd,ib_bdy) == ie_b(4)+1) .AND. &
786                       & nbjdta(ib,igrd,ib_bdy) == is_b(4) .AND. nbrdta(ib,igrd,ib_bdy) == ir) THEN
787                       com_north_b = 1 
788                     ENDIF
789                     IF( nbidta(ib,igrd,ib_bdy) >= iw_b(4) .AND. nbidta(ib,igrd,ib_bdy) <= ie_b(4) .AND.   &
790                       & nbjdta(ib,igrd,ib_bdy) >= is_b(4) .AND. nbjdta(ib,igrd,ib_bdy) <= in_b(4) .AND.   &
791                       & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
792                       ij = nbjdta(ib,igrd,ib_bdy)- is_b(4)+2
793                       if((com_north_b .ne. 1) .and. (ij .eq. 2)) then
794                          com_north_b = 1
795                       endif
796                     ENDIF
797                  ENDIF
798               ENDDO
799            ENDDO
800         ENDDO 
801         ! definition of the i- and j- direction local boundaries arrays
802         ! used for sending the boudaries
803         IF((com_east .eq. 1) .and. (com_west .eq. 1)) THEN
804            nbondi_bdy(ib_bdy) = 0
805         ELSEIF ((com_east .eq. 1) .and. (com_west .eq. 0)) THEN
806            nbondi_bdy(ib_bdy) = -1
807         ELSEIF ((com_east .eq. 0) .and. (com_west .eq. 1)) THEN
808            nbondi_bdy(ib_bdy) = 1
809         ENDIF
810
811         IF((com_north .eq. 1) .and. (com_south .eq. 1)) THEN
812            nbondj_bdy(ib_bdy) = 0
813         ELSEIF ((com_north .eq. 1) .and. (com_south .eq. 0)) THEN
814            nbondj_bdy(ib_bdy) = -1
815         ELSEIF ((com_north .eq. 0) .and. (com_south .eq. 1)) THEN
816            nbondj_bdy(ib_bdy) = 1
817         ENDIF
818
819         ! definition of the i- and j- direction local boundaries arrays
820         ! used for receiving the boudaries
821         IF((com_east_b .eq. 1) .and. (com_west_b .eq. 1)) THEN
822            nbondi_bdy_b(ib_bdy) = 0
823         ELSEIF ((com_east_b .eq. 1) .and. (com_west_b .eq. 0)) THEN
824            nbondi_bdy_b(ib_bdy) = -1
825         ELSEIF ((com_east_b .eq. 0) .and. (com_west_b .eq. 1)) THEN
826            nbondi_bdy_b(ib_bdy) = 1
827         ENDIF
828
829         IF((com_north_b .eq. 1) .and. (com_south_b .eq. 1)) THEN
830            nbondj_bdy_b(ib_bdy) = 0
831         ELSEIF ((com_north_b .eq. 1) .and. (com_south_b .eq. 0)) THEN
832            nbondj_bdy_b(ib_bdy) = -1
833         ELSEIF ((com_north_b .eq. 0) .and. (com_south_b .eq. 1)) THEN
834            nbondj_bdy_b(ib_bdy) = 1
835         ENDIF
836
837         ! Compute rim weights for FRS scheme
838         ! ----------------------------------
839         DO igrd = 1, jpbgrd
840            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
841               nbr => idx_bdy(ib_bdy)%nbr(ib,igrd)
842               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( FLOAT( nbr - 1 ) *0.5 )      ! tanh formulation
843!              idx_bdy(ib_bdy)%nbw(ib,igrd) = (FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth))**2      ! quadratic
844!              idx_bdy(ib_bdy)%nbw(ib,igrd) =  FLOAT(nn_rimwidth+1-nbr)/FLOAT(nn_rimwidth)          ! linear
845            END DO
846         END DO
847
848      ENDDO
849
850      ! ------------------------------------------------------
851      ! Initialise masks and find normal/tangential directions
852      ! ------------------------------------------------------
853
854      ! Read global 2D mask at T-points: bdytmask
855      ! -----------------------------------------
856      ! bdytmask = 1  on the computational domain AND on open boundaries
857      !          = 0  elsewhere   
858 
859      IF( cp_cfg == "eel" .AND. jp_cfg == 5 ) THEN          ! EEL configuration at 5km resolution
860         zmask(         :                ,:) = 0.e0
861         zmask(jpizoom+1:jpizoom+jpiglo-2,:) = 1.e0         
862      ELSE IF( ln_mask_file ) THEN
863         CALL iom_open( cn_mask_file, inum )
864         CALL iom_get ( inum, jpdom_data, 'bdy_msk', zmask(:,:) )
865         CALL iom_close( inum )
866      ELSE
867         zmask(:,:) = 1.e0
868      ENDIF
869
870      DO ij = 1, nlcj      ! Save mask over local domain     
871         DO ii = 1, nlci
872            bdytmask(ii,ij) = zmask( mig(ii), mjg(ij) )
873         END DO
874      END DO
875
876      ! Derive mask on U and V grid from mask on T grid
877      bdyumask(:,:) = 0.e0
878      bdyvmask(:,:) = 0.e0
879      DO ij=1, jpjm1
880         DO ii=1, jpim1
881            bdyumask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii+1, ij )
882            bdyvmask(ii,ij)=bdytmask(ii,ij)*bdytmask(ii  ,ij+1) 
883         END DO
884      END DO
885      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )      ! Lateral boundary cond.
886
887
888      ! Mask corrections
889      ! ----------------
890      DO ik = 1, jpkm1
891         DO ij = 1, jpj
892            DO ii = 1, jpi
893               tmask(ii,ij,ik) = tmask(ii,ij,ik) * bdytmask(ii,ij)
894               umask(ii,ij,ik) = umask(ii,ij,ik) * bdyumask(ii,ij)
895               vmask(ii,ij,ik) = vmask(ii,ij,ik) * bdyvmask(ii,ij)
896               bmask(ii,ij)    = bmask(ii,ij)    * bdytmask(ii,ij)
897            END DO     
898         END DO
899      END DO
900
901      DO ik = 1, jpkm1
902         DO ij = 2, jpjm1
903            DO ii = 2, jpim1
904               fmask(ii,ij,ik) = fmask(ii,ij,ik) * bdytmask(ii,ij  ) * bdytmask(ii+1,ij  )   &
905                  &                              * bdytmask(ii,ij+1) * bdytmask(ii+1,ij+1)
906            END DO     
907         END DO
908      END DO
909
910      tmask_i (:,:) = tmask(:,:,1) * tmask_i(:,:)             
911      bdytmask(:,:) = tmask(:,:,1)
912
913      ! bdy masks and bmask are now set to zero on boundary points:
914      igrd = 1       ! In the free surface case, bmask is at T-points
915      DO ib_bdy = 1, nb_bdy
916        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
917          bmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
918        ENDDO
919      ENDDO
920      !
921      igrd = 1
922      DO ib_bdy = 1, nb_bdy
923        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)     
924          bdytmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
925        ENDDO
926      ENDDO
927      !
928      igrd = 2
929      DO ib_bdy = 1, nb_bdy
930        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
931          bdyumask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
932        ENDDO
933      ENDDO
934      !
935      igrd = 3
936      DO ib_bdy = 1, nb_bdy
937        DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
938          bdyvmask(idx_bdy(ib_bdy)%nbi(ib,igrd), idx_bdy(ib_bdy)%nbj(ib,igrd)) = 0.e0
939        ENDDO
940      ENDDO
941
942      ! Lateral boundary conditions
943      CALL lbc_lnk( fmask        , 'F', 1. )   ;   CALL lbc_lnk( bdytmask(:,:), 'T', 1. )
944      CALL lbc_lnk( bdyumask(:,:), 'U', 1. )   ;   CALL lbc_lnk( bdyvmask(:,:), 'V', 1. )
945
946      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
947
948         idx_bdy(ib_bdy)%flagu(:) = 0.e0
949         idx_bdy(ib_bdy)%flagv(:) = 0.e0
950         icount = 0 
951
952         !flagu = -1 : u component is normal to the dynamical boundary but its direction is outward
953         !flagu =  0 : u is tangential
954         !flagu =  1 : u is normal to the boundary and is direction is inward
955 
956         igrd = 2      ! u-component
957         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
958            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
959            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
960            zefl = bdytmask(nbi  ,nbj)
961            zwfl = bdytmask(nbi+1,nbj)
962            IF( zefl + zwfl == 2 ) THEN
963               icount = icount + 1
964            ELSE
965               idx_bdy(ib_bdy)%flagu(ib)=-zefl+zwfl
966            ENDIF
967         END DO
968
969         !flagv = -1 : u component is normal to the dynamical boundary but its direction is outward
970         !flagv =  0 : u is tangential
971         !flagv =  1 : u is normal to the boundary and is direction is inward
972
973         igrd = 3      ! v-component
974         DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd) 
975            nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
976            nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
977            znfl = bdytmask(nbi,nbj  )
978            zsfl = bdytmask(nbi,nbj+1)
979            IF( znfl + zsfl == 2 ) THEN
980               icount = icount + 1
981            ELSE
982               idx_bdy(ib_bdy)%flagv(ib) = -znfl + zsfl
983            END IF
984         END DO
985 
986         IF( icount /= 0 ) THEN
987            IF(lwp) WRITE(numout,*)
988            IF(lwp) WRITE(numout,*) ' E R R O R : Some data velocity points,',   &
989               ' are not boundary points. Check nbi, nbj, indices for boundary set ',ib_bdy
990            IF(lwp) WRITE(numout,*) ' ========== '
991            IF(lwp) WRITE(numout,*)
992            nstop = nstop + 1
993         ENDIF
994   
995      ENDDO
996
997      ! Compute total lateral surface for volume correction:
998      ! ----------------------------------------------------
999      bdysurftot = 0.e0 
1000      IF( ln_vol ) THEN 
1001         igrd = 2      ! Lateral surface at U-points
1002         DO ib_bdy = 1, nb_bdy
1003            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1004               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1005               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1006               flagu => idx_bdy(ib_bdy)%flagu(ib)
1007               bdysurftot = bdysurftot + hu     (nbi  , nbj)                           &
1008                  &                    * e2u    (nbi  , nbj) * ABS( flagu ) &
1009                  &                    * tmask_i(nbi  , nbj)                           &
1010                  &                    * tmask_i(nbi+1, nbj)                   
1011            ENDDO
1012         ENDDO
1013
1014         igrd=3 ! Add lateral surface at V-points
1015         DO ib_bdy = 1, nb_bdy
1016            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
1017               nbi => idx_bdy(ib_bdy)%nbi(ib,igrd)
1018               nbj => idx_bdy(ib_bdy)%nbj(ib,igrd)
1019               flagv => idx_bdy(ib_bdy)%flagv(ib)
1020               bdysurftot = bdysurftot + hv     (nbi, nbj  )                           &
1021                  &                    * e1v    (nbi, nbj  ) * ABS( flagv ) &
1022                  &                    * tmask_i(nbi, nbj  )                           &
1023                  &                    * tmask_i(nbi, nbj+1)
1024            ENDDO
1025         ENDDO
1026         !
1027         IF( lk_mpp )   CALL mpp_sum( bdysurftot )      ! sum over the global domain
1028      END IF   
1029      !
1030      ! Tidy up
1031      !--------
1032      DEALLOCATE(nbidta, nbjdta, nbrdta)
1033
1034      IF( nn_timing == 1 ) CALL timing_stop('bdy_init')
1035
1036   END SUBROUTINE bdy_init
1037
1038#else
1039   !!---------------------------------------------------------------------------------
1040   !!   Dummy module                                   NO open boundaries
1041   !!---------------------------------------------------------------------------------
1042CONTAINS
1043   SUBROUTINE bdy_init      ! Dummy routine
1044   END SUBROUTINE bdy_init
1045#endif
1046
1047   !!=================================================================================
1048END MODULE bdyini
Note: See TracBrowser for help on using the repository browser.