Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90	(revision 5064)
@@ -21,5 +21,4 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fcor       !: coriolis coefficient
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   area       !: surface of grid cell 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tms, tmi   !: temperature mask, mask for stress
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tmu, tmv   !: mask at u and v velocity points
@@ -42,5 +41,5 @@
       !!-------------------------------------------------------------------
       !
-      ALLOCATE( fcor(jpi,jpj)   , area(jpi,jpj) ,      &
+      ALLOCATE( fcor(jpi,jpj)   ,                      &
          &      tms   (jpi,jpj) , tmi (jpi,jpj) ,      &
          &      tmu   (jpi,jpj) , tmv (jpi,jpj) ,      &
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5064)
@@ -109,7 +109,7 @@
    !! smv_i       |      -      |    Sea ice salt content         | ppt.m |
    !! oa_i        !      -      !    Sea ice areal age content    | day   |
-   !! e_i         !      -      !    Ice enthalpy                 | 10^9 J| 
+   !! e_i         !      -      !    Ice enthalpy                 | J/m2  | 
    !!      -      ! q_i_1d      !    Ice enthalpy per unit vol.   | J/m3  | 
-   !! e_s         !      -      !    Snow enthalpy                | 10^9 J| 
+   !! e_s         !      -      !    Snow enthalpy                | J/m2  | 
    !!      -      ! q_s_1d      !    Snow enthalpy per unit vol.  | J/m3  | 
    !!                                                                     |
@@ -224,5 +224,4 @@
 
    !                                     !!** define some parameters 
-   REAL(wp), PUBLIC, PARAMETER ::   unit_fac = 1.e+09_wp  !: conversion factor for ice / snow enthalpy
    REAL(wp), PUBLIC, PARAMETER ::   epsi06   = 1.e-06_wp  !: small number 
    REAL(wp), PUBLIC, PARAMETER ::   epsi10   = 1.e-10_wp  !: small number 
@@ -331,5 +330,5 @@
       
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_i        !: ice temperatures          [K]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i        !: ice thermal contents [Giga J]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i        !: ice thermal contents    [J/m2]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   s_i        !: ice salinities          [PSU]
 
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90	(revision 5064)
@@ -97,5 +97,5 @@
 
       !  Initialize volumes of boxes  (=area if adv_x first called, =psm otherwise)                                     
-      psm (:,:)  = MAX( pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
+      psm (:,:)  = MAX( pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20 )
 
       !  Calculate fluxes and moments between boxes i<-->i+1              
@@ -282,5 +282,5 @@
 
       !  Initialize volumes of boxes (=area if adv_x first called, =psm otherwise)
-      psm(:,:)  = MAX(  pcrh * area(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
+      psm(:,:)  = MAX(  pcrh * e12t(:,:) + ( 1.0 - pcrh ) * psm(:,:) , epsi20  )
 
       !  Calculate fluxes and moments between boxes j<-->j+1              
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 5064)
@@ -167,19 +167,23 @@
       REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft
       REAL(wp)                        :: zvmin, zamin, zamax 
+      REAL(wp)                        :: zconv
+
+      zconv = 1.e-9
 
       IF( icount == 0 ) THEN
 
-         zvi_b  = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) )
-         zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )
-         zei_b  = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) )
+         zvi_b  = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tms(:,:) )
+         zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tms(:,:) )
+         zei_b  = glob_sum( SUM( SUM(   e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv +  &
+            &               SUM( SUM(   e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv )
          zfw_b  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
             &                   wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    &
-            &             ) * area(:,:) * tms(:,:) )
+            &             ) * e12t(:,:) * tms(:,:) )
          zfs_b  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
             &                   sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &
-            &                 ) * area(:,:) * tms(:,:) )
+            &                 ) * e12t(:,:) * tms(:,:) )
          zft_b  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
             &                 - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
-            &                  ) * area(:,:) / unit_fac * tms(:,:) )
+            &                  ) * e12t(:,:) * tms(:,:) * zconv )
 
       ELSEIF( icount == 1 ) THEN
@@ -187,15 +191,16 @@
          zfs  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
             &                 sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
-            &                ) * area(:,:) * tms(:,:) ) - zfs_b
+            &                ) * e12t(:,:) * tms(:,:) ) - zfs_b
          zfw  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
             &                 wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    &
-            &                ) * area(:,:) * tms(:,:) ) - zfw_b
+            &                ) * e12t(:,:) * tms(:,:) ) - zfw_b
          zft  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
             &               - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
-            &                ) * area(:,:) / unit_fac * tms(:,:) ) - zft_b
+            &                ) * e12t(:,:) * tms(:,:) * zconv ) - zft_b
  
-         zvi  = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw 
-         zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic )
-         zei  =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zei_b * r1_rdtice + zft
+         zvi  = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw 
+         zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tms(:,:) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic )
+         zei  =   glob_sum( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv +   &
+            &               SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tms * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft
 
          zvmin = glob_min(v_i)
@@ -206,10 +211,10 @@
             IF ( ABS( zvi    ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (',cd_routine,') = ',(zvi * rday)
             IF ( ABS( zsmv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday)
-            IF ( ABS( zei    ) >  1.    ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (',cd_routine,') = ',(zei)
-            IF ( zvmin <  1.e-10        ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin)
-            IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > amax+1.e-10 ) THEN
+            IF ( ABS( zei    ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',(zei)
+            IF ( zvmin <  -epsi10       ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin)
+            IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > amax+epsi10 ) THEN
                                           WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
             ENDIF
-            IF ( zamin <  1.e-10        ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
+            IF ( zamin <  -epsi10       ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin
          ENDIF
 
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90	(revision 5064)
@@ -306,5 +306,5 @@
                WRITE(numout,*) ' - Cell values '
                WRITE(numout,*) '   ~~~~~~~~~~~ '
-               WRITE(numout,*) ' cell area     : ', area(ji,jj)
+               WRITE(numout,*) ' cell area     : ', e12t(ji,jj)
                WRITE(numout,*) ' at_i          : ', at_i(ji,jj)       
                WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)       
@@ -317,6 +317,6 @@
                   WRITE(numout,*) ' v_i           : ', v_i(ji,jj,jl)
                   WRITE(numout,*) ' v_s           : ', v_s(ji,jj,jl)
-                  WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1,jl)/1.0e9
-                  WRITE(numout,*) ' e_i           : ', e_i(ji,jj,1:nlay_i,jl)/1.0e9
+                  WRITE(numout,*) ' e_s           : ', e_s(ji,jj,1,jl)
+                  WRITE(numout,*) ' e_i           : ', e_i(ji,jj,1:nlay_i,jl)
                   WRITE(numout,*) ' t_su          : ', t_su(ji,jj,jl)
                   WRITE(numout,*) ' t_snow        : ', t_s(ji,jj,1,jl)
@@ -350,5 +350,5 @@
                WRITE(numout,*) ' - Cell values '
                WRITE(numout,*) '   ~~~~~~~~~~~ '
-               WRITE(numout,*) ' cell area     : ', area(ji,jj)
+               WRITE(numout,*) ' cell area     : ', e12t(ji,jj)
                WRITE(numout,*) ' at_i          : ', at_i(ji,jj)       
                WRITE(numout,*) ' vt_i          : ', vt_i(ji,jj)       
@@ -372,6 +372,6 @@
                   WRITE(numout,*) ' v_i        : ', v_i(ji,jj,jl)              , ' v_i_b      : ', v_i_b(ji,jj,jl)   
                   WRITE(numout,*) ' v_s        : ', v_s(ji,jj,jl)              , ' v_s_b      : ', v_s_b(ji,jj,jl)  
-                  WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)/1.0e9      , ' ei1        : ', e_i_b(ji,jj,1,jl)/1.0e9 
-                  WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)/1.0e9      , ' ei2_b      : ', e_i_b(ji,jj,2,jl)/1.0e9  
+                  WRITE(numout,*) ' e_i1       : ', e_i(ji,jj,1,jl)            , ' ei1        : ', e_i_b(ji,jj,1,jl) 
+                  WRITE(numout,*) ' e_i2       : ', e_i(ji,jj,2,jl)            , ' ei2_b      : ', e_i_b(ji,jj,2,jl)  
                   WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' e_snow_b   : ', e_s_b(ji,jj,1,jl) 
                   WRITE(numout,*) ' smv_i      : ', smv_i(ji,jj,jl)            , ' smv_i_b    : ', smv_i_b(ji,jj,jl)   
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 5064)
@@ -71,71 +71,71 @@
 
       ! 1/area
-      z1_area = 1._wp / MAX( glob_sum( area(:,:) * tms(:,:) ), epsi06 )
-
-      rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) )
+      z1_area = 1._wp / MAX( glob_sum( e12t(:,:) * tms(:,:) ), epsi06 )
+
+      rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( e12t(:,:) * tms(:,:) ) - epsi06 ) )
       ! ----------------------- !
       ! 1 -  Content variations !
       ! ----------------------- !
