Index: NEMO/branches/2020/tickets_icb_1900/src/OCE/LBC/mppini.F90
===================================================================
--- NEMO/branches/2020/tickets_icb_1900/src/OCE/LBC/mppini.F90	(revision 13277)
+++ NEMO/branches/2020/tickets_icb_1900/src/OCE/LBC/mppini.F90	(revision 13899)
@@ -8,6 +8,6 @@
    !!            8.0  !  1998-05  (M. Imbard, J. Escobar, L. Colombet )  SHMEM and MPI versions
    !!  NEMO      1.0  !  2004-01  (G. Madec, J.M Molines)  F90 : free form , north fold jpni > 1
-   !!            3.4  ! 2011-10  (A. C. Coward, NOCS & J. Donners, PRACE) add mpp_init_nfdcom
-   !!            3.   ! 2013-06  (I. Epicoco, S. Mocavero, CMCC) mpp_init_nfdcom: setup avoiding MPI communication 
+   !!            3.4  !  2011-10  (A. C. Coward, NOCS & J. Donners, PRACE)  add init_nfdcom
+   !!            3.   !  2013-06  (I. Epicoco, S. Mocavero, CMCC)  init_nfdcom: setup avoiding MPI communication 
    !!            4.0  !  2016-06  (G. Madec)  use domain configuration file instead of bathymetry file
    !!            4.0  !  2017-06  (J.M. Molines, T. Lovato) merge of mppini and mppini_2
@@ -15,15 +15,13 @@
 
    !!----------------------------------------------------------------------
-   !!  mpp_init          : Lay out the global domain over processors with/without land processor elimination
-   !!  mpp_init_mask     : Read global bathymetric information to facilitate land suppression
-   !!  mpp_init_ioipsl   : IOIPSL initialization in mpp 
-   !!  mpp_init_partition: Calculate MPP domain decomposition
-   !!  factorise         : Calculate the factors of the no. of MPI processes
-   !!  mpp_init_nfdcom   : Setup for north fold exchanges with explicit point-to-point messaging
+   !!  mpp_init       : Lay out the global domain over processors with/without land processor elimination
+   !!      init_ioipsl: IOIPSL initialization in mpp 
+   !!      init_nfdcom: Setup for north fold exchanges with explicit point-to-point messaging
+   !!      init_doloop: set the starting/ending indices of DO-loop used in do_loop_substitute 
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
    USE bdy_oce        ! open BounDarY  
    !
-   USE lbcnfd  , ONLY : isendto, nsndto, nfsloop, nfeloop   ! Setup of north fold exchanges 
+   USE lbcnfd  , ONLY : isendto, nsndto ! Setup of north fold exchanges 
    USE lib_mpp        ! distribued memory computing library
    USE iom            ! nemo I/O library 
@@ -34,8 +32,11 @@
    PRIVATE
 
-   PUBLIC mpp_init       ! called by opa.F90
-
-   INTEGER :: numbot = -1  ! 'bottom_level' local logical unit
-   INTEGER :: numbdy = -1  ! 'bdy_msk'      local logical unit
+   PUBLIC   mpp_init       ! called by nemogcm.F90
+   PUBLIC   mpp_getnum     ! called by prtctl
+   PUBLIC   mpp_basesplit  ! called by prtctl
+   PUBLIC   mpp_is_ocean   ! called by prtctl
+   
+   INTEGER ::   numbot = -1   ! 'bottom_level' local logical unit
+   INTEGER ::   numbdy = -1   ! 'bdy_msk'      local logical unit
    
    !!----------------------------------------------------------------------
@@ -61,4 +62,7 @@
       !!----------------------------------------------------------------------
       !
+      nn_hls = 1
+      jpiglo = Ni0glo + 2 * nn_hls
+      jpjglo = Nj0glo + 2 * nn_hls
       jpimax = jpiglo
       jpjmax = jpjglo
@@ -66,19 +70,13 @@
       jpj    = jpjglo
       jpk    = jpkglo
-      jpim1  = jpi-1                                            ! inner domain indices
-      jpjm1  = jpj-1                                            !   "           "
-      jpkm1  = MAX( 1, jpk-1 )                                  !   "           "
+      jpim1  = jpi-1                         ! inner domain indices
+      jpjm1  = jpj-1                         !   "           "
+      jpkm1  = MAX( 1, jpk-1 )               !   "           "
       jpij   = jpi*jpj
       jpni   = 1
       jpnj   = 1
       jpnij  = jpni*jpnj
-      nimpp  = 1           ! 
+      nimpp  = 1
       njmpp  = 1
-      nlci   = jpi
-      nlcj   = jpj
-      nldi   = 1
-      nldj   = 1
-      nlei   = jpi
-      nlej   = jpj
       nbondi = 2
       nbondj = 2
@@ -90,4 +88,6 @@
       l_Jperio = jpnj == 1 .AND. (jperio == 2 .OR. jperio == 7)
       !
+      CALL init_doloop                       ! set start/end indices or do-loop depending on the halo width value (nn_hls) 
+      !
       IF(lwp) THEN
          WRITE(numout,*)
@@ -98,8 +98,4 @@
       ENDIF
       !
-      IF(  jpni /= 1 .OR. jpnj /= 1 .OR. jpnij /= 1 )                                     &
-         CALL ctl_stop( 'mpp_init: equality  jpni = jpnj = jpnij = 1 is not satisfied',   &
-            &           'the domain is lay out for distributed memory computing!' )
-         !
 #if defined key_agrif
     IF (.NOT.agrif_root()) THEN
@@ -135,6 +131,4 @@
       !!                    njmpp     : latitudinal  index
       !!                    narea     : number for local area
-      !!                    nlci      : first dimension
-      !!                    nlcj      : second dimension
       !!                    nbondi    : mark for "east-west local boundary"
       !!                    nbondj    : mark for "north-south local boundary"
@@ -147,7 +141,6 @@
       INTEGER ::   ji, jj, jn, jproc, jarea   ! dummy loop indices
       INTEGER ::   inijmin
-      INTEGER ::   i2add
       INTEGER ::   inum                       ! local logical unit
-      INTEGER ::   idir, ifreq, icont         ! local integers
+      INTEGER ::   idir, ifreq                ! local integers
       INTEGER ::   ii, il1, ili, imil         !   -       -
       INTEGER ::   ij, il2, ilj, ijm1         !   -       -
@@ -162,8 +155,8 @@
       INTEGER, ALLOCATABLE, DIMENSION(:)     ::   iin, ii_nono, ii_noea          ! 1D workspace
       INTEGER, ALLOCATABLE, DIMENSION(:)     ::   ijn, ii_noso, ii_nowe          !  -     -
-      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ilci, ibondi, ipproc   ! 2D workspace
-      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ilcj, ibondj, ipolj    !  -     -
-      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ilei, ildi, iono, ioea         !  -     -
-      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ilej, ildj, ioso, iowe         !  -     -
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ijpi, ibondi, ipproc   ! 2D workspace
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ijpj, ibondj, ipolj    !  -     -
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iie0, iis0, iono, ioea         !  -     -
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ije0, ijs0, ioso, iowe         !  -     -
       LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   llisoce                        !  -     -
       NAMELIST/nambdy/ ln_bdy, nb_bdy, ln_coords_file, cn_coords_file,           &
@@ -173,5 +166,5 @@
            &             cn_ice, nn_ice_dta,                                     &
            &             ln_vol, nn_volctl, nn_rimwidth
-      NAMELIST/nammpp/ jpni, jpnj, ln_nnogather, ln_listonly
+      NAMELIST/nammpp/ jpni, jpnj, nn_hls, ln_nnogather, ln_listonly
       !!----------------------------------------------------------------------
       !
@@ -186,4 +179,5 @@
 902   IF( ios >  0 )   CALL ctl_nam ( ios , 'nammpp in configuration namelist' )   
       !
+      nn_hls = MAX(1, nn_hls)   ! nn_hls must be > 0
       IF(lwp) THEN
             WRITE(numout,*) '   Namelist nammpp'
@@ -195,8 +189,15 @@
          ENDIF
             WRITE(numout,*) '      avoid use of mpi_allgather at the north fold  ln_nnogather = ', ln_nnogather
+            WRITE(numout,*) '      halo width (applies to both rows and columns)       nn_hls = ', nn_hls
       ENDIF
       !
       IF(lwm)   WRITE( numond, nammpp )
