Index: trunk/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6472)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6478)
@@ -234,5 +234,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_oce, v_oce   !: surface ocean velocity used in ice dynamics
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ahiu , ahiv    !: hor. diffusivity coeff. at U- and V-points [m2/s]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   pahu , pahv    !: ice hor. eddy diffusivity coef. at U- and V-points
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ust2s, hicol   !: friction velocity, ice collection thickness accreted in leads
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   strp1, strp2   !: strength at previous time steps
@@ -303,5 +302,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness [W.m-2]
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice   
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   pahu3D , pahv3D
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_amax_2d  !: maximum ice concentration 2d array
 
@@ -429,5 +429,4 @@
       ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) ,                           &
          &      ahiu     (jpi,jpj) , ahiv     (jpi,jpj) ,                           &
-         &      pahu     (jpi,jpj) , pahv     (jpi,jpj) ,                           &
          &      ust2s    (jpi,jpj) , hicol    (jpi,jpj) ,                           &
          &      strp1    (jpi,jpj) , strp2    (jpi,jpj) , strength  (jpi,jpj) ,     &
@@ -442,7 +441,7 @@
          &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,     &
          &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) ,                        &
-         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl), qlead  (jpi,jpj) ,                     &
-         &      rn_amax_2d(jpi,jpj),                                                            &
-         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) ,                       &
+         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl), pahu3D(jpi,jpj,jpl+1), pahv3D(jpi,jpj,jpl+1),            &
+         &      rn_amax_2d (jpi,jpj), qlead (jpi, jpj),                                                           &
+         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,                                          &
          &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,    &
          &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) ,     & 
@@ -514,2 +513,3 @@
    !!======================================================================
 END MODULE ice
+
Index: trunk/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 6472)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 6478)
@@ -7,4 +7,5 @@
    !!             -   !  2001-05 (G. Madec, R. Hordoir) opa norm
    !!            1.0  !  2002-08 (C. Ethe)  F90, free form
+   !!            3.0  !  2015-08 (O. Tintó and M. Castrillo)  added lim_hdf (multiple)
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -27,5 +28,5 @@
    PRIVATE
 
-   PUBLIC   lim_hdf         ! called by lim_trp
+   PUBLIC   lim_hdf ! called by lim_trp
    PUBLIC   lim_hdf_init    ! called by sbc_lim_init
 
@@ -43,5 +44,5 @@
 CONTAINS
 
-   SUBROUTINE lim_hdf( ptab )
+   SUBROUTINE lim_hdf( ptab , ihdf_vars , jpl , nlay_i )
       !!-------------------------------------------------------------------
       !!                  ***  ROUTINE lim_hdf  ***
@@ -54,10 +55,13 @@
       !! ** Action  :    update ptab with the diffusive contribution
       !!-------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT( inout ) ::   ptab    ! Field on which the diffusion is applied
-      !
-      INTEGER                           ::  ji, jj                    ! dummy loop indices
+      INTEGER                           :: jpl, nlay_i, isize, ihdf_vars
+      REAL(wp),  DIMENSION(:,:,:), INTENT( inout ),TARGET ::   ptab    ! Field on which the diffusion is applied
+      !
+      INTEGER                           ::  ji, jj, jk, jl , jm               ! dummy loop indices
       INTEGER                           ::  iter, ierr           ! local integers
-      REAL(wp)                          ::  zrlxint, zconv     ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:) ::  zrlx, zflu, zflv, zdiv0, zdiv, ztab0
+      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
@@ -65,8 +69,29 @@
       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  ==!
+      ! +1 open water diffusion
+      isize = jpl*(ihdf_vars+nlay_i)+1
+      ALLOCATE( zconv (isize) )
+      ALLOCATE( pt2d_array(isize) , zrlx_array(isize) )
+      ALLOCATE( type_array(isize) )
+      ALLOCATE( psgn_array(isize) )
       
-      CALL wrk_alloc( jpi, jpj, zrlx, zflu, zflv, zdiv0, zdiv, ztab0 )
-
-      !                       !==  Initialisation  ==!
+      CALL wrk_alloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 )
+      CALL wrk_alloc( jpi, jpj, zflu, zflv, zdiv )
+
+      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)
@@ -74,5 +99,5 @@
          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 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_e1e2t(ji,jj)
@@ -83,23 +108,81 @@
       !                             ! Time integration parameters
       !
