Index: branches/2017/dev_r7832_HPC08_lbclnk_3rd_dim/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
===================================================================
--- branches/2017/dev_r7832_HPC08_lbclnk_3rd_dim/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 7844)
+++ branches/2017/dev_r7832_HPC08_lbclnk_3rd_dim/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 7897)
@@ -8,5 +8,5 @@
    !!            8.0  !  1998  (M. Imbard, J. Escobar, L. Colombet ) SHMEM and MPI
    !!                 !  1998  (J.M. Molines) Open boundary conditions
-   !!   NEMO     1.0  !  2003  (J.-M. Molines, G. Madec)  F90, free form
+   !!   NEMO     1.0  !  2003  (J.M. Molines, G. Madec)  F90, free form
    !!                 !  2003  (J.M. Molines) add mpp_ini_north(_3d,_2d)
    !!             -   !  2004  (R. Bourdalle Badie)  isend option in mpi
@@ -24,5 +24,6 @@
    !!            3.5  !  2013  ( C. Ethe, G. Madec ) message passing arrays as local variables 
    !!            3.5  !  2013  (S.Mocavero, I.Epicoco - CMCC) north fold optimizations
-   !!            3.6  !  2015 (O. Tintó and M. Castrillo - BSC) Added 'mpp_lnk_2d_multiple', 'mpp_lbc_north_2d_multiple', 'mpp_max_multiple' 
+   !!            3.6  !  2015  (O. Tintó and M. Castrillo - BSC) Added '_multiple' case for 2D lbc and max
+   !!            4.0  !  2017  (G. Madec) automatique allocation of array argument (use any 3rd dimension)
    !!----------------------------------------------------------------------
 
@@ -45,5 +46,5 @@
    !!   mpp_lnk_icb   : interface for message passing of 2d arrays with extra halo for icebergs (mpp_lnk_2d_icb)
    !!   mpprecv       :
-   !!   mppsend       :   SUBROUTINE mpp_ini_znl
+   !!   mppsend       :
    !!   mppscatter    :
    !!   mppgather     :
@@ -85,6 +86,7 @@
 
    TYPE arrayptr
-      REAL , DIMENSION (:,:),  POINTER :: pt2d
+      REAL(wp), DIMENSION (:,:),  POINTER ::   pt2d
    END TYPE arrayptr
+   !
    PUBLIC   arrayptr
    
@@ -101,5 +103,5 @@
    INTERFACE mpp_sum
       MODULE PROCEDURE mppsum_a_int, mppsum_int, mppsum_a_real, mppsum_real,   &
-                       mppsum_realdd, mppsum_a_realdd
+         &             mppsum_realdd, mppsum_a_realdd
    END INTERFACE
    INTERFACE mpp_lbc_north
@@ -112,5 +114,4 @@
       MODULE PROCEDURE mpp_maxloc2d ,mpp_maxloc3d
    END INTERFACE
-
    INTERFACE mpp_max_multiple
       MODULE PROCEDURE mppmax_real_multiple
@@ -138,39 +139,39 @@
    ! variables used in case of sea-ice
    INTEGER, PUBLIC ::   ncomm_ice       !: communicator made by the processors with sea-ice (public so that it can be freed in limthd)
-   INTEGER ::   ngrp_iworld     !  group ID for the world processors (for rheology)
-   INTEGER ::   ngrp_ice        !  group ID for the ice processors (for rheology)
-   INTEGER ::   ndim_rank_ice   !  number of 'ice' processors
-   INTEGER ::   n_ice_root      !  number (in the comm_ice) of proc 0 in the ice comm
+   INTEGER         ::   ngrp_iworld     !  group ID for the world processors (for rheology)
+   INTEGER         ::   ngrp_ice        !  group ID for the ice processors (for rheology)
+   INTEGER         ::   ndim_rank_ice   !  number of 'ice' processors
+   INTEGER         ::   n_ice_root      !  number (in the comm_ice) of proc 0 in the ice comm
    INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::   nrank_ice     ! dimension ndim_rank_ice
 
    ! variables used for zonal integration
    INTEGER, PUBLIC ::   ncomm_znl       !: communicator made by the processors on the same zonal average
-   LOGICAL, PUBLIC ::   l_znl_root      ! True on the 'left'most processor on the same row
-   INTEGER ::   ngrp_znl        ! group ID for the znl processors
-   INTEGER ::   ndim_rank_znl   ! number of processors on the same zonal average
+   LOGICAL, PUBLIC ::   l_znl_root      !  True on the 'left'most processor on the same row
+   INTEGER         ::   ngrp_znl        !  group ID for the znl processors
+   INTEGER         ::   ndim_rank_znl   !  number of processors on the same zonal average
    INTEGER, DIMENSION(:), ALLOCATABLE, SAVE ::   nrank_znl  ! dimension ndim_rank_znl, number of the procs into the same znl domain
 
    ! North fold condition in mpp_mpi with jpni > 1 (PUBLIC for TAM)
-   INTEGER, PUBLIC ::   ngrp_world        ! group ID for the world processors
-   INTEGER, PUBLIC ::   ngrp_opa          ! group ID for the opa processors
-   INTEGER, PUBLIC ::   ngrp_north        ! group ID for the northern processors (to be fold)
-   INTEGER, PUBLIC ::   ncomm_north       ! communicator made by the processors belonging to ngrp_north
-   INTEGER, PUBLIC ::   ndim_rank_north   ! number of 'sea' processor in the northern line (can be /= jpni !)
-   INTEGER, PUBLIC ::   njmppmax          ! value of njmpp for the processors of the northern line
-   INTEGER, PUBLIC ::   north_root        ! number (in the comm_opa) of proc 0 in the northern comm
-   INTEGER, DIMENSION(:), ALLOCATABLE, SAVE, PUBLIC ::   nrank_north   ! dimension ndim_rank_north
+   INTEGER, PUBLIC ::   ngrp_world        !: group ID for the world processors
+   INTEGER, PUBLIC ::   ngrp_opa          !: group ID for the opa processors
+   INTEGER, PUBLIC ::   ngrp_north        !: group ID for the northern processors (to be fold)
+   INTEGER, PUBLIC ::   ncomm_north       !: communicator made by the processors belonging to ngrp_north
+   INTEGER, PUBLIC ::   ndim_rank_north   !: number of 'sea' processor in the northern line (can be /= jpni !)
+   INTEGER, PUBLIC ::   njmppmax          !: value of njmpp for the processors of the northern line
+   INTEGER, PUBLIC ::   north_root        !: number (in the comm_opa) of proc 0 in the northern comm
+   INTEGER, PUBLIC, DIMENSION(:), ALLOCATABLE, SAVE ::   nrank_north   !: dimension ndim_rank_north
 
    ! Type of send : standard, buffered, immediate
-   CHARACTER(len=1), PUBLIC ::   cn_mpi_send   ! type od mpi send/recieve (S=standard, B=bsend, I=isend)
-   LOGICAL, PUBLIC          ::   l_isend = .FALSE.   ! isend use indicator (T if cn_mpi_send='I')
-   INTEGER, PUBLIC          ::   nn_buffer     ! size of the buffer in case of mpi_bsend
-
-   REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE :: tampon  ! buffer in case of bsend
-
-   LOGICAL, PUBLIC                                  ::   ln_nnogather       ! namelist control of northfold comms
-   LOGICAL, PUBLIC                                  ::   l_north_nogather = .FALSE.  ! internal control of northfold comms
-   INTEGER, PUBLIC                                  ::   ityp
+   CHARACTER(len=1), PUBLIC ::   cn_mpi_send         !: type od mpi send/recieve (S=standard, B=bsend, I=isend)
+   LOGICAL         , PUBLIC ::   l_isend = .FALSE.   !: isend use indicator (T if cn_mpi_send='I')
+   INTEGER         , PUBLIC ::   nn_buffer           !: size of the buffer in case of mpi_bsend
+
+   REAL(wp), DIMENSION(:), ALLOCATABLE, SAVE ::   tampon   ! buffer in case of bsend
+
+   LOGICAL, PUBLIC ::   ln_nnogather                !: namelist control of northfold comms
+   LOGICAL, PUBLIC ::   l_north_nogather = .FALSE.  !: internal control of northfold comms
+
    !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! NEMO/OPA 4.0 , NEMO Consortium (2017)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -178,6 +179,5 @@
 CONTAINS
 
-
-   FUNCTION mynode( ldtxt, ldname, kumnam_ref , kumnam_cfg , kumond , kstop, localComm )
+   FUNCTION mynode( ldtxt, ldname, kumnam_ref, kumnam_cfg, kumond, kstop, localComm )
       !!----------------------------------------------------------------------
       !!                  ***  routine mynode  ***
@@ -204,18 +204,17 @@
       WRITE(ldtxt(ii),*) '~~~~~~ '                                                        ;   ii = ii + 1
       !
-
       REWIND( kumnam_ref )              ! Namelist nammpp in reference namelist: mpi variables
       READ  ( kumnam_ref, nammpp, IOSTAT = ios, ERR = 901)
 901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammpp in reference namelist', lwp )
-
+      !
       REWIND( kumnam_cfg )              ! Namelist nammpp in configuration namelist: mpi variables
       READ  ( kumnam_cfg, nammpp, IOSTAT = ios, ERR = 902 )
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammpp in configuration namelist', lwp )
-
+      !
       !                              ! control print
       WRITE(ldtxt(ii),*) '   Namelist nammpp'                                             ;   ii = ii + 1
       WRITE(ldtxt(ii),*) '      mpi send type          cn_mpi_send = ', cn_mpi_send       ;   ii = ii + 1
       WRITE(ldtxt(ii),*) '      size exported buffer   nn_buffer   = ', nn_buffer,' bytes';   ii = ii + 1
-
+      !
 #if defined key_agrif
       IF( .NOT. Agrif_Root() ) THEN
@@ -225,12 +224,10 @@
       ENDIF
 #endif
-
-      IF(jpnij < 1)THEN
-         ! If jpnij is not specified in namelist then we calculate it - this
-         ! means there will be no land cutting out.
-         jpnij = jpni * jpnj
-      END IF
-
-      IF( (jpni < 1) .OR. (jpnj < 1) )THEN
+      !
+      IF( jpnij < 1 ) THEN         ! If jpnij is not specified in namelist then we calculate it
+         jpnij = jpni * jpnj       ! this means there will be no land cutting out.
+      ENDIF
+
+      IF( jpni < 1 .OR. jpnj < 1  ) THEN
          WRITE(ldtxt(ii),*) '      jpni, jpnj and jpnij will be calculated automatically' ;   ii = ii + 1
       ELSE
@@ -238,5 +235,5 @@
          WRITE(ldtxt(ii),*) '      processor grid extent in j         jpnj = ',jpnj       ;   ii = ii + 1
          WRITE(ldtxt(ii),*) '      number of local domains           jpnij = ',jpnij      ;   ii = ii + 1
-      END IF
+      ENDIF
 
       WRITE(ldtxt(ii),*) '      avoid use of mpi_allgather at the north fold  ln_nnogather = ', ln_nnogather  ; ii = ii + 1
@@ -268,5 +265,6 @@
             kstop = kstop + 1
          END SELECT
-      ELSE IF ( PRESENT(localComm) .and. .not. mpi_was_called ) THEN
+         !
+      ELSE IF ( PRESENT(localComm) .AND. .NOT. mpi_was_called ) THEN
          WRITE(ldtxt(ii),*) ' lib_mpp: You cannot provide a local communicator '          ;   ii = ii + 1
          WRITE(ldtxt(ii),*) '          without calling MPI_Init before ! '                ;   ii = ii + 1
@@ -309,5 +307,5 @@
 
 #if defined key_agrif
-      IF (Agrif_Root()) THEN
+      IF( Agrif_Root() ) THEN
          CALL Agrif_MPI_Init(mpi_comm_opa)
       ELSE
@@ -335,4 +333,213 @@
       !!
       !! ** Purpose :   Message passing manadgement
+      !!
+      !! ** Method  :   Use mppsend and mpprecv function for passing mask
+      !!              between processors following neighboring subdomains.
+      !!            domain parameters
+      !!                    nlci   : first dimension of the local subdomain
+      !!                    nlcj   : second dimension of the local subdomain
+      !!                    nbondi : mark for "east-west local boundary"
+      !!                    nbondj : mark for "north-south local boundary"
+      !!                    noea   : number for local neighboring processors
+      !!                    nowe   : number for local neighboring processors
+      !!                    noso   : number for local neighboring processors
+      !!                    nono   : number for local neighboring processors
+      !!
+      !! ** Action  :   ptab with update value at its periphery
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptab     ! 3D array on which the boundary condition is applied
+      CHARACTER(len=1)          , INTENT(in   ) ::   cd_type  ! nature of ptab array grid-points
+      REAL(wp)                  , INTENT(in   ) ::   psgn     ! sign used across the north fold boundary
+      CHARACTER(len=3), OPTIONAL, INTENT(in   ) ::   cd_mpp   ! fill the overlap area only
+      REAL(wp)        , OPTIONAL, INTENT(in   ) ::   pval     ! background value (used at closed boundaries)
+      !
+      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
+      INTEGER  ::   ipk                        ! 3rd dimension of the input array
+      INTEGER  ::   imigr, iihom, ijhom        ! temporary integers
+      INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
+      REAL(wp) ::   zland
+      INTEGER , DIMENSION(MPI_STATUS_SIZE)      ::   ml_stat        ! for key_mpi_isend
+      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   zt3ns, zt3sn   ! 3d for north-south & south-north
+      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   zt3ew, zt3we   ! 3d for east-west & west-east
+      !!----------------------------------------------------------------------
+      !
+      ipk = SIZE( ptab, 3 )
+      !
+      ALLOCATE( zt3ns(jpi,jprecj,ipk,2), zt3sn(jpi,jprecj,ipk,2),   &
+         &      zt3ew(jpj,jpreci,ipk,2), zt3we(jpj,jpreci,ipk,2)  )
+
+      !
+      IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value
+      ELSE                         ;   zland = 0._wp     ! zero by default
+      ENDIF
+
+      ! 1. standard boundary treatment
+      ! ------------------------------
+      IF( PRESENT( cd_mpp ) ) THEN      ! only fill added line/raw with existing values
+         !
+         ! WARNING ptab is defined only between nld and nle
+         DO jk = 1, ipk
+            DO jj = nlcj+1, jpj                 ! added line(s)   (inner only)
+               ptab(nldi  :nlei  , jj          ,jk) = ptab(nldi:nlei,     nlej,jk)
+               ptab(1     :nldi-1, jj          ,jk) = ptab(nldi     ,     nlej,jk)
+               ptab(nlei+1:nlci  , jj          ,jk) = ptab(     nlei,     nlej,jk)
+            END DO
+            DO ji = nlci+1, jpi                 ! added column(s) (full)
+               ptab(ji           ,nldj  :nlej  ,jk) = ptab(     nlei,nldj:nlej,jk)
+               ptab(ji           ,1     :nldj-1,jk) = ptab(     nlei,nldj     ,jk)
+               ptab(ji           ,nlej+1:jpj   ,jk) = ptab(     nlei,     nlej,jk)
+            END DO
+         END DO
+         !
+      ELSE                              ! standard close or cyclic treatment
+         !
+         !                                   ! East-West boundaries
+         !                                        !* Cyclic
+         IF( nbondi == 2 .AND. (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
+            ptab( 1 ,:,:) = ptab(jpim1,:,:)
+            ptab(jpi,:,:) = ptab(  2  ,:,:)
+         ELSE                                     !* closed
+            IF( .NOT. cd_type == 'F' )   ptab(     1       :jpreci,:,:) = zland    ! south except F-point
+                                         ptab(nlci-jpreci+1:jpi   ,:,:) = zland    ! north
+         ENDIF
+         !                                   ! North-South boundaries
+         !                                        !* cyclic (only with no mpp j-split)
+         IF( nbondj == 2 .AND. jperio == 7 ) THEN 
+            ptab(:,1 , :) = ptab(:, jpjm1,:)
+            ptab(:,jpj,:) = ptab(:,     2,:)
+         ELSE                                     !* closed
+            IF( .NOT. cd_type == 'F' )   ptab(:,     1       :jprecj,:) = zland       ! south except F-point
+                                         ptab(:,nlcj-jprecj+1:jpj   ,:) = zland       ! north
+         ENDIF
+         !
+      ENDIF
+
+      ! 2. East and west directions exchange
+      ! ------------------------------------
+      ! we play with the neigbours AND the row number because of the periodicity
+      !
+      SELECT CASE ( nbondi )      ! Read Dirichlet lateral conditions
+      CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case)
+         iihom = nlci-nreci
+         DO jl = 1, jpreci
+            zt3ew(:,jl,:,1) = ptab(jpreci+jl,:,:)
+            zt3we(:,jl,:,1) = ptab(iihom +jl,:,:)
+         END DO
+      END SELECT
+      !
+      !                           ! Migrations
+      imigr = jpreci * jpj * ipk
+      !
+      SELECT CASE ( nbondi )
+      CASE ( -1 )
+         CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req1 )
+         CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea )
+         IF(l_isend)   CALL mpi_wait(ml_req1, ml_stat, ml_err)
+      CASE ( 0 )
+         CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 )
+         CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req2 )
+         CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea )
+         CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe )
+         IF(l_isend)   CALL mpi_wait(ml_req1, ml_stat, ml_err)
+         IF(l_isend)   CALL mpi_wait(ml_req2, ml_stat, ml_err)
+      CASE ( 1 )
+         CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 )
+         CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe )
+         IF(l_isend)   CALL mpi_wait(ml_req1, ml_stat, ml_err)
+      END SELECT
+      !
+      !                           ! Write Dirichlet lateral conditions
+      iihom = nlci-jpreci
+      !
+      SELECT CASE ( nbondi )
+      CASE ( -1 )
+         DO jl = 1, jpreci
+            ptab(iihom+jl,:,:) = zt3ew(:,jl,:,2)
+         END DO
+      CASE ( 0 )
+         DO jl = 1, jpreci
+            ptab(jl      ,:,:) = zt3we(:,jl,:,2)
+            ptab(iihom+jl,:,:) = zt3ew(:,jl,:,2)
+         END DO
+      CASE ( 1 )
+         DO jl = 1, jpreci
+            ptab(jl      ,:,:) = zt3we(:,jl,:,2)
+         END DO
+      END SELECT
+
+      ! 3. North and south directions
+      ! -----------------------------
+      ! always closed : we play only with the neigbours
+      !
+      IF( nbondj /= 2 ) THEN      ! Read Dirichlet lateral conditions
+         ijhom = nlcj-nrecj
+         DO jl = 1, jprecj
+            zt3sn(:,jl,:,1) = ptab(:,ijhom +jl,:)
+            zt3ns(:,jl,:,1) = ptab(:,jprecj+jl,:)
+         END DO
+      ENDIF
+      !
+      !                           ! Migrations
+      imigr = jprecj * jpi * ipk
+      !
+      SELECT CASE ( nbondj )
+      CASE ( -1 )
+         CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req1 )
+         CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono )
+         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+      CASE ( 0 )
+         CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 )
+         CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req2 )
+         CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono )
+         CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso )
+         IF(l_isend)   CALL mpi_wait(ml_req1, ml_stat, ml_err)
+         IF(l_isend)   CALL mpi_wait(ml_req2, ml_stat, ml_err)
+      CASE ( 1 )
+         CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 )
+         CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso )
+         IF(l_isend)   CALL mpi_wait(ml_req1, ml_stat, ml_err )
+      END SELECT
+      !
+      !                           ! Write Dirichlet lateral conditions
+      ijhom = nlcj-jprecj
+      !
+      SELECT CASE ( nbondj )
+      CASE ( -1 )
+         DO jl = 1, jprecj
+            ptab(:,ijhom+jl,:) = zt3ns(:,jl,:,2)
+         END DO
+      CASE ( 0 )
+         DO jl = 1, jprecj
+            ptab(:,jl      ,:) = zt3sn(:,jl,:,2)
+            ptab(:,ijhom+jl,:) = zt3ns(:,jl,:,2)
+         END DO
+      CASE ( 1 )
+         DO jl = 1, jprecj
+            ptab(:,jl,:) = zt3sn(:,jl,:,2)
+         END DO
+      END SELECT
+
+      ! 4. north fold treatment
+      ! -----------------------
+      !
+      IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN
+         !
+         SELECT CASE ( jpni )
+         CASE ( 1 )     ;   CALL lbc_nfd      ( ptab, cd_type, psgn )   ! only 1 northern proc, no mpp
+         CASE DEFAULT   ;   CALL mpp_lbc_north( ptab, cd_type, psgn )   ! for all northern procs.
+         END SELECT
+         !
+      ENDIF
+      !
+      DEALLOCATE( zt3ns, zt3sn, zt3ew, zt3we )
+      !
+   END SUBROUTINE mpp_lnk_3d
+
+
+   SUBROUTINE mpp_lnk_2d_multiple( pt2d_array, type_array, psgn_array, kfld, cd_mpp, pval )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mpp_lnk_2d_multiple  ***
+      !!
+      !! ** Purpose :   Message passing management for multiple 2d arrays
       !!
       !! ** Method  :   Use mppsend and mpprecv function for passing mask
