Index: branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 7506)
+++ branches/2016/v3_6_CMIP6_ice_diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 7597)
@@ -325,114 +325,69 @@
          ! 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
+         IF( nn_ahi0 /= -1 ) THEN
+
+            ! --- Prepare diffusion for variables with categories --- !
+            !     mask eddy diffusivity coefficient at ocean U- and V-points
+            jm=1
+            DO jl = 1, jpl
+               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 (ihdf_vars must be increased)
+               !   zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1  
+               !   zhdfptab(:,:,jm) = variable_2 (:,:,1,jl); jm = jm + 1 
+               DO jk = 1, nlay_i
+                  zhdfptab(:,:,jm)=e_i(:,:,jk,jl); jm= jm+1
+               END DO
+            END DO
+            !
+            ! --- Prepare diffusion for open water area --- !
+            !     mask 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.
-                  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
-         !
-         !--------------------------------
-         !  diffusion of open water area
-         !--------------------------------
-         !                             ! 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
+                  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
+            !
+            zhdfptab(:,:,jm)= ato_i  (:,:);
+
+            ! --- Apply diffusion --- !
+            CALL lim_hdf( zhdfptab, ihdf_vars, jpl, nlay_i) 
+            
+            ! --- Recover properties --- !
+            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
+                  e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1 
+               END DO
+            END DO
+            ato_i  (:,:) = zhdfptab(:,:,jm)
+
+         ENDIF
          
-         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,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
-         !
-         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
-               e_i(:,:,jk,jl) = zhdfptab(:,:,jm);jm= jm + 1 
-            END DO
-         END DO
-
-         ato_i  (:,:) = zhdfptab(:,:,jm)
-
-         !------------------------------------------------------------------------------!
-         ! limit ice properties after transport                           
-         !------------------------------------------------------------------------------!
-!!gm & cr   :  MAX should not be active if adv scheme is positive !
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  v_s  (ji,jj,jl)   = MAX( 0._wp, v_s  (ji,jj,jl) )
-                  v_i  (ji,jj,jl)   = MAX( 0._wp, v_i  (ji,jj,jl) )
-                  smv_i(ji,jj,jl)   = MAX( 0._wp, smv_i(ji,jj,jl) )
-                  oa_i (ji,jj,jl)   = MAX( 0._wp, oa_i (ji,jj,jl) )
-                  a_i  (ji,jj,jl)   = MAX( 0._wp, a_i  (ji,jj,jl) )
-                  e_s  (ji,jj,1,jl) = MAX( 0._wp, e_s  (ji,jj,1,jl) )
-               END DO
-            END DO
-
-            DO jk = 1, nlay_i
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     e_i(ji,jj,jk,jl) = MAX( 0._wp, e_i(ji,jj,jk,jl) )
-                  END DO
-               END DO
-            END DO
-         END DO
-!!gm & cr 
-
          ! --- diags ---
          DO jj = 1, jpj
