Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90	(revision 11182)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90	(revision 11183)
@@ -34,4 +34,5 @@
 
    PUBLIC   bdy_init   ! routine called in nemo_init
+   PUBLIC   find_neib   ! routine called in bdy_nmn
 
    INTEGER, PARAMETER ::   jp_nseg = 100   ! 
@@ -131,9 +132,10 @@
       INTEGER  ::   jpbdtau, jpbdtas                       !   -       -
       INTEGER  ::   ib_bdy1, ib_bdy2, ib1, ib2             !   -       -
-      INTEGER  ::   i_offset, j_offset, inbdy, itreat      !   -       -
-      INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3, iibe, ijbe   !   -       -
-      INTEGER  ::   flagu, flagv                               ! short cuts
+      INTEGER  ::   i_offset, j_offset, inn                !   -       -
+      INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3           !   -       -
+      INTEGER  ::   iibe, ijbe, iibi, ijbi                 !   -       -
+      INTEGER  ::   flagu, flagv                           ! short cuts
       INTEGER , POINTER  ::  nbi, nbj, nbr                 ! short cuts
-      REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask    ! pointer to 2D mask fields
+      REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask   ! pointer to 2D mask fields
       REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
       INTEGER, DIMENSION (2)                  ::   kdimsz
@@ -147,5 +149,5 @@
       REAL(wp), TARGET, DIMENSION(jpi,jpj) ::   zfmask  ! temporary fmask array excluding coastal boundary condition (shlat)
       REAL(wp)        , DIMENSION(jpi,jpj) ::   ztmp
-      LOGICAL  ::   llnobdy, llsobdy, lleabdy, llwebdy     ! local logicals
+      LOGICAL  ::   llnon, llson, llean, llwen          ! local logicals
       !!
       CHARACTER(LEN=1)                     ::   ctypebdy   !     -        - 
@@ -1072,8 +1074,4 @@
       END DO
       
-      ! detect corner interior and its orientation index 1 to 4  depending on the orientation
-      ! detect corner exterior and its orientation index 5 to 8  depending on the orientation
-      ! detect geometries with 3 neighbours        index 9 to 12 depending on the orientation
-      ! else                                       index 0
       DO ib_bdy = 1, nb_bdy
          DO igrd = 1, jpbgrd
@@ -1087,43 +1085,56 @@
                ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE
