Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 6455)
@@ -253,20 +253,20 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhld        !: heat flux from the lead used for bottom melting
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw    !: snow-ocean mass exchange over 1 time step [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr    !: snow precipitation on ice over 1 time step [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub    !: snow sublimation over 1 time step [kg/m2]
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice    !: ice-ocean mass exchange over 1 time step [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sni    !: snow ice growth component of wfx_ice [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_opw    !: lateral ice growth component of wfx_ice [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bog    !: bottom ice growth component of wfx_ice [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_dyn    !: dynamical ice growth component of wfx_ice [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bom    !: bottom melt component of wfx_ice [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sum    !: surface melt component of wfx_ice [kg/m2]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res    !: residual component of wfx_ice [kg/m2]
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total) [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw    !: snow-ocean mass exchange   [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr    !: snow precipitation on ice  [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub    !: snow/ice sublimation       [kg.m-2.s-1]
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice    !: ice-ocean mass exchange                   [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sni    !: snow ice growth component of wfx_ice      [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_opw    !: lateral ice growth component of wfx_ice   [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bog    !: bottom ice growth component of wfx_ice    [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_dyn    !: dynamical ice growth component of wfx_ice [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bom    !: bottom melt component of wfx_ice          [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sum    !: surface melt component of wfx_ice         [kg.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res    !: residual component of wfx_ice             [kg.m-2.s-1]
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total)          [s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_thd     !: ice concentration tendency (thermodynamics) [s-1]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_dyn     !: ice concentration tendency (dynamics) [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_dyn     !: ice concentration tendency (dynamics)       [s-1]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s]
@@ -279,26 +279,30 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: residual salt flux due to correction of ice thickness [PSU/m2/s]
 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bom     !: total heat flux causing bottom ice melt 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sum     !: total heat flux causing surface ice melt 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_opw     !: total heat flux causing open water ice formation
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err     !: heat flux error after heat diffusion 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_out     !: heat flux remaining at the end of thermo transformations 
-
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sub     !: salt flux due to ice sublimation
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth        [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bom     !: total heat flux causing bottom ice melt          [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sum     !: total heat flux causing surface ice melt         [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_opw     !: total heat flux causing open water ice formation [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice   [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt                          [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err     !: heat flux error after heat diffusion             [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping             [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations   [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_out     !: heat flux remaining at the end of thermo transformations  [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_err_sub !: mass flux error after sublimation [kg.m-2.s-1]
+   
    ! heat flux associated with ice-atmosphere mass exchange
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation  [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation  [W.m-2]
 
    ! heat flux associated with ice-ocean mass exchange
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_thd     !: ice-ocean heat flux from thermo processes (limthd_dh) 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dyn     !: ice-ocean heat flux from mecanical processes (limitd_me) 
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_thd     !: ice-ocean heat flux from thermo processes (limthd_dh)  [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dyn     !: ice-ocean heat flux from mecanical processes (limitd_me)  [W.m-2]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness [W.m-2]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_amax_2d  !: maximum ice concentration 2d array
 
    !!--------------------------------------------------------------------------
@@ -372,11 +376,12 @@
    INTEGER          , PUBLIC ::   nlay_i          !: number of ice  layers 
    INTEGER          , PUBLIC ::   nlay_s          !: number of snow layers 
-   CHARACTER(len=32), PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input)
+   CHARACTER(len=80), PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input)
    CHARACTER(len=256), PUBLIC ::   cn_icerst_indir !: ice restart input directory
-   CHARACTER(len=32), PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output)
+   CHARACTER(len=80), PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output)
    CHARACTER(len=256), PUBLIC ::   cn_icerst_outdir!: ice restart output directory
    LOGICAL          , PUBLIC ::   ln_limdyn       !: flag for ice dynamics (T) or not (F)
    LOGICAL          , PUBLIC ::   ln_icectl       !: flag for sea-ice points output (T) or not (F)
-   REAL(wp)         , PUBLIC ::   rn_amax         !: maximum ice concentration
+   REAL(wp)         , PUBLIC ::   rn_amax_n       !: maximum ice concentration Northern hemisphere
+   REAL(wp)         , PUBLIC ::   rn_amax_s       !: maximum ice concentration Southern hemisphere
    INTEGER          , PUBLIC ::   iiceprt         !: debug i-point
    INTEGER          , PUBLIC ::   jiceprt         !: debug j-point
@@ -438,8 +443,9 @@
          &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) ,                        &
          &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl), qlead  (jpi,jpj) ,                     &
-         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,                        &
+         &      rn_amax_2d(jpi,jpj),                                                            &
+         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) ,                       &
          &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,    &
          &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) ,     & 
-         &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) ,                                   &
+         &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) , wfx_err_sub(jpi,jpj) ,       &
          &      hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) ,                           &
          &      hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , hfx_opw(jpi,jpj) ,    &
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 6455)
@@ -24,5 +24,6 @@
    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE sbc_oce , ONLY : sfx  ! Surface boundary condition: ocean fields
-
+   USE sbc_ice , ONLY : qevap_ice
+   
    IMPLICIT NONE
    PRIVATE
@@ -184,5 +185,5 @@
          ! salt flux
          zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &
+            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:)                   &
             &                ) *  e12t(:,:) * tmask(:,:,1) * zconv )
 
@@ -209,5 +210,5 @@
          ! salt flux
          zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
+            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:)                   & 
             &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zfs_b
 
@@ -256,5 +257,6 @@
             ENDIF
             IF (     zvmin   < -epsi10 ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',zvmin
-            IF (     zamax   > rn_amax+epsi10 .AND. cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' ) THEN
+            IF (     zamax   > MAX( rn_amax_n, rn_amax_s ) + epsi10 .AND. &
+               &                         cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' ) THEN
                                          WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax
             ENDIF
@@ -286,5 +288,6 @@
 #if ! defined key_bdy
       ! heat flux
-      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - hfx_sub ) * e12t * tmask(:,:,1) * zconv ) 
+      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - SUM( qevap_ice * a_i_b, dim=3 ) )  &
+         &              * e12t * tmask(:,:,1) * zconv ) 
       ! salt flux
       zsfx  = glob_sum( ( sfx + diag_smvi ) * e12t * tmask(:,:,1) * zconv ) * rday
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 6455)
@@ -56,5 +56,5 @@
       real(wp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc
       real(wp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni,   &
-      &               zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn 
+      &               zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn, zbg_sfx_sub 
       real(wp)   ::   zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn
       real(wp)   ::   zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res, zbg_vfx_spr, zbg_vfx_snw, zbg_vfx_sub  
@@ -111,4 +111,5 @@
       zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) )
       zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) )
+      zbg_sfx_sub = ztmp * glob_sum( sfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) )
 
       ! Heat budget
@@ -189,4 +190,5 @@
       CALL iom_put( 'ibgsfxbom' , zbg_sfx_bom                              )   ! salt flux bottom melt       -
       CALL iom_put( 'ibgsfxsum' , zbg_sfx_sum                              )   ! salt flux surface melt      -
+      CALL iom_put( 'ibgsfxsub' , zbg_sfx_sub                              )   ! salt flux sublimation      -
 
       CALL iom_put( 'ibghfxdhc' , zbg_hfx_dhc                              )   ! Heat content variation in snow and ice [W]
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 6455)
@@ -45,26 +45,17 @@
    REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
    REAL(wp), ALLOCATABLE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/
-   !                                                     ! closing associated w/ category n
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/closing associated w/ category n
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   hraft    ! thickness of rafted ice
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! mean ridge thickness/thickness of ridging ice 
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   krdg     ! thickness of ridging ice / mean ridge thickness
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   aridge   ! participating ice ridging
    REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   araft    ! participating ice rafting
 
    REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier
-   REAL(wp), PARAMETER ::   kraft   = 2.0_wp    ! rafting multipliyer
-   REAL(wp), PARAMETER ::   kamax   = 1.0_wp    ! max of ice area authorized (clem: scheme is not stable if kamax <= 0.99)
+   REAL(wp), PARAMETER ::   kraft   = 0.5_wp    ! rafting multipliyer
 
    REAL(wp) ::   Cp                             ! 
    !
-   !-----------------------------------------------------------------------
-   ! Ridging diagnostic arrays for history files
-   !-----------------------------------------------------------------------
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dardg1dt   ! rate of fractional area loss by ridging ice (1/s)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dardg2dt   ! rate of fractional area gain by new ridges (1/s)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dvirdgdt   ! rate of ice volume ridged (m/s)
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   opening    ! rate of opening due to divergence/shear (1/s)
    !
    !!----------------------------------------------------------------------
@@ -83,11 +74,6 @@
          &      asum (jpi,jpj)     , athorn(jpi,jpj,0:jpl)                    ,     &
          &      aksum(jpi,jpj)                                                ,     &
-         !
          &      hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) ,     &
-         &      hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) ,     &
-         !
-         !* Ridging diagnostic arrays for history files
-         &      dardg1dt(jpi,jpj)  , dardg2dt(jpi,jpj)                        ,     & 
-         &      dvirdgdt(jpi,jpj)  , opening(jpi,jpj)                         , STAT=lim_itd_me_alloc )
+         &      hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , STAT=lim_itd_me_alloc )
          !
       IF( lim_itd_me_alloc /= 0 )   CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' )
@@ -132,7 +118,4 @@
       REAL(wp), POINTER, DIMENSION(:,:)   ::   opning          ! rate of opening due to divergence/shear
       REAL(wp), POINTER, DIMENSION(:,:)   ::   closing_gross   ! rate at which area removed, not counting area of new ridges
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   msnow_mlt       ! mass of snow added to ocean (kg m-2)
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2)
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   vt_i_init, vt_i_final  !  ice volume summed over categories
       !
       INTEGER, PARAMETER ::   nitermax = 20    
@@ -142,5 +125,5 @@
       IF( nn_timing == 1 )  CALL timing_start('limitd_me')
 
-      CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
+      CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross )
 
       IF(ln_ctl) THEN
@@ -154,7 +137,4 @@
       IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-      CALL lim_var_zapsmall
-      CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting
-
       !-----------------------------------------------------------------------------!
       ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
@@ -164,14 +144,7 @@
       CALL lim_itd_me_ridgeprep                                    ! prepare ridging
       !
-      IF( con_i)   CALL lim_column_sum( jpl, v_i, vt_i_init )      ! conservation check
 
       DO jj = 1, jpj                                     ! Initialize arrays.
          DO ji = 1, jpi
-            msnow_mlt(ji,jj) = 0._wp
-            esnow_mlt(ji,jj) = 0._wp
-            dardg1dt (ji,jj) = 0._wp
-            dardg2dt (ji,jj) = 0._wp
-            dvirdgdt (ji,jj) = 0._wp
-            opening  (ji,jj) = 0._wp
 
             !-----------------------------------------------------------------------------!
@@ -204,6 +177,6 @@
             ! If divu_adv < 0, make sure the closing rate is large enough
             ! to give asum = 1.0 after ridging.
-
-            divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep
+            
+            divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep
 
             IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
@@ -224,25 +197,25 @@
       DO WHILE ( iterate_ridging > 0 .AND. niter < nitermax )
 
+         ! 3.2 closing_gross
+         !-----------------------------------------------------------------------------!
+         ! Based on the ITD of ridging and ridged ice, convert the net
+         !  closing rate to a gross closing rate.  
+         ! NOTE: 0 < aksum <= 1
+         closing_gross(:,:) = closing_net(:,:) / aksum(:,:)
+
+         ! correction to closing rate and opening if closing rate is excessive
+         !---------------------------------------------------------------------
+         ! Reduce the closing rate if more than 100% of the open water 
+         ! would be removed.  Reduce the opening rate proportionately.
          DO jj = 1, jpj
             DO ji = 1, jpi
-
-               ! 3.2 closing_gross
-               !-----------------------------------------------------------------------------!
-               ! Based on the ITD of ridging and ridged ice, convert the net
-               !  closing rate to a gross closing rate.  
-               ! NOTE: 0 < aksum <= 1
-               closing_gross(ji,jj) = closing_net(ji,jj) / aksum(ji,jj)
-
-               ! correction to closing rate and opening if closing rate is excessive
-               !---------------------------------------------------------------------
-               ! Reduce the closing rate if more than 100% of the open water 
-               ! would be removed.  Reduce the opening rate proportionately.
-               za   = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
-               IF( za > epsi20 ) THEN
-                  zfac = MIN( 1._wp, ato_i(ji,jj) / za )
-                  closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
-                  opning       (ji,jj) = opning       (ji,jj) * zfac
+               za   = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice
+               IF( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN  ! would lead to negative ato_i
+                  zfac = - ato_i(ji,jj) / za
+                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 
+               ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN  ! would lead to ato_i > asum
+                  zfac = ( asum(ji,jj) - ato_i(ji,jj) ) / za
+                  opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 
                ENDIF
-
             END DO
          END DO
@@ -256,8 +229,7 @@
                DO ji = 1, jpi
                   za = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
-                  IF( za  >  epsi20 ) THEN
-                     zfac = MIN( 1._wp, a_i(ji,jj,jl) / za )
+                  IF( za  >  a_i(ji,jj,jl) ) THEN
+                     zfac = a_i(ji,jj,jl) / za
                      closing_gross(ji,jj) = closing_gross(ji,jj) * zfac
-                     opning       (ji,jj) = opning       (ji,jj) * zfac
                   ENDIF
                END DO
@@ -268,13 +240,11 @@
          !-----------------------------------------------------------------------------!
 
-         CALL lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
-
+         CALL lim_itd_me_ridgeshift( opning, closing_gross )
+
+         
          ! 3.4 Compute total area of ice plus open water after ridging.
          !-----------------------------------------------------------------------------!
          ! This is in general not equal to one because of divergence during transport
-         asum(:,:) = ato_i(:,:)
-         DO jl = 1, jpl
-            asum(:,:) = asum(:,:) + a_i(:,:,jl)
-         END DO
+         asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
 
          ! 3.5 Do we keep on iterating ???
@@ -284,13 +254,12 @@
 
          iterate_ridging = 0
-
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN
+               IF( ABS( asum(ji,jj) - 1._wp ) < epsi10 ) THEN
                   closing_net(ji,jj) = 0._wp
                   opning     (ji,jj) = 0._wp
                ELSE
                   iterate_ridging    = 1
-                  divu_adv   (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice
+                  divu_adv   (ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice
                   closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
                   opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) )
@@ -309,60 +278,12 @@
 
          IF( iterate_ridging == 1 ) THEN
+            CALL lim_itd_me_ridgeprep
             IF( niter  >  nitermax ) THEN
                WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
                WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
             ENDIF
-            CALL lim_itd_me_ridgeprep
          ENDIF
 
       END DO !! on the do while over iter
-
-      !-----------------------------------------------------------------------------!
-      ! 4) Ridging diagnostics
-      !-----------------------------------------------------------------------------!
-      ! Convert ridging rate diagnostics to correct units.
-      ! Update fresh water and heat fluxes due to snow melt.
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-
-            dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice
-            dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice
-            dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice
-            opening (ji,jj) = opening (ji,jj) * r1_rdtice
-
-            !-----------------------------------------------------------------------------!
-            ! 5) Heat, salt and freshwater fluxes
-            !-----------------------------------------------------------------------------!
-            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) * r1_rdtice     ! heat sink for ocean (<0, W.m-2)
-
-         END DO
-      END DO
-
-      ! Check if there is a ridging error
-      IF( lwp ) THEN
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( ABS( asum(ji,jj) - kamax)  >  epsi10 ) THEN   ! there is a bug
-                  WRITE(numout,*) ' '
-                  WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
-                  WRITE(numout,*) ' limitd_me '
-                  WRITE(numout,*) ' POINT : ', ji, jj
-                  WRITE(numout,*) ' jpl, a_i, athorn '
-                  WRITE(numout,*) 0, ato_i(ji,jj), athorn(ji,jj,0)
-                  DO jl = 1, jpl
-                     WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
-                  END DO
-               ENDIF
-            END DO
-         END DO
-      END IF
-
-      ! Conservation check
-      IF ( con_i ) THEN
-         CALL lim_column_sum (jpl,   v_i, vt_i_final)
-         fieldid = ' v_i : limitd_me '
-         CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid) 
-      ENDIF
 
       CALL lim_var_agg( 1 ) 
