Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn2d.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn2d.F90	(revision 11257)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdydyn2d.F90	(revision 11258)
@@ -199,5 +199,5 @@
       ELSE                ;   ibeg = idx%nblenrim0(igrd)+1   ;   iend = idx%nblenrim(igrd)
       END IF
-      sshdta(:,:) = 0.0
+      !
       DO jb = ibeg, iend
          ii = idx%nbi(jb,igrd)
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 11257)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdyini.F90	(revision 11258)
@@ -116,4 +116,5 @@
          ! Open boundaries definition (arrays and masks)
          CALL bdy_def
+         IF( ln_meshmask )   CALL bdy_meshwri()
          !
          ! Open boundaries initialisation of external data arrays
@@ -146,5 +147,5 @@
       INTEGER  ::   ib_bdy, ii, ij, igrd, ib, ir, iseg     ! dummy loop indices
       INTEGER  ::   icount, icountr, icountr0, ibr_max     ! local integers
-      INTEGER  ::   ilen1, ibm1                            !   -       -
+      INTEGER  ::   ilen1                                  !   -       -
       INTEGER  ::   iwe, ies, iso, ino, inum, id_dummy     !   -       -
       INTEGER  ::   jpbdta                                 !   -       -
@@ -153,15 +154,14 @@
       INTEGER  ::   iibe, ijbe, iibi, ijbi                 !   -       -
       INTEGER  ::   flagu, flagv                           ! short cuts
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   zz_read      ! work space for 2D global boundary data
-      REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask   ! pointer to 2D mask fields
-      INTEGER, DIMENSION (4)                  ::   kdimsz
-      INTEGER, DIMENSION(jpbgrd,jp_bdy)       ::   nblendta         ! Length of index arrays 
-      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbidta, nbjdta   ! Index arrays: i and j indices of bdy dta
-      INTEGER, ALLOCATABLE, DIMENSION(:,:,:)  ::   nbrdta           ! Discrete distance from rim points
-      CHARACTER(LEN=1),DIMENSION(jpbgrd)      ::   cgrid
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                  ! temporary fmask array excluding coastal boundary condition (shlat)
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztmask, zumask, zvmask  ! temporary u/v mask array
-      !!
-      INTEGER                              ::   nbdyind, nbdybeg, nbdyend
+      INTEGER  ::   nbdyind, nbdybeg, nbdyend
+      INTEGER              , DIMENSION(4)             ::   kdimsz
+      INTEGER              , DIMENSION(jpbgrd,jp_bdy) ::   nblendta          ! Length of index arrays 
+      INTEGER,  ALLOCATABLE, DIMENSION(:,:,:)         ::   nbidta, nbjdta    ! Index arrays: i and j indices of bdy dta
+      INTEGER,  ALLOCATABLE, DIMENSION(:,:,:)         ::   nbrdta            ! Discrete distance from rim points
+      CHARACTER(LEN=1)     , DIMENSION(jpbgrd)        ::   cgrid
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:)     ::   zz_read                 ! work space for 2D global boundary data
+      REAL(wp), POINTER    , DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields
+      REAL(wp)             , DIMENSION(jpi,jpj) ::   zfmask   ! temporary fmask array excluding coastal boundary condition (shlat)
+      REAL(wp)             , DIMENSION(jpi,jpj) ::   ztmask, zumask, zvmask  ! temporary u/v mask array
       !!----------------------------------------------------------------------
       !
@@ -483,5 +483,5 @@
          END DO
       END DO
-
+      !
       ! Find lenght of boundaries and rim on local mpi domain
       !------------------------------------------------------
@@ -502,7 +502,6 @@
             DO ib = 1, nblendta(igrd,ib_bdy)
                ! check that data is in correct order in file
