Index: /branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
===================================================================
--- /branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 5578)
+++ /branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 5579)
@@ -28,4 +28,5 @@
 
    PUBLIC   lim_hdf         ! called by lim_trp
+   PUBLIC   lim_hdf_multiple ! called by lim_trp
    PUBLIC   lim_hdf_init    ! called by sbc_lim_init
 
@@ -124,5 +125,5 @@
          CALL lbc_lnk( zrlx, 'T', 1. )                   ! lateral boundary condition
          !
-         IF ( MOD( iter, nn_convfrq ) == 0 )  THEN    ! convergence test every nn_convfrq iterations (perf. optimization)
+         IF ( MOD( iter - 1 , nn_convfrq ) == 0 )  THEN    ! convergence test every nn_convfrq iterations (perf. optimization)
             zconv = 0._wp
             DO jj = 2, jpjm1
@@ -166,4 +167,198 @@
       !
    END SUBROUTINE lim_hdf
+   
+
+   SUBROUTINE lim_hdf_multiple( ptab , ihdf_vars , jpl , nlay_i )
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE lim_hdf  ***
+      !!
+      !! ** purpose :   Compute and add the diffusive trend on sea-ice variables
+      !!
+      !! ** method  :   Second order diffusive operator evaluated using a
+      !!              Cranck-Nicholson time Scheme.
+      !!
+      !! ** Action  :    update ptab with the diffusive contribution
+      !!-------------------------------------------------------------------
+      INTEGER                           :: jpl, nlay_i, isize, ihdf_vars
+      REAL(wp),  DIMENSION(:,:,:), INTENT( inout ),TARGET ::   ptab    ! Field on which the diffusion is applied
+      REAL(wp), POINTER, DIMENSION(:,:,:)        ::   pahu3D , pahv3D
+      !
+      INTEGER                           ::  ji, jj, jk, jl , jm               ! dummy loop indices
+      INTEGER                           ::  iter, ierr           ! local integers
+      REAL(wp)                          ::  zrlxint     ! local scalars
+      REAL(wp), POINTER , DIMENSION ( : )        :: zconv     ! local scalars
+      REAL(wp), POINTER , DIMENSION(:,:,:) ::  zrlx,zdiv0, ztab0
+      REAL(wp), POINTER , DIMENSION(:,:) ::  zflu, zflv, zdiv
+      CHARACTER(lc)                     ::  charout                   ! local character
+      REAL(wp), PARAMETER               ::  zrelax = 0.5_wp           ! relaxation constant for iterative procedure
+      REAL(wp), PARAMETER               ::  zalfa  = 0.5_wp           ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit
+      INTEGER , PARAMETER               ::  its    = 100              ! Maximum number of iteration
+      !!-------------------------------------------------------------------
+      TYPE(arrayptr)   , ALLOCATABLE, DIMENSION(:) ::   pt2d_array, zrlx_array
+      CHARACTER(len=1) , ALLOCATABLE, DIMENSION(:) ::   type_array ! define the nature of ptab array grid-points
+      !                                                            ! = T , U , V , F , W and I points
+      REAL(wp)        , ALLOCATABLE, DIMENSION(:)  ::   psgn_array    ! =-1 the sign change across the north fold boundary
+
+     !!--------------------------------------------------------------------- 
+
+      !                       !==  Initialisation  ==!
+         isize = jpl*(ihdf_vars+nlay_i)
+      ALLOCATE( zconv (isize) )
+      ALLOCATE( pt2d_array(isize) , zrlx_array(isize) )
+      ALLOCATE( type_array(isize) )
+      ALLOCATE( psgn_array(isize) )
+      
+      CALL wrk_alloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 )
+      CALL wrk_alloc( jpi, jpj, zflu, zflv, zdiv )
+      CALL wrk_alloc( jpi, jpj, jpl, pahu3D , pahv3D )
+
+
+      DO jl = 1 , jpl
+         jm = (jl-1)*(ihdf_vars+nlay_i)+1
+         DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
+            DO ji = 1 , fs_jpim1   ! vector opt.
+               pahu3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ptab(ji  ,jj,jm) ) ) )   &
+               &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ptab(ji+1, jj, jm ) ) ) ) * ahiu(ji,jj)
+               pahv3D(ji,jj,jl) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -ptab(ji, jj, jm ) ) ) )   &
+               &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- ptab(ji, jj+1, jm ) ) ) ) * ahiv(ji,jj)
+            END DO
+         END DO
+      END DO
+
+      DO jk= 1 , isize
+         pt2d_array(jk)%pt2d=>ptab(:,:,jk)
+         zrlx_array(jk)%pt2d=>zrlx(:,:,jk)
+         type_array(jk)='T'
+         psgn_array(jk)=1.
+      END DO
+
+      !
+      IF( linit ) THEN              ! Metric coefficient (compute at the first call and saved in efact)
+         ALLOCATE( efact(jpi,jpj) , STAT=ierr )
+         IF( lk_mpp    )   CALL mpp_sum( ierr )
+         IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'lim_hdf : unable to allocate arrays' )
+         DO jj = 2, jpjm1
+            DO ji = fs_2 , fs_jpim1   ! vector opt.
+               efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) * r1_e12t(ji,jj)
+            END DO
+         END DO
+         linit = .FALSE.
+      ENDIF
+      !                             ! Time integration parameters
+      !
+      zflu (jpi,: ) = 0._wp
+      zflv (jpi,: ) = 0._wp
+
+      DO jk=1 , isize
+         ztab0(:, : , jk ) = ptab(:,:,jk)      ! Arrays initialization
+         zdiv0(:, 1 , jk ) = 0._wp
+         zdiv0(:,jpj, jk ) = 0._wp
+         zdiv0(1,  :, jk ) = 0._wp
+         zdiv0(jpi,:, jk ) = 0._wp
+      END DO
+
+      zconv = 1._wp           !==  horizontal diffusion using a Crant-Nicholson scheme  ==!
+      iter  = 0
+      !
+      DO WHILE( MAXVAL(zconv(:)) > ( 2._wp * 1.e-04 ) .AND. iter <= its )   ! Sub-time step loop
+         !
+         iter = iter + 1                                 ! incrementation of the sub-time step number
+         !
+
+         DO jk = 1 , isize
+            jl = (jk-1) /( ihdf_vars+nlay_i)+1
+            IF (zconv(jk) > ( 2._wp * 1.e-04 )) THEN
+               DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
+                  DO ji = 1 , fs_jpim1   ! vector opt.
+                     zflu(ji,jj) = pahu3D(ji,jj,jl) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj,jk) - ptab(ji,jj,jk) )
+                     zflv(ji,jj) = pahv3D(ji,jj,jl) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1,jk) - ptab(ji,jj,jk) )
+                  END DO
+               END DO
+               !
+               DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
+                  DO ji = fs_2 , fs_jpim1   ! vector opt. 
+                     zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
+                  END DO
+               END DO
+               !
+               IF( iter == 1 )   zdiv0(:,:,jk) = zdiv(:,:)        ! save the 1st evaluation of the diffusive trend in zdiv0
+               !
+               DO jj = 2, jpjm1                                ! iterative evaluation
+                  DO ji = fs_2 , fs_jpim1   ! vector opt.
+                     zrlxint = (   ztab0(ji,jj,jk)    &
+                        &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj,jk) )   &
+                        &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj,jk) )                               &
+                        &      ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )
+                     zrlx(ji,jj,jk) = ptab(ji,jj,jk) + zrelax * ( zrlxint - ptab(ji,jj,jk) )
+                  END DO
+               END DO
+            END IF
+
+         END DO
+
+         CALL lbc_lnk_multi( zrlx_array, type_array , psgn_array , isize ) ! Multiple interchange of all the variables
+         !
+         
+         IF ( MOD( iter-1 , nn_convfrq ) == 0 )  THEN   !Convergence test every nn_convfrq iterations (perf. optimization ) 
+            DO jk=1,isize
+               zconv(jk) = 0._wp                                   ! convergence test
+               DO jj = 2, jpjm1
+                  DO ji = fs_2, fs_jpim1
+                     zconv(jk) = MAX( zconv(jk), ABS( zrlx(ji,jj,jk) - ptab(ji,jj,jk) )  )
+                  END DO
+               END DO
+            END DO
+            IF( lk_mpp ) CALL mpp_max_multiple( zconv , isize )            ! max over the global domain for all the variables
+         ENDIF
+         !
+         DO jk=1,isize
+            ptab(:,:,jk) = zrlx(:,:,jk)
+         END DO
+         !
+      END DO                                       ! end of sub-time step loop
+
+     ! -----------------------
+      !!! final step (clem) !!!
+      DO jk = 1, isize
+         jl = (jk-1) /( ihdf_vars+nlay_i)+1
+         DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
+            DO ji = 1 , fs_jpim1   ! vector opt.
+               zflu(ji,jj) = pahu3D(ji,jj,jl) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj,jk) - ptab(ji,jj,jk) )
+               zflv(ji,jj) = pahv3D(ji,jj,jl) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1,jk) - ptab(ji,jj,jk) )
+            END DO
+         END DO
+         !
+         DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
+            DO ji = fs_2 , fs_jpim1   ! vector opt. 
+               zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
+               ptab(ji,jj,jk) = ztab0(ji,jj,jk) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj,jk) )
+            END DO
+         END DO
+      END DO
+
+      CALL lbc_lnk_multi( pt2d_array, type_array , psgn_array , isize ) ! Multiple interchange of all the variables
+
+      !!! final step (clem) !!!
+      ! -----------------------
+
+      IF(ln_ctl)   THEN
+         DO jk = 1 , isize
+            zrlx(:,:,jk) = ptab(:,:,jk) - ztab0(:,:,jk)
+            WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter
+            CALL prt_ctl( tab2d_1=zrlx(:,:,jk), clinfo1=charout )
+         END DO
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 )
+      CALL wrk_dealloc( jpi, jpj, zflu, zflv, zdiv )
+      CALL wrk_dealloc( jpi, jpj, jpl, pahu3D , pahv3D )
+
+      DEALLOCATE( zconv )
+      DEALLOCATE( pt2d_array , zrlx_array )
+      DEALLOCATE( type_array )
+      DEALLOCATE( psgn_array )
+      !
+   END SUBROUTINE lim_hdf_multiple
+
 
    
@@ -179,5 +374,5 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicehdf/ nn_convfrq
+      NAMELIST/namicehdf/ nn_convfrq 
       !!-------------------------------------------------------------------
       !