@@ -410,9 +331,441 @@
       ENDIF  ! ln_limdyn=.true.
       !
-      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
+      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross )
       !
       IF( nn_timing == 1 )  CALL timing_stop('limitd_me')
    END SUBROUTINE lim_itd_me
 
+   SUBROUTINE lim_itd_me_ridgeprep
+      !!---------------------------------------------------------------------!
+      !!                ***  ROUTINE lim_itd_me_ridgeprep ***
+      !!
+      !! ** Purpose :   preparation for ridging and strength calculations
+      !!
+      !! ** Method  :   Compute the thickness distribution of the ice and open water 
+      !!              participating in ridging and of the resulting ridges.
+      !!---------------------------------------------------------------------!
+      INTEGER ::   ji,jj, jl    ! dummy loop indices
+      REAL(wp) ::   Gstari, astari, hrmean, zdummy   ! local scalar
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n
+      !------------------------------------------------------------------------------!
+
+      CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
+
+      Gstari     = 1.0/rn_gstar    
+      astari     = 1.0/rn_astar    
+      aksum(:,:)    = 0.0
+      athorn(:,:,:) = 0.0
+      aridge(:,:,:) = 0.0
+      araft (:,:,:) = 0.0
+
+      ! Zero out categories with very small areas
+      CALL lim_var_zapsmall
+
+      ! Ice thickness needed for rafting
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )
+               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
+            END DO
+         END DO
+      END DO
+
+      !------------------------------------------------------------------------------!
+      ! 1) Participation function 
+      !------------------------------------------------------------------------------!
+
+      ! Compute total area of ice plus open water.
+      ! This is in general not equal to one because of divergence during transport
+      asum(:,:) = ato_i(:,:) + SUM( a_i, dim=3 )
+
+      ! Compute cumulative thickness distribution function
+      ! Compute the cumulative thickness distribution function Gsum,
+      ! where Gsum(n) is the fractional area in categories 0 to n.
+      ! initial value (in h = 0) equals open water area
+      Gsum(:,:,-1) = 0._wp
+      Gsum(:,:,0 ) = ato_i(:,:)
+      ! for each value of h, you have to add ice concentration then
+      DO jl = 1, jpl
+         Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl)
+      END DO
+
+      ! Normalize the cumulative distribution to 1
+      DO jl = 0, jpl
+         Gsum(:,:,jl) = Gsum(:,:,jl) / asum(:,:)
+      END DO
+
+      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
+      !--------------------------------------------------------------------------------------------------
+      ! Compute the participation function athorn; this is analogous to
+      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
+      ! area lost from category n due to ridging/closing
+      ! athorn(n)   = total area lost due to ridging/closing
+      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar). 
+      !
+      ! The expressions for athorn are found by integrating b(h)g(h) between
+      ! the category boundaries.
+      ! athorn is always >= 0 and SUM(athorn(0:jpl))=1
+      !-----------------------------------------------------------------
+
+      IF( nn_partfun == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
+         DO jl = 0, jpl    
+            DO jj = 1, jpj 
+               DO ji = 1, jpi
+                  IF    ( Gsum(ji,jj,jl)   < rn_gstar ) THEN
+                     athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &
+                        &                        ( 2._wp - ( Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )
+                  ELSEIF( Gsum(ji,jj,jl-1) < rn_gstar ) THEN
+                     athorn(ji,jj,jl) = Gstari * ( rn_gstar       - Gsum(ji,jj,jl-1) ) *  &
+                        &                        ( 2._wp - ( Gsum(ji,jj,jl-1) + rn_gstar       ) * Gstari )
+                  ELSE
+                     athorn(ji,jj,jl) = 0._wp
+                  ENDIF
+               END DO
+            END DO
+         END DO
+
+      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007)
+         !                        
+         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array
+         DO jl = -1, jpl
+            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
+         END DO
+         DO jl = 0, jpl
+             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
+         END DO
+         !
+      ENDIF
+
+      IF( ln_rafting ) THEN      ! Ridging and rafting ice participation functions
+         !
+         DO jl = 1, jpl
+            DO jj = 1, jpj 
+               DO ji = 1, jpi
+                  zdummy           = TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) )
+                  aridge(ji,jj,jl) = ( 1._wp + zdummy ) * 0.5_wp * athorn(ji,jj,jl)
+                  araft (ji,jj,jl) = ( 1._wp - zdummy ) * 0.5_wp * athorn(ji,jj,jl)
+               END DO
+            END DO
+         END DO
+
+      ELSE
+         !
+         DO jl = 1, jpl
+            aridge(:,:,jl) = athorn(:,:,jl)
+         END DO
+         !
+      ENDIF
+
+      !-----------------------------------------------------------------
+      ! 2) Transfer function
+      !-----------------------------------------------------------------
+      ! Compute max and min ridged ice thickness for each ridging category.
+      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
+      ! 
+      ! This parameterization is a modified version of Hibler (1980).
+      ! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
+      !  and for very thick ridging ice must be >= krdgmin*hi
+      !
+      ! The minimum ridging thickness, hrmin, is equal to 2*hi 
+      !  (i.e., rafting) and for very thick ridging ice is
+      !  constrained by hrmin <= (hrmean + hi)/2.
+      ! 
+      ! The maximum ridging thickness, hrmax, is determined by
+      !  hrmean and hrmin.
+      !
+      ! These modifications have the effect of reducing the ice strength
+      ! (relative to the Hibler formulation) when very thick ice is
+      ! ridging.
+      !
+      ! aksum = net area removed/ total area removed
+      ! where total area removed = area of ice that ridges
+      !         net area removed = total area removed - area of new ridges
+      !-----------------------------------------------------------------
+
+      aksum(:,:) = athorn(:,:,0)
+      ! Transfer function
+      DO jl = 1, jpl !all categories have a specific transfer function
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               
+               IF( athorn(ji,jj,jl) > 0._wp ) THEN
+                  hrmean          = MAX( SQRT( rn_hstar * ht_i(ji,jj,jl) ), ht_i(ji,jj,jl) * krdgmin )
+                  hrmin(ji,jj,jl) = MIN( 2._wp * ht_i(ji,jj,jl), 0.5_wp * ( hrmean + ht_i(ji,jj,jl) ) )
+                  hrmax(ji,jj,jl) = 2._wp * hrmean - hrmin(ji,jj,jl)
+                  hraft(ji,jj,jl) = ht_i(ji,jj,jl) / kraft
+                  krdg(ji,jj,jl)  = ht_i(ji,jj,jl) / MAX( hrmean, epsi20 )
+
+                  ! Normalization factor : aksum, ensures mass conservation
+                  aksum(ji,jj) = aksum(ji,jj) + aridge(ji,jj,jl) * ( 1._wp - krdg(ji,jj,jl) )    &
+                     &                        + araft (ji,jj,jl) * ( 1._wp - kraft          )
+
+               ELSE
+                  hrmin(ji,jj,jl)  = 0._wp 
+                  hrmax(ji,jj,jl)  = 0._wp 
+                  hraft(ji,jj,jl)  = 0._wp 
+                  krdg (ji,jj,jl)  = 1._wp
+               ENDIF
+
+            END DO
+         END DO
+      END DO
+      !
+      CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
+      !
+   END SUBROUTINE lim_itd_me_ridgeprep
+
+
+   SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross )
+      !!----------------------------------------------------------------------
+      !!                ***  ROUTINE lim_itd_me_icestrength ***
+      !!
+      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
+      !!
+      !! ** Method  :   Remove area, volume, and energy from each ridging category
+      !!              and add to thicker ice categories.
+      !!----------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges
+      !
+      CHARACTER (len=80) ::   fieldid   ! field identifier
+      !
+      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
+      INTEGER ::   ij                ! horizontal index, combines i and j loops
+      INTEGER ::   icells            ! number of cells with a_i > puny
+      REAL(wp) ::   hL, hR, farea    ! left and right limits of integration
+
+      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices
+      REAL(wp), POINTER, DIMENSION(:) ::   zswitch, fvol   ! new ridge volume going to n2
+
+      REAL(wp), POINTER, DIMENSION(:) ::   afrac            ! fraction of category area ridged 
+      REAL(wp), POINTER, DIMENSION(:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges
+      REAL(wp), POINTER, DIMENSION(:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice
+      REAL(wp), POINTER, DIMENSION(:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
+
+      REAL(wp), POINTER, DIMENSION(:) ::   vrdg1   ! volume of ice ridged
+      REAL(wp), POINTER, DIMENSION(:) ::   vrdg2   ! volume of new ridges
+      REAL(wp), POINTER, DIMENSION(:) ::   vsw     ! volume of seawater trapped into ridges
+      REAL(wp), POINTER, DIMENSION(:) ::   srdg1   ! sal*volume of ice ridged
+      REAL(wp), POINTER, DIMENSION(:) ::   srdg2   ! sal*volume of new ridges
+      REAL(wp), POINTER, DIMENSION(:) ::   smsw    ! sal*volume of water trapped into ridges
+      REAL(wp), POINTER, DIMENSION(:) ::   oirdg1, oirdg2   ! ice age of ice ridged
+
+      REAL(wp), POINTER, DIMENSION(:) ::   afrft            ! fraction of category area rafted
+      REAL(wp), POINTER, DIMENSION(:) ::   arft1 , arft2    ! area of ice rafted and new rafted zone
+      REAL(wp), POINTER, DIMENSION(:) ::   virft , vsrft    ! ice & snow volume of rafting ice
+      REAL(wp), POINTER, DIMENSION(:) ::   esrft , smrft    ! snow energy & salinity of rafting ice
+      REAL(wp), POINTER, DIMENSION(:) ::   oirft1, oirft2   ! ice age of ice rafted
+
+      REAL(wp), POINTER, DIMENSION(:,:) ::   eirft      ! ice energy of rafting ice
+      REAL(wp), POINTER, DIMENSION(:,:) ::   erdg1      ! enth*volume of ice ridged
+      REAL(wp), POINTER, DIMENSION(:,:) ::   erdg2      ! enth*volume of new ridges
+      REAL(wp), POINTER, DIMENSION(:,:) ::   ersw       ! enth of water trapped into ridges
+      !!----------------------------------------------------------------------
+
+      CALL wrk_alloc( jpij,        indxi, indxj )
+      CALL wrk_alloc( jpij,        zswitch, fvol )
+      CALL wrk_alloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
+      CALL wrk_alloc( jpij,        vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 )
+      CALL wrk_alloc( jpij,        afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
+      CALL wrk_alloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw )
+
+      !-------------------------------------------------------------------------------
+      ! 1) Compute change in open water area due to closing and opening.
+      !-------------------------------------------------------------------------------
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ato_i(ji,jj) = MAX( 0._wp, ato_i(ji,jj) +  &
+               &                     ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice )
+         END DO
+      END DO
+
+      !-----------------------------------------------------------------
+      ! 3) Pump everything from ice which is being ridged / rafted
+      !-----------------------------------------------------------------
+      ! Compute the area, volume, and energy of ice ridging in each
+      ! category, along with the area of the resulting ridge.
+
+      DO jl1 = 1, jpl !jl1 describes the ridging category
+
+         !------------------------------------------------
+         ! 3.1) Identify grid cells with nonzero ridging
+         !------------------------------------------------
+         icells = 0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF( athorn(ji,jj,jl1) > 0._wp .AND. closing_gross(ji,jj) > 0._wp ) THEN
+                  icells = icells + 1
+                  indxi(icells) = ji
+                  indxj(icells) = jj
+               ENDIF
+            END DO
+         END DO
+
+         DO ij = 1, icells
+            ji = indxi(ij) ; jj = indxj(ij)
+
+            !--------------------------------------------------------------------
+            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
+            !--------------------------------------------------------------------
+            ardg1(ij) = aridge(ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
+            arft1(ij) = araft (ji,jj,jl1) * closing_gross(ji,jj) * rdt_ice
+
+            !---------------------------------------------------------------
+            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1 
+            !---------------------------------------------------------------
+            afrac(ij) = ardg1(ij) / a_i(ji,jj,jl1) !ridging
+            afrft(ij) = arft1(ij) / a_i(ji,jj,jl1) !rafting
+            ardg2(ij) = ardg1(ij) * krdg(ji,jj,jl1)
+            arft2(ij) = arft1(ij) * kraft
+
+            !--------------------------------------------------------------------------
+            ! 3.4) Subtract area, volume, and energy from ridging 
+            !     / rafting category n1.
+            !--------------------------------------------------------------------------
+            vrdg1(ij) = v_i(ji,jj,jl1) * afrac(ij)
+            vrdg2(ij) = vrdg1(ij) * ( 1. + rn_por_rdg )
+            vsw  (ij) = vrdg1(ij) * rn_por_rdg
+
+            vsrdg (ij) = v_s  (ji,jj,  jl1) * afrac(ij)
+            esrdg (ij) = e_s  (ji,jj,1,jl1) * afrac(ij)
+            srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij)
+            oirdg1(ij) = oa_i (ji,jj,  jl1) * afrac(ij)
+            oirdg2(ij) = oa_i (ji,jj,  jl1) * afrac(ij) * krdg(ji,jj,jl1) 
+
+            ! rafting volumes, heat contents ...
+            virft (ij) = v_i  (ji,jj,  jl1) * afrft(ij)
+            vsrft (ij) = v_s  (ji,jj,  jl1) * afrft(ij)
+            esrft (ij) = e_s  (ji,jj,1,jl1) * afrft(ij)
+            smrft (ij) = smv_i(ji,jj,  jl1) * afrft(ij) 
+            oirft1(ij) = oa_i (ji,jj,  jl1) * afrft(ij) 
+            oirft2(ij) = oa_i (ji,jj,  jl1) * afrft(ij) * kraft 
+
+            !-----------------------------------------------------------------
+            ! 3.5) Compute properties of new ridges
+            !-----------------------------------------------------------------
+            smsw(ij)  = vsw(ij) * sss_m(ji,jj)                   ! salt content of seawater frozen in voids
+            srdg2(ij) = srdg1(ij) + smsw(ij)                     ! salt content of new ridge
+            
+            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ij) * rhoic * r1_rdtice
+            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ij) * rhoic * r1_rdtice   ! increase in ice volume due to seawater frozen in voids
+            
+            !------------------------------------------            
+            ! 3.7 Put the snow somewhere in the ocean
+            !------------------------------------------            
+            !  Place part of the snow lost by ridging into the ocean. 
+            !  Note that esrdg > 0; the ocean must cool to melt snow.
+            !  If the ocean temp = Tf already, new ice must grow.
+            !  During the next time step, thermo_rates will determine whether
+            !  the ocean cools or new ice grows.
+            wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( rhosn * vsrdg(ij) * ( 1._wp - rn_fsnowrdg )   & 
+               &                              + rhosn * vsrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice  ! fresh water source for ocean
+
+            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ( - esrdg(ij) * ( 1._wp - rn_fsnowrdg )         & 
+               &                                - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
+
+            !-----------------------------------------------------------------
+            ! 3.8 Compute quantities used to apportion ice among categories
+            ! in the n2 loop below
+            !-----------------------------------------------------------------
+            dhr (ij) = 1._wp / ( hrmax(ji,jj,jl1)                    - hrmin(ji,jj,jl1)                    )
+            dhr2(ij) = 1._wp / ( hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1) )
+
+
+            ! update jl1 (removing ridged/rafted area)
+            a_i  (ji,jj,  jl1) = a_i  (ji,jj,  jl1) - ardg1 (ij) - arft1 (ij)
+            v_i  (ji,jj,  jl1) = v_i  (ji,jj,  jl1) - vrdg1 (ij) - virft (ij)
+            v_s  (ji,jj,  jl1) = v_s  (ji,jj,  jl1) - vsrdg (ij) - vsrft (ij)
+            e_s  (ji,jj,1,jl1) = e_s  (ji,jj,1,jl1) - esrdg (ij) - esrft (ij)
+            smv_i(ji,jj,  jl1) = smv_i(ji,jj,  jl1) - srdg1 (ij) - smrft (ij)
+            oa_i (ji,jj,  jl1) = oa_i (ji,jj,  jl1) - oirdg1(ij) - oirft1(ij)
+
+         END DO
+
+         !--------------------------------------------------------------------
+         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
+         !      compute ridged ice enthalpy 
+         !--------------------------------------------------------------------
+         DO jk = 1, nlay_i
+            DO ij = 1, icells
+               ji = indxi(ij) ; jj = indxj(ij)
+               ! heat content of ridged ice
+               erdg1(ij,jk) = e_i(ji,jj,jk,jl1) * afrac(ij) 
+               eirft(ij,jk) = e_i(ji,jj,jk,jl1) * afrft(ij)               
+               
+               ! enthalpy of the trapped seawater (J/m2, >0)
+               ! clem: if sst>0, then ersw <0 (is that possible?)
+               ersw(ij,jk)  = - rhoic * vsw(ij) * rcp * sst_m(ji,jj) * r1_nlay_i
+
+               ! heat flux to the ocean
+               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ij,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux 
+
+               ! it is added to sea ice because the sign convention is the opposite of the sign convention for the ocean
+               erdg2(ij,jk) = erdg1(ij,jk) + ersw(ij,jk)
+
+               ! update jl1
+               e_i  (ji,jj,jk,jl1) = e_i(ji,jj,jk,jl1) - erdg1(ij,jk) - eirft(ij,jk)
+
+            END DO
+         END DO
+
+         !-------------------------------------------------------------------------------
+         ! 4) Add area, volume, and energy of new ridge to each category jl2
+         !-------------------------------------------------------------------------------
+         DO jl2  = 1, jpl 
+            ! over categories to which ridged/rafted ice is transferred
+            DO ij = 1, icells
+               ji = indxi(ij) ; jj = indxj(ij)
+
+               ! Compute the fraction of ridged ice area and volume going to thickness category jl2.
+               IF( hrmin(ji,jj,jl1) <= hi_max(jl2) .AND. hrmax(ji,jj,jl1) > hi_max(jl2-1) ) THEN
+                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
+                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
+                  farea    = ( hR      - hL      ) * dhr(ij) 
+                  fvol(ij) = ( hR * hR - hL * hL ) * dhr2(ij)
+               ELSE
+                  farea    = 0._wp 
+                  fvol(ij) = 0._wp                  
+               ENDIF
+
+               ! Compute the fraction of rafted ice area and volume going to thickness category jl2
+               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
+                  zswitch(ij) = 1._wp
+               ELSE
+                  zswitch(ij) = 0._wp                  
+               ENDIF
+
+               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ( ardg2 (ij) * farea    + arft2 (ij) * zswitch(ij) )
+               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + ( oirdg2(ij) * farea    + oirft2(ij) * zswitch(ij) )
+               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + ( vrdg2 (ij) * fvol(ij) + virft (ij) * zswitch(ij) )
+               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + ( srdg2 (ij) * fvol(ij) + smrft (ij) * zswitch(ij) )
+               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + ( vsrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
+                  &                                        vsrft (ij) * rn_fsnowrft * zswitch(ij) )
+               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
+                  &                                        esrft (ij) * rn_fsnowrft * zswitch(ij) )
+
+            END DO
+
+            ! Transfer ice energy to category jl2 by ridging
+            DO jk = 1, nlay_i
+               DO ij = 1, icells
+                  ji = indxi(ij) ; jj = indxj(ij)
+                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + erdg2(ij,jk) * fvol(ij) + eirft(ij,jk) * zswitch(ij)                  
+               END DO
+            END DO
+            !
+         END DO ! jl2
+         
+      END DO ! jl1 (deforming categories)
+
+      !
+      CALL wrk_dealloc( jpij,        indxi, indxj )
+      CALL wrk_dealloc( jpij,        zswitch, fvol )
+      CALL wrk_dealloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
+      CALL wrk_dealloc( jpij,        vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 )
+      CALL wrk_dealloc( jpij,        afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
+      CALL wrk_dealloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw )
+      !
+   END SUBROUTINE lim_itd_me_ridgeshift
 
    SUBROUTINE lim_itd_me_icestrength( kstrngth )