-
+      !
+!!!------------------------------------
+!!!  nn_hls shloud be read in nammpp
+!!!------------------------------------
+      jpiglo = Ni0glo + 2 * nn_hls
+      jpjglo = Nj0glo + 2 * nn_hls
+      !
       ! do we need to take into account bdy_msk?
       READ  ( numnam_ref, nambdy, IOSTAT = ios, ERR = 903)
@@ -208,10 +209,10 @@
       IF( ln_bdy .AND. ln_mask_file ) CALL iom_open( cn_mask_file, numbdy )
       !
-      IF( ln_listonly )   CALL mpp_init_bestpartition( MAX(mppsize,jpni*jpnj), ldlist = .TRUE. )   ! must be done by all core
+      IF( ln_listonly )   CALL bestpartition( MAX(mppsize,jpni*jpnj), ldlist = .TRUE. )   ! must be done by all core
       !
       !  1. Dimension arrays for subdomains
       ! -----------------------------------
       !
-      ! If dimensions of processor grid weren't specified in the namelist file
+      ! If dimensions of processors grid weren't specified in the namelist file
       ! then we calculate them here now that we have our communicator size
       IF(lwp) THEN
@@ -221,14 +222,14 @@
       ENDIF
       IF( jpni < 1 .OR. jpnj < 1 ) THEN
-         CALL mpp_init_bestpartition( mppsize, jpni, jpnj )           ! best mpi decomposition for mppsize mpi processes
+         CALL bestpartition( mppsize, jpni, jpnj )           ! best mpi decomposition for mppsize mpi processes
          llauto = .TRUE.
          llbest = .TRUE.
       ELSE
          llauto = .FALSE.
-         CALL mpp_init_bestpartition( mppsize, inbi, inbj, icnt2 )    ! best mpi decomposition for mppsize mpi processes
+         CALL bestpartition( mppsize, inbi, inbj, icnt2 )    ! best mpi decomposition for mppsize mpi processes
          ! largest subdomain size for mpi decoposition jpni*jpnj given in the namelist
-         CALL mpp_basic_decomposition( jpni, jpnj, jpimax, jpjmax )
-         ! largest subdomain size for mpi decoposition inbi*inbj given by mpp_init_bestpartition
-         CALL mpp_basic_decomposition( inbi, inbj,  iimax,  ijmax )
+         CALL mpp_basesplit( jpiglo, jpjglo, nn_hls, jpni, jpnj, jpimax, jpjmax )
+         ! largest subdomain size for mpi decoposition inbi*inbj given by bestpartition
+         CALL mpp_basesplit( jpiglo, jpjglo, nn_hls, inbi, inbj,  iimax,  ijmax )
          icnt1 = jpni*jpnj - mppsize   ! number of land subdomains that should be removed to use mppsize mpi processes
          IF(lwp) THEN
@@ -261,5 +262,5 @@
       ! look for land mpi subdomains...
       ALLOCATE( llisoce(jpni,jpnj) )
-      CALL mpp_init_isoce( jpni, jpnj, llisoce )
+      CALL mpp_is_ocean( llisoce )
       inijmin = COUNT( llisoce )   ! number of oce subdomains
 
@@ -270,5 +271,5 @@
          WRITE(ctmp4,*) '   ==>>> There is the list of best domain decompositions you should use: '
          CALL ctl_stop( ctmp1, ctmp2, ctmp3, ' ', ctmp4, ' ' )
-         CALL mpp_init_bestpartition( mppsize, ldlist = .TRUE. )   ! must be done by all core
+         CALL bestpartition( mppsize, ldlist = .TRUE. )   ! must be done by all core
       ENDIF
 
@@ -294,5 +295,5 @@
             WRITE(numout,*)
          ENDIF
-         CALL mpp_init_bestpartition( mppsize, ldlist = .TRUE. )   ! must be done by all core
+         CALL bestpartition( mppsize, ldlist = .TRUE. )   ! must be done by all core
       ENDIF
 
@@ -319,17 +320,14 @@
 9003  FORMAT (a, i5)
 