-               ibm1 = MAX(1,ib-1)
-               IF(lwp) THEN         ! Since all procs read global data only need to do this check on one proc...
-                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ibm1,igrd,ib_bdy) ) THEN
+               IF( ib > 1 ) THEN
+                  IF( nbrdta(ib,igrd,ib_bdy) < nbrdta(ib-1,igrd,ib_bdy) ) THEN
                      CALL ctl_stop('bdy_segs : ERROR : boundary data in file must be defined ', &
                         &        ' in order of distance from edge nbr A utility for re-ordering ', &
@@ -636,8 +635,8 @@
       ! For the flagu/flagv calculation below we require a version of fmask without
       ! the land boundary condition (shlat) included:
-      DO ij = 2, jpjm1
-         DO ii = 2, jpim1
-            zfmask(ii,ij) = ztmask(ii,ij  ) * ztmask(ii+1,ij  )   &
-               &              * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
+      DO ij = 1, jpjm1
+         DO ii = 1, jpim1
+            zfmask(ii,ij) =  ztmask(ii,ij  ) * ztmask(ii+1,ij  )   &
+               &           * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
          END DO
       END DO
@@ -678,5 +677,6 @@
       ! handle rim1, do as if rim 0 was land
       ! ------------------------------------
-
+      
+      ! z[tuv]mask are now set to zero on rim 0 points:
       DO ib_bdy = 1, nb_bdy
          DO ib = 1, idx_bdy(ib_bdy)%nblenrim0(1)   ! extent of rim 0
@@ -692,8 +692,8 @@
 
       ! Recompute zfmask
-      DO ij = 2, jpjm1
-         DO ii = 2, jpim1
-            zfmask(ii,ij) = ztmask(ii,ij  ) * ztmask(ii+1,ij  )   &
-               &              * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
+      DO ij = 1, jpjm1
+         DO ii = 1, jpim1
+            zfmask(ii,ij) =  ztmask(ii,ij  ) * ztmask(ii+1,ij  )   &
+               &           * ztmask(ii,ij+1) * ztmask(ii+1,ij+1)
          END DO
       END DO
@@ -714,5 +714,4 @@
 
       CALL bdy_rim_treat( zumask, zvmask, zfmask, .false. )   ! compute flagu, flagv, ntreat on rim 1
-
       !
       ! Check which boundaries might need communication
@@ -727,9 +726,8 @@
          DO ib_bdy = 1, nb_bdy
             DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
+               IF( idx_bdy(ib_bdy)%ntreat(ib,igrd) == -1 ) CYCLE
                ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
                ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
-               IF( ib .LE. idx_bdy(ib_bdy)%nblenrim0(igrd) ) THEN   ;   ir = 0
-               ELSE                                                 ;   ir = 1
-               END IF
+               ir = idx_bdy(ib_bdy)%nbr(ib,igrd)
                flagu = NINT(idx_bdy(ib_bdy)%flagu(ib,igrd))
                flagv = NINT(idx_bdy(ib_bdy)%flagv(ib,igrd))
@@ -784,5 +782,5 @@
 
       DO ib_bdy = 1,nb_bdy
-         IF( cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. &
+         IF(  cn_dyn2d(ib_bdy) == 'orlanski' .OR. cn_dyn2d(ib_bdy) == 'orlanski_npo' .OR. &
             & cn_dyn3d(ib_bdy) == 'orlanski' .OR. cn_dyn3d(ib_bdy) == 'orlanski_npo' .OR. &
             & cn_tra(ib_bdy)   == 'orlanski' .OR. cn_tra(ib_bdy)   == 'orlanski_npo'      ) THEN
@@ -793,6 +791,5 @@
                   IF(  mig(ii) > 2 .AND. mig(ii) < jpiglo-2 .AND. mjg(ij) > 2 .AND. mjg(ij) < jpjglo-2  ) THEN
                      WRITE(ctmp1,*) ' Orlanski is not safe when the open boundaries are on the interior of the computational domain'
-                     WRITE(ctmp2,*) ' ========== '
-                     CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+                     CALL ctl_stop( ctmp1 )
                   END IF
                END DO
@@ -800,8 +797,5 @@
          END IF
       END DO
-
-      !
-      ! Tidy up
-      !--------
+      !
       DEALLOCATE( nbidta, nbjdta, nbrdta )
       !
@@ -809,5 +803,5 @@
 
 
-   SUBROUTINE bdy_rim_treat( zumask, zvmask, zfmask, lrim0 )
+   SUBROUTINE bdy_rim_treat( pumask, pvmask, pfmask, lrim0 )
       !!----------------------------------------------------------------------
       !!                 ***  ROUTINE bdy_rim_treat  ***
@@ -826,6 +820,6 @@
       !!                - and look at the ocean neighbours to compute ntreat
       !!----------------------------------------------------------------------
-      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: zfmask   ! temporary fmask excluding coastal boundary condition (shlat)
-      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: zumask, zvmask   ! temporary t/u/v mask array
+      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pfmask   ! temporary fmask excluding coastal boundary condition (shlat)
+      REAL(wp), TARGET, DIMENSION(jpi,jpj), INTENT (in   ) :: pumask, pvmask   ! temporary t/u/v mask array
       LOGICAL                             , INTENT (in   ) :: lrim0    ! .true. -> rim 0   .false. -> rim 1
       INTEGER  ::   ib_bdy, ii, ij, igrd, ib, icount       ! dummy loop indices
@@ -833,5 +827,5 @@
       INTEGER  ::   ibeg, iend                             ! local integer
       LOGICAL  ::   llnon, llson, llean, llwen             ! local logicals indicating the presence of a ocean neighbour
-      REAL(wp), POINTER, DIMENSION(:,:)       ::   pmask   ! pointer to 2D mask fields
+      REAL(wp), POINTER, DIMENSION(:,:)       ::   zmask   ! pointer to 2D mask fields
       REAL(wp) ::   zefl, zwfl, znfl, zsfl                 ! local scalars
       CHARACTER(LEN=1), DIMENSION(jpbgrd)     ::   cgrid
@@ -849,10 +843,10 @@
          DO igrd = 1, jpbgrd 
             SELECT CASE( igrd )
-               CASE( 1 )   ;   pmask => zumask     ;   i_offset = 0
-               CASE( 2 )   ;   pmask => bdytmask   ;   i_offset = 1
-               CASE( 3 )   ;   pmask => zfmask     ;   i_offset = 0
+               CASE( 1 )   ;   zmask => pumask     ;   i_offset = 0
+               CASE( 2 )   ;   zmask => bdytmask   ;   i_offset = 1
+               CASE( 3 )   ;   zmask => pfmask     ;   i_offset = 0
             END SELECT 
             icount = 0
-            ztmp(:,:) = 0._wp
+            ztmp(:,:) = -999._wp
             IF( lrim0 ) THEN   ! extent of rim 0
                ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
@@ -864,6 +858,6 @@
                ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
                IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE
-               zwfl = pmask(ii+i_offset-1,ij)
-               zefl = pmask(ii+i_offset  ,ij)
+               zwfl = zmask(ii+i_offset-1,ij)
+               zefl = zmask(ii+i_offset  ,ij)
                ! This error check only works if you are using the bdyXmask arrays
                IF( i_offset == 1 .and. zefl + zwfl == 2. ) THEN
@@ -875,10 +869,13 @@
             END DO
             IF( icount /= 0 ) THEN
-               WRITE(ctmp1,*) ' E R R O R : Some ',cgrid(igrd),' grid points,',   &
+               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   &
                   ' are not boundary points (flagu calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
-               WRITE(ctmp2,*) ' ========== '
-               CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+               CALL ctl_stop( ctmp1 )
             ENDIF 
-            CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
+            SELECT CASE( igrd )
+               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
+               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 
+               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 
+            END SELECT 
             DO ib = ibeg, iend
                ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -894,10 +891,10 @@
          DO igrd = 1, jpbgrd 
             SELECT CASE( igrd )
-               CASE( 1 )   ;   pmask => zvmask     ;   j_offset = 0
-               CASE( 2 )   ;   pmask => zfmask     ;   j_offset = 0
-               CASE( 3 )   ;   pmask => bdytmask   ;   j_offset = 1
+               CASE( 1 )   ;   zmask => pvmask     ;   j_offset = 0
+               CASE( 2 )   ;   zmask => pfmask     ;   j_offset = 0
+               CASE( 3 )   ;   zmask => bdytmask   ;   j_offset = 1
             END SELECT 
             icount = 0
-            ztmp(:,:) = 0._wp
+            ztmp(:,:) = -999._wp
             IF( lrim0 ) THEN   ! extent of rim 0
                ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
@@ -909,6 +906,6 @@
                ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
                IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )  CYCLE
-               zsfl = pmask(ii,ij+j_offset-1)
-               znfl = pmask(ii,ij+j_offset  )
+               zsfl = zmask(ii,ij+j_offset-1)
+               znfl = zmask(ii,ij+j_offset  )
                ! This error check only works if you are using the bdyXmask arrays
                IF( j_offset == 1 .and. znfl + zsfl == 2. ) THEN
@@ -920,10 +917,13 @@
             END DO
             IF( icount /= 0 ) THEN
-               WRITE(ctmp1,*) ' E R R O R : Some ',cgrid(igrd),' grid points,',   &
+               WRITE(ctmp1,*) 'Some ',cgrid(igrd),' grid points,',   &
                   ' are not boundary points (flagv calculation). Check nbi, nbj, indices for boundary set ',ib_bdy
-               WRITE(ctmp2,*) ' ========== '
-               CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+               CALL ctl_stop( ctmp1 )
             ENDIF
-            CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
+            SELECT CASE( igrd )
+               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
+               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 
+               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 
+            END SELECT 
             DO ib = ibeg, iend
                ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -938,9 +938,9 @@
          DO igrd = 1, jpbgrd
             SELECT CASE( igrd )
-               CASE( 1 )   ;   pmask => bdytmask 
-               CASE( 2 )   ;   pmask => bdyumask 
-               CASE( 3 )   ;   pmask => bdyvmask 
+               CASE( 1 )   ;   zmask => bdytmask 
+               CASE( 2 )   ;   zmask => bdyumask 
+               CASE( 3 )   ;   zmask => bdyvmask 
             END SELECT
-            ztmp(:,:) = 0._wp
+            ztmp(:,:) = -999._wp
             IF( lrim0 ) THEN   ! extent of rim 0
                ibeg = 1                                     ;   iend = idx_bdy(ib_bdy)%nblenrim0(igrd)
@@ -952,8 +952,8 @@
                ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
                IF( ii == 1 .OR. ii == jpi .OR. ij == 1 .OR. ij == jpj )   CYCLE
-               llnon = pmask(ii  ,ij+1) == 1.  
-               llson = pmask(ii  ,ij-1) == 1. 
-               llean = pmask(ii+1,ij  ) == 1. 
-               llwen = pmask(ii-1,ij  ) == 1. 
+               llnon = zmask(ii  ,ij+1) == 1.  
+               llson = zmask(ii  ,ij-1) == 1. 
+               llean = zmask(ii+1,ij  ) == 1. 
+               llwen = zmask(ii-1,ij  ) == 1. 
                inn  = COUNT( (/ llnon, llson, llean, llwen /) )
                IF( inn == 0 ) THEN   ! no neighbours -> interior of a corner  or  cluster of rim points
@@ -961,10 +961,10 @@
                   !  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,   &
+                  IF(     zmask(ii+1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 1.
+                  ELSEIF( zmask(ii-1,ij+1) == 1. ) THEN   ;   ztmp(ii,ij) = 2.
+                  ELSEIF( zmask(ii+1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 3.
+                  ELSEIF( zmask(ii-1,ij-1) == 1. ) THEN   ;   ztmp(ii,ij) = 4.
+                  ELSE                                    ;   ztmp(ii,ij) = -1.
+                     WRITE(ctmp1,*) 'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
                        ' on boundary set ', ib_bdy, ' has no free ocean neighbour'
                      IF( lrim0 ) THEN
@@ -973,6 +973,5 @@
                         WRITE(ctmp2,*) ' There seems to be a cluster of rim 1 points.'
                      END IF
-                     WRITE(ctmp3,*) ' ========== '
-                     CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' )
+                     CALL ctl_warn( ctmp1, ctmp2 )
                   END IF
                END IF
@@ -1005,11 +1004,14 @@
                END IF
                IF( inn == 4 ) THEN
-                  WRITE(ctmp1,*)  ' E R R O R : Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
+                  WRITE(ctmp1,*)  'Problem with  ',cgrid(igrd) ,' grid point', ii, ij,   &
                        ' on boundary set ', ib_bdy, ' have 4 neighbours'
-                  WRITE(ctmp2,*) ' ========== '
-                  CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+                  CALL ctl_stop( ctmp1 )
                END IF
             END DO
-            CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. )
+            SELECT CASE( igrd )
+               CASE( 1 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'T', 1. ) 
+               CASE( 2 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'U', 1. ) 
+               CASE( 3 )   ;   CALL lbc_lnk( 'bdyini', ztmp, 'V', 1. ) 
+            END SELECT 
             DO ib = ibeg, iend
                ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
@@ -1031,5 +1033,5 @@
       !!
       !! ** Method  :  use itreat input to select a case
-      !!               N.B. ntreat is defined for all bdy points in routine bdy_segs
+      !!               N.B. ntreat is defined for all bdy points in routine bdy_rim_treat
       !!
       !!----------------------------------------------------------------------
@@ -1243,15 +1245,15 @@
                      icorns(ib2,1) = npckgw(ib1)
                   ELSEIF ((jpisft(ib2)==jpiwob(ib1)).AND.(jpjwft(ib1)==jpjsob(ib2))) THEN
-                     WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
+                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
                         &                                     jpisft(ib2), jpjwft(ib1)
-                     WRITE(ctmp2,*) ' ==========  Not allowed yet'
-                     WRITE(ctmp3,*) '             Crossing problem with West segment: ',npckgw(ib1), & 
-                        &                                        ' and South segment: ',npckgs(ib2)
-                     CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' )
+                     WRITE(ctmp2,*) ' Not allowed yet'
+                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), & 
+                        &                            ' and South segment: ',npckgs(ib2)
+                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
                   ELSE
-                     WRITE(ctmp1,*) ' E R R O R : Check South and West Open boundary indices'
-                     WRITE(ctmp2,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1) , &
-                        &                                         ' and South segment: ',npckgs(ib2)
-                     CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+                     WRITE(ctmp1,*) ' Check South and West Open boundary indices'
+                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1) , &
+                        &                            ' and South segment: ',npckgs(ib2)
+                     CALL ctl_stop( ctmp1, ctmp2 )
                   END IF
                END IF
@@ -1275,15 +1277,15 @@
                      icorns(ib2,2) = npckge(ib1)
                   ELSEIF ((jpjeft(ib1)==jpjsob(ib2)).AND.(jpisdt(ib2)==jpieob(ib1)+1)) THEN
-                     WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
+                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
                         &                                     jpisdt(ib2), jpjeft(ib1)
-                     WRITE(ctmp2,*) ' ==========  Not allowed yet'
-                     WRITE(ctmp3,*) '             Crossing problem with East segment: ',npckge(ib1), &
-                        &                                        ' and South segment: ',npckgs(ib2)
-                     CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' )
+                     WRITE(ctmp2,*) ' Not allowed yet'
+                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
+                        &                            ' and South segment: ',npckgs(ib2)
+                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
                   ELSE
-                     WRITE(ctmp1,*) ' E R R O R : Check South and East Open boundary indices'
-                     WRITE(ctmp2,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), &
-                     &                                           ' and South segment: ',npckgs(ib2)
-                     CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+                     WRITE(ctmp1,*) ' Check South and East Open boundary indices'
+                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
+                     &                               ' and South segment: ',npckgs(ib2)
+                     CALL ctl_stop( ctmp1, ctmp2 )
                   END IF
                END IF
@@ -1307,15 +1309,15 @@
                      icornn(ib2,1) = npckgw(ib1)
                   ELSEIF ((jpjwdt(ib1)==jpjnob(ib2)+1).AND.(jpinft(ib2)==jpiwob(ib1))) THEN
-                     WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
+                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
                      &                                     jpinft(ib2), jpjwdt(ib1)
-                     WRITE(ctmp2,*) ' ==========  Not allowed yet'
-                     WRITE(ctmp3,*) '             Crossing problem with West segment: ',npckgw(ib1), &
-                     &                                                    ' and North segment: ',npckgn(ib2)
-                     CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' )
+                     WRITE(ctmp2,*) ' Not allowed yet'
+                     WRITE(ctmp3,*) ' Crossing problem with West segment: ',npckgw(ib1), &
+                     &                               ' and North segment: ',npckgn(ib2)
+                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
                   ELSE
-                     WRITE(ctmp1,*) ' E R R O R : Check North and West Open boundary indices'
-                     WRITE(ctmp2,*) ' ==========  Crossing problem with West segment: ',npckgw(ib1), &
-                     &                                                    ' and North segment: ',npckgn(ib2)
-                     CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+                     WRITE(ctmp1,*) ' Check North and West Open boundary indices'
+                     WRITE(ctmp2,*) ' Crossing problem with West segment: ',npckgw(ib1), &
+                     &                               ' and North segment: ',npckgn(ib2)
+                     CALL ctl_stop( ctmp1, ctmp2 )
                   END IF
                END IF
@@ -1339,15 +1341,15 @@
                      icornn(ib2,2) = npckge(ib1)
                   ELSEIF ((jpjedt(ib1)==jpjnob(ib2)+1).AND.(jpindt(ib2)==jpieob(ib1)+1)) THEN
-                     WRITE(ctmp1,*) ' E R R O R : Found an acute open boundary corner at point (i,j)= ', &
+                     WRITE(ctmp1,*) ' Found an acute open boundary corner at point (i,j)= ', &
                      &                                     jpindt(ib2), jpjedt(ib1)
-                     WRITE(ctmp2,*) ' ==========  Not allowed yet'
-                     WRITE(ctmp3,*) '             Crossing problem with East segment: ',npckge(ib1), &
-                     &                                           ' and North segment: ',npckgn(ib2)
-                     CALL ctl_stop( ' ', ctmp1, ctmp2, ctmp3, ' ' )
+                     WRITE(ctmp2,*) ' Not allowed yet'
+                     WRITE(ctmp3,*) ' Crossing problem with East segment: ',npckge(ib1), &
+                     &                               ' and North segment: ',npckgn(ib2)
+                     CALL ctl_stop( ctmp1, ctmp2, ctmp3 )
                   ELSE
-                     WRITE(ctmp1,*) ' E R R O R : Check North and East Open boundary indices'
-                     WRITE(ctmp2,*) ' ==========  Crossing problem with East segment: ',npckge(ib1), &
-                     &                                           ' and North segment: ',npckgn(ib2)
-                     CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+                     WRITE(ctmp1,*) ' Check North and East Open boundary indices'
+                     WRITE(ctmp2,*) ' Crossing problem with East segment: ',npckge(ib1), &
+                     &                               ' and North segment: ',npckgn(ib2)
+                     CALL ctl_stop( ctmp1, ctmp2 )
                   END IF
                END IF
@@ -1375,7 +1377,6 @@
          IF (ztestmask(1)==1) THEN 
             IF (icornw(ib,1)==0) THEN
-               WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgw(ib)
-               WRITE(ctmp2,*) ' ==========  does not start on land or on a corner'                                                  
-               CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
+               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
             ELSE
                ! This is a corner
@@ -1387,7 +1388,6 @@
          IF (ztestmask(2)==1) THEN
             IF (icornw(ib,2)==0) THEN
-               WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgw(ib)
-               WRITE(ctmp2,*) ' ==========  does not end on land or on a corner'                                                  
-               CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+               WRITE(ctmp1,*) ' Open boundary segment ', npckgw(ib)
+               CALL ctl_stop( ' ', ctmp1, ' does not end on land or on a corner' )
             ELSE
                ! This is a corner
@@ -1415,7 +1415,6 @@
          IF (ztestmask(1)==1) THEN
             IF (icorne(ib,1)==0) THEN
-               WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckge(ib)
-               WRITE(ctmp2,*) ' ==========  does not start on land or on a corner'                                                  
-               CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
+               CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
             ELSE
                ! This is a corner
@@ -1427,7 +1426,6 @@
          IF (ztestmask(2)==1) THEN
             IF (icorne(ib,2)==0) THEN
-               WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckge(ib)
-               WRITE(ctmp2,*) ' ==========  does not end on land or on a corner'                                                  
-               CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+               WRITE(ctmp1,*) ' Open boundary segment ', npckge(ib)
+               CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
             ELSE
                ! This is a corner
@@ -1454,12 +1452,10 @@
 
          IF ((ztestmask(1)==1).AND.(icorns(ib,1)==0)) THEN
-            WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgs(ib)
-            WRITE(ctmp2,*) ' ==========  does not start on land or on a corner'                                                  
-            CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
+            CALL ctl_stop( ctmp1, ' does not start on land or on a corner' )
          ENDIF
          IF ((ztestmask(2)==1).AND.(icorns(ib,2)==0)) THEN
-            WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgs(ib)
-            WRITE(ctmp2,*) ' ==========  does not end on land or on a corner'                                                  
-            CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+            WRITE(ctmp1,*) ' Open boundary segment ', npckgs(ib)
+            CALL ctl_stop( ctmp1, ' does not end on land or on a corner' )
          ENDIF
       END DO
@@ -1480,12 +1476,10 @@
 
          IF ((ztestmask(1)==1).AND.(icornn(ib,1)==0)) THEN
-            WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgn(ib)
-            WRITE(ctmp2,*) ' ==========  does not start on land'                                                  
-            CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
+            CALL ctl_stop( ctmp1, ' does not start on land' )
          ENDIF
          IF ((ztestmask(2)==1).AND.(icornn(ib,2)==0)) THEN
-            WRITE(ctmp1,*) ' E R R O R : Open boundary segment ', npckgn(ib)
-            WRITE(ctmp2,*) ' ==========  does not end on land'                                                  
-            CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+            WRITE(ctmp1,*) ' Open boundary segment ', npckgn(ib)
+            CALL ctl_stop( ctmp1, ' does not end on land' )
          ENDIF
       END DO
@@ -1727,11 +1721,79 @@
       !
       IF( itest>0 ) THEN
-         WRITE(ctmp1,*) ' E R R O R : Segments ', ib1, 'and ', ib2
-         WRITE(ctmp2,*) ' ==========  have different open bdy schemes'                                                  
-         CALL ctl_stop( ' ', ctmp1, ctmp2, ' ' )
+         WRITE(ctmp1,*) ' Segments ', ib1, 'and ', ib2
+         CALL ctl_stop( ctmp1, ' have different open bdy schemes' )
       ENDIF
       !
    END SUBROUTINE bdy_ctl_corn
 
+
+   SUBROUTINE bdy_meshwri()
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE bdy_meshwri  ***
+      !!         
+      !! ** Purpose :   write netcdf file with nbr, flagu, flagv, ntreat for T, U 
+      !!                and V points in 2D arrays for easier visualisation/control
+      !!
+      !! ** Method  :   use iom_rstput as in domwri.F
+      !!----------------------------------------------------------------------      
+      INTEGER  ::   ib_bdy, ii, ij, igrd, ib     ! dummy loop indices
+      INTEGER  ::   inum                                   !   -       -
+      REAL(wp), POINTER, DIMENSION(:,:)     ::   zmask                   ! pointer to 2D mask fields
+      REAL(wp)         , DIMENSION(jpi,jpj) ::   ztmp
+      CHARACTER(LEN=1) , DIMENSION(jpbgrd)  ::   cgrid
+      !!----------------------------------------------------------------------      
+      cgrid = (/'t','u','v'/)
+      CALL iom_open( 'bdy_mesh', inum, ldwrt = .TRUE. )
+      DO igrd = 1, jpbgrd
+         SELECT CASE( igrd )
+         CASE( 1 )   ;   zmask => tmask(:,:,1)
+         CASE( 2 )   ;   zmask => umask(:,:,1)
+         CASE( 3 )   ;   zmask => vmask(:,:,1)
+         END SELECT
+         ztmp(:,:) = zmask(:,:)
+         DO ib_bdy = 1, nb_bdy
+            DO ib = 1, idx_bdy(ib_bdy)%nblen(igrd)      ! nbr deined for all rims
+               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%nbr(ib,igrd), wp) + 10.
+               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
+            END DO
+         END DO
+         CALL iom_rstput( 0, 0, inum, 'bdy_nbr_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
+         ztmp(:,:) = zmask(:,:)
+         DO ib_bdy = 1, nb_bdy
+            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagu defined only for rims 0 and 1
+               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagu(ib,igrd), wp) + 10.
+               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
+            END DO
+         END DO
+         CALL iom_rstput( 0, 0, inum, 'flagu_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
+         ztmp(:,:) = zmask(:,:)
+         DO ib_bdy = 1, nb_bdy
+            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! flagv defined only for rims 0 and 1
+               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%flagv(ib,igrd), wp) + 10.
+               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
+            END DO
+         END DO
+         CALL iom_rstput( 0, 0, inum, 'flagv_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
+         ztmp(:,:) = zmask(:,:)
+         DO ib_bdy = 1, nb_bdy
+            DO ib = 1, idx_bdy(ib_bdy)%nblenrim(igrd)   ! ntreat defined only for rims 0 and 1
+               ii = idx_bdy(ib_bdy)%nbi(ib,igrd)
+               ij = idx_bdy(ib_bdy)%nbj(ib,igrd)
+               ztmp(ii,ij) = REAL(idx_bdy(ib_bdy)%ntreat(ib,igrd), wp) + 10.
+               IF( zmask(ii,ij) == 0. ) ztmp(ii,ij) = - ztmp(ii,ij)
+            END DO
+         END DO
+         CALL iom_rstput( 0, 0, inum, 'ntreat_'//cgrid(igrd), ztmp(:,:), ktype = jp_i4 )
+      END DO
+      CALL iom_close( inum )
+
+   END SUBROUTINE bdy_meshwri
+   
    !!=================================================================================
 END MODULE bdyini
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 11257)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/BDY/bdylib.F90	(revision 11258)
@@ -149,7 +149,7 @@
       REAL(wp) ::   zdy_1, zdy_2, zsign_ups
       REAL(wp), PARAMETER :: zepsilon = 1.e-30                 ! local small value
-      REAL(wp), POINTER, DIMENSION(:,:)          :: pmask      ! land/sea mask for field
-      REAL(wp), POINTER, DIMENSION(:,:)          :: pmask_xdif ! land/sea mask for x-derivatives
-      REAL(wp), POINTER, DIMENSION(:,:)          :: pmask_ydif ! land/sea mask for y-derivatives
+      REAL(wp), POINTER, DIMENSION(:,:)          :: zmask      ! land/sea mask for field
+      REAL(wp), POINTER, DIMENSION(:,:)          :: zmask_xdif ! land/sea mask for x-derivatives
+      REAL(wp), POINTER, DIMENSION(:,:)          :: zmask_ydif ! land/sea mask for y-derivatives
       REAL(wp), POINTER, DIMENSION(:,:)          :: pe_xdif    ! scale factors for x-derivatives
       REAL(wp), POINTER, DIMENSION(:,:)          :: pe_ydif    ! scale factors for y-derivatives
@@ -162,7 +162,7 @@
       SELECT CASE(igrd)
          CASE(1)
-            pmask      => tmask(:,:,1)
-            pmask_xdif => umask(:,:,1)
-            pmask_ydif => vmask(:,:,1)
+            zmask      => tmask(:,:,1)
+            zmask_xdif => umask(:,:,1)
+            zmask_ydif => vmask(:,:,1)
             pe_xdif    => e1u(:,:)
             pe_ydif    => e2v(:,:)
@@ -170,7 +170,7 @@
             ij_offset = 0
          CASE(2)
-            pmask      => umask(:,:,1)
-            pmask_xdif => tmask(:,:,1)
-            pmask_ydif => fmask(:,:,1)
+            zmask      => umask(:,:,1)
+            zmask_xdif => tmask(:,:,1)
+            zmask_ydif => fmask(:,:,1)
             pe_xdif    => e1t(:,:)
             pe_ydif    => e2f(:,:)
@@ -178,7 +178,7 @@
             ij_offset = 0
          CASE(3)
-            pmask      => vmask(:,:,1)
-            pmask_xdif => fmask(:,:,1)
-            pmask_ydif => tmask(:,:,1)
+            zmask      => vmask(:,:,1)
+            zmask_xdif => fmask(:,:,1)
+            zmask_ydif => tmask(:,:,1)
             pe_xdif    => e1f(:,:)
             pe_ydif    => e2t(:,:)
@@ -214,12 +214,12 @@
          !
          ! Calculate scale factors for calculation of spatial derivatives.        
-         zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1+ii_offset,ijbm1          )         &
-        &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1          ,ijbm1+ij_offset) ) 
-         zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2+ii_offset,ijbm2          )         &
-        &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2          ,ijbm2+ij_offset) ) 
-         zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )  & 
+         zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1   +ii_offset,ijbm1             )   &
+        &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1             ,ijbm1   +ij_offset) ) 
+         zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2   +ii_offset,ijbm2             )   &
+        &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2             ,ijbm2   +ij_offset) ) 
+         zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )   & 
         &      +  (ijbm1-ijbm1jm1) * pe_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) ) 
