Index: branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4744)
+++ branches/2014/dev_r4650_UKMO10_Tidally_Meaned_Diagnostics/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 5260)
@@ -22,5 +22,4 @@
    USE thd_ice        ! LIM thermodynamics
    USE dom_ice        ! LIM domain
-   USE par_ice        ! LIM parameters
    USE ice            ! LIM variables
    USE limtab         ! LIM 2D <==> 1D
@@ -29,5 +28,8 @@
    USE lib_mpp        ! MPP library
    USE wrk_nemo       ! work arrays
+   USE lbclnk         ! ocean lateral boundary conditions (or mpp link)
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+   USE limthd_ent
+   USE limvar
 
    IMPLICIT NONE
@@ -35,8 +37,4 @@
 
    PUBLIC lim_thd_lac     ! called by lim_thd
-
-   REAL(wp) ::   epsi10 = 1.e-10_wp   !
-   REAL(wp) ::   zzero  = 0._wp      !
-   REAL(wp) ::   zone   = 1._wp      !
 
    !!----------------------------------------------------------------------
@@ -71,16 +69,23 @@
       !!             - Computation of variation of ice volume and mass
       !!             - Computation of frldb after lateral accretion and 
-      !!               update ht_s_b, ht_i_b and tbif_1d(:,:)      
+      !!               update ht_s_1d, ht_i_1d and tbif_1d(:,:)      
       !!------------------------------------------------------------------------
-      INTEGER ::   ji,jj,jk,jl,jm   ! dummy loop indices
-      INTEGER ::   layer, nbpac     ! local integers 
-      INTEGER ::   ii, ij, iter   !   -       -
-      REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde  ! local scalars
+      INTEGER ::   ji,jj,jk,jl      ! dummy loop indices
+      INTEGER ::   nbpac            ! local integers 
+      INTEGER ::   ii, ij, iter     !   -       -
+      REAL(wp)  ::   ztmelts, zdv, zfrazb, zweight, zde                         ! local scalars
       REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      -
       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
       LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness
       CHARACTER (len = 15) :: fieldid
-      !
-      INTEGER , POINTER, DIMENSION(:) ::   zcatac      ! indexes of categories where new ice grows
+
+      REAL(wp) ::   zQm          ! enthalpy exchanged with the ocean (J/m2, >0 towards ocean)
+      REAL(wp) ::   zEi          ! sea ice specific enthalpy (J/kg)
+      REAL(wp) ::   zEw          ! seawater specific enthalpy (J/kg)
+      REAL(wp) ::   zfmdt        ! mass flux x time step (kg/m2, >0 towards ocean)
+     
+      REAL(wp) ::   zv_newfra
+  
+      INTEGER , POINTER, DIMENSION(:) ::   jcat        ! indexes of categories where new ice grows
       REAL(wp), POINTER, DIMENSION(:) ::   zswinew     ! switch for new ice or not
 
@@ -93,57 +98,30 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zdv_res     ! residual volume in case of excessive heat budget
       REAL(wp), POINTER, DIMENSION(:) ::   zda_res     ! residual area in case of excessive heat budget
-      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_ac    ! total ice fraction    
-      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_lev   ! total ice fraction for level ice only (type 1)   
-      REAL(wp), POINTER, DIMENSION(:) ::   zdh_frazb   ! accretion of frazil ice at the ice bottom
-      REAL(wp), POINTER, DIMENSION(:) ::   zvrel_ac    ! relative ice / frazil velocity (1D vector)
-
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zhice_old   ! previous ice thickness
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zdummy      ! dummy thickness of new ice 
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zdhicbot    ! thickness of new ice which is accreted vertically
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_old      ! old volume of ice in category jl
-      REAL(wp), POINTER, DIMENSION(:,:) ::   za_old      ! old area of ice in category jl
-      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_ac     ! 1-D version of a_i
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_ac     ! 1-D version of v_i
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_ac    ! 1-D version of oa_i
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_ac   ! 1-D version of smv_i
-
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_ac   !: 1-D version of e_i
-
-      REAL(wp), POINTER, DIMENSION(:) ::   zqbgow    ! heat budget of the open water (negative)
-      REAL(wp), POINTER, DIMENSION(:) ::   zdhex     ! excessively thick accreted sea ice (hlead-hice)
-
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqm0      ! old layer-system heat content
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zthick0   ! old ice thickness
-
-      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final   ! ice volume summed over categories
-      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   !  snow volume summed over categories
-      REAL(wp), POINTER, DIMENSION(:,:) ::   et_i_init, et_i_final   !  ice energy summed over categories
-      REAL(wp), POINTER, DIMENSION(:,:) ::   et_s_init               !  snow energy summed over categories
+      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_1d    ! total ice fraction    
+      REAL(wp), POINTER, DIMENSION(:) ::   zv_frazb    ! accretion of frazil ice at the ice bottom
+      REAL(wp), POINTER, DIMENSION(:) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector)
+
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_b      ! old volume of ice in category jl
+      REAL(wp), POINTER, DIMENSION(:,:) ::   za_b      ! old area of ice in category jl
+      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d   ! 1-D version of a_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d   ! 1-D version of v_i
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d ! 1-D version of smv_i
+
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_1d   !: 1-D version of e_i
+
       REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel                   ! relative ice / frazil velocity
