Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90	(revision 8622)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice1D.F90	(revision 8623)
@@ -89,5 +89,5 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d    !: <==> the 2D  sprecip
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_1d        !: <==> the 2D  at_i
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_1d       !: <==> the 2D  at_i
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   ato_i_1d      !: <==> the 2D  ato_i
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhtur_1d      !: <==> the 2D  fhtur
@@ -103,7 +103,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_i_1d        !: <==> the 2D  a_i
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ib_1d       !: <==> the 2D  a_i_b
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_i_1d       !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ib_1d      !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_s_1d       !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_i_1d        !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ib_1d       !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_s_1d        !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_su         !: Surface Conduction flux 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fc_bo_i       !: Bottom  Conduction flux 
@@ -112,22 +112,25 @@
    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_s_mlt      !: Snow melt [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_snowice    !: Snow ice formation             [m of ice]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_1d       !: Ice bulk salinity [ppt]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_1d        !: Ice bulk salinity [ppt]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_new       !: Salinity of new ice at the bottom
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_i_1d        !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_s_1d        !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sv_i_1d       !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   oa_i_1d       !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_1d       !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_ip_1d       !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_i_1d        !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_s_1d        !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sv_i_1d      !:
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   oa_i_1d       !:
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d   !: corresponding to the 2D var  t_s
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_i_1d   !: corresponding to the 2D var  t_i
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sz_i_1d  !: profiled ice salinity
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e_i_1d   !:    Ice  enthalpy per unit volume
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e_s_1d   !:    Snow enthalpy per unit volume
-
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   eh_i_old !: ice heat content (q*h, J.m-2)
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_old  !: ice thickness layer (m)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ip_1d       !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_frac_1d  !:
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_i_1d      !: corresponding to the 2D var  t_i
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sz_i_1d     !: profiled ice salinity
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e_i_1d      !:    Ice  enthalpy per unit volume
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   e_s_1d      !:    Snow enthalpy per unit volume
+
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   eh_i_old    !: ice heat content (q*h, J.m-2)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_old     !: ice thickness layer (m)
 
    ! Conduction flux diagnostics (SIMIP)
@@ -200,6 +203,7 @@
          &      h_i_1d  (jpij) , h_ib_1d  (jpij) , h_s_1d (jpij) , fc_su  (jpij) , fc_bo_i(jpij) ,  &    
          &      dh_s_tot (jpij) , dh_i_surf (jpij) , dh_i_sub(jpij) ,                                  &    
-         &      dh_i_bott(jpij) , dh_snowice(jpij) , s_i_1d (jpij) , s_i_new(jpij) , &
+         &      dh_i_bott(jpij) , dh_s_mlt(jpij), dh_snowice(jpij) , s_i_1d (jpij) , s_i_new(jpij) , &
          &      a_ip_1d  (jpij) , v_ip_1d   (jpij) , v_i_1d  (jpij) , v_s_1d (jpij) , &
+         &      h_ip_1d  (jpij) , a_ip_frac_1d(jpij) , &
          &      sv_i_1d (jpij) , oa_i_1d   (jpij) , STAT=ierr(ii) )
       !
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 8622)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90	(revision 8623)
@@ -76,7 +76,8 @@
       IF( icount == 0 ) THEN
          !                          ! water flux
-         pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     &
-            &                    wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
-            &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)    &
+         pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  &
+            &                    wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  &
+            &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  &
+            &                    wfx_ice_sub(:,:) + wfx_spr(:,:)  &
             &                  ) * e1e2t(:,:) ) * zconv
          !
@@ -96,12 +97,13 @@
 
          pdiag_t = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     &
-            &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * zconv
+            &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv
 
       ELSEIF( icount == 1 ) THEN
 
          ! water flux
-         zfv  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     &
-            &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  &
-            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        &
+         zfv  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  &
+            &                wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  &
+            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  &
+            &                wfx_ice_sub(:,:) + wfx_spr(:,:)  &
             &              ) * e1e2t(:,:) ) * zconv - pdiag_fv
 
@@ -117,8 +119,8 @@
  
          ! outputs
-         zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * zconv  &
+         zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv  &
             &     - pdiag_v ) * r1_rdtice - zfv ) * rday
 
-         zs = ( ( glob_sum( SUM( sv_i * rhoic            , dim=3 ) * e1e2t ) * zconv  &
+         zs = ( ( glob_sum( SUM( sv_i * rhoic             , dim=3 ) * e1e2t ) * zconv  &
             &     - pdiag_s ) * r1_rdtice + zfs ) * rday
 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8622)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8623)
