Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8499)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8500)
@@ -475,9 +475,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   t_si          !: Temperature at Snow-ice interface (K) 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   tm_si         !: mean temperature at the snow-ice interface (K) 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dmi_dyn  !: Change in ice mass due to ice dynamics (kg/m2/s)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_dms_dyn  !: Change in snow mass due to ice dynamics (kg/m2/s)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_bo    !: Bottom conduction flux (W/m2)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_fc_su    !: Surface conduction flux (W/m2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   diag_shear    !: Maximum shear of sea-ice velocity field
 
    !
@@ -597,7 +594,6 @@
       ii = ii + 1
       ALLOCATE( t_si (jpi,jpj,jpl)    , tm_si(jpi,jpj)        ,    & 
-                diag_dmi_dyn(jpi,jpj) , diag_dms_dyn(jpi,jpj) ,    &
                 diag_fc_bo(jpi,jpj)   , diag_fc_su(jpi,jpj)   ,    &
-                diag_shear(jpi,jpj)   , STAT = ierr(ii) )
+                STAT = ierr(ii) )
 
       ice_alloc = MAXVAL( ierr(:) )
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8499)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8500)
@@ -29,4 +29,5 @@
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE timing         ! Timing
+   USE iom            !
 
    IMPLICIT NONE
@@ -68,5 +69,5 @@
       !
       REAL(wp) ::    zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b
-      REAL(wp) ::    zdv, zda
+      REAL(wp) ::    zdv
       REAL(wp), DIMENSION(jpi,jpj)           ::   zatold, zeiold, zesold, zsmvold 
       REAL(wp), DIMENSION(jpi,jpj,jpl)       ::   zhimax, zviold, zvsold
@@ -80,5 +81,4 @@
       ! MV MP 2016
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:)   ::   z0ap , z0vp
-      REAL(wp) ::   za_old
       ! END MV MP 2016
       REAL(wp), ALLOCATABLE, DIMENSION(:,:,:,:) ::   z0ei
@@ -111,13 +111,12 @@
       !--- Thickness correction init. --- !
       zatold(:,:) = at_i
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
-               ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
-               ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
-            END DO
-         END DO
-      END DO
+      WHERE( a_i(:,:,:) >= epsi20 )
+         ht_i(:,:,:) = v_i(:,:,:) / a_i(:,:,:)
+         ht_s(:,:,:) = v_s(:,:,:) / a_i(:,:,:)
+      ELSEWHERE
+         ht_i(:,:,:) = 0._wp
+         ht_s(:,:,:) = 0._wp         
+      END WHERE
+         
       ! --- Record max of the surrounding ice thicknesses for correction in case advection creates ice too thick --- !
       zhimax(:,:,:) = ht_i(:,:,:) + ht_s(:,:,:)
@@ -126,5 +125,5 @@
             DO ji = 2, jpim1
 !!gm use of MAXVAL here is very probably less efficient than expending the 9 values
-               zhimax(ji,jj,jl) = MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) )
+               zhimax(ji,jj,jl) = MAX( epsi20, MAXVAL( ht_i(ji-1:ji+1,jj-1:jj+1,jl) + ht_s(ji-1:ji+1,jj-1:jj+1,jl) ) )
             END DO
          END DO
@@ -399,6 +398,11 @@
          END DO
       END DO
-      
-      IF( nn_limdyn == 2) THEN
+      IF( iom_use('icetrp') )   CALL iom_put( "icetrp" , diag_trp_vi * rday  )         ! ice volume transport
+      IF( iom_use('snwtrp') )   CALL iom_put( "snwtrp" , diag_trp_vs * rday  )         ! snw volume transport
+      IF( iom_use('saltrp') )   CALL iom_put( "saltrp" , diag_trp_smv * rday * rhoic ) ! salt content transport
+      IF( iom_use('deitrp') )   CALL iom_put( "deitrp" , diag_trp_ei         )         ! advected ice enthalpy (W/m2)
+      IF( iom_use('destrp') )   CALL iom_put( "destrp" , diag_trp_es         )         ! advected snw enthalpy (W/m2)
+      
+      IF( nn_limdyn == 2 ) THEN
          !
          CALL ice_var_zapsmall      !--- zap small areas ---!
@@ -407,57 +411,29 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! 
                   IF ( v_i(ji,jj,jl) > 0._wp ) THEN
                      !