-      IF( numbot /= -1 )   CALL iom_close( numbot )
-      IF( numbdy /= -1 )   CALL iom_close( numbdy )
-    
-      ALLOCATE(  nfiimpp(jpni,jpnj), nfipproc(jpni,jpnj), nfilcit(jpni,jpnj) ,    &
-         &       nimppt(jpnij) , ibonit(jpnij) , nlcit(jpnij) , nlcjt(jpnij) ,    &
-         &       njmppt(jpnij) , ibonjt(jpnij) , nldit(jpnij) , nldjt(jpnij) ,    &
-         &                                       nleit(jpnij) , nlejt(jpnij) ,    &
+      ALLOCATE(  nfimpp(jpni ) , nfproc(jpni ) ,   nfjpi(jpni ) ,                     &
+         &       nimppt(jpnij) , ibonit(jpnij) ,  jpiall(jpnij) ,  jpjall(jpnij) ,    &
+         &       njmppt(jpnij) , ibonjt(jpnij) , nis0all(jpnij) , njs0all(jpnij) ,    &
+         &                                       nie0all(jpnij) , nje0all(jpnij) ,    &
          &       iin(jpnij), ii_nono(jpnij), ii_noea(jpnij),   &
          &       ijn(jpnij), ii_noso(jpnij), ii_nowe(jpnij),   &
-         &       iimppt(jpni,jpnj), ilci(jpni,jpnj), ibondi(jpni,jpnj), ipproc(jpni,jpnj),   &
-         &       ijmppt(jpni,jpnj), ilcj(jpni,jpnj), ibondj(jpni,jpnj), ipolj(jpni,jpnj),   &
-         &       ilei(jpni,jpnj), ildi(jpni,jpnj), iono(jpni,jpnj), ioea(jpni,jpnj),   &
-         &       ilej(jpni,jpnj), ildj(jpni,jpnj), ioso(jpni,jpnj), iowe(jpni,jpnj),   &
+         &       iimppt(jpni,jpnj), ijpi(jpni,jpnj), ibondi(jpni,jpnj), ipproc(jpni,jpnj),   &
+         &       ijmppt(jpni,jpnj), ijpj(jpni,jpnj), ibondj(jpni,jpnj),  ipolj(jpni,jpnj),   &
+         &         iie0(jpni,jpnj), iis0(jpni,jpnj),   iono(jpni,jpnj),   ioea(jpni,jpnj),   &
+         &         ije0(jpni,jpnj), ijs0(jpni,jpnj),   ioso(jpni,jpnj),   iowe(jpni,jpnj),   &
          &       STAT=ierr )
       CALL mpp_sum( 'mppini', ierr )
@@ -345,9 +343,18 @@
       ! -----------------------------------
       !
-      nreci = 2 * nn_hls
-      nrecj = 2 * nn_hls
-      CALL mpp_basic_decomposition( jpni, jpnj, jpimax, jpjmax, iimppt, ijmppt, ilci, ilcj )
-      nfiimpp(:,:) = iimppt(:,:)
-      nfilcit(:,:) = ilci(:,:)
+      CALL mpp_basesplit( jpiglo, jpjglo, nn_hls, jpni, jpnj, jpimax, jpjmax, iimppt, ijmppt, ijpi, ijpj )
+      CALL mpp_getnum( llisoce, ipproc, iin, ijn )
+      !
+      !DO jn = 1, jpni
+      !   jproc = ipproc(jn,jpnj)
+      !   ii = iin(jproc+1)
+      !   ij = ijn(jproc+1)
+      !   nfproc(jn) = jproc
+      !   nfimpp(jn) = iimppt(ii,ij)
+      !   nfjpi (jn) =   ijpi(ii,ij)
+      !END DO
+      nfproc(:) = ipproc(:,jpnj) 
+      nfimpp(:) = iimppt(:,jpnj) 
+      nfjpi (:) =   ijpi(:,jpnj)
       !
       IF(lwp) THEN
@@ -358,7 +365,8 @@
          WRITE(numout,*) '      jpni = ', jpni  
          WRITE(numout,*) '      jpnj = ', jpnj
+         WRITE(numout,*) '     jpnij = ', jpnij
          WRITE(numout,*)
-         WRITE(numout,*) '      sum ilci(i,1) = ', sum(ilci(:,1)), ' jpiglo = ', jpiglo
-         WRITE(numout,*) '      sum ilcj(1,j) = ', sum(ilcj(1,:)), ' jpjglo = ', jpjglo
+         WRITE(numout,*) '      sum ijpi(i,1) = ', sum(ijpi(:,1)), ' jpiglo = ', jpiglo
+         WRITE(numout,*) '      sum ijpj(1,j) = ', sum(ijpj(1,:)), ' jpjglo = ', jpjglo
       ENDIF
      
@@ -375,6 +383,6 @@
          ii = 1 + MOD(iarea0,jpni)
          ij = 1 +     iarea0/jpni
-         ili = ilci(ii,ij)
-         ilj = ilcj(ii,ij)
+         ili = ijpi(ii,ij)
+         ilj = ijpj(ii,ij)
          ibondi(ii,ij) = 0                         ! default: has e-w neighbours
          IF( ii   ==    1 )   ibondi(ii,ij) = -1   ! first column, has only e neighbour
@@ -391,8 +399,8 @@
          ioea(ii,ij) = iarea0 + 1
          iono(ii,ij) = iarea0 + jpni
-         ildi(ii,ij) =  1  + nn_hls
-         ilei(ii,ij) = ili - nn_hls
-         ildj(ii,ij) =  1  + nn_hls
-         ilej(ii,ij) = ilj - nn_hls
+         iis0(ii,ij) =  1  + nn_hls
+         iie0(ii,ij) = ili - nn_hls
+         ijs0(ii,ij) =  1  + nn_hls
+         ije0(ii,ij) = ilj - nn_hls
 
          ! East-West periodicity: change ibondi, ioea, iowe
@@ -432,34 +440,4 @@
       ! ----------------------------
       !
-      ! specify which subdomains are oce subdomains; other are land subdomains
-      ipproc(:,:) = -1
-      icont = -1
-      DO jarea = 1, jpni*jpnj
-         iarea0 = jarea - 1
-         ii = 1 + MOD(iarea0,jpni)
-         ij = 1 +     iarea0/jpni
-         IF( llisoce(ii,ij) ) THEN
-            icont = icont + 1
-            ipproc(ii,ij) = icont
-            iin(icont+1) = ii
-            ijn(icont+1) = ij
-         ENDIF
-      END DO
-      ! if needed add some land subdomains to reach jpnij active subdomains
-      i2add = jpnij - inijmin
-      DO jarea = 1, jpni*jpnj
-         iarea0 = jarea - 1
-         ii = 1 + MOD(iarea0,jpni)
-         ij = 1 +     iarea0/jpni
-         IF( .NOT. llisoce(ii,ij) .AND. i2add > 0 ) THEN
-            icont = icont + 1
-            ipproc(ii,ij) = icont
-            iin(icont+1) = ii
-            ijn(icont+1) = ij
-            i2add = i2add - 1
-         ENDIF
-      END DO
-      nfipproc(:,:) = ipproc(:,:)
-
       ! neighbour treatment: change ibondi, ibondj if next to a land zone
       DO jarea = 1, jpni*jpnj
@@ -500,15 +478,4 @@
          ENDIF
       END DO
-
-      ! Update il[de][ij] according to modified ibond[ij]
-      ! ----------------------
-      DO jproc = 1, jpnij
-         ii = iin(jproc)
-         ij = ijn(jproc)
-         IF( ibondi(ii,ij) == -1 .OR. ibondi(ii,ij) == 2 ) ildi(ii,ij) =  1
-         IF( ibondi(ii,ij) ==  1 .OR. ibondi(ii,ij) == 2 ) ilei(ii,ij) = ilci(ii,ij)
-         IF( ibondj(ii,ij) == -1 .OR. ibondj(ii,ij) == 2 ) ildj(ii,ij) =  1
-         IF( ibondj(ii,ij) ==  1 .OR. ibondj(ii,ij) == 2 ) ilej(ii,ij) = ilcj(ii,ij)
-      END DO
       
       ! 5. Subdomain print
@@ -523,5 +490,5 @@
             DO jj = jpnj, 1, -1
                WRITE(numout,9403) ('   ',ji=il1,il2-1)
-               WRITE(numout,9402) jj, (ilci(ji,jj),ilcj(ji,jj),ji=il1,il2)
+               WRITE(numout,9402) jj, (ijpi(ji,jj),ijpj(ji,jj),ji=il1,il2)
                WRITE(numout,9404) (ipproc(ji,jj),ji=il1,il2)
                WRITE(numout,9403) ('   ',ji=il1,il2-1)
@@ -580,36 +547,31 @@
       noea = ii_noea(narea)
       nono = ii_nono(narea)
-      nlci = ilci(ii,ij)  
-      nldi = ildi(ii,ij)
-      nlei = ilei(ii,ij)
-      nlcj = ilcj(ii,ij)  
-      nldj = ildj(ii,ij)
-      nlej = ilej(ii,ij)
+      jpi    = ijpi(ii,ij)  
+!!$      Nis0  = iis0(ii,ij)
+!!$      Nie0  = iie0(ii,ij)
+      jpj    = ijpj(ii,ij)  
+!!$      Njs0  = ijs0(ii,ij)
+!!$      Nje0  = ije0(ii,ij)
       nbondi = ibondi(ii,ij)
       nbondj = ibondj(ii,ij)
       nimpp = iimppt(ii,ij)  
       njmpp = ijmppt(ii,ij)
-      jpi = nlci
-      jpj = nlcj
-      jpk = jpkglo                                             ! third dim
-#if defined key_agrif
-      ! simple trick to use same vertical grid as parent but different number of levels: 
-      ! Save maximum number of levels in jpkglo, then define all vertical grids with this number.
-      ! Suppress once vertical online interpolation is ok
-!!$      IF(.NOT.Agrif_Root())   jpkglo = Agrif_Parent( jpkglo )
-#endif
-      jpim1 = jpi-1                                            ! inner domain indices
-      jpjm1 = jpj-1                                            !   "           "
-      jpkm1 = MAX( 1, jpk-1 )                                  !   "           "
-      jpij  = jpi*jpj                                          !  jpi x j
+      jpk = jpkglo                              ! third dim
+      !
+      CALL init_doloop                          ! set start/end indices of do-loop, depending on the halo width value (nn_hls) 
+      !
+      jpim1 = jpi-1                             ! inner domain indices
+      jpjm1 = jpj-1                             !   "           "
+      jpkm1 = MAX( 1, jpk-1 )                   !   "           "
+      jpij  = jpi*jpj                           !  jpi x j
       DO jproc = 1, jpnij
          ii = iin(jproc)
          ij = ijn(jproc)
-         nlcit(jproc) = ilci(ii,ij)
-         nldit(jproc) = ildi(ii,ij)
-         nleit(jproc) = ilei(ii,ij)
-         nlcjt(jproc) = ilcj(ii,ij)
-         nldjt(jproc) = ildj(ii,ij)
-         nlejt(jproc) = ilej(ii,ij)
+         jpiall (jproc) = ijpi(ii,ij)
+         nis0all(jproc) = iis0(ii,ij)
+         nie0all(jproc) = iie0(ii,ij)
+         jpjall (jproc) = ijpj(ii,ij)
+         njs0all(jproc) = ijs0(ii,ij)
+         nje0all(jproc) = ije0(ii,ij)
          ibonit(jproc) = ibondi(ii,ij)
          ibonjt(jproc) = ibondj(ii,ij)
@@ -625,10 +587,10 @@
          WRITE(inum,'(6i8,a,3i8,a)') jpnij,jpimax,jpjmax,jpk,jpiglo,jpjglo,&
    &           ' ( local: ',narea,jpi,jpj,' )'
-         WRITE(inum,'(a)') 'nproc nlci nlcj nldi nldj nlei nlej nimp njmp nono noso nowe noea nbondi nbondj '
+         WRITE(inum,'(a)') 'nproc   jpi  jpj Nis0 Njs0 Nie0 Nje0 nimp njmp nono noso nowe noea nbondi nbondj '
 
          DO jproc = 1, jpnij
-            WRITE(inum,'(13i5,2i7)')   jproc-1, nlcit  (jproc), nlcjt  (jproc),   &
-               &                                nldit  (jproc), nldjt  (jproc),   &
-               &                                nleit  (jproc), nlejt  (jproc),   &
+            WRITE(inum,'(13i5,2i7)')   jproc-1,  jpiall(jproc),  jpjall(jproc),   &
+               &                                nis0all(jproc), njs0all(jproc),   &
+               &                                nie0all(jproc), nje0all(jproc),   &
                &                                nimppt (jproc), njmppt (jproc),   & 
                &                                ii_nono(jproc), ii_noso(jproc),   &
@@ -664,11 +626,6 @@
          WRITE(numout,*) '    l_Iperio = ', l_Iperio
          WRITE(numout,*) '    l_Jperio = ', l_Jperio
-         WRITE(numout,*) '      nlci   = ', nlci
-         WRITE(numout,*) '      nlcj   = ', nlcj
          WRITE(numout,*) '      nimpp  = ', nimpp
          WRITE(numout,*) '      njmpp  = ', njmpp
-         WRITE(numout,*) '      nreci  = ', nreci  
-         WRITE(numout,*) '      nrecj  = ', nrecj  
-         WRITE(numout,*) '      nn_hls = ', nn_hls 
       ENDIF
 
@@ -692,14 +649,12 @@
       ENDIF
       !
-      CALL mpp_init_ioipsl       ! Prepare NetCDF output file (if necessary)
+      CALL init_ioipsl       ! Prepare NetCDF output file (if necessary)
       !      
       IF (( jperio >= 3 .AND. jperio <= 6 .AND. jpni > 1 ).AND.( ln_nnogather )) THEN
-         CALL mpp_init_nfdcom     ! northfold neighbour lists
+         CALL init_nfdcom     ! northfold neighbour lists
          IF (llwrtlay) THEN
             WRITE(inum,*)
             WRITE(inum,*)
             WRITE(inum,*) 'north fold exchanges with explicit point-to-point messaging :'
-            WRITE(inum,*) 'nfsloop : ', nfsloop
-            WRITE(inum,*) 'nfeloop : ', nfeloop
             WRITE(inum,*) 'nsndto : ', nsndto
             WRITE(inum,*) 'isendto : ', isendto
@@ -711,13 +666,14 @@
       DEALLOCATE(iin, ijn, ii_nono, ii_noea, ii_noso, ii_nowe,    &
          &       iimppt, ijmppt, ibondi, ibondj, ipproc, ipolj,   &
-         &       ilci, ilcj, ilei, ilej, ildi, ildj,              &
+         &       ijpi, ijpj, iie0, ije0, iis0, ijs0,              &
          &       iono, ioea, ioso, iowe, llisoce)
       !
     END SUBROUTINE mpp_init
 
-
-    SUBROUTINE mpp_basic_decomposition( knbi, knbj, kimax, kjmax, kimppt, kjmppt, klci, klcj)
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE mpp_basic_decomposition  ***
+#endif
+
+    SUBROUTINE mpp_basesplit( kiglo, kjglo, khls, knbi, knbj, kimax, kjmax, kimppt, kjmppt, klci, klcj)
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE mpp_basesplit  ***
       !!                    
       !! ** Purpose :   Lay out the global domain over processors.
@@ -731,4 +687,6 @@
       !!                    klcj       : second dimension
       !!----------------------------------------------------------------------
+      INTEGER,                                 INTENT(in   ) ::   kiglo, kjglo
+      INTEGER,                                 INTENT(in   ) ::   khls
       INTEGER,                                 INTENT(in   ) ::   knbi, knbj
       INTEGER,                                 INTENT(  out) ::   kimax, kjmax
@@ -737,14 +695,15 @@
       !
       INTEGER ::   ji, jj
+      INTEGER ::   i2hls 
       INTEGER ::   iresti, irestj, irm, ijpjmin
-      INTEGER ::   ireci, irecj
-      !!----------------------------------------------------------------------
+      !!----------------------------------------------------------------------
+      i2hls = 2*khls
       !
 #if defined key_nemocice_decomp
-      kimax = ( nx_global+2-2*nn_hls + (knbi-1) ) / knbi + 2*nn_hls    ! first  dim.
-      kjmax = ( ny_global+2-2*nn_hls + (knbj-1) ) / knbj + 2*nn_hls    ! second dim. 
+      kimax = ( nx_global+2-i2hls + (knbi-1) ) / knbi + i2hls    ! first  dim.
+      kjmax = ( ny_global+2-i2hls + (knbj-1) ) / knbj + i2hls    ! second dim. 
 #else
-      kimax = ( jpiglo - 2*nn_hls + (knbi-1) ) / knbi + 2*nn_hls    ! first  dim.
-      kjmax = ( jpjglo - 2*nn_hls + (knbj-1) ) / knbj + 2*nn_hls    ! second dim.
+      kimax = ( kiglo - i2hls + (knbi-1) ) / knbi + i2hls    ! first  dim.
+      kjmax = ( kjglo - i2hls + (knbj-1) ) / knbj + i2hls    ! second dim.
 #endif
       IF( .NOT. PRESENT(kimppt) ) RETURN
@@ -753,25 +712,23 @@
       ! -----------------------------------
       !  Computation of local domain sizes klci() klcj()
-      !  These dimensions depend on global sizes knbi,knbj and jpiglo,jpjglo
+      !  These dimensions depend on global sizes knbi,knbj and kiglo,kjglo
       !  The subdomains are squares lesser than or equal to the global
       !  dimensions divided by the number of processors minus the overlap array.
       !
-      ireci = 2 * nn_hls
-      irecj = 2 * nn_hls
-      iresti = 1 + MOD( jpiglo - ireci -1 , knbi )
-      irestj = 1 + MOD( jpjglo - irecj -1 , knbj )
+      iresti = 1 + MOD( kiglo - i2hls - 1 , knbi )
+      irestj = 1 + MOD( kjglo - i2hls - 1 , knbj )
       !
       !  Need to use kimax and kjmax here since jpi and jpj not yet defined
 #if defined key_nemocice_decomp
       ! Change padding to be consistent with CICE
-      klci(1:knbi-1      ,:) = kimax
-      klci(knbi          ,:) = jpiglo - (knbi - 1) * (kimax - nreci)
-      klcj(:,      1:knbj-1) = kjmax
-      klcj(:,          knbj) = jpjglo - (knbj - 1) * (kjmax - nrecj)
+      klci(1:knbi-1,:       ) = kimax
+      klci(  knbi  ,:       ) = kiglo - (knbi - 1) * (kimax - i2hls)
+      klcj(:       ,1:knbj-1) = kjmax
+      klcj(:       ,  knbj  ) = kjglo - (knbj - 1) * (kjmax - i2hls)
 #else
       klci(1:iresti      ,:) = kimax
       klci(iresti+1:knbi ,:) = kimax-1
-      IF( MINVAL(klci) < 3 ) THEN
-         WRITE(ctmp1,*) '   mpp_basic_decomposition: minimum value of jpi must be >= 3'
+      IF( MINVAL(klci) < 2*i2hls ) THEN
+         WRITE(ctmp1,*) '   mpp_basesplit: minimum value of jpi must be >= ', 2*i2hls
          WRITE(ctmp2,*) '   We have ', MINVAL(klci)
         CALL ctl_stop( 'STOP', ctmp1, ctmp2 )
@@ -779,19 +736,17 @@
       IF( jperio == 3 .OR. jperio == 4 .OR. jperio == 5 .OR. jperio == 6 ) THEN
          ! minimize the size of the last row to compensate for the north pole folding coast
-         IF( jperio == 3 .OR. jperio == 4 )   ijpjmin = 5   ! V and F folding involves line jpj-3 that must not be south boundary
-         IF( jperio == 5 .OR. jperio == 6 )   ijpjmin = 4   ! V and F folding involves line jpj-2 that must not be south boundary
-         irm = knbj - irestj                                    ! total number of lines to be removed
-         klcj(:,            knbj) = MAX( ijpjmin, kjmax-irm )   ! we must have jpj >= ijpjmin in the last row
-         irm = irm - ( kjmax - klcj(1,knbj) )                   ! remaining number of lines to remove 
-         irestj = knbj - 1 - irm                        
-         klcj(:,        1:irestj) = kjmax
+         IF( jperio == 3 .OR. jperio == 4 )   ijpjmin = 2+3*khls   ! V and F folding must be outside of southern halos
+         IF( jperio == 5 .OR. jperio == 6 )   ijpjmin = 1+3*khls   ! V and F folding must be outside of southern halos
+         irm = knbj - irestj                                       ! total number of lines to be removed
+         klcj(:,knbj) = MAX( ijpjmin, kjmax-irm )                  ! we must have jpj >= ijpjmin in the last row
+         irm = irm - ( kjmax - klcj(1,knbj) )                      ! remaining number of lines to remove 
+         irestj = knbj - 1 - irm
          klcj(:, irestj+1:knbj-1) = kjmax-1
       ELSE
-         ijpjmin = 3
-         klcj(:,      1:irestj) = kjmax
-         klcj(:, irestj+1:knbj) = kjmax-1
-      ENDIF
-      IF( MINVAL(klcj) < ijpjmin ) THEN
-         WRITE(ctmp1,*) '   mpp_basic_decomposition: minimum value of jpj must be >= ', ijpjmin
+         klcj(:, irestj+1:knbj  ) = kjmax-1
+      ENDIF
+      klcj(:,1:irestj) = kjmax
+      IF( MINVAL(klcj) < 2*i2hls ) THEN
+         WRITE(ctmp1,*) '   mpp_basesplit: minimum value of jpj must be >= ', 2*i2hls
          WRITE(ctmp2,*) '   We have ', MINVAL(klcj)
          CALL ctl_stop( 'STOP', ctmp1, ctmp2 )
@@ -807,5 +762,5 @@
          DO jj = 1, knbj
             DO ji = 2, knbi
-               kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - ireci
+               kimppt(ji,jj) = kimppt(ji-1,jj) + klci(ji-1,jj) - i2hls
             END DO
          END DO
@@ -815,15 +770,15 @@
          DO jj = 2, knbj
             DO ji = 1, knbi
-               kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - irecj
+               kjmppt(ji,jj) = kjmppt(ji,jj-1) + klcj(ji,jj-1) - i2hls
             END DO
          END DO
       ENDIF
       
-   END SUBROUTINE mpp_basic_decomposition
-
-
-   SUBROUTINE mpp_init_bestpartition( knbij, knbi, knbj, knbcnt, ldlist )
-      !!----------------------------------------------------------------------
-      !!                 ***  ROUTINE mpp_init_bestpartition  ***
+   END SUBROUTINE mpp_basesplit
+
+
+   SUBROUTINE bestpartition( knbij, knbi, knbj, knbcnt, ldlist )
+      !!----------------------------------------------------------------------
+      !!                 ***  ROUTINE bestpartition  ***
       !!
       !! ** Purpose :
@@ -831,5 +786,5 @@
       !! ** Method  :
       !!----------------------------------------------------------------------
-      INTEGER,           INTENT(in   ) ::   knbij         ! total number if subdomains               (knbi*knbj)
+      INTEGER,           INTENT(in   ) ::   knbij         ! total number of subdomains (knbi*knbj)
       INTEGER, OPTIONAL, INTENT(  out) ::   knbi, knbj    ! number if subdomains along i and j (knbi and knbj)
       INTEGER, OPTIONAL, INTENT(  out) ::   knbcnt        ! number of land subdomains
@@ -839,4 +794,5 @@
       INTEGER :: iszitst, iszjtst
       INTEGER :: isziref, iszjref
+      INTEGER :: iszimin, iszjmin
       INTEGER :: inbij, iszij
       INTEGER :: inbimax, inbjmax, inbijmax, inbijold
@@ -867,6 +823,11 @@
       inbimax = 0
       inbjmax = 0
-      isziref = jpiglo*jpjglo+1
+      isziref = jpiglo*jpjglo+1   ! define a value that is larger than the largest possible
       iszjref = jpiglo*jpjglo+1
+      !
+      iszimin = 4*nn_hls          ! minimum size of the MPI subdomain so halos are always adressing neighbor inner domain
+      iszjmin = 4*nn_hls
+      IF( jperio == 3 .OR. jperio == 4 )   iszjmin = MAX(iszjmin, 2+3*nn_hls)   ! V and F folding must be outside of southern halos
+      IF( jperio == 5 .OR. jperio == 6 )   iszjmin = MAX(iszjmin, 1+3*nn_hls)   ! V and F folding must be outside of southern halos
       !
       ! get the list of knbi that gives a smaller jpimax than knbi-1
@@ -876,7 +837,7 @@
          iszitst = ( nx_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim.
 #else
-         iszitst = ( jpiglo - 2*nn_hls + (ji-1) ) / ji + 2*nn_hls
+         iszitst = ( Ni0glo + (ji-1) ) / ji + 2*nn_hls   ! max subdomain i-size
 #endif
-         IF( iszitst < isziref ) THEN
+         IF( iszitst < isziref .AND. iszitst >= iszimin ) THEN
             isziref = iszitst
             inbimax = inbimax + 1
@@ -887,7 +848,7 @@
          iszjtst = ( ny_global+2-2*nn_hls + (ji-1) ) / ji + 2*nn_hls    ! first  dim.
 #else
-         iszjtst = ( jpjglo - 2*nn_hls + (ji-1) ) / ji + 2*nn_hls
+         iszjtst = ( Nj0glo + (ji-1) ) / ji + 2*nn_hls   ! max subdomain j-size
 #endif
-         IF( iszjtst < iszjref ) THEN
+         IF( iszjtst < iszjref .AND. iszjtst >= iszjmin ) THEN
             iszjref = iszjtst
             inbjmax = inbjmax + 1
@@ -927,5 +888,5 @@
       iszij1(:) = iszi1(:) * iszj1(:)
 
-      ! if therr is no land and no print
+      ! if there is no land and no print
       IF( .NOT. llist .AND. numbot == -1 .AND. numbdy == -1 ) THEN
          ! get the smaller partition which gives the smallest subdomain size
@@ -946,4 +907,5 @@
          ii = MINLOC(iszij1, mask = inbij1 == inbij, dim = 1)   ! warning: send back the first occurence if multiple results
          IF ( iszij1(ii) < iszij ) THEN
+            ii = MINLOC( iszi1+iszj1, mask = iszij1 == iszij1(ii) .AND. inbij1 == inbij, dim = 1)  ! select the smaller perimeter if multiple min
             isz0 = isz0 + 1
             indexok(isz0) = ii
@@ -975,10 +937,10 @@
          ji = isz0   ! initialization with the largest value
          ALLOCATE( llisoce(inbi0(ji), inbj0(ji)) )
-         CALL mpp_init_isoce( inbi0(ji), inbj0(ji), llisoce ) ! Warning: must be call by all cores (call mpp_sum)
+         CALL mpp_is_ocean( llisoce )   ! Warning: must be call by all cores (call mpp_sum)
          inbijold = COUNT(llisoce)
          DEALLOCATE( llisoce )
          DO ji =isz0-1,1,-1
             ALLOCATE( llisoce(inbi0(ji), inbj0(ji)) )
-            CALL mpp_init_isoce( inbi0(ji), inbj0(ji), llisoce ) ! Warning: must be call by all cores (call mpp_sum)
+            CALL mpp_is_ocean( llisoce )   ! Warning: must be call by all cores (call mpp_sum)
             inbij = COUNT(llisoce)
             DEALLOCATE( llisoce )
@@ -1006,5 +968,5 @@
          ii = ii -1 
          ALLOCATE( llisoce(inbi0(ii), inbj0(ii)) )
-         CALL mpp_init_isoce( inbi0(ii), inbj0(ii), llisoce )            ! must be done by all core
+         CALL mpp_is_ocean( llisoce )            ! must be done by all core
          inbij = COUNT(llisoce)
          DEALLOCATE( llisoce )
@@ -1015,5 +977,5 @@
       DEALLOCATE( inbi0, inbj0 )
       !
-   END SUBROUTINE mpp_init_bestpartition
+   END SUBROUTINE bestpartition
    
    
@@ -1024,5 +986,5 @@
       !! ** Purpose : the the proportion of land points in the surface land-sea mask
       !!
-      !! ** Method  : read iproc strips (of length jpiglo) of the land-sea mask
+      !! ** Method  : read iproc strips (of length Ni0glo) of the land-sea mask
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(  out) :: propland    ! proportion of land points in the global domain (between 0 and 1)
@@ -1041,5 +1003,5 @@
 
       ! number of processes reading the bathymetry file 
-      iproc = MINVAL( (/mppsize, jpjglo/2, 100/) )  ! read a least 2 lines, no more that 100 processes reading at the same time
+      iproc = MINVAL( (/mppsize, Nj0glo/2, 100/) )  ! read a least 2 lines, no more that 100 processes reading at the same time
       
       ! we want to read iproc strips of the land-sea mask. -> pick up iproc processes every idiv processes starting at 1
@@ -1051,10 +1013,10 @@
       IF( MOD( narea-1, idiv ) == 0 .AND. iarea < iproc ) THEN   ! beware idiv can be = to 1
          !
-         ijsz = jpjglo / iproc                                               ! width of the stripe to read
-         IF( iarea < MOD(jpjglo,iproc) ) ijsz = ijsz + 1
-         ijstr = iarea*(jpjglo/iproc) + MIN(iarea, MOD(jpjglo,iproc)) + 1    ! starting j position of the reading
+         ijsz = Nj0glo / iproc                                               ! width of the stripe to read
+         IF( iarea < MOD(Nj0glo,iproc) ) ijsz = ijsz + 1
+         ijstr = iarea*(Nj0glo/iproc) + MIN(iarea, MOD(Nj0glo,iproc)) + 1    ! starting j position of the reading
          !
-         ALLOCATE( lloce(jpiglo, ijsz) )                                     ! allocate the strip
-         CALL mpp_init_readbot_strip( ijstr, ijsz, lloce )
+         ALLOCATE( lloce(Ni0glo, ijsz) )                                     ! allocate the strip
+         CALL readbot_strip( ijstr, ijsz, lloce )
          inboce = COUNT(lloce)                                               ! number of ocean point in the stripe
          DEALLOCATE(lloce)
@@ -1065,28 +1027,34 @@
       CALL mpp_sum( 'mppini', inboce )   ! total number of ocean points over the global domain
       !
-      propland = REAL( jpiglo*jpjglo - inboce, wp ) / REAL( jpiglo*jpjglo, wp ) 
+      propland = REAL( Ni0glo*Nj0glo - inboce, wp ) / REAL( Ni0glo*Nj0glo, wp ) 
       !
    END SUBROUTINE mpp_init_landprop
    
    
-   SUBROUTINE mpp_init_isoce( knbi, knbj, ldisoce )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE mpp_init_nboce  ***
-      !!
-      !! ** Purpose : check for a mpi domain decomposition knbi x knbj which
-      !!              subdomains contain at least 1 ocean point
-      !!
-      !! ** Method  : read knbj strips (of length jpiglo) of the land-sea mask
-      !!----------------------------------------------------------------------
-      INTEGER,                       INTENT(in   ) ::   knbi, knbj     ! domain decomposition
-      LOGICAL, DIMENSION(knbi,knbj), INTENT(  out) ::   ldisoce        ! .true. if a sub domain constains 1 ocean point 
-      !
-      INTEGER, DIMENSION(knbi,knbj) ::   inboce                        ! number oce oce pint in each mpi subdomain
-      INTEGER, DIMENSION(knbi*knbj) ::   inboce_1d
+   SUBROUTINE mpp_is_ocean( ldisoce )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE mpp_is_ocean  ***
+      !!
+      !! ** Purpose : Check for a mpi domain decomposition inbi x inbj which
+      !!              subdomains, including 1 halo (even if nn_hls>1), contain
+      !!              at least 1 ocean point.
+      !!              We must indeed ensure that each subdomain that is a neighbour
+      !!              of a land subdomain as only land points on its boundary
+      !!              (inside the inner subdomain) with the land subdomain.
+      !!              This is needed to get the proper bondary conditions on
+      !!              a subdomain with a closed boundary.
+      !!
+      !! ** Method  : read inbj strips (of length Ni0glo) of the land-sea mask
+      !!----------------------------------------------------------------------
+      LOGICAL, DIMENSION(:,:), INTENT(  out) ::   ldisoce        ! .true. if a sub domain constains 1 ocean point 
+      !
       INTEGER :: idiv, iimax, ijmax, iarea
+      INTEGER :: inbi, inbj, inx, iny, inry, isty
       INTEGER :: ji, jn
-      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce                  ! lloce(i,j) = .true. if the point (i,j) is ocean 
-      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ilci
-      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ilcj
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   inboce           ! number oce oce pint in each mpi subdomain
+      INTEGER, ALLOCATABLE, DIMENSION(:  ) ::   inboce_1d
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   iimppt, ijpi
+      INTEGER, ALLOCATABLE, DIMENSION(:,:) ::   ijmppt, ijpj
+      LOGICAL, ALLOCATABLE, DIMENSION(:,:) ::   lloce            ! lloce(i,j) = .true. if the point (i,j) is ocean 
       !!----------------------------------------------------------------------
       ! do nothing if there is no land-sea mask
@@ -1095,43 +1063,88 @@
          RETURN
       ENDIF
-
-      ! we want to read knbj strips of the land-sea mask. -> pick up knbj processes every idiv processes starting at 1
-      IF           ( knbj == 1 ) THEN   ;   idiv = mppsize
-      ELSE IF ( mppsize < knbj ) THEN   ;   idiv = 1
-      ELSE                              ;   idiv = ( mppsize - 1 ) / ( knbj - 1 )
-      ENDIF
+      !
+      inbi = SIZE( ldisoce, dim = 1 )
+      inbj = SIZE( ldisoce, dim = 2 )
+      !
+      ! we want to read inbj strips of the land-sea mask. -> pick up inbj processes every idiv processes starting at 1
+      IF           ( inbj == 1 ) THEN   ;   idiv = mppsize
+      ELSE IF ( mppsize < inbj ) THEN   ;   idiv = 1
+      ELSE                              ;   idiv = ( mppsize - 1 ) / ( inbj - 1 )
+      ENDIF
+      !
+      ALLOCATE( inboce(inbi,inbj), inboce_1d(inbi*inbj) )
       inboce(:,:) = 0          ! default no ocean point found
-
-      DO jn = 0, (knbj-1)/mppsize   ! if mppsize < knbj : more strips than mpi processes (because of potential land domains)
+      !
+      DO jn = 0, (inbj-1)/mppsize   ! if mppsize < inbj : more strips than mpi processes (because of potential land domains)
          !
-         iarea = (narea-1)/idiv + jn * mppsize   ! involed process number (starting counting at 0)
-         IF( MOD( narea-1, idiv ) == 0 .AND. iarea < knbj ) THEN   ! beware idiv can be = to 1
+         iarea = (narea-1)/idiv + jn * mppsize + 1                     ! involed process number (starting counting at 1)
+         IF( MOD( narea-1, idiv ) == 0 .AND. iarea <= inbj ) THEN      ! beware idiv can be = to 1
             !
-            ALLOCATE( iimppt(knbi,knbj), ijmppt(knbi,knbj), ilci(knbi,knbj), ilcj(knbi,knbj) )
-            CALL mpp_basic_decomposition( knbi, knbj, iimax, ijmax, iimppt, ijmppt, ilci, ilcj )
+            ALLOCATE( iimppt(inbi,inbj), ijmppt(inbi,inbj), ijpi(inbi,inbj), ijpj(inbi,inbj) )
+            CALL mpp_basesplit( Ni0glo, Nj0glo, 0, inbi, inbj, iimax, ijmax, iimppt, ijmppt, ijpi, ijpj )
             !
-            ALLOCATE( lloce(jpiglo, ilcj(1,iarea+1)) )                                         ! allocate the strip
-            CALL mpp_init_readbot_strip( ijmppt(1,iarea+1), ilcj(1,iarea+1), lloce )           ! read the strip
-            DO  ji = 1, knbi
-               inboce(ji,iarea+1) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ilci(ji,1)-1,:) )   ! number of ocean point in subdomain
+            inx = Ni0glo + 2   ;   iny = ijpj(1,iarea) + 2             ! strip size + 1 halo on each direction (even if nn_hls>1)
+            ALLOCATE( lloce(inx, iny) )                                ! allocate the strip
+            inry = iny - COUNT( (/ iarea == 1, iarea == inbj /) )      ! number of point to read in y-direction
+            isty = 1 + COUNT( (/ iarea == 1 /) )                       ! read from the first or the second line?
+            CALL readbot_strip( ijmppt(1,iarea) - 2 + isty, inry, lloce(2:inx-1, isty:inry+isty-1) )   ! read the strip
+            ! 
+            IF( iarea == 1    ) THEN                                   ! the first line was not read
+               IF( jperio == 2 .OR. jperio == 7 ) THEN                 !   north-south periodocity
+                  CALL readbot_strip( Nj0glo, 1, lloce(2:inx-1, 1) )   !   read the last line -> first line of lloce
+               ELSE
+                  lloce(2:inx-1,  1) = .FALSE.                         !   closed boundary
+               ENDIF
+            ENDIF
+            IF( iarea == inbj ) THEN                                   ! the last line was not read
+               IF( jperio == 2 .OR. jperio == 7 ) THEN                 !   north-south periodocity
+                  CALL readbot_strip( 1, 1, lloce(2:inx-1,iny) )       !      read the first line -> last line of lloce
+               ELSEIF( jperio == 3 .OR. jperio == 4 ) THEN             !   north-pole folding T-pivot, T-point 
+                  lloce(2,iny) = lloce(2,iny-2)                        !      here we have 1 halo (even if nn_hls>1)
+                  DO ji = 3,inx-1
+                     lloce(ji,iny  ) = lloce(inx-ji+2,iny-2)           !      ok, we have at least 3 lines
+                  END DO
+                  DO ji = inx/2+2,inx-1
+                     lloce(ji,iny-1) = lloce(inx-ji+2,iny-1)
+                  END DO
+               ELSEIF( jperio == 5 .OR. jperio == 6 ) THEN             !   north-pole folding F-pivot, T-point, 1 halo
+                  lloce(inx/2+1,iny-1) = lloce(inx/2,iny-1)            !      here we have 1 halo (even if nn_hls>1)
+                  lloce(inx  -1,iny-1) = lloce(2    ,iny-1)
+                  DO ji = 2,inx-1
+                     lloce(ji,iny) = lloce(inx-ji+1,iny-1)
+                  END DO
+               ELSE                                                    !   closed boundary
+                  lloce(2:inx-1,iny) = .FALSE.
+               ENDIF
+            ENDIF
+            !                                                          ! first and last column were not read
+            IF( jperio == 1 .OR. jperio == 4 .OR. jperio == 6 .OR. jperio == 7 ) THEN
+               lloce(1,:) = lloce(inx-1,:)   ;   lloce(inx,:) = lloce(2,:)   ! east-west periodocity
+            ELSE
+               lloce(1,:) = .FALSE.          ;   lloce(inx,:) = .FALSE.      ! closed boundary
+            ENDIF
+            !
+            DO  ji = 1, inbi
+               inboce(ji,iarea) = COUNT( lloce(iimppt(ji,1):iimppt(ji,1)+ijpi(ji,1)+1,:) )   ! lloce as 2 points more than Ni0glo
             END DO
             !
             DEALLOCATE(lloce)