-      zbg_ivo = glob_sum( vt_i(:,:) * area(:,:) * tms(:,:) ) ! volume ice 
-      zbg_svo = glob_sum( vt_s(:,:) * area(:,:) * tms(:,:) ) ! volume snow
-      zbg_are = glob_sum( at_i(:,:) * area(:,:) * tms(:,:) ) ! area
-      zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )       ! mean salt content
-      zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * area(:,:) * tms(:,:) )  ! mean temp content
-
-      !zbg_ihc = glob_sum( et_i(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content
-      !zbg_shc = glob_sum( et_s(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_svo,epsi06 ) ! snow heat content
+      zbg_ivo = glob_sum( vt_i(:,:) * e12t(:,:) * tms(:,:) ) ! volume ice 
+      zbg_svo = glob_sum( vt_s(:,:) * e12t(:,:) * tms(:,:) ) ! volume snow
+      zbg_are = glob_sum( at_i(:,:) * e12t(:,:) * tms(:,:) ) ! area
+      zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tms(:,:) )       ! mean salt content
+      zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * e12t(:,:) * tms(:,:) )  ! mean temp content
+
+      !zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * tms(:,:) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content
+      !zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * tms(:,:) ) / MAX( zbg_svo,epsi06 ) ! snow heat content
 
       ! Volume
       ztmp = rswitch * z1_area * r1_rau0 * rday
-      zbg_vfx     = ztmp * glob_sum(     emp(:,:) * area(:,:) * tms(:,:) )
-      zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * area(:,:) * tms(:,:) )
-      zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * area(:,:) * tms(:,:) )
-      zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * area(:,:) * tms(:,:) )
-      zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * area(:,:) * tms(:,:) )
-      zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * area(:,:) * tms(:,:) )
-      zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * area(:,:) * tms(:,:) )
-      zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * area(:,:) * tms(:,:) )
-      zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * area(:,:) * tms(:,:) )
-      zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * area(:,:) * tms(:,:) )
-      zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * area(:,:) * tms(:,:) )
+      zbg_vfx     = ztmp * glob_sum(     emp(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * e12t(:,:) * tms(:,:) )
 
       ! Salt
-      zbg_sfx     = ztmp * glob_sum(     sfx(:,:) * area(:,:) * tms(:,:) )
-      zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * area(:,:) * tms(:,:) )
-      zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * area(:,:) * tms(:,:) )
-      zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * area(:,:) * tms(:,:) )
-
-      zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * area(:,:) * tms(:,:) )
-      zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * area(:,:) * tms(:,:) )
-      zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * area(:,:) * tms(:,:) )
-      zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * area(:,:) * tms(:,:) )
-      zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * area(:,:) * tms(:,:) )
+      zbg_sfx     = ztmp * glob_sum(     sfx(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * e12t(:,:) * tms(:,:) )
+
+      zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tms(:,:) )
+      zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tms(:,:) )
 
       ! Heat budget
       zbg_ihc      = glob_sum( et_i(:,:) * 1.e-20 ) ! ice heat content  [1.e-20 J]
       zbg_shc      = glob_sum( et_s(:,:) * 1.e-20 ) ! snow heat content [1.e-20 J]