Index: /branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- /branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5578)
+++ /branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5579)
@@ -63,5 +63,5 @@
       INTEGER, INTENT(in) ::   kt           ! number of iteration
       !
-      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices
+      INTEGER  ::   ji, jj, jk, jm , jl, jt      ! dummy loop indices
       INTEGER  ::   initad                  ! number of sub-timestep for the advection
       REAL(wp) ::   zcfl , zusnit           !   -      -
@@ -75,6 +75,11 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhimax                   ! old ice thickness
       REAL(wp), POINTER, DIMENSION(:,:)      ::   zatold, zeiold, zesold   ! old concentration, enthalpies
+      REAL(wp), POINTER, DIMENSION(:,:,:)             ::   zhdfptab
       REAL(wp) ::    zdv, zvi, zvs, zsmv, zes, zei
       REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
+      !!---------------------------------------------------------------------
+      INTEGER                                ::  ihdf_vars  = 6  !!Number of variables in which we apply horizontal diffusion
+                                                                   !!  inside limtrp for each ice category , not counting the 
+                                                                   !!  variables corresponding to ice_layers 
       !!---------------------------------------------------------------------
       IF( nn_timing == 1 )  CALL timing_start('limtrp')
@@ -85,4 +90,5 @@
       CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
       CALL wrk_alloc( jpi,jpj,jpl,        zhimax, zviold, zvsold, zsmvold )