-            DEALLOCATE(iimppt, ijmppt, ilci, ilcj)
+            DEALLOCATE(iimppt, ijmppt, ijpi, ijpj)
             !
          ENDIF
       END DO
    
-      inboce_1d = RESHAPE(inboce, (/ knbi*knbj /))
+      inboce_1d = RESHAPE(inboce, (/ inbi*inbj /))
       CALL mpp_sum( 'mppini', inboce_1d )
-      inboce = RESHAPE(inboce_1d, (/knbi, knbj/))
+      inboce = RESHAPE(inboce_1d, (/inbi, inbj/))
       ldisoce(:,:) = inboce(:,:) /= 0
-      !
-   END SUBROUTINE mpp_init_isoce
+      DEALLOCATE(inboce, inboce_1d)
+      !
+   END SUBROUTINE mpp_is_ocean
    
    
-   SUBROUTINE mpp_init_readbot_strip( kjstr, kjcnt, ldoce )
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE mpp_init_readbot_strip  ***
+   SUBROUTINE readbot_strip( kjstr, kjcnt, ldoce )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE readbot_strip  ***
       !!
       !! ** Purpose : Read relevant bathymetric information in order to
@@ -1139,36 +1152,89 @@
       !!              of land domains, in an mpp computation.
       !!
-      !! ** Method  : read stipe of size (jpiglo,...)
-      !!----------------------------------------------------------------------
-      INTEGER                         , INTENT(in   ) :: kjstr       ! starting j position of the reading
-      INTEGER                         , INTENT(in   ) :: kjcnt       ! number of lines to read
-      LOGICAL, DIMENSION(jpiglo,kjcnt), INTENT(  out) :: ldoce       ! ldoce(i,j) = .true. if the point (i,j) is ocean 
+      !! ** Method  : read stipe of size (Ni0glo,...)
+      !!----------------------------------------------------------------------
+      INTEGER                         , INTENT(in   ) ::   kjstr       ! starting j position of the reading
+      INTEGER                         , INTENT(in   ) ::   kjcnt       ! number of lines to read
+      LOGICAL, DIMENSION(Ni0glo,kjcnt), INTENT(  out) ::   ldoce       ! ldoce(i,j) = .true. if the point (i,j) is ocean 
       !
       INTEGER                           ::   inumsave                ! local logical unit
