Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90	(revision 8558)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90	(revision 8559)
@@ -100,5 +100,4 @@
    !                                                     ! to reintegrate longwave flux inside the ice thermodynamics
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hicol_1d      !: Ice collection thickness accumulated in leads
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_su_1d       !: <==> the 2D  t_su
@@ -197,5 +196,5 @@
          &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  &
          &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  &
-         &      sfx_lam_1d (jpij) , sfx_dyn_1d(jpij), hicol_1d   (jpij)    , STAT=ierr(ii) )
+         &      sfx_lam_1d (jpij) , sfx_dyn_1d(jpij)  , STAT=ierr(ii) )
       !
       ii = ii + 1
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 8558)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 8559)
@@ -293,7 +293,7 @@
       DO jk = 1, nlay_i
          DO ji = 1, nidx
-            ztmelts           = - tmut * s_i_1d(ji,jk) + rt0          ! Melting point of layer k [K]
-            
-            IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
+            ztmelts           = - tmut * s_i_1d(ji,jk)          ! Melting point of layer k [C]
+            
+            IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting
 
                zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]       
@@ -320,5 +320,5 @@
                
                zEi            = - e_i_1d(ji,jk) * r1_rhoic            ! Specific enthalpy of layer k [J/kg, <0]
-               zEw            =    rcp * ( ztmelts - rt0 )            ! Specific enthalpy of resulting meltwater [J/kg, <0]
+               zEw            =    rcp * ztmelts                      ! Specific enthalpy of resulting meltwater [J/kg, <0]
                zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
                
@@ -430,11 +430,11 @@
                                   + ( 1. - zswitch_sal ) * sm_i_1d(ji) 
                ! New ice growth
-               ztmelts            = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
+               ztmelts            = - tmut * s_i_new(ji)                 ! New ice melting point (C)
 
                zt_i_new           = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
                
-               zEi                = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)      
-                  &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
-                  &               + rcp  * ( ztmelts-rt0 )          
+               zEi                = cpic * ( zt_i_new - (ztmelts+rt0) ) &     ! Specific enthalpy of forming ice (J/kg, <0)      
+                  &               - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) )   &
+                  &               + rcp  * ztmelts          
 
                zEw                = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
@@ -450,11 +450,11 @@
             zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0)
 
-            ztmelts        = - tmut * s_i_new(ji) + rt0          ! New ice melting point (K)
+            ztmelts        = - tmut * s_i_new(ji)                ! New ice melting point (C)
             
             zt_i_new       = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
             
-            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)      
-               &               - lfus * ( 1.0 - ( ztmelts - rt0 ) / ( zt_i_new - rt0 ) )   &
-               &               + rcp  * ( ztmelts-rt0 )          
+            zEi            = cpic * ( zt_i_new - (ztmelts+rt0) ) &     ! Specific enthalpy of forming ice (J/kg, <0)      
+               &               - lfus * ( 1.0 - ztmelts / ( zt_i_new - rt0 ) )   &
+               &               + rcp  * ztmelts          
             
             zEw            = rcp  * ( t_bo_1d(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0)
@@ -490,7 +490,7 @@
             IF(  zf_tt(ji)  >  0._wp  .AND. jk > icount(ji,jk) ) THEN   ! do not calculate where layer has already disappeared by surface melting 
 
-               ztmelts = - tmut * s_i_1d(ji,jk) + rt0  ! Melting point of layer jk (K)
-
-               IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
+               ztmelts = - tmut * s_i_1d(ji,jk)  ! Melting point of layer jk (C)
+
+               IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !!! Internal melting
 
                   zEi               = - e_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
@@ -520,5 +520,5 @@
 
                   zEi             = - e_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
-                  zEw             = rcp * ( ztmelts - rt0 )    ! Specific enthalpy of meltwater (J/kg, <0)
+                  zEw             = rcp * ztmelts              ! Specific enthalpy of meltwater (J/kg, <0)
                   zdE             = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)
 
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90	(revision 8558)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_do.F90	(revision 8559)
@@ -80,5 +80,5 @@
       INTEGER  ::   ji,jj,jk,jl      ! dummy loop indices
       INTEGER  ::   iter     !   -       -