@@ -347,28 +554,23 @@
       !!                    noso   : number for local neighboring processors
       !!                    nono   : number for local neighboring processors
-      !!
-      !! ** Action  :   ptab with update value at its periphery
-      !!
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptab     ! 3D array on which the boundary condition is applied
-      CHARACTER(len=1)                , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points
-      !                                                             ! = T , U , V , F , W points
-      REAL(wp)                        , INTENT(in   ) ::   psgn     ! =-1 the sign change across the north fold boundary
-      !                                                             ! =  1. , the sign is kept
-      CHARACTER(len=3), OPTIONAL      , INTENT(in   ) ::   cd_mpp   ! fill the overlap area only
-      REAL(wp)        , OPTIONAL      , INTENT(in   ) ::   pval     ! background value (used at closed boundaries)
-      !
-      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
+      !!----------------------------------------------------------------------
+      TYPE( arrayptr ), DIMENSION(:), INTENT(inout) ::   pt2d_array   ! pointer array of 2D fields 
+      CHARACTER(len=1), DIMENSION(:), INTENT(in   ) ::   type_array   ! nature of pt2d_array grid-points
+      REAL(wp)        , DIMENSION(:), INTENT(in   ) ::   psgn_array   ! sign used across the north fold boundary
+      INTEGER                       , INTENT(in   ) ::   kfld         ! number of pt2d arrays
+      CHARACTER(len=3), OPTIONAL    , INTENT(in   ) ::   cd_mpp       ! fill the overlap area only
+      REAL(wp)        , OPTIONAL    , INTENT(in   ) ::   pval         ! background value (used at closed boundaries)
+      !
+      INTEGER  ::   ji, jj, jl, jf   ! dummy loop indices
       INTEGER  ::   imigr, iihom, ijhom        ! temporary integers
       INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
       REAL(wp) ::   zland
-      INTEGER , DIMENSION(MPI_STATUS_SIZE)      ::   ml_stat        ! for key_mpi_isend
-      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   zt3ns, zt3sn   ! 3d for north-south & south-north
-      REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   zt3ew, zt3we   ! 3d for east-west & west-east
-      !!----------------------------------------------------------------------
-      
-      ALLOCATE( zt3ns(jpi,jprecj,jpk,2), zt3sn(jpi,jprecj,jpk,2),   &
-         &      zt3ew(jpj,jpreci,jpk,2), zt3we(jpj,jpreci,jpk,2)  )
-
+      INTEGER , DIMENSION(MPI_STATUS_SIZE)    ::   ml_stat       ! for key_mpi_isend
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ns, zt2sn   ! 2d for north-south & south-north
+      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ew, zt2we   ! 2d for east-west & west-east
+      !!----------------------------------------------------------------------
+      !
+      ALLOCATE( zt2ns(jpi,jprecj,2*kfld), zt2sn(jpi,jprecj,2*kfld),  &
+         &      zt2ew(jpj,jpreci,2*kfld), zt2we(jpj,jpreci,2*kfld)   )
       !
       IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value
@@ -378,41 +580,43 @@
       ! 1. standard boundary treatment
       ! ------------------------------
-      IF( PRESENT( cd_mpp ) ) THEN      ! only fill added line/raw with existing values
-         !
-         ! WARNING ptab is defined only between nld and nle
-         DO jk = 1, jpk
+      !
+      !First Array
+      DO jf = 1 , kfld
+         IF( PRESENT( cd_mpp ) ) THEN      ! only fill added line/raw with existing values
+            !
+            ! WARNING pt2d is defined only between nld and nle
             DO jj = nlcj+1, jpj                 ! added line(s)   (inner only)
-               ptab(nldi  :nlei  , jj          ,jk) = ptab(nldi:nlei,     nlej,jk)
-               ptab(1     :nldi-1, jj          ,jk) = ptab(nldi     ,     nlej,jk)
-               ptab(nlei+1:nlci  , jj          ,jk) = ptab(     nlei,     nlej,jk)
+               pt2d_array(jf)%pt2d(nldi  :nlei  , jj) = pt2d_array(jf)%pt2d(nldi:nlei, nlej)
+               pt2d_array(jf)%pt2d(1     :nldi-1, jj) = pt2d_array(jf)%pt2d(nldi     , nlej)
+               pt2d_array(jf)%pt2d(nlei+1:nlci  , jj) = pt2d_array(jf)%pt2d(     nlei, nlej) 
             END DO
             DO ji = nlci+1, jpi                 ! added column(s) (full)
-               ptab(ji           ,nldj  :nlej  ,jk) = ptab(     nlei,nldj:nlej,jk)
-               ptab(ji           ,1     :nldj-1,jk) = ptab(     nlei,nldj     ,jk)
-               ptab(ji           ,nlej+1:jpj   ,jk) = ptab(     nlei,     nlej,jk)
+               pt2d_array(jf)%pt2d(ji, nldj  :nlej  ) = pt2d_array(jf)%pt2d(nlei, nldj:nlej)
+               pt2d_array(jf)%pt2d(ji, 1     :nldj-1) = pt2d_array(jf)%pt2d(nlei, nldj     )
+               pt2d_array(jf)%pt2d(ji, nlej+1:jpj   ) = pt2d_array(jf)%pt2d(nlei,      nlej)
             END DO
-         END DO
-         !
-      ELSE                              ! standard close or cyclic treatment
-         !
-         !                                   ! East-West boundaries
-         !                                        !* Cyclic east-west
-         IF( nbondi == 2 .AND. (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
-            ptab( 1 ,:,:) = ptab(jpim1,:,:)
-            ptab(jpi,:,:) = ptab(  2  ,:,:)
-         ELSE                                     !* closed
-            IF( .NOT. cd_type == 'F' )   ptab(     1       :jpreci,:,:) = zland    ! south except F-point
-                                         ptab(nlci-jpreci+1:jpi   ,:,:) = zland    ! north
+            !
+         ELSE                              ! standard close or cyclic treatment
+            !
+            !                                   ! East-West boundaries
+            IF( nbondi == 2 .AND.   &                !* Cyclic
+               &    (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
+               pt2d_array(jf)%pt2d(  1  , : ) = pt2d_array(jf)%pt2d( jpim1, : )                             ! west
+               pt2d_array(jf)%pt2d( jpi , : ) = pt2d_array(jf)%pt2d(   2  , : )                             ! east
+            ELSE                                     !* Closed
+               IF( .NOT. type_array(jf) == 'F' )   pt2d_array(jf)%pt2d(            1 : jpreci,:) = zland    ! south except F-point
+                                                   pt2d_array(jf)%pt2d(nlci-jpreci+1 : jpi   ,:) = zland    ! north
+            ENDIF
+            !                                   ! North-South boundaries
+            !                                        !* Cyclic
+            IF( nbondj == 2 .AND. jperio == 7 ) THEN
+               pt2d_array(jf)%pt2d(:,  1  ) =   pt2d_array(jf)%pt2d(:, jpjm1 )
+               pt2d_array(jf)%pt2d(:, jpj ) =   pt2d_array(jf)%pt2d(:,   2   )          
+            ELSE                                     !* Closed             
+               IF( .NOT. type_array(jf) == 'F' )   pt2d_array(jf)%pt2d(:,             1:jprecj ) = zland    ! south except F-point
+                                                   pt2d_array(jf)%pt2d(:, nlcj-jprecj+1:jpj    ) = zland    ! north
+            ENDIF
          ENDIF
-                                          ! North-south cyclic
-         IF ( nbondj == 2 .AND. jperio == 7 )    THEN !* cyclic north south only with no mpp split in latitude
-            ptab(:,1 , :) = ptab(:, jpjm1,:)
-            ptab(:,jpj,:) = ptab(:,     2,:)
-         ELSE   !                                   ! North-South boundaries (closed)
-            IF( .NOT. cd_type == 'F' )   ptab(:,     1       :jprecj,:) = zland       ! south except F-point
-                                         ptab(:,nlcj-jprecj+1:jpj   ,:) = zland       ! north
-         ENDIF
-         !
-      ENDIF
+      END DO
 
       ! 2. East and west directions exchange
@@ -420,127 +624,202 @@
       ! we play with the neigbours AND the row number because of the periodicity
       !
-      SELECT CASE ( nbondi )      ! Read Dirichlet lateral conditions
-      CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case)
-         iihom = nlci-nreci
-         DO jl = 1, jpreci
-            zt3ew(:,jl,:,1) = ptab(jpreci+jl,:,:)
-            zt3we(:,jl,:,1) = ptab(iihom +jl,:,:)
-         END DO
-      END SELECT
+      DO jf = 1 , kfld
+         SELECT CASE ( nbondi )      ! Read Dirichlet lateral conditions
+         CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case)
+            iihom = nlci-nreci
+            DO jl = 1, jpreci
+               zt2ew( : , jl , jf ) = pt2d_array(jf)%pt2d( jpreci+jl , : )
+               zt2we( : , jl , jf ) = pt2d_array(jf)%pt2d( iihom +jl , : )
+            END DO
+         END SELECT
+      END DO
       !
       !                           ! Migrations
-      imigr = jpreci * jpj * jpk
+      imigr = jpreci * jpj
       !
       SELECT CASE ( nbondi )
       CASE ( -1 )
-         CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req1 )
-         CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea )
-         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
-      CASE ( 0 )
-         CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 )
-         CALL mppsend( 2, zt3we(1,1,1,1), imigr, noea, ml_req2 )
-         CALL mpprecv( 1, zt3ew(1,1,1,2), imigr, noea )
-         CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe )
-         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
-         IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err)
-      CASE ( 1 )
-         CALL mppsend( 1, zt3ew(1,1,1,1), imigr, nowe, ml_req1 )
-         CALL mpprecv( 2, zt3we(1,1,1,2), imigr, nowe )
-         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+         CALL mppsend( 2, zt2we(1,1,1), kfld*imigr, noea, ml_req1 )
+         CALL mpprecv( 1, zt2ew(1,1,kfld+1), kfld*imigr, noea )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+      CASE ( 0 )
+         CALL mppsend( 1, zt2ew(1,1,1), kfld*imigr, nowe, ml_req1 )
+         CALL mppsend( 2, zt2we(1,1,1), kfld*imigr, noea, ml_req2 )
+         CALL mpprecv( 1, zt2ew(1,1,kfld+1), kfld*imigr, noea )
+         CALL mpprecv( 2, zt2we(1,1,kfld+1), kfld*imigr, nowe )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
+         IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
+      CASE ( 1 )
+         CALL mppsend( 1, zt2ew(1,1,1), kfld*imigr, nowe, ml_req1 )
+         CALL mpprecv( 2, zt2we(1,1,kfld+1), kfld*imigr, nowe )
+         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       END SELECT
       !
       !                           ! Write Dirichlet lateral conditions
-      iihom = nlci-jpreci
-      !
-      SELECT CASE ( nbondi )
-      CASE ( -1 )
-         DO jl = 1, jpreci
-            ptab(iihom+jl,:,:) = zt3ew(:,jl,:,2)
-         END DO
-      CASE ( 0 )
-         DO jl = 1, jpreci
-            ptab(jl      ,:,:) = zt3we(:,jl,:,2)
-            ptab(iihom+jl,:,:) = zt3ew(:,jl,:,2)
-         END DO
-      CASE ( 1 )
-         DO jl = 1, jpreci
-            ptab(jl      ,:,:) = zt3we(:,jl,:,2)
-         END DO
-      END SELECT
-
+      iihom = nlci - jpreci
+      !
+
+      DO jf = 1 , kfld
+         SELECT CASE ( nbondi )
+         CASE ( -1 )
+            DO jl = 1, jpreci
+               pt2d_array(jf)%pt2d( iihom+jl ,:) = zt2ew(:,jl,kfld+jf)
+            END DO
+         CASE ( 0 )
+            DO jl = 1, jpreci
+               pt2d_array(jf)%pt2d(       jl ,:) = zt2we(:,jl,kfld+jf)
+               pt2d_array(jf)%pt2d( iihom+jl ,:) = zt2ew(:,jl,kfld+jf)
+            END DO
+         CASE ( 1 )
+            DO jl = 1, jpreci
+               pt2d_array(jf)%pt2d( jl ,:)= zt2we(:,jl,kfld+jf)
+            END DO
+         END SELECT
+      END DO
+      
       ! 3. North and south directions
       ! -----------------------------
       ! always closed : we play only with the neigbours
       !