-      REAL(wp), DIMENSION(jpiglo,kjcnt) ::   zbot, zbdy 
+      REAL(wp), DIMENSION(Ni0glo,kjcnt) ::   zbot, zbdy 
       !!----------------------------------------------------------------------
       !
       inumsave = numout   ;   numout = numnul   !   redirect all print to /dev/null
       !
-      IF( numbot /= -1 ) THEN
-         CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
+      IF( numbot /= -1 ) THEN   
+         CALL iom_get( numbot, jpdom_unknown, 'bottom_level', zbot, kstart = (/1,kjstr/), kcount = (/Ni0glo, kjcnt/) )
       ELSE
-         zbot(:,:) = 1.                         ! put a non-null value
-      ENDIF
-
-       IF( numbdy /= -1 ) THEN                  ! Adjust with bdy_msk if it exists    
-         CALL iom_get ( numbdy, jpdom_unknown, 'bdy_msk', zbdy, kstart = (/1,kjstr/), kcount = (/jpiglo, kjcnt/) )
+         zbot(:,:) = 1._wp                      ! put a non-null value
+      ENDIF
+      !
+      IF( numbdy /= -1 ) THEN                   ! Adjust with bdy_msk if it exists    
+         CALL iom_get ( numbdy, jpdom_unknown,     'bdy_msk', zbdy, kstart = (/1,kjstr/), kcount = (/Ni0glo, kjcnt/) )
          zbot(:,:) = zbot(:,:) * zbdy(:,:)
       ENDIF
       !
