Index: trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 6483)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 6490)
@@ -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
+
