Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8498)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8500)
@@ -422,52 +422,57 @@
             !
             jl1 = kdonor(ji,jl)
-            IF ( jl1 == jl  ) THEN   ;   jl2 = jl1+1
-            ELSE                     ;   jl2 = jl 
-            ENDIF
-            !
-            IF( v_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
-            ELSE                                  ;   zworkv(ji) = 0._wp
-            ENDIF
-            IF( a_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1)
-            ELSE                                  ;   zworka(ji) = 0._wp
-            ENDIF
-            !
-            a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas
-            a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
-            !
-            v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl)       ! Ice volumes
-            v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
-            !
-            ztrans         = v_s_2d(ji,jl1) * zworkv(ji)          ! Snow volumes
-            v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
-            v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 
-            !          
-            !                                                     ! Ice age 
-            ztrans             = oa_i_2d(ji,jl1) * pdaice(ji,jl)  !!clem: should be * zworka(ji) but it does not work ????
-            oa_i_2d(ji,jl1)    = oa_i_2d(ji,jl1) - ztrans
-            oa_i_2d(ji,jl2)    = oa_i_2d(ji,jl2) + ztrans
-            !
-            ztrans             = smv_i_2d(ji,jl1) * zworkv(ji)    ! Ice salinity
-
-            smv_i_2d(ji,jl1)   = smv_i_2d(ji,jl1) - ztrans
-            smv_i_2d(ji,jl2)   = smv_i_2d(ji,jl2) + ztrans
-            !
-            !                                                     ! Surface temperature
-            ztrans          = t_su_2d(ji,jl1) * pdaice(ji,jl)     !!clem: should be * zworka(ji) but it does not work ????
-            zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans
-            zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
-            !  
-            ! MV MP 2016 
-            IF ( nn_pnd_scheme > 0 ) THEN
-               !                                                  ! Pond fraction
-               ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
-               a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
-               a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
-               !                                                  ! Pond volume (also proportional to da/a)
-               ztrans          = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
-               v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
-               v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
-            ENDIF
-            ! END MV MP 2016 
+            !
+            IF( jl1 > 0 ) THEN
+               !
+               IF ( jl1 == jl  ) THEN   ;   jl2 = jl1+1
+               ELSE                     ;   jl2 = jl 
+               ENDIF
+               !
+               IF( v_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworkv(ji) = pdvice(ji,jl) / v_i_2d(ji,jl1)
+               ELSE                                  ;   zworkv(ji) = 0._wp
+               ENDIF
+               IF( a_i_2d(ji,jl1) >= epsi10 ) THEN   ;   zworka(ji) = pdaice(ji,jl) / a_i_2d(ji,jl1)
+               ELSE                                  ;   zworka(ji) = 0._wp
+               ENDIF
+               !
+               a_i_2d(ji,jl1) = a_i_2d(ji,jl1) - pdaice(ji,jl)       ! Ice areas
+               a_i_2d(ji,jl2) = a_i_2d(ji,jl2) + pdaice(ji,jl)
+               !
+               v_i_2d(ji,jl1) = v_i_2d(ji,jl1) - pdvice(ji,jl)       ! Ice volumes
+               v_i_2d(ji,jl2) = v_i_2d(ji,jl2) + pdvice(ji,jl)
+               !
+               ztrans         = v_s_2d(ji,jl1) * zworkv(ji)          ! Snow volumes
+               v_s_2d(ji,jl1) = v_s_2d(ji,jl1) - ztrans
+               v_s_2d(ji,jl2) = v_s_2d(ji,jl2) + ztrans 
+               !          
+               !                                                     ! Ice age 
+               ztrans             = oa_i_2d(ji,jl1) * pdaice(ji,jl)  !!clem: should be * zworka(ji) but it does not work ????
+               oa_i_2d(ji,jl1)    = oa_i_2d(ji,jl1) - ztrans
+               oa_i_2d(ji,jl2)    = oa_i_2d(ji,jl2) + ztrans
+               !
+               ztrans             = smv_i_2d(ji,jl1) * zworkv(ji)    ! Ice salinity
+               !
+               smv_i_2d(ji,jl1)   = smv_i_2d(ji,jl1) - ztrans
+               smv_i_2d(ji,jl2)   = smv_i_2d(ji,jl2) + ztrans
+               !
+               !                                                     ! Surface temperature
+               ztrans          = t_su_2d(ji,jl1) * pdaice(ji,jl)     !!clem: should be * zworka(ji) but it does not work ????
+               zaTsfn(ji,jl1)  = zaTsfn(ji,jl1) - ztrans
+               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
+               !  
+               ! MV MP 2016 
+               IF ( nn_pnd_scheme > 0 ) THEN
+                  !                                                  ! Pond fraction
+                  ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
+                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans
+                  a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans
+                  !                                                  ! Pond volume (also proportional to da/a)
+                  ztrans          = v_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
+                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
+                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
+               ENDIF
+               ! END MV MP 2016
+               !
+            ENDIF   ! jl1 >0
          END DO
          !
@@ -479,11 +484,14 @@
                !
                jl1 = kdonor(ji,jl)
-               IF(jl1 == jl) THEN  ;  jl2 = jl+1
-               ELSE                ;  jl2 = jl
-               ENDIF
-               !
-               ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji)
-               e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
-               e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans 
+               !
+               IF( jl1 > 0 ) THEN
+                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
+                  ELSE                ;  jl2 = jl
+                  ENDIF
+                  !
+                  ztrans            = e_s(ii,ij,jk,jl1) * zworkv(ji)
+                  e_s(ii,ij,jk,jl1) = e_s(ii,ij,jk,jl1) - ztrans
+                  e_s(ii,ij,jk,jl2) = e_s(ii,ij,jk,jl2) + ztrans
+               ENDIF
             END DO
          END DO