-         zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1+ii_offset,ijbm1)                  &
-        &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1          ,ijbm1+ij_offset) ) 
+         zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1   +ii_offset,ijbm1             )   &
+        &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1             ,ijbm1   +ij_offset) ) 
          ! make sure scale factors are nonzero
          if( zey1 .lt. rsmall ) zey1 = zey2
@@ -228,11 +228,11 @@
          zey1 = max(zey1,rsmall); zey2 = max(zey2,rsmall); 
          !
-         ! Calculate masks for calculation of spatial derivatives.        
-         zmask_x = ( abs(iibm1-iibm2) * pmask_xdif(iibm2+ii_offset,ijbm2          )         &
-        &          + abs(ijbm1-ijbm2) * pmask_ydif(iibm2          ,ijbm2+ij_offset) ) 
-         zmask_y1 = ( (iibm1-iibm1jm1) * pmask_xdif(iibm1jm1+ii_offset,ijbm1jm1          )  & 
-        &          +  (ijbm1-ijbm1jm1) * pmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) ) 
-         zmask_y2 = ( (iibm1jp1-iibm1) * pmask_xdif(iibm1+ii_offset,ijbm1)                  &
-        &          +  (ijbm1jp1-ijbm1) * pmask_ydif(iibm1          ,ijbm1+ij_offset) ) 
+         ! Calculate masks for calculation of spatial derivatives.
+         zmask_x  = ( abs(iibm1-iibm2) * zmask_xdif(iibm2   +ii_offset,ijbm2               )   &
+        &           + abs(ijbm1-ijbm2) * zmask_ydif(iibm2             ,ijbm2   +ij_offset) ) 
+         zmask_y1 = ( (iibm1-iibm1jm1) * zmask_xdif(iibm1jm1+ii_offset,ijbm1jm1            )   & 
+        &          +  (ijbm1-ijbm1jm1) * zmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) ) 
+         zmask_y2 = ( (iibm1jp1-iibm1) * zmask_xdif(iibm1   +ii_offset,ijbm1               )   &
+        &          +  (ijbm1jp1-ijbm1) * zmask_ydif(iibm1             ,ijbm1   +ij_offset) ) 
 
          ! Calculation of terms required for both versions of the scheme. 