-      REAL(wp) ::   ztmelts, zdv, zfrazb, zweight, zde                          ! local scalars
+      REAL(wp) ::   ztmelts, zfrazb, zweight, zde                          ! local scalars
       REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      -
       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
@@ -102,5 +102,4 @@
       REAL(wp), DIMENSION(jpij) ::   zdv_res     ! residual volume in case of excessive heat budget
       REAL(wp), DIMENSION(jpij) ::   zda_res     ! residual area in case of excessive heat budget
-      REAL(wp), DIMENSION(jpij) ::   zat_i_1d    ! total ice fraction    
       REAL(wp), DIMENSION(jpij) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
       REAL(wp), DIMENSION(jpij) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector)
@@ -108,9 +107,6 @@
       REAL(wp), DIMENSION(jpij,jpl) ::   zv_b      ! old volume of ice in category jl
       REAL(wp), DIMENSION(jpij,jpl) ::   za_b      ! old area of ice in category jl
-      REAL(wp), DIMENSION(jpij,jpl) ::   za_i_1d   ! 1-D version of a_i
-      REAL(wp), DIMENSION(jpij,jpl) ::   zv_i_1d   ! 1-D version of v_i
-      REAL(wp), DIMENSION(jpij,jpl) ::   zsmv_i_1d ! 1-D version of smv_i
-
-      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_1d !: 1-D version of e_i
+
+      REAL(wp), DIMENSION(jpij,nlay_i,jpl) ::   ze_i_2d !: 1-D version of e_i
 
       REAL(wp), DIMENSION(jpi,jpj) ::   zvrel     ! relative ice / frazil velocity
@@ -243,25 +239,24 @@
       IF ( nidx > 0 ) THEN
 
-         CALL tab_2d_1d( nidx, idxice(1:nidx), zat_i_1d  (1:nidx)     , at_i           )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), at_i_1d (1:nidx)      , at_i         )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i  (:,:,:) )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i  (:,:,:) )
+         CALL tab_3d_2d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) )
          DO jl = 1, jpl
-            CALL tab_2d_1d( nidx, idxice(1:nidx), za_i_1d  (1:nidx,jl), a_i  (:,:,jl)  )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), zv_i_1d  (1:nidx,jl), v_i  (:,:,jl)  )
-            CALL tab_2d_1d( nidx, idxice(1:nidx), zsmv_i_1d(1:nidx,jl), smv_i(:,:,jl)  )
             DO jk = 1, nlay_i
-               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl)   )
-            END DO
-         END DO
-
-         CALL tab_2d_1d( nidx, idxice(1:nidx), qlead_1d  (1:nidx)     , qlead       )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d   (1:nidx)     , t_bo        )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx)     , sfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx)     , wfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hicol_1d  (1:nidx)     , hicol       )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d  (1:nidx)     , zvrel       )
-
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx)     , hfx_thd     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx)     , hfx_opw     )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), rn_amax_1d(1:nidx)     , rn_amax_2d  )
-         CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d    (1:nidx)     , sss_m       )
+               CALL tab_2d_1d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )
+            END DO
+         END DO
+         CALL tab_2d_1d( nidx, idxice(1:nidx), qlead_1d  (1:nidx) , qlead       )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), t_bo_1d   (1:nidx) , t_bo        )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx) , sfx_opw     )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx) , wfx_opw     )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), zh_newice (1:nidx) , hicol       )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), zvrel_1d  (1:nidx) , zvrel       )
+
+         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx) , hfx_thd     )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx) , hfx_opw     )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), rn_amax_1d(1:nidx) , rn_amax_2d  )
+         CALL tab_2d_1d( nidx, idxice(1:nidx), sss_1d    (1:nidx) , sss_m       )
 
          !------------------------------------------------------------------------------|
@@ -269,8 +264,10 @@
          !------------------------------------------------------------------------------|
          DO jl = 1, jpl