-      zbg_hfx_dhc  = glob_sum( diag_heat_dhc(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-
-      zbg_hfx_thd  = glob_sum( hfx_thd(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_dyn  = glob_sum( hfx_dyn(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_res  = glob_sum( hfx_res(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_sub  = glob_sum( hfx_sub(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_snw  = glob_sum( hfx_snw(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_sum  = glob_sum( hfx_sum(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_bom  = glob_sum( hfx_bom(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_bog  = glob_sum( hfx_bog(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_dif  = glob_sum( hfx_dif(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_opw  = glob_sum( hfx_opw(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_out  = glob_sum( hfx_out(:,:) * area(:,:) * tms(:,:) ) ! [in W]
-      zbg_hfx_in   = glob_sum(  hfx_in(:,:) * area(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_dhc  = glob_sum( diag_heat_dhc(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+
+      zbg_hfx_thd  = glob_sum( hfx_thd(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_dyn  = glob_sum( hfx_dyn(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_res  = glob_sum( hfx_res(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_sub  = glob_sum( hfx_sub(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_snw  = glob_sum( hfx_snw(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_sum  = glob_sum( hfx_sum(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_bom  = glob_sum( hfx_bom(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_bog  = glob_sum( hfx_bog(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_dif  = glob_sum( hfx_dif(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_opw  = glob_sum( hfx_opw(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_out  = glob_sum( hfx_out(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
+      zbg_hfx_in   = glob_sum(  hfx_in(:,:) * e12t(:,:) * tms(:,:) ) ! [in W]
     
       ! --------------------------------------------- !
       ! 2 - Trends due to forcing and ice growth/melt !
       ! --------------------------------------------- !
-      z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * area(:,:) * tms(:,:) ) ! volume fluxes
-      z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * area(:,:) * tms(:,:) ) ! salt fluxes
+      z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * e12t(:,:) * tms(:,:) ) ! volume fluxes
+      z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * e12t(:,:) * tms(:,:) ) ! salt fluxes
       z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + &
-                          &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * area(:,:) * tms(:,:) ) ! volume fluxes
+                          &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * e12t(:,:) * tms(:,:) ) ! volume fluxes
       !
       frc_vol  = frc_vol  + z_frc_vol  * rdt_ice
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 5064)
@@ -191,5 +191,5 @@
          CALL prt_ctl(tab2d_1=delta_i   , clinfo1=' lim_dyn  : delta_i   :')
          CALL prt_ctl(tab2d_1=strength  , clinfo1=' lim_dyn  : strength  :')
-         CALL prt_ctl(tab2d_1=area      , clinfo1=' lim_dyn  : cell area :')
+         CALL prt_ctl(tab2d_1=e12t      , clinfo1=' lim_dyn  : cell area :')
          CALL prt_ctl(tab2d_1=at_i      , clinfo1=' lim_dyn  : at_i      :')
          CALL prt_ctl(tab2d_1=vt_i      , clinfo1=' lim_dyn  : vt_i      :')
@@ -274,10 +274,10 @@
       rhoco  = rau0  * cw
       !
-      !  Diffusion coefficients.
+      !  Diffusion coefficients
       zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
       IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
            
       za00 = ahi0 / ( 1.e05_wp )              ! 1.e05 = 100km = max grid space at 60° latitude in orca2
-                                              !                 (i.e. 60° = min latitude for ice cover)  
+                                              !                    (60° = min latitude for ice cover)  
       DO jj = 1, jpj
          DO ji = 1, jpi
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 5064)
@@ -345,8 +345,7 @@
                    ! Snow energy of melting
                    e_s(ji,jj,jk,jl) = zswitch(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
-                   ! Change dimensions
-                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
-                   ! Multiply by volume, so that heat content in Joules
-                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s
+
+                   ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
+                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) / nlay_s
                END DO ! ji
             END DO ! jj
@@ -368,9 +367,6 @@
                       -   rcp     * ( ztmelts - rtt ) )
 
-                   ! Correct dimensions to avoid big values
-                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
-
-                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J
-                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i
+                   ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
+                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) / nlay_i
                END DO ! ji
             END DO ! jj
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 5064)
@@ -339,5 +339,5 @@
             !-----------------------------------------------------------------------------!
             wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean
-            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * unit_fac / area(ji,jj) * r1_rdtice  ! heat sink for ocean (<0, W.m-2)
+            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean (<0, W.m-2)
 
          END DO
@@ -381,5 +381,5 @@
          CALL prt_ctl_info(' - Cell values : ')
          CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
-         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_me  : cell area :')
+         CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_me  : cell area :')
          CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me  : at_i      :')
          CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me  : vt_i      :')
@@ -1093,5 +1093,5 @@
                &                                + rhosn*vsrft(ji,jj)*(1.0-fsnowrft)
 
-            ! in 1e-9 Joules (same as e_s)
+            ! in J/m2 (same as e_s)
             esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-fsnowrdg)         &   !rafting included
                &                                - esrft(ji,jj)*(1.0-fsnowrft)          
@@ -1133,12 +1133,5 @@
                hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux 
 
-               ! Correct dimensions to avoid big values
-               ersw(ji,jj,jk)   = ersw(ji,jj,jk) / unit_fac
-
-               ! Mutliply by ice volume, and divide by number of layers to get heat content in 1.e9 J
-               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean 
-               !! MV HC 2014
-               ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj)
-
+               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
                erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk)
 
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 5064)
@@ -129,7 +129,7 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               rswitch             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes
+               rswitch           = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i(ji,jj,jl) + epsi10 ) )     !0 if no ice and 1 if yes
                ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi10 ) * rswitch
-               rswitch             = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes
+               rswitch           = 1.0 - MAX( 0.0, SIGN( 1.0, - a_i_b(ji,jj,jl) + epsi10) ) !0 if no ice and 1 if yes
                zht_i_b(ji,jj,jl) = v_i_b(ji,jj,jl) / MAX( a_i_b(ji,jj,jl), epsi10 ) * rswitch
                IF( a_i(ji,jj,jl) > epsi10 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_b(ji,jj,jl) 
@@ -167,7 +167,4 @@
       !-----------------------------------------------------------------------------------------------
       !- 4.1 Compute category boundaries
-      ! Tricky trick see limitd_me.F90
-      ! will be soon removed, CT
-      ! hi_max(kubnd) = 99.
       zhbnew(:,:,:) = 0._wp
 
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90	(revision 5064)
@@ -114,6 +114,4 @@
       CALL lbc_lnk( tmf(:,:), 'F', 1. )           ! lateral boundary conditions
 
-      !                           !==  unmasked and masked area of T-grid cell
-      area(:,:) = e1t(:,:) * e2t(:,:)
       !
    END SUBROUTINE lim_msh
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 5064)
@@ -116,15 +116,15 @@
       CHARACTER (len=50) ::   charout
       REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         !
-      REAL(wp) ::   za, zstms, zmask   ! local scalars
-      REAL(wp) ::   zc1, zc2, zc3             ! ice mass
-
-      REAL(wp) ::   dtevp              ! time step for subcycling
-      REAL(wp) ::   dtotel, ecc2, ecci ! square of yield ellipse eccenticity
-      REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars
-      REAL(wp) ::   zu_ice2, zv_ice1   !
-      REAL(wp) ::   zddc, zdtc         ! delta on corners and on centre
-      REAL(wp) ::   zdst               ! shear at the center of the grid point
-      REAL(wp) ::   zdsshx, zdsshy     ! term for the gradient of ocean surface
-      REAL(wp) ::   sigma1, sigma2     ! internal ice stress
+      REAL(wp) ::   za, zstms          ! local scalars
+      REAL(wp) ::   zc1, zc2, zc3      ! ice mass
+
+      REAL(wp) ::   dtevp , z1_dtevp              ! time step for subcycling
+      REAL(wp) ::   dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity
+      REAL(wp) ::   z0, zr, zcca, zccb            ! temporary scalars
+      REAL(wp) ::   zu_ice2, zv_ice1              !
+      REAL(wp) ::   zddc, zdtc                    ! delta on corners and on centre
+      REAL(wp) ::   zdst                          ! shear at the center of the grid point
+      REAL(wp) ::   zdsshx, zdsshy                ! term for the gradient of ocean surface
+      REAL(wp) ::   sigma1, sigma2                ! internal ice stress
 
       REAL(wp) ::   zresm         ! Maximal error on ice velocity
@@ -137,6 +137,6 @@
       REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points
       REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays
-      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce1, v_oce1   ! ocean u/v component on U points                           
-      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2, v_oce2   ! ocean u/v component on V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1           ! ocean u/v component on U points                           
+      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2           ! ocean u/v component on V points
       REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point
       REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses
@@ -156,5 +156,5 @@
 
       CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
-      CALL wrk_alloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                )
+      CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1                )
       CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
       CALL wrk_alloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
@@ -228,8 +228,6 @@
       !  zcorl2: Coriolis parameter on V-points                            
       !  (ztagnx,ztagny): wind stress on U/V points                       
-      !  u_oce1: ocean u component on u points                           
       !  v_oce1: ocean v component on u points                          
       !  u_oce2: ocean u component on v points                         
-      !  v_oce2: ocean v component on v points                        
 
       IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
@@ -266,22 +264,19 @@
 
             ! Mass, coriolis coeff. and currents
-            zmass1(ji,jj) = ( zt12*zc1 + zt11*zc2 ) / (zt11+zt12+zepsi)
-            zmass2(ji,jj) = ( zt22*zc1 + zt21*zc3 ) / (zt21+zt22+zepsi)
-            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) )   &
+            zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi )
+            zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi )
+            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) )   &
                &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi )
-            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) )   &
+            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) )   &
                &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi )
             !
-            u_oce1(ji,jj)  = u_oce(ji,jj)
-            v_oce2(ji,jj)  = v_oce(ji,jj)
-
             ! Ocean has no slip boundary condition
-            v_oce1(ji,jj)  = 0.5*( (v_oce(ji,jj)+v_oce(ji,jj-1))*e1t(ji,jj)    &
-               &                 +(v_oce(ji+1,jj)+v_oce(ji+1,jj-1))*e1t(ji+1,jj)) &
-               &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)  
-
-            u_oce2(ji,jj)  = 0.5*((u_oce(ji,jj)+u_oce(ji-1,jj))*e2t(ji,jj)     &
-               &                 +(u_oce(ji,jj+1)+u_oce(ji-1,jj+1))*e2t(ji,jj+1)) &
-               &                / (e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
+            v_oce1(ji,jj)  = 0.5 * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji,jj)      &
+               &                   + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) )  &
+               &                   / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)  
+
+            u_oce2(ji,jj)  = 0.5 * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj)      &
+               &                   + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) )  &
+               &                   / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)
 
             ! Wind stress at U,V-point
@@ -316,4 +311,6 @@
       dtotel = dtevp / ( 2._wp * telast )
 #endif
+      z1_dtotel = 1._wp / ( 1._wp + dtotel )
+      z1_dtevp  = 1._wp / dtevp
       !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
       ecc2 = ecc * ecc
@@ -334,5 +331,5 @@
 
          DO jj = k_j1+1, k_jpj-1
-            DO ji = fs_2, jpim1   !RB bug no vect opt due to tmi
+            DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to tmi
 
                !  
@@ -357,28 +354,27 @@
                divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
-                  &            ) / area(ji,jj)
+                  &            ) * r1_e12t(ji,jj)
 
                zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
                   &         - ( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
-                  &         ) / area(ji,jj)
+                  &         ) * r1_e12t(ji,jj)
 
                !
                zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
                   &         + ( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-                  &         ) / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2._wp - tmf(ji,jj) )   &
+                  &         ) * r1_e12f(ji,jj) * ( 2._wp - tmf(ji,jj) )   &
                   &         * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1)
 
 
-               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji,jj)   + v_ice(ji,jj-1)   ) * e1t(ji+1,jj)   &
-                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji,jj) )   &
+               v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
+                  &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   &
                   &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj) 
 
