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 NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/BDY – NEMO

source: NEMO/branches/2020/dev_r12563_ASINTER-06_ABL_improvement/src/OCE/BDY/bdyini.F90 @ 14325

Last change on this file since 14325 was 13159, checked in by gsamson, 4 years ago

merge trunk@r13136 into ASINTER-06 branch; pass all SETTE tests; results identical to trunk@r13136; ticket #2419

  • Property svn:keywords set to Id
File size: 89.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.4  !  2012     (J. Chanut) straight open boundary case update
14   !!            3.5  !  2012     (S. Mocavero, I. Epicoco) optimization of BDY communications
15   !!            3.7  !  2016     (T. Lovato) Remove bdy macro, call here init for dta and tides
16   !!----------------------------------------------------------------------
17   !!   bdy_init      : Initialization of unstructured open boundaries
18   !!----------------------------------------------------------------------
19   USE oce            ! ocean dynamics and tracers variables
20   USE dom_oce        ! ocean space and time domain
21   USE sbc_oce , ONLY: nn_ice
22   USE bdy_oce        ! unstructured open boundary conditions
23   USE bdydta         ! open boundary cond. setting   (bdy_dta_init routine)
24   USE bdytides       ! open boundary cond. setting   (bdytide_init routine)
25   USE tide_mod, ONLY: ln_tide ! tidal forcing
26   USE phycst  , ONLY: rday
27   !
28   USE in_out_manager ! I/O units
29   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
30   USE lib_mpp        ! for mpp_sum 
31   USE iom            ! I/O
32
33   IMPLICIT NONE
34   PRIVATE
35
36   PUBLIC   bdy_init    ! routine called in nemo_init
37   PUBLIC   find_neib   ! routine called in bdy_nmn
38
39   INTEGER, PARAMETER ::   jp_nseg = 100   !
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/OCE 4.0 , NEMO Consortium (2018)
48   !! $Id$
49   !! Software governed by the CeCILL license (see ./LICENSE)
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/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,         &
66         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     &
67         &             cn_dyn3d, nn_dyn3d_dta, cn_tra, nn_tra_dta,             &
68         &             ln_tra_dmp, ln_dyn3d_dmp, rn_time_dmp, rn_time_dmp_out, &
69         &             cn_ice, nn_ice_dta,                                     &
70         &             ln_vol, nn_volctl, nn_rimwidth
71         !
72      INTEGER  ::   ios                 ! Local integer output status for namelist read
73      !!----------------------------------------------------------------------
74
75      ! ------------------------
76      ! Read namelist parameters
77      ! ------------------------
78      READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 901)
79901   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy in reference namelist' )
80      ! make sur that all elements of the namelist variables have a default definition from namelist_ref
81      ln_coords_file (2:jp_bdy) = ln_coords_file (1)
82      cn_coords_file (2:jp_bdy) = cn_coords_file (1)
83      cn_dyn2d       (2:jp_bdy) = cn_dyn2d       (1)
84      nn_dyn2d_dta   (2:jp_bdy) = nn_dyn2d_dta   (1)
85      cn_dyn3d       (2:jp_bdy) = cn_dyn3d       (1)
86      nn_dyn3d_dta   (2:jp_bdy) = nn_dyn3d_dta   (1)
87      cn_tra         (2:jp_bdy) = cn_tra         (1)
88      nn_tra_dta     (2:jp_bdy) = nn_tra_dta     (1)   
89      ln_tra_dmp     (2:jp_bdy) = ln_tra_dmp     (1)
90      ln_dyn3d_dmp   (2:jp_bdy) = ln_dyn3d_dmp   (1)
91      rn_time_dmp    (2:jp_bdy) = rn_time_dmp    (1)
92      rn_time_dmp_out(2:jp_bdy) = rn_time_dmp_out(1)
93      cn_ice         (2:jp_bdy) = cn_ice         (1)
94      nn_ice_dta     (2:jp_bdy) = nn_ice_dta     (1)
95      READ  ( numnam_cfg, nambdy, IOSTAT = ios, ERR = 902 )
96902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nambdy in configuration namelist' )
97      IF(lwm) WRITE ( numond, nambdy )
98
99      IF( .NOT. Agrif_Root() ) ln_bdy = .FALSE.   ! forced for Agrif children
100
101      IF( nb_bdy == 0 ) ln_bdy = .FALSE.
102     
103      ! -----------------------------------------
104      ! unstructured open boundaries use control
105      ! -----------------------------------------
106      IF ( ln_bdy ) THEN
107         IF(lwp) WRITE(numout,*)
108         IF(lwp) WRITE(numout,*) 'bdy_init : initialization of open boundaries'
109         IF(lwp) WRITE(numout,*) '~~~~~~~~'
110         !
111         ! Open boundaries definition (arrays and masks)
112         CALL bdy_def
113         IF( ln_meshmask )   CALL bdy_meshwri()
114         !
115         ! Open boundaries initialisation of external data arrays
116         CALL bdy_dta_init
117         !
118         ! Open boundaries initialisation of tidal harmonic forcing
119         IF( ln_tide ) CALL bdytide_init
120         !
121      ELSE
122         IF(lwp) WRITE(numout,*)
123         IF(lwp) WRITE(numout,*) 'bdy_init : open boundaries not used (ln_bdy = F)'
124         IF(lwp) WRITE(numout,*) '~~~~~~~~'
125         !
126      ENDIF
127      !
128   END SUBROUTINE bdy_init
129
130
131   SUBROUTINE bdy_def
132      !!----------------------------------------------------------------------
133      !!                 ***  ROUTINE bdy_init  ***
134      !!         
135      !! ** Purpose :   Definition of unstructured open boundaries.
136      !!
137      !! ** Method  :   Read initialization arrays (mask, indices) to identify
138      !!              an unstructured open boundary
139      !!
140      !! ** Input   :  bdy_init.nc, input file for unstructured open boundaries
141      !!----------------------------------------------------------------------     
142      INTEGER  ::   ib_bdy, ii, ij, igrd, ib, ir, iseg     ! dummy loop indices
143      INTEGER  ::   icount, icountr, icountr0, ibr_max     ! local integers
144      INTEGER  ::   ilen1                                  !   -       -
145      INTEGER  ::   iwe, ies, iso, ino, inum, id_dummy     !   -       -
146      INTEGER  ::   jpbdta                                 !   -       -
147      INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       -
148      INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3           !   -       -
149      INTEGER  ::   iibe, ijbe, iibi, ijbi                 !   -       -
150      INTEGER  ::   flagu, flagv                           ! short cuts
151      INTEGER  ::   nbdyind, nbdybeg, nbdyend
152      INTEGER              , DIMENSION(4)             ::   kdimsz
153      INTEGER              , DIMENSION(jpbgrd,jp_bdy) ::   nblendta          ! Length of index arrays
154      INTEGER,  ALLOCATABLE, DIMENSION(:,:,:)         ::   nbidta, nbjdta    ! Index arrays: i and j indices of bdy dta
155      INTEGER,  ALLOCATABLE, DIMENSION(:,:,:)         ::   nbrdta            ! Discrete distance from rim points
156      CHARACTER(LEN=1)     , DIMENSION(jpbgrd)        ::   cgrid
157      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zz_read                 ! work space for 2D global boundary data
158      REAL(wp), POINTER    , DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields
159      REAL(wp)             , DIMENSION(jpi,jpj) ::   zfmask   ! temporary fmask array excluding coastal boundary condition (shlat)
160      REAL(wp)             , DIMENSION(jpi,jpj) ::   ztmask, zumask, zvmask  ! temporary u/v mask array
161      !!----------------------------------------------------------------------
162      !
163      cgrid = (/'t','u','v'/)
164
165      ! -----------------------------------------
166      ! Check and write out namelist parameters
167      ! -----------------------------------------
168      IF( jperio /= 0 )   CALL ctl_stop( 'bdy_segs: Cyclic or symmetric,',   &
169         &                               ' and general open boundary condition are not compatible' )
170
171      IF(lwp) WRITE(numout,*) 'Number of open boundary sets : ', nb_bdy
172
173      DO ib_bdy = 1,nb_bdy
174
175         IF(lwp) THEN
176            WRITE(numout,*) ' ' 
177            WRITE(numout,*) '------ Open boundary data set ',ib_bdy,' ------' 
178            IF( ln_coords_file(ib_bdy) ) THEN
179               WRITE(numout,*) 'Boundary definition read from file '//TRIM(cn_coords_file(ib_bdy))
180            ELSE
181               WRITE(numout,*) 'Boundary defined in namelist.'
182            ENDIF
183            WRITE(numout,*)
184         ENDIF
185
186         ! barotropic bdy
187         !----------------
188         IF(lwp) THEN
189            WRITE(numout,*) 'Boundary conditions for barotropic solution:  '
190            SELECT CASE( cn_dyn2d(ib_bdy) )                 
191            CASE( 'none' )           ;   WRITE(numout,*) '      no open boundary condition'       
192            CASE( 'frs' )            ;   WRITE(numout,*) '      Flow Relaxation Scheme'
193            CASE( 'flather' )        ;   WRITE(numout,*) '      Flather radiation condition'
194            CASE( 'orlanski' )       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
195            CASE( 'orlanski_npo' )   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
196            CASE DEFAULT             ;   CALL ctl_stop( 'unrecognised value for cn_dyn2d' )
197            END SELECT
198         ENDIF
199
200         dta_bdy(ib_bdy)%lneed_ssh   = cn_dyn2d(ib_bdy) == 'flather'
201         dta_bdy(ib_bdy)%lneed_dyn2d = cn_dyn2d(ib_bdy) /= 'none'
202
203         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn2d ) THEN
204            SELECT CASE( nn_dyn2d_dta(ib_bdy) )                   !
205            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'       
206            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file'
207            CASE( 2 )      ;   WRITE(numout,*) '      tidal harmonic forcing taken from file'
208            CASE( 3 )      ;   WRITE(numout,*) '      boundary data AND tidal harmonic forcing taken from files'
209            CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn2d_dta must be between 0 and 3' )
210            END SELECT
211         ENDIF
212         IF ( dta_bdy(ib_bdy)%lneed_dyn2d .AND. nn_dyn2d_dta(ib_bdy) .GE. 2  .AND. .NOT.ln_tide ) THEN
213            CALL ctl_stop( 'You must activate with ln_tide to add tidal forcing at open boundaries' )
214         ENDIF
215         IF(lwp) WRITE(numout,*)
216
217         ! baroclinic bdy
218         !----------------
219         IF(lwp) THEN
220            WRITE(numout,*) 'Boundary conditions for baroclinic velocities:  '
221            SELECT CASE( cn_dyn3d(ib_bdy) )                 
222            CASE('none')           ;   WRITE(numout,*) '      no open boundary condition'       
223            CASE('frs')            ;   WRITE(numout,*) '      Flow Relaxation Scheme'
224            CASE('specified')      ;   WRITE(numout,*) '      Specified value'
225            CASE('neumann')        ;   WRITE(numout,*) '      Neumann conditions'
226            CASE('zerograd')       ;   WRITE(numout,*) '      Zero gradient for baroclinic velocities'
227            CASE('zero')           ;   WRITE(numout,*) '      Zero baroclinic velocities (runoff case)'
228            CASE('orlanski')       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
229            CASE('orlanski_npo')   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
230            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_dyn3d' )
231            END SELECT
232         ENDIF
233
234         dta_bdy(ib_bdy)%lneed_dyn3d = cn_dyn3d(ib_bdy) == 'frs'      .OR. cn_dyn3d(ib_bdy) == 'specified'   &
235            &                     .OR. cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo'
236
237         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_dyn3d ) THEN
238            SELECT CASE( nn_dyn3d_dta(ib_bdy) )                   !
239            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'       
240            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file'
241            CASE DEFAULT   ;   CALL ctl_stop( 'nn_dyn3d_dta must be 0 or 1' )
242            END SELECT
243         END IF
244
245         IF ( ln_dyn3d_dmp(ib_bdy) ) THEN
246            IF ( cn_dyn3d(ib_bdy) == 'none' ) THEN
247               IF(lwp) WRITE(numout,*) 'No open boundary condition for baroclinic velocities: ln_dyn3d_dmp is set to .false.'
248               ln_dyn3d_dmp(ib_bdy) = .false.
249            ELSEIF ( cn_dyn3d(ib_bdy) == 'frs' ) THEN
250               CALL ctl_stop( 'Use FRS OR relaxation' )
251            ELSE
252               IF(lwp) WRITE(numout,*) '      + baroclinic velocities relaxation zone'
253               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
254               IF(rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
255               dta_bdy(ib_bdy)%lneed_dyn3d = .TRUE.
256            ENDIF
257         ELSE
258            IF(lwp) WRITE(numout,*) '      NO relaxation on baroclinic velocities'
259         ENDIF
260         IF(lwp) WRITE(numout,*)
261
262         !    tra bdy
263         !----------------
264         IF(lwp) THEN
265            WRITE(numout,*) 'Boundary conditions for temperature and salinity:  '
266            SELECT CASE( cn_tra(ib_bdy) )                 
267            CASE('none')           ;   WRITE(numout,*) '      no open boundary condition'       
268            CASE('frs')            ;   WRITE(numout,*) '      Flow Relaxation Scheme'
269            CASE('specified')      ;   WRITE(numout,*) '      Specified value'
270            CASE('neumann')        ;   WRITE(numout,*) '      Neumann conditions'
271            CASE('runoff')         ;   WRITE(numout,*) '      Runoff conditions : Neumann for T and specified to 0.1 for salinity'
272            CASE('orlanski')       ;   WRITE(numout,*) '      Orlanski (fully oblique) radiation condition with adaptive nudging'
273            CASE('orlanski_npo')   ;   WRITE(numout,*) '      Orlanski (NPO) radiation condition with adaptive nudging'
274            CASE DEFAULT           ;   CALL ctl_stop( 'unrecognised value for cn_tra' )
275            END SELECT
276         ENDIF
277
278         dta_bdy(ib_bdy)%lneed_tra = cn_tra(ib_bdy) == 'frs'       .OR. cn_tra(ib_bdy) == 'specified'   &
279            &                   .OR. cn_tra(ib_bdy) == 'orlanski'  .OR. cn_tra(ib_bdy) == 'orlanski_npo' 
280
281         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_tra ) THEN
282            SELECT CASE( nn_tra_dta(ib_bdy) )                   !
283            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'       
284            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file'
285            CASE DEFAULT   ;   CALL ctl_stop( 'nn_tra_dta must be 0 or 1' )
286            END SELECT
287         ENDIF
288
289         IF ( ln_tra_dmp(ib_bdy) ) THEN
290            IF ( cn_tra(ib_bdy) == 'none' ) THEN
291               IF(lwp) WRITE(numout,*) 'No open boundary condition for tracers: ln_tra_dmp is set to .false.'
292               ln_tra_dmp(ib_bdy) = .false.
293            ELSEIF ( cn_tra(ib_bdy) == 'frs' ) THEN
294               CALL ctl_stop( 'Use FRS OR relaxation' )
295            ELSE
296               IF(lwp) WRITE(numout,*) '      + T/S relaxation zone'
297               IF(lwp) WRITE(numout,*) '      Damping time scale: ',rn_time_dmp(ib_bdy),' days'
298               IF(lwp) WRITE(numout,*) '      Outflow damping time scale: ',rn_time_dmp_out(ib_bdy),' days'
299               IF(lwp.AND.rn_time_dmp(ib_bdy)<0) CALL ctl_stop( 'Time scale must be positive' )
300               dta_bdy(ib_bdy)%lneed_tra = .TRUE.
301            ENDIF
302         ELSE
303            IF(lwp) WRITE(numout,*) '      NO T/S relaxation'
304         ENDIF
305         IF(lwp) WRITE(numout,*)
306
307#if defined key_si3
308         IF(lwp) THEN
309            WRITE(numout,*) 'Boundary conditions for sea ice:  '
310            SELECT CASE( cn_ice(ib_bdy) )                 
311            CASE('none')   ;   WRITE(numout,*) '      no open boundary condition'       
312            CASE('frs')    ;   WRITE(numout,*) '      Flow Relaxation Scheme'
313            CASE DEFAULT   ;   CALL ctl_stop( 'unrecognised value for cn_ice' )
314            END SELECT
315         ENDIF
316
317         dta_bdy(ib_bdy)%lneed_ice = cn_ice(ib_bdy) /= 'none'
318
319         IF( dta_bdy(ib_bdy)%lneed_ice .AND. nn_ice /= 2 ) THEN
320            WRITE(ctmp1,*) 'bdy number ', ib_bdy,', needs ice model but nn_ice = ', nn_ice
321            CALL ctl_stop( ctmp1 )
322         ENDIF
323
324         IF( lwp .AND. dta_bdy(ib_bdy)%lneed_ice ) THEN
325            SELECT CASE( nn_ice_dta(ib_bdy) )                   !
326            CASE( 0 )      ;   WRITE(numout,*) '      initial state used for bdy data'       
327            CASE( 1 )      ;   WRITE(numout,*) '      boundary data taken from file'
328            CASE DEFAULT   ;   CALL ctl_stop( 'nn_ice_dta must be 0 or 1' )
329            END SELECT
330         ENDIF
331#else
332         dta_bdy(ib_bdy)%lneed_ice = .FALSE.
333#endif
334         !
335         IF(lwp) WRITE(numout,*)
336         IF(lwp) WRITE(numout,*) '      Width of relaxation zone = ', nn_rimwidth(ib_bdy)
337         IF(lwp) WRITE(numout,*)
338         !
339      END DO   ! nb_bdy
340
341      IF( lwp ) THEN
342         IF( ln_vol ) THEN                     ! check volume conservation (nn_volctl value)
343            WRITE(numout,*) 'Volume correction applied at open boundaries'
344            WRITE(numout,*)
345            SELECT CASE ( nn_volctl )
346            CASE( 1 )      ;   WRITE(numout,*) '      The total volume will be constant'
347            CASE( 0 )      ;   WRITE(numout,*) '      The total volume will vary according to the surface E-P flux'
348            CASE DEFAULT   ;   CALL ctl_stop( 'nn_volctl must be 0 or 1' )
349            END SELECT
350            WRITE(numout,*)
351            !
352            ! sanity check if used with tides       
353            IF( ln_tide ) THEN
354               WRITE(numout,*) ' The total volume correction is not working with tides. '
355               WRITE(numout,*) ' Set ln_vol to .FALSE. '
356               WRITE(numout,*) ' or '
357               WRITE(numout,*) ' equilibriate your bdy input files '
358               CALL ctl_stop( 'The total volume correction is not working with tides.' )
359            END IF
360         ELSE
361            WRITE(numout,*) 'No volume correction applied at open boundaries'
362            WRITE(numout,*)
363         ENDIF
364      ENDIF
365
366      ! -------------------------------------------------
367      ! Initialise indices arrays for open boundaries
368      ! -------------------------------------------------
369
370      nblendta(:,:) = 0
371      nbdysege = 0
372      nbdysegw = 0
373      nbdysegn = 0
374      nbdysegs = 0
375
376      ! Define all boundaries
377      ! ---------------------
378      DO ib_bdy = 1, nb_bdy
379         !
380         IF( .NOT. ln_coords_file(ib_bdy) ) THEN     ! build bdy coordinates with segments defined in namelist
381
382            CALL bdy_read_seg( ib_bdy, nblendta(:,ib_bdy) )
383
384         ELSE                                        ! Read size of arrays in boundary coordinates file.
385           
386            CALL iom_open( cn_coords_file(ib_bdy), inum )
387            DO igrd = 1, jpbgrd
388               id_dummy = iom_varid( inum, 'nbi'//cgrid(igrd), kdimsz=kdimsz ) 
389               nblendta(igrd,ib_bdy) = MAXVAL(kdimsz)
390            END DO
391            CALL iom_close( inum )
392         ENDIF
393         !
394      END DO ! ib_bdy
395
396      ! Now look for crossings in user (namelist) defined open boundary segments:
397      IF( nbdysege > 0 .OR. nbdysegw > 0 .OR. nbdysegn > 0 .OR. nbdysegs > 0)   CALL bdy_ctl_seg
398     
399      ! Allocate arrays
400      !---------------
401      jpbdta = MAXVAL(nblendta(1:jpbgrd,1:nb_bdy))
402      ALLOCATE( nbidta(jpbdta, jpbgrd, nb_bdy), nbjdta(jpbdta, jpbgrd, nb_bdy), nbrdta(jpbdta, jpbgrd, nb_bdy) )
403      nbrdta(:,:,:) = 0   ! initialize nbrdta as it may not be completely defined for each bdy
404     
405      ! Calculate global boundary index arrays or read in from file
406      !------------------------------------------------------------               
407      ! 1. Read global index arrays from boundary coordinates file.
408      DO ib_bdy = 1, nb_bdy
409         !
410         IF( ln_coords_file(ib_bdy) ) THEN
411            !
412            ALLOCATE( zz_read( MAXVAL(nblendta), 1 ) )         
413            CALL iom_open( cn_coords_file(ib_bdy), inum )
414            !
415            DO igrd = 1, jpbgrd
416               CALL iom_get( inum, jpdom_unknown, 'nbi'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
417               DO ii = 1,nblendta(igrd,ib_bdy)
418                  nbidta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
419               END DO
420               CALL iom_get( inum, jpdom_unknown, 'nbj'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
421               DO ii = 1,nblendta(igrd,ib_bdy)
422                  nbjdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
423               END DO
424               CALL iom_get( inum, jpdom_unknown, 'nbr'//cgrid(igrd), zz_read(1:nblendta(igrd,ib_bdy),:) )
425               DO ii = 1,nblendta(igrd,ib_bdy)
426                  nbrdta(ii,igrd,ib_bdy) = NINT( zz_read(ii,1) )
427               END DO
428               !
429               ibr_max = MAXVAL( nbrdta(:,igrd,ib_bdy) )
430               IF(lwp) WRITE(numout,*)
431               IF(lwp) WRITE(numout,*) ' Maximum rimwidth in file is ', ibr_max
432               IF(lwp) WRITE(numout,*) ' nn_rimwidth from namelist is ', nn_rimwidth(ib_bdy)
433               IF (ibr_max < nn_rimwidth(ib_bdy))   &
434                  CALL ctl_stop( 'nn_rimwidth is larger than maximum rimwidth in file',cn_coords_file(ib_bdy) )
435            END DO
436            !
437            CALL iom_close( inum )
438            DEALLOCATE( zz_read )
439            !
440         ENDIF
441         !
442      END DO
443
444      ! 2. Now fill indices corresponding to straight open boundary arrays:
445      CALL bdy_coords_seg( nbidta, nbjdta, nbrdta )
446
447      !  Deal with duplicated points
448      !-----------------------------
449      ! We assign negative indices to duplicated points (to remove them from bdy points to be updated)
450      ! if their distance to the bdy is greater than the other
451      ! If their distance are the same, just keep only one to avoid updating a point twice
452      DO igrd = 1, jpbgrd
453         DO ib_bdy1 = 1, nb_bdy
454            DO ib_bdy2 = 1, nb_bdy
455               IF (ib_bdy1/=ib_bdy2) THEN
456                  DO ib1 = 1, nblendta(igrd,ib_bdy1)
457                     DO ib2 = 1, nblendta(igrd,ib_bdy2)
458                        IF ((nbidta(ib1, igrd, ib_bdy1)==nbidta(ib2, igrd, ib_bdy2)).AND. &
459                           &   (nbjdta(ib1, igrd, ib_bdy1)==nbjdta(ib2, igrd, ib_bdy2))) THEN
460                           !                           IF ((lwp).AND.(igrd==1)) WRITE(numout,*) ' found coincident point ji, jj:', &
461                           !                                                       &              nbidta(ib1, igrd, ib_bdy1),      &
462                           !                                                       &              nbjdta(ib2, igrd, ib_bdy2)
463                           ! keep only points with the lowest distance to boundary:
464                           IF (nbrdta(ib1, igrd, ib_bdy1)<nbrdta(ib2, igrd, ib_bdy2)) THEN
465                              nbidta(ib2, igrd, ib_bdy2) =-ib_bdy2
466                              nbjdta(ib2, igrd, ib_bdy2) =-ib_bdy2
467                           ELSEIF (nbrdta(ib1, igrd, ib_bdy1)>nbrdta(ib2, igrd, ib_bdy2)) THEN
468                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
469                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
470                              ! Arbitrary choice if distances are the same:
471                           ELSE
472                              nbidta(ib1, igrd, ib_bdy1) =-ib_bdy1
473                              nbjdta(ib1, igrd, ib_bdy1) =-ib_bdy1
474                           ENDIF
475                        END IF
476                     END DO
477                  END DO
478               ENDIF
479            END DO
480         END DO
481      END DO
482      !
483      ! Find lenght of boundaries and rim on local mpi domain
484      !------------------------------------------------------
485      !
486      iwe = mig(1)
487      ies = mig(jpi)
488      iso = mjg(1) 
489      ino = mjg(jpj) 
490      !
491      DO ib_bdy = 1, nb_bdy
492         DO igrd = 1, jpbgrd
493            icount   = 0   ! initialization of local bdy length
494            icountr  = 0   ! initialization of local rim 0 and rim 1 bdy length
495            icountr0 = 0   ! initialization of local rim 0 bdy length
496            idx_bdy(ib_bdy)%nblen(igrd)     = 0
497            idx_bdy(ib_bdy)%nblenrim(igrd)  = 0
498            idx_bdy(ib_bdy)%nblenrim0(igrd) = 0
499            DO ib = 1, nblendta(igrd,ib_bdy)
500               ! check that data is in correct order in file
501               IF( ib > 1 ) THEN
502                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ib-1,igrd,ib_bdy) ) THEN
503                     CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', &
504                        &        ' in order of distance from edge nbr A utility for re-ordering ', &
505                        &        ' boundary coordinates and data files exists in the TOOLS/OBC directory')
506                  ENDIF
507               ENDIF
508               ! check if point is in local domain
509               IF(  nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND.   &
510                  & nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino      ) THEN
511                  !
512                  icount = icount + 1
513                  IF( nbrdta(ib,igrd,ib_bdy) == 1 .OR. nbrdta(ib,igrd,ib_bdy) == 0 )   icountr = icountr + 1
514                  IF( nbrdta(ib,igrd,ib_bdy) == 0 )   icountr0 = icountr0 + 1
515               ENDIF
516            END DO
517            idx_bdy(ib_bdy)%nblen    (igrd) = icount   !: length of boundary data on each proc
518            idx_bdy(ib_bdy)%nblenrim (igrd) = icountr  !: length of rim 0 and rim 1 boundary data on each proc   
519            idx_bdy(ib_bdy)%nblenrim0(igrd) = icountr0 !: length of rim 0 boundary data on each proc     
520         END DO   ! igrd
521
522         ! Allocate index arrays for this boundary set
523         !--------------------------------------------
524         ilen1 = MAXVAL( idx_bdy(ib_bdy)%nblen(:) )
525         ALLOCATE( idx_bdy(ib_bdy)%nbi   (ilen1,jpbgrd) ,   &
526            &      idx_bdy(ib_bdy)%nbj   (ilen1,jpbgrd) ,   &
527            &      idx_bdy(ib_bdy)%nbr   (ilen1,jpbgrd) ,   &
528            &      idx_bdy(ib_bdy)%nbd   (ilen1,jpbgrd) ,   &
529            &      idx_bdy(ib_bdy)%nbdout(ilen1,jpbgrd) ,   &
530            &      idx_bdy(ib_bdy)%ntreat(ilen1,jpbgrd) ,   &
531            &      idx_bdy(ib_bdy)%nbmap (ilen1,jpbgrd) ,   &
532            &      idx_bdy(ib_bdy)%nbw   (ilen1,jpbgrd) ,   &
533            &      idx_bdy(ib_bdy)%flagu (ilen1,jpbgrd) ,   &
534            &      idx_bdy(ib_bdy)%flagv (ilen1,jpbgrd) )
535
536         ! Dispatch mapping indices and discrete distances on each processor
537         ! -----------------------------------------------------------------
538         DO igrd = 1, jpbgrd
539            icount  = 0
540            ! Outer loop on rimwidth to ensure outermost points come first in the local arrays.
541            DO ir = 0, nn_rimwidth(ib_bdy)
542               DO ib = 1, nblendta(igrd,ib_bdy)
543                  ! check if point is in local domain and equals ir
544                  IF(  nbidta(ib,igrd,ib_bdy) >= iwe .AND. nbidta(ib,igrd,ib_bdy) <= ies .AND.   &
545                     & nbjdta(ib,igrd,ib_bdy) >= iso .AND. nbjdta(ib,igrd,ib_bdy) <= ino .AND.   &
546                     & nbrdta(ib,igrd,ib_bdy) == ir  ) THEN
547                     !
548                     icount = icount  + 1
549                     idx_bdy(ib_bdy)%nbi(icount,igrd)   = nbidta(ib,igrd,ib_bdy)- mig(1)+1   ! global to local indexes
550                     idx_bdy(ib_bdy)%nbj(icount,igrd)   = nbjdta(ib,igrd,ib_bdy)- mjg(1)+1   ! global to local indexes
551                     idx_bdy(ib_bdy)%nbr(icount,igrd)   = nbrdta(ib,igrd,ib_bdy)
552                     idx_bdy(ib_bdy)%nbmap(icount,igrd) = ib
553                  ENDIF
554               END DO
555            END DO
556         END DO   ! igrd
557
558      END DO   ! ib_bdy
559
560      ! Initialize array indicating communications in bdy
561      ! -------------------------------------------------
562      ALLOCATE( lsend_bdy(nb_bdy,jpbgrd,4,0:1), lrecv_bdy(nb_bdy,jpbgrd,4,0:1) )
563      lsend_bdy(:,:,:,:) = .false.
564      lrecv_bdy(:,:,:,:) = .false. 
565
566      DO ib_bdy = 1, nb_bdy
567         DO igrd = 1, jpbgrd
568            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! only the rim triggers communications, see bdy routines
569               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
570               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
571               IF( ib .LE. idx_bdy(ib_bdy)%nblenrim0(igrd) ) THEN   ;   ir = 0
572               ELSE                                                 ;   ir = 1
573               END IF
574               !
575               ! check if point has to be sent     to   a neighbour
576               ! W neighbour and on the inner left  side
577               IF( ii == 2     .and. (nbondi == 0 .or. nbondi ==  1) )   lsend_bdy(ib_bdy,igrd,1,ir) = .true.
578               ! E neighbour and on the inner right side
579               IF( ii == jpi-1 .and. (nbondi == 0 .or. nbondi == -1) )   lsend_bdy(ib_bdy,igrd,2,ir) = .true.
580               ! S neighbour and on the inner down side
581               IF( ij == 2     .and. (nbondj == 0 .or. nbondj ==  1) )   lsend_bdy(ib_bdy,igrd,3,ir) = .true.
582               ! N neighbour and on the inner up   side
583               IF( ij == jpj-1 .and. (nbondj == 0 .or. nbondj == -1) )   lsend_bdy(ib_bdy,igrd,4,ir) = .true.
584               !
585               ! check if point has to be received from a neighbour
586               ! W neighbour and on the outter left  side
587               IF( ii == 1     .and. (nbondi == 0 .or. nbondi ==  1) )   lrecv_bdy(ib_bdy,igrd,1,ir) = .true.
588               ! E neighbour and on the outter right side
589               IF( ii == jpi   .and. (nbondi == 0 .or. nbondi == -1) )   lrecv_bdy(ib_bdy,igrd,2,ir) = .true.
590               ! S neighbour and on the outter down side
591               IF( ij == 1     .and. (nbondj == 0 .or. nbondj ==  1) )   lrecv_bdy(ib_bdy,igrd,3,ir) = .true.
592               ! N neighbour and on the outter up   side
593               IF( ij == jpj   .and. (nbondj == 0 .or. nbondj == -1) )   lrecv_bdy(ib_bdy,igrd,4,ir) = .true.
594               !
595            END DO
596         END DO  ! igrd
597
598         ! Compute rim weights for FRS scheme
599         ! ----------------------------------
600         DO igrd = 1, jpbgrd
601            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
602               ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same weights
603               idx_bdy(ib_bdy)%nbw(ib,igrd) = 1.- TANH( REAL( ir - 1 ) *0.5 )      ! tanh formulation
604               !               idx_bdy(ib_bdy)%nbw(ib,igrd) = (REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.  ! quadratic
605               !               idx_bdy(ib_bdy)%nbw(ib,igrd) =  REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy))       ! linear
606            END DO
607         END DO
608
609         ! Compute damping coefficients
610         ! ----------------------------
611         DO igrd = 1, jpbgrd
612            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)
613               ir = MAX( 1, idx_bdy(ib_bdy)%nbr(ib,igrd) )   ! both rim 0 and rim 1 have the same damping coefficients
614               idx_bdy(ib_bdy)%nbd(ib,igrd) = 1. / ( rn_time_dmp(ib_bdy) * rday ) & 
615                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
616               idx_bdy(ib_bdy)%nbdout(ib,igrd) = 1. / ( rn_time_dmp_out(ib_bdy) * rday ) & 
617                  & *(REAL(nn_rimwidth(ib_bdy)+1-ir)/REAL(nn_rimwidth(ib_bdy)))**2.   ! quadratic
618            END DO
619         END DO
620
621      END DO ! ib_bdy
622
623      ! ------------------------------------------------------
624      ! Initialise masks and find normal/tangential directions
625      ! ------------------------------------------------------
626
627      ! ------------------------------------------
628      ! handle rim0, do as if rim 1 was free ocean
629      ! ------------------------------------------
630
631      ztmask(:,:) = tmask(:,:,1)   ;   zumask(:,:) = umask(:,:,1)   ;   zvmask(:,:) = vmask(:,:,1)
632      ! For the flagu/flagv calculation below we require a version of fmask without
633      ! the land boundary condition (shlat) included:
634      DO ij = 1, jpjm1
635         DO ii = 1, jpim1
636            zfmask(ii,ij) =  ztmask(ii,ij  ) * ztmask(ii+1,ij  )   &
637               &           * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
638         END DO
639      END DO
640      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. )
641
642      ! Read global 2D mask at T-points: bdytmask
643      ! -----------------------------------------
644      ! bdytmask = 1  on the computational domain AND on open boundaries
645      !          = 0  elsewhere   
646
647      bdytmask(:,:) = ssmask(:,:)
648
649      ! Derive mask on U and V grid from mask on T grid
650      DO ij = 1, jpjm1
651         DO ii = 1, jpim1
652            bdyumask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii+1,ij  )
653            bdyvmask(ii,ij) = bdytmask(ii,ij) * bdytmask(ii  ,ij+1) 
654         END DO
655      END DO
656      CALL lbc_lnk_multi( 'bdyini', bdyumask, 'U', 1., bdyvmask, 'V', 1. )   ! Lateral boundary cond.
657
658      ! bdy masks are now set to zero on rim 0 points:
659      DO ib_bdy = 1, nb_bdy
660         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1)   ! extent of rim 0
661            bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
662         END DO
663         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2)   ! extent of rim 0
664            bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
665         END DO
666         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3)   ! extent of rim 0
667            bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
668         END DO
669      END DO
670
671      CALL bdy_rim_treat( zumask, zvmask, zfmask, .true. )   ! compute flagu, flagv, ntreat on rim 0
672
673      ! ------------------------------------
674      ! handle rim1, do as if rim 0 was land
675      ! ------------------------------------
676     
677      ! z[tuv]mask are now set to zero on rim 0 points:
678      DO ib_bdy = 1, nb_bdy
679         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1)   ! extent of rim 0
680            ztmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
681         END DO
682         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(2)   ! extent of rim 0
683            zumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
684         END DO
685         DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(3)   ! extent of rim 0
686            zvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
687         END DO
688      END DO
689
690      ! Recompute zfmask
691      DO ij = 1, jpjm1
692         DO ii = 1, jpim1
693            zfmask(ii,ij) =  ztmask(ii,ij  ) * ztmask(ii+1,ij  )   &
694               &           * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
695         END DO
696      END DO
697      CALL lbc_lnk( 'bdyini', zfmask, 'F', 1. )
698
699      ! bdy masks are now set to zero on rim1 points:
700      DO ib_bdy = 1, nb_bdy
701         DO ib = idx_bdy(ib_bdy)%nblenrim0(1) + 1,  idx_bdy(ib_bdy)%nblenrim(1)   ! extent of rim 1
702            bdytmask(idx_bdy(ib_bdy)%nbi(ib,1), idx_bdy(ib_bdy)%nbj(ib,1)) = 0._wp
703         END DO
704         DO ib = idx_bdy(ib_bdy)%nblenrim0(2) + 1,  idx_bdy(ib_bdy)%nblenrim(2)   ! extent of rim 1
705            bdyumask(idx_bdy(ib_bdy)%nbi(ib,2), idx_bdy(ib_bdy)%nbj(ib,2)) = 0._wp
706         END DO
707         DO ib = idx_bdy(ib_bdy)%nblenrim0(3) + 1,  idx_bdy(ib_bdy)%nblenrim(3)   ! extent of rim 1
708            bdyvmask(idx_bdy(ib_bdy)%nbi(ib,3), idx_bdy(ib_bdy)%nbj(ib,3)) = 0._wp
709         END DO
710      END DO
711
712      CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. )   ! compute flagu, flagv, ntreat on rim 1
713      !
714      ! Check which boundaries might need communication
715      ALLOCATE( lsend_bdyint(nb_bdy,jpbgrd,4,0:1), lrecv_bdyint(nb_bdy,jpbgrd,4,0:1) )
716      lsend_bdyint(:,:,:,:) = .false.
717      lrecv_bdyint(:,:,:,:) = .false. 
718      ALLOCATE( lsend_bdyext(nb_bdy,jpbgrd,4,0:1), lrecv_bdyext(nb_bdy,jpbgrd,4,0:1) )
719      lsend_bdyext(:,:,:,:) = .false.
720      lrecv_bdyext(:,:,:,:) = .false.
721      !
722      DO igrd = 1, jpbgrd
723         DO ib_bdy = 1, nb_bdy
724            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
725               IF( idx_bdy(ib_bdy)%ntreat(ib,igrd) == -1 ) CYCLE
726               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
727               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
728               ir = idx_bdy(ib_bdy)%nbr(ib,igrd)
729               flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd))
730               flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd))
731               iibe = ii - flagu   ! neighbouring point towards the exterior of the computational domain
732               ijbe = ij - flagv
733               iibi = ii + flagu   ! neighbouring point towards the interior of the computational domain
734               ijbi = ij + flagv
735               CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 )   ! free ocean neighbours
736               !
737               ! search neighbour in the  west/east  direction
738               ! Rim is on the halo and computed ocean is towards exterior of mpi domain 
739               !      <--    (o exterior)     --> 
740               ! (1)  o|x         OR    (2)   x|o
741               !       |___                 ___|
742               IF( iibi == 0     .OR. ii1 == 0     .OR. ii2 == 0     .OR. ii3 == 0     )   lrecv_bdyint(ib_bdy,igrd,1,ir) = .true.
743               IF( iibi == jpi+1 .OR. ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 )   lrecv_bdyint(ib_bdy,igrd,2,ir) = .true. 
744               IF( iibe == 0                                                           )   lrecv_bdyext(ib_bdy,igrd,1,ir) = .true.
745               IF( iibe == jpi+1                                                       )   lrecv_bdyext(ib_bdy,igrd,2,ir) = .true. 
746               ! Check if neighbour has its rim parallel to its mpi subdomain border and located next to its halo
747               ! :¨¨¨¨¨|¨¨-->    |                                             |    <--¨¨|¨¨¨¨¨:
748               ! :     |  x:o    |    neighbour limited by ... would need o    |    o:x  |     :
749               ! :.....|_._:_____|   (1) W neighbour         E neighbour (2)   |_____:_._|.....:
750               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. &
751                  & ( iibi == 3     .OR. ii1 == 3     .OR. ii2 == 3     .OR. ii3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,1,ir)=.true.
752               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. &
753                  & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) )   lsend_bdyint(ib_bdy,igrd,2,ir)=.true.
754               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. iibe == 3     )   lsend_bdyext(ib_bdy,igrd,1,ir)=.true.
755               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2 )   lsend_bdyext(ib_bdy,igrd,2,ir)=.true.
756               !
757               ! search neighbour in the north/south direction   
758               ! Rim is on the halo and computed ocean is towards exterior of mpi domain
759               !(3)   |       |         ^   ___o___     
760               !  |   |___x___|   OR    |  |   x   |
761               !  v       o           (4)  |       |
762               IF( ijbi == 0     .OR. ij1 == 0     .OR. ij2 == 0     .OR. ij3 == 0     )   lrecv_bdyint(ib_bdy,igrd,3,ir) = .true.
763               IF( ijbi == jpj+1 .OR. ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 )   lrecv_bdyint(ib_bdy,igrd,4,ir) = .true.
764               IF( ijbe == 0                                                           )   lrecv_bdyext(ib_bdy,igrd,3,ir) = .true.
765               IF( ijbe == jpj+1                                                       )   lrecv_bdyext(ib_bdy,igrd,4,ir) = .true.
766               ! Check if neighbour has its rim parallel to its mpi subdomain     _________  border and next to its halo
767               !   ^  |    o    |                                                :         :
768               !   |  |¨¨¨¨x¨¨¨¨|   neighbour limited by ... would need o     |  |....x....|
769               !      :_________:  (3) S neighbour          N neighbour (4)   v  |    o    |   
770               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. &
771                  & ( ijbi == 3     .OR. ij1 == 3     .OR. ij2 == 3     .OR. ij3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,3,ir)=.true.
772               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. &
773                  & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) )   lsend_bdyint(ib_bdy,igrd,4,ir)=.true.
774               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. ijbe == 3     )   lsend_bdyext(ib_bdy,igrd,3,ir)=.true.
775               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2 )   lsend_bdyext(ib_bdy,igrd,4,ir)=.true.
776            END DO
777         END DO
778      END DO
779
780      DO ib_bdy = 1,nb_bdy
781         IF(  cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. &
782            & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. &
783            & cn_tra(ib_bdy)   == 'orlanski' .OR. cn_tra(ib_bdy)   == 'orlanski_npo'      ) THEN
784            DO igrd = 1, jpbgrd
785               DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
786                  ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
787                  ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
788                  IF(  mig(ii) > 2 .AND. mig(ii) < jpiglo-2 .AND. mjg(ij) > 2 .AND. mjg(ij) < jpjglo-2  ) THEN
789                     WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain'
790                     CALL ctl_stop( ctmp1 )
791                  END IF
792               END DO
793            END DO
794         END IF
795      END DO
796      !
797      DEALLOCATE( nbidta, nbjdta, nbrdta )
798      !
799   END SUBROUTINE bdy_def
800
801
802   SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 )
803      !!----------------------------------------------------------------------
804      !!                 ***  ROUTINE bdy_rim_treat  ***
805      !!
806      !! ** Purpose :   Initialize structures ( flagu, flagv, ntreat ) indicating how rim points
807      !!                  are to be handled in the boundary condition treatment
808      !!
809      !! ** Method  :   - to handle rim 0 zmasks must indicate ocean points      (set at one on rim 0 and rim 1 and interior)
810      !!                            and bdymasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
811      !!                    (as if rim 1 was free ocean)
812      !!                - to handle rim 1 zmasks must be set at 0 on rim 0       (set at one on rim 1 and interior)
813      !!                            and bdymasks must indicate free ocean points (set at one on interior)
814      !!                    (as if rim 0 was land)
815      !!                - we can then check in which direction the interior of the computational domain is with the difference
816      !!                         mask array values on both sides to compute flagu and flagv
817      !!                - and look at the ocean neighbours to compute ntreat
818      !!----------------------------------------------------------------------
819      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pfmask   ! temporary fmask excluding coastal boundary condition (shlat)
820      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pumask, pvmask   ! temporary t/u/v mask array
821      LOGICAL                             , INTENT (in   ) :: lrim0    ! .true. -> rim 0   .false. -> rim 1
822      INTEGER  ::   ib_bdy, ii, ij, igrd, ib, icount       ! dummy loop indices
823      INTEGER  ::   i_offset, j_offset, inn                ! local integer
824      INTEGER  ::   ibeg, iend                             ! local integer
825      LOGICAL  ::   llnon, llson, llean, llwen             ! local logicals indicating the presence of a ocean neighbour
826      REAL(wp), POINTER, DIMENSION(:,:)       ::   zmask   ! pointer to 2D mask fields
827      REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
828      CHARACTER(LEN=1), DIMENSION(jpbgrd)     ::   cgrid
829      REAL(wp)        , DIMENSION(jpi,jpj)    ::   ztmp
830      !!----------------------------------------------------------------------
831
832      cgrid = (/'t','u','v'/)
833
834      DO ib_bdy = 1, nb_bdy       ! Indices and directions of rim velocity components
835
836         ! Calculate relationship of U direction to the local orientation of the boundary
837         ! flagu = -1 : u component is normal to the dynamical boundary and its direction is outward
838         ! flagu =  0 : u is tangential
839         ! flagu =  1 : u is normal to the boundary and is direction is inward
840         DO igrd = 1, jpbgrd 
841            SELECT CASE( igrd )
842               CASE( 1 )   ;   zmask => pumask     ;   i_offset = 0
843               CASE( 2 )   ;   zmask => bdytmask   ;   i_offset = 1
844               CASE( 3 )   ;   zmask => pfmask     ;   i_offset = 0
845            END SELECT
846            icount = 0
847            ztmp(:,:) = -999._wp
848            IF( lrim0 ) THEN   ! extent of rim 0
849               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
850            ELSE               ! extent of rim 1
851               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd)
852            END IF
853            DO ib = ibeg, iend 
854               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
855               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
856               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE
857               zwfl = zmask(ii+i_offset-1,ij)
858               zefl = zmask(ii+i_offset  ,ij)
859               ! This error check only works if you are using the bdyXmask arrays
860               IF( i_offset == 1 .and. zefl + zwfl == 2. ) THEN
861                  icount = icount + 1
862                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
863               ELSE
864                  ztmp(ii,ij) = -zwfl + zefl
865               ENDIF
866            END DO
867            IF( icount /= 0 ) THEN
868               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   &
869                  ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
870               CALL ctl_stop( ctmp1 )
871            ENDIF
872            SELECT CASE( igrd )
873               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
874               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 
875               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 
876            END SELECT
877            DO ib = ibeg, iend
878               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
879               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
880               idx_bdy(ib_bdy)%flagu(ib,igrd) = ztmp(ii,ij)
881            END DO
882         END DO
883
884         ! Calculate relationship of V direction to the local orientation of the boundary
885         ! flagv = -1 : v component is normal to the dynamical boundary but its direction is outward
886         ! flagv =  0 : v is tangential
887         ! flagv =  1 : v is normal to the boundary and is direction is inward
888         DO igrd = 1, jpbgrd 
889            SELECT CASE( igrd )
890               CASE( 1 )   ;   zmask => pvmask     ;   j_offset = 0
891               CASE( 2 )   ;   zmask => pfmask     ;   j_offset = 0
892               CASE( 3 )   ;   zmask => bdytmask   ;   j_offset = 1
893            END SELECT
894            icount = 0
895            ztmp(:,:) = -999._wp
896            IF( lrim0 ) THEN   ! extent of rim 0
897               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
898            ELSE               ! extent of rim 1
899               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd)
900            END IF
901            DO ib = ibeg, iend
902               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
903               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
904               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE
905               zsfl = zmask(ii,ij+j_offset-1)
906               znfl = zmask(ii,ij+j_offset  )
907               ! This error check only works if you are using the bdyXmask arrays
908               IF( j_offset == 1 .and. znfl + zsfl == 2. ) THEN
909                  IF(lwp) WRITE(numout,*) 'Problem with igrd = ',igrd,' at (global) nbi, nbj : ',mig(ii),mjg(ij)
910                  icount = icount + 1
911               ELSE
912                  ztmp(ii,ij) = -zsfl + znfl
913               END IF
914            END DO
915            IF( icount /= 0 ) THEN
916               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   &
917                  ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
918               CALL ctl_stop( ctmp1 )
919            ENDIF
920            SELECT CASE( igrd )
921               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
922               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 
923               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 
924            END SELECT
925            DO ib = ibeg, iend
926               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
927               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
928               idx_bdy(ib_bdy)%flagv(ib,igrd) = ztmp(ii,ij)
929            END DO
930         END DO
931         !
932      END DO ! ib_bdy
933     
934      DO ib_bdy = 1, nb_bdy
935         DO igrd = 1, jpbgrd
936            SELECT CASE( igrd )
937               CASE( 1 )   ;   zmask => bdytmask 
938               CASE( 2 )   ;   zmask => bdyumask 
939               CASE( 3 )   ;   zmask => bdyvmask 
940            END SELECT
941            ztmp(:,:) = -999._wp
942            IF( lrim0 ) THEN   ! extent of rim 0
943               ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
944            ELSE               ! extent of rim 1
945               ibeg = idx_bdy(ib_bdy)%nblenrim0(igrd) + 1   ;   iend = idx_bdy(ib_bdy)%nblenrim(igrd)
946            END IF
947            DO ib = ibeg, iend
948               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
949               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
950               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE
951               llnon = zmask(ii  ,ij+1) == 1. 
952               llson = zmask(ii  ,ij-1) == 1. 
953               llean = zmask(ii+1,ij  ) == 1. 
954               llwen = zmask(ii-1,ij  ) == 1. 
955               inn  = COUNT( (/ llnon, llson, llean, llwen /) )
956               IF( inn == 0 ) THEN   ! no neighbours -> interior of a corner  or  cluster of rim points
957                  !               !              !     _____     !     _____    !    __     __
958                  !  1 |   o      !  2  o   |    !  3 | x        !  4     x |   !      |   |   -> error
959                  !    |_x_ _     !    _ _x_|    !    |   o      !      o   |   !      |x_x|
960                  IF(     zmask(ii+1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 1.
961                  ELSEIF( zmask(ii-1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 2.
962                  ELSEIF( zmask(ii+1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 3.
963                  ELSEIF( zmask(ii-1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 4.
964                  ELSE                                    ;   ztmp(ii,ij) = -1.
965                     WRITE(ctmp1,*) 'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
966                       ' on boundary set ', ib_bdy, ' has no free ocean neighbour'
967                     IF( lrim0 ) THEN
968                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 0 points.'
969                     ELSE
970                        WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.'
971                     END IF
972                     CALL ctl_warn( ctmp1, ctmp2 )
973                  END IF
974               END IF
975               IF( inn == 1 ) THEN   ! middle of linear bdy  or incomplete corner  ! ___ o
976                  !    |         !         |   !      o     !    ______            !    |x___
977                  ! 5  | x o     ! 6   o x |   ! 7  __x__   ! 8    x
978                  !    |         !         |   !            !      o
979                  IF( llean )   ztmp(ii,ij) = 5.
980                  IF( llwen )   ztmp(ii,ij) = 6.
981                  IF( llnon )   ztmp(ii,ij) = 7.
982                  IF( llson )   ztmp(ii,ij) = 8.
983               END IF
984               IF( inn == 2 ) THEN   ! exterior of a corner
985                  !        o      !        o      !    _____|       !       |_____ 
986                  !  9 ____x o    ! 10   o x___   ! 11     x o      ! 12   o x     
987                  !         |     !       |       !        o        !        o
988                  IF( llnon .AND. llean )   ztmp(ii,ij) =  9.
989                  IF( llnon .AND. llwen )   ztmp(ii,ij) = 10.
990                  IF( llson .AND. llean )   ztmp(ii,ij) = 11.
991                  IF( llson .AND. llwen )   ztmp(ii,ij) = 12.
992               END IF
993               IF( inn == 3 ) THEN   ! 3 neighbours     __   __
994                  !    |_  o      !        o  _|  !       |_|     !       o         
995                  ! 13  _| x o    ! 14   o x |_   ! 15   o x o    ! 16  o x o       
996                  !    |   o      !        o   |  !        o      !    __|¨|__   
997                  IF( llnon .AND. llean .AND. llson )   ztmp(ii,ij) = 13.
998                  IF( llnon .AND. llwen .AND. llson )   ztmp(ii,ij) = 14.
999                  IF( llwen .AND. llson .AND. llean )   ztmp(ii,ij) = 15.
1000                  IF( llwen .AND. llnon .AND. llean )   ztmp(ii,ij) = 16.
1001               END IF
1002               IF( inn == 4 ) THEN
1003                  WRITE(ctmp1,*)  'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
1004                       ' on boundary set ', ib_bdy, ' have 4 neighbours'
1005                  CALL ctl_stop( ctmp1 )
1006               END IF
1007            END DO
1008            SELECT CASE( igrd )
1009               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
1010               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 
1011               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 
1012            END SELECT
1013            DO ib = ibeg, iend
1014               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1015               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1016               idx_bdy(ib_bdy)%ntreat(ib,igrd) = NINT(ztmp(ii,ij))
1017            END DO
1018         END DO
1019      END DO
1020
1021    END SUBROUTINE bdy_rim_treat
1022
1023   
1024    SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 )
1025      !!----------------------------------------------------------------------
1026      !!                 ***  ROUTINE find_neib  ***
1027      !!
1028      !! ** Purpose :   get ii1, ij1, ii2, ij2, ii3, ij3, the indices of
1029      !!               the free ocean neighbours of (ii,ij) for bdy treatment
1030      !!
1031      !! ** Method  :  use itreat input to select a case
1032      !!               N.B. ntreat is defined for all bdy points in routine bdy_rim_treat
1033      !!
1034      !!----------------------------------------------------------------------
1035      INTEGER, INTENT(in   )      ::   ii, ij, itreat
1036      INTEGER, INTENT(  out)      ::   ii1, ij1, ii2, ij2, ii3, ij3
1037      !!----------------------------------------------------------------------
1038      SELECT CASE( itreat )   ! points that will be used by bdy routines, -1 will be discarded
1039         !               !               !     _____     !     _____     
1040         !  1 |   o      !  2  o   |     !  3 | x        !  4     x |   
1041         !    |_x_ _     !    _ _x_|     !    |   o      !      o   |
1042      CASE( 1 )    ;   ii1 = ii+1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1043      CASE( 2 )    ;   ii1 = ii-1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1044      CASE( 3 )    ;   ii1 = ii+1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1045      CASE( 4 )    ;   ii1 = ii-1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1046         !    |          !         |     !      o        !    ______                   ! or incomplete corner
1047         ! 5  | x o      ! 6   o x |     ! 7  __x__      ! 8    x                      !  7  ____ o
1048         !    |          !         |     !               !      o                      !         |x___
1049      CASE( 5 )    ;   ii1 = ii+1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1050      CASE( 6 )    ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1051      CASE( 7 )    ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1052      CASE( 8 )    ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
1053         !        o      !        o      !    _____|     !       |_____ 
1054         !  9 ____x o    ! 10   o x___   ! 11     x o    ! 12   o x     
1055         !         |     !       |       !        o      !        o     
1056      CASE( 9  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
1057      CASE( 10 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1058      CASE( 11 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1059      CASE( 12 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
1060         !    |_  o      !        o  _|  !     ¨¨|_|¨¨   !       o         
1061         ! 13  _| x o    !  14  o x |_   !  15  o x o    ! 16  o x o       
1062         !    |   o      !        o   |  !        o      !    __|¨|__
1063      CASE( 13 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1   
1064      CASE( 14 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1 
1065      CASE( 15 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij-1   ;   ii3 = ii+1   ;   ij3 = ij   
1066      CASE( 16 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij
1067      END SELECT
1068   END SUBROUTINE find_neib
1069   
1070
1071   SUBROUTINE bdy_read_seg( kb_bdy, knblendta ) 
1072      !!----------------------------------------------------------------------
1073      !!                 ***  ROUTINE bdy_coords_seg  ***
1074      !!
1075      !! ** Purpose :  build bdy coordinates with segments defined in namelist
1076      !!
1077      !! ** Method  :  read namelist nambdy_index blocks
1078      !!
1079      !!----------------------------------------------------------------------
1080      INTEGER                   , INTENT (in   ) ::   kb_bdy           ! bdy number
1081      INTEGER, DIMENSION(jpbgrd), INTENT (  out) ::   knblendta        ! length of index arrays
1082      !!
1083      INTEGER          ::   ios                 ! Local integer output status for namelist read
1084      INTEGER          ::   nbdyind, nbdybeg, nbdyend
1085      INTEGER          ::   nbdy_count, nbdy_rdstart, nbdy_loc
1086      CHARACTER(LEN=1) ::   ctypebdy   !     -        -
1087      CHARACTER(LEN=50)::   cerrmsg    !     -        -
1088      NAMELIST/nambdy_index/ ctypebdy, nbdyind, nbdybeg, nbdyend
1089      !!----------------------------------------------------------------------
1090      ! Need to support possibility of reading more than one nambdy_index from
1091      ! the namelist_cfg internal file.
1092      ! Do this by finding the kb_bdy'th occurence of nambdy_index in the
1093      ! character buffer as the starting point.
1094      nbdy_rdstart = 1
1095      DO nbdy_count = 1, kb_bdy
1096       nbdy_loc = INDEX( numnam_cfg( nbdy_rdstart: ), 'nambdy_index' )
1097       IF( nbdy_loc .GT. 0 ) THEN
1098          nbdy_rdstart = nbdy_rdstart + nbdy_loc
1099       ELSE
1100          WRITE(cerrmsg,'(A,I4,A)') 'Error: entry number ',kb_bdy,' of nambdy_index not found'
1101          ios = -1
1102          CALL ctl_nam ( ios , cerrmsg )
1103       ENDIF
1104      END DO
1105      nbdy_rdstart = MAX( 1, nbdy_rdstart - 2 )
1106      READ  ( numnam_cfg( nbdy_rdstart: ), nambdy_index, IOSTAT = ios, ERR = 904)
1107904   IF( ios /= 0 )   CALL ctl_nam ( ios , 'nambdy_index in configuration namelist' )
1108      IF(lwm) WRITE ( numond, nambdy_index )
1109     
1110      SELECT CASE ( TRIM(ctypebdy) )
1111      CASE( 'N' )
1112         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1113            nbdyind  = jpjglo - 2  ! set boundary to whole side of model domain.
1114            nbdybeg  = 2
1115            nbdyend  = jpiglo - 1
1116         ENDIF
1117         nbdysegn = nbdysegn + 1
1118         npckgn(nbdysegn) = kb_bdy ! Save bdy package number
1119         jpjnob(nbdysegn) = nbdyind
1120         jpindt(nbdysegn) = nbdybeg
1121         jpinft(nbdysegn) = nbdyend
1122         !
1123      CASE( 'S' )
1124         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1125            nbdyind  = 2           ! set boundary to whole side of model domain.
1126            nbdybeg  = 2
1127            nbdyend  = jpiglo - 1
1128         ENDIF
1129         nbdysegs = nbdysegs + 1
1130         npckgs(nbdysegs) = kb_bdy ! Save bdy package number
1131         jpjsob(nbdysegs) = nbdyind
1132         jpisdt(nbdysegs) = nbdybeg
1133         jpisft(nbdysegs) = nbdyend
1134         !
1135      CASE( 'E' )
1136         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1137            nbdyind  = jpiglo - 2  ! set boundary to whole side of model domain.
1138            nbdybeg  = 2
1139            nbdyend  = jpjglo - 1
1140         ENDIF
1141         nbdysege = nbdysege + 1 
1142         npckge(nbdysege) = kb_bdy ! Save bdy package number
1143         jpieob(nbdysege) = nbdyind
1144         jpjedt(nbdysege) = nbdybeg
1145         jpjeft(nbdysege) = nbdyend
1146         !
1147      CASE( 'W' )
1148         IF( nbdyind == -1 ) THEN  ! Automatic boundary definition: if nbdysegX = -1
1149            nbdyind  = 2           ! set boundary to whole side of model domain.
1150            nbdybeg  = 2
1151            nbdyend  = jpjglo - 1
1152         ENDIF
1153         nbdysegw = nbdysegw + 1
1154         npckgw(nbdysegw) = kb_bdy ! Save bdy package number
1155         jpiwob(nbdysegw) = nbdyind
1156         jpjwdt(nbdysegw) = nbdybeg
1157         jpjwft(nbdysegw) = nbdyend
1158         !
1159      CASE DEFAULT   ;   CALL ctl_stop( 'ctypebdy must be N, S, E or W' )
1160      END SELECT
1161     
1162      ! For simplicity we assume that in case of straight bdy, arrays have the same length
1163      ! (even if it is true that last tangential velocity points
1164      ! are useless). This simplifies a little bit boundary data format (and agrees with format
1165      ! used so far in obc package)
1166     
1167      knblendta(1:jpbgrd) =  (nbdyend - nbdybeg + 1) * nn_rimwidth(kb_bdy)
1168     
1169   END SUBROUTINE bdy_read_seg
1170
1171   
1172   SUBROUTINE bdy_ctl_seg
1173      !!----------------------------------------------------------------------
1174      !!                 ***  ROUTINE bdy_ctl_seg  ***
1175      !!
1176      !! ** Purpose :   Check straight open boundary segments location
1177      !!
1178      !! ** Method  :   - Look for open boundary corners
1179      !!                - Check that segments start or end on land
1180      !!----------------------------------------------------------------------
1181      INTEGER  ::   ib, ib1, ib2, ji ,jj, itest 
1182      INTEGER, DIMENSION(jp_nseg,2) :: icorne, icornw, icornn, icorns 
1183      REAL(wp), DIMENSION(2) ::   ztestmask
1184      !!----------------------------------------------------------------------
1185      !
1186      IF (lwp) WRITE(numout,*) ' '
1187      IF (lwp) WRITE(numout,*) 'bdy_ctl_seg: Check analytical segments'
1188      IF (lwp) WRITE(numout,*) '~~~~~~~~~~~~'
1189      !
1190      IF(lwp) WRITE(numout,*) 'Number of east  segments     : ', nbdysege
1191      IF(lwp) WRITE(numout,*) 'Number of west  segments     : ', nbdysegw
1192      IF(lwp) WRITE(numout,*) 'Number of north segments     : ', nbdysegn
1193      IF(lwp) WRITE(numout,*) 'Number of south segments     : ', nbdysegs
1194      ! 1. Check bounds
1195      !----------------
1196      DO ib = 1, nbdysegn
1197         IF (lwp) WRITE(numout,*) '**check north seg bounds pckg: ', npckgn(ib)
1198         IF ((jpjnob(ib).ge.jpjglo-1).or.& 
1199            &(jpjnob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1200         IF (jpindt(ib).ge.jpinft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1201         IF (jpindt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1202         IF (jpinft(ib).gt.jpiglo)     CALL ctl_stop( 'End index out of domain' )
1203      END DO
1204      !
1205      DO ib = 1, nbdysegs
1206         IF (lwp) WRITE(numout,*) '**check south seg bounds pckg: ', npckgs(ib)
1207         IF ((jpjsob(ib).ge.jpjglo-1).or.& 
1208            &(jpjsob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1209         IF (jpisdt(ib).ge.jpisft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1210         IF (jpisdt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1211         IF (jpisft(ib).gt.jpiglo)     CALL ctl_stop( 'End index out of domain' )
1212      END DO
1213      !
1214      DO ib = 1, nbdysege
1215         IF (lwp) WRITE(numout,*) '**check east  seg bounds pckg: ', npckge(ib)
1216         IF ((jpieob(ib).ge.jpiglo-1).or.& 
1217            &(jpieob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1218         IF (jpjedt(ib).ge.jpjeft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1219         IF (jpjedt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1220         IF (jpjeft(ib).gt.jpjglo)     CALL ctl_stop( 'End index out of domain' )
1221      END DO
1222      !
1223      DO ib = 1, nbdysegw
1224         IF (lwp) WRITE(numout,*) '**check west  seg bounds pckg: ', npckgw(ib)
1225         IF ((jpiwob(ib).ge.jpiglo-1).or.& 
1226            &(jpiwob(ib).le.1))        CALL ctl_stop( 'nbdyind out of domain' )
1227         IF (jpjwdt(ib).ge.jpjwft(ib)) CALL ctl_stop( 'Bdy start index is greater than end index' )
1228         IF (jpjwdt(ib).lt.1     )     CALL ctl_stop( 'Start index out of domain' )
1229         IF (jpjwft(ib).gt.jpjglo)     CALL ctl_stop( 'End index out of domain' )
1230      ENDDO
1231      !
1232      !     
1233      ! 2. Look for segment crossings
1234      !------------------------------
1235      IF (lwp) WRITE(numout,*) '**Look for segments corners  :'
1236      !
1237      itest = 0 ! corner number
1238      !
1239      ! flag to detect if start or end of open boundary belongs to a corner
1240      ! if not (=0), it must be on land.
1241      ! if a corner is detected, save bdy package number for further tests
1242      icorne(:,:)=0. ; icornw(:,:)=0. ; icornn(:,:)=0. ; icorns(:,:)=0.
1243      ! South/West crossings
1244      IF ((nbdysegw > 0).AND.(nbdysegs > 0)) THEN
1245         DO ib1 = 1, nbdysegw       
1246            DO ib2 = 1, nbdysegs
1247               IF (( jpisdt(ib2)<=jpiwob(ib1)).AND. &
1248                &  ( jpisft(ib2)>=jpiwob(ib1)).AND. &
1249                &  ( jpjwdt(ib1)<=jpjsob(ib2)).AND. &
1250                &  ( jpjwft(ib1)>=jpjsob(ib2))) THEN
1251                  IF ((jpjwdt(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpiwob(ib1))) THEN 
1252                     ! We have a possible South-West corner                     
1253!                     WRITE(numout,*) ' Found a South-West corner at (i,j): ', jpisdt(ib2), jpjwdt(ib1)
1254!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgs(ib2)
1255                     icornw(ib1,1) = npckgs(ib2)
1256                     icorns(ib2,1) = npckgw(ib1)
1257                  ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
1258                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1259                        &                                     jpisft(ib2), jpjwft(ib1)
1260                     WRITE(ctmp2,*) ' Not allowed yet'
1261                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), & 
1262                        &                            ' and South segment: ',npckgs(ib2)
1263                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1264                  ELSE
1265                     WRITE(ctmp1,*) ' Check South and West Open boundary indices'
1266                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , &
1267                        &                            ' and South segment: ',npckgs(ib2)
1268                     CALL ctl_stop( ctmp1, ctmp2 )
1269                  END IF
1270               END IF
1271            END DO
1272         END DO
1273      END IF
1274      !
1275      ! South/East crossings
1276      IF ((nbdysege > 0).AND.(nbdysegs > 0)) THEN
1277         DO ib1 = 1, nbdysege
1278            DO ib2 = 1, nbdysegs
1279               IF (( jpisdt(ib2)<=jpieob(ib1)+1).AND. &
1280                &  ( jpisft(ib2)>=jpieob(ib1)+1).AND. &
1281                &  ( jpjedt(ib1)<=jpjsob(ib2)  ).AND. &
1282                &  ( jpjeft(ib1)>=jpjsob(ib2)  )) THEN
1283                  IF ((jpjedt(ib1)==jpjsob(ib2)).AND.(jpisft(ib2)==jpieob(ib1)+1)) THEN
1284                     ! We have a possible South-East corner
1285!                     WRITE(numout,*) ' Found a South-East corner at (i,j): ', jpisft(ib2), jpjedt(ib1)
1286!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgs(ib2)
1287                     icorne(ib1,1) = npckgs(ib2)
1288                     icorns(ib2,2) = npckge(ib1)
1289                  ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
1290                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1291                        &                                     jpisdt(ib2), jpjeft(ib1)
1292                     WRITE(ctmp2,*) ' Not allowed yet'
1293                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1294                        &                            ' and South segment: ',npckgs(ib2)
1295                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1296                  ELSE
1297                     WRITE(ctmp1,*) ' Check South and East Open boundary indices'
1298                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1299                     &                               ' and South segment: ',npckgs(ib2)
1300                     CALL ctl_stop( ctmp1, ctmp2 )
1301                  END IF
1302               END IF
1303            END DO
1304         END DO
1305      END IF
1306      !
1307      ! North/West crossings
1308      IF ((nbdysegn > 0).AND.(nbdysegw > 0)) THEN
1309         DO ib1 = 1, nbdysegw       
1310            DO ib2 = 1, nbdysegn
1311               IF (( jpindt(ib2)<=jpiwob(ib1)  ).AND. &
1312                &  ( jpinft(ib2)>=jpiwob(ib1)  ).AND. &
1313                &  ( jpjwdt(ib1)<=jpjnob(ib2)+1).AND. &
1314                &  ( jpjwft(ib1)>=jpjnob(ib2)+1)) THEN
1315                  IF ((jpjwft(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpiwob(ib1))) THEN
1316                     ! We have a possible North-West corner
1317!                     WRITE(numout,*) ' Found a North-West corner at (i,j): ', jpindt(ib2), jpjwft(ib1)
1318!                     WRITE(numout,*) ' between segments: ', npckgw(ib1), npckgn(ib2)
1319                     icornw(ib1,2) = npckgn(ib2)
1320                     icornn(ib2,1) = npckgw(ib1)
1321                  ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
1322                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1323                     &                                     jpinft(ib2), jpjwdt(ib1)
1324                     WRITE(ctmp2,*) ' Not allowed yet'
1325                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1326                     &                               ' and North segment: ',npckgn(ib2)
1327                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1328                  ELSE
1329                     WRITE(ctmp1,*) ' Check North and West Open boundary indices'
1330                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), &
1331                     &                               ' and North segment: ',npckgn(ib2)
1332                     CALL ctl_stop( ctmp1, ctmp2 )
1333                  END IF
1334               END IF
1335            END DO
1336         END DO
1337      END IF
1338      !
1339      ! North/East crossings
1340      IF ((nbdysegn > 0).AND.(nbdysege > 0)) THEN
1341         DO ib1 = 1, nbdysege       
1342            DO ib2 = 1, nbdysegn
1343               IF (( jpindt(ib2)<=jpieob(ib1)+1).AND. &
1344                &  ( jpinft(ib2)>=jpieob(ib1)+1).AND. &
1345                &  ( jpjedt(ib1)<=jpjnob(ib2)+1).AND. &
1346                &  ( jpjeft(ib1)>=jpjnob(ib2)+1)) THEN
1347                  IF ((jpjeft(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpieob(ib1)+1)) THEN
1348                     ! We have a possible North-East corner
1349!                     WRITE(numout,*) ' Found a North-East corner at (i,j): ', jpinft(ib2), jpjeft(ib1)
1350!                     WRITE(numout,*) ' between segments: ', npckge(ib1), npckgn(ib2)
1351                     icorne(ib1,2) = npckgn(ib2)
1352                     icornn(ib2,2) = npckge(ib1)
1353                  ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
1354                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
1355                     &                                     jpindt(ib2), jpjedt(ib1)
1356                     WRITE(ctmp2,*) ' Not allowed yet'
1357                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
1358                     &                               ' and North segment: ',npckgn(ib2)
1359                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
1360                  ELSE
1361                     WRITE(ctmp1,*) ' Check North and East Open boundary indices'
1362                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
1363                     &                               ' and North segment: ',npckgn(ib2)
1364                     CALL ctl_stop( ctmp1, ctmp2 )
1365                  END IF
1366               END IF
1367            END DO
1368         END DO
1369      END IF
1370      !
1371      ! 3. Check if segment extremities are on land
1372      !--------------------------------------------
1373      !
1374      ! West segments
1375      DO ib = 1, nbdysegw
1376         ! get mask at boundary extremities:
1377         ztestmask(1:2)=0.
1378         DO ji = 1, jpi
1379            DO jj = 1, jpj             
1380              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
1381               &  ((jj + njmpp - 1) == jpjwdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1382              IF (((ji + nimpp - 1) == jpiwob(ib)).AND. & 
1383               &  ((jj + njmpp - 1) == jpjwft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1384            END DO
1385         END DO
1386         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1387
1388         IF (ztestmask(1)==1) THEN
1389            IF (icornw(ib,1)==0) THEN
1390               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1391               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1392            ELSE
1393               ! This is a corner
1394               IF(lwp) WRITE(numout,*) 'Found a South-West corner at (i,j): ', jpiwob(ib), jpjwdt(ib)
1395               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,1))
1396               itest=itest+1
1397            ENDIF
1398         ENDIF
1399         IF (ztestmask(2)==1) THEN
1400            IF (icornw(ib,2)==0) THEN
1401               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
1402               CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' )
1403            ELSE
1404               ! This is a corner
1405               IF(lwp) WRITE(numout,*) 'Found a North-West corner at (i,j): ', jpiwob(ib), jpjwft(ib)
1406               CALL bdy_ctl_corn(npckgw(ib), icornw(ib,2))
1407               itest=itest+1
1408            ENDIF
1409         ENDIF
1410      END DO
1411      !
1412      ! East segments
1413      DO ib = 1, nbdysege
1414         ! get mask at boundary extremities:
1415         ztestmask(1:2)=0.
1416         DO ji = 1, jpi
1417            DO jj = 1, jpj             
1418              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
1419               &  ((jj + njmpp - 1) == jpjedt(ib))) ztestmask(1)=tmask(ji,jj,1)
1420              IF (((ji + nimpp - 1) == jpieob(ib)+1).AND. & 
1421               &  ((jj + njmpp - 1) == jpjeft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1422            END DO
1423         END DO
1424         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1425
1426         IF (ztestmask(1)==1) THEN
1427            IF (icorne(ib,1)==0) THEN
1428               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1429               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1430            ELSE
1431               ! This is a corner
1432               IF(lwp) WRITE(numout,*) 'Found a South-East corner at (i,j): ', jpieob(ib)+1, jpjedt(ib)
1433               CALL bdy_ctl_corn(npckge(ib), icorne(ib,1))
1434               itest=itest+1
1435            ENDIF
1436         ENDIF
1437         IF (ztestmask(2)==1) THEN
1438            IF (icorne(ib,2)==0) THEN
1439               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
1440               CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1441            ELSE
1442               ! This is a corner
1443               IF(lwp) WRITE(numout,*) 'Found a North-East corner at (i,j): ', jpieob(ib)+1, jpjeft(ib)
1444               CALL bdy_ctl_corn(npckge(ib), icorne(ib,2))
1445               itest=itest+1
1446            ENDIF
1447         ENDIF
1448      END DO
1449      !
1450      ! South segments
1451      DO ib = 1, nbdysegs
1452         ! get mask at boundary extremities:
1453         ztestmask(1:2)=0.
1454         DO ji = 1, jpi
1455            DO jj = 1, jpj             
1456              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
1457               &  ((ji + nimpp - 1) == jpisdt(ib))) ztestmask(1)=tmask(ji,jj,1)
1458              IF (((jj + njmpp - 1) == jpjsob(ib)).AND. & 
1459               &  ((ji + nimpp - 1) == jpisft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1460            END DO
1461         END DO
1462         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1463
1464         IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
1465            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1466            CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
1467         ENDIF
1468         IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
1469            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
1470            CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
1471         ENDIF
1472      END DO
1473      !
1474      ! North segments
1475      DO ib = 1, nbdysegn
1476         ! get mask at boundary extremities:
1477         ztestmask(1:2)=0.
1478         DO ji = 1, jpi
1479            DO jj = 1, jpj             
1480              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
1481               &  ((ji + nimpp - 1) == jpindt(ib))) ztestmask(1)=tmask(ji,jj,1)
1482              IF (((jj + njmpp - 1) == jpjnob(ib)+1).AND. & 
1483               &  ((ji + nimpp - 1) == jpinft(ib))) ztestmask(2)=tmask(ji,jj,1) 
1484            END DO
1485         END DO
1486         CALL mpp_sum( 'bdyini', ztestmask, 2 )   ! sum over the global domain
1487
1488         IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
1489            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1490            CALL ctl_stop( ctmp1, ' does not start on land' )
1491         ENDIF
1492         IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
1493            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
1494            CALL ctl_stop( ctmp1, ' does not end on land' )
1495         ENDIF
1496      END DO
1497      !
1498      IF ((itest==0).AND.(lwp)) WRITE(numout,*) 'NO open boundary corner found'
1499      !
1500      ! Other tests TBD:
1501      ! segments completly on land
1502      ! optimized open boundary array length according to landmask
1503      ! Nudging layers that overlap with interior domain
1504      !
1505   END SUBROUTINE bdy_ctl_seg
1506
1507   
1508   SUBROUTINE bdy_coords_seg( nbidta, nbjdta, nbrdta ) 
1509      !!----------------------------------------------------------------------
1510      !!                 ***  ROUTINE bdy_coords_seg  ***
1511      !!
1512      !! ** Purpose :  build nbidta, nbidta, nbrdta for bdy built with segments
1513      !!
1514      !! ** Method  : 
1515      !!
1516      !!----------------------------------------------------------------------
1517      INTEGER, DIMENSION(:,:,:), intent(  out)  ::   nbidta, nbjdta, nbrdta   ! Index arrays: i and j indices of bdy dta
1518      !!
1519      INTEGER  ::   ii, ij, ir, iseg
1520      INTEGER  ::   igrd         ! grid type (t=1, u=2, v=3)
1521      INTEGER  ::   icount       !
1522      INTEGER  ::   ib_bdy       ! bdy number
1523      !!----------------------------------------------------------------------
1524
1525      ! East
1526      !-----
1527      DO iseg = 1, nbdysege
1528         ib_bdy = npckge(iseg)
1529         !
1530         ! ------------ T points -------------
1531         igrd=1
1532         icount=0
1533         DO ir = 1, nn_rimwidth(ib_bdy)
1534            DO ij = jpjedt(iseg), jpjeft(iseg)
1535               icount = icount + 1
1536               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
1537               nbjdta(icount, igrd, ib_bdy) = ij
1538               nbrdta(icount, igrd, ib_bdy) = ir
1539            ENDDO
1540         ENDDO
1541         !
1542         ! ------------ U points -------------
1543         igrd=2
1544         icount=0
1545         DO ir = 1, nn_rimwidth(ib_bdy)
1546            DO ij = jpjedt(iseg), jpjeft(iseg)
1547               icount = icount + 1
1548               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 1 - ir
1549               nbjdta(icount, igrd, ib_bdy) = ij
1550               nbrdta(icount, igrd, ib_bdy) = ir
1551            ENDDO
1552         ENDDO
1553         !
1554         ! ------------ V points -------------
1555         igrd=3
1556         icount=0
1557         DO ir = 1, nn_rimwidth(ib_bdy)
1558            !            DO ij = jpjedt(iseg), jpjeft(iseg) - 1
1559            DO ij = jpjedt(iseg), jpjeft(iseg)
1560               icount = icount + 1
1561               nbidta(icount, igrd, ib_bdy) = jpieob(iseg) + 2 - ir
1562               nbjdta(icount, igrd, ib_bdy) = ij
1563               nbrdta(icount, igrd, ib_bdy) = ir
1564            ENDDO
1565            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1566            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1567         ENDDO
1568      ENDDO
1569      !
1570      ! West
1571      !-----
1572      DO iseg = 1, nbdysegw
1573         ib_bdy = npckgw(iseg)
1574         !
1575         ! ------------ T points -------------
1576         igrd=1
1577         icount=0
1578         DO ir = 1, nn_rimwidth(ib_bdy)
1579            DO ij = jpjwdt(iseg), jpjwft(iseg)
1580               icount = icount + 1
1581               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
1582               nbjdta(icount, igrd, ib_bdy) = ij
1583               nbrdta(icount, igrd, ib_bdy) = ir
1584            ENDDO
1585         ENDDO
1586         !
1587         ! ------------ U points -------------
1588         igrd=2
1589         icount=0
1590         DO ir = 1, nn_rimwidth(ib_bdy)
1591            DO ij = jpjwdt(iseg), jpjwft(iseg)
1592               icount = icount + 1
1593               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
1594               nbjdta(icount, igrd, ib_bdy) = ij
1595               nbrdta(icount, igrd, ib_bdy) = ir
1596            ENDDO
1597         ENDDO
1598         !
1599         ! ------------ V points -------------
1600         igrd=3
1601         icount=0
1602         DO ir = 1, nn_rimwidth(ib_bdy)
1603            !            DO ij = jpjwdt(iseg), jpjwft(iseg) - 1
1604            DO ij = jpjwdt(iseg), jpjwft(iseg)
1605               icount = icount + 1
1606               nbidta(icount, igrd, ib_bdy) = jpiwob(iseg) + ir - 1
1607               nbjdta(icount, igrd, ib_bdy) = ij
1608               nbrdta(icount, igrd, ib_bdy) = ir
1609            ENDDO
1610            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1611            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1612         ENDDO
1613      ENDDO
1614      !
1615      ! North
1616      !-----
1617      DO iseg = 1, nbdysegn
1618         ib_bdy = npckgn(iseg)
1619         !
1620         ! ------------ T points -------------
1621         igrd=1
1622         icount=0
1623         DO ir = 1, nn_rimwidth(ib_bdy)
1624            DO ii = jpindt(iseg), jpinft(iseg)
1625               icount = icount + 1
1626               nbidta(icount, igrd, ib_bdy) = ii
1627               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir 
1628               nbrdta(icount, igrd, ib_bdy) = ir
1629            ENDDO
1630         ENDDO
1631         !
1632         ! ------------ U points -------------
1633         igrd=2
1634         icount=0
1635         DO ir = 1, nn_rimwidth(ib_bdy)
1636            !            DO ii = jpindt(iseg), jpinft(iseg) - 1
1637            DO ii = jpindt(iseg), jpinft(iseg)
1638               icount = icount + 1
1639               nbidta(icount, igrd, ib_bdy) = ii
1640               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 2 - ir
1641               nbrdta(icount, igrd, ib_bdy) = ir
1642            ENDDO
1643            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1644            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1645         ENDDO
1646         !
1647         ! ------------ V points -------------
1648         igrd=3
1649         icount=0
1650         DO ir = 1, nn_rimwidth(ib_bdy)
1651            DO ii = jpindt(iseg), jpinft(iseg)
1652               icount = icount + 1
1653               nbidta(icount, igrd, ib_bdy) = ii
1654               nbjdta(icount, igrd, ib_bdy) = jpjnob(iseg) + 1 - ir
1655               nbrdta(icount, igrd, ib_bdy) = ir
1656            ENDDO
1657         ENDDO
1658      ENDDO
1659      !
1660      ! South
1661      !-----
1662      DO iseg = 1, nbdysegs
1663         ib_bdy = npckgs(iseg)
1664         !
1665         ! ------------ T points -------------
1666         igrd=1
1667         icount=0
1668         DO ir = 1, nn_rimwidth(ib_bdy)
1669            DO ii = jpisdt(iseg), jpisft(iseg)
1670               icount = icount + 1
1671               nbidta(icount, igrd, ib_bdy) = ii
1672               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
1673               nbrdta(icount, igrd, ib_bdy) = ir
1674            ENDDO
1675         ENDDO
1676         !
1677         ! ------------ U points -------------
1678         igrd=2
1679         icount=0
1680         DO ir = 1, nn_rimwidth(ib_bdy)
1681            !            DO ii = jpisdt(iseg), jpisft(iseg) - 1
1682            DO ii = jpisdt(iseg), jpisft(iseg)
1683               icount = icount + 1
1684               nbidta(icount, igrd, ib_bdy) = ii
1685               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
1686               nbrdta(icount, igrd, ib_bdy) = ir
1687            ENDDO
1688            nbidta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1689            nbjdta(icount, igrd, ib_bdy) = -ib_bdy ! Discount this point
1690         ENDDO
1691         !
1692         ! ------------ V points -------------
1693         igrd=3
1694         icount=0
1695         DO ir = 1, nn_rimwidth(ib_bdy)
1696            DO ii = jpisdt(iseg), jpisft(iseg)
1697               icount = icount + 1
1698               nbidta(icount, igrd, ib_bdy) = ii
1699               nbjdta(icount, igrd, ib_bdy) = jpjsob(iseg) + ir - 1
1700               nbrdta(icount, igrd, ib_bdy) = ir
1701            ENDDO
1702         ENDDO
1703      ENDDO
1704
1705     
1706   END SUBROUTINE bdy_coords_seg
1707   
1708   
1709   SUBROUTINE bdy_ctl_corn( ib1, ib2 )
1710      !!----------------------------------------------------------------------
1711      !!                 ***  ROUTINE bdy_ctl_corn  ***
1712      !!
1713      !! ** Purpose :   Check numerical schemes consistency between
1714      !!                segments having a common corner
1715      !!
1716      !! ** Method  :   
1717      !!----------------------------------------------------------------------
1718      INTEGER, INTENT(in)  ::   ib1, ib2
1719      INTEGER :: itest
1720      !!----------------------------------------------------------------------
1721      itest = 0
1722
1723      IF( cn_dyn2d(ib1) /= cn_dyn2d(ib2) )   itest = itest + 1
1724      IF( cn_dyn3d(ib1) /= cn_dyn3d(ib2) )   itest = itest + 1
1725      IF( cn_tra  (ib1) /= cn_tra  (ib2) )   itest = itest + 1
1726      !
1727      IF( nn_dyn2d_dta(ib1) /= nn_dyn2d_dta(ib2) )   itest = itest + 1
1728      IF( nn_dyn3d_dta(ib1) /= nn_dyn3d_dta(ib2) )   itest = itest + 1
1729      IF( nn_tra_dta  (ib1) /= nn_tra_dta  (ib2) )   itest = itest + 1
1730      !
1731      IF( nn_rimwidth(ib1) /= nn_rimwidth(ib2) )   itest = itest + 1   
1732      !
1733      IF( itest>0 ) THEN
1734         WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2
1735         CALL ctl_stop( ctmp1, ' have different open bdy schemes' )
1736      ENDIF
1737      !
1738   END SUBROUTINE bdy_ctl_corn
1739
1740
1741   SUBROUTINE bdy_meshwri()
1742      !!----------------------------------------------------------------------
1743      !!                 ***  ROUTINE bdy_meshwri  ***
1744      !!         
1745      !! ** Purpose :   write netcdf file with nbr, flagu, flagv, ntreat for T, U
1746      !!                and V points in 2D arrays for easier visualisation/control
1747      !!
1748      !! ** Method  :   use iom_rstput as in domwri.F
1749      !!----------------------------------------------------------------------     
1750      INTEGER  ::   ib_bdy, ii, ij, igrd, ib     ! dummy loop indices
1751      INTEGER  ::   inum                                   !   -       -
1752      REAL(wp), POINTER, DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields
1753      REAL(wp)         , DIMENSION(jpi,jpj) ::   ztmp
1754      CHARACTER(LEN=1) , DIMENSION(jpbgrd)  ::   cgrid
1755      !!----------------------------------------------------------------------     
1756      cgrid = (/'t','u','v'/)
1757      CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. )
1758      DO igrd = 1, jpbgrd
1759         SELECT CASE( igrd )
1760         CASE( 1 )   ;   zmask => tmask(:,:,1)
1761         CASE( 2 )   ;   zmask => umask(:,:,1)
1762         CASE( 3 )   ;   zmask => vmask(:,:,1)
1763         END SELECT
1764         ztmp(:,:) = zmask(:,:)
1765         DO ib_bdy = 1, nb_bdy
1766            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)      ! nbr deined for all rims
1767               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1768               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1769               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10.
1770               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1771            END DO
1772         END DO
1773         CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1774         ztmp(:,:) = zmask(:,:)
1775         DO ib_bdy = 1, nb_bdy
1776            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagu defined only for rims 0 and 1
1777               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1778               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1779               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10.
1780               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1781            END DO
1782         END DO
1783         CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1784         ztmp(:,:) = zmask(:,:)
1785         DO ib_bdy = 1, nb_bdy
1786            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagv defined only for rims 0 and 1
1787               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1788               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1789               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10.
1790               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1791            END DO
1792         END DO
1793         CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1794         ztmp(:,:) = zmask(:,:)
1795         DO ib_bdy = 1, nb_bdy
1796            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! ntreat defined only for rims 0 and 1
1797               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
1798               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
1799               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10.
1800               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
1801            END DO
1802         END DO
1803         CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
1804      END DO
1805      CALL iom_close( inum )
1806
1807   END SUBROUTINE bdy_meshwri
1808   
1809   !!=================================================================================
1810END MODULE bdyini
Note: See TracBrowser for help on using the repository browser.