Index: NEMO/trunk/src/ICE/iceistate.F90
===================================================================
--- NEMO/trunk/src/ICE/iceistate.F90	(revision 10994)
+++ NEMO/trunk/src/ICE/iceistate.F90	(revision 11229)
@@ -504,6 +504,5 @@
       !!
       !!-----------------------------------------------------------------------------
-      INTEGER ::   ji, jj
-      INTEGER ::   ios, ierr, inum_ice   ! Local integer output status for namelist read
+      INTEGER ::   ios   ! Local integer output status for namelist read
       INTEGER ::   ifpr, ierror
       !
Index: NEMO/trunk/src/ICE/icethd_do.F90
===================================================================
--- NEMO/trunk/src/ICE/icethd_do.F90	(revision 10994)
+++ NEMO/trunk/src/ICE/icethd_do.F90	(revision 11229)
@@ -128,6 +128,6 @@
 
       ! Default new ice thickness
-      WHERE( qlead(:,:) < 0._wp )   ;   ht_i_new(:,:) = rn_hinew
-      ELSEWHERE                     ;   ht_i_new(:,:) = 0._wp
+      WHERE( qlead(:,:) < 0._wp  .AND. tau_icebfr(:,:) == 0._wp )   ;   ht_i_new(:,:) = rn_hinew ! if cooling and no landfast
+      ELSEWHERE                                                     ;   ht_i_new(:,:) = 0._wp
       END WHERE
 
@@ -182,4 +182,7 @@
                   END DO
                   !
+                  ! bound ht_i_new (though I don't see why it should be necessary)
+                  ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) )
+                  !
                ENDIF
                !
Index: NEMO/trunk/src/ICE/icevar.F90
===================================================================
--- NEMO/trunk/src/ICE/icevar.F90	(revision 10994)
+++ NEMO/trunk/src/ICE/icevar.F90	(revision 11229)
@@ -73,9 +73,12 @@
    PUBLIC   ice_var_zapneg
    PUBLIC   ice_var_roundoff
-   PUBLIC   ice_var_itd
-   PUBLIC   ice_var_itd2
    PUBLIC   ice_var_bv           
    PUBLIC   ice_var_enthalpy           
    PUBLIC   ice_var_sshdyn
+   PUBLIC   ice_var_itd
+
+   INTERFACE ice_var_itd
+      MODULE PROCEDURE ice_var_itd_1c1c, ice_var_itd_Nc1c, ice_var_itd_1cMc, ice_var_itd_NcMc
+   END INTERFACE
 
    !!----------------------------------------------------------------------
@@ -656,10 +659,154 @@
    END SUBROUTINE ice_var_roundoff
    
+
+   SUBROUTINE ice_var_bv
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE ice_var_bv ***
+      !!
+      !! ** Purpose :   computes mean brine volume (%) in sea ice
+      !!
+      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
+      !!
+      !! References : Vancoppenolle et al., JGR, 2007
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
+      !!-------------------------------------------------------------------
+      !
+!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
+!!   instead of setting everything to zero as just below
+      bv_i (:,:,:) = 0._wp
+      DO jl = 1, jpl
+         DO jk = 1, nlay_i
+            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
+               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
+            END WHERE
+         END DO
+      END DO
+      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
+      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
+      END WHERE
+      !
+   END SUBROUTINE ice_var_bv
+
+
+   SUBROUTINE ice_var_enthalpy
+      !!-------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_var_enthalpy *** 
+      !!                 
+      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
+      !!
+      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ji, jk   ! dummy loop indices
+      REAL(wp) ::   ztmelts  ! local scalar 
+      !!-------------------------------------------------------------------
+      !
+      DO jk = 1, nlay_i             ! Sea ice energy of melting
+         DO ji = 1, npti
+            ztmelts      = - rTmlt  * sz_i_1d(ji,jk)
+            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue
+                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
+            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
+               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
+               &                   - rcp   * ztmelts )
+         END DO
+      END DO
+      DO jk = 1, nlay_s             ! Snow energy of melting
+         DO ji = 1, npti
+            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
+         END DO
+      END DO
+      !
+   END SUBROUTINE ice_var_enthalpy
+
    