-            DO jk = 1, nlay_i
-               DO ji = 1, nidx
-                  IF( zv_i_1d(ji,jl) > 0._wp )   ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) / zv_i_1d(ji,jl) * REAL( nlay_i )
-               END DO
+            DO jk = 1, nlay_i               
+               WHERE( v_i_2d(1:nidx,jl) > 0._wp )
+                  ze_i_2d(1:nidx,jk,jl) = ze_i_2d(1:nidx,jk,jl) / v_i_2d(1:nidx,jl) * REAL( nlay_i )
+               ELSEWHERE
+                  ze_i_2d(1:nidx,jk,jl) = 0._wp
+               END WHERE
             END DO
          END DO
@@ -282,11 +279,6 @@
          ! Keep old ice areas and volume in memory
          !-----------------------------------------
-         zv_b(1:nidx,:) = zv_i_1d(1:nidx,:) 
-         za_b(1:nidx,:) = za_i_1d(1:nidx,:)
-
-         !----------------------
-         ! Thickness of new ice
-         !----------------------
-         zh_newice(1:nidx) = hicol_1d(1:nidx)
+         zv_b(1:nidx,:) = v_i_2d(1:nidx,:) 
+         za_b(1:nidx,:) = a_i_2d(1:nidx,:)
 
          !----------------------
@@ -309,8 +301,8 @@
          ! We assume that new ice is formed at the seawater freezing point
          DO ji = 1, nidx
-            ztmelts       = - tmut * zs_newice(ji) + rt0                  ! Melting point (K)
-            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_1d(ji) )                                         &
-               &                       + lfus * ( 1.0 - ( ztmelts - rt0 ) / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
-               &                       - rcp  *         ( ztmelts - rt0 )  )
+            ztmelts       = - tmut * zs_newice(ji)                  ! Melting point (C)
+            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - ( t_bo_1d(ji) - rt0 ) )                     &
+               &                       + lfus * ( 1.0 - ztmelts / MIN( t_bo_1d(ji) - rt0, -epsi10 ) )   &
+               &                       - rcp  *         ztmelts )
          END DO
 
@@ -318,7 +310,5 @@
          ! Age of new ice
          !----------------
-         DO ji = 1, nidx
-            zo_newice(ji) = 0._wp
-         END DO
+         zo_newice(1:nidx) = 0._wp
 
          !-------------------
@@ -350,9 +340,9 @@
          END DO
          
-         zv_frazb(:) = 0._wp
+         zv_frazb(1:nidx) = 0._wp
          IF( ln_frazil ) THEN
             ! A fraction zfrazb of frazil ice is accreted at the ice bottom
             DO ji = 1, nidx
-               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
+               rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - at_i_1d(ji) ) )
                zfrazb        = rswitch * ( TANH( rn_Cfraz * ( zvrel_1d(ji) - rn_vfraz ) ) + 1.0 ) * 0.5 * rn_maxfraz
                zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
@@ -378,6 +368,6 @@
          ! we keep the excessive volume in memory and attribute it later to bottom accretion
          DO ji = 1, nidx
-            IF ( za_newice(ji) >  ( rn_amax_1d(ji) - zat_i_1d(ji) ) ) THEN
-               zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - zat_i_1d(ji) )
+            IF ( za_newice(ji) >  ( rn_amax_1d(ji) - at_i_1d(ji) ) ) THEN
+               zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - at_i_1d(ji) )
                zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
                za_newice(ji) = za_newice(ji) - zda_res  (ji)
@@ -390,15 +380,14 @@
 
          ! find which category to fill
-         zat_i_1d(:) = 0._wp
          DO jl = 1, jpl
             DO ji = 1, nidx
                IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN
-                  za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji)
-                  zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji)
-                  jcat    (ji)    = jl
+                  a_i_2d(ji,jl) = a_i_2d(ji,jl) + za_newice(ji)
+                  v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newice(ji)
+                  jcat(ji) = jl
                ENDIF
-               zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d  (ji,jl)
-            END DO
-         END DO
+            END DO
+         END DO
+         at_i_1d(1:nidx) = SUM( a_i_2d(1:nidx,:), dim=2 )
 
          ! Heat content
@@ -411,8 +400,8 @@
             DO ji = 1, nidx
                jl = jcat(ji)