@@ -242,8 +242,8 @@
          ! Note no rdt factor in expression for zdt because it cancels in the expressions for 
          ! zrx and zry.
-         zdt = phia(iibm1,ijbm1) - phib(iibm1,ijbm1)
-         zdx = ( ( phia(iibm1,ijbm1) - phia(iibm2,ijbm2) ) / zex2 ) * zmask_x 
+         zdt   =     phia(iibm1   ,ijbm1   ) - phib(iibm1   ,ijbm1   )
+         zdx   = ( ( phia(iibm1   ,ijbm1   ) - phia(iibm2   ,ijbm2   ) ) / zex2 ) * zmask_x 
          zdy_1 = ( ( phib(iibm1   ,ijbm1   ) - phib(iibm1jm1,ijbm1jm1) ) / zey1 ) * zmask_y1    
-         zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1)    ) / zey2 ) * zmask_y2 
+         zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1) - phib(iibm1   ,ijbm1   ) ) / zey2 ) * zmask_y2 
          zdy_centred = 0.5 * ( zdy_1 + zdy_2 )
 !!$         zdy_centred = phib(iibm1jp1,ijbm1jp1) - phib(iibm1jm1,ijbm1jm1)
@@ -276,5 +276,5 @@
            &                    + zwgt * ( phi_ext(jb) - phib(ii,ij) ) ) / ( 1. + zrx ) 
          end if
