source: NEMO/trunk/src/OCE/BDY/bdyini.F90 @ 11536

Last change on this file since 11536 was 11536, checked in by smasson, 13 months ago

trunk: merge dev_r10984_HPC-13 into the trunk

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