@@ -231,9 +231,10 @@
             s_i_new   (1:npti) = 0._wp ; dh_s_tot (1:npti) = 0._wp  ! --- some init --- !  (important to have them here) 
             dh_i_surf (1:npti) = 0._wp ; dh_i_bott(1:npti) = 0._wp
-            dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp
+            dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp
             !
             IF( ln_icedH ) THEN                                     ! --- growing/melting --- !
                               CALL ice_thd_zdf                             ! Ice/Snow Temperature profile
                               CALL ice_thd_dh                              ! Ice/Snow thickness   
+                              CALL ice_thd_pnd                             ! Melt ponds formation
                               CALL ice_thd_ent( e_i_1d(1:npti,:) )         ! Ice enthalpy remapping
             ENDIF
@@ -273,6 +274,4 @@
       !
       IF( ln_icedO )       CALL ice_thd_do                 ! --- frazil ice growing in leads --- !
-      !
-                           CALL ice_thd_pnd( kt )          ! --- Melt ponds --- !
       !
       ! controls
@@ -365,17 +364,20 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
          CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     )
-         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,kl)     )
-         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d(1:npti), h_s(:,:,kl)     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl)     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl)     )
          CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     )
-         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,kl)     )
+         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl)     )
          DO jk = 1, nlay_s
-            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)   )
-            CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)   )
+            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)    )
+            CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)    )
          END DO
          DO jk = 1, nlay_i
-            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl)   )
-            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl)   )
-            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)   )
-         END DO
+            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl)  )
+            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl)  )
+            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  )
+         END DO
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
          !
          CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d(1:npti), qprec_ice        )
@@ -406,4 +408,5 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr          )
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd          )
          !
          CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog          )
@@ -457,23 +460,27 @@
          !
          ! Change thickness to volume (replaces routine ice_var_eqv2glo)
-         v_i_1d(1:npti)  = h_i_1d(1:npti) * a_i_1d(1:npti)
-         v_s_1d(1:npti)  = h_s_1d(1:npti) * a_i_1d(1:npti)
-         sv_i_1d(1:npti) = s_i_1d(1:npti) * v_i_1d(1:npti)
+         v_i_1d (1:npti) = h_i_1d (1:npti) * a_i_1d (1:npti)
+         v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti)
+         sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti)
+         v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
          
          CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             )
          CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     )
-         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,kl)     )
-         CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d(1:npti), h_s(:,:,kl)     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl)     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl)     )
          CALL tab_1d_2d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     )
-         CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,kl)     )
+         CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl)     )
          DO jk = 1, nlay_s
-            CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) )
-            CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) )
+            CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)    )
+            CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)    )
          END DO
          DO jk = 1, nlay_i
-            CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl) )
-            CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl) )
-            CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) )
-         END DO
+            CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl)  )
+            CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl)  )
+            CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  )
+         END DO
+         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
          !
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
@@ -491,4 +498,5 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr        )
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam        )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        )
          !
          CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog        )
@@ -526,4 +534,5 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) )
          !
       END SELECT
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 8622)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90	(revision 8623)
@@ -148,6 +148,7 @@
             ! Contribution to mass flux
             wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * h_s_1d(ji) * a_i_1d(ji) * r1_rdtice
+            dh_s_mlt(ji)       = dh_s_mlt(ji) - h_s_1d(ji)
             ! updates
-            h_s_1d(ji)   = 0._wp
+            h_s_1d(ji)    = 0._wp
             e_s_1d (ji,1) = 0._wp
             t_s_1d (ji,1) = rt0
@@ -211,4 +212,5 @@
          ! snow melting only = water into the ocean (then without snow precip), >0
          wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
+         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
          ! updates available heat + precipitations after melting
          zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )      
@@ -233,7 +235,8 @@
             hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice 
             ! snow melting only = water into the ocean (then without snow precip)
-            wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
+            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice
+            dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,jk)
             ! updates available heat + thickness
-            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
+            zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
             h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
          END DO
@@ -571,5 +574,5 @@
          zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting
          dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1)
-         h_s_1d   (ji)  = h_s_1d(ji)   + zdeltah(ji,1)
+         h_s_1d   (ji)   = h_s_1d(ji)   + zdeltah(ji,1)
         
          zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                ! update available heat (J.m-2)
@@ -577,5 +580,6 @@
          hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)
          ! Contribution to mass flux