-               rswitch = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) )
-               ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    &
-                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_b(ji,jl) )  &
-                  &        * rswitch / MAX( zv_i_1d(ji,jl), epsi20 )
+               rswitch = MAX( 0._wp, SIGN( 1._wp , v_i_2d(ji,jl) - epsi20 ) )
+               ze_i_2d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                    &
+                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_2d(ji,jk,jl) * zv_b(ji,jl) )  &
+                  &        * rswitch / MAX( v_i_2d(ji,jl), epsi20 )
             END DO
          END DO
@@ -428,6 +417,6 @@
             DO jk = 1, nlay_i
                DO ji = 1, nidx
-                  h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i
-                  eh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk)
+                  h_i_old (ji,jk) = v_i_2d(ji,jl) * r1_nlay_i
+                  eh_i_old(ji,jk) = ze_i_2d(ji,jk,jl) * h_i_old(ji,jk)
                END DO
             END DO
@@ -435,8 +424,8 @@
             ! new volumes including lateral/bottom accretion + residual
             DO ji = 1, nidx
-               rswitch        = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) )
-               zv_newfra      = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 )
-               za_i_1d(ji,jl) = rswitch * za_i_1d(ji,jl)               
-               zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra
+               rswitch        = MAX( 0._wp, SIGN( 1._wp , at_i_1d(ji) - epsi20 ) )
+               zv_newfra     = rswitch * ( zdv_res(ji) + zv_frazb(ji) ) * a_i_2d(ji,jl) / MAX( at_i_1d(ji) , epsi20 )
+               a_i_2d(ji,jl) = rswitch * a_i_2d(ji,jl)               
+               v_i_2d(ji,jl) = v_i_2d(ji,jl) + zv_newfra
                ! for remapping
                h_i_old (ji,nlay_i+1) = zv_newfra
@@ -444,5 +433,5 @@
             ENDDO
             ! --- Ice enthalpy remapping --- !
-            CALL ice_thd_ent( ze_i_1d(1:nidx,:,jl) ) 
+            CALL ice_thd_ent( ze_i_2d(1:nidx,:,jl) ) 
          ENDDO
 
@@ -452,6 +441,5 @@
          DO jl = 1, jpl
             DO ji = 1, nidx
-               zdv   = zv_i_1d(ji,jl) - zv_b(ji,jl)
-               zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji)
+               smv_i_2d(ji,jl) = smv_i_2d(ji,jl) + zs_newice(ji) * ( v_i_2d(ji,jl) - zv_b(ji,jl) )
             END DO
          END DO
@@ -462,7 +450,5 @@
          DO jl = 1, jpl
             DO jk = 1, nlay_i
-               DO ji = 1, nidx
-                  ze_i_1d(ji,jk,jl) = ze_i_1d(ji,jk,jl) * zv_i_1d(ji,jl) * r1_nlay_i 
-               END DO
+               ze_i_2d(1:nidx,jk,jl) = ze_i_2d(1:nidx,jk,jl) * v_i_2d(1:nidx,jl) * r1_nlay_i 
             END DO
          END DO
@@ -470,16 +456,16 @@
          ! 7) Change 2D vectors to 1D vectors 
          !------------------------------------------------------------------------------!
-         DO jl = 1, jpl
-            CALL tab_1d_2d( nidx, idxice(1:nidx), za_i_1d (1:nidx,jl), a_i (:,:,jl)  )
-            CALL tab_1d_2d( nidx, idxice(1:nidx), zv_i_1d (1:nidx,jl), v_i (:,:,jl)  )
-            CALL tab_1d_2d( nidx, idxice(1:nidx), zsmv_i_1d(1:nidx,jl), smv_i (:,:,jl)   )
+         CALL tab_2d_3d( nidx, idxice(1:nidx), a_i_2d  (1:nidx,1:jpl), a_i  (:,:,:) )
+         CALL tab_2d_3d( nidx, idxice(1:nidx), v_i_2d  (1:nidx,1:jpl), v_i  (:,:,:) )
+         CALL tab_2d_3d( nidx, idxice(1:nidx), smv_i_2d(1:nidx,1:jpl), smv_i(:,:,:) )
+          DO jl = 1, jpl
             DO jk = 1, nlay_i