@@ -434,5 +787,5 @@
       INTEGER             ::   ksmooth     ! smoothing the resistance to deformation
       INTEGER             ::   numts_rm    ! number of time steps for the P smoothing
-      REAL(wp)            ::   zhi, zp, z1_3  ! local scalars
+      REAL(wp)            ::   zp, z1_3    ! local scalars
       REAL(wp), POINTER, DIMENSION(:,:) ::   zworka   ! temporary array used here
       !!----------------------------------------------------------------------
@@ -459,21 +812,22 @@
                DO ji = 1, jpi
                   !
-                  IF( a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0._wp ) THEN
-                     zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
+                  IF( athorn(ji,jj,jl) > 0._wp ) THEN
                      !----------------------------
                      ! PE loss from deforming ice
                      !----------------------------
-                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * zhi * zhi
+                     strength(ji,jj) = strength(ji,jj) - athorn(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
 
                      !--------------------------
                      ! PE gain from rafting ice
                      !--------------------------
-                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * zhi * zhi
+                     strength(ji,jj) = strength(ji,jj) + 2._wp * araft(ji,jj,jl) * ht_i(ji,jj,jl) * ht_i(ji,jj,jl)
 
                      !----------------------------
                      ! PE gain from ridging ice
                      !----------------------------
-                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) / krdg(ji,jj,jl)     &
-                        * z1_3 * ( hrmax(ji,jj,jl)**2 + hrmin(ji,jj,jl)**2 +  hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )  
+                     strength(ji,jj) = strength(ji,jj) + aridge(ji,jj,jl) * krdg(ji,jj,jl) * z1_3 *  &
+                        &                              ( hrmax(ji,jj,jl) * hrmax(ji,jj,jl) +         &
+                        &                                hrmin(ji,jj,jl) * hrmin(ji,jj,jl) +         &
+                        &                                hrmax(ji,jj,jl) * hrmin(ji,jj,jl) )  
                         !!(a**3-b**3)/(a-b) = a*a+ab+b*b                      
                   ENDIF