-      ldoce(:,:) = zbot(:,:) > 0.
+      ldoce(:,:) = zbot(:,:) > 0._wp
       numout = inumsave
       !
-   END SUBROUTINE mpp_init_readbot_strip
-
-
-   SUBROUTINE mpp_init_ioipsl
-      !!----------------------------------------------------------------------
-      !!                  ***  ROUTINE mpp_init_ioipsl  ***
+   END SUBROUTINE readbot_strip
+
+
+   SUBROUTINE mpp_getnum( ldisoce, kproc, kipos, kjpos )
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE mpp_getnum  ***
+      !!
+      !! ** Purpose : give a number to each MPI subdomains (starting at 0)
+      !!
+      !! ** Method  : start from bottom left. First skip land subdomain, and finally use them if needed
+      !!----------------------------------------------------------------------
+      LOGICAL, DIMENSION(:,:), INTENT(in   ) ::   ldisoce     ! F if land process
+      INTEGER, DIMENSION(:,:), INTENT(  out) ::   kproc       ! subdomain number (-1 if supressed, starting at 0)
+      INTEGER, DIMENSION(  :), INTENT(  out) ::   kipos       ! i-position of the subdomain (from 1 to jpni)
+      INTEGER, DIMENSION(  :), INTENT(  out) ::   kjpos       ! j-position of the subdomain (from 1 to jpnj)
+      !
+      INTEGER :: ii, ij, jarea, iarea0
+      INTEGER :: icont, i2add , ini, inj, inij
+      !!----------------------------------------------------------------------
+      !
+      ini = SIZE(ldisoce, dim = 1)
+      inj = SIZE(ldisoce, dim = 2)
+      inij = SIZE(kipos)
+      !
+      ! specify which subdomains are oce subdomains; other are land subdomains
+      kproc(:,:) = -1
+      icont = -1
+      DO jarea = 1, ini*inj
+         iarea0 = jarea - 1
+         ii = 1 + MOD(iarea0,ini)
+         ij = 1 +     iarea0/ini
+         IF( ldisoce(ii,ij) ) THEN
+            icont = icont + 1
+            kproc(ii,ij) = icont
+            kipos(icont+1) = ii
+            kjpos(icont+1) = ij
+         ENDIF
+      END DO
+      ! if needed add some land subdomains to reach inij active subdomains
+      i2add = inij - COUNT( ldisoce )
+      DO jarea = 1, ini*inj
+         iarea0 = jarea - 1
+         ii = 1 + MOD(iarea0,ini)
+         ij = 1 +     iarea0/ini
+         IF( .NOT. ldisoce(ii,ij) .AND. i2add > 0 ) THEN
+            icont = icont + 1
+            kproc(ii,ij) = icont
+            kipos(icont+1) = ii
+            kjpos(icont+1) = ij
+            i2add = i2add - 1
+         ENDIF
+      END DO
+      !
+   END SUBROUTINE mpp_getnum
+
+
+   SUBROUTINE init_ioipsl
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE init_ioipsl  ***
       !!
       !! ** Purpose :   