-      ztab0(:, : ) = ptab(:,:)      ! Arrays initialization
-      zdiv0(:, 1 ) = 0._wp
-      zdiv0(:,jpj) = 0._wp
-      zflu (jpi,:) = 0._wp   
-      zflv (jpi,:) = 0._wp
-      zdiv0(1,  :) = 0._wp
-      zdiv0(jpi,:) = 0._wp
+      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( zconv > ( 2._wp * 1.e-04 ) .AND. iter <= its )   ! Sub-time step loop
+      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_e1e2t(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) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
-               zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
+               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
@@ -108,62 +191,32 @@
             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_e1e2t(ji,jj)
-            END DO
-         END DO
-         !
-         IF( iter == 1 )   zdiv0(:,:) = 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)    &
-                  &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) )   &
-                  &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )                               & 
-                  &      ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )
-               zrlx(ji,jj) = ptab(ji,jj) + zrelax * ( zrlxint - ptab(ji,jj) )
-            END DO
-         END DO
-         CALL lbc_lnk( zrlx, 'T', 1. )                   ! lateral boundary condition
-         !
-         IF ( MOD( iter, nn_convfrq ) == 0 )  THEN    ! convergence test every nn_convfrq iterations (perf. optimization)
-            zconv = 0._wp
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                  zconv = MAX( zconv, ABS( zrlx(ji,jj) - ptab(ji,jj) )  )
-               END DO
-            END DO
-            IF( lk_mpp )   CALL mpp_max( zconv )      ! max over the global domain
-         ENDIF
-         !
-         ptab(:,:) = zrlx(:,:)
-         !
-      END DO                                       ! end of sub-time step loop
-
-      ! -----------------------
-      !!! final step (clem) !!!
-      DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
-         DO ji = 1 , fs_jpim1   ! vector opt.
-            zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
-            zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(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
-      !
-      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_e1e2t(ji,jj)
-            ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) )
-         END DO
-      END DO
-      CALL lbc_lnk( ptab, 'T', 1. )                   ! lateral boundary condition
+
+      CALL lbc_lnk_multi( pt2d_array, type_array , psgn_array , isize ) ! Multiple interchange of all the variables
+
       !!! final step (clem) !!!
       ! -----------------------
 
       IF(ln_ctl)   THEN
-         zrlx(:,:) = ptab(:,:) - ztab0(:,:)
-         WRITE(charout,FMT="(' lim_hdf  : zconv =',D23.16, ' iter =',I4,2X)") zconv, iter
-         CALL prt_ctl( tab2d_1=zrlx, clinfo1=charout )
-      ENDIF
-      !
-      CALL wrk_dealloc( jpi, jpj, zrlx, zflu, zflv, zdiv0, zdiv, ztab0 )
+         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 )
+
+      DEALLOCATE( zconv )
+      DEALLOCATE( pt2d_array , zrlx_array )
+      DEALLOCATE( type_array )
+      DEALLOCATE( psgn_array )
       !
    END SUBROUTINE lim_hdf
+
 
    
@@ -179,5 +232,5 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicehdf/ nn_convfrq
+      NAMELIST/namicehdf/ nn_convfrq 
       !!-------------------------------------------------------------------
       !
@@ -212,2 +265,3 @@
    !!======================================================================
 END MODULE limhdf
+
Index: trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 6472)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 6478)
@@ -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)+1,zhdfptab)
 
       IF( numit == nstart .AND. lwp ) THEN
@@ -170,5 +176,5 @@
             z0oi (:,:,jl)   = oa_i (:,:,  jl) * e1e2t(:,:)  ! Age content
             z0es (:,:,jl)   = e_s  (:,:,1,jl) * e1e2t(:,:)  ! Snow heat content
-            DO jk = 1, nlay_i
+           DO jk = 1, nlay_i
                z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
             END DO
@@ -284,5 +290,43 @@
          ! Diffusion of Ice fields                  
          !------------------------------------------------------------------------------!
-
+         !------------------------------------
+         !  Diffusion of other ice variables
+         !------------------------------------
+         jm=1
+         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
+            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, -a_i(ji  ,jj,  jl ) ) ) )   &
+                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -a_i(ji+1,jj,  jl ) ) ) ) * ahiu(ji,jj)
+                  pahv3D(ji,jj,jl) = ( 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
+
+            zhdfptab(:,:,jm)= a_i  (:,:,  jl); jm = jm + 1    
+            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 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
+              zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1
+            END DO
+         END DO
          !
          !--------------------------------
@@ -290,39 +334,43 @@
          !--------------------------------
          !                             ! 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, -at_i(ji  ,jj) ) ) )   &