@@ -497,5 +851,4 @@
          !
       ENDIF                     ! kstrngth
-
       !
       !------------------------------------------------------------------------------!
@@ -503,7 +856,5 @@
       !------------------------------------------------------------------------------!
       ! CAN BE REMOVED
-      !
       IF( ln_icestr_bvf ) THEN
-
          DO jj = 1, jpj
             DO ji = 1, jpi
@@ -511,7 +862,5 @@
             END DO
          END DO
-
       ENDIF
-
       !
       !------------------------------------------------------------------------------!
@@ -558,5 +907,4 @@
       IF ( ksmooth == 2 ) THEN
 
-
          CALL lbc_lnk( strength, 'T', 1. )
 
@@ -565,6 +913,6 @@
                IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN 
                   numts_rm = 1 ! number of time steps for the running mean
-                  IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
-                  IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
+                  IF ( strp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
+                  IF ( strp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1
                   zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
                   strp2(ji,jj) = strp1(ji,jj)
@@ -583,616 +931,4 @@
       !
    END SUBROUTINE lim_itd_me_icestrength
-
-
-   SUBROUTINE lim_itd_me_ridgeprep
-      !!---------------------------------------------------------------------!
-      !!                ***  ROUTINE lim_itd_me_ridgeprep ***
-      !!
-      !! ** Purpose :   preparation for ridging and strength calculations
-      !!
-      !! ** Method  :   Compute the thickness distribution of the ice and open water 
-      !!              participating in ridging and of the resulting ridges.
-      !!---------------------------------------------------------------------!
-      INTEGER ::   ji,jj, jl    ! dummy loop indices
-      REAL(wp) ::   Gstari, astari, zhi, hrmean, zdummy   ! local scalar
-      REAL(wp), POINTER, DIMENSION(:,:)   ::   zworka    ! temporary array used here
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   Gsum      ! Gsum(n) = sum of areas in categories 0 to n
-      !------------------------------------------------------------------------------!
-
-      CALL wrk_alloc( jpi,jpj, zworka )
-      CALL wrk_alloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
-
-      Gstari     = 1.0/rn_gstar    
-      astari     = 1.0/rn_astar    
-      aksum(:,:)    = 0.0
-      athorn(:,:,:) = 0.0
-      aridge(:,:,:) = 0.0
-      araft (:,:,:) = 0.0
-      hrmin(:,:,:)  = 0.0 
-      hrmax(:,:,:)  = 0.0 
-      hraft(:,:,:)  = 0.0 
-      krdg (:,:,:)  = 1.0
-
-      !     ! Zero out categories with very small areas
-      CALL lim_var_zapsmall
-
-      !------------------------------------------------------------------------------!
-      ! 1) Participation function 
-      !------------------------------------------------------------------------------!
-
-      ! Compute total area of ice plus open water.
-      ! This is in general not equal to one because of divergence during transport
-      asum(:,:) = ato_i(:,:)
-      DO jl = 1, jpl
-         asum(:,:) = asum(:,:) + a_i(:,:,jl)
-      END DO
-
-      ! Compute cumulative thickness distribution function
-      ! Compute the cumulative thickness distribution function Gsum,
-      ! where Gsum(n) is the fractional area in categories 0 to n.
-      ! initial value (in h = 0) equals open water area
-
-      Gsum(:,:,-1) = 0._wp
-      Gsum(:,:,0 ) = ato_i(:,:)
-
-      ! for each value of h, you have to add ice concentration then
-      DO jl = 1, jpl
-         Gsum(:,:,jl) = Gsum(:,:,jl-1) + a_i(:,:,jl)
-      END DO
-
-      ! Normalize the cumulative distribution to 1
-      zworka(:,:) = 1._wp / Gsum(:,:,jpl)
-      DO jl = 0, jpl
-         Gsum(:,:,jl) = Gsum(:,:,jl) * zworka(:,:)
-      END DO
-
-      ! 1.3 Compute participation function a(h) = b(h).g(h) (athorn)
-      !--------------------------------------------------------------------------------------------------
-      ! Compute the participation function athorn; this is analogous to
-      ! a(h) = b(h)g(h) as defined in Thorndike et al. (1975).
-      ! area lost from category n due to ridging/closing
-      ! athorn(n)   = total area lost due to ridging/closing
-      ! assume b(h) = (2/Gstar) * (1 - G(h)/Gstar). 
-      !
-      ! The expressions for athorn are found by integrating b(h)g(h) between
-      ! the category boundaries.
-      !-----------------------------------------------------------------
-
-      IF( nn_partfun == 0 ) THEN       !--- Linear formulation (Thorndike et al., 1975)
-         DO jl = 0, jpl    
-            DO jj = 1, jpj 
-               DO ji = 1, jpi
-                  IF( Gsum(ji,jj,jl) < rn_gstar) THEN
-                     athorn(ji,jj,jl) = Gstari * ( Gsum(ji,jj,jl) - Gsum(ji,jj,jl-1) ) * &
-                        &                        ( 2.0 - (Gsum(ji,jj,jl-1) + Gsum(ji,jj,jl) ) * Gstari )
-                  ELSEIF (Gsum(ji,jj,jl-1) < rn_gstar) THEN
-                     athorn(ji,jj,jl) = Gstari * ( rn_gstar - Gsum(ji,jj,jl-1) ) *  &
-                        &                        ( 2.0 - ( Gsum(ji,jj,jl-1) + rn_gstar ) * Gstari )
-                  ELSE
-                     athorn(ji,jj,jl) = 0.0
-                  ENDIF
-               END DO
-            END DO
-         END DO
-
-      ELSE                             !--- Exponential, more stable formulation (Lipscomb et al, 2007)
-         !                        
-         zdummy = 1._wp / ( 1._wp - EXP(-astari) )        ! precompute exponential terms using Gsum as a work array
-         DO jl = -1, jpl
-            Gsum(:,:,jl) = EXP( -Gsum(:,:,jl) * astari ) * zdummy
-         END DO
-         DO jl = 0, jpl
-             athorn(:,:,jl) = Gsum(:,:,jl-1) - Gsum(:,:,jl)
-         END DO
-         !
-      ENDIF
-
-      IF( ln_rafting ) THEN      ! Ridging and rafting ice participation functions
-         !
-         DO jl = 1, jpl
-            DO jj = 1, jpj 
-               DO ji = 1, jpi
-                  IF ( athorn(ji,jj,jl) > 0._wp ) THEN
-!!gm  TANH( -X ) = - TANH( X )  so can be computed only 1 time....
-                     aridge(ji,jj,jl) = ( TANH (  rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
-                     araft (ji,jj,jl) = ( TANH ( -rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
-                     IF ( araft(ji,jj,jl) < epsi06 )   araft(ji,jj,jl)  = 0._wp
-                     aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )
-                  ENDIF
-               END DO
-            END DO
-         END DO
-
-      ELSE
-         !
-         DO jl = 1, jpl
-            aridge(:,:,jl) = athorn(:,:,jl)
-         END DO
-         !
-      ENDIF
-
-      IF( ln_rafting ) THEN
-
-         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 .AND. lwp ) THEN
-            DO jl = 1, jpl
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN
-                        WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
-                        WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
-                        WRITE(numout,*) ' lat, lon   : ', gphit(ji,jj), glamt(ji,jj)
-                        WRITE(numout,*) ' aridge     : ', aridge(ji,jj,1:jpl)
-                        WRITE(numout,*) ' araft      : ', araft(ji,jj,1:jpl)
-                        WRITE(numout,*) ' athorn     : ', athorn(ji,jj,1:jpl)
-                     ENDIF
-                  END DO
-               END DO
-            END DO
-         ENDIF
-
-      ENDIF
-
-      !-----------------------------------------------------------------
-      ! 2) Transfer function
-      !-----------------------------------------------------------------
-      ! Compute max and min ridged ice thickness for each ridging category.
-      ! Assume ridged ice is uniformly distributed between hrmin and hrmax.
-      ! 
-      ! This parameterization is a modified version of Hibler (1980).
-      ! The mean ridging thickness, hrmean, is proportional to hi^(0.5)
-      !  and for very thick ridging ice must be >= krdgmin*hi
-      !
-      ! The minimum ridging thickness, hrmin, is equal to 2*hi 
-      !  (i.e., rafting) and for very thick ridging ice is
-      !  constrained by hrmin <= (hrmean + hi)/2.
-      ! 
-      ! The maximum ridging thickness, hrmax, is determined by
-      !  hrmean and hrmin.
-      !
-      ! These modifications have the effect of reducing the ice strength
-      ! (relative to the Hibler formulation) when very thick ice is
-      ! ridging.
-      !
-      ! aksum = net area removed/ total area removed
-      ! where total area removed = area of ice that ridges
-      !         net area removed = total area removed - area of new ridges
-      !-----------------------------------------------------------------
-
-      ! Transfer function
-      DO jl = 1, jpl !all categories have a specific transfer function
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-
-               IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN
-                  zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
-                  hrmean          = MAX(SQRT(rn_hstar*zhi), zhi*krdgmin)
-                  hrmin(ji,jj,jl) = MIN(2.0*zhi, 0.5*(hrmean + zhi))
-                  hrmax(ji,jj,jl) = 2.0*hrmean - hrmin(ji,jj,jl)
-                  hraft(ji,jj,jl) = kraft*zhi
-                  krdg(ji,jj,jl)  = hrmean / zhi
-               ELSE
-                  hraft(ji,jj,jl) = 0.0
-                  hrmin(ji,jj,jl) = 0.0 
-                  hrmax(ji,jj,jl) = 0.0 
-                  krdg (ji,jj,jl) = 1.0
-               ENDIF
-
-            END DO
-         END DO
-      END DO
-
-      ! Normalization factor : aksum, ensures mass conservation
-      aksum(:,:) = athorn(:,:,0)
-      DO jl = 1, jpl 
-         aksum(:,:)    = aksum(:,:) + aridge(:,:,jl) * ( 1._wp - 1._wp / krdg(:,:,jl) )    &
-            &                       + araft (:,:,jl) * ( 1._wp - 1._wp / kraft        )
-      END DO
-      !
-      CALL wrk_dealloc( jpi,jpj, zworka )
-      CALL wrk_dealloc( jpi,jpj,jpl+2, Gsum, kkstart = -1 )
-      !
-   END SUBROUTINE lim_itd_me_ridgeprep
-
-
-   SUBROUTINE lim_itd_me_ridgeshift( opning, closing_gross, msnow_mlt, esnow_mlt )
-      !!----------------------------------------------------------------------
-      !!                ***  ROUTINE lim_itd_me_icestrength ***
-      !!
-      !! ** Purpose :   shift ridging ice among thickness categories of ice thickness
-      !!
-      !! ** Method  :   Remove area, volume, and energy from each ridging category
-      !!              and add to thicker ice categories.
-      !!----------------------------------------------------------------------
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   opning         ! rate of opening due to divergence/shear
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   closing_gross  ! rate at which area removed, excluding area of new ridges
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   msnow_mlt      ! mass of snow added to ocean (kg m-2)
-      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   esnow_mlt      ! energy needed to melt snow in ocean (J m-2)
-      !
-      CHARACTER (len=80) ::   fieldid   ! field identifier
-      LOGICAL, PARAMETER ::   l_conservation_check = .true.  ! if true, check conservation (useful for debugging)
-      !
-      INTEGER ::   ji, jj, jl, jl1, jl2, jk   ! dummy loop indices
-      INTEGER ::   ij                ! horizontal index, combines i and j loops
-      INTEGER ::   icells            ! number of cells with aicen > puny
-      REAL(wp) ::   hL, hR, farea, ztmelts    ! left and right limits of integration
-
-      INTEGER , POINTER, DIMENSION(:) ::   indxi, indxj   ! compressed indices
-
-      REAL(wp), POINTER, DIMENSION(:,:) ::   vice_init, vice_final   ! ice volume summed over categories
-      REAL(wp), POINTER, DIMENSION(:,:) ::   eice_init, eice_final   ! ice energy summed over layers
-
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   aicen_init, vicen_init   ! ice  area    & volume before ridging
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   vsnwn_init, esnwn_init   ! snow volume  & energy before ridging
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   smv_i_init, oa_i_init    ! ice salinity & age    before ridging
-
-      REAL(wp), POINTER, DIMENSION(:,:,:,:) ::   eicen_init        ! ice energy before ridging
-
-      REAL(wp), POINTER, DIMENSION(:,:) ::   afrac , fvol     ! fraction of category area ridged & new ridge volume going to n2
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges
-      REAL(wp), POINTER, DIMENSION(:,:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice
-      REAL(wp), POINTER, DIMENSION(:,:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
-
-      REAL(wp), POINTER, DIMENSION(:,:) ::   vrdg1   ! volume of ice ridged
-      REAL(wp), POINTER, DIMENSION(:,:) ::   vrdg2   ! volume of new ridges
-      REAL(wp), POINTER, DIMENSION(:,:) ::   vsw     ! volume of seawater trapped into ridges
-      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg1   ! sal*volume of ice ridged
-      REAL(wp), POINTER, DIMENSION(:,:) ::   srdg2   ! sal*volume of new ridges
-      REAL(wp), POINTER, DIMENSION(:,:) ::   smsw    ! sal*volume of water trapped into ridges
-      REAL(wp), POINTER, DIMENSION(:,:) ::   oirdg1, oirdg2   ! ice age of ice ridged
-
-      REAL(wp), POINTER, DIMENSION(:,:) ::   afrft            ! fraction of category area rafted
-      REAL(wp), POINTER, DIMENSION(:,:) ::   arft1 , arft2    ! area of ice rafted and new rafted zone
-      REAL(wp), POINTER, DIMENSION(:,:) ::   virft , vsrft    ! ice & snow volume of rafting ice
-      REAL(wp), POINTER, DIMENSION(:,:) ::   esrft , smrft    ! snow energy & salinity of rafting ice
-      REAL(wp), POINTER, DIMENSION(:,:) ::   oirft1, oirft2   ! ice age of ice rafted
-
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   eirft      ! ice energy of rafting ice
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   erdg1      ! enth*volume of ice ridged
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   erdg2      ! enth*volume of new ridges
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ersw       ! enth of water trapped into ridges
-      !!----------------------------------------------------------------------
-
-      CALL wrk_alloc( (jpi+1)*(jpj+1),       indxi, indxj )
-      CALL wrk_alloc( jpi, jpj,              vice_init, vice_final, eice_init, eice_final )
-      CALL wrk_alloc( jpi, jpj,              afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
-      CALL wrk_alloc( jpi, jpj,              vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 )
-      CALL wrk_alloc( jpi, jpj,              afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
-      CALL wrk_alloc( jpi, jpj, jpl,         aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
-      CALL wrk_alloc( jpi, jpj, nlay_i,      eirft, erdg1, erdg2, ersw )
-      CALL wrk_alloc( jpi, jpj, nlay_i, jpl, eicen_init )
-
-      ! Conservation check
-      eice_init(:,:) = 0._wp
-
-      IF( con_i ) THEN
-         CALL lim_column_sum        (jpl,    v_i,       vice_init )
-         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_init )
-         DO ji = mi0(iiceprt), mi1(iiceprt)
-            DO jj = mj0(jiceprt), mj1(jiceprt)
-               WRITE(numout,*) ' vice_init  : ', vice_init(ji,jj)
-               WRITE(numout,*) ' eice_init  : ', eice_init(ji,jj)
-            END DO
-         END DO
-      ENDIF
-
-      !-------------------------------------------------------------------------------
-      ! 1) Compute change in open water area due to closing and opening.
-      !-------------------------------------------------------------------------------
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ato_i(ji,jj) = ato_i(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice        &
-               &                        + opning(ji,jj)                          * rdt_ice
-            IF    ( ato_i(ji,jj) < -epsi10 ) THEN    ! there is a bug
-               IF(lwp)   WRITE(numout,*) 'Ridging error: ato_i < 0 -- ato_i : ',ato_i(ji,jj)
-            ELSEIF( ato_i(ji,jj) < 0._wp   ) THEN    ! roundoff error
-               ato_i(ji,jj) = 0._wp
-            ENDIF
-         END DO
-      END DO
-
-      !-----------------------------------------------------------------
-      ! 2) Save initial state variables
-      !-----------------------------------------------------------------
-      aicen_init(:,:,:)   = a_i  (:,:,:)
-      vicen_init(:,:,:)   = v_i  (:,:,:)
-      vsnwn_init(:,:,:)   = v_s  (:,:,:)
-      smv_i_init(:,:,:)   = smv_i(:,:,:)
-      esnwn_init(:,:,:)   = e_s  (:,:,1,:)
-      eicen_init(:,:,:,:) = e_i  (:,:,:,:)
-      oa_i_init (:,:,:)   = oa_i (:,:,:)
-
-      !
-      !-----------------------------------------------------------------
-      ! 3) Pump everything from ice which is being ridged / rafted
-      !-----------------------------------------------------------------
-      ! Compute the area, volume, and energy of ice ridging in each
-      ! category, along with the area of the resulting ridge.
-
-      DO jl1 = 1, jpl !jl1 describes the ridging category
-
-         !------------------------------------------------
-         ! 3.1) Identify grid cells with nonzero ridging
-         !------------------------------------------------
-
-         icells = 0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF( aicen_init(ji,jj,jl1) > epsi10 .AND. athorn(ji,jj,jl1) > 0._wp  &
-                  &   .AND. closing_gross(ji,jj) > 0._wp ) THEN
-                  icells = icells + 1
-                  indxi(icells) = ji
-                  indxj(icells) = jj
-               ENDIF
-            END DO
-         END DO
-
-         DO ij = 1, icells
-            ji = indxi(ij)
-            jj = indxj(ij)
-
-            !--------------------------------------------------------------------
-            ! 3.2) Compute area of ridging ice (ardg1) and of new ridge (ardg2)
-            !--------------------------------------------------------------------
-
-            ardg1(ji,jj) = aridge(ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
-            arft1(ji,jj) = araft (ji,jj,jl1)*closing_gross(ji,jj)*rdt_ice
-            ardg2(ji,jj) = ardg1(ji,jj) / krdg(ji,jj,jl1)
-            arft2(ji,jj) = arft1(ji,jj) / kraft
-
-            !---------------------------------------------------------------
-            ! 3.3) Compute ridging /rafting fractions, make sure afrac <=1 
-            !---------------------------------------------------------------
-
-            afrac(ji,jj) = ardg1(ji,jj) / aicen_init(ji,jj,jl1) !ridging
-            afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting
-
-            IF( afrac(ji,jj) > kamax + epsi10 ) THEN  ! there is a bug
-               IF(lwp)   WRITE(numout,*) ' ardg > a_i -- ardg, aicen_init : ', ardg1(ji,jj), aicen_init(ji,jj,jl1)
-            ELSEIF( afrac(ji,jj) > kamax ) THEN       ! roundoff error
-               afrac(ji,jj) = kamax
-            ENDIF
-
-            IF( afrft(ji,jj) > kamax + epsi10 ) THEN ! there is a bug
-               IF(lwp)   WRITE(numout,*) ' arft > a_i -- arft, aicen_init : ', arft1(ji,jj), aicen_init(ji,jj,jl1) 
-            ELSEIF( afrft(ji,jj) > kamax) THEN       ! roundoff error
-               afrft(ji,jj) = kamax
-            ENDIF
-
-            !--------------------------------------------------------------------------
-            ! 3.4) Subtract area, volume, and energy from ridging 
-            !     / rafting category n1.
-            !--------------------------------------------------------------------------
-            vrdg1(ji,jj) = vicen_init(ji,jj,jl1) * afrac(ji,jj)
-            vrdg2(ji,jj) = vrdg1(ji,jj) * ( 1. + rn_por_rdg )
-            vsw  (ji,jj) = vrdg1(ji,jj) * rn_por_rdg
-
-            vsrdg (ji,jj) = vsnwn_init(ji,jj,jl1) * afrac(ji,jj)
-            esrdg (ji,jj) = esnwn_init(ji,jj,jl1) * afrac(ji,jj)
-            srdg1 (ji,jj) = smv_i_init(ji,jj,jl1) * afrac(ji,jj)
-            oirdg1(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj)
-            oirdg2(ji,jj) = oa_i_init (ji,jj,jl1) * afrac(ji,jj) / krdg(ji,jj,jl1) 
-
-            ! rafting volumes, heat contents ...
-            virft (ji,jj) = vicen_init(ji,jj,jl1) * afrft(ji,jj)
-            vsrft (ji,jj) = vsnwn_init(ji,jj,jl1) * afrft(ji,jj)
-            esrft (ji,jj) = esnwn_init(ji,jj,jl1) * afrft(ji,jj)
-            smrft (ji,jj) = smv_i_init(ji,jj,jl1) * afrft(ji,jj) 
-            oirft1(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj) 
-            oirft2(ji,jj) = oa_i_init (ji,jj,jl1) * afrft(ji,jj) / kraft 
-
-            ! substract everything
-            a_i(ji,jj,jl1)   = a_i(ji,jj,jl1)   - ardg1 (ji,jj) - arft1 (ji,jj)
-            v_i(ji,jj,jl1)   = v_i(ji,jj,jl1)   - vrdg1 (ji,jj) - virft (ji,jj)
-            v_s(ji,jj,jl1)   = v_s(ji,jj,jl1)   - vsrdg (ji,jj) - vsrft (ji,jj)
-            e_s(ji,jj,1,jl1) = e_s(ji,jj,1,jl1) - esrdg (ji,jj) - esrft (ji,jj)
-            smv_i(ji,jj,jl1) = smv_i(ji,jj,jl1) - srdg1 (ji,jj) - smrft (ji,jj)
-            oa_i(ji,jj,jl1)  = oa_i(ji,jj,jl1)  - oirdg1(ji,jj) - oirft1(ji,jj)
-
-            !-----------------------------------------------------------------
-            ! 3.5) Compute properties of new ridges
-            !-----------------------------------------------------------------
-            !---------
-            ! Salinity
-            !---------
-            smsw(ji,jj)  = vsw(ji,jj) * sss_m(ji,jj)                      ! salt content of seawater frozen in voids !! MV HC2014
-            srdg2(ji,jj) = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge
-
-            !srdg2(ji,jj) = MIN( rn_simax * vrdg2(ji,jj) , zsrdg2 )         ! impose a maximum salinity
-            
-            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice
-            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice   ! increase in ice volume du to seawater frozen in voids             
-
-            !------------------------------------            
-            ! 3.6 Increment ridging diagnostics
-            !------------------------------------            
-
-            !        jl1 looping 1-jpl
-            !           ij looping 1-icells
-
-            dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)
-            dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)
-            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice
-
-            IF( con_i )   vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)
-
-            !------------------------------------------            
-            ! 3.7 Put the snow somewhere in the ocean
-            !------------------------------------------            
-            !  Place part of the snow lost by ridging into the ocean. 
-            !  Note that esnow_mlt < 0; the ocean must cool to melt snow.
-            !  If the ocean temp = Tf already, new ice must grow.
-            !  During the next time step, thermo_rates will determine whether
-            !  the ocean cools or new ice grows.
-            !        jl1 looping 1-jpl
-            !           ij looping 1-icells
-
-            msnow_mlt(ji,jj) = msnow_mlt(ji,jj) + rhosn*vsrdg(ji,jj)*(1.0-rn_fsnowrdg)   &   ! rafting included
-               &                                + rhosn*vsrft(ji,jj)*(1.0-rn_fsnowrft)
-
-            ! in J/m2 (same as e_s)
-            esnow_mlt(ji,jj) = esnow_mlt(ji,jj) - esrdg(ji,jj)*(1.0-rn_fsnowrdg)         &   !rafting included
-               &                                - esrft(ji,jj)*(1.0-rn_fsnowrft)          
-
-            !-----------------------------------------------------------------
-            ! 3.8 Compute quantities used to apportion ice among categories
-            ! in the n2 loop below
-            !-----------------------------------------------------------------
-
-            !        jl1 looping 1-jpl
-            !           ij looping 1-icells
-
-            dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)
-            dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)
-
-         END DO
-
-         !--------------------------------------------------------------------
-         ! 3.9 Compute ridging ice enthalpy, remove it from ridging ice and
-         !      compute ridged ice enthalpy 
-         !--------------------------------------------------------------------
-         DO jk = 1, nlay_i
-            DO ij = 1, icells
-               ji = indxi(ij)
-               jj = indxj(ij)
-               ! heat content of ridged ice
-               erdg1(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrac(ji,jj) 
-               eirft(ji,jj,jk)      = eicen_init(ji,jj,jk,jl1) * afrft(ji,jj)
-               e_i  (ji,jj,jk,jl1)  = e_i(ji,jj,jk,jl1) - erdg1(ji,jj,jk) - eirft(ji,jj,jk)
-               
-               
-               ! enthalpy of the trapped seawater (J/m2, >0)
-               ! clem: if sst>0, then ersw <0 (is that possible?)
-               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * sst_m(ji,jj) * r1_nlay_i
-
-               ! heat flux to the ocean
-               hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + ersw(ji,jj,jk) * r1_rdtice  ! > 0 [W.m-2] ocean->ice flux 
-
-               ! 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)
-
-            END DO
-         END DO
-
-
-         IF( con_i ) THEN
-            DO jk = 1, nlay_i
-               DO ij = 1, icells
-                  ji = indxi(ij)
-                  jj = indxj(ij)
-                  eice_init(ji,jj) = eice_init(ji,jj) + erdg2(ji,jj,jk) - erdg1(ji,jj,jk)
-               END DO
-            END DO
-         ENDIF
-
-         !-------------------------------------------------------------------------------
-         ! 4) Add area, volume, and energy of new ridge to each category jl2
-         !-------------------------------------------------------------------------------
-         !        jl1 looping 1-jpl
-         DO jl2  = 1, jpl 
-            ! over categories to which ridged ice is transferred
-            DO ij = 1, icells
-               ji = indxi(ij)
-               jj = indxj(ij)
-
-               ! Compute the fraction of ridged ice area and volume going to 
-               ! thickness category jl2.
-               ! Transfer area, volume, and energy accordingly.
-
-               IF( hrmin(ji,jj,jl1) >= hi_max(jl2) .OR. hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN
-                  hL = 0._wp
-                  hR = 0._wp
-               ELSE
-                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
-                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
-               ENDIF
-
-               ! fraction of ridged ice area and volume going to n2
-               farea = ( hR - hL ) / dhr(ji,jj) 
-               fvol(ji,jj) = ( hR*hR - hL*hL ) / dhr2(ji,jj)
-
-               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ardg2 (ji,jj) * farea
-               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj)
-               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg
-               e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * rn_fsnowrdg
-               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + srdg2 (ji,jj) * fvol(ji,jj)
-               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirdg2(ji,jj) * farea
-
-            END DO
-
-            ! Transfer ice energy to category jl2 by ridging
-            DO jk = 1, nlay_i
-               DO ij = 1, icells
-                  ji = indxi(ij)
-                  jj = indxj(ij)
-                  e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + fvol(ji,jj) * erdg2(ji,jj,jk)
-               END DO
-            END DO
-            !
-         END DO                 ! jl2 (new ridges)            
-
-         DO jl2 = 1, jpl 
-
-            DO ij = 1, icells
-               ji = indxi(ij)
-               jj = indxj(ij)
-               ! Compute the fraction of rafted ice area and volume going to 
-               ! thickness category jl2, transfer area, volume, and energy accordingly.
-               !
-               IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
-                  a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + arft2 (ji,jj)
-                  v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + virft (ji,jj)
-                  v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrft (ji,jj) * rn_fsnowrft
-                  e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * rn_fsnowrft
-                  smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)    
-                  oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj)
-               ENDIF
-               !
-            END DO
-
-            ! Transfer rafted ice energy to category jl2 
-            DO jk = 1, nlay_i
-               DO ij = 1, icells
-                  ji = indxi(ij)
-                  jj = indxj(ij)
-                  IF( hraft(ji,jj,jl1) <= hi_max(jl2) .AND. hraft(ji,jj,jl1) > hi_max(jl2-1)  ) THEN
-                     e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)
-                  ENDIF
-               END DO
-            END DO
-
-         END DO
-
-      END DO ! jl1 (deforming categories)
-
-      ! Conservation check
-      IF ( con_i ) THEN
-         CALL lim_column_sum (jpl,   v_i, vice_final)
-         fieldid = ' v_i : limitd_me '
-         CALL lim_cons_check (vice_init, vice_final, 1.0e-6, fieldid) 
-
-         CALL lim_column_sum_energy (jpl, nlay_i,  e_i, eice_final )
-         fieldid = ' e_i : limitd_me '
-         CALL lim_cons_check (eice_init, eice_final, 1.0e-2, fieldid) 
-
-         DO ji = mi0(iiceprt), mi1(iiceprt)
-            DO jj = mj0(jiceprt), mj1(jiceprt)
-               WRITE(numout,*) ' vice_init  : ', vice_init (ji,jj)
-               WRITE(numout,*) ' vice_final : ', vice_final(ji,jj)
-               WRITE(numout,*) ' eice_init  : ', eice_init (ji,jj)
-               WRITE(numout,*) ' eice_final : ', eice_final(ji,jj)
-            END DO
-         END DO
-      ENDIF
-      !
-      CALL wrk_dealloc( (jpi+1)*(jpj+1),        indxi, indxj )
-      CALL wrk_dealloc( jpi, jpj,               vice_init, vice_final, eice_init, eice_final )
-      CALL wrk_dealloc( jpi, jpj,               afrac, fvol , ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
-      CALL wrk_dealloc( jpi, jpj,               vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 )
-      CALL wrk_dealloc( jpi, jpj,               afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
-      CALL wrk_dealloc( jpi, jpj, jpl,          aicen_init, vicen_init, vsnwn_init, esnwn_init, smv_i_init, oa_i_init )
-      CALL wrk_dealloc( jpi, jpj, nlay_i,       eirft, erdg1, erdg2, ersw )
-      CALL wrk_dealloc( jpi, jpj, nlay_i, jpl,  eicen_init )
-      !
-   END SUBROUTINE lim_itd_me_ridgeshift
 
    SUBROUTINE lim_itd_me_init
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 6455)
@@ -94,5 +94,5 @@
       !!              - fr_i    : ice fraction
       !!              - tn_ice  : sea-ice surface temperature
-      !!              - alb_ice : sea-ice albedo (only useful in coupled mode)
+      !!              - alb_ice : sea-ice albedo (recomputed only for coupled mode)
       !!
       !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
@@ -106,8 +106,10 @@
       REAL(wp) ::   zqsr                                           ! New solar flux received by the ocean
       !
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_cs, zalb_os     ! 2D/3D workspace
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb_cs, zalb_os     ! 3D workspace
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zalb                 ! 2D workspace
       !!---------------------------------------------------------------------
 
       ! make calls for heat fluxes before it is modified
+      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
       IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface
       IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface
@@ -118,8 +120,24 @@
       IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   &
          &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) )