-         phia(ii,ij) = phia(ii,ij) * pmask(ii,ij)
+         phia(ii,ij) = phia(ii,ij) * zmask(ii,ij)
       END DO
       !
@@ -314,7 +314,7 @@
       REAL(wp) ::   zdy_1, zdy_2,  zsign_ups
       REAL(wp), PARAMETER :: zepsilon = 1.e-30                 ! local small value
-      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask      ! land/sea mask for field
-      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask_xdif ! land/sea mask for x-derivatives
-      REAL(wp), POINTER, DIMENSION(:,:,:)        :: pmask_ydif ! land/sea mask for y-derivatives
+      REAL(wp), POINTER, DIMENSION(:,:,:)        :: zmask      ! land/sea mask for field
+      REAL(wp), POINTER, DIMENSION(:,:,:)        :: zmask_xdif ! land/sea mask for x-derivatives
+      REAL(wp), POINTER, DIMENSION(:,:,:)        :: zmask_ydif ! land/sea mask for y-derivatives
       REAL(wp), POINTER, DIMENSION(:,:)          :: pe_xdif    ! scale factors for x-derivatives
       REAL(wp), POINTER, DIMENSION(:,:)          :: pe_ydif    ! scale factors for y-derivatives
@@ -327,7 +327,7 @@
       SELECT CASE(igrd)
          CASE(1)
