Index: branches/UKMO/dev_r5518_GO6_missing_STABLE_revs/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
===================================================================
--- branches/UKMO/dev_r5518_GO6_missing_STABLE_revs/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 7939)
+++ branches/UKMO/dev_r5518_GO6_missing_STABLE_revs/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 7965)
@@ -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_e12t(ji,jj)
@@ -83,23 +108,82 @@
       !                             ! 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_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) = 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 +192,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_e12t(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_e12t(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 +233,5 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicehdf/ nn_convfrq
+      NAMELIST/namicehdf/ nn_convfrq 
       !!-------------------------------------------------------------------
       !
@@ -212,2 +266,3 @@
    !!======================================================================
 END MODULE limhdf
+