-      IF( nbondj /= 2 ) THEN      ! Read Dirichlet lateral conditions
-         ijhom = nlcj-nrecj
-         DO jl = 1, jprecj
-            zt3sn(:,jl,:,1) = ptab(:,ijhom +jl,:)
-            zt3ns(:,jl,:,1) = ptab(:,jprecj+jl,:)
-         END DO
-      ENDIF
+      !First Array
+      DO jf = 1 , kfld
+         IF( nbondj /= 2 ) THEN      ! Read Dirichlet lateral conditions
+            ijhom = nlcj-nrecj
+            DO jl = 1, jprecj
+               zt2sn(:,jl,jf) = pt2d_array(jf)%pt2d(:, ijhom +jl )
+               zt2ns(:,jl,jf) = pt2d_array(jf)%pt2d(:, jprecj+jl )
+            END DO
+         ENDIF
+      END DO
       !
       !                           ! Migrations
-      imigr = jprecj * jpi * jpk
+      imigr = jprecj * jpi
       !
       SELECT CASE ( nbondj )
       CASE ( -1 )
-         CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req1 )
-         CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono )
-         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
-      CASE ( 0 )
-         CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 )
-         CALL mppsend( 4, zt3sn(1,1,1,1), imigr, nono, ml_req2 )
-         CALL mpprecv( 3, zt3ns(1,1,1,2), imigr, nono )
-         CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso )
-         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
-         IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err)
-      CASE ( 1 )
-         CALL mppsend( 3, zt3ns(1,1,1,1), imigr, noso, ml_req1 )
-         CALL mpprecv( 4, zt3sn(1,1,1,2), imigr, noso )
-         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+         CALL mppsend( 4, zt2sn(1,1,     1), kfld*imigr, nono, ml_req1 )
+         CALL mpprecv( 3, zt2ns(1,1,kfld+1), kfld*imigr, nono )
+         IF(l_isend)   CALL mpi_wait( ml_req1, ml_stat, ml_err )
+      CASE ( 0 )
+         CALL mppsend( 3, zt2ns(1,1,     1), kfld*imigr, noso, ml_req1 )
+         CALL mppsend( 4, zt2sn(1,1,     1), kfld*imigr, nono, ml_req2 )
+         CALL mpprecv( 3, zt2ns(1,1,kfld+1), kfld*imigr, nono )
+         CALL mpprecv( 4, zt2sn(1,1,kfld+1), kfld*imigr, noso )
+         IF(l_isend)   CALL mpi_wait(ml_req1,ml_stat,ml_err)
+         IF(l_isend)   CALL mpi_wait(ml_req2,ml_stat,ml_err)
+      CASE ( 1 )
+         CALL mppsend( 3, zt2ns(1,1,     1), kfld*imigr, noso, ml_req1 )
+         CALL mpprecv( 4, zt2sn(1,1,kfld+1), kfld*imigr, noso )
+         IF(l_isend)   CALL mpi_wait(ml_req1,ml_stat,ml_err)
       END SELECT
       !
       !                           ! Write Dirichlet lateral conditions
-      ijhom = nlcj-jprecj
-      !
-      SELECT CASE ( nbondj )
-      CASE ( -1 )
-         DO jl = 1, jprecj
-            ptab(:,ijhom+jl,:) = zt3ns(:,jl,:,2)
-         END DO
-      CASE ( 0 )
-         DO jl = 1, jprecj
-            ptab(:,jl      ,:) = zt3sn(:,jl,:,2)
-            ptab(:,ijhom+jl,:) = zt3ns(:,jl,:,2)
-         END DO
-      CASE ( 1 )
-         DO jl = 1, jprecj
-            ptab(:,jl,:) = zt3sn(:,jl,:,2)
-         END DO
-      END SELECT
-
+      ijhom = nlcj - jprecj
+      !
+      DO jf = 1 , kfld
+         SELECT CASE ( nbondj )
+         CASE ( -1 )
+            DO jl = 1, jprecj
+               pt2d_array(jf)%pt2d(:, ijhom+jl ) = zt2ns(:,jl, kfld+jf )
+            END DO
+         CASE ( 0 )
+            DO jl = 1, jprecj
+               pt2d_array(jf)%pt2d(:,       jl ) = zt2sn(:,jl, kfld+jf )
+               pt2d_array(jf)%pt2d(:, ijhom+jl ) = zt2ns(:,jl, kfld+jf )
+            END DO
+         CASE ( 1 )
+            DO jl = 1, jprecj
+               pt2d_array(jf)%pt2d(:,       jl ) = zt2sn(:,jl, kfld+jf )
+            END DO
+         END SELECT
+      END DO
+      
       ! 4. north fold treatment
       ! -----------------------
       !
-      IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN
+      IF( npolj /= 0 .AND. .NOT.PRESENT(cd_mpp) ) THEN
          !
          SELECT CASE ( jpni )
-         CASE ( 1 )     ;   CALL lbc_nfd      ( ptab, cd_type, psgn )   ! only 1 northern proc, no mpp
-         CASE DEFAULT   ;   CALL mpp_lbc_north( ptab, cd_type, psgn )   ! for all northern procs.
+         CASE ( 1 )  
+            DO jf = 1, kfld  
+               CALL lbc_nfd( pt2d_array(jf)%pt2d(:,:), type_array(jf), psgn_array(jf) )  ! only 1 northern proc, no mpp
+            END DO
+         CASE DEFAULT   
+            CALL mpp_lbc_north_2d_multiple( pt2d_array, type_array, psgn_array, kfld )   ! for all northern procs.
          END SELECT
          !
       ENDIF
       !
-      DEALLOCATE( zt3ns, zt3sn, zt3ew, zt3we )
-      !
-   END SUBROUTINE mpp_lnk_3d
-
-
-   SUBROUTINE mpp_lnk_2d_multiple( pt2d_array , type_array , psgn_array , num_fields , cd_mpp, pval )
-      !!----------------------------------------------------------------------
-      !!                  ***  routine mpp_lnk_2d_multiple  ***
-      !!
-      !! ** Purpose :   Message passing management for multiple 2d arrays
+      DEALLOCATE( zt2ns, zt2sn, zt2ew, zt2we )
+      !
+   END SUBROUTINE mpp_lnk_2d_multiple
+
+   
+   SUBROUTINE load_array( pt2d, cd_type, psgn, pt2d_array, type_array, psgn_array, kfld )
+      !!---------------------------------------------------------------------
+      REAL(wp)        , DIMENSION(:,:), TARGET, INTENT(inout) ::   pt2d         ! 2D array on which the boundary condition is applied
+      CHARACTER(len=1)                        , INTENT(in   ) ::   cd_type      ! nature of ptab array grid-points
+      REAL(wp)                                , INTENT(in   ) ::   psgn         ! sign used across the north fold boundary
+      TYPE(arrayptr)  , DIMENSION(:)          , INTENT(inout) ::   pt2d_array   ! 
+      CHARACTER(len=1), DIMENSION(:)          , INTENT(inout) ::   type_array   ! nature of pt2d_array array grid-points
+      REAL(wp)        , DIMENSION(:)          , INTENT(inout) ::   psgn_array   ! sign used across the north fold boundary
+      INTEGER                                 , INTENT(inout) ::   kfld         !
+      !!---------------------------------------------------------------------
+      !
+      kfld                  =  kfld + 1
+      pt2d_array(kfld)%pt2d => pt2d
+      type_array(kfld)      =  cd_type
+      psgn_array(kfld)      =  psgn
+      !
+   END SUBROUTINE load_array
+   
+   
+   SUBROUTINE mpp_lnk_2d_9( pt2dA, cd_typeA, psgnA, pt2dB, cd_typeB, psgnB, pt2dC, cd_typeC, psgnC   &
+      &                   , pt2dD, cd_typeD, psgnD, pt2dE, cd_typeE, psgnE, pt2dF, cd_typeF, psgnF   &
+      &                   , pt2dG, cd_typeG, psgnG, pt2dH, cd_typeH, psgnH, pt2dI, cd_typeI, psgnI, cd_mpp, pval)
+      !!---------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), TARGET          , INTENT(inout) ::   pt2dA    ! 2D arrays on which the lbc is applied
+      REAL(wp), DIMENSION(jpi,jpj), TARGET, OPTIONAL, INTENT(inout) ::   pt2dB , pt2dC , pt2dD , pt2dE
+      REAL(wp), DIMENSION(jpi,jpj), TARGET, OPTIONAL, INTENT(inout) ::   pt2dF , pt2dG , pt2dH , pt2dI 
+      CHARACTER(len=1)                              , INTENT(in   ) ::   cd_typeA ! nature of pt2D. array grid-points
+      CHARACTER(len=1)                    , OPTIONAL, INTENT(in   ) ::   cd_typeB , cd_typeC , cd_typeD , cd_typeE
+      CHARACTER(len=1)                    , OPTIONAL, INTENT(in   ) ::   cd_typeF , cd_typeG , cd_typeH , cd_typeI
+      REAL(wp)                                      , INTENT(in   ) ::   psgnA    ! sign used across the north fold
+      REAL(wp)                            , OPTIONAL, INTENT(in   ) ::   psgnB , psgnC , psgnD , psgnE
+      REAL(wp)                            , OPTIONAL, INTENT(in   ) ::   psgnF , psgnG , psgnH , psgnI   
+      CHARACTER(len=3)                    , OPTIONAL, INTENT(in   ) ::   cd_mpp   ! fill the overlap area only
+      REAL(wp)                            , OPTIONAL, INTENT(in   ) ::   pval     ! background value (used at closed boundaries)
+      !!
+      INTEGER :: kfld
+      TYPE(arrayptr)   , DIMENSION(9) ::   pt2d_array 
+      CHARACTER(len=1) , DIMENSION(9) ::   type_array    ! define the nature of ptab array grid-points
+      REAL(wp)         , DIMENSION(9) ::   psgn_array    ! sign used across the north fold boundary
+      !!---------------------------------------------------------------------
+      !
+      kfld = 0
+      !
+      !                 ! Load the first array
+      CALL load_array( pt2dA, cd_typeA, psgnA, pt2d_array, type_array, psgn_array, kfld )
+      !
+      !                 ! Look if more arrays are added
+      IF( PRESENT(psgnB) )   CALL load_array( pt2dB, cd_typeB, psgnB, pt2d_array, type_array, psgn_array, kfld )
+      IF( PRESENT(psgnC) )   CALL load_array( pt2dC, cd_typeC, psgnC, pt2d_array, type_array, psgn_array, kfld )
+      IF( PRESENT(psgnD) )   CALL load_array( pt2dD, cd_typeD, psgnD, pt2d_array, type_array, psgn_array, kfld )
+      IF( PRESENT(psgnE) )   CALL load_array( pt2dE, cd_typeE, psgnE, pt2d_array, type_array, psgn_array, kfld )
+      IF( PRESENT(psgnF) )   CALL load_array( pt2dF, cd_typeF, psgnF, pt2d_array, type_array, psgn_array, kfld )
+      IF( PRESENT(psgnG) )   CALL load_array( pt2dG, cd_typeG, psgnG, pt2d_array, type_array, psgn_array, kfld )
+      IF( PRESENT(psgnH) )   CALL load_array( pt2dH, cd_typeH, psgnH, pt2d_array, type_array, psgn_array, kfld )
+      IF( PRESENT(psgnI) )   CALL load_array( pt2dI, cd_typeI, psgnI, pt2d_array, type_array, psgn_array, kfld )
+      !
+      CALL mpp_lnk_2d_multiple( pt2d_array, type_array, psgn_array, kfld, cd_mpp,pval )
+      !
+   END SUBROUTINE mpp_lnk_2d_9
+
+
+   SUBROUTINE mpp_lnk_2d( pt2d, cd_type, psgn, cd_mpp, pval )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mpp_lnk_2d  ***
+      !!
+      !! ** Purpose :   Message passing manadgement for 2d array
       !!
       !! ** Method  :   Use mppsend and mpprecv function for passing mask
@@ -555,27 +834,23 @@
       !!                    noso   : number for local neighboring processors
       !!                    nono   : number for local neighboring processors
-      !!----------------------------------------------------------------------
-      CHARACTER(len=1), DIMENSION(:), INTENT(in   ) ::   type_array   ! define the nature of ptab array grid-points
-      !                                                               ! = T , U , V , F , W and I points
-      REAL(wp)        , DIMENSION(:), INTENT(in   ) ::   psgn_array   ! =-1 the sign change across the north fold boundary
-      !                                                               ! =  1. , the sign is kept
-      CHARACTER(len=3), OPTIONAL    , INTENT(in   ) ::   cd_mpp       ! fill the overlap area only
-      REAL(wp)        , OPTIONAL    , INTENT(in   ) ::   pval         ! background value (used at closed boundaries)
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pt2d     ! 2D array on which the boundary condition is applied
+      CHARACTER(len=1)            , INTENT(in   ) ::   cd_type  ! nature of ptab array grid-points
+      REAL(wp)                    , INTENT(in   ) ::   psgn     ! sign used across the north fold boundary
+      CHARACTER(len=3), OPTIONAL  , INTENT(in   ) ::   cd_mpp   ! fill the overlap area only
+      REAL(wp)        , OPTIONAL  , INTENT(in   ) ::   pval     ! background value (used at closed boundaries)
       !!
       INTEGER  ::   ji, jj, jl   ! dummy loop indices
-      INTEGER  ::   ii    !!MULTI SEND DUMMY LOOP INDICES
       INTEGER  ::   imigr, iihom, ijhom        ! temporary integers
       INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
-      INTEGER :: num_fields
-      TYPE( arrayptr ), DIMENSION(:) :: pt2d_array
       REAL(wp) ::   zland
-      INTEGER , DIMENSION(MPI_STATUS_SIZE)    ::   ml_stat       ! for key_mpi_isend
+      INTEGER, DIMENSION(MPI_STATUS_SIZE)     ::   ml_stat       ! for key_mpi_isend
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ns, zt2sn   ! 2d for north-south & south-north
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ew, zt2we   ! 2d for east-west & west-east
-
-      !!----------------------------------------------------------------------
-      !
-      ALLOCATE( zt2ns(jpi,jprecj,2*num_fields), zt2sn(jpi,jprecj,2*num_fields),  &
-         &      zt2ew(jpj,jpreci,2*num_fields), zt2we(jpj,jpreci,2*num_fields)   )
+      !!----------------------------------------------------------------------
+      !
+      ALLOCATE( zt2ns(jpi,jprecj,2), zt2sn(jpi,jprecj,2),  &
+         &      zt2ew(jpj,jpreci,2), zt2we(jpj,jpreci,2)   )
       !
       IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value
@@ -586,43 +861,39 @@
       ! ------------------------------
       !