-               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_1d(1:nidx,jk,jl), e_i(:,:,jk,jl)  )
-            END DO
-         END DO
-         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx), sfx_opw  )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx), wfx_opw  )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd  )
-         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx), hfx_opw  )
+               CALL tab_1d_2d( nidx, idxice(1:nidx), ze_i_2d(1:nidx,jk,jl), e_i(:,:,jk,jl) )
+            END DO
+         END DO
+         CALL tab_1d_2d( nidx, idxice(1:nidx), sfx_opw_1d(1:nidx), sfx_opw )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), wfx_opw_1d(1:nidx), wfx_opw )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_thd_1d(1:nidx), hfx_thd )
+         CALL tab_1d_2d( nidx, idxice(1:nidx), hfx_opw_1d(1:nidx), hfx_opw )
          !
       ENDIF ! nidx > 0
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8558)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8559)
@@ -59,4 +59,5 @@
       REAL(wp) ::   iflush, igravdr              ! local scalars
       REAL(wp) ::   zs_sni, zsm_i_gd, zsm_i_fl, zsm_i_si, zsm_i_bg   ! local scalars
+      REAL(wp) ::   z1_time_gd, z1_time_fl
       !!---------------------------------------------------------------------
 
@@ -66,4 +67,7 @@
       CASE( 2 )       !  time varying salinity with linear profile  !
       !               !---------------------------------------------!
+         z1_time_gd = 1._wp / rn_time_gd * rdt_ice
+         z1_time_fl = 1._wp / rn_time_fl * rdt_ice
+         !
          DO ji = 1, nidx
 
@@ -71,11 +75,11 @@
             !  Update ice salinity from snow-ice and bottom growth
             !---------------------------------------------------------
-            zs_sni   = sss_1d(ji) * ( rhoic - rhosn ) * r1_rhoic   ! Salinity of snow ice
-            rswitch  = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
-            zsm_i_si = ( zs_sni      - sm_i_1d(ji) ) *             dh_snowice(ji)  / MAX( ht_i_1d(ji), epsi20 ) * rswitch ! snow-ice    
-            zsm_i_bg = ( s_i_new(ji) - sm_i_1d(ji) ) * MAX( 0._wp, dh_i_bott(ji) ) / MAX( ht_i_1d(ji), epsi20 ) * rswitch ! bottom growth
-
-            ! Update salinity (nb: salt flux already included in icethd_dh)
-            sm_i_1d(ji) = sm_i_1d(ji) + zsm_i_bg + zsm_i_si
+            IF( ht_i_1d(ji) > 0._wp ) THEN
+               zs_sni   = sss_1d(ji) * ( rhoic - rhosn ) * r1_rhoic                                 ! Salinity of snow ice
+               zsm_i_si = ( zs_sni      - sm_i_1d(ji) ) *             dh_snowice(ji)  / ht_i_1d(ji) ! snow-ice    
+               zsm_i_bg = ( s_i_new(ji) - sm_i_1d(ji) ) * MAX( 0._wp, dh_i_bott(ji) ) / ht_i_1d(ji) ! bottom growth
+               ! Update salinity (nb: salt flux already included in icethd_dh)
+               sm_i_1d(ji) = sm_i_1d(ji) + zsm_i_bg + zsm_i_si
+            ENDIF
 
             IF( ld_sal ) THEN
@@ -85,6 +89,7 @@
                iflush   = MAX( 0._wp , SIGN( 1._wp , t_su_1d(ji) - rt0         ) )  ! =1 if summer 
                igravdr  = MAX( 0._wp , SIGN( 1._wp , t_bo_1d(ji) - t_su_1d(ji) ) )  ! =1 if t_su < t_bo