+      CALL wrk_alloc( jpi,jpj,jpl*(ihdf_vars + nlay_i),zhdfptab)
 
       IF( numit == nstart .AND. lwp ) THEN
@@ -170,5 +176,5 @@
             z0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
             z0es (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
-            DO jk = 1, nlay_i
+           DO jk = 1, nlay_i
                z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
             END DO
@@ -305,24 +311,41 @@
          !------------------------------------
          DO jl = 1, jpl
-         !                             ! Masked eddy diffusivity coefficient at ocean U- and V-points
-            DO jj = 1, jpjm1                 ! NB: has not to be defined on jpj line and jpi row
-               DO ji = 1 , fs_jpim1   ! vector opt.
-                  pahu(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji  ,jj,jl) ) ) )   &
-                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,jl) ) ) ) * ahiu(ji,jj)
-                  pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji,jj  ,jl) ) ) )   &
-                     &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- a_i(ji,jj+1,jl) ) ) ) * ahiv(ji,jj)
-               END DO
-            END DO
-
-            CALL lim_hdf( v_i  (:,:,  jl) )
-            CALL lim_hdf( v_s  (:,:,  jl) )
-            CALL lim_hdf( smv_i(:,:,  jl) )
-            CALL lim_hdf( oa_i (:,:,  jl) )
-            CALL lim_hdf( a_i  (:,:,  jl) )
-            CALL lim_hdf( e_s  (:,:,1,jl) )
+            jm=(jl-1)*(ihdf_vars+nlay_i)+1
+            zhdfptab(:,:,jm)= a_i  (:,:,  jl); jm = jm + 1   ! IMPORTANT  a_i must go in the first position because 
+                                                             ! it is needed at this position inside lim_hdf_multiple.
+            zhdfptab(:,:,jm)= v_i  (:,:,  jl); jm = jm + 1
+            zhdfptab(:,:,jm)= v_s  (:,:,  jl); jm = jm + 1 
+            zhdfptab(:,:,jm)= smv_i(:,:,  jl); jm = jm + 1
+            zhdfptab(:,:,jm)= oa_i (:,:,  jl); jm = jm + 1
+            zhdfptab(:,:,jm)= e_s  (:,:,1,jl); jm = jm + 1
+         ! Sample of adding more variables to apply lim_hdf using lim_hdf_multiple optimization---
+         !   zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1  
+         !   zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1 
+         !
+         ! and in this example the parameter ihdf_vars musb be changed to 8 (necessary for allocation)
+         !----------------------------------------------------------------------------------------
             DO jk = 1, nlay_i
