New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
bdyini.F90 in NEMO/trunk/src/OCE/BDY – NEMO

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

Last change on this file since 12377 was 12377, checked in by acc, 4 years ago

The big one. Merging all 2019 developments from the option 1 branch back onto the trunk.

This changeset reproduces 2019/dev_r11943_MERGE_2019 on the trunk using a 2-URL merge
onto a working copy of the trunk. I.e.:

svn merge --ignore-ancestry \

svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/trunk \
svn+ssh://acc@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/NEMO/branches/2019/dev_r11943_MERGE_2019 ./

The --ignore-ancestry flag avoids problems that may otherwise arise from the fact that
the merge history been trunk and branch may have been applied in a different order but
care has been taken before this step to ensure that all applicable fixes and updates
are present in the merge branch.

The trunk state just before this step has been branched to releases/release-4.0-HEAD
and that branch has been immediately tagged as releases/release-4.0.2. Any fixes
or additions in response to tickets on 4.0, 4.0.1 or 4.0.2 should be done on
releases/release-4.0-HEAD. From now on future 'point' releases (e.g. 4.0.2) will
remain unchanged with periodic releases as needs demand. Note release-4.0-HEAD is a
transitional naming convention. Future full releases, say 4.2, will have a release-4.2
branch which fulfills this role and the first point release (e.g. 4.2.0) will be made
immediately following the release branch creation.

2020 developments can be started from any trunk revision later than this one.

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