-               zsm_i_gd = - igravdr * MAX( sm_i_1d(ji) - rn_sal_gd , 0._wp ) / rn_time_gd * rdt_ice  ! gravity drainage 
-               zsm_i_fl = - iflush  * MAX( sm_i_1d(ji) - rn_sal_fl , 0._wp ) / rn_time_fl * rdt_ice  ! flushing
+
+               zsm_i_gd = - igravdr * MAX( sm_i_1d(ji) - rn_sal_gd , 0._wp ) * z1_time_gd  ! gravity drainage 
+               zsm_i_fl = - iflush  * MAX( sm_i_1d(ji) - rn_sal_fl , 0._wp ) * z1_time_fl  ! flushing
                
                ! Update salinity   
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 8558)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8559)
@@ -42,8 +42,8 @@
    !!   ice_var_eqv2glo   : transform from VEQV to VGLO
    !!   ice_var_salprof   : salinity profile in the ice
-   !!   ice_var_bv        : brine volume
    !!   ice_var_salprof1d : salinity profile in the ice 1D
    !!   ice_var_zapsmall  : remove very small area and volume
    !!   ice_var_itd       : convert 1-cat to multiple cat
+   !!   ice_var_bv        : brine volume
    !!----------------------------------------------------------------------
    USE dom_oce        ! ocean space and time domain
@@ -64,8 +64,8 @@
    PUBLIC   ice_var_eqv2glo      
    PUBLIC   ice_var_salprof      
-   PUBLIC   ice_var_bv           
    PUBLIC   ice_var_salprof1d    
    PUBLIC   ice_var_zapsmall
    PUBLIC   ice_var_itd
+   PUBLIC   ice_var_bv           
 
    !!----------------------------------------------------------------------
@@ -170,4 +170,8 @@
       END WHERE
       !
+      WHERE( v_i(:,:,:) > epsi20 )   ;   z1_v_i(:,:,:) = 1._wp / v_i(:,:,:)
+      ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp
+      END WHERE
+      !
       ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)    !--- ice thickness
 
@@ -177,5 +181,5 @@
          ht_i  (:,:,jpl) = zhmax
          a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 
-         z1_a_i(:,:,jpl) = zhmax / v_i(:,:,jpl)          ! NB: v_i always /=0 as ht_i > hi_max
+         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)          ! NB: v_i always /=0 as ht_i > hi_max
       END WHERE
 
@@ -185,5 +189,5 @@
 
       IF( nn_icesal == 2 ) THEN                    !--- salinity (with a minimum value imposed everywhere)
-         WHERE( v_i(:,:,:) > epsi20 )   ;   sm_i(:,:,:) = MAX( rn_simin , smv_i(:,:,:) / v_i(:,:,:) )
+         WHERE( v_i(:,:,:) > epsi20 )   ;   sm_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, smv_i(:,:,:) * z1_v_i(:,:,:) ) )
          ELSEWHERE                      ;   sm_i(:,:,:) = rn_simin
          END WHERE
@@ -202,5 +206,5 @@
                   IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area 
                      !
-                     ze_i             =   e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3]
+                     ze_i             =   e_i(ji,jj,jk,jl) * z1_v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3]
                      ztmelts          = - s_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
                      ! Conversion q(S,T) -> T (second order equation)
@@ -295,6 +299,8 @@
          ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
          !
-         DO jk = 1, nlay_i
-            s_i(:,:,jk,:)  = sm_i(:,:,:)
+         DO jl = 1, jpl
+            DO jk = 1, nlay_i
+               s_i(:,:,jk,jl)  = sm_i(:,:,jl)
+            END DO
          END DO
          !                                      ! Slope of the linear profile 
@@ -353,35 +359,4 @@
    END SUBROUTINE ice_var_salprof
 
-
-   SUBROUTINE ice_var_bv
-      !!-------------------------------------------------------------------
-      !!                ***  ROUTINE ice_var_bv ***
-      !!
-      !! ** Purpose :   computes mean brine volume (%) in sea ice
-      !!
-      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
-      !!
-      !! References : Vancoppenolle et al., JGR, 2007
-      !!-------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      !!-------------------------------------------------------------------
-      !
-!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
-!!   instead of setting everything to zero as just below
-      bv_i (:,:,:) = 0._wp
-      DO jl = 1, jpl
-         DO jk = 1, nlay_i
-            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
-               bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
-            END WHERE
-         END DO
-      END DO
-      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
-      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
-      END WHERE
-      !
-   END SUBROUTINE ice_var_bv
-
-
    SUBROUTINE ice_var_salprof1d
       !!-------------------------------------------------------------------