-               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj)   + u_ice(ji-1,jj)   ) * e2t(ji,jj+1)   &
-                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj) )   &
+               u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
+                  &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   &
                   &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)
-
-            END DO
-         END DO
-         CALL lbc_lnk( v_ice1, 'U', -1. )   ;   CALL lbc_lnk( u_ice2, 'V', -1. )      ! lateral boundary cond.
+            END DO
+         END DO
+         CALL lbc_lnk( v_ice1  , 'U', -1. )   ;   CALL lbc_lnk( u_ice2  , 'V', -1. )      ! lateral boundary cond.
 
          DO jj = k_j1+1, k_jpj-1
@@ -386,55 +382,35 @@
 
                !- Calculate Delta at centre of grid cells
-               zdst      = ( e2u(ji, jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   &
-                  &        + e1v(ji, jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   &
-                  &        ) / area(ji,jj)
-
-               delta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) * usecc2 )  
+               zdst          = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   &
+                  &            + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   &
+                  &            ) * r1_e12t(ji,jj)
+
+               delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )  
                delta_i(ji,jj) = delta + creepl
-               !-Calculate stress tensor components zs1 and zs2 
-               !-at centre of grid cells (see section 3.5 of CICE user's guide).
-               zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( divu_i(ji,jj) / delta_i(ji,jj) - delta / delta_i(ji,jj) )  &
-                  &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
-               zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / delta_i(ji,jj) * zpresh(ji,jj) ) )  &
-                  &         / ( 1._wp + dtotel )
-
-            END DO
-         END DO
-
-         CALL lbc_lnk( zs1(:,:), 'T', 1. )
-         CALL lbc_lnk( zs2(:,:), 'T', 1. )
-
-         DO jj = k_j1+1, k_jpj-1
-            DO ji = fs_2, fs_jpim1
+
                !- Calculate Delta on corners
-               zddc  =      ( ( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
-                  &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
-                  &           )*e1f(ji,jj)*e1f(ji,jj)                       &
-                  &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                &
-                  &            -u_ice2(ji,jj)/e2v(ji,jj)                    &
-                  &           )*e2f(ji,jj)*e2f(ji,jj)                       &
-                  &         )                                               &
-                  &        / ( e1f(ji,jj) * e2f(ji,jj) )
-
-               zdtc  =      (-( v_ice1(ji,jj+1)/e1u(ji,jj+1)                &
-                  &            -v_ice1(ji,jj)/e1u(ji,jj)                    &
-                  &           )*e1f(ji,jj)*e1f(ji,jj)                       &
-                  &          +( u_ice2(ji+1,jj)/e2v(ji+1,jj)                &
-                  &            -u_ice2(ji,jj)/e2v(ji,jj)                    &
-                  &           )*e2f(ji,jj)*e2f(ji,jj)                       &
-                  &         )                                               &
-                  &        / ( e1f(ji,jj) * e2f(ji,jj) )
-
-               zddc = SQRT(zddc**2+(zdtc**2+zds(ji,jj)**2)*usecc2) + creepl
-
-               !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
-               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
-                  &          ( ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj) ) )  &
-                  &          / ( 1.0 + dtotel ) 
-
-            END DO ! ji
-         END DO ! jj
-
-         CALL lbc_lnk( zs12(:,:), 'F', 1. )
+               zddc  = (  ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
+                  &     + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
+                  &    ) * r1_e12f(ji,jj)
+
+               zdtc  = (- ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
+                  &     + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
+                  &    ) * r1_e12f(ji,jj)
+
+               zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + creepl
+
+               !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide).
+               zs1(ji,jj)  = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj)   * zpresh(ji,jj)   &
+                  &          ) * z1_dtotel
+               zs2(ji,jj)  = ( zs2 (ji,jj) + dtotel *         ecci * zdt(ji,jj) / delta_i(ji,jj)   * zpresh(ji,jj)   &
+                  &          ) * z1_dtotel
+               !-Calculate stress tensor component zs12 at corners
+               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *         ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj)  &
+                  &          ) * z1_dtotel 
+
+            END DO
+         END DO
+         CALL lbc_lnk( zs1 , 'T', 1. )   ;   CALL lbc_lnk( zs2, 'T', 1. )
+         CALL lbc_lnk( zs12, 'F', 1. )
 
          ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)
@@ -442,14 +418,13 @@
             DO ji = fs_2, fs_jpim1
                !- contribution of zs1, zs2 and zs12 to zf1
-               zf1(ji,jj) = 0.5*( (zs1(ji+1,jj)-zs1(ji,jj))*e2u(ji,jj) &
-                  &              +(zs2(ji+1,jj)*e2t(ji+1,jj)**2-zs2(ji,jj)*e2t(ji,jj)**2)/e2u(ji,jj) &
-                  &              +2.0*(zs12(ji,jj)*e1f(ji,jj)**2-zs12(ji,jj-1)*e1f(ji,jj-1)**2)/e1u(ji,jj) &
-                  &             ) / ( e1u(ji,jj)*e2u(ji,jj) )
+               zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  &
+                  &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) / e2u(ji,jj)          &
+                  &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) / e1u(ji,jj)  &
+                  &                ) * r1_e12u(ji,jj)
                ! contribution of zs1, zs2 and zs12 to zf2
