source: NEMO/branches/2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles/src/OCE/BDY/bdyini.F90 @ 11641

Last change on this file since 11641 was 11641, checked in by acc, 19 months ago

Branch 2019/dev_r11613_ENHANCE-04_namelists_as_internalfiles. Further substantive changes to the BDY routines that read multiple copies of some namelists from namelist_cfg. Special treatment is required to replicate the original behaviour with internal files. These changes emulate previous behaviour but removal of the reliance on multiple, same-named namelist blocks would be preferable.

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