@@ -393,7 +368,7 @@
       INTEGER  ::   ji, jk    ! dummy loop indices
       REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars
-      REAL(wp) ::   zalpha, zs, zs0              !   -      -
-      !
-      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s   !
+      REAL(wp) ::   zs, zs0              !   -      -
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s, zalpha   !
       REAL(wp), PARAMETER :: zsi0 = 3.5_wp
       REAL(wp), PARAMETER :: zsi1 = 4.5_wp
@@ -405,5 +380,5 @@
       CASE( 1 )       !  constant salinity in time and space  !
          !            !---------------------------------------!
-         s_i_1d(:,:) = rn_icesal
+         s_i_1d(1:nidx,:) = rn_icesal
          !
          !            !---------------------------------------------!
@@ -411,24 +386,29 @@
          !            !---------------------------------------------!
          !
-         ALLOCATE( z_slope_s(jpij) )
-         !
-         DO ji = 1, nidx          ! Slope of the linear profile
-            rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
-            z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )
-         END DO
-
+         ALLOCATE( z_slope_s(jpij), zalpha(jpij) )
+         !
+         !                                      ! Slope of the linear profile 
+         WHERE( ht_i_1d(1:nidx) > epsi20 )   ;   z_slope_s(1:nidx) = 2._wp * sm_i_1d(1:nidx) / ht_i_1d(1:nidx)
+         ELSEWHERE                           ;   z_slope_s(1:nidx) = 0._wp
+         END WHERE
+         
          z1_dS = 1._wp / ( zsi1 - zsi0 )
+         DO ji = 1, nidx
+            zalpha(ji) = MAX(  0._wp , MIN(  ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp  )  )
+            !                             ! force a constant profile when SSS too low (Baltic Sea)
+            IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) )   zalpha(ji) = 0._wp
+         END DO
+         !
+         ! Computation of the profile
          DO jk = 1, nlay_i
             DO ji = 1, nidx
-               zalpha = MAX(  0._wp , MIN(  ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp  )  )
-               IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) )   zalpha = 0._wp               ! cst profile when SSS too low (Baltic Sea)
-               !
-               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i   ! linear profile with 0 surface value
-               zs  = zalpha * zs0 + ( 1._wp - zalpha ) * sm_i_1d(ji)                      ! weighting the profile
-               s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )                     ! bounding salinity
-            END DO
-         END DO
-         !
-         DEALLOCATE( z_slope_s )
+               !                          ! linear profile with 0 surface value
+               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
+               zs  = zalpha(ji) * zs0 + ( 1._wp - zalpha(ji) ) * sm_i_1d(ji)
+               s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )
+            END DO
+         END DO
+         !
+         DEALLOCATE( z_slope_s, zalpha )
 
          !            !-------------------------------------------!
@@ -436,5 +416,5 @@
          !            !-------------------------------------------!                                   (mean = 2.30)
          !
-         sm_i_1d(:) = 2.30_wp
+         sm_i_1d(1:nidx) = 2.30_wp
          !
 !!gm cf remark in ice_var_salprof routine, CASE( 3 )
@@ -459,5 +439,5 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
-      REAL(wp) ::   zsal, zvi, zvs, zei, zes, zvp
+      REAL(wp), DIMENSION(jpi,jpj) ::   zswitch
       !!-------------------------------------------------------------------
       !
@@ -467,17 +447,19 @@
          ! Zap ice energy and use ocean heat to melt ice
          !-----------------------------------------------------------------
+         WHERE( a_i(:,:,jl) > epsi20 )   ;   ht_i(:,:,jl) = v_i(:,:,jl) / a_i(:,:,jl)
+         ELSEWHERE                       ;   ht_i(:,:,jl) = 0._wp
+         END WHERE
+         !
+         WHERE( a_i(:,:,jl) < epsi10 .OR. v_i(:,:,jl) < epsi10 .OR. ht_i(:,:,jl) < epsi10 )   ;   zswitch(:,:) = 0._wp
+         ELSEWHERE                                                                            ;   zswitch(:,:) = 1._wp
+         END WHERE
+         
          DO jk = 1, nlay_i
             DO jj = 1 , jpj
                DO ji = 1 , jpi