-      IF( iom_use('qemp_oce' ) )   CALL iom_put( "qemp_oce"  , qemp_oce(:,:) )  
-      IF( iom_use('qemp_ice' ) )   CALL iom_put( "qemp_ice"  , qemp_ice(:,:) )  
-
-      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
+      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )  
+      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )  
+      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce(:,:) )   ! emp over ocean (taking into account the snow blown away from the ice)
+      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice(:,:) )   ! emp over ice   (taking into account the snow blown away from the ice)
+
+      ! albedo output
+      CALL wrk_alloc( jpi,jpj, zalb )    
+
+      zalb(:,:) = 0._wp
+      WHERE     ( SUM( a_i_b, dim=3 ) <= epsi06 )  ;  zalb(:,:) = 0.066_wp
+      ELSEWHERE                                    ;  zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / SUM( a_i_b, dim=3 )
+      END WHERE
+      IF( iom_use('alb_ice' ) )  CALL iom_put( "alb_ice"  , zalb(:,:) )           ! ice albedo output
+
+      zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + 0.066_wp * ( 1._wp - SUM( a_i_b, dim=3 ) )      
+      IF( iom_use('albedo'  ) )  CALL iom_put( "albedo"  , zalb(:,:) )           ! ice albedo output
+
+      CALL wrk_dealloc( jpi,jpj, zalb )    
+      !
+      
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -140,10 +158,11 @@
             hfx_out(ji,jj) = hfx_out(ji,jj) + zqmass + zqsr
 