@@ -495,11 +503,14 @@
                !
                jl1 = kdonor(ji,jl)
-               IF(jl1 == jl) THEN  ;  jl2 = jl+1
-               ELSE                ;  jl2 = jl
-               ENDIF
-               !
-               ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji)
-               e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
-               e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans 
+               !
+               IF( jl1 > 0 ) THEN
+                  IF(jl1 == jl) THEN  ;  jl2 = jl+1
+                  ELSE                ;  jl2 = jl
+                  ENDIF
+                  !
+                  ztrans            = e_i(ii,ij,jk,jl1) * zworkv(ji)
+                  e_i(ii,ij,jk,jl1) = e_i(ii,ij,jk,jl1) - ztrans
+                  e_i(ii,ij,jk,jl2) = e_i(ii,ij,jk,jl2) + ztrans
+               ENDIF
             END DO
          END DO
@@ -510,15 +521,11 @@
       ! 3) Update ice thickness and temperature
       !-------------------------------------------------------------------------------
-      DO jl = 1, jpl
-         DO ji = 1, nidx
-            IF ( a_i_2d(ji,jl) > epsi10 ) THEN 
-               ht_i_2d(ji,jl)  =  v_i_2d(ji,jl) / a_i_2d(ji,jl) 
-               t_su_2d(ji,jl)  =  zaTsfn(ji,jl) / a_i_2d(ji,jl) 
-            ELSE
-               ht_i_2d(ji,jl)  = 0._wp
-               t_su_2d(ji,jl)  = rt0
-            ENDIF
-         END DO
-      END DO
+      WHERE( a_i_2d(1:nidx,:) >= epsi20 )
+         ht_i_2d(1:nidx,:)  =  v_i_2d(1:nidx,:) / a_i_2d(1:nidx,:) 
+         t_su_2d(1:nidx,:)  =  zaTsfn(1:nidx,:) / a_i_2d(1:nidx,:) 
+      ELSEWHERE
+         ht_i_2d(1:nidx,:)  = 0._wp
+         t_su_2d(1:nidx,:)  = rt0
+      END WHERE
       !
       CALL tab_2d_3d( nidx, idxice(1:nidx), ht_i_2d (1:nidx,1:jpl), ht_i  )