-   SUBROUTINE ice_var_itd( zhti, zhts, zati, zh_i, zh_s, za_i )
-      !!-------------------------------------------------------------------
-      !!                ***  ROUTINE ice_var_itd   ***
-      !!
-      !! ** Purpose :  converting 1-cat ice to multiple ice categories
+   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_var_sshdyn  ***
+      !!                     
+      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
+      !!
+      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0
+      !!
+      !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
+      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*, 
+      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008
+      !!----------------------------------------------------------------------
+      !
+      ! input
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh            !: ssh [m]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass    !: mass of snow and ice at current  ice time step [Kg/m2]
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b  !: mass of snow and ice at previous ice time step [Kg/m2]
+      !
+      ! output
+      REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn  ! equivalent ssh in lead [m]
+      !
+      ! temporary
+      REAL(wp) :: zintn, zintb                     ! time interpolation weights []
+      REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m]
+      !
+      ! compute ice load used to define the equivalent ssh in lead
+      IF( ln_ice_embd ) THEN
+         !                                            
+         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
+         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1}
+         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
+         !
+         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
+         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
+         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
+         !
+         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
+         !
+      ELSE
+         zsnwiceload(:,:) = 0.0_wp
+      ENDIF
+      ! compute equivalent ssh in lead
+      ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
+      !
+   END FUNCTION ice_var_sshdyn
+
+   
+   !!-------------------------------------------------------------------
+   !!                ***  INTERFACE ice_var_itd   ***
+   !!
+   !! ** Purpose :  converting N-cat ice to jpl ice categories
+   !!-------------------------------------------------------------------
+   SUBROUTINE ice_var_itd_1c1c( zhti, zhts, zati, zh_i, zh_s, za_i )
+      !!-------------------------------------------------------------------
+      !! ** Purpose :  converting 1-cat ice to 1 ice category
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
+      REAL(wp), DIMENSION(:), INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
+      !!-------------------------------------------------------------------
+      zh_i(:) = zhti(:)
+      zh_s(:) = zhts(:)
+      za_i(:) = zati(:)
+   END SUBROUTINE ice_var_itd_1c1c
+
+   SUBROUTINE ice_var_itd_Nc1c( zhti, zhts, zati, zh_i, zh_s, za_i )
+      !!-------------------------------------------------------------------
+      !! ** Purpose :  converting N-cat ice to 1 ice category
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(:,:), INTENT(in)    ::   zhti, zhts, zati    ! input ice/snow variables
+      REAL(wp), DIMENSION(:)  , INTENT(inout) ::   zh_i, zh_s, za_i    ! output ice/snow variables
+      !!-------------------------------------------------------------------
+      !
+      za_i(:) = SUM( zati(:,:), dim=2 )
+      !
+      WHERE( za_i(:) /= 0._wp )
+         zh_i(:) = SUM( zhti(:,:) * zati(:,:), dim=2 ) / za_i(:)
+         zh_s(:) = SUM( zhts(:,:) * zati(:,:), dim=2 ) / za_i(:)
+      ELSEWHERE
+         zh_i(:) = 0._wp
+         zh_s(:) = 0._wp
+      END WHERE
+      !
+   END SUBROUTINE ice_var_itd_Nc1c
+   
+   SUBROUTINE ice_var_itd_1cMc( zhti, zhts, zati, zh_i, zh_s, za_i )
+      !!-------------------------------------------------------------------
+      !!
+      !! ** Purpose :  converting 1-cat ice to jpl ice categories
       !!
       !!                  ice thickness distribution follows a gaussian law
@@ -801,10 +948,8 @@
       END DO
       !
-   END SUBROUTINE ice_var_itd
-
-
-   SUBROUTINE ice_var_itd2( zhti, zhts, zati, zh_i, zh_s, za_i )
-      !!-------------------------------------------------------------------
-      !!                ***  ROUTINE ice_var_itd2   ***
+   END SUBROUTINE ice_var_itd_1cMc
+
+   SUBROUTINE ice_var_itd_NcMc( zhti, zhts, zati, zh_i, zh_s, za_i )
+      !!-------------------------------------------------------------------
       !!
       !! ** Purpose :  converting N-cat ice to jpl ice categories
@@ -845,197 +990,108 @@
       idim = SIZE( zhti, 1 )
       icat = SIZE( zhti, 2 )