+
+      REAL(wp) :: zcai = 1.4e-3_wp
       !!-----------------------------------------------------------------------!
 
-      CALL wrk_alloc( jpij, zcatac )   ! integer
+      CALL wrk_alloc( jpij, jcat )   ! integer
       CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
-      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
-      CALL wrk_alloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )
-      CALL wrk_alloc( jpij,jkmax,jpl, ze_i_ac )
-      CALL wrk_alloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
-      CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel )
-
-      et_i_init(:,:) = 0._wp
-      et_s_init(:,:) = 0._wp
-      vt_i_init(:,:) = 0._wp
-      vt_s_init(:,:) = 0._wp
-
-      !------------------------------------------------------------------------------!
-      ! 1) Conservation check and changes in each ice category
-      !------------------------------------------------------------------------------!
-      IF( con_i ) THEN
-         CALL lim_column_sum        ( jpl, v_i          , vt_i_init)
-         CALL lim_column_sum        ( jpl, v_s          , vt_s_init)
-         CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init)
-         CALL lim_column_sum        ( jpl, e_s(:,:,1,:) , et_s_init)
-      ENDIF
-
+      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d )
+      CALL wrk_alloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d )
+      CALL wrk_alloc( jpij,nlay_i,jpl, ze_i_1d )
+      CALL wrk_alloc( jpi,jpj, zvrel )
+
+      CALL lim_var_agg(1)
+      CALL lim_var_glo2eqv
       !------------------------------------------------------------------------------|
       ! 2) Convert units for ice internal energy
@@ -154,7 +132,6 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * REAL( nlay_i )
-                  zindb = 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) = e_i(ji,jj,jk,jl) * unit_fac * zindb
+                  rswitch          = 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
@@ -179,7 +156,7 @@
 
       ! Default new ice thickness 
-      hicol(:,:) = hiccrit(1)
-
-      IF( fraz_swi == 1._wp ) THEN
+      hicol(:,:) = rn_hnewice
+
+      IF( ln_frazil ) THEN
 
          !--------------------
@@ -190,27 +167,27 @@
          zhicrit = 0.04 ! frazil ice thickness
          ztwogp  = 2. * rau0 / ( grav * 0.3 * ( rau0 - rhoic ) ) ! reduced grav
-         zsqcd   = 1.0 / SQRT( 1.3 * cai ) ! 1/SQRT(airdensity*drag)
+         zsqcd   = 1.0 / SQRT( 1.3 * zcai ) ! 1/SQRT(airdensity*drag)
          zgamafr = 0.03
 
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
-               IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) ) > 0.e0 ) THEN
+         DO jj = 2, jpj
+            DO ji = 2, jpi
+               IF ( qlead(ji,jj) < 0._wp ) THEN
                   !-------------
                   ! Wind stress
                   !-------------
                   ! C-grid wind stress components
-                  ztaux         = ( utau_ice(ji-1,jj  ) * tmu(ji-1,jj  )   &
-                     &          +   utau_ice(ji  ,jj  ) * tmu(ji  ,jj  ) ) * 0.5_wp
-                  ztauy         = ( vtau_ice(ji  ,jj-1) * tmv(ji  ,jj-1)   &
-                     &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) * 0.5_wp
+                  ztaux         = ( utau_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)   &
+                     &          +   utau_ice(ji  ,jj  ) * umask(ji  ,jj  ,1) ) * 0.5_wp
+                  ztauy         = ( vtau_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)   &
+                     &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp
                   ! Square root of wind stress