-               llnobdy = pmask(ii  ,ij+1) == 1.  
-               llsobdy = pmask(ii  ,ij-1) == 1. 
-               lleabdy = pmask(ii+1,ij  ) == 1. 
-               llwebdy = pmask(ii-1,ij  ) == 1. 
-               inbdy  = COUNT( (/ llnobdy, llsobdy, lleabdy, llwebdy /) )
-               IF( inbdy == 0 )   THEN   ! no neighbours -> interior of a corner
-                  !               !              !     _____     !     _____     
-                  !  1 |   o      !  2  o   |    !  3 | x        !  4     x |    
-                  !    |_x_ _     !    _ _x_|    !    |   o      !      o   |
-                  IF( pmask(ii+1,ij+1) == 1. )   ztmp(ii,ij) = 1.
-                  IF( pmask(ii-1,ij+1) == 1. )   ztmp(ii,ij) = 2.
-                  IF( pmask(ii+1,ij-1) == 1. )   ztmp(ii,ij) = 3.
-                  IF( pmask(ii-1,ij-1) == 1. )   ztmp(ii,ij) = 4.
+               IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE 
+               llean = pmask(ii+1,ij  ) == 1.
+               llwen = pmask(ii-1,ij  ) == 1.
+               llnon = pmask(ii  ,ij+1) == 1.
+               llson = pmask(ii  ,ij-1) == 1.
+               inn  = COUNT( (/ llnon, llson, llean, llwen /) )
+               IF( inn == 0 ) THEN   ! no neighbours -> interior of a corner  or  cluster of rim points
+                  !               !              !     _____     !     _____    !    __     __
+                  !  1 |   o      !  2  o   |    !  3 | x        !  4     x |   !      |   |   -> error
+                  !    |_x_ _     !    _ _x_|    !    |   o      !      o   |   !      |x_x|
+                  IF(     pmask(ii+1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 1.
+                  ELSEIF( pmask(ii-1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 2.
+                  ELSEIF( pmask(ii+1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 3.
+                  ELSEIF( pmask(ii-1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 4.
+                  ELSE
+                     WRITE(ctmp1,*) ' E R R O R : Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
+                          ' on boundary set ', ib_bdy, ' has no free ocean neighbour'
+                     WRITE(ctmp2,*) ' There seems to be a cluster of rim points.'
+                     WRITE(ctmp3,*) ' ========== '
+                     CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' )
+                  END IF
                END IF
-               IF( inbdy == 1 )   THEN   ! middle of linear bdy
-                  ztmp(ii,ij) = 0.   ! regular treatment with flags
+               IF( inn == 1 ) THEN   ! middle of linear bdy  or incomplete corner  ! ___ o
+                  !    |         !         |   !      o     !    ______            !    |x___
+                  ! 5  | x o     ! 6   o x |   ! 7  __x__   ! 8    x
+                  !    |         !         |   !            !      o
+                  IF( llean )   ztmp(ii,ij) = 5.
+                  IF( llwen )   ztmp(ii,ij) = 6.
+                  IF( llnon )   ztmp(ii,ij) = 7.
+                  IF( llson )   ztmp(ii,ij) = 8.
                END IF
-               IF( inbdy == 2 )   THEN   ! exterior of a corner
+               IF( inn == 2 ) THEN   ! exterior of a corner
                   !        o      !        o      !    _____|       !       |_____  
-                  !  5 ____x o    !  6   o x___   ! 7      x o      !  8   o x      
+                  !  9 ____x o    ! 10   o x___   ! 11     x o      ! 12   o x      
                   !         |     !       |       !        o        !        o 
-                  IF( llnobdy .AND. lleabdy )   ztmp(ii,ij) = 5.
-                  IF( llnobdy .AND. llwebdy )   ztmp(ii,ij) = 6.
-                  IF( llsobdy .AND. lleabdy )   ztmp(ii,ij) = 7.
-                  IF( llsobdy .AND. llwebdy )   ztmp(ii,ij) = 8.
+                  IF( llnon .AND. llean )   ztmp(ii,ij) =  9.
+                  IF( llnon .AND. llwen )   ztmp(ii,ij) = 10.
+                  IF( llson .AND. llean )   ztmp(ii,ij) = 11.
+                  IF( llson .AND. llwen )   ztmp(ii,ij) = 12.
                END IF
-               IF( inbdy == 3 )   THEN   ! 3 neighbours __   __
+               IF( inn == 3 ) THEN   ! 3 neighbours     __   __
                   !    |_  o      !        o  _|  !       |_|     !       o         
-                  !  9  _| x o    ! 10   o x |_   ! 11   o x o    ! 12  o x o       
+                  ! 13  _| x o    ! 14   o x |_   ! 15   o x o    ! 16  o x o       
                   !    |   o      !        o   |  !        o      !    __|¨|__    
-                  IF( llnobdy .AND. lleabdy .AND. llsobdy )   ztmp(ii,ij) = 9.
-                  IF( llnobdy .AND. llwebdy .AND. llsobdy )   ztmp(ii,ij) = 10.
-                  IF( llwebdy .AND. llsobdy .AND. lleabdy )   ztmp(ii,ij) = 11.
-                  IF( llwebdy .AND. llnobdy .AND. lleabdy )   ztmp(ii,ij) = 12.
+                  IF( llnon .AND. llean .AND. llson )   ztmp(ii,ij) = 13.
+                  IF( llnon .AND. llwen .AND. llson )   ztmp(ii,ij) = 14.
+                  IF( llwen .AND. llson .AND. llean )   ztmp(ii,ij) = 15.
+                  IF( llwen .AND. llnon .AND. llean )   ztmp(ii,ij) = 16.
                END IF
-               IF( inbdy == 4 )   THEN
-                  WRITE(ctmp1,*) ' E R R O R : Problem with  ',cgrid(igrd) ,' grid points,',   &
-                       '  some points on boundary set ', ib_bdy, ' have 4 neighbours'
+               IF( inn == 4 ) THEN
+                  WRITE(ctmp1,*)  ' E R R O R : Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
+                       ' on boundary set ', ib_bdy, ' have 4 neighbours'
                   WRITE(ctmp2,*) ' ========== '
                   CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
@@ -1131,5 +1142,5 @@
             END DO
             CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )
-            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)  
+            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
                ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
@@ -1156,21 +1167,8 @@
                flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd))
                iibe = ii - flagu   ! neighbouring point towards the exterior of the computational domain
-               ijbe = ij - flagv   ! neighbouring point towards the exterior of the computational domain
-               SELECT CASE( idx_bdy(ib_bdy)%ntreat(ib,igrd) )   ! points that will be used by bdy routines, -1 will be discarded
-               CASE( 0 )    ;   ii1 = ii + flagu   ;   ii2 = -1   ;   ii3 = -1
-                                ij1 = ij + flagv   ;   ij2 = -1   ;   ij3 = -1
-               CASE( 1  )   ;   ii1 = ii+1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
-               CASE( 2  )   ;   ii1 = ii-1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
-               CASE( 3  )   ;   ii1 = ii+1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
-               CASE( 4  )   ;   ii1 = ii-1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
-               CASE( 5  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
-               CASE( 6  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
-               CASE( 7  )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
-               CASE( 8  )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
-               CASE( 9  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1   
-               CASE( 10 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1 
-               CASE( 11 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij-1   ;   ii3 = ii+1   ;   ij3 = ij   
-               CASE( 12 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij
-               END SELECT
+               ijbe = ij - flagv
+               iibi = ii + flagu   ! neighbouring point towards the interior of the computational domain
+               ijbi = ij + flagv
+               CALL find_neib( ii, ij, idx_bdy(ib_bdy)%ntreat(ib,igrd), ii1, ij1, ii2, ij2, ii3, ij3 )
                !
                ! search neighbour in the  west/east  direction
@@ -1179,6 +1177,6 @@
                ! (1)  o|x         OR    (2)   x|o
                !       |___                 ___| 
-               IF( ii1 == 0     .OR. ii2 == 0     .OR. ii3 == 0     )   lrecv_bdyint(ib_bdy,igrd,1) = .true.
-               IF( ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 )   lrecv_bdyint(ib_bdy,igrd,2) = .true.  
+               IF( iibi == 0     .OR. ii1 == 0     .OR. ii2 == 0     .OR. ii3 == 0     )   lrecv_bdyint(ib_bdy,igrd,1) = .true.
+               IF( iibi == jpi+1 .OR. ii1 == jpi+1 .OR. ii2 == jpi+1 .OR. ii3 == jpi+1 )   lrecv_bdyint(ib_bdy,igrd,2) = .true.  
                IF( iibe == 0     )   lrecv_bdyext(ib_bdy,igrd,1) = .true.
                IF( iibe == jpi+1 )   lrecv_bdyext(ib_bdy,igrd,2) = .true.  
@@ -1187,8 +1185,8 @@
                ! :     |  x:o    |    neighbour limited by ... would need o    |    o:x  |     :
                ! :.....|_._:_____|   (1) W neighbour         E neighbour (2)   |_____:_._|.....:
-               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) &
-                    &          .AND. ( ii1 == 3     .OR. ii2 == 3     .OR. ii3 == 3     ) )   lsend_bdyint(ib_bdy,igrd,1) = .true.
-               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) &
-                    &          .AND. ( ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2 ) )   lsend_bdyint(ib_bdy,igrd,2) = .true.
+               IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. &
+                    & ( iibi == 3     .OR. ii1 == 3     .OR. ii2 == 3     .OR. ii3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,1) = .true.
+               IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. &
+                    & ( iibi == jpi-2 .OR. ii1 == jpi-2 .OR. ii2 == jpi-2 .OR. ii3 == jpi-2) )   lsend_bdyint(ib_bdy,igrd,2) = .true.
                IF( ii == 2     .AND. ( nbondi ==  1 .OR. nbondi == 0 ) .AND. iibe == 3     )   lsend_bdyext(ib_bdy,igrd,1) = .true.
                IF( ii == jpi-1 .AND. ( nbondi == -1 .OR. nbondi == 0 ) .AND. iibe == jpi-2 )   lsend_bdyext(ib_bdy,igrd,2) = .true.
@@ -1199,6 +1197,6 @@
                !  |   |___x___|   OR    |  |   x   |
                !  v       o           (4)  |       |
-               IF( ij1 == 0     .OR. ij2 == 0     .OR. ij3 == 0     )   lrecv_bdyint(ib_bdy,igrd,3) = .true.
-               IF( ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 )   lrecv_bdyint(ib_bdy,igrd,4) = .true.
+               IF( ijbi == 0     .OR. ij1 == 0     .OR. ij2 == 0     .OR. ij3 == 0     )   lrecv_bdyint(ib_bdy,igrd,3) = .true.
+               IF( ijbi == jpj+1 .OR. ij1 == jpj+1 .OR. ij2 == jpj+1 .OR. ij3 == jpj+1 )   lrecv_bdyint(ib_bdy,igrd,4) = .true.
                IF( ijbe == 0     )   lrecv_bdyext(ib_bdy,igrd,3) = .true.
                IF( ijbe == jpj+1 )   lrecv_bdyext(ib_bdy,igrd,4) = .true.
@@ -1207,8 +1205,8 @@
                !   |  |¨¨¨¨x¨¨¨¨|   neighbour limited by ... would need o     |  |....x....|
                !      :_________:  (3) S neighbour          N neighbour (4)   v  |    o    |   
-               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) &
-                    &          .AND. ( ij1 == 3     .OR. ij2 == 3     .OR. ij3 == 3     ) )   lsend_bdyint(ib_bdy,igrd,3) = .true.
-               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) &
-                    &          .AND. ( ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2 ) )   lsend_bdyint(ib_bdy,igrd,4) = .true.
+               IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. &
+                    & ( ijbi == 3     .OR. ij1 == 3     .OR. ij2 == 3     .OR. ij3 == 3    ) )   lsend_bdyint(ib_bdy,igrd,3) = .true.
+               IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. &
+                    & ( ijbi == jpj-2 .OR. ij1 == jpj-2 .OR. ij2 == jpj-2 .OR. ij3 == jpj-2) )   lsend_bdyint(ib_bdy,igrd,4) = .true.
                IF( ij == 2     .AND. ( nbondj ==  1 .OR. nbondj == 0 ) .AND. ijbe == 3     )   lsend_bdyext(ib_bdy,igrd,3) = .true.
                IF( ij == jpj-1 .AND. ( nbondj == -1 .OR. nbondj == 0 ) .AND. ijbe == jpj-2 )   lsend_bdyext(ib_bdy,igrd,4) = .true.