-            ! Add the residual from heat diffusion equation (W.m-2)
-            !-------------------------------------------------------
-            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj)
+            ! Add the residual from heat diffusion equation and sublimation (W.m-2)
+            !----------------------------------------------------------------------
+            hfx_out(ji,jj) = hfx_out(ji,jj) + hfx_err_dif(ji,jj) +   &
+               &           ( hfx_sub(ji,jj) - SUM( qevap_ice(ji,jj,:) * a_i_b(ji,jj,:) ) )
 
             ! New qsr and qns used to compute the oceanic heat flux at the next time step
-            !---------------------------------------------------
+            !----------------------------------------------------------------------------
             qsr(ji,jj) = zqsr                                      
             qns(ji,jj) = hfx_out(ji,jj) - zqsr              
@@ -165,7 +184,6 @@
 
             ! mass flux at the ocean/ice interface
-            fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) ) * r1_rdtice  ! F/M mass flux save at least for biogeochemical model
-            emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange)
-            
+            fmmflx(ji,jj) = - ( wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_err_sub(ji,jj) )              ! F/M mass flux save at least for biogeochemical model
+            emp(ji,jj)    = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_err_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange)
          END DO
       END DO
@@ -175,5 +193,5 @@
       !------------------------------------------!
       sfx(:,:) = sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:)   &
-         &     + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:)
+         &     + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:) + sfx_sub(:,:)
 
       !-------------------------------------------------------------!
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 6455)
@@ -461,5 +461,5 @@
 
       DO ji = kideb, kiut