-                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
+                  ztenagm       =  SQRT( SQRT( ztaux**2 + ztauy**2 ) )
 
                   !---------------------
                   ! Frazil ice velocity
                   !---------------------
-                  zvfrx         = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10)
-                  zvfry         = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10)
+                  rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )
+                  zvfrx   = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )
+                  zvfry   = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )
 
                   !-------------------
@@ -218,9 +195,7 @@
                   !-------------------
                   ! C-grid ice velocity
-                  zindb = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  )
-                  zvgx  = zindb * (  u_ice(ji-1,jj  ) * tmu(ji-1,jj  )    &
-                     &             + u_ice(ji,jj    ) * tmu(ji  ,jj  )  ) * 0.5_wp
-                  zvgy  = zindb * (  v_ice(ji  ,jj-1) * tmv(ji  ,jj-1)    &
-                     &             + v_ice(ji,jj    ) * tmv(ji  ,jj  )  ) * 0.5_wp
+                  rswitch = MAX(  0._wp, SIGN( 1._wp , at_i(ji,jj) )  )
+                  zvgx    = rswitch * ( u_ice(ji-1,jj  ) * umask(ji-1,jj  ,1)  + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp
+                  zvgy    = rswitch * ( v_ice(ji  ,jj-1) * vmask(ji  ,jj-1,1)  + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp
 
                   !-----------------------------------
@@ -248,5 +223,5 @@
                   iterate_frazil = .true.
 
-                  DO WHILE ( iter .LT. 100 .AND. iterate_frazil ) 
+                  DO WHILE ( iter < 100 .AND. iterate_frazil ) 
                      zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) &
                         - hicol(ji,jj) * zhicrit * ztwogp * zvrel2
@@ -264,4 +239,7 @@
             END DO ! loop on ji ends
          END DO ! loop on jj ends
+      ! 
+      CALL lbc_lnk( zvrel(:,:), 'T', 1. )
+      CALL lbc_lnk( hicol(:,:), 'T', 1. )
 
       ENDIF ! End of computation of frazil ice collection thickness
@@ -276,7 +254,9 @@
       ! This occurs if open water energy budget is negative
       nbpac = 0
+      npac(:) = 0
+      !
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) )  >  0._wp ) THEN
+            IF ( qlead(ji,jj)  <  0._wp ) THEN
                nbpac = nbpac + 1
                npac( nbpac ) = (jj - 1) * jpi + ji
@@ -287,8 +267,8 @@
       ! debug point to follow
       jiindex_1d = 0
-      IF( ln_nicep ) THEN
-         DO ji = mi0(jiindx), mi1(jiindx)
-            DO jj = mj0(jjindx), mj1(jjindx)
-               IF ( tms(ji,jj) * ( qcmif(ji,jj) - qldif(ji,jj) )  >  0._wp ) THEN
+      IF( ln_icectl ) THEN
+         DO ji = mi0(iiceprt), mi1(iiceprt)
+            DO jj = mj0(jiceprt), mj1(jiceprt)
+               IF ( qlead(ji,jj)  <  0._wp ) THEN
                   jiindex_1d = (jj - 1) * jpi + ji
                ENDIF
@@ -297,5 +277,5 @@
       ENDIF
    
-      IF( ln_nicep ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
+      IF( ln_icectl ) WRITE(numout,*) 'lim_thd_lac : nbpac = ', nbpac
 
       !------------------------------
@@ -307,22 +287,23 @@
       IF ( nbpac > 0 ) THEN
 
-         CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, zat_i_1d  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) )
          DO jl = 1, jpl
-            CALL tab_2d_1d( nbpac, za_i_ac  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
-            CALL tab_2d_1d( nbpac, zv_i_ac  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
-            CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) )
-            CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
+            CALL tab_2d_1d( nbpac, za_i_1d  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
+            CALL tab_2d_1d( nbpac, zv_i_1d  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) )
+            CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
             DO jk = 1, nlay_i