-               zf2(ji,jj) = 0.5*( (zs1(ji,jj+1)-zs1(ji,jj))*e1v(ji,jj) &
-                  &              -(zs2(ji,jj+1)*e1t(ji,jj+1)**2 - zs2(ji,jj)*e1t(ji,jj)**2)/e1v(ji,jj) &
-                  &              + 2.0*(zs12(ji,jj)*e2f(ji,jj)**2 -    &
-                  zs12(ji-1,jj)*e2f(ji-1,jj)**2)/e2v(ji,jj) &
-                  &             ) / ( e1v(ji,jj)*e2v(ji,jj) )
+               zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  &
+                  &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) / e1v(ji,jj)          &
+                  &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) / e2v(ji,jj)  &
+                  &               )  * r1_e12v(ji,jj)
             END DO
          END DO
@@ -463,19 +438,17 @@
             DO jj = k_j1+1, k_jpj-1
                DO ji = fs_2, fs_jpim1
-                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)
-                  z0           = zmass1(ji,jj)/dtevp
+                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * tmu(ji,jj)
+                  z0           = zmass1(ji,jj) * z1_dtevp
 
                   ! SB modif because ocean has no slip boundary condition
-                  zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)         &
-                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   &
-                     &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
-                  za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
-                     (zv_ice1-v_oce1(ji,jj))**2) * (1.0-zfrld1(ji,jj))
-                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
-                     za*(u_oce1(ji,jj))
-                  zcca         = z0+za
+                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji  ,jj)     &
+                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
+                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)
+                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
+                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
+                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
+                  zcca         = z0 + za
                   zccb         = zcorl1(ji,jj)
-                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+zepsi)*zmask 
-
+                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
                END DO
             END DO
@@ -492,18 +465,16 @@
                DO ji = fs_2, fs_jpim1
 
-                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)
-                  z0           = zmass2(ji,jj)/dtevp
+                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * tmv(ji,jj)
+                  z0           = zmass2(ji,jj) * z1_dtevp
                   ! SB modif because ocean has no slip boundary condition
-                  zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)     &
-                     &                 + (u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   &
-                     &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)
-                  za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + & 
-                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
-                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
-                     za2ct(ji,jj) + za*(v_oce2(ji,jj))
-                  zcca         = z0+za
+                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
+                     &                 + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
+                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)
+                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  & 
+                     &                         ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) )
+                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
+                  zcca         = z0 + za
                   zccb         = zcorl2(ji,jj)
-                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+zepsi)*zmask
-
+                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
                END DO
             END DO
@@ -520,19 +491,17 @@
             DO jj = k_j1+1, k_jpj-1
                DO ji = fs_2, fs_jpim1
-                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass2(ji,jj))))*tmv(ji,jj)
-                  z0           = zmass2(ji,jj)/dtevp
+                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * tmv(ji,jj)
+                  z0           = zmass2(ji,jj) * z1_dtevp
                   ! SB modif because ocean has no slip boundary condition
-                  zu_ice2       = 0.5*( (u_ice(ji,jj)+u_ice(ji-1,jj))*e2t(ji,jj)      &
-                     &                 +(u_ice(ji,jj+1)+u_ice(ji-1,jj+1))*e2t(ji,jj+1))   &
-                     &               /(e2t(ji,jj+1)+e2t(ji,jj)) * tmv(ji,jj)   
-
-                  za           = rhoco*SQRT((zu_ice2-u_oce2(ji,jj))**2 + &
-                     (v_ice(ji,jj)-v_oce2(ji,jj))**2)*(1.0-zfrld2(ji,jj))
-                  zr           = z0*v_ice(ji,jj) + zf2(ji,jj) + &
-                     za2ct(ji,jj) + za*(v_oce2(ji,jj))
-                  zcca         = z0+za
+                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
+                     &                  +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
+                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * tmv(ji,jj)   
+
+                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  &
+                     &                         ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) )
+                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
+                  zcca         = z0 + za
                   zccb         = zcorl2(ji,jj)
-                  v_ice(ji,jj) = (zr-zccb*zu_ice2)/(zcca+zepsi)*zmask
-
+                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
                END DO
             END DO
@@ -548,19 +517,18 @@
             DO jj = k_j1+1, k_jpj-1
                DO ji = fs_2, fs_jpim1
-                  zmask        = (1.0-MAX(0._wp,SIGN(1._wp,-zmass1(ji,jj))))*tmu(ji,jj)
-                  z0           = zmass1(ji,jj)/dtevp
-                  zv_ice1       = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji,jj)      &
-                     &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji+1,jj))   &
-                     &               /(e1t(ji+1,jj)+e1t(ji,jj)) * tmu(ji,jj)
-
-                  za           = rhoco*SQRT((u_ice(ji,jj)-u_oce1(ji,jj))**2 + &
-                     (zv_ice1-v_oce1(ji,jj))**2)*(1.0-zfrld1(ji,jj))
-                  zr           = z0*u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + &
-                     za*(u_oce1(ji,jj))
-                  zcca         = z0+za
+                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * tmu(ji,jj)
+                  z0           = zmass1(ji,jj) * z1_dtevp
+                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji,jj)       &
+                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
+                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * tmu(ji,jj)
+
+                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
+                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
+                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
+                  zcca         = z0 + za
                   zccb         = zcorl1(ji,jj)
-                  u_ice(ji,jj) = (zr+zccb*zv_ice1)/(zcca+zepsi)*zmask 
-               END DO ! ji
-            END DO ! jj
+                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
+               END DO
+            END DO
 
             CALL lbc_lnk( u_ice(:,:), 'U', -1. )
@@ -577,8 +545,7 @@
             !---  Convergence test.
             DO jj = k_j1+1 , k_jpj-1
-               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) ,           &
-                  ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
-            END DO
-            zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
+               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
+            END DO
+            zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) )
             IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
          ENDIF
@@ -637,43 +604,27 @@
             IF ( vt_i(ji,jj) <= zvmin ) THEN
 
-               divu_i(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      &
-                  &             -e2u(ji-1,jj)*u_ice(ji-1,jj)                  &
-                  &             +e1v(ji,jj)*v_ice(ji,jj)                      &
-                  &             -e1v(ji,jj-1)*v_ice(ji,jj-1)                  &
-                  &            )                                              &
-                  &            / area(ji,jj)
-
-               zdt(ji,jj) = ( ( u_ice(ji,jj)/e2u(ji,jj)                    &
-                  &            -u_ice(ji-1,jj)/e2u(ji-1,jj)                &
-                  &           )*e2t(ji,jj)*e2t(ji,jj)                      &
-                  &          -( v_ice(ji,jj)/e1v(ji,jj)                    &
-                  &            -v_ice(ji,jj-1)/e1v(ji,jj-1)                &
-                  &           )*e1t(ji,jj)*e1t(ji,jj)                      &
-                  &         )                                              &
-                  &        / area(ji,jj)
+               divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj  ) * u_ice(ji-1,jj  )   &
+                  &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji  ,jj-1) * v_ice(ji  ,jj-1)   &
+                  &            ) * r1_e12t(ji,jj)
+
+               zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  &
+                  &          -( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  &
+                  &         ) * r1_e12t(ji,jj)
                !
                ! SB modif because ocean has no slip boundary condition 