-            pmask      => tmask(:,:,:)
-            pmask_xdif => umask(:,:,:)
-            pmask_ydif => vmask(:,:,:)
+            zmask      => tmask(:,:,:)
+            zmask_xdif => umask(:,:,:)
+            zmask_ydif => vmask(:,:,:)
             pe_xdif    => e1u(:,:)
             pe_ydif    => e2v(:,:)
@@ -335,7 +335,7 @@
             ij_offset = 0
          CASE(2)
-            pmask      => umask(:,:,:)
-            pmask_xdif => tmask(:,:,:)
-            pmask_ydif => fmask(:,:,:)
+            zmask      => umask(:,:,:)
+            zmask_xdif => tmask(:,:,:)
+            zmask_ydif => fmask(:,:,:)
             pe_xdif    => e1t(:,:)
             pe_ydif    => e2f(:,:)
@@ -343,7 +343,7 @@
             ij_offset = 0
          CASE(3)
-            pmask      => vmask(:,:,:)
-            pmask_xdif => fmask(:,:,:)
-            pmask_ydif => tmask(:,:,:)
+            zmask      => vmask(:,:,:)
+            zmask_xdif => fmask(:,:,:)
+            zmask_ydif => tmask(:,:,:)
             pe_xdif    => e1f(:,:)
             pe_ydif    => e2t(:,:)
@@ -381,12 +381,12 @@
             !
             ! Calculate scale factors for calculation of spatial derivatives.        
-            zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1+ii_offset,ijbm1          )         &
-           &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1          ,ijbm1+ij_offset) ) 
-            zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2+ii_offset,ijbm2          )         &
-           &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2          ,ijbm2+ij_offset) ) 
-            zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )  & 
+            zex1 = ( abs(iibm1-iibm2) * pe_xdif(iibm1   +ii_offset,ijbm1             )   &
+           &       + abs(ijbm1-ijbm2) * pe_ydif(iibm1             ,ijbm1+ij_offset   ) ) 
+            zex2 = ( abs(iibm1-iibm2) * pe_xdif(iibm2   +ii_offset,ijbm2             )   &
+           &       + abs(ijbm1-ijbm2) * pe_ydif(iibm2             ,ijbm2+ij_offset   ) ) 
+            zey1 = ( (iibm1-iibm1jm1) * pe_xdif(iibm1jm1+ii_offset,ijbm1jm1          )   & 
            &      +  (ijbm1-ijbm1jm1) * pe_ydif(iibm1jm1          ,ijbm1jm1+ij_offset) ) 
-            zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1+ii_offset,ijbm1)                  &
-           &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1          ,ijbm1+ij_offset) ) 
+            zey2 = ( (iibm1jp1-iibm1) * pe_xdif(iibm1   +ii_offset,ijbm1             )   &
+           &      +  (ijbm1jp1-ijbm1) * pe_ydif(iibm1             ,ijbm1+ij_offset   ) ) 
             ! make sure scale factors are nonzero
             if( zey1 .lt. rsmall ) zey1 = zey2
@@ -396,10 +396,10 @@
             !
             ! Calculate masks for calculation of spatial derivatives.        