-               CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
-            END DO ! jk
-         END DO ! jl
-
-         CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif  , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif  , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, sfx_thd_1d(1:nbpac)     , sfx_thd, jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, rdm_ice_1d(1:nbpac)     , rdm_ice, jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) )
+               CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
+            END DO
+         END DO
+
+         CALL tab_2d_1d( nbpac, qlead_1d  (1:nbpac)     , qlead  , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, t_bo_1d   (1:nbpac)     , t_bo   , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, sfx_opw_1d(1:nbpac)     , sfx_opw, jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw, jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, hicol_1d  (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, zvrel_1d  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) )
+
+         CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac)     , hfx_thd, jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac)     , hfx_opw, jpi, jpj, npac(1:nbpac) )
 
          !------------------------------------------------------------------------------!
@@ -330,29 +311,32 @@
          !------------------------------------------------------------------------------!
 
+         !-----------------------------------------
+         ! Keep old ice areas and volume in memory
+         !-----------------------------------------
+         zv_b(1:nbpac,:) = zv_i_1d(1:nbpac,:) 
+         za_b(1:nbpac,:) = za_i_1d(1:nbpac,:)
          !----------------------
          ! Thickness of new ice
          !----------------------
          DO ji = 1, nbpac
-            zh_newice(ji) = hiccrit(1)
-         END DO
-         IF( fraz_swi == 1.0 )   zh_newice(:) = hicol_b(:)
+            zh_newice(ji) = rn_hnewice
+         END DO
+         IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac)
 
          !----------------------
          ! Salinity of new ice 
          !----------------------
-
-         SELECT CASE ( num_sal )
+         SELECT CASE ( nn_icesal )
          CASE ( 1 )                    ! Sice = constant 
-            zs_newice(:) = bulk_sal
+            zs_newice(1:nbpac) = rn_icesal
          CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
             DO ji = 1, nbpac
                ii =   MOD( npac(ji) - 1 , jpi ) + 1
                ij =      ( npac(ji) - 1 ) / jpi + 1
-               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij)  )
+               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , rn_simax , 0.5 * sss_m(ii,ij)  )
             END DO
          CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
-            zs_newice(:) =   2.3
+            zs_newice(1:nbpac) =   2.3
          END SELECT
-
 
          !-------------------------
@@ -361,11 +345,10 @@
          ! We assume that new ice is formed at the seawater freezing point
          DO ji = 1, nbpac
-            ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K)
-            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             &
-               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   &
-               &                       - rcp  *         ( ztmelts - rtt )  )
-            ze_newice(ji) =   MAX( ze_newice(ji) , 0._wp )    &
-               &          +   MAX(  0.0 , SIGN( 1.0 , - ze_newice(ji) )  ) * rhoic * lfus
-         END DO ! ji
+            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 )  )
+         END DO
+
          !----------------
          ! Age of new ice
@@ -373,12 +356,5 @@
          DO ji = 1, nbpac
             zo_newice(ji) = 0._wp
-         END DO ! ji
-
-         !--------------------------
-         ! Open water energy budget 
-         !--------------------------
-         DO ji = 1, nbpac
-            zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji)     !<0
-         END DO ! ji
+         END DO
 
          !-------------------
@@ -386,28 +362,33 @@
          !-------------------
          DO ji = 1, nbpac
-            zv_newice(ji) = - zqbgow(ji) / ze_newice(ji)
+
+            zEi           = - ze_newice(ji) * r1_rhoic             ! specific enthalpy of forming ice [J/kg]
+
+            zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg]
+                                                                   ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 
+                                                                   
+            zdE           = zEi - zEw                              ! specific enthalpy difference [J/kg]
+                                              
+            zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0) 
+                                                                   ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point   
+            zv_newice(ji) = - zfmdt * r1_rhoic
+
+            zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux  
+
+            ! Contribution to heat flux to the ocean [W.m-2], >0  
+            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice
+            ! Total heat flux used in this process [W.m-2]  
+            hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice
+            ! mass flux
+            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice
+            ! salt flux
+            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice
 
             ! A fraction zfrazb of frazil ice is accreted at the ice bottom
-            zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
-            zdh_frazb(ji) =         zfrazb   * zv_newice(ji)
+            rswitch       = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) )
+            zfrazb        = rswitch * ( TANH ( rn_Cfrazb * ( zvrel_1d(ji) - rn_vfrazb ) ) + 1.0 ) * 0.5 * rn_maxfrazb
+            zv_frazb(ji)  =         zfrazb   * zv_newice(ji)
             zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
          END DO