-               zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1)              &
-                  &           - u_ice(ji,jj)   / e1u(ji,jj) )              &
-                  &           * e1f(ji,jj) * e1f(ji,jj)                    &
-                  &          + ( v_ice(ji+1,jj) / e2v(ji+1,jj)             &
-                  &            - v_ice(ji,jj)  / e2v(ji,jj) )              &
-                  &           * e2f(ji,jj) * e2f(ji,jj) )                  &
-                  &        / ( e1f(ji,jj) * e2f(ji,jj) ) * ( 2.0 - tmf(ji,jj) ) &
-                  &        * tmi(ji,jj) * tmi(ji,jj+1)                     &
-                  &        * tmi(ji+1,jj) * tmi(ji+1,jj+1)
-
-               zdst = (  e2u( ji  , jj   ) * v_ice1(ji  ,jj  )    &
-                  &           - e2u( ji-1, jj   ) * v_ice1(ji-1,jj  )    &
-                  &           + e1v( ji  , jj   ) * u_ice2(ji  ,jj  )    &
-                  &           - e1v( ji  , jj-1 ) * u_ice2(ji  ,jj-1)  ) / area(ji,jj)
-
-               delta = SQRT( divu_i(ji,jj)*divu_i(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )  
+               zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
+                  &          +( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
+                  &         ) * r1_e12f(ji,jj) * ( 2.0 - tmf(ji,jj) )                                     &
+                  &         * tmi(ji,jj) * tmi(ji,jj+1) * tmi(ji+1,jj) * tmi(ji+1,jj+1)
+
+               zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )    &
+                  &   + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1) ) * r1_e12t(ji,jj)
+
+               delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )  
                delta_i(ji,jj) = delta + creepl
             
             ENDIF
-
-         END DO !jj
-      END DO !ji
+         END DO
+      END DO
       !
       !------------------------------------------------------------------------------!
@@ -684,8 +635,6 @@
       DO jj = k_j1+1, k_jpj-1
          DO ji = fs_2, fs_jpim1
-            zdst= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &   
-               &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &   
-               &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &   
-               &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj) 
+            zdst           = (  e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)  &   
+               &              + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj) 
             shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst )
          END DO
@@ -741,5 +690,5 @@
       !
       CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
-      CALL wrk_dealloc( jpi,jpj, u_oce1, u_oce2, u_ice2, v_oce1 , v_oce2, v_ice1                )
+      CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1                )
       CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
       CALL wrk_dealloc( jpi,jpj, zdt   , zds   , zs1   , zs2   , zs12   , zresr , zpice                 )
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 5064)
@@ -26,5 +26,5 @@
    USE phycst           ! physical constants
    USE dom_oce          ! ocean domain
-   USE dom_ice,    ONLY : tms, area
+   USE dom_ice,    ONLY : tms
    USE ice              ! LIM sea-ice variables
    USE sbc_ice          ! Surface boundary condition: sea-ice fields
@@ -99,11 +99,9 @@
       !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
       !!              These refs are now obsolete since everything has been revised
-      !!              The ref should be Rousset et al., 2015?
+      !!              The ref should be Rousset et al., 2015
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt                                   ! number of iteration
-      !
       INTEGER  ::   ji, jj, jl, jk                                  ! dummy loop indices
-      !
-      REAL(wp) ::   zemp                                            !  local scalars
+      REAL(wp) ::   zemp                                            ! local scalars
       REAL(wp) ::   zf_mass                                         ! Heat flux associated with mass exchange ice->ocean (W.m-2)
       REAL(wp) ::   zfcm1                                           ! New solar flux received by the ocean
@@ -321,6 +319,4 @@
       !! ** input   : Namelist namicedia
       !!-------------------------------------------------------------------
-      REAL(wp) :: zsum, zarea
-      !
       INTEGER  ::   ji, jj, jk                       ! dummy loop indices
       REAL(wp) ::   zcoefu, zcoefv, zcoeff          ! local scalar
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5064)
@@ -109,5 +109,5 @@
       ! 1.2) Heat content    
       !--------------------
-      ! Change the units of heat content; from global units to J.m3
+      ! Change the units of heat content; from J/m2 to J/m3
       DO jl = 1, jpl
          DO jk = 1, nlay_i
@@ -115,9 +115,7 @@
                DO ji = 1, jpi
                   !0 if no ice and 1 if yes
-                  rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 )  )
+                  rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 )  )
                   !Energy of melting q(S,T) [J.m-3]
-                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi10 ) ) * REAL( nlay_i )
-                  !convert units ! very important that this line is here        
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac 
+                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL( nlay_i )
                END DO
             END DO
@@ -127,9 +125,7 @@
                DO ji = 1, jpi
                   !0 if no ice and 1 if yes
-                  rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi10 )  )
+                  rswitch = 1.0 - MAX(  0.0 , SIGN( 1.0 , - v_s(ji,jj,jl) + epsi20 )  )
                   !Energy of melting q(S,T) [J.m-3]
-                  e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL( nlay_s )
-                  !convert units ! very important that this line is here
-                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac 
+                  e_s(ji,jj,jk,jl) = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL( nlay_s )
                END DO
             END DO
@@ -453,8 +449,8 @@
       ! Ice heat content              
       !------------------------
-      ! Enthalpies are global variables we have to readjust the units (heat content in Joules)
+      ! Enthalpies are global variables we have to readjust the units (heat content in J/m2)
       DO jl = 1, jpl
          DO jk = 1, nlay_i
-            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_i(:,:,jl) / ( unit_fac * REAL( nlay_i ) )
+            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) / REAL( nlay_i )
          END DO
       END DO
@@ -463,8 +459,8 @@
       ! Snow heat content              
       !------------------------
-      ! Enthalpies are global variables we have to readjust the units (heat content in Joules)
+      ! Enthalpies are global variables we have to readjust the units (heat content in J/m2)
       DO jl = 1, jpl
          DO jk = 1, nlay_s
-            e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * area(:,:) * a_i(:,:,jl) * ht_s(:,:,jl) / ( unit_fac * REAL( nlay_s ) )
+            e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) / REAL( nlay_s )
          END DO
       END DO
@@ -489,5 +485,5 @@
          CALL prt_ctl_info(' - Cell values : ')
          CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
-         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_thd  : cell area :')
+         CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_thd  : cell area :')
          CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_thd  : at_i      :')
          CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_thd  : vt_i      :')
@@ -520,9 +516,11 @@
       CALL wrk_dealloc( jpi, jpj, zqsr, zqns )
 
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !------------------------------------------------------------------------------|
       !  1) Transport of ice between thickness categories.                           |
       !------------------------------------------------------------------------------|
-      ! Given thermodynamic growth rates, transport ice between
-      ! thickness categories.
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
+      ! Given thermodynamic growth rates, transport ice between thickness categories.
       IF( jpl > 1 )   CALL lim_itd_th_rem( 1, jpl, kt )
       !
@@ -530,4 +528,6 @@
       CALL lim_var_agg(1)
 
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th_rem', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !------------------------------------------------------------------------------|
       !  3) Add frazil ice growing in leads.
@@ -536,9 +536,12 @@
       CALL lim_var_glo2eqv    ! only for info
       
+      ! conservation test
+      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
       IF(ln_ctl) THEN   ! Control print
          CALL prt_ctl_info(' ')
          CALL prt_ctl_info(' - Cell values : ')
          CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
-         CALL prt_ctl(tab2d_1=area , clinfo1=' lim_itd_th  : cell area :')
+         CALL prt_ctl(tab2d_1=e12t , clinfo1=' lim_itd_th  : cell area :')
          CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_th  : at_i      :')
          CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_th  : vt_i      :')
@@ -568,7 +571,4 @@
       ENDIF
       !
-      ! conservation test
-      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
       IF( nn_timing == 1 )  CALL timing_stop('limthd')
 