-                     rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
-                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
-                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
-                     !
-                     zdv  = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl)  
+                     ht_i  (ji,jj,jl) = v_i (ji,jj,jl) / a_i(ji,jj,jl)
+                     ht_s  (ji,jj,jl) = v_s (ji,jj,jl) / a_i(ji,jj,jl)
+                     zdv = v_i(ji,jj,jl) + v_s(ji,jj,jl) - zviold(ji,jj,jl) - zvsold(ji,jj,jl)  
                      !
                      IF ( ( zdv >  0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) .AND. zatold(ji,jj) < 0.80 ) .OR. &
                         & ( zdv <= 0.0 .AND. (ht_i(ji,jj,jl)+ht_s(ji,jj,jl)) > zhimax(ji,jj,jl) ) ) THEN
-                        !
-                        rswitch        = MAX( 0._wp, SIGN( 1._wp, zhimax(ji,jj,jl) - epsi20 ) )
-                        a_i(ji,jj,jl)  = rswitch * ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / MAX( zhimax(ji,jj,jl), epsi20 )
-                        !
-                        ! small correction due to *rswitch for a_i
-                        v_i  (ji,jj,jl)        = rswitch * v_i  (ji,jj,jl)
-                        v_s  (ji,jj,jl)        = rswitch * v_s  (ji,jj,jl)
-                        smv_i(ji,jj,jl)        = rswitch * smv_i(ji,jj,jl)
-                        e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl)
-                        e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
-
-                        ! MV MP 2016
-                        IF ( nn_pnd_scheme > 0 ) THEN
-                           a_ip (ji,jj,jl)        = rswitch * a_ip (ji,jj,jl)
-                           v_ip (ji,jj,jl)        = rswitch * v_ip (ji,jj,jl)
-                        ENDIF
-                        ! END MV MP 2016
-                        !
+                        a_i (ji,jj,jl) = ( v_i(ji,jj,jl) + v_s(ji,jj,jl) ) / zhimax(ji,jj,jl)
+                        ht_i(ji,jj,jl) =   v_i(ji,jj,jl) / a_i(ji,jj,jl)
                      ENDIF
                      !
                   ENDIF
-                  !
                END DO
             END DO
          END DO
-         
-         DO jj = 1, jpj             !--- bound ht_i to hi_max (99 m).
-            DO ji = 1, jpi
-               ! MV MP 2016
-               za_old = a_i(ji,jj,jpl)
-               ! END MV MP 2016
-               rswitch         = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) )
-               ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) )
-               a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch
-               ! MV MP 2016
-               IF ( nn_pnd_scheme > 0 ) THEN
-                  ! correct pond fraction to avoid a_ip > a_i
-                  a_ip(ji,jj,jpl) = a_ip(ji,jj,jpl) * a_i(ji,jj,jpl) / MAX( za_old , epsi20 ) * rswitch
-               ENDIF
-               ! END MP 2016
-            END DO
-         END DO
+                  
+         WHERE( ht_i(:,:,jpl) > hi_max(jpl) )             !--- bound ht_i to hi_max (99 m)
+            ht_i(:,:,jpl) = hi_max(jpl)
+            a_i (:,:,jpl) = v_i(:,:,jpl) / hi_max(jpl)
+         END WHERE
+
+         IF ( nn_pnd_scheme > 0 ) THEN                    !--- correct pond fraction to avoid a_ip > a_i
+            WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
+         ENDIF
          !
       ENDIF
@@ -467,24 +443,11 @@
       !------------------------------------------------------------
       !
-!!gm remplace the test by, l_piling a logical compute one for all in icestp.F90 (and its declaration in ice.F90
-!!gm      IF ( nn_limdyn == 1 .OR. ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) ) THEN ! simple conservative piling, comparable with LIM2
-      IF( l_piling ) THEN
+      IF( l_piling ) THEN ! simple conservative piling, comparable with 1-cat models
          at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
          DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch = MAX( 0._wp, SIGN( 1._wp, at_i(ji,jj) - epsi20 ) )