-      !
-      ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
-      ALLOCATE( jlmin(idim), jlmax(idim) )
-
-      ! --- initialize output fields to 0 --- !
-      zh_i(1:idim,1:jpl) = 0._wp
-      zh_s(1:idim,1:jpl) = 0._wp
-      za_i(1:idim,1:jpl) = 0._wp
-      !
-      ! --- fill the categories --- !
-      !     find where cat-input = cat-output and fill cat-output fields  
-      jlmax(:) = 0
-      jlmin(:) = 999
-      jlfil(:,:) = 0
-      DO jl1 = 1, jpl
-         DO jl2 = 1, icat
+      !                                 ! ---------------------- !
+      IF( icat == jpl ) THEN            ! input cat = output cat !
+         !                              ! ---------------------- !
+         zh_i(:,:) = zhti(:,:)
+         zh_s(:,:) = zhts(:,:)
+         za_i(:,:) = zati(:,:)
+         !                              ! ---------------------- !
+      ELSEIF( icat == 1 ) THEN          ! input cat = 1          !
+         !                              ! ---------------------- !
+         CALL  ice_var_itd_1cMc( zhti(:,1), zhts(:,1), zati(:,1), zh_i(:,:), zh_s(:,:), za_i(:,:) )
+         !                              ! ---------------------- !
+      ELSEIF( jpl == 1 ) THEN           ! output cat = 1         !
+         !                              ! ---------------------- !
+         CALL  ice_var_itd_Nc1c( zhti(:,:), zhts(:,:), zati(:,:), zh_i(:,1), zh_s(:,1), za_i(:,1) )         
+         !                              ! ----------------------- !
+      ELSE                              ! input cat /= output cat !
+         !                              ! ----------------------- !
+         
+         ALLOCATE( jlfil(idim,jpl), jlfil2(idim,jpl) )       ! allocate arrays
+         ALLOCATE( jlmin(idim), jlmax(idim) )
+
+         ! --- initialize output fields to 0 --- !
+         zh_i(1:idim,1:jpl) = 0._wp
+         zh_s(1:idim,1:jpl) = 0._wp
+         za_i(1:idim,1:jpl) = 0._wp
+         !
+         ! --- fill the categories --- !
+         !     find where cat-input = cat-output and fill cat-output fields  
+         jlmax(:) = 0
+         jlmin(:) = 999
+         jlfil(:,:) = 0
+         DO jl1 = 1, jpl
+            DO jl2 = 1, icat
+               DO ji = 1, idim
+                  IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
+                     ! fill the right category
+                     zh_i(ji,jl1) = zhti(ji,jl2)
+                     zh_s(ji,jl1) = zhts(ji,jl2)
+                     za_i(ji,jl1) = zati(ji,jl2)
+                     ! record categories that are filled
+                     jlmax(ji) = MAX( jlmax(ji), jl1 )
+                     jlmin(ji) = MIN( jlmin(ji), jl1 )
+                     jlfil(ji,jl1) = jl1
+                  ENDIF
+               END DO
+            END DO
+         END DO
+         !
+         ! --- fill the gaps between categories --- !  
+         !     transfer from categories filled at the previous step to the empty ones in between
+         DO ji = 1, idim
+            jl1 = jlmin(ji)
+            jl2 = jlmax(ji)
+            IF( jl1 > 1 ) THEN
+               ! fill the lower cat (jl1-1)
+               za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
+               zh_i(ji,jl1-1) = hi_mean(jl1-1)
+               ! remove from cat jl1
+               za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
+            ENDIF
+            IF( jl2 < jpl ) THEN
+               ! fill the upper cat (jl2+1)
+               za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
+               zh_i(ji,jl2+1) = hi_mean(jl2+1)
+               ! remove from cat jl2
+               za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
+            ENDIF
+         END DO
+         !
+         jlfil2(:,:) = jlfil(:,:) 
+         ! fill categories from low to high
+         DO jl = 2, jpl-1
             DO ji = 1, idim