-         zdh_mel = MIN( 0._wp, 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) + dh_i_sub(ji) )
          IF( zdh_mel < 0._wp .AND. a_i_1d(ji) > 0._wp )  THEN
             zvi          = a_i_1d(ji) * ht_i_1d(ji)
@@ -470,5 +470,5 @@
             zda_mel     = rswitch * a_i_1d(ji) * zdh_mel / ( 2._wp * MAX( zhi_bef, epsi20 ) )
             a_i_1d(ji)  = MAX( epsi20, a_i_1d(ji) + zda_mel ) 
-             ! adjust thickness
+            ! adjust thickness
             ht_i_1d(ji) = zvi / a_i_1d(ji)            
             ht_s_1d(ji) = zvs / a_i_1d(ji)            
@@ -514,4 +514,5 @@
          
          CALL tab_2d_1d( nbpb, qprec_ice_1d(1:nbpb), qprec_ice(:,:) , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, qevap_ice_1d(1:nbpb), qevap_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, qsr_ice_1d (1:nbpb), qsr_ice(:,:,jl) , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb), fr1_i0          , jpi, jpj, npb(1:nbpb) )
@@ -543,5 +544,6 @@
          CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res         , jpi, jpj, npb(1:nbpb) )
-         
+         CALL tab_2d_1d( nbpb, sfx_sub_1d (1:nbpb), sfx_sub         , jpi, jpj,npb(1:nbpb) )
+ 
          CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr         , jpi, jpj, npb(1:nbpb) )
@@ -593,5 +595,6 @@
          CALL tab_1d_2d( nbpb, sfx_res       , npb, sfx_res_1d(1:nbpb)   , jpi, jpj )
          CALL tab_1d_2d( nbpb, sfx_bri       , npb, sfx_bri_1d(1:nbpb)   , jpi, jpj )
-         
+         CALL tab_1d_2d( nbpb, sfx_sub       , npb, sfx_sub_1d(1:nbpb)   , jpi, jpj )        
+ 
          CALL tab_1d_2d( nbpb, hfx_thd       , npb, hfx_thd_1d(1:nbpb)   , jpi, jpj )
          CALL tab_1d_2d( nbpb, hfx_spr       , npb, hfx_spr_1d(1:nbpb)   , jpi, jpj )
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 6455)
@@ -74,5 +74,5 @@
 
       REAL(wp) ::   ztmelts             ! local scalar
-      REAL(wp) ::   zfdum       
+      REAL(wp) ::   zdum       
       REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
       REAL(wp) ::   zs_snic      ! snow-ice salinity
@@ -95,4 +95,5 @@
       REAL(wp), POINTER, DIMENSION(:) ::   zq_rema     ! remaining heat at the end of the routine    (J.m-2)
       REAL(wp), POINTER, DIMENSION(:) ::   zf_tt       ! Heat budget to determine melting or freezing(W.m-2)
+      REAL(wp), POINTER, DIMENSION(:) ::   zevap_rema  ! remaining mass flux from sublimation        (kg.m-2)
 
       REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_mel   ! snow melt 
@@ -105,6 +106,4 @@
 
       REAL(wp), POINTER, DIMENSION(:) ::   zqh_i       ! total ice heat content  (J.m-2)
-      REAL(wp), POINTER, DIMENSION(:) ::   zqh_s       ! total snow heat content (J.m-2)
-      REAL(wp), POINTER, DIMENSION(:) ::   zq_s        ! total snow enthalpy     (J.m-3)
       REAL(wp), POINTER, DIMENSION(:) ::   zsnw        ! distribution of snow after wind blowing
 
@@ -122,16 +121,15 @@
       END SELECT
 
-      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw )
-      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
+      CALL wrk_alloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema )
+      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i )
       CALL wrk_alloc( jpij, nlay_i, zdeltah, zh_i )
       CALL wrk_alloc( jpij, nlay_i, icount )
        
-      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp
+      dh_i_surf  (:) = 0._wp ; dh_i_bott  (:) = 0._wp ; dh_snowice(:) = 0._wp ; dh_i_sub(:) = 0._wp
       dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp   
 
       zqprec   (:) = 0._wp ; zq_su    (:) = 0._wp ; zq_bo    (:) = 0._wp ; zf_tt(:) = 0._wp
-      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp
+      zq_rema  (:) = 0._wp ; zsnw     (:) = 0._wp ; zevap_rema(:) = 0._wp ;
       zdh_s_mel(:) = 0._wp ; zdh_s_pre(:) = 0._wp ; zdh_s_sub(:) = 0._wp ; zqh_i(:) = 0._wp
-      zqh_s    (:) = 0._wp ; zq_s     (:) = 0._wp     
 
       zdeltah(:,:) = 0._wp ; zh_i(:,:) = 0._wp       
@@ -159,8 +157,8 @@
       !
       DO ji = kideb, kiut
-         zfdum      = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
+         zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji) 
          zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji) 
 
-         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
+         zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
          zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
       END DO
@@ -187,10 +185,4 @@
       !  2) Computing layer thicknesses and enthalpies.            !
       !------------------------------------------------------------!
-      !
-      DO jk = 1, nlay_s
-         DO ji = kideb, kiut
-            zqh_s(ji) =  zqh_s(ji) + q_s_1d(ji,jk) * ht_s_1d(ji) * r1_nlay_s
-         END DO
-      END DO
       !
       DO jk = 1, nlay_i
@@ -275,17 +267,15 @@
       END DO
 
-      !----------------------
-      ! 3.2 Snow sublimation 
-      !----------------------
+      !------------------------------
+      ! 3.2 Sublimation (part1: snow) 
+      !------------------------------
       ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
       ! clem comment: not counted in mass/heat exchange in limsbc since this is an exchange with atm. (not ocean)
-      ! clem comment: ice should also sublimate
       zdeltah(:,:) = 0._wp
-      ! coupled mode: sublimation is set to 0 (evap_ice = 0) until further notice
-      ! forced  mode: snow thickness change due to sublimation
-      DO ji = kideb, kiut
-         zdh_s_sub(ji)  =  MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
-         ! Heat flux by sublimation [W.m-2], < 0
-         !      sublimate first snow that had fallen, then pre-existing snow
+      DO ji = kideb, kiut
+         zdh_s_sub(ji)  = MAX( - ht_s_1d(ji) , - evap_ice_1d(ji) * r1_rhosn * rdt_ice )
+         ! remaining evap in kg.m-2 (used for ice melting later on)
+         zevap_rema(ji)  = evap_ice_1d(ji) * rdt_ice + zdh_s_sub(ji) * rhosn
+         ! Heat flux by sublimation [W.m-2], < 0 (sublimate first snow that had fallen, then pre-existing snow)
          zdeltah(ji,1)  = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
          hfx_sub_1d(ji) = hfx_sub_1d(ji) + ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * q_s_1d(ji,1)  &
@@ -309,12 +299,10 @@
       !-------------------------------------------
       ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
-      zq_s(:) = 0._wp 
       DO jk = 1, nlay_s
          DO ji = kideb,kiut
-            rswitch       = MAX(  0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 )  )
-            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *                          &
-              &            ( (   zdh_s_pre(ji)             ) * zqprec(ji) +  &
-              &              (   ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
-            zq_s(ji)     =  zq_s(ji) + q_s_1d(ji,jk)
+            rswitch       = MAX( 0._wp , SIGN( 1._wp, ht_s_1d(ji) - epsi20 ) )
+            q_s_1d(ji,jk) = rswitch / MAX( ht_s_1d(ji), epsi20 ) *           &
+              &            ( ( zdh_s_pre(ji)               ) * zqprec(ji) +  &
+              &              ( ht_s_1d(ji) - zdh_s_pre(ji) ) * rhosn * ( cpic * ( rt0 - t_s_1d(ji,jk) ) + lfus ) )
          END DO
       END DO
@@ -370,5 +358,5 @@
                zQm            = zfmdt * zEw                           ! Energy of the melt water sent to the ocean [J/m2, <0]
                
-               ! Contribution to salt flux (clem: using sm_i_1d and not s_i_1d(jk) is ok)
+               ! Contribution to salt flux >0 (clem: using sm_i_1d and not s_i_1d(jk) is ok)
                sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoic * a_i_1d(ji) * zdeltah(ji,jk) * sm_i_1d(ji) * r1_rdtice
                
@@ -383,4 +371,20 @@
                
             END IF
+            ! ----------------------
+            ! Sublimation part2: ice
+            ! ----------------------
+            zdum      = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoic )
+            zdeltah(ji,jk) = zdeltah(ji,jk) + zdum
+            dh_i_sub(ji)  = dh_i_sub(ji) + zdum
+            ! Salt flux > 0 (clem2016: flux is sent to the ocean for simplicity but salt should remain in the ice except if all ice is melted.
+            !                          It must be corrected at some point)
+            sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * sm_i_1d(ji) * r1_rdtice
+            ! Heat flux [W.m-2], < 0
+            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * q_i_1d(ji,jk) * a_i_1d(ji) * r1_rdtice
+            ! Mass flux > 0
+            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhoic * a_i_1d(ji) * zdum * r1_rdtice
+            ! update remaining mass flux
+            zevap_rema(ji)  = zevap_rema(ji) + zdum * rhoic
+            
             ! record which layers have disappeared (for bottom melting) 
             !    => icount=0 : no layer has vanished
@@ -389,5 +393,5 @@
             icount(ji,jk) = NINT( rswitch )
             zh_i(ji,jk)   = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
-
+                        
             ! update heat content (J.m-2) and layer thickness
             qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_1d(ji,jk)
@@ -397,5 +401,11 @@
       ! update ice thickness
       DO ji = kideb, kiut
-         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) )
+         ht_i_1d(ji) =  MAX( 0._wp , ht_i_1d(ji) + dh_i_surf(ji) + dh_i_sub(ji) )
+      END DO
+
+      ! remaining "potential" evap is sent to ocean
+      DO ji = kideb, kiut
+         ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
+         wfx_err_sub(ii,ij) = wfx_err_sub(ii,ij) - zevap_rema(ji) * a_i_1d(ji) * r1_rdtice  ! <=0 (net evap for the ocean in kg.m-2.s-1)
       END DO
 
@@ -686,6 +696,6 @@
       WHERE( ht_i_1d == 0._wp ) a_i_1d = 0._wp
       
-      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw )
-      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
+      CALL wrk_dealloc( jpij, zqprec, zq_su, zq_bo, zf_tt, zq_rema, zsnw, zevap_rema )
+      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i )
       CALL wrk_dealloc( jpij, nlay_i, zdeltah, zh_i )
       CALL wrk_dealloc( jpij, nlay_i, icount )
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 6455)
@@ -75,7 +75,6 @@
       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) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf                     !   -      -
       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
-      LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness
       CHARACTER (len = 15) :: fieldid
 
@@ -108,9 +107,9 @@
       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
+      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                     ! ice-air drag (clem: should be dependent on coupling/forcing used)
       !!-----------------------------------------------------------------------!
 
@@ -143,6 +142,5 @@
       !------------------------------------------------------------------------------!    
       ! hicol is the thickness of new ice formed in open water
-      ! hicol can be either prescribed (frazswi = 0)
-      ! or computed (frazswi = 1)
+      ! hicol can be either prescribed (frazswi = 0) or computed (frazswi = 1)
       ! Frazil ice forms in open water, is transported by wind
       ! accumulates at the edge of the consolidated ice edge
@@ -155,6 +153,8 @@
       zvrel(:,:) = 0._wp
 
-      ! Default new ice thickness 
-      hicol(:,:) = rn_hnewice
+      ! Default new ice thickness
+      WHERE( qlead(:,:) < 0._wp ) ; hicol = rn_hnewice
+      ELSEWHERE                   ; hicol = 0._wp
+      END WHERE
 
       IF( ln_frazil ) THEN
@@ -182,5 +182,5 @@
                      &          +   vtau_ice(ji  ,jj  ) * vmask(ji  ,jj  ,1) ) * 0.5_wp
                   ! Square root of wind stress
-                  ztenagm       =  SQRT( SQRT( ztaux**2 + ztauy**2 ) )
+                  ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )
 
                   !---------------------