-         wfx_snw_sum_1d(ji)  =  wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
+         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice
+         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1)
          !    
          ! Remaining heat flux (W.m-2) is sent to the ocean heat budget
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90	(revision 8622)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90	(revision 8623)
@@ -6,22 +6,23 @@
    !! history :       ! Original code by Daniela Flocco and Adrian Turner
    !!            1.0  ! 2012    (O. Lecomte) Adaptation for scientific tests (NEMO3.1)
-   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6
+   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO4
    !!----------------------------------------------------------------------
 #if defined key_lim3
    !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                 LIM3 sea-ice model
+   !!   'key_lim3' :                                     ESIM sea-ice model
    !!----------------------------------------------------------------------
    !!   ice_thd_pnd_init      : some initialization and namelist read
    !!   ice_thd_pnd           : main calling routine
    !!----------------------------------------------------------------------
-   USE phycst           ! physical constants
-   USE dom_oce          ! ocean space and time domain
-   USE ice              ! LIM-3 variables
+   USE phycst         ! physical constants
+   USE dom_oce        ! ocean space and time domain
+   USE ice            ! sea-ice: variables
+   USE ice1D          ! sea-ice: thermodynamics variables
+   USE icetab         ! sea-ice: 1D <==> 2D transformation
    !
-   USE lbclnk           ! lateral boundary conditions - MPP exchanges
-   USE lib_mpp          ! MPP library
-   USE in_out_manager   ! I/O manager
-   USE lib_fortran      ! glob_sum
-   USE timing           ! Timing
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
+   USE timing         ! Timing
 
    IMPLICIT NONE
@@ -40,11 +41,11 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
-   !! $Id: limdyn.F90 6994 2016-10-05 13:07:10Z clem $
+   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
+   !! $Id: icethd_pnd.F90 8420 2017-10-05 13:07:10Z clem $
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
 
-   SUBROUTINE ice_thd_pnd( kt )
+   SUBROUTINE ice_thd_pnd
       !!-------------------------------------------------------------------
       !!               ***  ROUTINE ice_thd_pnd   ***
@@ -52,11 +53,8 @@
       !! ** Purpose :   change melt pond fraction
       !!                
-      !! ** Method  :   brutal force
+      !! ** Method  :   brut force
       !!
       !! ** Action  : - 
       !!              - 
-      !!-------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! number of iteration
-      INTEGER ::   ji, jj, jl     ! dummy loop indices
       !!-------------------------------------------------------------------
 
@@ -96,16 +94,19 @@
       !!    
       !!-------------------------------------------------------------------
-      
-      WHERE ( a_i(:,:,:) > 0._wp .AND. t_su(:,:,:) >= rt0 ) 
-         a_ip_frac = rn_apnd
-         h_ip      = rn_hpnd    
-         v_ip      = a_ip_frac * a_i * h_ip 
-         a_ip      = a_ip_frac * a_i
-      ELSEWHERE
-         a_ip      = 0._wp
-         h_ip      = 0._wp
-         v_ip      = 0._wp
-         a_ip_frac = 0._wp
-      END WHERE
+      INTEGER  ::   ji        ! loop indices
+      !!-------------------------------------------------------------------
+      DO ji = 1, npti
+         
+         IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
+            a_ip_frac_1d(ji) = rn_apnd
+            h_ip_1d(ji)      = rn_hpnd    
+            a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji)
+         ELSE
+            a_ip_frac_1d(ji) = 0._wp
+            h_ip_1d(ji)      = 0._wp    
+            a_ip_1d(ji)      = 0._wp
+         ENDIF
+         
+      END DO
       
    END SUBROUTINE pnd_CST
@@ -117,7 +118,11 @@
       !! ** Purpose    : Compute melt pond evolution
       !!
-      !! ** Method     : Empirical method. A fraction of meltwater is accumulated 
-      !!                 in pond volume. It is then released exponentially when
-      !!                 surface is freezing.
+      !! ** Method     : Empirical method. A fraction of meltwater is accumulated in ponds 
+      !!                 and sent to ocean when surface is freezing
+      !!
+      !!                 pond growth:      Vp = Vp + dVmelt
+      !!                    with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i
+      !!                 pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)
+      !!                    with Tp = -2degC
       !!
       !! ** Tunable parameters : (no real expertise yet, ideas?)
@@ -131,130 +136,91 @@
       !!    
       !!-------------------------------------------------------------------
-      
-      INTEGER, DIMENSION(jpij)         ::   indxi             ! compressed indices for cells with ice melting
-      INTEGER, DIMENSION(jpij)         ::   indxj             !
-
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zwfx_mlw          ! available meltwater for melt ponding
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zrfrac            ! fraction of available meltwater retained for melt ponding
-
       REAL(wp), PARAMETER ::   zrmin  = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
       REAL(wp), PARAMETER ::   zrmax  = 0.70_wp  ! maximum   ''           ''       ''        ''            ''