@@ -1187,34 +1253,28 @@
       ! Set idompar values equivalent to the jpdom_local_noextra definition
       ! used in IOM. This works even if jpnij .ne. jpni*jpnj.
-      iglo(1) = jpiglo
-      iglo(2) = jpjglo
-      iloc(1) = nlci
-      iloc(2) = nlcj
-      iabsf(1) = nimppt(narea)
-      iabsf(2) = njmppt(narea)
+      iglo( :) = (/ Ni0glo, Nj0glo /)
+      iloc( :) = (/ Ni_0  , Nj_0   /)
+      iabsf(:) = (/ Nis0  , Njs0   /) + (/ nimpp, njmpp /) - 1 - nn_hls   ! corresponds to mig0(Nis0) but mig0 is not yet defined!
       iabsl(:) = iabsf(:) + iloc(:) - 1
-      ihals(1) = nldi - 1
-      ihals(2) = nldj - 1
-      ihale(1) = nlci - nlei
-      ihale(2) = nlcj - nlej
-      idid(1) = 1
-      idid(2) = 2
+      ihals(:) = (/ 0     , 0      /)
+      ihale(:) = (/ 0     , 0      /)
+      idid( :) = (/ 1     , 2      /)
 
       IF(lwp) THEN
           WRITE(numout,*)