-               CALL lim_hdf( e_i(:,:,jk,jl) )
-            END DO
-         END DO
+              zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1
+            END DO
+         END DO
+         CALL lim_hdf_multiple( zhdfptab, ihdf_vars, jpl, nlay_i) 
+
+         DO jl = 1, jpl
+            jm=(jl-1)*(ihdf_vars+nlay_i)+1
+            a_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1      
+            v_i  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
+            v_s  (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
+            smv_i(:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
+            oa_i (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1 
+            e_s  (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1 
+         ! Sample of adding more variables to apply lim_hdf---------
+         !   variable_1  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1 
+         !   variable_2  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1
+         !-----------------------------------------------------------
+            DO jk = 1, nlay_i
+               e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1 
+            END DO
+         END DO
+
 
          !------------------------------------------------------------------------------!
@@ -464,4 +487,5 @@
       CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
       CALL wrk_dealloc( jpi,jpj,jpl,        zviold, zvsold, zhimax, zsmvold )
+      CALL wrk_dealloc( jpi,jpj,jpl*(ihdf_vars+nlay_i),zhdfptab)
       !
       IF( nn_timing == 1 )  CALL timing_stop('limtrp')
Index: /branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90
===================================================================
--- /branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90	(revision 5578)
+++ /branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90	(revision 5579)
@@ -24,5 +24,5 @@
 
    INTERFACE lbc_lnk_multi
-      MODULE PROCEDURE mpp_lnk_2d_9
+      MODULE PROCEDURE mpp_lnk_2d_9 , mpp_lnk_2d_multiple
    END INTERFACE
 
Index: /branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
===================================================================
--- /branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 5578)
+++ /branches/2015/dev_r5546_CNRS19_HPC_scalability/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 5579)
@@ -62,4 +62,5 @@
    USE lbcnfd         ! north fold treatment
    USE in_out_manager ! I/O manager
+   USE wrk_nemo       ! work arrays
 
    IMPLICIT NONE
@@ -70,6 +71,7 @@
    PUBLIC   mpp_ini_north, mpp_lbc_north, mpp_lbc_north_e
    PUBLIC   mpp_min, mpp_max, mpp_sum, mpp_minloc, mpp_maxloc
+   PUBLIC   mpp_max_multiple
    PUBLIC   mpp_lnk_3d, mpp_lnk_3d_gather, mpp_lnk_2d, mpp_lnk_2d_e
-   PUBLIC   mpp_lnk_2d_9 
+   PUBLIC   mpp_lnk_2d_9 , mpp_lnk_2d_multiple 
    PUBLIC   mppscatter, mppgather
    PUBLIC   mpp_ini_ice, mpp_ini_znl
@@ -78,8 +80,10 @@
    PUBLIC   mpp_lnk_bdy_2d, mpp_lnk_bdy_3d
    PUBLIC   mpp_lbc_north_icb, mpp_lnk_2d_icb
+   PUBLIC   mpprank
 
    TYPE arrayptr
       REAL , DIMENSION (:,:),  POINTER :: pt2d
    END TYPE arrayptr
+   PUBLIC   arrayptr
    
    !! * Interfaces
@@ -105,4 +109,8 @@
    INTERFACE mpp_maxloc
       MODULE PROCEDURE mpp_maxloc2d ,mpp_maxloc3d
+   END INTERFACE
+
+   INTERFACE mpp_max_multiple
+      MODULE PROCEDURE mppmax_real_multiple
    END INTERFACE
 
@@ -724,16 +732,17 @@
       ! -----------------------
       !
-      DO ii = 1 , num_fields
          !First Array
-         IF( npolj /= 0 .AND. .NOT. PRESENT(cd_mpp) ) THEN
-            !
-            SELECT CASE ( jpni )
-            CASE ( 1 )     ;   CALL lbc_nfd      ( pt2d_array(ii)%pt2d( : , : ), type_array(ii) , psgn_array(ii) )   ! only 1 northern proc, no mpp
-            CASE DEFAULT   ;   CALL mpp_lbc_north( pt2d_array(ii)%pt2d( : , : ), type_array(ii), psgn_array(ii) )   ! for all northern procs.
-            END SELECT
-            !
-         ENDIF
-         !
-      END DO
+      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.
+         END SELECT
+         !
+      ENDIF
+        !
       
       DEALLOCATE( zt2ns, zt2sn, zt2ew, zt2we )
@@ -1681,4 +1690,29 @@
    END SUBROUTINE mppmax_real
 
+   SUBROUTINE mppmax_real_multiple( ptab, NUM , kcom  )
+      !!----------------------------------------------------------------------
+      !!                  ***  routine mppmax_real  ***
+      !!
+      !! ** Purpose :   Maximum
+      !!
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(:) ,  INTENT(inout)           ::   ptab   ! ???
+      INTEGER , INTENT(in   )           ::   NUM
+      INTEGER , INTENT(in   ), OPTIONAL ::   kcom   ! ???
+      !!
+      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)
+      !
+   END SUBROUTINE mppmax_real_multiple
+
 
    SUBROUTINE mppmin_a_real( ptab, kdim, kcom )
@@ -2575,4 +2609,165 @@
    END SUBROUTINE mpp_lbc_north_2d
 
+   SUBROUTINE mpp_lbc_north_2d_multiple( pt2d_array, cd_type, psgn, num_fields)
+      !!---------------------------------------------------------------------
+      !!                   ***  routine mpp_lbc_north_2d  ***
+      !!
+      !! ** Purpose :   Ensure proper north fold horizontal bondary condition
+      !!              in mpp configuration in case of jpn1 > 1
+      !!              (for multiple 2d arrays )
+      !!
+      !! ** Method  :   North fold condition and mpp with more than one proc
+      !!              in i-direction require a specific treatment. We gather
+      !!              the 4 northern lines of the global domain on 1 processor
+      !!              and apply lbc north-fold on this sub array. Then we
+      !!              scatter the north fold array back to the processors.
+      !!
+      !!----------------------------------------------------------------------
+      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
+      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
+      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) )
+      !
+      ijpj   = 4
+      ijpjm1 = 3
+      !
+      
+      DO jk = 1, num_fields
+         DO jj = nlcj-ijpj+1, nlcj             ! put in znorthloc the last 4 jlines of pt2d (for every variable)
+            ij = jj - nlcj + ijpj
+            znorthloc(:,ij,jk) = pt2d_array(jk)%pt2d(:,jj)
+         END DO
+      END DO
+      !                                     ! Build in procs of ncomm_north the znorthgloio
+      itaille = jpi * ijpj
+                                                                  
+      IF ( l_north_nogather ) THEN
+         !
+         ! Avoid the use of mpi_allgather by exchanging only with the processes already identified 
+         ! (in nemo_northcomms) as being  involved in this process' northern boundary exchange
+         !
+         ztabr(:,:,:) = 0
+         ztabl(:,:,:) = 0
+
+         DO jk = 1, num_fields
+            DO jj = nlcj-ijpj+1, nlcj          ! First put local values into the global array
+               ij = jj - nlcj + ijpj
+               DO ji = nfsloop, nfeloop
+                  ztabl(ji,ij,jk) = pt2d_array(jk)%pt2d(ji,jj)
+               END DO
+            END DO
+         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
+            ENDIF
+         END DO
+         DO jr = 1,nsndto
+            iproc = nfipproc(isendto(jr),jpnj)
+            IF(iproc .ne. -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
+                 DO jj = 1, ijpj
+                    DO ji = ildi, ilei
+                       ztabr(iilb+ji,jj,jk) = zfoldwk(ji,jj,jk)       ! Modified to 3D
+                    END DO
+                 END DO
+              END DO
+            ELSE IF (iproc .eq. (narea-1)) THEN
+              DO jk = 1, num_fields
+                 DO jj = 1, ijpj
+                    DO ji = ildi, ilei
+                          ztabr(iilb+ji,jj,jk) = pt2d_array(jk)%pt2d(ji,nlcj-ijpj+jj)       ! Modified to 3D
+                    END DO
+                 END DO
+              END DO
+            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
+                  CALL mpi_wait(ml_req_nf(jr), ml_stat, ml_err)
+               ENDIF
+            END DO
+         ENDIF
+         !
+         DO ji = 1, num_fields     ! 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 jj = nlcj-ijpj+1, nlcj             ! Scatter back to pt2d
+               ij = jj - nlcj + ijpj
+               DO ji = 1, nlci
+                  pt2d_array(jk)%pt2d(ji,jj) = ztabl(ji,ij,jk)       ! Modified to 3D
+               END DO
+            END DO
+         END DO
+         
+         !
+      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
+            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)
+               DO jj = 1, ijpj
+                  DO ji = ildi, ilei
+                     ztab(ji+iilb-1,jj,jk) = znorthgloio(ji,jj,jk,jr)
+                  END DO
+               END DO
+            END DO
+         END DO
+         
+         DO ji = 1, num_fields
+            CALL lbc_nfd( ztab(:,:,ji), cd_type(ji), psgn(ji) )   ! North fold boundary condition
+         END DO
+         !
+         DO jk = 1, num_fields
+            DO jj = nlcj-ijpj+1, nlcj             ! Scatter back to pt2d
+               ij = jj - nlcj + ijpj
+               DO ji = 1, nlci
+                  pt2d_array(jk)%pt2d(ji,jj) = ztab(ji+nimpp-1,ij,jk)
+               END DO
+            END DO
+         END DO
+         !
+         !
+      ENDIF
+      DEALLOCATE( ztab, znorthloc, zfoldwk, znorthgloio )
+      DEALLOCATE( ztabl, ztabr )
+      !
+   END SUBROUTINE mpp_lbc_north_2d_multiple
 
    SUBROUTINE mpp_lbc_north_e( pt2d, cd_type, psgn)