-
-         !------------------------------------
-         ! Diags for energy conservation test
-         !------------------------------------
-         DO ji = 1, nbpac
-            ii = MOD( npac(ji) - 1 , jpi ) + 1
-            ij =    ( npac(ji) - 1 ) / jpi + 1
-            !
-            zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji)
-            !
-            vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji)             ! volume
-            et_i_init(ii,ij) = et_i_init(ii,ij) + zde                       ! Energy
-
-         END DO
-
-         ! keep new ice volume in memory
-         CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj )
 
          !-----------------
@@ -415,9 +396,6 @@
          !-----------------
          DO ji = 1, nbpac
-            ii = MOD( npac(ji) - 1 , jpi ) + 1
-            ij =    ( npac(ji) - 1 ) / jpi + 1
             za_newice(ji) = zv_newice(ji) / zh_newice(ji)
-            diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem
-         END DO !ji
+         END DO
 
          !------------------------------------------------------------------------------!
@@ -425,18 +403,12 @@
          !------------------------------------------------------------------------------!
 
-         !-----------------------------------------
-         ! Keep old ice areas and volume in memory
-         !-----------------------------------------
-         zv_old(:,:) = zv_i_ac(:,:) 
-         za_old(:,:) = za_i_ac(:,:)
-
-         !-------------------------------------------
-         ! Compute excessive new ice area and volume
-         !-------------------------------------------
+         !------------------------
+         ! 6.1) lateral ice growth
+         !------------------------
          ! If lateral ice growth gives an ice concentration gt 1, then
          ! we keep the excessive volume in memory and attribute it later to bottom accretion
          DO ji = 1, nbpac
-            IF ( za_newice(ji) >  ( amax - zat_i_ac(ji) ) ) THEN
-               zda_res(ji)   = za_newice(ji) - ( amax - zat_i_ac(ji) )
+            IF ( za_newice(ji) >  ( rn_amax - zat_i_1d(ji) ) ) THEN
+               zda_res(ji)   = za_newice(ji) - ( rn_amax - zat_i_1d(ji) )
                zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
                za_newice(ji) = za_newice(ji) - zda_res  (ji)
@@ -446,232 +418,114 @@
                zdv_res(ji) = 0._wp
             ENDIF
-         END DO ! ji
-
-         !------------------------------------------------
-         ! Laterally redistribute new ice volume and area
-         !------------------------------------------------
-         zat_i_ac(:) = 0._wp
+         END DO
+
+         ! find which category to fill
+         zat_i_1d(:) = 0._wp
          DO jl = 1, jpl
             DO ji = 1, nbpac
-               IF(  hi_max   (jl-1)  <   zh_newice(ji)   .AND.   &
-                  & zh_newice(ji)    <=  hi_max   (jl)         ) THEN
-                  za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji)
-                  zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji)
-                  zat_i_ac(ji)    = zat_i_ac(ji)    + za_i_ac  (ji,jl)
-                  zcatac  (ji)    = jl
+               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
                ENDIF
-            END DO
-         END DO
-
-         !----------------------------------
-         ! Heat content - lateral accretion
-         !----------------------------------
-         DO ji = 1, nbpac
-            jl = zcatac(ji)                                                           ! categroy in which new ice is put
-            zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) )             ! zindb=1 if ice =0 otherwise
-            zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb  ! old ice thickness
-            zdhex    (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) )           ! difference in thickness
-            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) )   ! ice totally new in jl category
+               zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d  (ji,jl)
+            END DO
+         END DO
+
+         ! Heat content
+         DO ji = 1, nbpac
+            jl = jcat(ji)                                                    ! categroy in which new ice is put
+            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_b(ji,jl) ) )   ! 0 if old ice
          END DO
 
          DO jk = 1, nlay_i
             DO ji = 1, nbpac