-      !First Array
-      DO ii = 1 , num_fields
-         IF( PRESENT( cd_mpp ) ) THEN      ! only fill added line/raw with existing values
-            !
-            ! WARNING pt2d is defined only between nld and nle
-            DO jj = nlcj+1, jpj                 ! added line(s)   (inner only)
-               pt2d_array(ii)%pt2d(nldi  :nlei  , jj) = pt2d_array(ii)%pt2d(nldi:nlei, nlej)
-               pt2d_array(ii)%pt2d(1     :nldi-1, jj) = pt2d_array(ii)%pt2d(nldi     , nlej)
-               pt2d_array(ii)%pt2d(nlei+1:nlci  , jj) = pt2d_array(ii)%pt2d(     nlei, nlej) 
-            END DO
-            DO ji = nlci+1, jpi                 ! added column(s) (full)
-               pt2d_array(ii)%pt2d(ji, nldj  :nlej  ) = pt2d_array(ii)%pt2d(nlei, nldj:nlej)
-               pt2d_array(ii)%pt2d(ji, 1     :nldj-1) = pt2d_array(ii)%pt2d(nlei, nldj     )
-               pt2d_array(ii)%pt2d(ji, nlej+1:jpj   ) = pt2d_array(ii)%pt2d(nlei,      nlej)
-            END DO
-            !
-         ELSE                              ! standard close or cyclic treatment
-            !
-            !                                   ! East-West boundaries
-            IF( nbondi == 2 .AND.   &                ! Cyclic east-west
-               &    (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
-               pt2d_array(ii)%pt2d(  1  , : ) = pt2d_array(ii)%pt2d( jpim1, : )                                    ! west
-               pt2d_array(ii)%pt2d( jpi , : ) = pt2d_array(ii)%pt2d(   2  , : )                                    ! east
-            ELSE                                     ! closed
-               IF( .NOT. type_array(ii) == 'F' )   pt2d_array(ii)%pt2d(            1 : jpreci,:) = zland    ! south except F-point
-                                                   pt2d_array(ii)%pt2d(nlci-jpreci+1 : jpi   ,:) = zland    ! north
-            ENDIF
-                                                ! Noth-South boundaries
-            IF ( nbondj == 2 .AND. jperio == 7 )    THEN !* cyclic north south
-               pt2d_array(ii)%pt2d(:, 1   ) =   pt2d_array(ii)%pt2d(:, jpjm1 )
-               pt2d_array(ii)%pt2d(:, jpj ) =   pt2d_array(ii)%pt2d(:, 2 )          
-            ELSE   !             
-               !                                   ! North-South boundaries (closed)
-               IF( .NOT. type_array(ii) == 'F' )   pt2d_array(ii)%pt2d(:,             1:jprecj ) = zland    ! south except F-point
-                                                   pt2d_array(ii)%pt2d(:, nlcj-jprecj+1:jpj    ) = zland    ! north
-            !
-            ENDIF
-          ENDIF
-      END DO
+      IF( PRESENT( cd_mpp ) ) THEN      ! only fill added line/raw with existing values
+         !
+         ! WARNING pt2d is defined only between nld and nle
+         DO jj = nlcj+1, jpj                 ! added line(s)   (inner only)
+            pt2d(nldi  :nlei  , jj          ) = pt2d(nldi:nlei,     nlej)
+            pt2d(1     :nldi-1, jj          ) = pt2d(nldi     ,     nlej)
+            pt2d(nlei+1:nlci  , jj          ) = pt2d(     nlei,     nlej)
+         END DO
+         DO ji = nlci+1, jpi                 ! added column(s) (full)
+            pt2d(ji           ,nldj  :nlej  ) = pt2d(     nlei,nldj:nlej)
+            pt2d(ji           ,1     :nldj-1) = pt2d(     nlei,nldj     )
+            pt2d(ji           ,nlej+1:jpj   ) = pt2d(     nlei,     nlej)
+         END DO
+         !
+      ELSE                              ! standard close or cyclic treatment
+         !
+         !                                   ! East-West boundaries
+         IF( nbondi == 2 .AND.   &                !* cyclic
+            &    (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
+            pt2d( 1 ,:) = pt2d(jpim1,:)                                          ! west
+            pt2d(jpi,:) = pt2d(  2  ,:)                                          ! east
+         ELSE                                     !* closed
+            IF( .NOT. cd_type == 'F' )   pt2d(     1       :jpreci,:) = zland    ! south except F-point
+                                         pt2d(nlci-jpreci+1:jpi   ,:) = zland    ! north
+         ENDIF
+         !                                   ! North-South boundaries
+         !                                        !* cyclic
+         IF( nbondj == 2 .AND. jperio == 7 ) THEN
+            pt2d(:,  1 ) = pt2d(:,jpjm1)
+            pt2d(:, jpj) = pt2d(:,    2)
+         ELSE                                     !* closed
+            IF( .NOT. cd_type == 'F' )   pt2d(:,     1       :jprecj) = zland    !south except F-point
+                                         pt2d(:,nlcj-jprecj+1:jpj   ) = zland    ! north
+         ENDIF
+      ENDIF
 
       ! 2. East and west directions exchange
@@ -630,14 +901,12 @@
       ! we play with the neigbours AND the row number because of the periodicity
       !
-      DO ii = 1 , num_fields
-         SELECT CASE ( nbondi )      ! Read Dirichlet lateral conditions
-         CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case)
-            iihom = nlci-nreci
-            DO jl = 1, jpreci
-               zt2ew( : , jl , ii ) = pt2d_array(ii)%pt2d( jpreci+jl , : )
-               zt2we( : , jl , ii ) = pt2d_array(ii)%pt2d( iihom +jl , : )
-            END DO
-         END SELECT
-      END DO
+      SELECT CASE ( nbondi )      ! Read Dirichlet lateral conditions
+      CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case)
+         iihom = nlci-nreci
+         DO jl = 1, jpreci
+            zt2ew(:,jl,1) = pt2d(jpreci+jl,:)
+            zt2we(:,jl,1) = pt2d(iihom +jl,:)
+         END DO
+      END SELECT
       !
       !                           ! Migrations
@@ -646,17 +915,17 @@
       SELECT CASE ( nbondi )
       CASE ( -1 )
-         CALL mppsend( 2, zt2we(1,1,1), num_fields*imigr, noea, ml_req1 )
-         CALL mpprecv( 1, zt2ew(1,1,num_fields+1), num_fields*imigr, noea )
+         CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req1 )
+         CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       CASE ( 0 )
-         CALL mppsend( 1, zt2ew(1,1,1), num_fields*imigr, nowe, ml_req1 )
-         CALL mppsend( 2, zt2we(1,1,1), num_fields*imigr, noea, ml_req2 )
-         CALL mpprecv( 1, zt2ew(1,1,num_fields+1), num_fields*imigr, noea )
-         CALL mpprecv( 2, zt2we(1,1,num_fields+1), num_fields*imigr, nowe )
+         CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 )
+         CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req2 )
+         CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea )
+         CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
          IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
       CASE ( 1 )
-         CALL mppsend( 1, zt2ew(1,1,1), num_fields*imigr, nowe, ml_req1 )
-         CALL mpprecv( 2, zt2we(1,1,num_fields+1), num_fields*imigr, nowe )
+         CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 )
+         CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       END SELECT
@@ -665,37 +934,31 @@
       iihom = nlci - jpreci
       !
-
-      DO ii = 1 , num_fields
-         SELECT CASE ( nbondi )
-         CASE ( -1 )
-            DO jl = 1, jpreci
-               pt2d_array(ii)%pt2d( iihom+jl , : ) = zt2ew(:,jl,num_fields+ii)
-            END DO
-         CASE ( 0 )
-            DO jl = 1, jpreci
-               pt2d_array(ii)%pt2d( jl , : ) = zt2we(:,jl,num_fields+ii)
-               pt2d_array(ii)%pt2d( iihom+jl , : ) = zt2ew(:,jl,num_fields+ii)
-            END DO
-         CASE ( 1 )
-            DO jl = 1, jpreci
-               pt2d_array(ii)%pt2d( jl , : )= zt2we(:,jl,num_fields+ii)
-            END DO
-         END SELECT
-      END DO
-      
+      SELECT CASE ( nbondi )
+      CASE ( -1 )
+         DO jl = 1, jpreci
+            pt2d(iihom+jl,:) = zt2ew(:,jl,2)
+         END DO
+      CASE ( 0 )
+         DO jl = 1, jpreci
+            pt2d(jl      ,:) = zt2we(:,jl,2)
+            pt2d(iihom+jl,:) = zt2ew(:,jl,2)
+         END DO
+      CASE ( 1 )
+         DO jl = 1, jpreci
+            pt2d(jl      ,:) = zt2we(:,jl,2)
+         END DO
+      END SELECT
+
       ! 3. North and south directions
       ! -----------------------------
       ! always closed : we play only with the neigbours
       !
-      !First Array
-      DO ii = 1 , num_fields
-         IF( nbondj /= 2 ) THEN      ! Read Dirichlet lateral conditions
-            ijhom = nlcj-nrecj
-            DO jl = 1, jprecj
-               zt2sn(:,jl , ii) = pt2d_array(ii)%pt2d( : , ijhom +jl )
-               zt2ns(:,jl , ii) = pt2d_array(ii)%pt2d( : , jprecj+jl )
-            END DO
-         ENDIF
-      END DO
+      IF( nbondj /= 2 ) THEN      ! Read Dirichlet lateral conditions
+         ijhom = nlcj-nrecj
+         DO jl = 1, jprecj
+            zt2sn(:,jl,1) = pt2d(:,ijhom +jl)
+            zt2ns(:,jl,1) = pt2d(:,jprecj+jl)
+         END DO
+      ENDIF
       !
       !                           ! Migrations
@@ -704,17 +967,17 @@
       SELECT CASE ( nbondj )
       CASE ( -1 )
-         CALL mppsend( 4, zt2sn(1,1,1), num_fields*imigr, nono, ml_req1 )
-         CALL mpprecv( 3, zt2ns(1,1,num_fields+1), num_fields*imigr, nono )
+         CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req1 )
+         CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       CASE ( 0 )
-         CALL mppsend( 3, zt2ns(1,1,1), num_fields*imigr, noso, ml_req1 )
-         CALL mppsend( 4, zt2sn(1,1,1), num_fields*imigr, nono, ml_req2 )
-         CALL mpprecv( 3, zt2ns(1,1,num_fields+1), num_fields*imigr, nono )
-         CALL mpprecv( 4, zt2sn(1,1,num_fields+1), num_fields*imigr, noso )
+         CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 )
+         CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req2 )
+         CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono )
+         CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
          IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
       CASE ( 1 )
-         CALL mppsend( 3, zt2ns(1,1,1), num_fields*imigr, noso, ml_req1 )
-         CALL mpprecv( 4, zt2sn(1,1,num_fields+1), num_fields*imigr, noso )
+         CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 )
+         CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso )
          IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
       END SELECT
@@ -723,114 +986,42 @@
       ijhom = nlcj - jprecj
       !
-
-      DO ii = 1 , num_fields
-         !First Array
-         SELECT CASE ( nbondj )
-         CASE ( -1 )
-            DO jl = 1, jprecj
-               pt2d_array(ii)%pt2d( : , ijhom+jl ) = zt2ns( : , jl , num_fields+ii )
-            END DO
-         CASE ( 0 )
-            DO jl = 1, jprecj
-               pt2d_array(ii)%pt2d( : , jl ) = zt2sn( : , jl , num_fields + ii)
-               pt2d_array(ii)%pt2d( : , ijhom + jl ) = zt2ns( : , jl , num_fields + ii )
-            END DO
-         CASE ( 1 )
-            DO jl = 1, jprecj
-               pt2d_array(ii)%pt2d( : , jl ) = zt2sn( : , jl , num_fields + ii )
-            END DO
-         END SELECT
-      END DO
-      
+      SELECT CASE ( nbondj )
+      CASE ( -1 )
+         DO jl = 1, jprecj
+            pt2d(:,ijhom+jl) = zt2ns(:,jl,2)
+         END DO
+      CASE ( 0 )
+         DO jl = 1, jprecj
+            pt2d(:,jl      ) = zt2sn(:,jl,2)
+            pt2d(:,ijhom+jl) = zt2ns(:,jl,2)
+         END DO
+      CASE ( 1 )
+         DO jl = 1, jprecj
+            pt2d(:,jl      ) = zt2sn(:,jl,2)
+         END DO
+      END SELECT
+
       ! 4. north fold treatment
       ! -----------------------
       !
-         !First Array
       IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN
          !
          SELECT CASE ( jpni )
-         CASE ( 1 )     ;   
-             DO ii = 1 , num_fields  
-                       CALL lbc_nfd      ( pt2d_array(ii)%pt2d( : , : ), type_array(ii) , psgn_array(ii) )   ! only 1 northern proc, no mpp
-             END DO
-         CASE DEFAULT   ;   CALL mpp_lbc_north_2d_multiple( pt2d_array, type_array, psgn_array, num_fields )   ! for all northern procs.
+         CASE ( 1 )     ;   CALL lbc_nfd      ( pt2d, cd_type, psgn )   ! only 1 northern proc, no mpp
+         CASE DEFAULT   ;   CALL mpp_lbc_north( pt2d, cd_type, psgn )   ! for all northern procs.
          END SELECT
          !
       ENDIF
-        !
       !
       DEALLOCATE( zt2ns, zt2sn, zt2ew, zt2we )
       !
-   END SUBROUTINE mpp_lnk_2d_multiple
-
-   
-   SUBROUTINE load_array( pt2d, cd_type, psgn, pt2d_array, type_array, psgn_array, num_fields )
-      !!---------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), TARGET, INTENT(inout) ::   pt2d    ! Second 2D array on which the boundary condition is applied
-      CHARACTER(len=1)                    , INTENT(in   ) ::   cd_type ! define the nature of ptab array grid-points
-      REAL(wp)                            , INTENT(in   ) ::   psgn    ! =-1 the sign change across the north fold boundary
-      TYPE(arrayptr)   , DIMENSION(9) ::   pt2d_array
-      CHARACTER(len=1) , DIMENSION(9) ::   type_array    ! define the nature of ptab array grid-points
-      REAL(wp)         , DIMENSION(9) ::   psgn_array    ! =-1 the sign change across the north fold boundary
-      INTEGER                            , INTENT (inout) :: num_fields 
-      !!---------------------------------------------------------------------
-      num_fields = num_fields + 1
-      pt2d_array(num_fields)%pt2d => pt2d
-      type_array(num_fields)      =  cd_type
-      psgn_array(num_fields)      =  psgn
-   END SUBROUTINE load_array
-   
-   
-   SUBROUTINE mpp_lnk_2d_9( pt2dA, cd_typeA, psgnA, pt2dB, cd_typeB, psgnB, pt2dC, cd_typeC, psgnC   &
-      &                   , pt2dD, cd_typeD, psgnD, pt2dE, cd_typeE, psgnE, pt2dF, cd_typeF, psgnF   &
-      &                   , pt2dG, cd_typeG, psgnG, pt2dH, cd_typeH, psgnH, pt2dI, cd_typeI, psgnI, cd_mpp, pval)
-      !!---------------------------------------------------------------------
-      ! Second 2D array on which the boundary condition is applied
-      REAL(wp), DIMENSION(jpi,jpj), TARGET          , INTENT(inout) ::   pt2dA    
-      REAL(wp), DIMENSION(jpi,jpj), TARGET, OPTIONAL, INTENT(inout) ::   pt2dB , pt2dC , pt2dD , pt2dE
-      REAL(wp), DIMENSION(jpi,jpj), TARGET, OPTIONAL, INTENT(inout) ::   pt2dF , pt2dG , pt2dH , pt2dI 
-      ! define the nature of ptab array grid-points
-      CHARACTER(len=1)                              , INTENT(in   ) ::   cd_typeA
-      CHARACTER(len=1)                    , OPTIONAL, INTENT(in   ) ::   cd_typeB , cd_typeC , cd_typeD , cd_typeE
-      CHARACTER(len=1)                    , OPTIONAL, INTENT(in   ) ::   cd_typeF , cd_typeG , cd_typeH , cd_typeI
-      ! =-1 the sign change across the north fold boundary
-      REAL(wp)                                      , INTENT(in   ) ::   psgnA    
-      REAL(wp)                            , OPTIONAL, INTENT(in   ) ::   psgnB , psgnC , psgnD , psgnE
-      REAL(wp)                            , OPTIONAL, INTENT(in   ) ::   psgnF , psgnG , psgnH , psgnI   
-      CHARACTER(len=3)                    , OPTIONAL, INTENT(in   ) ::   cd_mpp   ! fill the overlap area only
-      REAL(wp)                            , OPTIONAL, INTENT(in   ) ::   pval     ! background value (used at closed boundaries)
-      !!
-      TYPE(arrayptr)   , DIMENSION(9) ::   pt2d_array 
-      CHARACTER(len=1) , DIMENSION(9) ::   type_array    ! define the nature of ptab array grid-points
-      !                                                         ! = T , U , V , F , W and I points
-      REAL(wp)         , DIMENSION(9) ::   psgn_array    ! =-1 the sign change across the north fold boundary
-      INTEGER :: num_fields
-      !!---------------------------------------------------------------------
-      !
-      num_fields = 0
-      !
-      ! Load the first array
-      CALL load_array( pt2dA, cd_typeA, psgnA, pt2d_array, type_array, psgn_array, num_fields )
-      !
-      ! Look if more arrays are added
-      IF( PRESENT(psgnB) )   CALL load_array(pt2dB,cd_typeB,psgnB,pt2d_array, type_array, psgn_array,num_fields)
-      IF( PRESENT(psgnC) )   CALL load_array(pt2dC,cd_typeC,psgnC,pt2d_array, type_array, psgn_array,num_fields)
-      IF( PRESENT(psgnD) )   CALL load_array(pt2dD,cd_typeD,psgnD,pt2d_array, type_array, psgn_array,num_fields)
-      IF( PRESENT(psgnE) )   CALL load_array(pt2dE,cd_typeE,psgnE,pt2d_array, type_array, psgn_array,num_fields)
-      IF( PRESENT(psgnF) )   CALL load_array(pt2dF,cd_typeF,psgnF,pt2d_array, type_array, psgn_array,num_fields)
-      IF( PRESENT(psgnG) )   CALL load_array(pt2dG,cd_typeG,psgnG,pt2d_array, type_array, psgn_array,num_fields)
-      IF( PRESENT(psgnH) )   CALL load_array(pt2dH,cd_typeH,psgnH,pt2d_array, type_array, psgn_array,num_fields)
-      IF( PRESENT(psgnI) )   CALL load_array(pt2dI,cd_typeI,psgnI,pt2d_array, type_array, psgn_array,num_fields)
-      !
-      CALL mpp_lnk_2d_multiple( pt2d_array, type_array, psgn_array, num_fields, cd_mpp,pval )
-      !
-   END SUBROUTINE mpp_lnk_2d_9
-
-
-   SUBROUTINE mpp_lnk_2d( pt2d, cd_type, psgn, cd_mpp, pval )
-      !!----------------------------------------------------------------------
-      !!                  ***  routine mpp_lnk_2d  ***
-      !!
-      !! ** Purpose :   Message passing manadgement for 2d array
+   END SUBROUTINE mpp_lnk_2d
+
+
+   SUBROUTINE mpp_lnk_3d_gather( ptab1, cd_type1, ptab2, cd_type2, psgn )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mpp_lnk_3d_gather  ***
+      !!
+      !! ** Purpose :   Message passing manadgement for two 3D arrays
       !!
       !! ** Method  :   Use mppsend and mpprecv function for passing mask