-      REAL(wp), PARAMETER ::   zrexp  = 0.01_wp  ! rate constant to refreeze melt ponds
       REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp ! pond aspect ratio
-
-      REAL(wp) ::   zhi               ! dummy ice thickness
-      REAL(wp) ::   zhs               ! dummy snow depth
-      REAL(wp) ::   zTp, z1_Tp        ! reference temperature
-      REAL(wp) ::   zdTs              ! dummy temperature difference
+      REAL(wp), PARAMETER ::   zTp = -2._wp ! 
+
+      REAL(wp) ::   zfr_mlt            ! fraction of available meltwater retained for melt ponding
+      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding
+      REAL(wp) ::   z1_Tp        ! reference temperature
       REAL(wp) ::   z1_rhofw          ! inverse freshwater density
       REAL(wp) ::   z1_zpnd_aspect    ! inverse pond aspect ratio
-      REAL(wp) ::   zvpold            ! dummy pond volume
-
-      INTEGER  ::   ji, jj, jl, ij    ! loop indices
-      INTEGER  ::   icells            ! size of dummy array
+      REAL(wp) ::   zvpnd            ! dummy pond volume
+      REAL(wp) ::   zfac, zdum
+
+      INTEGER  ::   ji        ! loop indices
       !!-------------------------------------------------------------------
       z1_rhofw       = 1. / rhofw 
       z1_zpnd_aspect = 1. / zpnd_aspect
-      zTp            = -2. 
       z1_Tp          = 1._wp / zTp 
 
-      a_ip_frac(:,:,:) = 0._wp
-      h_ip     (:,:,:) = 0._wp
-
-      !------------------------------------------------------------------
-      ! Available melt water for melt ponding and corresponding fraction
-      !------------------------------------------------------------------
-
-      zwfx_mlw(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:), 0._wp )        ! available meltwater for melt ponding
-
-      ! NB: zwfx_mlw can be slightly negative for very small values (why?)
-      ! This can in some occasions give negative
-      ! v_ip in the first category, which then gives crazy pond
-      ! fractions and crashes the code as soon as the melt-pond
-      ! radiative coupling is activated
-      ! if we understand and remove why wfx_sum or wfx_snw could be
-      ! negative, then, we can remove the MAX
-      ! NB: I now changed to wfx_snw_sum, this may fix the problem. 
-      ! We should check
-
-      zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:)  
-
-      DO jl = 1, jpl   
-
-         !------------------------------------------------------------------------------
-         ! Identify grid cells where ponds should be updated (can probably be improved)
-         !------------------------------------------------------------------------------
-         indxi(:) = 0
-         indxj(:) = 0
-         icells   = 0
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF ( a_i(ji,jj,jl) > epsi10 ) 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)
-
-            zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
-            zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)
-
-            IF ( zhi < rn_himin) THEN   !--- Remove ponds on thin ice if ice is too thin
-
-               a_ip(ji,jj,jl)      = 0._wp                               !--- Dump ponds
-               v_ip(ji,jj,jl)      = 0._wp
-               a_ip_frac(ji,jj,jl) = 0._wp
-               h_ip(ji,jj,jl)      = 0._wp
-
-               IF( ln_pnd_fwb )   wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + v_ip(ji,jj,jl)  !--- Give freshwater to the ocean
-
-            ELSE                        !--- Update pond characteristics
-
-               !--- Add retained melt water to melt ponds
-               ! v_ip should never be positive, otherwise code crashes
-               ! MV: as far as I saw, UM5 can create very small negative v_ip values
-               ! hence I added the max, which was not required with Prather (1 yr run)
-               v_ip(ji,jj,jl) = MAX( v_ip(ji,jj,jl), 0._wp ) + zrfrac(ji,jj,jl) * z1_rhofw * zwfx_mlw(ji,jj) * a_i(ji,jj,jl) * rdt_ice
-
-               !--- Shrink pond due to refreezing
-               zdTs           = MAX ( zTp - t_su(ji,jj,jl) + rt0 , 0. )
-
-               zvpold         = v_ip(ji,jj,jl)
-
-               v_ip(ji,jj,jl) = v_ip(ji,jj,jl) * EXP( zrexp * zdTs * z1_Tp )
-
-               !--- Dump meltwater due to refreezing ( of course this is wrong
-               !--- but this parameterization is too simple )
-               IF ( ln_pnd_fwb ) THEN
-                  wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice
-               ENDIF
-               
-               a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) )
-               !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative
-
-               h_ip(ji,jj,jl)      = zpnd_aspect * a_ip_frac(ji,jj,jl)
-
-               a_ip(ji,jj,jl)      = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl)
-
+      DO ji = 1, npti
+         !                                                        !--------------------------------!
+         IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  !
+            !                                                     !--------------------------------!
+            !--- Remove ponds on thin ice and send water into the ocean
+            IF( ln_pnd_fwb ) THEN
+               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + h_ip_1d(ji) * a_ip_1d(ji) * rhofw * r1_rdtice
             ENDIF