-                  zda     = rswitch * MIN( rn_amax_2d(ji,jj) - at_i(ji,jj), 0._wp ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj), epsi20 )
-                  a_i(ji,jj,jl) = a_i(ji,jj,jl) + zda
-               END DO
-            END DO
-         END DO
-!!gm better and faster coding?
-!         DO jl = 1, jpl
-!            WHERE( at_i(:,:) > epsi20 )
-!               a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
-!            END WHERE
-!         END DO
-!!gm end
+            WHERE( at_i(:,:) > epsi20 )
+               a_i(:,:,jl) = a_i(:,:,jl) * (  1._wp + MIN( rn_amax_2d(:,:) - at_i(:,:) , 0._wp ) / at_i(:,:)  )
+            END WHERE
+         END DO
       ENDIF
       
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 8499)
+++ /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  )
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8499)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8500)
@@ -612,6 +612,4 @@
                &                                      + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
 
-!!!            wfx_snw_dyn(ji,jj) = wfx_snw_dyn(ji,jj) + zwfx_snw
-
             hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
                &                                - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
@@ -734,8 +732,4 @@
          !
       END DO ! jl1 (deforming categories)
-
-      ! SIMIP diagnostics
-      diag_dmi_dyn(:,:) = - wfx_dyn(:,:)     + rhoic * diag_trp_vi(:,:)
-      diag_dms_dyn(:,:) = - wfx_snw_dyn(:,:) + rhosn * diag_trp_vs(:,:)      
       !
    END SUBROUTINE ice_rdgrft_ridgeshift
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90	(revision 8499)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90	(revision 8500)
@@ -153,5 +153,6 @@
       REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
       !! --- diags
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zswi, zsig1, zsig2, zsig3
+      REAL(wp), DIMENSION(jpi,jpj) ::   zswi
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
       !! --- SIMIP diags
       REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_sig1      ! Average normal stress in sea ice   
@@ -692,15 +693,20 @@
       ! 5) diagnostics
       !------------------------------------------------------------------------------!
-      ! --- ellipse --- !
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
+         END DO
+      END DO
+
+      ! --- divergence, shear and strength --- !
+      IF( iom_use('idive')  )   CALL iom_put( "idive"  , divu_i(:,:)   * zswi(:,:) )   ! divergence
+      IF( iom_use('ishear') )   CALL iom_put( "ishear" , shear_i(:,:)  * zswi(:,:) )   ! shear
+      IF( iom_use('icestr') )   CALL iom_put( "icestr" , strength(:,:) * zswi(:,:) )   ! Ice strength
+
+      ! --- charge ellipse --- !
       IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') ) THEN
          !
-         ALLOCATE( zswi(jpi,jpj) , zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
-         !
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
-            END DO
-         END DO
-         
+         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
+         !         
          DO jj = 2, jpjm1
             DO ji = 2, jpim1
@@ -728,5 +734,5 @@
          IF( iom_use('isig3') )   CALL iom_put( "isig3" , zsig3 )
          !
-         DEALLOCATE( zswi , zsig1 , zsig2 , zsig3 )
+         DEALLOCATE( zsig1 , zsig2 , zsig3 )
       ENDIF
       
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8499)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8500)
@@ -571,5 +571,4 @@
             
             ! SIMIP diagnostics
-            diag_dms_dyn(ji,jj)  = 0._wp ; diag_dmi_dyn(ji,jj)  = 0._wp
             diag_fc_bo(ji,jj)    = 0._wp ; diag_fc_su(ji,jj)    = 0._wp
             
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90	(revision 8499)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dif.F90	(revision 8500)
@@ -738,6 +738,8 @@
       ! MV SIMIP 2016
       !--- Conduction fluxes (positive downwards)
-      diag_fc_bo_1d(:)   = diag_fc_bo_1d(:) + fc_bo_i(:) * a_i_1d(:) / at_i_1d(:)
-      diag_fc_su_1d(:)   = diag_fc_su_1d(:) + fc_su(:)   * a_i_1d(:) / at_i_1d(:)
+      DO ji = 1, nidx
+         diag_fc_bo_1d(ji)   = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji)
+         diag_fc_su_1d(ji)   = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji)
+      END DO
       ! END MV SIMIP 2016
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8499)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8500)
@@ -116,8 +116,8 @@
          !
          ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) )