@@ -846,220 +1037,15 @@
       !!                    nono   : number for local neighboring processors
       !!
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pt2d     ! 2D array on which the boundary condition is applied
-      CHARACTER(len=1)            , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points
-      !                                                         ! = T , U , V , F , W and I points
-      REAL(wp)                    , INTENT(in   ) ::   psgn     ! =-1 the sign change across the north fold boundary
-      !                                                         ! =  1. , the sign is kept
-      CHARACTER(len=3), OPTIONAL  , INTENT(in   ) ::   cd_mpp   ! fill the overlap area only
-      REAL(wp)        , OPTIONAL  , INTENT(in   ) ::   pval     ! background value (used at closed boundaries)
-      !!
-      INTEGER  ::   ji, jj, jl   ! dummy loop indices
-      INTEGER  ::   imigr, iihom, ijhom        ! temporary integers
-      INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
-      REAL(wp) ::   zland
-      INTEGER, DIMENSION(MPI_STATUS_SIZE)     ::   ml_stat       ! for key_mpi_isend
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ns, zt2sn   ! 2d for north-south & south-north
-      REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ew, zt2we   ! 2d for east-west & west-east
-      !!----------------------------------------------------------------------
-      !
-      ALLOCATE( zt2ns(jpi,jprecj,2), zt2sn(jpi,jprecj,2),  &
-         &      zt2ew(jpj,jpreci,2), zt2we(jpj,jpreci,2)   )
-      !
-      IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value
-      ELSE                         ;   zland = 0._wp     ! zero by default
-      ENDIF
-
-      ! 1. standard boundary treatment
-      ! ------------------------------
-      !
-      IF( PRESENT( cd_mpp ) ) THEN      ! only fill added line/raw with existing values
-         !
-         ! WARNING pt2d is defined only between nld and nle
-         DO jj = nlcj+1, jpj                 ! added line(s)   (inner only)
-            pt2d(nldi  :nlei  , jj          ) = pt2d(nldi:nlei,     nlej)
-            pt2d(1     :nldi-1, jj          ) = pt2d(nldi     ,     nlej)
-            pt2d(nlei+1:nlci  , jj          ) = pt2d(     nlei,     nlej)
-         END DO
-         DO ji = nlci+1, jpi                 ! added column(s) (full)
-            pt2d(ji           ,nldj  :nlej  ) = pt2d(     nlei,nldj:nlej)
-            pt2d(ji           ,1     :nldj-1) = pt2d(     nlei,nldj     )
-            pt2d(ji           ,nlej+1:jpj   ) = pt2d(     nlei,     nlej)
-         END DO
-         !
-      ELSE                              ! standard close or cyclic treatment
-         !
-         !                                   ! East-West boundaries
-         IF( nbondi == 2 .AND.   &                ! Cyclic east-west
-            &    (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
-            pt2d( 1 ,:) = pt2d(jpim1,:)                                    ! west
-            pt2d(jpi,:) = pt2d(  2  ,:)                                    ! east
-         ELSE                                     ! closed
-            IF( .NOT. cd_type == 'F' )   pt2d(     1       :jpreci,:) = zland    ! south except F-point
-                                         pt2d(nlci-jpreci+1:jpi   ,:) = zland    ! north
-         ENDIF
-                                            ! North-South boudaries
-         IF ( nbondj == 2 .AND. jperio == 7 )    THEN !* cyclic north south
-            pt2d(:,  1 ) = pt2d(:,jpjm1)
-            pt2d(:, jpj) = pt2d(:,    2)
-         ELSE    
-         !                                   ! North-South boundaries (closed)
-            IF( .NOT. cd_type == 'F' )   pt2d(:,     1       :jprecj) = zland    !south except F-point
-                                         pt2d(:,nlcj-jprecj+1:jpj   ) = zland    ! north
-         ENDIF     
-      ENDIF
-
-      ! 2. East and west directions exchange
-      ! ------------------------------------
-      ! we play with the neigbours AND the row number because of the periodicity
-      !
-      SELECT CASE ( nbondi )      ! Read Dirichlet lateral conditions
-      CASE ( -1, 0, 1 )                ! all exept 2 (i.e. close case)
-         iihom = nlci-nreci
-         DO jl = 1, jpreci
-            zt2ew(:,jl,1) = pt2d(jpreci+jl,:)
-            zt2we(:,jl,1) = pt2d(iihom +jl,:)
-         END DO
-      END SELECT
-      !
-      !                           ! Migrations
-      imigr = jpreci * jpj
-      !
-      SELECT CASE ( nbondi )
-      CASE ( -1 )
-         CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req1 )
-         CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea )
-         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
-      CASE ( 0 )
-         CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 )
-         CALL mppsend( 2, zt2we(1,1,1), imigr, noea, ml_req2 )
-         CALL mpprecv( 1, zt2ew(1,1,2), imigr, noea )
-         CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe )
-         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
-         IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
-      CASE ( 1 )
-         CALL mppsend( 1, zt2ew(1,1,1), imigr, nowe, ml_req1 )
-         CALL mpprecv( 2, zt2we(1,1,2), imigr, nowe )
-         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
-      END SELECT
-      !
-      !                           ! Write Dirichlet lateral conditions
-      iihom = nlci - jpreci
-      !
-      SELECT CASE ( nbondi )
-      CASE ( -1 )
-         DO jl = 1, jpreci
-            pt2d(iihom+jl,:) = zt2ew(:,jl,2)
-         END DO
-      CASE ( 0 )
-         DO jl = 1, jpreci
-            pt2d(jl      ,:) = zt2we(:,jl,2)
-            pt2d(iihom+jl,:) = zt2ew(:,jl,2)
-         END DO
-      CASE ( 1 )
-         DO jl = 1, jpreci
-            pt2d(jl      ,:) = zt2we(:,jl,2)
-         END DO
-      END SELECT
-
-
-      ! 3. North and south directions
-      ! -----------------------------
-      ! always closed : we play only with the neigbours
-      !
-      IF( nbondj /= 2 ) THEN      ! Read Dirichlet lateral conditions
-         ijhom = nlcj-nrecj
-         DO jl = 1, jprecj
-            zt2sn(:,jl,1) = pt2d(:,ijhom +jl)
-            zt2ns(:,jl,1) = pt2d(:,jprecj+jl)
-         END DO
-      ENDIF
-      !
-      !                           ! Migrations
-      imigr = jprecj * jpi
-      !
-      SELECT CASE ( nbondj )
-      CASE ( -1 )
-         CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req1 )
-         CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono )
-         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
-      CASE ( 0 )
-         CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 )
-         CALL mppsend( 4, zt2sn(1,1,1), imigr, nono, ml_req2 )
-         CALL mpprecv( 3, zt2ns(1,1,2), imigr, nono )
-         CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso )
-         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
-         IF(l_isend) CALL mpi_wait(ml_req2,ml_stat,ml_err)
-      CASE ( 1 )
-         CALL mppsend( 3, zt2ns(1,1,1), imigr, noso, ml_req1 )
-         CALL mpprecv( 4, zt2sn(1,1,2), imigr, noso )
-         IF(l_isend) CALL mpi_wait(ml_req1,ml_stat,ml_err)
-      END SELECT
-      !
-      !                           ! Write Dirichlet lateral conditions
-      ijhom = nlcj - jprecj
-      !
-      SELECT CASE ( nbondj )
-      CASE ( -1 )
-         DO jl = 1, jprecj
-            pt2d(:,ijhom+jl) = zt2ns(:,jl,2)
-         END DO
-      CASE ( 0 )
-         DO jl = 1, jprecj
-            pt2d(:,jl      ) = zt2sn(:,jl,2)
-            pt2d(:,ijhom+jl) = zt2ns(:,jl,2)
-         END DO
-      CASE ( 1 )
-         DO jl = 1, jprecj
-            pt2d(:,jl      ) = zt2sn(:,jl,2)
-         END DO
-      END SELECT
-
-
-      ! 4. north fold treatment
-      ! -----------------------
-      !
-      IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN
-         !
-         SELECT CASE ( jpni )
-         CASE ( 1 )     ;   CALL lbc_nfd      ( pt2d, cd_type, psgn )   ! only 1 northern proc, no mpp
-         CASE DEFAULT   ;   CALL mpp_lbc_north( pt2d, cd_type, psgn )   ! for all northern procs.
-         END SELECT
-         !
-      ENDIF
-      !
-      DEALLOCATE( zt2ns, zt2sn, zt2ew, zt2we )
-      !
-   END SUBROUTINE mpp_lnk_2d
-
-
-   SUBROUTINE mpp_lnk_3d_gather( ptab1, cd_type1, ptab2, cd_type2, psgn )
-      !!----------------------------------------------------------------------
-      !!                  ***  routine mpp_lnk_3d_gather  ***
-      !!
-      !! ** Purpose :   Message passing manadgement for two 3D arrays
-      !!
-      !! ** Method  :   Use mppsend and mpprecv function for passing mask
-      !!      between processors following neighboring subdomains.
-      !!            domain parameters
-      !!                    nlci   : first dimension of the local subdomain
-      !!                    nlcj   : second dimension of the local subdomain
-      !!                    nbondi : mark for "east-west local boundary"
-      !!                    nbondj : mark for "north-south local boundary"
-      !!                    noea   : number for local neighboring processors
-      !!                    nowe   : number for local neighboring processors
-      !!                    noso   : number for local neighboring processors
-      !!                    nono   : number for local neighboring processors
-      !!
       !! ** Action  :   ptab1 and ptab2  with update value at its periphery
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptab1     ! first and second 3D array on which
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptab2     ! the boundary condition is applied
-      CHARACTER(len=1)                , INTENT(in   ) ::   cd_type1  ! nature of ptab1 and ptab2 arrays
-      CHARACTER(len=1)                , INTENT(in   ) ::   cd_type2  ! i.e. grid-points = T , U , V , F or W points
-      REAL(wp)                        , INTENT(in   ) ::   psgn      ! =-1 the sign change across the north fold boundary
-      !!                                                             ! =  1. , the sign is kept
-      INTEGER  ::   jl   ! dummy loop indices
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptab1     ! 1st 3D array on which the boundary condition is applied
+      CHARACTER(len=1)          , INTENT(in   ) ::   cd_type1  ! nature of ptab1 arrays
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptab2     ! 3nd 3D array on which the boundary condition is applied
+      CHARACTER(len=1)          , INTENT(in   ) ::   cd_type2  ! nature of ptab2 arrays
+      REAL(wp)                  , INTENT(in   ) ::   psgn      ! sign used across the north fold boundary
+      !
+      INTEGER  ::   jl                         ! dummy loop indices
+      INTEGER  ::   ipk                        ! 3rd dimension of the input array
       INTEGER  ::   imigr, iihom, ijhom        ! temporary integers
       INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
@@ -1069,11 +1055,13 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( zt4ns(jpi,jprecj,jpk,2,2), zt4sn(jpi,jprecj,jpk,2,2) ,    &
-         &      zt4ew(jpj,jpreci,jpk,2,2), zt4we(jpj,jpreci,jpk,2,2) )
-      !
+      ipk = SIZE( ptab1, 3 )
+      !
+      ALLOCATE( zt4ns(jpi,jprecj,ipk,2,2), zt4sn(jpi,jprecj,ipk,2,2) ,    &
+         &      zt4ew(jpj,jpreci,ipk,2,2), zt4we(jpj,jpreci,ipk,2,2) )
+
       ! 1. standard boundary treatment
       ! ------------------------------
       !                                      ! East-West boundaries