-!!gm I think we can do better/faster  :  to be discussed
-                  rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
-                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
-                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
-                     &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
-                  zei              = e_i(ji,jj,jk,jl)
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
-                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
                   ! update exchanges with ocean
-                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
+                  hfx_res(ji,jj)   = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_i(ji,jj,jk,jl) * r1_rdtice ! W.m-2 <0
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * zswitch(ji,jj)
+                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
                END DO
             END DO
@@ -486,48 +468,43 @@
          DO jj = 1 , jpj
             DO ji = 1 , jpi
-               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
-               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
-               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
-                  &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
-               zsal = smv_i(ji,jj,  jl)
-               zvi  = v_i  (ji,jj,  jl)
-               zvs  = v_s  (ji,jj,  jl)
-               zes  = e_s  (ji,jj,1,jl)
-               IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw )   zvp  = v_ip (ji,jj  ,jl)
+               ! update exchanges with ocean
+               sfx_res(ji,jj)  = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * smv_i(ji,jj,jl)   * rhoic * r1_rdtice
+               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i  (ji,jj,jl)   * rhoic * r1_rdtice
+               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s  (ji,jj,jl)   * rhosn * r1_rdtice
+               hfx_res(ji,jj)  = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s  (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0
                !-----------------------------------------------------------------
                ! Zap snow energy 
                !-----------------------------------------------------------------
-               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
-               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
+               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * zswitch(ji,jj) + rt0 * ( 1._wp - zswitch(ji,jj) )
+               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * zswitch(ji,jj)
 
                !-----------------------------------------------------------------
                ! zap ice and snow volume, add water and salt to ocean
                !-----------------------------------------------------------------
-               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
-               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch
-               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch
-               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch
-               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
-               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
-               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
-
-               ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * rswitch
-               ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * rswitch
-
-               ! MV MP 2016
-               IF ( ln_pnd ) THEN 
-                  a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch
-                  v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch
-                  IF ( ln_pnd_fw )   wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_ip(ji,jj,jl)  - zvp  ) * rhofw * r1_rdtice
-               ENDIF
-               ! END MV MP 2016
-
-               ! update exchanges with ocean
-               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
-               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
-               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
-               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
-            END DO
-         END DO
+               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - zswitch(ji,jj) ) + ato_i(ji,jj)
+               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * zswitch(ji,jj)
+               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * zswitch(ji,jj)
+               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * zswitch(ji,jj)
+               t_su (ji,jj,jl) = t_su (ji,jj,jl) * zswitch(ji,jj) + t_bo(ji,jj) * ( 1._wp - zswitch(ji,jj) )
+               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * zswitch(ji,jj)
+               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * zswitch(ji,jj)
+
+               ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * zswitch(ji,jj)
+               ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * zswitch(ji,jj)
+
+            END DO
+         END DO
+
+         IF( ln_pnd ) THEN 
+            DO jj = 1 , jpj
+               DO ji = 1 , jpi
+                  IF( ln_pnd_fw )   &
+                     &   wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
+                  a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
+                  v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
+               END DO
+            END DO
+         ENDIF
+         
       END DO 
 
@@ -688,4 +665,35 @@
     END SUBROUTINE ice_var_itd
 
+
+    SUBROUTINE ice_var_bv
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE ice_var_bv ***
+      !!
+      !! ** Purpose :   computes mean brine volume (%) in sea ice
+      !!
+      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
+      !!
+      !! References : Vancoppenolle et al., JGR, 2007
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
+      !!-------------------------------------------------------------------
+      !
+!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
+!!   instead of setting everything to zero as just below
+      bv_i (:,:,:) = 0._wp
+      DO jl = 1, jpl
+         DO jk = 1, nlay_i
+            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
+               bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
+            END WHERE
+         END DO
+      END DO
+      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
+      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
+      END WHERE
+      !
+   END SUBROUTINE ice_var_bv
+
+
 #else
    !!----------------------------------------------------------------------