-               jl = zcatac(ji)
-               zqold   = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]
-               zalphai = MIN( zhice_old(ji,jl) * REAL( jk )     / REAL( nlay_i ), zh_newice(ji) )   &
-                  &    - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) )
-               ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji)                                     &
-                  + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl)  * zqold * zhice_old(ji,jl) / REAL( nlay_i )  &
-                  + za_newice(ji)  * ze_newice(ji) * zalphai                                       &
-                  + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) )
-            END DO
-         END DO
-
-         !-----------------------------------------------
-         ! Add excessive volume of new ice at the bottom
-         !-----------------------------------------------
-         ! If the ice concentration exceeds 1, the remaining volume of new ice
-         ! is equally redistributed among all ice categories in which there is
-         ! ice
-
-         ! Fraction of level ice
-         jm = 1
-         zat_i_lev(:) = 0._wp
-
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-            DO ji = 1, nbpac
-               zat_i_lev(ji) = zat_i_lev(ji) + za_i_ac(ji,jl) 
-            END DO
-         END DO
-
-         IF( ln_nicep .AND. jiindex_1d > 0 ) WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-            DO ji = 1, nbpac
-               zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) )
-               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi10 ) )  ! clem
-               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi10 )
-            END DO
-         END DO
-         IF( ln_nicep .AND. jiindex_1d > 0 )   WRITE(numout,*) ' zv_i_ac : ', zv_i_ac(jiindex_1d, 1:jpl)
-
-         !---------------------------------
-         ! Heat content - bottom accretion
-         !---------------------------------
-         jm = 1
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-            DO ji = 1, nbpac
-               zindb =  1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) )       ! zindb=1 if ice =0 otherwise
-               zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
-               zdhicbot (ji,jl) = zdv_res(ji)    / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb    &
-                  &             +  zindb * zdh_frazb(ji)                               ! frazil ice may coalesce
-               zdummy(ji,jl)    = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb      ! thickness of residual ice
-            END DO
-         END DO
-
-         ! old layers thicknesses and enthalpies
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
+               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 )
+            END DO
+         END DO
+
+         !------------------------------------------------
+         ! 6.2) bottom ice growth + ice enthalpy remapping
+         !------------------------------------------------
+         DO jl = 1, jpl
+
+            ! for remapping
+            h_i_old (1:nbpac,0:nlay_i+1) = 0._wp
+            qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp
             DO jk = 1, nlay_i
                DO ji = 1, nbpac
-                  zthick0(ji,jk,jl) =  zhice_old(ji,jl) / REAL( nlay_i )
-                  zqm0   (ji,jk,jl) =  ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)
+                  h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i
+                  qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk)
                END DO
             END DO
-         END DO
-!!gm ???  why the previous do loop  if ocerwriten by the following one ?
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
+
+            ! new volumes including lateral/bottom accretion + residual
             DO ji = 1, nbpac
-               zthick0(ji,nlay_i+1,jl) =  zdhicbot(ji,jl)
-               zqm0   (ji,nlay_i+1,jl) =  ze_newice(ji) * zdhicbot(ji,jl)
-            END DO ! ji
-         END DO ! jl
-
-         ! Redistributing energy on the new grid
-         ze_i_ac(:,:,:) = 0._wp
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-            DO jk = 1, nlay_i
-               DO layer = 1, nlay_i + 1
-                  DO ji = 1, nbpac
-                     zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) 
-                     ! Redistributing energy on the new grid
-                     zweight = MAX (  MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) )   &
-                        &    - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp )   &
-                        &    /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb
-                     ze_i_ac(ji,jk,jl) =  ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl)  
-                  END DO ! ji
-               END DO ! layer
-            END DO ! jk
-         END DO ! jl
-
-         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-            DO jk = 1, nlay_i
-               DO ji = 1, nbpac
-                  zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 
-                  ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl)   &
-                     &              / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb
-               END DO
-            END DO
-         END DO
-
-         !------------
-         ! Update age 
-         !------------
-         DO jl = 1, jpl
-            DO ji = 1, nbpac
-               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
-               zoa_i_ac(ji,jl)  = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb   
-            END DO 
-         END DO   
+               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
+               ! for remapping
+               h_i_old (ji,nlay_i+1) = zv_newfra
+               qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra
+            ENDDO
+            ! --- Ice enthalpy remapping --- !
+            CALL lim_thd_ent( 1, nbpac, ze_i_1d(1:nbpac,:,jl) ) 
+         ENDDO
 
          !-----------------
          ! Update salinity
          !-----------------