-
-         END DO
-
+            a_ip_1d(ji)      = 0._wp
+            a_ip_frac_1d(ji) = 0._wp
+            h_ip_1d(ji)      = 0._wp
+            !                                                     !--------------------------------!
+         ELSE                                                     ! Case ice thickness >= rn_himin !
+            !                                                     !--------------------------------!
+            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
+            zvpnd       = v_ip_1d(ji)
+
+            ! available meltwater for melt ponding [m, >0] and fraction
+            zdv_mlt = -( dh_i_surf(ji)*rhoic + dh_s_mlt(ji)*rhosn ) * z1_rhofw * a_i_1d(ji)
+            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  !!clem CICE doc
+!!!            zfr_mlt = zrmin + zrmax * a_i_1d(ji)  !!clem Holland paper 
+
+            !--- Pond gowth
+            ! v_ip should never be negative, otherwise code crashes
+            ! MV: as far as I saw, UM5 can create very small negative v_ip values (not Prather)
+            v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt )
+
+            IF( ln_pnd_fwb .AND. zdv_mlt > 0._wp ) THEN
+               ! melt ponds mass flux (<0)
+               zfac = zfr_mlt * zdv_mlt * rhofw * r1_rdtice
+               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
+
+               ! adjust ice/snow melting (>0) to take into account ponding
+               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )
+               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
+               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum)
+            ENDIF
+
+            !--- Pond contraction (due to refreezing)
+            v_ip_1d(ji)      = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp )
+
+            ! --- Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
+            !        h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i
+            a_ip_1d(ji)      = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) )
+            !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative
+            a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
+            h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
+
+            !! clem: there is no clever way to conserve mass here...
+            
+!            IF( ln_pnd_fwb ) THEN
+!               ! melt ponds mass flux (>0 when ponds shrink)
+!               zfac = ( v_ip_1d(ji) - zvpnd ) * rhofw * r1_rdtice
+!               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac!!
+!
+!               ! adjust ice/snow
+!               zfac = ( v_ip_1d(ji) - zvpnd ) * rhofw / ( a_i_1d(ji)*h_s_1d(ji)*rhosn + a_i_1d(ji)*h_i_1d(ji)*rhoic )!
+!
+!               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + h_s_1d(ji) * zfac * a_i_1d(ji) * rhosn * r1_rdtice
+!               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     + h_i_1d(ji) * zfac * a_i_1d(ji) * rhoic * r1_rdtice
+!               
+!               !!h_s_1d(ji) = h_s_1d(ji) * ( 1._wp + zfac )
+!               !!h_i_1d(ji) = h_i_1d(ji) * ( 1._wp + zfac )!
+!
+!            ENDIF
+
+         ENDIF
       END DO
-
-      !--- Remove retained meltwater from surface fluxes 
-
-      IF ( ln_pnd_fwb ) THEN
-
-         wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) ) 
-         wfx_sum(:,:)     = wfx_sum(:,:)     *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) )
-
-      ENDIF
 
    END SUBROUTINE pnd_H12
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8622)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8623)
@@ -164,7 +164,5 @@
             ! mass flux from ice/ocean
             wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   &
-               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 
-
-            IF ( ln_pnd_fwb )   wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj)
+               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj)
 
             ! add the snow melt water to snow mass flux to the ocean
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8622)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8623)
@@ -255,7 +255,8 @@
       !!-------------------------------------------------------------------
       !
-      v_i (:,:,:) = h_i(:,:,:) * a_i(:,:,:)
-      v_s (:,:,:) = h_s(:,:,:) * a_i(:,:,:)
-      sv_i(:,:,:) = s_i(:,:,:) * v_i(:,:,:)
+      v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:)
+      v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:)
+      sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:)
+      v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
       !
    END SUBROUTINE ice_var_eqv2glo
@@ -512,4 +513,5 @@
       ! open water = 1 if at_i=0
       WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp
+
       !
    END SUBROUTINE ice_var_zapsmall