@@ -1234,5 +1232,5 @@
          END IF
       END DO
-      
+
       !
       ! Tidy up
@@ -1241,4 +1239,48 @@
       !
    END SUBROUTINE bdy_segs
+
+
+   SUBROUTINE find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 )
+     !!----------------------------------------------------------------------
+     !!                 ***  ROUTINE find_neib  ***
+     !!
+     !! ** Purpose :   get ii1, ij1, ii2, ij2, ii3, ij3, the indices of
+     !!               the free ocean neighbours of (ii,ij) for bdy treatment
+     !!
+     !!----------------------------------------------------------------------
+     INTEGER, INTENT(in   )      ::   ii, ij, itreat
+     INTEGER, INTENT(  out)      ::   ii1, ij1, ii2, ij2, ii3, ij3
+     !!----------------------------------------------------------------------
+     SELECT CASE( itreat )   ! points that will be used by bdy routines, -1 will be discarded
+        !               !               !     _____     !     _____     
+        !  1 |   o      !  2  o   |     !  3 | x        !  4     x |    
+        !    |_x_ _     !    _ _x_|     !    |   o      !      o   |
+     CASE( 1 )    ;   ii1 = ii+1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
+     CASE( 2 )    ;   ii1 = ii-1   ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
+     CASE( 3 )    ;   ii1 = ii+1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
+     CASE( 4 )    ;   ii1 = ii-1   ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
+        !    |          !         |     !      o        !    ______                   ! or incomplete corner
+        ! 5  | x o      ! 6   o x |     ! 7  __x__      ! 8    x                      !  7  ____ o
+        !    |          !         |     !               !      o                      !         |x___
+     CASE( 5 )    ;   ii1 = ii+1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
+     CASE( 6 )    ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
+     CASE( 7 )    ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
+     CASE( 8 )    ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = -1     ;   ij2 = -1     ;   ii3 = -1     ;   ij3 = -1
+        !        o      !        o      !    _____|     !       |_____  
+        !  9 ____x o    ! 10   o x___   ! 11     x o    ! 12   o x      
+        !         |     !       |       !        o      !        o      
+     CASE( 9  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1 
+     CASE( 10 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
+     CASE( 11 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
+     CASE( 12 )   ;   ii1 = ii     ;   ij1 = ij-1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = -1     ;   ij3 = -1
+        !    |_  o      !        o  _|  !     ¨¨|_|¨¨   !       o         
+        ! 13  _| x o    !  14  o x |_   !  15  o x o    ! 16  o x o       
+        !    |   o      !        o   |  !        o      !    __|¨|__ 
+     CASE( 13 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1   
+     CASE( 14 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1 
+     CASE( 15 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij-1   ;   ii3 = ii+1   ;   ij3 = ij   
+     CASE( 16 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij
+     END SELECT
+   END SUBROUTINE find_neib
 
 
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdylib.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdylib.F90	(revision 11182)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdylib.F90	(revision 11183)
@@ -15,4 +15,5 @@
    USE bdy_oce        ! ocean open boundary conditions
    USE phycst         ! physical constants
+   USE bdyini
    !
    USE in_out_manager !
@@ -452,6 +453,5 @@
       INTEGER  ::   ii, ij   ! 2D addresses
       INTEGER  ::   ipkm1    ! size of phia third dimension minus 1
-      INTEGER  ::   flagu, flagv                        ! short cuts
-      INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3
+      INTEGER  ::   ii1, ii2, ii3, ij1, ij2, ij3, itreat
       !!----------------------------------------------------------------------
       !
@@ -468,56 +468,28 @@
          ii = idx%nbi(ib,igrd)
          ij = idx%nbj(ib,igrd)
-         !
-         SELECT CASE( idx%ntreat(ib,igrd) )   ! select free ocean neighbours
-            !     o 
-            !  ___x___         ! either flagu or flagv = 0
-            CASE( 0 )   ;   flagu = NINT( idx%flagu(ib,igrd) )   ;   ii1 = ii+flagu
-                            flagv = NINT( idx%flagv(ib,igrd) )   ;   ij1 = ij+flagv
-            !               !              !     _____     !     _____     
-            !  1 |   o      !  2  o   |    !  3 | x        !  4     x |    
-            !    |_x_ _     !    _ _x_|    !    |   o      !      o   |
-            CASE( 1 )   ;   ii1 = ii+1   ;   ij1 = ij+1
-            CASE( 2 )   ;   ii1 = ii-1   ;   ij1 = ij+1
-            CASE( 3 )   ;   ii1 = ii+1   ;   ij1 = ij-1
-            CASE( 4 )   ;   ii1 = ii-1   ;   ij1 = ij-1
-            !        o      !        o      !    _____|       !       |_____  
-            !  5 ____x o    !  6   o x___   ! 7      x o      !  8   o x      
-            !         |     !       |       !        o        !        o      
-            CASE( 5 )   ;   ii1 = ii   ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij   
-            CASE( 6 )   ;   ii1 = ii   ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij
-            CASE( 7 )   ;   ii1 = ii   ;   ij1 = ij-1   ;   ii2 = ii+1   ;   ij2 = ij  
-            CASE( 8 )   ;   ii1 = ii   ;   ij1 = ij-1   ;   ii2 = ii-1   ;   ij2 = ij  
-            !    |_  o      !        o  _|  !     ¨¨|_|¨¨   !       o         
-            !  9  _| x o    !  10  o x |_   !  11  o x o    ! 12  o x o       
-            !    |   o      !        o   |  !        o      !    __|¨|__ 
-            CASE( 9  )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii+1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1   
-            CASE( 10 )   ;   ii1 = ii     ;   ij1 = ij+1   ;   ii2 = ii-1   ;   ij2 = ij     ;   ii3 = ii     ;   ij3 = ij-1 
-            CASE( 11 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij-1   ;   ii3 = ii+1   ;   ij3 = ij   
-            CASE( 12 )   ;   ii1 = ii-1   ;   ij1 = ij     ;   ii2 = ii     ;   ij2 = ij+1   ;   ii3 = ii+1   ;   ij3 = ij
+         itreat = idx%ntreat(ib,igrd)
+         CALL find_neib( ii, ij, itreat, ii1, ij1, ii2, ij2, ii3, ij3 )
+         SELECT CASE( itreat )
+         CASE( 1:8 )
+            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE
+            DO ik = 1, ipkm1
+               IF( pmask(ii1,ij1,ik) /= 0. )   phia(ii,ij,ik) = phia(ii1,ij1,ik)  
+            END DO
+         CASE( 9:12 )
+            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE
+            IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE
+            DO ik = 1, ipkm1
+               zweight = pmask(ii1,ij1,ik) + pmask(ii2,ij2,ik)
+               IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) ) / zweight
+            END DO
+         CASE( 13:16 )
+            IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE
+            IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE
+            IF( ii3 < 1 .OR. ii3 > jpi .OR. ij3 < 1 .OR. ij3 > jpj )   CYCLE
+            DO ik = 1, ipkm1
+               zweight = pmask(ii1,ij1,ik) + pmask(ii2,ij2,ik) + pmask(ii3,ij3,ik)
+               IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) + phia(ii3,ij3,ik) ) / zweight
+            END DO
          END SELECT
-         !
-         SELECT CASE( idx%ntreat(ib,igrd) )
-            CASE( 0:4 )
-               IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE
-               DO ik = 1, ipkm1
-                  IF( pmask(ii1,ij1,ik) /= 0. )   phia(ii,ij,ik) = phia(ii1,ij1,ik)  
-               END DO
-            CASE( 5:8 )
-               IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE
-               IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE
-               DO ik = 1, ipkm1
-                  zweight = pmask(ii1,ij1,ik) + pmask(ii2,ij2,ik)
-                  IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia (ii2,ij2,ik) ) / zweight
-               END DO
-            CASE( 9:12 )
-               IF( ii1 < 1 .OR. ii1 > jpi .OR. ij1 < 1 .OR. ij1 > jpj )   CYCLE
-               IF( ii2 < 1 .OR. ii2 > jpi .OR. ij2 < 1 .OR. ij2 > jpj )   CYCLE
-               IF( ii3 < 1 .OR. ii3 > jpi .OR. ij3 < 1 .OR. ij3 > jpj )   CYCLE
-               DO ik = 1, ipkm1
-                  zweight = pmask(ii1,ij1,ik) + pmask(ii2,ij2,ik) + pmask(ii3,ij3,ik)
-                  IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia (ii2,ij2,ik) + phia (ii3,ij3,ik) ) / zweight
-               END DO
-         END SELECT
-         !
       END DO
       !