-         WHERE( at_i(:,:) > epsi10 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
+         WHERE( at_i(:,:) > epsi20 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
          ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
          END WHERE
-         WHERE( vt_i(:,:) > epsi10 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
+         WHERE( vt_i(:,:) > epsi20 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
          ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp
          END WHERE
@@ -163,6 +163,4 @@
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
       !!------------------------------------------------------------------
-
-!!gm Question 1:  here use epsi20 , in ice_var_agg it is epsi10 which is used...  why ??
 
 !!gm Question 2:  It is possible to define existence of sea-ice in a common way between 
@@ -549,8 +547,8 @@
       ! to be sure that at_i is the sum of a_i(jl)
       at_i (:,:) = a_i(:,:,1)
-	  vt_i (:,:) = v_i(:,:,1)
+      vt_i (:,:) = v_i(:,:,1)
       DO jl = 2, jpl
          at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-	     vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
+         vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
       END DO
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90	(revision 8499)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90	(revision 8500)
@@ -145,14 +145,5 @@
       CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity
       CALL iom_put( "snowvol"     , vt_s    * zswi      )        ! snow volume
-
-      CALL iom_put( "idive"       , divu_i(:,:)  * zswi(:,:) + zmiss(:,:) )   ! divergence
-      CALL iom_put( "ishear"      , shear_i(:,:) * zswi(:,:) + zmiss(:,:) )   ! shear
       
-      CALL iom_put( "icetrp"      , diag_trp_vi * rday  )        ! ice volume transport
-      CALL iom_put( "snwtrp"      , diag_trp_vs * rday  )        ! snw volume transport
-      CALL iom_put( "saltrp"      , diag_trp_smv * rday * rhoic ) ! salt content transport
-      CALL iom_put( "deitrp"      , diag_trp_ei         )        ! advected ice enthalpy (W/m2)
-      CALL iom_put( "destrp"      , diag_trp_es         )        ! advected snw enthalpy (W/m2)
-
       CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from bottom growth
       CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melting
@@ -275,5 +266,5 @@
       IF ( iom_use( "dmithd"   ) ) CALL iom_put( "dmithd"      , - wfx_bog - wfx_bom - wfx_sum   &                       ! Sea-ice mass change from thermodynamics
               &                                                   - wfx_sni - wfx_opw - wfx_res )
-      IF ( iom_use( "dmidyn"   ) ) CALL iom_put( "dmidyn"      ,   diag_dmi_dyn             )                            ! Sea-ice mass change from dynamics
+      IF ( iom_use( "dmidyn"   ) ) CALL iom_put( "dmidyn"      , - wfx_dyn(:,:) + rhoic * diag_trp_vi(:,:) )             ! Sea-ice mass change from dynamics(kg/m2/s)
       IF ( iom_use( "dmiopw"   ) ) CALL iom_put( "dmiopw"      , - wfx_opw                  )                            ! Sea-ice mass change through growth in open water
       IF ( iom_use( "dmibog"   ) ) CALL iom_put( "dmibog"      , - wfx_bog                  )                            ! Sea-ice mass change through basal growth
@@ -290,5 +281,5 @@
 
       IF ( iom_use( "dmsmel"   ) ) CALL iom_put( "dmsmel"      , - wfx_snw_sum              )                            ! Snow mass change through melt
-      IF ( iom_use( "dmsdyn"   ) ) CALL iom_put( "dmsdyn"      ,   diag_dms_dyn             )                            ! Snow mass change through dynamics
+      IF ( iom_use( "dmsdyn"   ) ) CALL iom_put( "dmsdyn"      , - wfx_snw_dyn(:,:) + rhosn * diag_trp_vs(:,:) )         ! Snow mass change through dynamics(kg/m2/s)
 
       IF ( iom_use( "hfxsenso" ) ) CALL iom_put( "hfxsenso"    ,   -fhtur(:,:)              * zswi(:,:) + zmiss(:,:) )   ! Sensible oceanic heat flux
@@ -305,7 +296,4 @@
       IF ( iom_use( "utau_ice" ) ) CALL iom_put( "utau_ice"     ,  utau_ice(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Wind stress term in force balance (x)
       IF ( iom_use( "vtau_ice" ) ) CALL iom_put( "vtau_ice"     ,  vtau_ice(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Wind stress term in force balance (y)
-
-
-      IF ( iom_use( "icestr"   ) ) CALL iom_put( "icestr"      ,   strength(:,:)            * zswi(:,:) + zmiss(:,:) )   ! Ice strength
 
       !--------------------------------