-               IF( hi_max(jl1-1) <= zhti(ji,jl2) .AND. hi_max(jl1) > zhti(ji,jl2) ) THEN
-                  ! fill the right category
-                  zh_i(ji,jl1) = zhti(ji,jl2)
-                  zh_s(ji,jl1) = zhts(ji,jl2)
-                  za_i(ji,jl1) = zati(ji,jl2)
-                  ! record categories that are filled
-                  jlmax(ji) = MAX( jlmax(ji), jl1 )
-                  jlmin(ji) = MIN( jlmin(ji), jl1 )
-                  jlfil(ji,jl1) = jl1
+               IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
+                  ! fill high
+                  za_i(ji,jl) = ztrans * za_i(ji,jl-1)
+                  zh_i(ji,jl) = hi_mean(jl)
+                  jlfil(ji,jl) = jl
+                  ! remove low
+                  za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
                ENDIF
             END DO
          END DO
-      END DO
-      !
-      ! --- fill the gaps between categories --- !  
-      !     transfer from categories filled at the previous step to the empty ones in between
-      DO ji = 1, idim
-         jl1 = jlmin(ji)
-         jl2 = jlmax(ji)
-         IF( jl1 > 1 ) THEN
-            ! fill the lower cat (jl1-1)
-            za_i(ji,jl1-1) = ztrans * za_i(ji,jl1)
-            zh_i(ji,jl1-1) = hi_mean(jl1-1)
-            ! remove from cat jl1
-            za_i(ji,jl1  ) = ( 1._wp - ztrans ) * za_i(ji,jl1)
-         ENDIF
-         IF( jl2 < jpl ) THEN
-            ! fill the upper cat (jl2+1)
-            za_i(ji,jl2+1) = ztrans * za_i(ji,jl2)
-            zh_i(ji,jl2+1) = hi_mean(jl2+1)
-            ! remove from cat jl2
-            za_i(ji,jl2  ) = ( 1._wp - ztrans ) * za_i(ji,jl2)
-         ENDIF
-      END DO
-      !
-      jlfil2(:,:) = jlfil(:,:) 
-      ! fill categories from low to high
-      DO jl = 2, jpl-1
-         DO ji = 1, idim
-            IF( jlfil(ji,jl-1) /= 0 .AND. jlfil(ji,jl) == 0 ) THEN
-               ! fill high
-               za_i(ji,jl) = ztrans * za_i(ji,jl-1)
-               zh_i(ji,jl) = hi_mean(jl)
-               jlfil(ji,jl) = jl
-               ! remove low
-               za_i(ji,jl-1) = ( 1._wp - ztrans ) * za_i(ji,jl-1)
-            ENDIF
-         END DO
-      END DO
-      !
-      ! fill categories from high to low
-      DO jl = jpl-1, 2, -1
-         DO ji = 1, idim
-            IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
-               ! fill low
-               za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
-               zh_i(ji,jl) = hi_mean(jl) 
-               jlfil2(ji,jl) = jl
-               ! remove high
-               za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
-            ENDIF
-         END DO
-      END DO
-      !
-      DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
-      DEALLOCATE( jlmin, jlmax )
-      !
-   END SUBROUTINE ice_var_itd2
-
-
-   SUBROUTINE ice_var_bv
-      !!-------------------------------------------------------------------
-      !!                ***  ROUTINE ice_var_bv ***
-      !!
-      !! ** Purpose :   computes mean brine volume (%) in sea ice
-      !!
-      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
-      !!
-      !! References : Vancoppenolle et al., JGR, 2007
-      !!-------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      !!-------------------------------------------------------------------
-      !
-!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
-!!   instead of setting everything to zero as just below
-      bv_i (:,:,:) = 0._wp
-      DO jl = 1, jpl
-         DO jk = 1, nlay_i
-            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
-               bv_i(:,:,jl) = bv_i(:,:,jl) - rTmlt * sz_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
-            END WHERE
-         END DO
-      END DO
-      WHERE( vt_i(:,:) > epsi20 )   ;   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
-      ELSEWHERE                     ;   bvm_i(:,:) = 0._wp
-      END WHERE
-      !
-   END SUBROUTINE ice_var_bv
-
-
-   SUBROUTINE ice_var_enthalpy
-      !!-------------------------------------------------------------------
-      !!                   ***  ROUTINE ice_var_enthalpy *** 
-      !!                 
-      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature
-      !!
-      !! ** Method  :   Formula (Bitz and Lipscomb, 1999)
-      !!-------------------------------------------------------------------
-      INTEGER  ::   ji, jk   ! dummy loop indices
-      REAL(wp) ::   ztmelts  ! local scalar 
-      !!-------------------------------------------------------------------
-      !
-      DO jk = 1, nlay_i             ! Sea ice energy of melting
-         DO ji = 1, npti
-            ztmelts      = - rTmlt  * sz_i_1d(ji,jk)
-            t_i_1d(ji,jk) = MIN( t_i_1d(ji,jk), ztmelts + rt0 ) ! Force t_i_1d to be lower than melting point => likely conservation issue
-                                                                !   (sometimes zdf scheme produces abnormally high temperatures)   
-            e_i_1d(ji,jk) = rhoi * ( rcpi  * ( ztmelts - ( t_i_1d(ji,jk) - rt0 ) )           &
-               &                   + rLfus * ( 1._wp - ztmelts / ( t_i_1d(ji,jk) - rt0 ) )   &
-               &                   - rcp   * ztmelts )
-         END DO
-      END DO
-      DO jk = 1, nlay_s             ! Snow energy of melting
-         DO ji = 1, npti
-            e_s_1d(ji,jk) = rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus )
-         END DO
-      END DO
-      !
-   END SUBROUTINE ice_var_enthalpy
-
-   FUNCTION ice_var_sshdyn(pssh, psnwice_mass, psnwice_mass_b)
-      !!---------------------------------------------------------------------
-      !!                   ***  ROUTINE ice_var_sshdyn  ***
-      !!                     
-      !! ** Purpose :  compute the equivalent ssh in lead when sea ice is embedded
-      !!
-      !! ** Method  :  ssh_lead = ssh + (Mice + Msnow) / rau0
-      !!
-      !! ** Reference : Jean-Michel Campin, John Marshall, David Ferreira,
-      !!                Sea ice-ocean coupling using a rescaled vertical coordinate z*, 
-      !!                Ocean Modelling, Volume 24, Issues 1-2, 2008
-      !!----------------------------------------------------------------------
-      !
-      ! input
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: pssh            !: ssh [m]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass    !: mass of snow and ice at current  ice time step [Kg/m2]
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) :: psnwice_mass_b  !: mass of snow and ice at previous ice time step [Kg/m2]
-      !
-      ! output
-      REAL(wp), DIMENSION(jpi,jpj) :: ice_var_sshdyn  ! equivalent ssh in lead [m]
-      !
-      ! temporary
-      REAL(wp) :: zintn, zintb                     ! time interpolation weights []
-      REAL(wp), DIMENSION(jpi,jpj) :: zsnwiceload  ! snow and ice load [m]
-      !
-      ! compute ice load used to define the equivalent ssh in lead
-      IF( ln_ice_embd ) THEN
-         !                                            
-         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
-         !                                               = (1/nn_fsbc)^2 * {SUM[n]        , n=0,nn_fsbc-1}
-         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp
-         !
-         ! average interpolation coeff as used in dynspg = (1/nn_fsbc)   *    {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
-         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
-         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
-         !
-         zsnwiceload(:,:) = ( zintn * psnwice_mass(:,:) + zintb * psnwice_mass_b(:,:) ) * r1_rau0
-         !
-      ELSE
-         zsnwiceload(:,:) = 0.0_wp
+         !
+         ! fill categories from high to low
+         DO jl = jpl-1, 2, -1
+            DO ji = 1, idim
+               IF( jlfil2(ji,jl+1) /= 0 .AND. jlfil2(ji,jl) == 0 ) THEN
+                  ! fill low
+                  za_i(ji,jl) = za_i(ji,jl) + ztrans * za_i(ji,jl+1)
+                  zh_i(ji,jl) = hi_mean(jl) 
+                  jlfil2(ji,jl) = jl
+                  ! remove high
+                  za_i(ji,jl+1) = ( 1._wp - ztrans ) * za_i(ji,jl+1)
+               ENDIF
+            END DO
+         END DO
+         !
+         DEALLOCATE( jlfil, jlfil2 )      ! deallocate arrays
+         DEALLOCATE( jlmin, jlmax )
+         !
       ENDIF
-      ! compute equivalent ssh in lead
-      ice_var_sshdyn(:,:) = pssh(:,:) + zsnwiceload(:,:)
-      !
-   END FUNCTION ice_var_sshdyn
-
+      !
+   END SUBROUTINE ice_var_itd_NcMc
 
 #else