-         !clem IF(  num_sal == 2  ) THEN
-            DO jl = 1, jpl
-               DO ji = 1, nbpac
-                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
-                  zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
-                  zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif
-               END DO
-            END DO   
-         !clem ENDIF
-
-         !--------------------------------
-         ! Update mass/salt fluxes (clem)
-         !--------------------------------
          DO jl = 1, jpl
             DO ji = 1, nbpac
-               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
-               zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
-               rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb
-               sfx_thd_1d(ji)   =   sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb
-           END DO
+               zdv   = zv_i_1d(ji,jl) - zv_b(ji,jl)
+               zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji)
+            END DO
          END DO
 
          !------------------------------------------------------------------------------!
-         ! 8) Change 2D vectors to 1D vectors 
+         ! 7) Change 2D vectors to 1D vectors 
          !------------------------------------------------------------------------------!
          DO jl = 1, jpl
-            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj )
-            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj )
-            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj )
-            !clem IF (  num_sal == 2  )   &
-               CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
+            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj )
+            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj )
+            CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj )
             DO jk = 1, nlay_i
-               CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj )
-            END DO
-         END DO
-         CALL tab_1d_2d( nbpac, sfx_thd, npac(1:nbpac), sfx_thd_1d(1:nbpac), jpi, jpj )
-         CALL tab_1d_2d( nbpac, rdm_ice, npac(1:nbpac), rdm_ice_1d(1:nbpac), jpi, jpj )
+               CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_1d(1:nbpac,jk,jl), jpi, jpj )
+            END DO
+         END DO
+         CALL tab_1d_2d( nbpac, sfx_opw, npac(1:nbpac), sfx_opw_1d(1:nbpac), jpi, jpj )
+         CALL tab_1d_2d( nbpac, wfx_opw, npac(1:nbpac), wfx_opw_1d(1:nbpac), jpi, jpj )
+
+         CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj )
+         CALL tab_1d_2d( nbpac, hfx_opw, npac(1:nbpac), hfx_opw_1d(1:nbpac), jpi, jpj )
          !
       ENDIF ! nbpac > 0
 
       !------------------------------------------------------------------------------!
-      ! 9) Change units for e_i
+      ! 8) Change units for e_i
       !------------------------------------------------------------------------------!    
       DO jl = 1, jpl
-         DO jk = 1, nlay_i          ! heat content in 10^9 Joules
-            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i )  / unit_fac 
+         DO jk = 1, nlay_i
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  ! heat content in J/m2
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 
+               END DO
+            END DO
          END DO
       END DO
 
-      !------------------------------------------------------------------------------|
-      ! 10) Conservation check and changes in each ice category
-      !------------------------------------------------------------------------------|
-      IF( con_i ) THEN 
-         CALL lim_column_sum (jpl,   v_i, vt_i_final)
-         fieldid = 'v_i, limthd_lac'
-         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
-         !
-         CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final)
-         fieldid = 'e_i, limthd_lac'
-         CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid) 
-         !
-         CALL lim_column_sum (jpl,   v_s, vt_s_final)
-         fieldid = 'v_s, limthd_lac'
-         CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid) 
-         !
-         !     CALL lim_column_sum (jpl,   e_s(:,:,1,:) , et_s_init)
-         !     fieldid = 'e_s, limthd_lac'
-         !     CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid) 
-         IF( ln_nicep ) THEN
-            DO ji = mi0(jiindx), mi1(jiindx)
-               DO jj = mj0(jjindx), mj1(jjindx)
-                  WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj)
-                  WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj)
-                  WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj)
-                  WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj)
-               END DO
-            END DO
-         ENDIF
-         !
-      ENDIF
       !
-      CALL wrk_dealloc( jpij, zcatac )   ! integer
+      CALL wrk_dealloc( jpij, jcat )   ! integer
       CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice )
-      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zdh_frazb, zvrel_ac, zqbgow, zdhex )
-      CALL wrk_dealloc( jpij,jpl, zhice_old, zdummy, zdhicbot, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac )
-      CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac )
-      CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 )
-      CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel )
+      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zv_frazb, zvrel_1d )
+      CALL wrk_dealloc( jpij,jpl, zv_b, za_b, za_i_1d, zv_i_1d, zsmv_i_1d )
+      CALL wrk_dealloc( jpij,nlay_i,jpl, ze_i_1d )
+      CALL wrk_dealloc( jpi,jpj, zvrel )
       !
    END SUBROUTINE lim_thd_lac