-      !                                           !* Cyclic east-west
+      !                                           !* Cyclic 
       IF( nbondi == 2 .AND. (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
          ptab1( 1 ,:,:) = ptab1(jpim1,:,:)
@@ -1082,22 +1070,23 @@
          ptab2(jpi,:,:) = ptab2(  2  ,:,:)
       ELSE                                        !* closed
-         IF( .NOT. cd_type1 == 'F' )   ptab1(     1       :jpreci,:,:) = 0.e0    ! south except at F-point
-         IF( .NOT. cd_type2 == 'F' )   ptab2(     1       :jpreci,:,:) = 0.e0
-                                       ptab1(nlci-jpreci+1:jpi   ,:,:) = 0.e0    ! north
-                                       ptab2(nlci-jpreci+1:jpi   ,:,:) = 0.e0
-      ENDIF
-                                            ! North-South boundaries
-      IF ( nbondj == 2 .AND. jperio == 7 )    THEN !* cyclic north south
-         ptab1(:,     1       ,:) = ptab1(: ,  jpjm1 , :)
-         ptab1(:,   jpj       ,:) = ptab1(: ,      2 , :)
-         ptab2(:,     1       ,:) = ptab2(: ,  jpjm1 , :)
-         ptab2(:,   jpj       ,:) = ptab2(: ,      2 , :)
+         IF( .NOT. cd_type1 == 'F' )   ptab1(     1       :jpreci,:,:) = 0._wp   ! south except at F-point
+         IF( .NOT. cd_type2 == 'F' )   ptab2(     1       :jpreci,:,:) = 0._wp
+                                       ptab1(nlci-jpreci+1:jpi   ,:,:) = 0._wp   ! north
+                                       ptab2(nlci-jpreci+1:jpi   ,:,:) = 0._wp
+      ENDIF
+      !                                     ! North-South boundaries
+      !                                           !* cyclic
+      IF( nbondj == 2 .AND. jperio == 7 ) THEN
+         ptab1(:,  1  ,:) = ptab1(:, jpjm1 , :)
+         ptab1(:, jpj ,:) = ptab1(:,   2   , :)
+         ptab2(:,  1  ,:) = ptab2(:, jpjm1 , :)
+         ptab2(:, jpj ,:) = ptab2(:,   2   , :)
       ELSE     
-      !                                      ! North-South boundaries closed
-      IF( .NOT. cd_type1 == 'F' )   ptab1(:,     1       :jprecj,:) = 0.e0    ! south except at F-point
-      IF( .NOT. cd_type2 == 'F' )   ptab2(:,     1       :jprecj,:) = 0.e0
-                                    ptab1(:,nlcj-jprecj+1:jpj   ,:) = 0.e0    ! north
-                                    ptab2(:,nlcj-jprecj+1:jpj   ,:) = 0.e0
-      ENDIF     
+         !                                        !* closed
+         IF( .NOT. cd_type1 == 'F' )   ptab1(:,     1       :jprecj,:) = 0._wp   ! south except at F-point
+         IF( .NOT. cd_type2 == 'F' )   ptab2(:,     1       :jprecj,:) = 0._wp
+                                       ptab1(:,nlcj-jprecj+1:jpj   ,:) = 0._wp   ! north
+                                       ptab2(:,nlcj-jprecj+1:jpj   ,:) = 0._wp
+      ENDIF
 
       ! 2. East and west directions exchange
@@ -1117,5 +1106,5 @@
       !
       !                           ! Migrations
-      imigr = jpreci * jpj * jpk *2
+      imigr = jpreci * jpj * ipk *2
       !
       SELECT CASE ( nbondi )
@@ -1159,5 +1148,4 @@
          END DO
       END SELECT
-
 
       ! 3. North and south directions
@@ -1176,5 +1164,5 @@
       !
       !                           ! Migrations
-      imigr = jprecj * jpi * jpk * 2
+      imigr = jprecj * jpi * ipk * 2
       !
       SELECT CASE ( nbondj )
@@ -1218,5 +1206,4 @@
          END DO
       END SELECT
-
 
       ! 4. north fold treatment
@@ -1284,38 +1271,32 @@
 
 
-      ! 1. standard boundary treatment
+      ! 1. standard boundary treatment   (CAUTION: the order matters Here !!!! )
       ! ------------------------------
-      ! Order matters Here !!!!
-      !
-                                           ! North-South cyclic
-      IF ( nbondj == 2 .AND. jperio == 7 )    THEN !* cyclic north south
-         pt2d(:, 1-jprj:  1     ) = pt2d ( :, jpjm1-jprj:jpjm1)
+      !                                !== North-South boundaries
+      !                                      !* cyclic
+      IF( nbondj == 2 .AND. jperio == 7 ) THEN
+         pt2d(:, 1-jprj:  1     ) = pt2d ( :, jpjm1-jprj:jpjm1 )
          pt2d(:, jpj   :jpj+jprj) = pt2d ( :, 2         :2+jprj)
-      ELSE
-        
-      !                                      !* North-South boundaries (closed)
-      IF( .NOT. cd_type == 'F' )   pt2d(:,  1-jprj   :  jprecj  ) = 0.e0    ! south except at F-point
-                                   pt2d(:,nlcj-jprecj+1:jpj+jprj) = 0.e0    ! north
-      ENDIF
-                                
-      !                                      ! East-West boundaries
-      !                                           !* Cyclic east-west
+      ELSE                                   !* closed
+         IF( .NOT. cd_type == 'F' )   pt2d(:,  1-jprj   :  jprecj  ) = 0._wp     ! south except at F-point
+                                      pt2d(:,nlcj-jprecj+1:jpj+jprj) = 0._wp     ! north
+      ENDIF
+      !                                !== East-West boundaries
+      !                                      !* Cyclic east-west
       IF( nbondi == 2 .AND. (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) ) THEN
-         pt2d(1-jpri:     1    ,:) = pt2d(jpim1-jpri:  jpim1 ,:)       ! east
-         pt2d(   jpi  :jpi+jpri,:) = pt2d(     2      :2+jpri,:)       ! west
-         !
-      ELSE                                        !* closed
-         IF( .NOT. cd_type == 'F' )   pt2d(  1-jpri   :jpreci    ,:) = 0.e0    ! south except at F-point
-                                      pt2d(nlci-jpreci+1:jpi+jpri,:) = 0.e0    ! north
-      ENDIF
-      !
-
+         pt2d(1-jpri:     1    ,:) = pt2d(jpim1-jpri:  jpim1 ,:)              ! east
+         pt2d(   jpi  :jpi+jpri,:) = pt2d(     2      :2+jpri,:)              ! west
+      ELSE                                   !* closed
+         IF( .NOT. cd_type == 'F' )   pt2d(  1-jpri   :jpreci    ,:) = 0._wp  ! south except at F-point
+                                      pt2d(nlci-jpreci+1:jpi+jpri,:) = 0._wp  ! north
+      ENDIF
+      !
       ! north fold treatment
-      ! -----------------------
+      ! --------------------
       IF( npolj /= 0 ) THEN
          !
          SELECT CASE ( jpni )
          CASE ( 1 )     ;   CALL lbc_nfd        ( pt2d(1:jpi,1:jpj+jprj), cd_type, psgn, pr2dj=jprj )
-         CASE DEFAULT   ;   CALL mpp_lbc_north_e( pt2d                    , cd_type, psgn               )
+         CASE DEFAULT   ;   CALL mpp_lbc_north_e( pt2d                  , cd_type, psgn             )
          END SELECT
          !
@@ -1375,5 +1356,4 @@
       END SELECT
 
-
       ! 3. North and south directions
       ! -----------------------------
@@ -1429,4 +1409,5 @@
       !
    END SUBROUTINE mpp_lnk_2d_e
+
 
    SUBROUTINE mpp_lnk_sum_3d( ptab, cd_type, psgn, cd_mpp, pval )
@@ -1452,11 +1433,9 @@
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptab     ! 3D array on which the boundary condition is applied
-      CHARACTER(len=1)                , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points
-      !                                                             ! = T , U , V , F , W points
-      REAL(wp)                        , INTENT(in   ) ::   psgn     ! =-1 the sign change across the north fold boundary
-      !                                                             ! =  1. , the sign is kept
+      CHARACTER(len=1)                , INTENT(in   ) ::   cd_type  !  nature of ptab array grid-points
+      REAL(wp)                        , INTENT(in   ) ::   psgn     ! sign used across the north fold boundary
       CHARACTER(len=3), OPTIONAL      , INTENT(in   ) ::   cd_mpp   ! fill the overlap area only
       REAL(wp)        , OPTIONAL      , INTENT(in   ) ::   pval     ! background value (used at closed boundaries)
-      !!
+      !
       INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
       INTEGER  ::   imigr, iihom, ijhom        ! temporary integers
@@ -1467,13 +1446,11 @@
       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   zt3ns, zt3sn   ! 3d for north-south & south-north
       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE ::   zt3ew, zt3we   ! 3d for east-west & west-east
-
-      !!----------------------------------------------------------------------
-      
+      !!----------------------------------------------------------------------
+      !
       ALLOCATE( zt3ns(jpi,jprecj,jpk,2), zt3sn(jpi,jprecj,jpk,2),   &
          &      zt3ew(jpj,jpreci,jpk,2), zt3we(jpj,jpreci,jpk,2)  )
-
       !
       IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value
-      ELSE                         ;   zland = 0.e0      ! zero by default
+      ELSE                         ;   zland = 0._wp     ! zero by default
       ENDIF
 
@@ -1488,6 +1465,6 @@
       iihom = nlci-jpreci
          DO jl = 1, jpreci
-            zt3ew(:,jl,:,1) = ptab(jl      ,:,:) ; ptab(jl      ,:,:) = 0.0_wp
-            zt3we(:,jl,:,1) = ptab(iihom+jl,:,:) ; ptab(iihom+jl,:,:) = 0.0_wp 
+            zt3ew(:,jl,:,1) = ptab(jl      ,:,:) ; ptab(jl      ,:,:) = 0._wp
+            zt3we(:,jl,:,1) = ptab(iihom+jl,:,:) ; ptab(iihom+jl,:,:) = 0._wp 
          END DO
       END SELECT
@@ -1520,5 +1497,5 @@
       CASE ( -1 )
          DO jl = 1, jpreci
-            ptab(iihom+jl,:,:) = ptab(iihom+jl,:,:) + zt3ew(:,jl,:,2)
+            ptab(iihom +jl,:,:) = ptab(iihom +jl,:,:) + zt3ew(:,jl,:,2)
          END DO
       CASE ( 0 )
@@ -1533,5 +1510,4 @@
       END SELECT
 
-
       ! 3. North and south directions
       ! -----------------------------
@@ -1541,6 +1517,6 @@
          ijhom = nlcj-jprecj
          DO jl = 1, jprecj
-            zt3sn(:,jl,:,1) = ptab(:,ijhom+jl,:) ; ptab(:,ijhom+jl,:) = 0.0_wp
-            zt3ns(:,jl,:,1) = ptab(:,jl      ,:) ; ptab(:,jl      ,:) = 0.0_wp
+            zt3sn(:,jl,:,1) = ptab(:,ijhom+jl,:)   ;   ptab(:,ijhom+jl,:) = 0._wp
+            zt3ns(:,jl,:,1) = ptab(:,jl      ,:)   ;   ptab(:,jl      ,:) = 0._wp
          END DO
       ENDIF
@@ -1586,5 +1562,4 @@
       END SELECT
 
-
       ! 4. north fold treatment
       ! -----------------------
@@ -1602,4 +1577,5 @@
       !
    END SUBROUTINE mpp_lnk_sum_3d
+
 
    SUBROUTINE mpp_lnk_sum_2d( pt2d, cd_type, psgn, cd_mpp, pval )
@@ -1620,11 +1596,8 @@
       !!                    noso   : number for local neighboring processors
       !!                    nono   : number for local neighboring processors
-      !!
       !!----------------------------------------------------------------------
       REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   pt2d     ! 2D array on which the boundary condition is applied
-      CHARACTER(len=1)            , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points
-      !                                                         ! = T , U , V , F , W and I points
-      REAL(wp)                    , INTENT(in   ) ::   psgn     ! =-1 the sign change across the north fold boundary
-      !                                                         ! =  1. , the sign is kept
+      CHARACTER(len=1)            , INTENT(in   ) ::   cd_type  ! nature of pt2d array grid-points
+      REAL(wp)                    , INTENT(in   ) ::   psgn     ! sign used across the north fold boundary
       CHARACTER(len=3), OPTIONAL  , INTENT(in   ) ::   cd_mpp   ! fill the overlap area only
       REAL(wp)        , OPTIONAL  , INTENT(in   ) ::   pval     ! background value (used at closed boundaries)
@@ -1638,13 +1611,11 @@
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ns, zt2sn   ! 2d for north-south & south-north
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::  zt2ew, zt2we   ! 2d for east-west & west-east
-
-      !!----------------------------------------------------------------------
-
+      !!----------------------------------------------------------------------
+      !
       ALLOCATE( zt2ns(jpi,jprecj,2), zt2sn(jpi,jprecj,2),  &
          &      zt2ew(jpj,jpreci,2), zt2we(jpj,jpreci,2)   )
-
       !
       IF( PRESENT( pval ) ) THEN   ;   zland = pval      ! set land value
-      ELSE                         ;   zland = 0.e0      ! zero by default
+      ELSE                         ;   zland = 0._wp     ! zero by default
       ENDIF
 
@@ -1757,5 +1728,4 @@
       END SELECT
 
-
       ! 4. north fold treatment
       ! -----------------------
@@ -1773,4 +1743,5 @@
       !
    END SUBROUTINE mpp_lnk_sum_2d
+
 
    SUBROUTINE mppsend( ktyp, pmess, kbytes, kdest, md_req )
@@ -2015,5 +1986,5 @@
       !!                 ***  routine mppmax_a_real  ***
       !!
-      !! ** Purpose :   Maximum
+      !! ** Purpose :   Maximum of a 1D array
       !!
       !!----------------------------------------------------------------------
@@ -2039,5 +2010,5 @@
       !!                  ***  routine mppmax_real  ***
       !!
-      !! ** Purpose :   Maximum
+      !! ** Purpose :   Maximum for each element of a 1D array
       !!
       !!----------------------------------------------------------------------
@@ -2057,5 +2028,6 @@
    END SUBROUTINE mppmax_real
 
-   SUBROUTINE mppmax_real_multiple( ptab, NUM , kcom  )
+
+   SUBROUTINE mppmax_real_multiple( pt1d, kdim, kcom  )
       !!----------------------------------------------------------------------
       !!                  ***  routine mppmax_real  ***
@@ -2064,19 +2036,16 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(:) ,  INTENT(inout)           ::   ptab   ! ???
-      INTEGER , INTENT(in   )           ::   NUM
-      INTEGER , INTENT(in   ), OPTIONAL ::   kcom   ! ???
+      REAL(wp), DIMENSION(kdim), INTENT(inout) ::   pt1d   ! 1D arrays
+      INTEGER                  , INTENT(in   ) ::   kdim
+      INTEGER , OPTIONAL       , INTENT(in   ) ::   kcom   ! local communicator
       !!
       INTEGER  ::   ierror, localcomm
-      REAL(wp) , POINTER , DIMENSION(:) ::   zwork
-      !!----------------------------------------------------------------------
-      !
-      CALL wrk_alloc(NUM , zwork)
+      !!----------------------------------------------------------------------
+      !
       localcomm = mpi_comm_opa
       IF( PRESENT(kcom) )   localcomm = kcom
       !
-      CALL mpi_allreduce( ptab, zwork, NUM, mpi_double_precision, mpi_max, localcomm, ierror )
-      ptab = zwork
-      CALL wrk_dealloc(NUM , zwork)
+      CALL mpi_allreduce( ptab, zwork, kdim, mpi_double_precision, mpi_max, localcomm, ierror )
+      ptab(:) = zwork(:)
       !
    END SUBROUTINE mppmax_real_multiple
@@ -2251,6 +2220,6 @@
       !!-----------------------------------------------------------------------
       !
-      zmin  = MINVAL( ptab(:,:) , mask= pmask == 1.e0 )
-      ilocs = MINLOC( ptab(:,:) , mask= pmask == 1.e0 )
+      zmin  = MINVAL( ptab(:,:) , mask= pmask == 1._wp )
+      ilocs = MINLOC( ptab(:,:) , mask= pmask == 1._wp )
       !
       ki = ilocs(1) + nimpp - 1
@@ -2279,9 +2248,9 @@
       !!
       !!--------------------------------------------------------------------------
-      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   ptab         ! Local 2D array
-      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pmask        ! Local mask
-      REAL(wp)                         , INTENT(  out) ::   pmin         ! Global minimum of ptab
-      INTEGER                          , INTENT(  out) ::   ki, kj, kk   ! index of minimum in global frame
-      !!
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   ptab         ! Local 2D array
+      REAL(wp), DIMENSION(:,:,:), INTENT(in   ) ::   pmask        ! Local mask
+      REAL(wp)                  , INTENT(  out) ::   pmin         ! Global minimum of ptab
+      INTEGER                   , INTENT(  out) ::   ki, kj, kk   ! index of minimum in global frame
+      !
       INTEGER  ::   ierror
       REAL(wp) ::   zmin     ! local minimum
@@ -2290,6 +2259,6 @@
       !!-----------------------------------------------------------------------
       !
-      zmin  = MINVAL( ptab(:,:,:) , mask= pmask == 1.e0 )
-      ilocs = MINLOC( ptab(:,:,:) , mask= pmask == 1.e0 )
+      zmin  = MINVAL( ptab(:,:,:) , mask= pmask == 1._wp )
+      ilocs = MINLOC( ptab(:,:,:) , mask= pmask == 1._wp )
       !
       ki = ilocs(1) + nimpp - 1
@@ -2297,6 +2266,6 @@
       kk = ilocs(3)
       !
-      zain(1,:)=zmin
-      zain(2,:)=ki+10000.*kj+100000000.*kk
+      zain(1,:) = zmin
+      zain(2,:) = ki + 10000.*kj + 100000000.*kk
       !
       CALL MPI_ALLREDUCE( zain,zaout, 1, MPI_2DOUBLE_PRECISION,MPI_MINLOC,MPI_COMM_OPA,ierror)
@@ -2331,6 +2300,6 @@
       !!-----------------------------------------------------------------------
       !
-      zmax  = MAXVAL( ptab(:,:) , mask= pmask == 1.e0 )
-      ilocs = MAXLOC( ptab(:,:) , mask= pmask == 1.e0 )
+      zmax  = MAXVAL( ptab(:,:) , mask= pmask == 1._wp )
+      ilocs = MAXLOC( ptab(:,:) , mask= pmask == 1._wp )
       !
       ki = ilocs(1) + nimpp - 1
@@ -2359,17 +2328,17 @@
       !!
       !!--------------------------------------------------------------------------
-      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   ptab         ! Local 2D array
-      REAL(wp), DIMENSION (jpi,jpj,jpk), INTENT(in   ) ::   pmask        ! Local mask
-      REAL(wp)                         , INTENT(  out) ::   pmax         ! Global maximum of ptab
-      INTEGER                          , INTENT(  out) ::   ki, kj, kk   ! index of maximum in global frame
-      !!
-      REAL(wp) :: zmax   ! local maximum
+      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   ptab         ! Local 2D array
+      REAL(wp), DIMENSION (:,:,:), INTENT(in   ) ::   pmask        ! Local mask
+      REAL(wp)                   , INTENT(  out) ::   pmax         ! Global maximum of ptab
+      INTEGER                    , INTENT(  out) ::   ki, kj, kk   ! index of maximum in global frame
+      !
+      INTEGER  ::   ierror   ! local integer
+      REAL(wp) ::   zmax     ! local maximum
       REAL(wp), DIMENSION(2,1) ::   zain, zaout
       INTEGER , DIMENSION(3)   ::   ilocs
-      INTEGER :: ierror
       !!-----------------------------------------------------------------------
       !
-      zmax  = MAXVAL( ptab(:,:,:) , mask= pmask == 1.e0 )
-      ilocs = MAXLOC( ptab(:,:,:) , mask= pmask == 1.e0 )
+      zmax  = MAXVAL( ptab(:,:,:) , mask= pmask == 1._wp )
+      ilocs = MAXLOC( ptab(:,:,:) , mask= pmask == 1._wp )
       !
       ki = ilocs(1) + nimpp - 1
@@ -2377,6 +2346,6 @@
       kk = ilocs(3)
       !
-      zain(1,:)=zmax
-      zain(2,:)=ki+10000.*kj+100000000.*kk
+      zain(1,:) = zmax
+      zain(2,:) = ki + 10000.*kj + 100000000.*kk
       !
       CALL MPI_ALLREDUCE( zain,zaout, 1, MPI_2DOUBLE_PRECISION,MPI_MAXLOC,MPI_COMM_OPA,ierror)
@@ -2422,5 +2391,4 @@
 
    SUBROUTINE mpp_comm_free( kcom )
-      !!----------------------------------------------------------------------
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kcom
@@ -2692,12 +2660,11 @@
       !!              and apply lbc north-fold on this sub array. Then we
       !!              scatter the north fold array back to the processors.
-      !!
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   pt3d      ! 3D array on which the b.c. is applied
-      CHARACTER(len=1)                , INTENT(in   ) ::   cd_type   ! nature of pt3d grid-points
-      !                                                              !   = T ,  U , V , F or W  gridpoints
-      REAL(wp)                        , INTENT(in   ) ::   psgn      ! = -1. the sign change across the north fold 
-      !!                                                             ! =  1. , the sign is kept
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pt3d      ! 3D array on which the b.c. is applied
+      CHARACTER(len=1)          , INTENT(in   ) ::   cd_type   ! nature of pt3d grid-points
+      REAL(wp)                  , INTENT(in   ) ::   psgn      ! sign used across the north fold
+      !
       INTEGER ::   ji, jj, jr, jk
+      INTEGER ::   ipk                  ! 3rd dimension of the input array
       INTEGER ::   ierr, itaille, ildi, ilei, iilb
       INTEGER ::   ijpj, ijpjm1, ij, iproc
@@ -2715,12 +2682,14 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( ztab(jpiglo,4,jpk) , znorthloc(jpi,4,jpk), zfoldwk(jpi,4,jpk), znorthgloio(jpi,4,jpk,jpni) )
-      ALLOCATE( ztabl(jpi,4,jpk), ztabr(jpi*jpmaxngh, 4, jpk) ) 
+      ipk = SIZE( pt3d, 3 )
+      !
+      ALLOCATE( ztab (jpiglo,4,ipk), znorthloc(jpi,4,ipk), zfoldwk(jpi,4,ipk), znorthgloio(jpi,4,ipk,jpni) )
+      ALLOCATE( ztabl(jpi   ,4,ipk), ztabr(jpi*jpmaxngh,4,ipk)   ) 
 
       ijpj   = 4
       ijpjm1 = 3
       !
-      znorthloc(:,:,:) = 0
-      DO jk = 1, jpk
+      znorthloc(:,:,:) = 0._wp
+      DO jk = 1, ipk
          DO jj = nlcj - ijpj +1, nlcj          ! put in xnorthloc the last 4 jlines of pt3d
             ij = jj - nlcj + ijpj
@@ -2730,12 +2699,12 @@
       !
       !                                     ! Build in procs of ncomm_north the znorthgloio
-      itaille = jpi * jpk * ijpj
+      itaille = jpi * ipk * ijpj
 
       IF ( l_north_nogather ) THEN
          !
-        ztabr(:,:,:) = 0
-        ztabl(:,:,:) = 0
-
-        DO jk = 1, jpk
+        ztabr(:,:,:) = 0._wp
+        ztabl(:,:,:) = 0._wp
+
+        DO jk = 1, ipk
            DO jj = nlcj-ijpj+1, nlcj          ! First put local values into the global array
               ij = jj - nlcj + ijpj
@@ -2747,5 +2716,5 @@
 
          DO jr = 1,nsndto
-            IF ((nfipproc(isendto(jr),jpnj) .ne. (narea-1)) .and. (nfipproc(isendto(jr),jpnj) .ne. -1)) THEN
+            IF ((nfipproc(isendto(jr),jpnj) /= (narea-1)) .and. (nfipproc(isendto(jr),jpnj) /= -1)) THEN
               CALL mppsend( 5, znorthloc, itaille, nfipproc(isendto(jr),jpnj), ml_req_nf(jr) )
             ENDIF
@@ -2753,12 +2722,12 @@
          DO jr = 1,nsndto
             iproc = nfipproc(isendto(jr),jpnj)
-            IF(iproc .ne. -1) THEN
+            IF(iproc /= -1) THEN
                ilei = nleit (iproc+1)
                ildi = nldit (iproc+1)
                iilb = nfiimpp(isendto(jr),jpnj) - nfiimpp(isendto(1),jpnj)
             ENDIF
-            IF((iproc .ne. (narea-1)) .and. (iproc .ne. -1)) THEN
+            IF((iproc /= (narea-1)) .and. (iproc /= -1)) THEN
               CALL mpprecv(5, zfoldwk, itaille, iproc)
-              DO jk = 1, jpk
+              DO jk = 1, ipk
                  DO jj = 1, ijpj
                     DO ji = ildi, ilei
@@ -2767,6 +2736,6 @@
                  END DO
               END DO
-           ELSE IF (iproc .eq. (narea-1)) THEN
-              DO jk = 1, jpk
+           ELSE IF( iproc == narea-1 ) THEN
+              DO jk = 1, ipk
                  DO jj = 1, ijpj
                     DO ji = ildi, ilei
@@ -2779,5 +2748,5 @@
          IF (l_isend) THEN
             DO jr = 1,nsndto
-               IF ((nfipproc(isendto(jr),jpnj) .ne. (narea-1)) .and. (nfipproc(isendto(jr),jpnj) .ne. -1)) THEN
+               IF ((nfipproc(isendto(jr),jpnj) /= (narea-1)) .and. (nfipproc(isendto(jr),jpnj) /= -1)) THEN
                   CALL mpi_wait(ml_req_nf(jr), ml_stat, ml_err)
                ENDIF    
@@ -2785,5 +2754,5 @@
          ENDIF
          CALL mpp_lbc_nfd( ztabl, ztabr, cd_type, psgn )   ! North fold boundary condition
-         DO jk = 1, jpk
+         DO jk = 1, ipk
             DO jj = nlcj-ijpj+1, nlcj             ! Scatter back to pt3d
                ij = jj - nlcj + ijpj
@@ -2794,10 +2763,9 @@
          END DO
          !
-
       ELSE
          CALL MPI_ALLGATHER( znorthloc  , itaille, MPI_DOUBLE_PRECISION,                &
             &                znorthgloio, itaille, MPI_DOUBLE_PRECISION, ncomm_north, ierr )
          !
-         ztab(:,:,:) = 0.e0
+         ztab(:,:,:) = 0._wp
          DO jr = 1, ndim_rank_north         ! recover the global north array
             iproc = nrank_north(jr) + 1
@@ -2805,5 +2773,5 @@
             ilei  = nleit (iproc)
             iilb  = nimppt(iproc)
-            DO jk = 1, jpk
+            DO jk = 1, ipk
                DO jj = 1, ijpj
                   DO ji = ildi, ilei
@@ -2815,5 +2783,5 @@
          CALL lbc_nfd( ztab, cd_type, psgn )   ! North fold boundary condition
          !
-         DO jk = 1, jpk
+         DO jk = 1, ipk
             DO jj = nlcj-ijpj+1, nlcj             ! Scatter back to pt3d
                ij = jj - nlcj + ijpj
@@ -2902,5 +2870,5 @@
 
          DO jr = 1,nsndto
-            IF ((nfipproc(isendto(jr),jpnj) .ne. (narea-1)) .and. (nfipproc(isendto(jr),jpnj) .ne. -1)) THEN
+            IF ((nfipproc(isendto(jr),jpnj) /= (narea-1)) .and. (nfipproc(isendto(jr),jpnj) /= -1)) THEN
                CALL mppsend(5, znorthloc, itaille, nfipproc(isendto(jr),jpnj), ml_req_nf(jr))
             ENDIF
@@ -2908,10 +2876,10 @@
          DO jr = 1,nsndto
             iproc = nfipproc(isendto(jr),jpnj)
-            IF(iproc .ne. -1) THEN
+            IF(iproc /= -1) THEN
                ilei = nleit (iproc+1)
                ildi = nldit (iproc+1)
                iilb = nfiimpp(isendto(jr),jpnj) - nfiimpp(isendto(1),jpnj)
             ENDIF
-            IF((iproc .ne. (narea-1)) .and. (iproc .ne. -1)) THEN
+            IF( iproc /= narea-1 .AND. iproc /= -1 ) THEN
               CALL mpprecv(5, zfoldwk, itaille, iproc)
               DO jj = 1, ijpj
@@ -2920,5 +2888,5 @@
                  END DO
               END DO
-            ELSE IF (iproc .eq. (narea-1)) THEN
+            ELSE IF( iproc == narea-1 ) THEN
               DO jj = 1, ijpj
                  DO ji = ildi, ilei
@@ -2930,5 +2898,5 @@
          IF (l_isend) THEN
             DO jr = 1,nsndto
-               IF ((nfipproc(isendto(jr),jpnj) .ne. (narea-1)) .and. (nfipproc(isendto(jr),jpnj) .ne. -1)) THEN
+               IF ((nfipproc(isendto(jr),jpnj) /= (narea-1)) .and. (nfipproc(isendto(jr),jpnj) /= -1)) THEN
                   CALL mpi_wait(ml_req_nf(jr), ml_stat, ml_err)
                ENDIF
@@ -2948,5 +2916,5 @@
             &                znorthgloio, itaille, MPI_DOUBLE_PRECISION, ncomm_north, ierr )
          !
-         ztab(:,:) = 0.e0
+         ztab(:,:) = 0._wp
          DO jr = 1, ndim_rank_north            ! recover the global north array
             iproc = nrank_north(jr) + 1
@@ -2975,5 +2943,6 @@
    END SUBROUTINE mpp_lbc_north_2d
 
-   SUBROUTINE mpp_lbc_north_2d_multiple( pt2d_array, cd_type, psgn, num_fields)
+
+   SUBROUTINE mpp_lbc_north_2d_multiple( pt2d_array, cd_type, psgn, kfld )
       !!---------------------------------------------------------------------
       !!                   ***  routine mpp_lbc_north_2d  ***
@@ -2990,28 +2959,26 @@
       !!
       !!----------------------------------------------------------------------
-      INTEGER ,  INTENT (in   ) ::   num_fields  ! number of variables contained in pt2d
-      TYPE( arrayptr ), DIMENSION(:) :: pt2d_array
-      CHARACTER(len=1), DIMENSION(:), INTENT(in   ) ::   cd_type   ! nature of pt2d grid-points
-      !                                                          !   = T ,  U , V , F or W  gridpoints
-      REAL(wp), DIMENSION(:), INTENT(in   ) ::   psgn      ! = -1. the sign change across the north fold 
-      !!                                                             ! =  1. , the sign is kept
+      TYPE( arrayptr ), DIMENSION(:), INTENT(inout) ::   pt2d_array   ! pointer array of 2D fields
+      CHARACTER(len=1), DIMENSION(:), INTENT(in   ) ::   cd_type      ! nature of pt2d grid-points
+      REAL(wp)        , DIMENSION(:), INTENT(in   ) ::   psgn         ! sign used across the north fold 
+      INTEGER                       , INTENT(in   ) ::   kfld         ! number of variables contained in pt2d
+      !
       INTEGER ::   ji, jj, jr, jk
       INTEGER ::   ierr, itaille, ildi, ilei, iilb
-      INTEGER ::   ijpj, ijpjm1, ij, iproc
-      INTEGER, DIMENSION (jpmaxngh)      ::   ml_req_nf          !for mpi_isend when avoiding mpi_allgather
-      INTEGER                            ::   ml_err             ! for mpi_isend when avoiding mpi_allgather
-      INTEGER, DIMENSION(MPI_STATUS_SIZE)::   ml_stat            ! for mpi_isend when avoiding mpi_allgather
-      !                                                              ! Workspace for message transfers avoiding mpi_allgather
+      INTEGER ::   ijpj, ijpjm1, ij, iproc, iflag
+      INTEGER, DIMENSION (jpmaxngh)      ::   ml_req_nf   ! for mpi_isend when avoiding mpi_allgather
+      INTEGER                            ::   ml_err      ! for mpi_isend when avoiding mpi_allgather
+      INTEGER, DIMENSION(MPI_STATUS_SIZE)::   ml_stat     ! for mpi_isend when avoiding mpi_allgather
+      !                                                   ! Workspace for message transfers avoiding mpi_allgather
+      INTEGER :: istatus(mpi_status_size)
       REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE   :: ztab
       REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE   :: znorthloc, zfoldwk
       REAL(wp), DIMENSION(:,:,:,:), ALLOCATABLE   :: znorthgloio
       REAL(wp), DIMENSION(:,:,:)  , ALLOCATABLE   :: ztabl, ztabr
-      INTEGER :: istatus(mpi_status_size)
-      INTEGER :: iflag
-      !!----------------------------------------------------------------------
-      !
-      ALLOCATE( ztab(jpiglo,4,num_fields), znorthloc(jpi,4,num_fields), zfoldwk(jpi,4,num_fields),   &
-            &   znorthgloio(jpi,4,num_fields,jpni) )   ! expanded to 3 dimensions
-      ALLOCATE( ztabl(jpi,4,num_fields), ztabr(jpi*jpmaxngh, 4,num_fields) )
+      !!----------------------------------------------------------------------
+      !
+      ALLOCATE( ztab(jpiglo,4,kfld), znorthloc  (jpi,4,kfld),        &
+         &      zfoldwk(jpi,4,kfld), znorthgloio(jpi,4,kfld,jpni),   &
+         &      ztabl  (jpi,4,kfld), ztabr(jpi*jpmaxngh, 4,kfld)   )
       !
       ijpj   = 4
@@ -3019,5 +2986,5 @@
       !
       
-      DO jk = 1, num_fields
+      DO jk = 1, kfld
          DO jj = nlcj-ijpj+1, nlcj             ! put in znorthloc the last 4 jlines of pt2d (for every variable)
             ij = jj - nlcj + ijpj
@@ -3033,8 +3000,8 @@
          ! (in nemo_northcomms) as being  involved in this process' northern boundary exchange
          !
-         ztabr(:,:,:) = 0
-         ztabl(:,:,:) = 0
-
-         DO jk = 1, num_fields
+         ztabr(:,:,:) = 0._wp
+         ztabl(:,:,:) = 0._wp
+
+         DO jk = 1, kfld
             DO jj = nlcj-ijpj+1, nlcj          ! First put local values into the global array
                ij = jj - nlcj + ijpj
@@ -3045,19 +3012,19 @@
          END DO
 
-         DO jr = 1,nsndto
-            IF ((nfipproc(isendto(jr),jpnj) .ne. (narea-1)) .and. (nfipproc(isendto(jr),jpnj) .ne. -1)) THEN
-               CALL mppsend(5, znorthloc, itaille*num_fields, nfipproc(isendto(jr),jpnj), ml_req_nf(jr)) ! Buffer expanded "num_fields" times
+         DO jr = 1, nsndto
+            IF( nfipproc(isendto(jr),jpnj) /= narea-1 .AND. nfipproc(isendto(jr),jpnj) /= -1 ) THEN
+               CALL mppsend(5, znorthloc, itaille*kfld, nfipproc(isendto(jr),jpnj), ml_req_nf(jr)) ! Buffer expanded "kfld" times
             ENDIF
          END DO
-         DO jr = 1,nsndto
+         DO jr = 1, nsndto
             iproc = nfipproc(isendto(jr),jpnj)
-            IF(iproc .ne. -1) THEN
+            IF(iproc /= -1) THEN
                ilei = nleit (iproc+1)
                ildi = nldit (iproc+1)
                iilb = nfiimpp(isendto(jr),jpnj) - nfiimpp(isendto(1),jpnj)
             ENDIF
-            IF((iproc .ne. (narea-1)) .and. (iproc .ne. -1)) THEN
-              CALL mpprecv(5, zfoldwk, itaille*num_fields, iproc) ! Buffer expanded "num_fields" times
-              DO jk = 1 , num_fields
+            IF( iproc /= narea-1 .AND. iproc /= -1 ) THEN
+              CALL mpprecv(5, zfoldwk, itaille*kfld, iproc) ! Buffer expanded "kfld" times
+              DO jk = 1 , kfld
                  DO jj = 1, ijpj
                     DO ji = ildi, ilei
@@ -3066,6 +3033,6 @@
                  END DO
               END DO
-            ELSE IF (iproc .eq. (narea-1)) THEN
-              DO jk = 1, num_fields
+            ELSEIF ( iproc == narea-1 ) THEN
+              DO jk = 1, kfld
                  DO jj = 1, ijpj
                     DO ji = ildi, ilei
@@ -3076,7 +3043,7 @@
             ENDIF
          END DO
-         IF (l_isend) THEN
-            DO jr = 1,nsndto
-               IF ((nfipproc(isendto(jr),jpnj) .ne. (narea-1)) .and. (nfipproc(isendto(jr),jpnj) .ne. -1)) THEN
+         IF( l_isend ) THEN
+            DO jr = 1, nsndto
+               IF( nfipproc(isendto(jr),jpnj) /= narea-1 .AND. nfipproc(isendto(jr),jpnj) /= -1 ) THEN
                   CALL mpi_wait(ml_req_nf(jr), ml_stat, ml_err)
                ENDIF
@@ -3084,9 +3051,9 @@
          ENDIF
          !
-         DO ji = 1, num_fields     ! Loop to manage 3D variables
+         DO ji = 1, kfld     ! Loop to manage 3D variables
             CALL mpp_lbc_nfd( ztabl(:,:,ji), ztabr(:,:,ji), cd_type(ji), psgn(ji) )  ! North fold boundary condition
          END DO
          !
-         DO jk = 1, num_fields
+         DO jk = 1, kfld
             DO jj = nlcj-ijpj+1, nlcj             ! Scatter back to pt2d
                ij = jj - nlcj + ijpj
@@ -3100,9 +3067,9 @@
       ELSE
          !
-         CALL MPI_ALLGATHER( znorthloc  , itaille*num_fields, MPI_DOUBLE_PRECISION,        &
-            &                znorthgloio, itaille*num_fields, MPI_DOUBLE_PRECISION, ncomm_north, ierr )
-         !
-         ztab(:,:,:) = 0.e0
-         DO jk = 1, num_fields
+         CALL MPI_ALLGATHER( znorthloc  , itaille*kfld, MPI_DOUBLE_PRECISION,        &
+            &                znorthgloio, itaille*kfld, MPI_DOUBLE_PRECISION, ncomm_north, ierr )
+         !
+         ztab(:,:,:) = 0._wp
+         DO jk = 1, kfld
             DO jr = 1, ndim_rank_north            ! recover the global north array
                iproc = nrank_north(jr) + 1
@@ -3118,9 +3085,9 @@
          END DO
          
-         DO ji = 1, num_fields
+         DO ji = 1, kfld
             CALL lbc_nfd( ztab(:,:,ji), cd_type(ji), psgn(ji) )   ! North fold boundary condition
          END DO
          !
-         DO jk = 1, num_fields
+         DO jk = 1, kfld
             DO jj = nlcj-ijpj+1, nlcj             ! Scatter back to pt2d
                ij = jj - nlcj + ijpj
@@ -3138,4 +3105,5 @@
    END SUBROUTINE mpp_lbc_north_2d_multiple
 
+
    SUBROUTINE mpp_lbc_north_e( pt2d, cd_type, psgn)
       !!---------------------------------------------------------------------
@@ -3155,40 +3123,36 @@
       REAL(wp), DIMENSION(1-jpr2di:jpi+jpr2di,1-jpr2dj:jpj+jpr2dj), INTENT(inout) ::   pt2d     ! 2D array with extra halo
       CHARACTER(len=1)                                            , INTENT(in   ) ::   cd_type  ! nature of pt3d grid-points
-      !                                                                                         !   = T ,  U , V , F or W -points
-      REAL(wp)                                                    , INTENT(in   ) ::   psgn     ! = -1. the sign change across the
-      !!                                                                                        ! north fold, =  1. otherwise
+      REAL(wp)                                                    , INTENT(in   ) ::   psgn     ! sign used across the north fold
+      !
       INTEGER ::   ji, jj, jr
       INTEGER ::   ierr, itaille, ildi, ilei, iilb
       INTEGER ::   ijpj, ij, iproc
-      !
       REAL(wp), DIMENSION(:,:)  , ALLOCATABLE  ::  ztab_e, znorthloc_e
       REAL(wp), DIMENSION(:,:,:), ALLOCATABLE  ::  znorthgloio_e
-
       !!----------------------------------------------------------------------
       !
       ALLOCATE( ztab_e(jpiglo,4+2*jpr2dj), znorthloc_e(jpi,4+2*jpr2dj), znorthgloio_e(jpi,4+2*jpr2dj,jpni) )
-
       !
       ijpj=4
-      ztab_e(:,:) = 0.e0
-
-      ij=0
+      ztab_e(:,:) = 0._wp
+
+      ij = 0
       ! put in znorthloc_e the last 4 jlines of pt2d
       DO jj = nlcj - ijpj + 1 - jpr2dj, nlcj +jpr2dj
          ij = ij + 1
          DO ji = 1, jpi
-            znorthloc_e(ji,ij)=pt2d(ji,jj)
+            znorthloc_e(ji,ij) = pt2d(ji,jj)
          END DO
       END DO
       !
       itaille = jpi * ( ijpj + 2 * jpr2dj )
-      CALL MPI_ALLGATHER( znorthloc_e(1,1)  , itaille, MPI_DOUBLE_PRECISION,    &
+      CALL MPI_ALLGATHER( znorthloc_e(1,1)    , itaille, MPI_DOUBLE_PRECISION,    &
          &                znorthgloio_e(1,1,1), itaille, MPI_DOUBLE_PRECISION, ncomm_north, ierr )
       !
       DO jr = 1, ndim_rank_north            ! recover the global north array
          iproc = nrank_north(jr) + 1
-         ildi = nldit (iproc)
-         ilei = nleit (iproc)
-         iilb = nimppt(iproc)
+         ildi  = nldit (iproc)
+         ilei  = nleit (iproc)
+         iilb  = nimppt(iproc)
          DO jj = 1, ijpj+2*jpr2dj
             DO ji = ildi, ilei
@@ -3197,5 +3161,4 @@
          END DO
       END DO
-
 
       ! 2. North-Fold boundary conditions
@@ -3238,12 +3201,11 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   ptab     ! 3D array on which the boundary condition is applied
-      CHARACTER(len=1)                , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points
-      !                                                             ! = T , U , V , F , W points
-      REAL(wp)                        , INTENT(in   ) ::   psgn     ! =-1 the sign change across the north fold boundary
-      !                                                             ! =  1. , the sign is kept
-      INTEGER                         , INTENT(in   ) ::   ib_bdy   ! BDY boundary set
+      REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   ptab     ! 3D array on which the boundary condition is applied
+      CHARACTER(len=1)          , INTENT(in   ) ::   cd_type  ! nature of ptab grid point
+      REAL(wp)                  , INTENT(in   ) ::   psgn     ! sign used across the north fold boundary
+      INTEGER                   , INTENT(in   ) ::   ib_bdy   ! BDY boundary set
       !
       INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
+      INTEGER  ::   ipk                        ! 3rd dimension of the input array
       INTEGER  ::   imigr, iihom, ijhom        ! local integers
       INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
@@ -3255,6 +3217,8 @@
       !!----------------------------------------------------------------------
       !
-      ALLOCATE( zt3ns(jpi,jprecj,jpk,2), zt3sn(jpi,jprecj,jpk,2),   &
-         &      zt3ew(jpj,jpreci,jpk,2), zt3we(jpj,jpreci,jpk,2)  )
+      ipk = SIZE( ptab, 3 )
+      !      
+      ALLOCATE( zt3ns(jpi,jprecj,ipk,2), zt3sn(jpi,jprecj,ipk,2),   &
+         &      zt3ew(jpj,jpreci,ipk,2), zt3we(jpj,jpreci,ipk,2)  )
 
       zland = 0._wp
@@ -3263,5 +3227,5 @@
       ! ------------------------------
       !                                   ! East-West boundaries
-      !                                        !* Cyclic east-west
+      !                                        !* Cyclic
       IF( nbondi == 2) THEN
          IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
@@ -3273,5 +3237,5 @@
          ENDIF
       ELSEIF(nbondi == -1) THEN
-         IF( .NOT. cd_type == 'F' )   ptab(1:jpreci,:,:) = zland    ! south except F-point
+         IF( .NOT. cd_type == 'F' )   ptab(1:jpreci,:,:) = zland       ! south except F-point
       ELSEIF(nbondi == 1) THEN
          ptab(nlci-jpreci+1:jpi,:,:) = zland    ! north
@@ -3298,5 +3262,5 @@
       !
       !                           ! Migrations
-      imigr = jpreci * jpj * jpk
+      imigr = jpreci * jpj * ipk
       !
       SELECT CASE ( nbondi_bdy(ib_bdy) )
@@ -3348,5 +3312,4 @@
          END DO
       END SELECT
-
 
       ! 3. North and south directions
@@ -3363,5 +3326,5 @@
       !
       !                           ! Migrations
-      imigr = jprecj * jpi * jpk
+      imigr = jprecj * jpi * ipk
       !
       SELECT CASE ( nbondj_bdy(ib_bdy) )
@@ -3413,5 +3376,4 @@
          END DO
       END SELECT
-
 
       ! 4. north fold treatment
@@ -3453,12 +3415,10 @@
       !!
       !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   ptab     ! 3D array on which the boundary condition is applied
-      CHARACTER(len=1)            , INTENT(in   ) ::   cd_type  ! define the nature of ptab array grid-points
-      !                                                         ! = T , U , V , F , W points
-      REAL(wp)                    , INTENT(in   ) ::   psgn     ! =-1 the sign change across the north fold boundary
-      !                                                         ! =  1. , the sign is kept
-      INTEGER                     , INTENT(in   ) ::   ib_bdy   ! BDY boundary set
-      !
-      INTEGER  ::   ji, jj, jl             ! dummy loop indices
+      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(inout) ::   ptab     ! 3D array on which the boundary condition is applied
+      CHARACTER(len=1)                , INTENT(in   ) ::   cd_type  ! nature of ptab array grid-points
+      REAL(wp)                        , INTENT(in   ) ::   psgn     ! sign used across the north fold boundary
+      INTEGER                         , INTENT(in   ) ::   ib_bdy   ! BDY boundary set
+      !
+      INTEGER  ::   ji, jj, jl                 ! dummy loop indices
       INTEGER  ::   imigr, iihom, ijhom        ! local integers
       INTEGER  ::   ml_req1, ml_req2, ml_err   ! for key_mpi_isend
@@ -3478,17 +3438,17 @@
       ! ------------------------------
       !                                   ! East-West boundaries
-      !                                      !* Cyclic east-west
+      !                                         !* Cyclic
       IF( nbondi == 2 ) THEN
-         IF (nperio == 1 .OR. nperio == 4 .OR. nperio == 6) THEN
+         IF( nperio == 1 .OR. nperio == 4 .OR. nperio == 6 ) THEN
             ptab( 1 ,:) = ptab(jpim1,:)
             ptab(jpi,:) = ptab(  2  ,:)
          ELSE
-            IF(.NOT.cd_type == 'F' )  ptab(     1       :jpreci,:) = zland    ! south except F-point
-                                      ptab(nlci-jpreci+1:jpi   ,:) = zland    ! north
+            IF(.NOT.cd_type == 'F' )   ptab(     1       :jpreci,:) = zland    ! south except F-point
+                                       ptab(nlci-jpreci+1:jpi   ,:) = zland    ! north
          ENDIF
       ELSEIF(nbondi == -1) THEN
-         IF( .NOT.cd_type == 'F' )    ptab(     1       :jpreci,:) = zland    ! south except F-point
+         IF(.NOT.cd_type == 'F' )      ptab(     1       :jpreci,:) = zland    ! south except F-point
       ELSEIF(nbondi == 1) THEN
-                                      ptab(nlci-jpreci+1:jpi   ,:) = zland    ! north
+                                       ptab(nlci-jpreci+1:jpi   ,:) = zland    ! north
       ENDIF
       !                                      !* closed
@@ -3537,10 +3497,10 @@
       SELECT CASE ( nbondi_bdy(ib_bdy) )
       CASE ( -1 )
-         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
-      CASE ( 0 )
-         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
-         IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err)
-      CASE ( 1 )
-         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err)
+         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err )
+      CASE ( 0 )
+         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err )
+         IF(l_isend) CALL mpi_wait(ml_req2, ml_stat, ml_err )
+      CASE ( 1 )
+         IF(l_isend) CALL mpi_wait(ml_req1, ml_stat, ml_err )
       END SELECT
       !