-          WRITE(numout,*) 'mpp_init_ioipsl :   iloc  = ', iloc (1), iloc (2)
-          WRITE(numout,*) '~~~~~~~~~~~~~~~     iabsf = ', iabsf(1), iabsf(2)
-          WRITE(numout,*) '                    ihals = ', ihals(1), ihals(2)
-          WRITE(numout,*) '                    ihale = ', ihale(1), ihale(2)
+          WRITE(numout,*) 'mpp init_ioipsl :   iloc  = ', iloc
+          WRITE(numout,*) '~~~~~~~~~~~~~~~     iabsf = ', iabsf
+          WRITE(numout,*) '                    ihals = ', ihals
+          WRITE(numout,*) '                    ihale = ', ihale
       ENDIF
       !
       CALL flio_dom_set ( jpnij, nproc, idid, iglo, iloc, iabsf, iabsl, ihals, ihale, 'BOX', nidom)
       !
-   END SUBROUTINE mpp_init_ioipsl  
-
-
-   SUBROUTINE mpp_init_nfdcom
-      !!----------------------------------------------------------------------
-      !!                     ***  ROUTINE  mpp_init_nfdcom  ***
+   END SUBROUTINE init_ioipsl  
+
+
+   SUBROUTINE init_nfdcom
+      !!----------------------------------------------------------------------
+      !!                     ***  ROUTINE  init_nfdcom  ***
       !! ** Purpose :   Setup for north fold exchanges with explicit 
       !!                point-to-point messaging
@@ -1226,8 +1286,5 @@
       !!----------------------------------------------------------------------
       INTEGER  ::   sxM, dxM, sxT, dxT, jn
-      INTEGER  ::   njmppmax
-      !!----------------------------------------------------------------------
-      !
-      njmppmax = MAXVAL( njmppt )
+      !!----------------------------------------------------------------------
       !
       !initializes the north-fold communication variables
@@ -1235,8 +1292,8 @@
       nsndto     = 0
       !
-      IF ( njmpp == njmppmax ) THEN      ! if I am a process in the north
+      IF ( njmpp == MAXVAL( njmppt ) ) THEN      ! if I am a process in the north
          !
          !sxM is the first point (in the global domain) needed to compute the north-fold for the current process
-         sxM = jpiglo - nimppt(narea) - nlcit(narea) + 1
+         sxM = jpiglo - nimppt(narea) - jpiall(narea) + 1
          !dxM is the last point (in the global domain) needed to compute the north-fold for the current process
          dxM = jpiglo - nimppt(narea) + 2
@@ -1247,6 +1304,6 @@
          DO jn = 1, jpni
             !
-            sxT = nfiimpp(jn, jpnj)                            ! sxT = 1st  point (in the global domain) of the jn process
-            dxT = nfiimpp(jn, jpnj) + nfilcit(jn, jpnj) - 1    ! dxT = last point (in the global domain) of the jn process
+            sxT = nfimpp(jn)                    ! sxT = 1st  point (in the global domain) of the jn process
+            dxT = nfimpp(jn) + nfjpi(jn) - 1    ! dxT = last point (in the global domain) of the jn process
             !
             IF    ( sxT < sxM  .AND.  sxM < dxT ) THEN
@@ -1262,21 +1319,46 @@
             !
          END DO
-         nfsloop = 1
-         nfeloop = nlci
-         DO jn = 2,jpni-1
-            IF( nfipproc(jn,jpnj) == (narea - 1) ) THEN
-               IF( nfipproc(jn-1,jpnj) == -1 )   nfsloop = nldi
-               IF( nfipproc(jn+1,jpnj) == -1 )   nfeloop = nlei
-            ENDIF
-         END DO
          !
       ENDIF
       l_north_nogather = .TRUE.
       !
-   END SUBROUTINE mpp_init_nfdcom
-
-
-#endif
-
+   END SUBROUTINE init_nfdcom
+
+
+   SUBROUTINE init_doloop
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE init_doloop  ***
+      !!
+      !! ** Purpose :   set the starting/ending indices of DO-loop
+      !!              These indices are used in do_loop_substitute.h90
+      !!----------------------------------------------------------------------
+      !
+      Nis0 =   1+nn_hls   ;   Nis1 = Nis0-1   ;   Nis2 = MAX(  1, Nis0-2)
+      Njs0 =   1+nn_hls   ;   Njs1 = Njs0-1   ;   Njs2 = MAX(  1, Njs0-2)  
+      !                                                 
+      Nie0 = jpi-nn_hls   ;   Nie1 = Nie0+1   ;   Nie2 = MIN(jpi, Nie0+2)
+      Nje0 = jpj-nn_hls   ;   Nje1 = Nje0+1   ;   Nje2 = MIN(jpj, Nje0+2)
+      !
+      IF( nn_hls == 1 ) THEN          !* halo size of 1
+         !
+         Nis1nxt2 = Nis0   ;   Njs1nxt2 = Njs0
+         Nie1nxt2 = Nie0   ;   Nje1nxt2 = Nje0
+         !
+      ELSE                            !* larger halo size... 
+         !
+         Nis1nxt2 = Nis1   ;   Njs1nxt2 = Njs1
+         Nie1nxt2 = Nie1   ;   Nje1nxt2 = Nje1
+         !
+      ENDIF
+      !
+      Ni_0 = Nie0 - Nis0 + 1
+      Nj_0 = Nje0 - Njs0 + 1
+      Ni_1 = Nie1 - Nis1 + 1
+      Nj_1 = Nje1 - Njs1 + 1
+      Ni_2 = Nie2 - Nis2 + 1
+      Nj_2 = Nje2 - Njs2 + 1
+      !
+   END SUBROUTINE init_doloop
+   
    !!======================================================================
 END MODULE mppini