+         !         &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
+         !      pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
+         !         &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
+         !   END DO
+         !END DO
+         
          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, -at_i(ji  ,jj) ) ) )   &
-                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
-               pahv(ji,jj) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
-                  &        * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
+               pahu3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji  ,jj) ) ) )   &
+                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji+1,jj) ) ) ) * ahiu(ji,jj)
+               pahv3D(ji,jj,jpl+1) = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, -at_i(ji,jj  ) ) ) )   &
+                  &                * ( 1._wp - MAX( 0._wp, SIGN( 1._wp,- at_i(ji,jj+1) ) ) ) * ahiv(ji,jj)
             END DO
          END DO
          !
-         CALL lim_hdf( ato_i (:,:) )
-
-         !------------------------------------
-         !  Diffusion of other ice variables
-         !------------------------------------
-         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) )
+         zhdfptab(:,:,jm)= ato_i  (:,:);
+         CALL lim_hdf( zhdfptab, ihdf_vars, jpl, nlay_i) 
+
+         jm=1
+         DO jl = 1, jpl
+            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
-               CALL lim_hdf( e_i(:,:,jk,jl) )
-            END DO
-         END DO
+               e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1 
+            END DO
+         END DO
+
+         ato_i  (:,:) = zhdfptab(:,:,jm)
 
          !------------------------------------------------------------------------------!
@@ -464,4 +512,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)+1,zhdfptab)
       !
       IF( nn_timing == 1 )  CALL timing_stop('limtrp')
@@ -479,2 +528,3 @@
    !!======================================================================
 END MODULE limtrp
+
Index: trunk/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90	(revision 6472)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/LBC/lbclnk.F90	(revision 6478)
@@ -9,4 +9,5 @@
    !!            3.5  ! 2012     (S.Mocavero, I. Epicoco) optimization of BDY comm. via lbc_bdy_lnk and lbc_obc_lnk
    !!            3.4  ! 2012-12  (R. Bourdalle-Badie, G. Reffray)  add a C1D case  
+   !!            3.6  ! 2015-06  (O. Tintó and M. Castrillo)  add lbc_lnk_multi  
    !!----------------------------------------------------------------------
 #if defined key_mpp_mpi
@@ -22,5 +23,5 @@
 
    INTERFACE lbc_lnk_multi
-      MODULE PROCEDURE mpp_lnk_2d_9
+      MODULE PROCEDURE mpp_lnk_2d_9, mpp_lnk_2d_multiple
    END INTERFACE
    !
@@ -90,4 +91,8 @@
    END INTERFACE
    !
+   INTERFACE lbc_lnk_multi
+      MODULE PROCEDURE lbc_lnk_2d_9, lbc_lnk_2d_multiple
+   END INTERFACE
+
    INTERFACE lbc_bdy_lnk
       MODULE PROCEDURE lbc_bdy_lnk_2d, lbc_bdy_lnk_3d
@@ -97,7 +102,13 @@
       MODULE PROCEDURE lbc_lnk_2d_e
    END INTERFACE
+   
+   TYPE arrayptr
+      REAL , DIMENSION (:,:),  POINTER :: pt2d
+   END TYPE arrayptr
+   PUBLIC   arrayptr
 
    PUBLIC   lbc_lnk       ! ocean/ice  lateral boundary conditions
    PUBLIC   lbc_lnk_e     !
+   PUBLIC   lbc_lnk_multi ! modified ocean lateral boundary conditions
    PUBLIC   lbc_bdy_lnk   ! ocean lateral BDY boundary conditions
    PUBLIC   lbc_lnk_icb   !
@@ -181,4 +192,61 @@
       !
    END SUBROUTINE lbc_lnk_2d