@@ -614,16 +614,25 @@
       !!                          ( dA = A/2h dh )
       !!-----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop
-      INTEGER             ::   ji            ! dummy loop indices
-      REAL(wp)            ::   zhi_bef       ! ice thickness before thermo
-      REAL(wp)            ::   zdh_mel       ! net melting
+      INTEGER, INTENT(in) ::   kideb, kiut        ! bounds for the spatial loop
+      INTEGER             ::   ji                 ! dummy loop indices
+      REAL(wp)            ::   zhi_bef            ! ice thickness before thermo
+      REAL(wp)            ::   zdh_mel, zda_mel   ! net melting
+      REAL(wp)            ::   zv                 ! ice volume 
 
       DO ji = kideb, kiut
-         zhi_bef = ht_i_1d(ji) - ( dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) )
-         zdh_mel = MIN( 0._wp,     dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) )
-         IF( zdh_mel < 0._wp )   &
-            &   a_i_1d(ji) = MAX( 0._wp, a_i_1d(ji) + a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) ) ) 
+         zdh_mel = MIN( 0._wp, dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji) )
+         IF( zdh_mel < 0._wp )  THEN
+            zv         = a_i_1d(ji) * ht_i_1d(ji)
+            ! lateral melting = concentration change
+            zhi_bef     = ht_i_1d(ji) - zdh_mel
+            zda_mel     =  a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi10 ) )
+            a_i_1d(ji)  = MAX( 0._wp, a_i_1d(ji) + zda_mel ) 
+            ! adjust thickness
+            rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - a_i_1d(ji) + epsi20 ) )
+            ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 )
+            ! adjust thickness
+            at_i_1d(ji) = a_i_1d(ji)
+         END IF
       END DO
-      at_i_1d(:) = a_i_1d(:)
       
    END SUBROUTINE lim_thd_lam
@@ -665,5 +674,5 @@
       IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN 
          nn_monocat = 0
-         WRITE(numout, *) ' nn_monocat must be 0 in multi-category case '
+         IF(lwp) WRITE(numout, *) '   nn_monocat must be 0 in multi-category case '
       ENDIF
 
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5064)
@@ -157,5 +157,4 @@
       REAL(wp)                            :: zh_thres      ! thickness thres. for G(h) computation
       REAL(wp)                            :: zhe           ! dummy factor
-      REAL(wp)                            :: zswitch       ! dummy switch
       REAL(wp)                            :: zkimean       ! mean sea ice thermal conductivity
       REAL(wp)                            :: zfac          ! dummy factor
@@ -363,11 +362,7 @@
          ! Fichefet and Morales Maqueda, JGR 1997
 
-         zghe(:) = 0._wp
+         zghe(:) = 1._wp
 
          SELECT CASE ( nn_monocat )
-
-         CASE (0,2,4)
-
-            zghe(kideb:kiut) = 1._wp
 
          CASE (1,3) ! LIM3
@@ -388,9 +383,9 @@
 
                ! G(he)
-               zswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if >
-               zghe(ji) = ( 1.0 - zswitch ) + zswitch * ( 0.5 + 0.5 * LOG( 2.*zhe / zepsilon ) )
+               rswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if >
+               zghe(ji) = ( 1._wp - rswitch ) + rswitch * 0.5_wp * ( 1._wp + LOG( 2.*zhe / zepsilon ) )
    
                ! Impose G(he) < 2.
-               zghe(ji) = MIN( zghe(ji), 2.0 )
+               zghe(ji) = MIN( zghe(ji), 2._wp )
 
             END DO
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 5064)
@@ -132,6 +132,6 @@
       DO jk1 = 1, nlay_i
          DO ji = kideb, kiut
-            rswitch      = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi10 ) ) 
-            qnew(ji,jk1) = rswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi10 )
+            rswitch      = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi20 ) ) 
+            qnew(ji,jk1) = rswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 )
          ENDDO
       ENDDO
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 5064)
@@ -130,8 +130,6 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  rswitch          = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi10 )  )   !0 if no ice and 1 if yes
-                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) &
-                      &   / ( area(ji,jj) * MAX( v_i(ji,jj,jl) ,  epsi10 ) ) * REAL( nlay_i, wp )
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac
+                  rswitch          = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) + epsi20 )  )   !0 if no ice
+                  e_i(ji,jj,jk,jl) = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl), epsi20 ) * REAL( nlay_i, wp )
                END DO
             END DO
@@ -526,6 +524,6 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! heat content in Joules
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ,wp ) * unit_fac ) 
+                  ! heat content in J/m2
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) / REAL( nlay_i ,wp ) 
                END DO
             END DO
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5064)
@@ -95,5 +95,5 @@
       ENDIF
 
-      zsm(:,:) = area(:,:)
+      zsm(:,:) = e12t(:,:)
       
       !                             !-------------------------------------!
@@ -153,13 +153,15 @@
          ! transported fields                                        
          !-------------------------
-         zs0ow(:,:,1) = ato_i(:,:) * area(:,:)              ! Open water area 
-         DO jl = 1, jpl
-            zs0sn (:,:,jl)   = v_s  (:,:,jl) * area(:,:)    ! Snow volume
-            zs0ice(:,:,jl)   = v_i  (:,:,jl) * area(:,:)    ! Ice  volume
-            zs0a  (:,:,jl)   = a_i  (:,:,jl) * area(:,:)    ! Ice area
-            zs0sm (:,:,jl)   = smv_i(:,:,jl) * area(:,:)    ! Salt content
-            zs0oi (:,:,jl)   = oa_i (:,:,jl) * area(:,:)    ! Age content
-            zs0c0 (:,:,jl)   = e_s  (:,:,1,jl)              ! Snow heat content
-            zs0e  (:,:,:,jl) = e_i  (:,:,:,jl)              ! Ice  heat content
+         zs0ow(:,:,1) = ato_i(:,:) * e12t(:,:)              ! Open water area 
+         DO jl = 1, jpl
+            zs0sn (:,:,jl)   = v_s  (:,:,jl) * e12t(:,:)    ! Snow volume
+            zs0ice(:,:,jl)   = v_i  (:,:,jl) * e12t(:,:)    ! Ice  volume
+            zs0a  (:,:,jl)   = a_i  (:,:,jl) * e12t(:,:)    ! Ice area
+            zs0sm (:,:,jl)   = smv_i(:,:,jl) * e12t(:,:)    ! Salt content
+            zs0oi (:,:,jl)   = oa_i (:,:,jl) * e12t(:,:)    ! Age content
+            zs0c0 (:,:,jl)   = e_s  (:,:,1,jl) * e12t(:,:)  ! Snow heat content
+            DO jk = 1, nlay_i
+               zs0e  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
+            END DO
          END DO
 
@@ -253,13 +255,15 @@
          ! Recover the properties from their contents
          !-------------------------------------------