@@ -3628,5 +3588,4 @@
          END DO
       END SELECT
-
 
       ! 4. north fold treatment
@@ -3713,10 +3672,11 @@
       !!   This subroutine computes yddb(i) = ydda(i)+yddb(i)
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in)                         :: ilen, itype
-      COMPLEX(wp), DIMENSION(ilen), INTENT(in)     :: ydda
-      COMPLEX(wp), DIMENSION(ilen), INTENT(inout)  :: yddb
+      INTEGER                     , INTENT(in)    ::   ilen, itype
+      COMPLEX(wp), DIMENSION(ilen), INTENT(in)    ::   ydda
+      COMPLEX(wp), DIMENSION(ilen), INTENT(inout) ::   yddb
       !
       REAL(wp) :: zerr, zt1, zt2    ! local work variables
-      INTEGER :: ji, ztmp           ! local scalar
+      INTEGER  :: ji, ztmp           ! local scalar
+      !!---------------------------------------------------------------------
 
       ztmp = itype   ! avoid compilation warning
@@ -3841,11 +3801,10 @@
       !!                    nono   : number for local neighboring processors
       !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(1-jpri:jpi+jpri,1-jprj:jpj+jprj), INTENT(inout) ::   pt2d     ! 2D array with extra halo
+      CHARACTER(len=1)                                    , INTENT(in   ) ::   cd_type  ! nature of ptab array grid-points
+      REAL(wp)                                            , INTENT(in   ) ::   psgn     ! sign used across the north fold
       INTEGER                                             , INTENT(in   ) ::   jpri
       INTEGER                                             , INTENT(in   ) ::   jprj