+   
+   SUBROUTINE lbc_lnk_2d_multiple( pt2d_array , type_array , psgn_array , num_fields )
+      !!
+      INTEGER :: num_fields
+      TYPE( arrayptr ), DIMENSION(:) :: pt2d_array
+      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
+      !
+      INTEGER  ::   ii    !!MULTI SEND DUMMY LOOP INDICES
+      !
+      DO ii = 1, num_fields
+        CALL lbc_lnk_2d( pt2d_array(ii)%pt2d, type_array(ii), psgn_array(ii) )
+      END DO     
+      !
+   END SUBROUTINE lbc_lnk_2d_multiple
+
+   SUBROUTINE lbc_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)
+      !!
+      !!---------------------------------------------------------------------
+
+      !!The first array
+      CALL lbc_lnk( pt2dA, cd_typeA, psgnA ) 
+
+      !! Look if more arrays to process
+      IF(PRESENT (psgnB) )CALL lbc_lnk( pt2dA, cd_typeA, psgnA ) 
+      IF(PRESENT (psgnC) )CALL lbc_lnk( pt2dC, cd_typeC, psgnC ) 
+      IF(PRESENT (psgnD) )CALL lbc_lnk( pt2dD, cd_typeD, psgnD ) 
+      IF(PRESENT (psgnE) )CALL lbc_lnk( pt2dE, cd_typeE, psgnE ) 
+      IF(PRESENT (psgnF) )CALL lbc_lnk( pt2dF, cd_typeF, psgnF ) 
+      IF(PRESENT (psgnG) )CALL lbc_lnk( pt2dG, cd_typeG, psgnG ) 
+      IF(PRESENT (psgnH) )CALL lbc_lnk( pt2dH, cd_typeH, psgnH ) 
+      IF(PRESENT (psgnI) )CALL lbc_lnk( pt2dI, cd_typeI, psgnI ) 
+
+   END SUBROUTINE lbc_lnk_2d_9
+
+
+
+
 
 #else
@@ -379,4 +447,58 @@
       !    
    END SUBROUTINE lbc_lnk_2d
+   
+   SUBROUTINE lbc_lnk_2d_multiple( pt2d_array , type_array , psgn_array , num_fields )
+      !!
+      INTEGER :: num_fields
+      TYPE( arrayptr ), DIMENSION(:) :: pt2d_array
+      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
+      !
+      INTEGER  ::   ii    !!MULTI SEND DUMMY LOOP INDICES
+      !
+      DO ii = 1, num_fields
+        CALL lbc_lnk_2d( pt2d_array(ii)%pt2d, type_array(ii), psgn_array(ii) )
+      END DO     
+      !
+   END SUBROUTINE lbc_lnk_2d_multiple
+
+   SUBROUTINE lbc_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)
+      !!
+      !!---------------------------------------------------------------------
+
+      !!The first array
+      CALL lbc_lnk( pt2dA, cd_typeA, psgnA ) 
+
+      !! Look if more arrays to process
+      IF(PRESENT (psgnB) )CALL lbc_lnk( pt2dA, cd_typeA, psgnA ) 
+      IF(PRESENT (psgnC) )CALL lbc_lnk( pt2dC, cd_typeC, psgnC ) 
+      IF(PRESENT (psgnD) )CALL lbc_lnk( pt2dD, cd_typeD, psgnD ) 
+      IF(PRESENT (psgnE) )CALL lbc_lnk( pt2dE, cd_typeE, psgnE ) 
+      IF(PRESENT (psgnF) )CALL lbc_lnk( pt2dF, cd_typeF, psgnF ) 
+      IF(PRESENT (psgnG) )CALL lbc_lnk( pt2dG, cd_typeG, psgnG ) 
+      IF(PRESENT (psgnH) )CALL lbc_lnk( pt2dH, cd_typeH, psgnH ) 
+      IF(PRESENT (psgnI) )CALL lbc_lnk( pt2dI, cd_typeI, psgnI ) 
+
+   END SUBROUTINE lbc_lnk_2d_9
+
 
 #endif
@@ -448,2 +570,3 @@
    !!======================================================================
 END MODULE lbclnk
+
Index: trunk/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90
===================================================================
--- trunk/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 6472)
+++ trunk/NEMOGCM/NEMO/OPA_SRC/LBC/lib_mpp.F90	(revision 6478)
@@ -24,4 +24,5 @@
    !!            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' 
    !!----------------------------------------------------------------------
 
@@ -62,4 +63,5 @@
    USE lbcnfd         ! north fold treatment
    USE in_out_manager ! I/O manager
+   USE wrk_nemo       ! work arrays
 
    IMPLICIT NONE
@@ -70,6 +72,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   mpp_lnk_sum_3d, mpp_lnk_sum_2d
    PUBLIC   mppscatter, mppgather
@@ -79,8 +82,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
@@ -106,4 +111,8 @@
    INTERFACE mpp_maxloc
       MODULE PROCEDURE mpp_maxloc2d ,mpp_maxloc3d
+   END INTERFACE
+
+   INTERFACE mpp_max_multiple
+      MODULE PROCEDURE mppmax_real_multiple
    END INTERFACE
 
@@ -726,16 +735,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 )
@@ -2019,4 +2029,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 )
@@ -2912,4 +2947,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)