@@ -205,41 +205,29 @@
                   zvrel2 = MAX(  ( zvfrx - zvgx ) * ( zvfrx - zvgx )   &
                      &         + ( zvfry - zvgy ) * ( zvfry - zvgy ) , 0.15 * 0.15 )
-                  zvrel(ji,jj)  = SQRT( zvrel2 )
+                  zvrel(ji,jj) = SQRT( zvrel2 )
 
                   !---------------------
                   ! Iterative procedure
                   !---------------------
-                  hicol(ji,jj) = zhicrit + 0.1 
-                  hicol(ji,jj) = zhicrit +   hicol(ji,jj)    &
-                     &                   / ( hicol(ji,jj) * hicol(ji,jj) -  zhicrit * zhicrit ) * ztwogp * zvrel2
-
-!!gm better coding: above: hicol(ji,jj) * hicol(ji,jj) = (zhicrit + 0.1)*(zhicrit + 0.1)
-!!gm                                                   = zhicrit**2 + 0.2*zhicrit +0.01
-!!gm                therefore the 2 lines with hicol can be replaced by 1 line:
-!!gm              hicol(ji,jj) = zhicrit + (zhicrit + 0.1) / ( 0.2 * zhicrit + 0.01 ) * ztwogp * zvrel2
-!!gm further more (zhicrit + 0.1)/(0.2 * zhicrit + 0.01 )*ztwogp can be computed one for all outside the DO loop
+                  hicol(ji,jj) = zhicrit +   ( zhicrit + 0.1 )    &
+                     &                   / ( ( zhicrit + 0.1 ) * ( zhicrit + 0.1 ) -  zhicrit * zhicrit ) * ztwogp * zvrel2
 
                   iter = 1
-                  iterate_frazil = .true.
-
-                  DO WHILE ( iter < 100 .AND. iterate_frazil ) 
-                     zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) &
-                        - hicol(ji,jj) * zhicrit * ztwogp * zvrel2
-                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0*hicol(ji,jj) + zhicrit ) &
-                        - zhicrit * ztwogp * zvrel2
-                     zhicol_new = hicol(ji,jj) - zf/zfp
-                     hicol(ji,jj)   = zhicol_new
-
+                  DO WHILE ( iter < 20 ) 
+                     zf  = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj) * hicol(ji,jj) - zhicrit * zhicrit ) -   &
+                        &    hicol(ji,jj) * zhicrit * ztwogp * zvrel2
+                     zfp = ( hicol(ji,jj) - zhicrit ) * ( 3.0 * hicol(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2
+
+                     hicol(ji,jj) = hicol(ji,jj) - zf/zfp
                      iter = iter + 1
-
-                  END DO ! do while
+                  END DO
 
                ENDIF ! end of selection of pixels where ice forms
 
-            END DO ! loop on ji ends
-         END DO ! loop on jj ends
-      ! 
-      CALL lbc_lnk( zvrel(:,:), 'T', 1. )
-      CALL lbc_lnk( hicol(:,:), 'T', 1. )
+            END DO 
+         END DO 
+         ! 
+         CALL lbc_lnk( zvrel(:,:), 'T', 1. )
+         CALL lbc_lnk( hicol(:,:), 'T', 1. )
 
       ENDIF ! End of computation of frazil ice collection thickness
@@ -282,6 +270,5 @@
       ! Move from 2-D to 1-D vectors
       !------------------------------
-      ! If ocean gains heat do nothing 
-      ! 0therwise compute new ice formation
+      ! If ocean gains heat do nothing. Otherwise compute new ice formation
 
       IF ( nbpac > 0 ) THEN
@@ -297,13 +284,14 @@
          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) )
+         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) )
+         CALL tab_2d_1d( nbpac, rn_amax_1d(1:nbpac)     , rn_amax_2d, jpi, jpj, npac(1:nbpac) )
 
          !------------------------------------------------------------------------------!
@@ -316,11 +304,9 @@
          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) = rn_hnewice
-         END DO
-         IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac)
+         zh_newice(1:nbpac) = hicol_1d(1:nbpac)
 
          !----------------------
@@ -384,12 +370,17 @@
             ! salt flux
             sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice
-
+         END DO
+         
+         zv_frazb(:) = 0._wp
+         IF( ln_frazil ) THEN
             ! A fraction zfrazb of frazil ice is accreted at the ice bottom
-            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
-
+            DO ji = 1, nbpac
+               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
+         END IF
+         
          !-----------------
          ! Area of new ice
@@ -409,6 +400,6 @@
          ! we keep the excessive volume in memory and attribute it later to bottom accretion
          DO ji = 1, nbpac
-            IF ( za_newice(ji) >  ( rn_amax - zat_i_1d(ji) ) ) THEN
-               zda_res(ji)   = za_newice(ji) - ( rn_amax - zat_i_1d(ji) )
+            IF ( za_newice(ji) >  ( rn_amax_1d(ji) - zat_i_1d(ji) ) ) THEN
+               zda_res(ji)   = za_newice(ji) - ( rn_amax_1d(ji) - zat_i_1d(ji) )
                zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
                za_newice(ji) = za_newice(ji) - zda_res  (ji)
@@ -443,5 +434,5 @@
                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) +                                                      &
+               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 )
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 6455)
@@ -422,5 +422,5 @@
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  a_i(ji,jj,1)  = MIN( a_i(ji,jj,1), rn_amax )
+                  a_i(ji,jj,1)  = MIN( a_i(ji,jj,1), rn_amax_2d(ji,jj) )
                END DO
             END DO
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 6455)
@@ -80,7 +80,7 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN
-                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) )
-                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) )
+               IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN
+                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
+                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
                ENDIF
             END DO
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 6455)
@@ -94,7 +94,7 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF( at_i(ji,jj) > rn_amax .AND. a_i(ji,jj,jl) > 0._wp ) THEN
-                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) )
-                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax / at_i(ji,jj) ) )
+               IF( at_i(ji,jj) > rn_amax_2d(ji,jj) .AND. a_i(ji,jj,jl) > 0._wp ) THEN
+                  a_i (ji,jj,jl) = a_i (ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
+                  oa_i(ji,jj,jl) = oa_i(ji,jj,jl) * ( 1._wp - ( 1._wp - rn_amax_2d(ji,jj) / at_i(ji,jj) ) )
                ENDIF
             END DO
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 6455)
@@ -163,4 +163,20 @@
                rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes
                ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
+            END DO
+         END DO
+      END DO
+      ! Force the upper limit of ht_i to always be < hi_max (99 m).
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            rswitch = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) )
+            ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) )
+            a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch
+         END DO
+      END DO
+
+      DO jl = 1, jpl
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes
                ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
                o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
@@ -168,5 +184,5 @@
          END DO
       END DO
-
+      
       IF(  nn_icesal == 2  )THEN
          DO jl = 1, jpl
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 6455)
@@ -157,13 +157,5 @@
       ENDIF
 
-      IF ( iom_use( "icecolf" ) ) THEN 
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) ) )
-               z2d(ji,jj) = hicol(ji,jj) * rswitch
-            END DO
-         END DO
-         CALL iom_put( "icecolf"     , z2d              )        ! frazil ice collection thickness
-      ENDIF
+      IF ( iom_use( "icecolf" ) )   CALL iom_put( "icecolf", hicol )  ! frazil ice collection thickness
 
       CALL iom_put( "isst"        , sst_m               )        ! sea surface temperature
@@ -190,12 +182,13 @@
       CALL iom_put( "destrp"      , diag_trp_es         )        ! advected snw enthalpy (W/m2)
 
-      CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from brines
-      CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from brines
-      CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from brines
-      CALL iom_put( "sfxsni"      , sfx_sni * rday      )        ! salt flux from brines
-      CALL iom_put( "sfxopw"      , sfx_opw * rday      )        ! salt flux from brines
+      CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from bottom growth
+      CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melt
+      CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melt
+      CALL iom_put( "sfxsni"      , sfx_sni * rday      )        ! salt flux from snow ice formation
+      CALL iom_put( "sfxopw"      , sfx_opw * rday      )        ! salt flux from open water formation
       CALL iom_put( "sfxdyn"      , sfx_dyn * rday      )        ! salt flux from ridging rafting
-      CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from limupdate (resultant)
+      CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from residual
       CALL iom_put( "sfxbri"      , sfx_bri * rday      )        ! salt flux from brines
+      CALL iom_put( "sfxsub"      , sfx_sub * rday      )        ! salt flux from sublimation
       CALL iom_put( "sfx"         , sfx     * rday      )        ! total salt flux
 
@@ -235,4 +228,17 @@
       CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice 
       CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip 
+
+
+      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations  
+         DO jj = 1, jpj 
+            DO ji = 1, jpi
+               z2d(ji,jj)  = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) ! mean ice thickness
+            END DO
+         END DO
+         WHERE( z2d(:,:) < 0.2 .AND. z2d(:,:) > 0. ) ; z2da = wfx_bog
+         ELSEWHERE                                   ; z2da = 0._wp
+         END WHERE
+         CALL iom_put( "vfxthin", ( wfx_opw + z2da ) * ztmp )
+      ENDIF
       
       !--------------------------------
@@ -311,5 +317,5 @@
       !!
       !! History :
-      !!   4.1  !  2013-06  (C. Rousset)
+      !!   4.0  !  2013-06  (C. Rousset)
       !!----------------------------------------------------------------------
       INTEGER, INTENT( in ) ::   kt               ! ocean time-step index)
Index: branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 6452)
+++ branches/CNRS/dev_r6270_PISCES_QUOTA/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 6455)
@@ -45,4 +45,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qns_ice_1d  
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_bo_1d     
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rn_amax_1d
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sum_1d
@@ -83,4 +84,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_res_1d  
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_sub_1d
+
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d    !: <==> the 2D  sprecip
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   frld_1d       !: <==> the 2D  frld
@@ -91,4 +94,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   evap_ice_1d   !: <==> the 2D  evap_ice
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qprec_ice_1d  !: <==> the 2D  qprec_ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qevap_ice_1d  !: <==> the 3D  qevap_ice
    !                                                     ! to reintegrate longwave flux inside the ice thermodynamics
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   i0            !: fraction of radiation transmitted to the ice
@@ -107,4 +111,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_tot      !: Snow accretion/ablation        [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_surf     !: Ice surface accretion/ablation [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_sub      !: Ice surface sublimation [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bott     !: Ice bottom accretion/ablation  [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_snowice    !: Snow ice formation             [m of ice]
@@ -144,4 +149,5 @@
          &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,    & 
          &      hfx_dif_1d(jpij) , hfx_opw_1d(jpij) ,                      &
+         &      rn_amax_1d(jpij) ,                                         &
          &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) ,                      &
          &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) ,   &
@@ -153,7 +159,7 @@
          &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) ,  &
          &      dqns_ice_1d(jpij) , evap_ice_1d (jpij),                                         &
-         &      qprec_ice_1d(jpij), i0         (jpij) ,                                         &  
+         &      qprec_ice_1d(jpij), qevap_ice_1d(jpij), i0         (jpij) ,                     &  
          &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  &
-         &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) ,                     &
+         &      sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  &
          &      dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) ,                     &     
          &      dsm_i_si_1d(jpij) , hicol_1d    (jpij)                     , STAT=ierr(2) )
@@ -161,6 +167,6 @@
       ALLOCATE( t_su_1d   (jpij) , a_i_1d   (jpij) , ht_i_1d  (jpij) ,    &   
          &      ht_s_1d   (jpij) , fc_su    (jpij) , fc_bo_i  (jpij) ,    &    
-         &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) ,    &    
-         &      dh_snowice(jpij) , sm_i_1d  (jpij) , s_i_new  (jpij) ,    &
+         &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_sub (jpij) ,    &    
+         &      dh_i_bott (jpij) ,dh_snowice(jpij) , sm_i_1d  (jpij) , s_i_new  (jpij) ,    &
          &      t_s_1d(jpij,nlay_s) , t_i_1d(jpij,nlay_i) , s_i_1d(jpij,nlay_i) ,  &            
          &      q_i_1d(jpij,nlay_i+1) , q_s_1d(jpij,nlay_s) ,                        &