-      REAL(wp), DIMENSION(1-jpri:jpi+jpri,1-jprj:jpj+jprj), INTENT(inout) ::   pt2d     ! 2D array with extra halo
-      CHARACTER(len=1)                                    , INTENT(in   ) ::   cd_type  ! nature of ptab array grid-points
-      !                                                                                 ! = T , U , V , F , W and I points
-      REAL(wp)                                            , INTENT(in   ) ::   psgn     ! =-1 the sign change across the
-      !!                                                                                ! north boundary, =  1. otherwise
+      !
       INTEGER  ::   jl   ! dummy loop indices
       INTEGER  ::   imigr, iihom, ijhom        ! temporary integers
@@ -3875,6 +3834,6 @@
          !
       ELSE                                        !* closed
-         IF( .NOT. cd_type == 'F' )   pt2d(  1-jpri   :jpreci    ,:) = 0.e0    ! south except at F-point
-                                      pt2d(nlci-jpreci+1:jpi+jpri,:) = 0.e0    ! north
+         IF( .NOT. cd_type == 'F' )   pt2d(  1-jpri   :jpreci    ,:) = 0._wp    ! south except at F-point
+                                      pt2d(nlci-jpreci+1:jpi+jpri,:) = 0._wp    ! north
       ENDIF
       !
@@ -3996,5 +3955,5 @@
          END DO
       END SELECT
-
+      !
    END SUBROUTINE mpp_lnk_2d_icb
    
@@ -4020,4 +3979,7 @@
       MODULE PROCEDURE mpp_maxloc2d ,mpp_maxloc3d
    END INTERFACE
+   INTERFACE mpp_max_multiple
+      MODULE PROCEDURE mppmax_real_multiple
+   END INTERFACE
 
    LOGICAL, PUBLIC, PARAMETER ::   lk_mpp = .FALSE.      !: mpp flag
@@ -4191,4 +4153,12 @@
       WRITE(*,*) 'mpp_comm_free: You should not have seen this print! error?', kcom
    END SUBROUTINE mpp_comm_free
+   
+   SUBROUTINE mppmax_real_multiple( ptab, kdim , kcom  )
+      REAL, DIMENSION(:) ::   ptab   ! 
+      INTEGER            ::   kdim   ! 
+      INTEGER, OPTIONAL  ::   kcom   ! 
+      WRITE(*,*) 'mppmax_real_multiple: You should not have seen this print! error?', ptab(1), kdim
+   END SUBROUTINE mppmax_real_multiple
+
 #endif
 
@@ -4225,5 +4195,5 @@
                                CALL FLUSH(numout    )
       IF( numstp     /= -1 )   CALL FLUSH(numstp    )
-      IF( numsol     /= -1 )   CALL FLUSH(numsol    )
+      IF( numrun     /= -1 )   CALL FLUSH(numrun    )
       IF( numevo_ice /= -1 )   CALL FLUSH(numevo_ice)
       !
@@ -4332,5 +4302,5 @@
             WRITE(kout,*)
          ENDIF
-         CALL FLUSH(kout) 
+         CALL FLUSH( kout ) 
          STOP 'ctl_opn bad opening'
       ENDIF