-         ato_i(:,:) = zs0ow(:,:,1) / area(:,:)
-         DO jl = 1, jpl
-            v_i  (:,:,jl)   = zs0ice(:,:,jl) / area(:,:)
-            v_s  (:,:,jl)   = zs0sn (:,:,jl) / area(:,:)
-            smv_i(:,:,jl)   = zs0sm (:,:,jl) / area(:,:)
-            oa_i (:,:,jl)   = zs0oi (:,:,jl) / area(:,:)
-            a_i  (:,:,jl)   = zs0a  (:,:,jl) / area(:,:)
-            e_s  (:,:,1,jl) = zs0c0 (:,:,jl)              
-            e_i  (:,:,:,jl) = zs0e  (:,:,:,jl)           
+         ato_i(:,:) = zs0ow(:,:,1) * r1_e12t(:,:)
+         DO jl = 1, jpl
+            v_i  (:,:,jl)   = zs0ice(:,:,jl) * r1_e12t(:,:)
+            v_s  (:,:,jl)   = zs0sn (:,:,jl) * r1_e12t(:,:)
+            smv_i(:,:,jl)   = zs0sm (:,:,jl) * r1_e12t(:,:)
+            oa_i (:,:,jl)   = zs0oi (:,:,jl) * r1_e12t(:,:)
+            a_i  (:,:,jl)   = zs0a  (:,:,jl) * r1_e12t(:,:)
+            e_s  (:,:,1,jl) = zs0c0 (:,:,jl) * r1_e12t(:,:)
+            DO jk = 1, nlay_i
+               e_i  (:,:,jk,jl) = zs0e  (:,:,jk,jl) * r1_e12t(:,:)
+            END DO
          END DO
 
@@ -380,6 +384,6 @@
                      wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice
                      sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice 
-                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
-                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
+                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
+                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( SUM( e_i(ji,jj,1:nlay_i,jl) ) - zei ) * r1_rdtice ! W.m-2 <0
                   ENDIF
 
@@ -404,6 +408,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice
-               diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) / area(ji,jj) * unit_fac * r1_rdtice
+               diag_trp_ei(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) ) - zeiold(ji,jj) ) * r1_rdtice
+               diag_trp_es(ji,jj) = ( SUM( e_s(ji,jj,1:nlay_s,:) ) - zesold(ji,jj) ) * r1_rdtice
 
                diag_trp_vi(ji,jj) = SUM( v_i(ji,jj,:) - zviold(ji,jj,:) ) * r1_rdtice
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 5064)
@@ -16,4 +16,5 @@
    USE sbc_ice         ! Surface boundary condition: ice fields
    USE dom_ice
+   USE dom_oce
    USE phycst          ! physical constants
    USE ice
@@ -148,5 +149,5 @@
             diag_heat_dhc(ji,jj) = ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
                &                     SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
-               &                   ) * unit_fac * r1_rdtice / area(ji,jj)   
+               &                   ) * r1_rdtice   
          END DO
       END DO
@@ -159,5 +160,5 @@
          CALL prt_ctl_info(' - Cell values : ')
          CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
-         CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update1  : cell area   :')
+         CALL prt_ctl(tab2d_1=e12t       , clinfo1=' lim_update1  : cell area   :')
          CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update1  : at_i        :')
          CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update1  : vt_i        :')
@@ -183,8 +184,8 @@
             CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update1  : v_s         : ')
             CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update1  : v_s_b       : ')
-            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1        : ')
-            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)/1.0e9, clinfo1= ' lim_update1  : e_i1_b      : ')
-            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2        : ')
-            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)/1.0e9, clinfo1= ' lim_update1  : e_i2_b      : ')
+            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_i1        : ')
+            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_i1_b      : ')
+            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)      , clinfo1= ' lim_update1  : e_i2        : ')
+            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)      , clinfo1= ' lim_update1  : e_i2_b      : ')
             CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow      : ')
             CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update1  : e_snow_b    : ')
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 5064)
@@ -16,4 +16,5 @@
    USE sbc_ice         ! Surface boundary condition: ice fields
    USE dom_ice
+   USE dom_oce
    USE phycst          ! physical constants
    USE ice
@@ -190,5 +191,5 @@
                &                   ( SUM( e_i(ji,jj,1:nlay_i,:) - e_i_b(ji,jj,1:nlay_i,:) ) +  & 
                &                     SUM( e_s(ji,jj,1:nlay_s,:) - e_s_b(ji,jj,1:nlay_s,:) )    &
-               &                   ) * unit_fac * r1_rdtice / area(ji,jj)   
+               &                   ) * r1_rdtice   
          END DO
       END DO
@@ -203,5 +204,5 @@
          CALL prt_ctl_info(' - Cell values : ')
          CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
-         CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update2  : cell area   :')
+         CALL prt_ctl(tab2d_1=e12t       , clinfo1=' lim_update2  : cell area   :')
          CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update2  : at_i        :')
          CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update2  : vt_i        :')
@@ -227,8 +228,8 @@
             CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update2  : v_s         : ')
             CALL prt_ctl(tab2d_1=v_s_b      (:,:,jl)        , clinfo1= ' lim_update2  : v_s_b       : ')
-            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1        : ')
-            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)/1.0e9, clinfo1= ' lim_update2  : e_i1_b      : ')
-            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2        : ')
-            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)/1.0e9, clinfo1= ' lim_update2  : e_i2_b      : ')
+            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_i1        : ')
+            CALL prt_ctl(tab2d_1=e_i_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_i1_b      : ')
+            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)      , clinfo1= ' lim_update2  : e_i2        : ')
+            CALL prt_ctl(tab2d_1=e_i_b      (:,:,2,jl)      , clinfo1= ' lim_update2  : e_i2_b      : ')
             CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow      : ')
             CALL prt_ctl(tab2d_1=e_s_b      (:,:,1,jl)      , clinfo1= ' lim_update2  : e_snow_b    : ')
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 5064)
@@ -190,7 +190,6 @@
                DO ji = 1, jpi
                   !                                                              ! Energy of melting q(S,T) [J.m-3]
-                  rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes
-                  zq_i    = rswitch * e_i(ji,jj,jk,jl) / area(ji,jj) / MAX( v_i(ji,jj,jl) , epsi10 ) * REAL(nlay_i,wp) 
-                  zq_i    = zq_i * unit_fac                             ! convert units
+                  rswitch   = 1.0 - MAX( 0.0 , SIGN( 1.0 , - v_i(ji,jj,jl) + epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes
+                  zq_i    = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp) 
                   ztmelts = -tmut * s_i(ji,jj,jk,jl) + rtt              ! Ice layer melt temperature
                   !
@@ -216,7 +215,6 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi10 ) )     ! rswitch = 0 if no ice and 1 if yes
-                  zq_s  = rswitch * e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi10 ) ) * REAL(nlay_s,wp)
-                  zq_s  = zq_s * unit_fac                                    ! convert units
+                  rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - v_s(ji,jj,jl) + epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes
+                  zq_s  = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp)
                   !
                   t_s(ji,jj,jk,jl) = rtt + rswitch * ( - zfac1 * zq_s + zfac2 )
@@ -541,5 +539,5 @@
                   t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rtt * ( 1._wp - rswitch )
                   ! update exchanges with ocean
-                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
+                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
                END DO
             END DO
@@ -579,5 +577,5 @@
                wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
                wfx_snw(ji,jj)  = wfx_snw(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 ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0
+               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
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 5060)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 5064)
@@ -229,8 +229,6 @@
                ! Snow energy of melting
                e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
-               ! Change dimensions
-               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
-               ! Multiply by volume, so that heat content in 10^9 Joules
-               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s
+               ! Multiply by volume, so that heat content in J/m2
+               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) / nlay_s
             END DO
             DO jk = 1, nlay_i
@@ -241,8 +239,6 @@
                   +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) &
                   - rcp      * ( ztmelts - rtt ) )
-               ! Correct dimensions to avoid big values
-               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
-               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
-               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i
+               ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
+               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / nlay_i
             END DO
 
