Index: branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/file_def_nemo-lim.xml
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/file_def_nemo-lim.xml	(revision 8900)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/file_def_nemo-lim.xml	(revision 8906)
@@ -14,8 +14,8 @@
    <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE.">  <!-- 5d files -->   	
      <file id="file21" name_suffix="_icemod" description="ice variables" enabled=".true." >
-       <field field_ref="snothic"         name="snthic" />
+       <field field_ref="snwthic"         name="snthic" />
        <field field_ref="icethic"         name="sithic" />
        <field field_ref="icevolu"         name="sivolu" />
-       <field field_ref="snowvol"         name="snvolu" />
+       <field field_ref="snwvolu"         name="snvolu" />
        <field field_ref="iceconc"         name="siconc" />
      </file>
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90	(revision 8900)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90	(revision 8906)
@@ -669,8 +669,10 @@
                ! Put the melt pond water into the ocean
                !------------------------------------------            
-               IF ( ln_pnd_fwb ) THEN
-                  wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean
-                     &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice
-               ENDIF
+               ! clem: I think the following lines must be commented since there
+               !       is no net mass flux between melt ponds and the ocean (see icethd_pnd.F90 for ex.)
+               !IF ( ln_pnd_fwb ) THEN
+               !   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean
+               !      &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice
+               !ENDIF
 
                ! virtual salt flux to keep salinity constant
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8900)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8906)
@@ -138,14 +138,14 @@
       !
       CASE( jp_usr )                !--- user defined formulation
-                                CALL usrdef_sbc_ice_flx( kt )
+                                  CALL usrdef_sbc_ice_flx( kt )
          !
       CASE( jp_blk )                !--- bulk formulation
-                                CALL blk_ice_flx( t_su, h_s, h_i, alb_ice )    ! 
-         IF( ln_mixcpl       )  CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
-         IF( nn_flxdist /= 2 )  CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
+                                  CALL blk_ice_flx( t_su, h_s, h_i, alb_ice )    ! 
+         IF( ln_mixcpl        )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
+         IF( nn_flxdist /= -1 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
          !
       CASE ( jp_purecpl )           !--- coupled formulation
-                                CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
-         IF( nn_flxdist == 2 )  CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
+                                  CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )
+         IF( nn_flxdist /= -1 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist )
          !
       END SELECT
@@ -294,13 +294,12 @@
       CASE( -1  )
          IF(lwp) WRITE(numout,*) '   ESIM: use per-category fluxes (nn_flxdist = -1) '
-         IF( ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_flxdist for ESIM in coupled mode must be 0 or 2' )
       CASE(  0  )
          IF(lwp) WRITE(numout,*) '   ESIM: use average per-category fluxes (nn_flxdist = 0) '
       CASE(  1  )
          IF(lwp) WRITE(numout,*) '   ESIM: use average then redistribute per-category fluxes (nn_flxdist = 1) '
-         IF( ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_flxdist for ESIM in coupled mode must be 0 or 2' )
+         IF( ln_cpl )         CALL ctl_stop( 'ice_thd_init: the chosen nn_flxdist for ESIM in coupled mode must be /=1' )
       CASE(  2  )
          IF(lwp) WRITE(numout,*) '   ESIM: Redistribute a single flux over categories (nn_flxdist = 2) '
-         IF( .NOT. ln_cpl )   CALL ctl_stop( 'ice_thd_init : the chosen nn_flxdist for ESIM in forced mode cannot be 2' )
+         IF( .NOT. ln_cpl )   CALL ctl_stop( 'ice_thd_init: the chosen nn_flxdist for ESIM in forced mode must be /=2' )
       CASE DEFAULT
          CALL ctl_stop( 'ice_thd_init: ESIM option, nn_flxdist, should be between -1 and 2' )
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8900)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8906)
@@ -296,19 +296,19 @@
          ! 7) Make sure h_i >= minimum ice thickness hi_min
          !----------------------------------------------------------------------------------------------
-         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,1)   )
-         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,1)    )
-         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,1)   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,1) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,1) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,1) )
          !
          DO ji = 1, npti
             IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN
                a_i_1d (ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 
-               IF ( ln_pnd_H12 )    a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
+               IF( ln_pnd_H12 )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
                h_i_1d(ji) = rn_himin
             ENDIF
          END DO
          !
-         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i(:,:,1)   )
-         CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d  (1:npti), a_i(:,:,1)    )
-         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,1)   )
+         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,1) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,1) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,1) )
          !
       ENDIF
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90	(revision 8900)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90	(revision 8906)
@@ -125,5 +125,5 @@
       !!                 pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp)
       !!                    with Tp = -2degC
-      !!
+      !!  
       !! ** Tunable parameters : (no real expertise yet, ideas?)
       !! 
@@ -136,18 +136,17 @@
       !!    
       !!-------------------------------------------------------------------
-      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 ::   zpnd_aspect = 0.8_wp ! pond aspect ratio
-      REAL(wp), PARAMETER ::   zTp = -2._wp ! 
-
-      REAL(wp) ::   zfr_mlt            ! 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 ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio
+      REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature
+
+      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) ::   zvpnd            ! dummy pond volume
+      REAL(wp) ::   z1_Tp            ! inverse reference temperature
+      REAL(wp) ::   z1_rhofw         ! inverse freshwater density
+      REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio
       REAL(wp) ::   zfac, zdum
 
-      INTEGER  ::   ji        ! loop indices
+      INTEGER  ::   ji   ! loop indices
       !!-------------------------------------------------------------------
       z1_rhofw       = 1. / rhofw 
@@ -159,8 +158,5 @@
          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
+            !--- Remove ponds on thin ice
             a_ip_1d(ji)      = 0._wp
             a_ip_frac_1d(ji) = 0._wp
@@ -170,22 +166,21 @@
             !                                                     !--------------------------------!
             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
+            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc
+            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from 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 )
 
+            ! melt pond mass flux (<0)
             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
+               ! adjust ice/snow melting flux to balance melt pond flux (>0)
                zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )
                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
@@ -193,31 +188,12 @@
             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
+            !--- 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
Index: branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8900)
+++ branches/2017/dev_CNRS_2017/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8906)
@@ -479,7 +479,4 @@
                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice
                hfx_res(ji,jj)  = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0
-               IF( ln_pnd_fwb ) THEN 
-                  wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
-               ENDIF
                !-----------------------------------------------------------------
                ! Zap snow energy 