-            zmask_x = ( abs(iibm1-iibm2) * pmask_xdif(iibm2+ii_offset,ijbm2          ,jk)          &
-           &          + abs(ijbm1-ijbm2) * pmask_ydif(iibm2          ,ijbm2+ij_offset,jk) ) 
-            zmask_y1 = ( (iibm1-iibm1jm1) * pmask_xdif(iibm1jm1+ii_offset,ijbm1jm1          ,jk)   & 
-           &          +  (ijbm1-ijbm1jm1) * pmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset,jk) ) 
-            zmask_y2 = ( (iibm1jp1-iibm1) * pmask_xdif(iibm1+ii_offset,ijbm1          ,jk)         &
-           &          +  (ijbm1jp1-ijbm1) * pmask_ydif(iibm1          ,ijbm1+ij_offset,jk) ) 
+            zmask_x  = ( abs(iibm1-iibm2) * zmask_xdif(iibm2   +ii_offset,ijbm2             ,jk)   &
+           &           + abs(ijbm1-ijbm2) * zmask_ydif(iibm2             ,ijbm2   +ij_offset,jk) ) 
+            zmask_y1 = ( (iibm1-iibm1jm1) * zmask_xdif(iibm1jm1+ii_offset,ijbm1jm1          ,jk)   & 
+           &          +  (ijbm1-ijbm1jm1) * zmask_ydif(iibm1jm1          ,ijbm1jm1+ij_offset,jk) ) 
+            zmask_y2 = ( (iibm1jp1-iibm1) * zmask_xdif(iibm1   +ii_offset,ijbm1             ,jk)   &
+           &          +  (ijbm1jp1-ijbm1) * zmask_ydif(iibm1             ,ijbm1   +ij_offset,jk) ) 
             !
             ! Calculate normal (zrx) and tangential (zry) components of radiation velocities.
@@ -407,6 +407,6 @@
             ! Centred derivative is calculated as average of "left" and "right" derivatives for 
             ! this reason. 
-            zdt = phia(iibm1,ijbm1,jk) - phib(iibm1,ijbm1,jk)
-            zdx = ( ( phia(iibm1,ijbm1,jk) - phia(iibm2,ijbm2,jk) ) / zex2 ) * zmask_x                  
+            zdt   =     phia(iibm1   ,ijbm1   ,jk) - phib(iibm1   ,ijbm1   ,jk)
+            zdx   = ( ( phia(iibm1   ,ijbm1   ,jk) - phia(iibm2   ,ijbm2   ,jk) ) / zex2 ) * zmask_x                  
             zdy_1 = ( ( phib(iibm1   ,ijbm1   ,jk) - phib(iibm1jm1,ijbm1jm1,jk) ) / zey1 ) * zmask_y1  
             zdy_2 = ( ( phib(iibm1jp1,ijbm1jp1,jk) - phib(iibm1   ,ijbm1   ,jk) ) / zey2 ) * zmask_y2      
@@ -442,5 +442,5 @@
               &                       + zwgt * ( phi_ext(jb,jk) - phib(ii,ij,jk) ) ) / ( 1. + zrx ) 
             end if
-            phia(ii,ij,jk) = phia(ii,ij,jk) * pmask(ii,ij,jk)
+            phia(ii,ij,jk) = phia(ii,ij,jk) * zmask(ii,ij,jk)
          END DO
          !
@@ -469,5 +469,5 @@
       !! 
       REAL(wp) ::   zweight
-      REAL(wp), POINTER, DIMENSION(:,:,:)      :: pmask         ! land/sea mask for field
+      REAL(wp), POINTER, DIMENSION(:,:,:)      :: zmask         ! land/sea mask for field
       INTEGER  ::   ib, ik   ! dummy loop indices
       INTEGER  ::   ii, ij   ! 2D addresses
@@ -480,7 +480,7 @@
       !
       SELECT CASE(igrd)
-         CASE(1)   ;   pmask => tmask(:,:,:)
-         CASE(2)   ;   pmask => umask(:,:,:)
-         CASE(3)   ;   pmask => vmask(:,:,:)
+         CASE(1)   ;   zmask => tmask(:,:,:)
+         CASE(2)   ;   zmask => umask(:,:,:)
+         CASE(3)   ;   zmask => vmask(:,:,:)
          CASE DEFAULT ;   CALL ctl_stop( 'unrecognised value for igrd in bdy_nmn' )
       END SELECT
@@ -502,5 +502,5 @@
             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)  
+               IF( zmask(ii1,ij1,ik) /= 0. )   phia(ii,ij,ik) = phia(ii1,ij1,ik)  
             END DO
          CASE( 9:12 )
@@ -508,5 +508,5 @@
             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)
+               zweight = zmask(ii1,ij1,ik) + zmask(ii2,ij2,ik)
                IF( zweight /= 0. )   phia(ii,ij,ik) = ( phia(ii1,ij1,ik) + phia(ii2,ij2,ik) ) / zweight
             END DO
@@ -516,5 +516,5 @@
             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)
+               zweight = zmask(ii1,ij1,ik) + zmask(ii2,ij2,ik) + zmask(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
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DOM/domain.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DOM/domain.F90	(revision 11257)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DOM/domain.F90	(revision 11258)
@@ -101,5 +101,5 @@
          CASE( 0 )   ;   WRITE(numout,*) '         (i.e. closed)'
          CASE( 1 )   ;   WRITE(numout,*) '         (i.e. cyclic east-west)'
-         CASE( 2 )   ;   WRITE(numout,*) '         (i.e. equatorial symmetric)'
+         CASE( 2 )   ;   WRITE(numout,*) '         (i.e. cyclic north-south)'
          CASE( 3 )   ;   WRITE(numout,*) '         (i.e. north fold with T-point pivot)'
          CASE( 4 )   ;   WRITE(numout,*) '         (i.e. cyclic east-west and north fold with T-point pivot)'
@@ -527,5 +527,6 @@
       INTEGER ::   inum, ii   ! local integer
       REAL(wp) ::   zorca_res                     ! local scalars
-      REAL(wp) ::   ziglo, zjglo, zkglo, zperio   !   -      -
+      REAL(wp) ::   zperio                        !   -      -
+      INTEGER, DIMENSION(4) ::   idvar, idimsz    ! size   of dimensions
       !!----------------------------------------------------------------------
       !
@@ -559,8 +560,9 @@
          !
       ENDIF
-      !
-      CALL iom_get( inum, 'jpiglo', ziglo  )   ;   kpi = NINT( ziglo )
-      CALL iom_get( inum, 'jpjglo', zjglo  )   ;   kpj = NINT( zjglo )
-      CALL iom_get( inum, 'jpkglo', zkglo  )   ;   kpk = NINT( zkglo )
+       !
+      idvar = iom_varid( inum, 'e3t_0', kdimsz = idimsz )   ! use e3t_0, that must exist, to get jp(ijk)glo
+      kpi = idimsz(1)
+      kpj = idimsz(2)
+      kpk = idimsz(3)
       CALL iom_get( inum, 'jperio', zperio )   ;   kperio = NINT( zperio )
       CALL iom_close( inum )
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynkeg.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynkeg.F90	(revision 11257)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/DYN/dynkeg.F90	(revision 11258)
@@ -74,11 +74,8 @@
       INTEGER, INTENT( in ) ::   kscheme   ! =0/1   type of KEG scheme 
       !
-      INTEGER  ::   ji, jj, jk, jb           ! dummy loop indices
-      INTEGER  ::   igrd, ib_bdy             ! local integers
+      INTEGER  ::   ji, jj, jk             ! dummy loop indices
       REAL(wp) ::   zu, zv                   ! local scalars
       REAL(wp), DIMENSION(jpi,jpj,jpk)        ::   zhke
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   ztrdu, ztrdv 
-      REAL(wp)  :: zweightu, zweightv
-      LOGICAL, DIMENSION(4) :: llsend1, llrecv1  ! indicate how bdy communications are to be carried out
       !!----------------------------------------------------------------------
       !
@@ -113,42 +110,4 @@
             END DO
          END DO
-         !
-         IF (ln_bdy) THEN
-            ! Maria Luneva & Fred Wobus: July-2016
-            ! compensate for lack of turbulent kinetic energy on liquid bdy points
-            DO ib_bdy = 1, nb_bdy
-               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
-                  igrd = 1           ! compensating null velocity on the bdy
-                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
-                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 1 to jpi
-                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 1 to jpj
-                     IF( ji == 1 .OR. jj == 1 )  CYCLE
-                     DO jk = 1, jpkm1
-                        zhke(ji,jj,jk) = 0._wp
-                        zweightu = umask(ji-1,jj  ,jk) + umask(ji,jj,jk)
-                        zweightv = vmask(ji  ,jj-1,jk) + vmask(ji,jj,jk)
-                        zu = un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)  +  un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk)
-                        zv = vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)  +  vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk)
-                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zu / (2._wp * zweightu) 
-                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) + zv / (2._wp * zweightv) 
-                     END DO
-                  END DO
-               END IF
-            END DO
-            ! send jpi-1, jpj-1 and receive 1 used in the computation of the speed tendencies
-            llsend1(:) = .false.
-            llrecv1(:) = .false.
-            DO ib_bdy = 1, nb_bdy
-               llsend1(2) = llsend1(2) .OR. ANY(lsend_bdy(ib_bdy,igrd,2,:))   ! send east
-               llsend1(4) = llsend1(4) .OR. ANY(lsend_bdy(ib_bdy,igrd,4,:))   ! send north
-               llrecv1(1) = llrecv1(1) .OR. ANY(lrecv_bdy(ib_bdy,igrd,1,:))   ! receive west 
-               llrecv1(3) = llrecv1(3) .OR. ANY(lrecv_bdy(ib_bdy,igrd,3,:))   ! receive south
-            END DO
-   
-            IF( ANY(llsend1) .OR. ANY(llrecv1) ) THEN   ! if need to send/recv in at least one direction
-               CALL lbc_lnk( 'bdydyn2d', zhke, 'T',  1., kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1 )
-            END IF
-         END IF
-         !
       CASE ( nkeg_HW )                          !--  Hollingsworth scheme  --!
          DO jk = 1, jpkm1
@@ -168,35 +127,4 @@
             END DO
          END DO
-         IF (ln_bdy) THEN
-            ! Maria Luneva & Fred Wobus: July-2016
-            ! compensate for lack of turbulent kinetic energy on liquid bdy points
-            DO ib_bdy = 1, nb_bdy
-               IF( cn_dyn3d(ib_bdy) /= 'none' ) THEN
-                  igrd = 1           ! compensation null velocity on land at the bdy
-                  DO jb = 1, idx_bdy(ib_bdy)%nblenrim(igrd)
-                     ji   = idx_bdy(ib_bdy)%nbi(jb,igrd)   ! maximum extent : from 1 to jpi
-                     jj   = idx_bdy(ib_bdy)%nbj(jb,igrd)   ! maximum extent : from 1 to jpj
-                     IF( ji == 1 .OR. ji == jpi .OR. jj == 1 .OR. jj == jpj )  CYCLE
-                     DO jk = 1, jpkm1
-                        zhke(ji,jj,jk) = 0._wp
-                        zweightu = 8._wp * ( umask(ji-1,jj  ,jk) + umask(ji  ,jj  ,jk) ) &
-                             &   + 2._wp * ( umask(ji-1,jj-1,jk) + umask(ji-1,jj+1,jk) + umask(ji  ,jj-1,jk) + umask(ji  ,jj+1,jk) ) 
-                        zweightv = 8._wp * ( vmask(ji  ,jj-1,jk) + vmask(ji  ,jj-1,jk) ) &
-                             &   + 2._wp * ( vmask(ji-1,jj-1,jk) + vmask(ji+1,jj-1,jk) + vmask(ji-1,jj  ,jk) + vmask(ji+1,jj  ,jk) )
-                        zu = 8._wp * ( un(ji-1,jj  ,jk) * un(ji-1,jj  ,jk)    &
-                           &         + un(ji  ,jj  ,jk) * un(ji  ,jj  ,jk) )  &
-                           &   +     ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) ) * ( un(ji-1,jj-1,jk) + un(ji-1,jj+1,jk) )   &
-                           &   +     ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) ) * ( un(ji  ,jj-1,jk) + un(ji  ,jj+1,jk) )
-                        zv = 8._wp * ( vn(ji  ,jj-1,jk) * vn(ji  ,jj-1,jk)    &
-                           &         + vn(ji  ,jj  ,jk) * vn(ji  ,jj  ,jk) )  &
-                           &  +      ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) ) * ( vn(ji-1,jj-1,jk) + vn(ji+1,jj-1,jk) )   &
-                           &  +      ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) ) * ( vn(ji-1,jj  ,jk) + vn(ji+1,jj  ,jk) )
-                        IF( zweightu > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zu / ( 2._wp * zweightu )
-                        IF( zweightv > 0._wp )   zhke(ji,jj,jk) =  zhke(ji,jj,jk) +  zv / ( 2._wp * zweightv )
-                     END DO
-                  END DO
-               END IF
-            END DO
-         END IF
          CALL lbc_lnk( 'dynkeg', zhke, 'T', 1. )
          !
Index: /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/fldread.F90
===================================================================
--- /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/fldread.F90	(revision 11257)
+++ /NEMO/branches/2019/dev_r10984_HPC-13_IRRMANN_BDY_optimization/src/OCE/SBC/fldread.F90	(revision 11258)
@@ -465,5 +465,5 @@
       ELSE                                      ;   it_offset = 0
       ENDIF
-      IF( PRESENT(kt_offset) )   it_offset = kt_offset
+      IF( PRESENT(kt_offset) )      it_offset = kt_offset
       IF( PRESENT(kit) ) THEN   ;   it_offset = ( kit + it_offset ) * NINT( rdt/REAL(nn_baro,wp) )
       ELSE                      ;   it_offset =         it_offset   * NINT(       rdt            )
@@ -486,6 +486,6 @@
             !       forcing record :    1 
             !                            
-            ztmp = REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 &
-           &       + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday )
+            ztmp =  REAL( nsec_year, wp ) / ( REAL( nyear_len(1), wp ) * rday ) + 0.5 &
+               &  + REAL( it_offset, wp ) / ( REAL( nyear_len(1), wp ) * rday )
             sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
             ! swap at the middle of the year
@@ -517,6 +517,6 @@
             !       forcing record :  nmonth 
             !                            
-            ztmp = REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 &
-           &       + REAL( it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday )
+            ztmp =  REAL( nsec_month, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday ) + 0.5 &
+           &      + REAL(  it_offset, wp ) / ( REAL( nmonth_len(nmonth), wp ) * rday )
             imth = nmonth + INT( ztmp ) - COUNT((/llbefore/))
             IF( sdjf%cltype == 'monthly' ) THEN   ;   sdjf%nrec_a(1) = 1 + INT( ztmp ) - COUNT((/llbefore/))
