Index: /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8178)
+++ /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8179)
@@ -323,4 +323,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw     !: snow-ocean mass exchange   [kg.m-2.s-1]
    ! MV MP 2016
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sum !: surface melt component of wfx_snw [kg.m-2.s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_pnd     !: melt pond-ocean mass exchange   [kg.m-2.s-1]
    ! END MV MP 2016
@@ -508,5 +509,5 @@
          &      wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) , wfx_lam(jpi,jpj) ,     &
          ! MV MP 2016
-         &      wfx_pnd(jpi,jpj) ,                                                              &
+         &      wfx_pnd(jpi,jpj) , wfx_snw_sum(jpi,jpj) ,                                       &
          ! END MV MP 2016
          &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     &
Index: /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90
===================================================================
--- /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8178)
+++ /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8179)
@@ -129,4 +129,12 @@
       ENDIF
 
+      IF ( ln_pnd .AND. ( nn_pnd_scheme == 0 ) .AND. ( ln_pnd_fw ) ) THEN
+         WRITE(numout,*) ' Prescribed melt ponds do not conserve fresh water mass, hence ln_pnd_fw must be set to false ' 
+         ln_pnd_fw     = .FALSE.
+         IF(lwp) THEN                     ! control print
+            WRITE(numout,*) '    Active melt ponds freshwater coupling                       ln_pnd_fw     = ', ln_pnd_fw
+         ENDIF
+      ENDIF
+
       IF ( ln_pnd .AND. ( nn_pnd_scheme == 2 ) .AND. ( jpl == 1 ) ) THEN
          WRITE(numout,*) ' Topographic melt ponds are incompatible with jpl = 1 '
@@ -280,5 +288,5 @@
         !------------------------------------------------------------------
  
-        zwfx_mlw(:,:) = MAX( wfx_sum(:,:) + wfx_snw(:,:), 0._wp )        ! available meltwater for melt ponding
+        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?)
@@ -289,4 +297,6 @@
                 ! 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(:,:,:)  
@@ -352,4 +362,5 @@
  
                 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)
@@ -374,8 +385,8 @@
         IF ( ln_pnd_fw ) THEN
  
-            wfx_snw(:,:) = wfx_snw(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
- 
-            wfx_sum(:,:) = wfx_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
- 
+            wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) ) 
+
+            wfx_sum(:,:)     = wfx_sum(:,:) * ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
+
         ENDIF
 
@@ -496,5 +507,5 @@
        kcells(:)   = 0
 
-       zvolp(:,:) = wfx_sum(:,:) + wfx_snw(:,:) + vt_ip(:,:) ! Total available melt water, to be distributed as melt ponds
+       zvolp(:,:) = wfx_sum(:,:) + wfx_snw_sum(:,:) + vt_ip(:,:) ! Total available melt water, to be distributed as melt ponds
        zTsfcn(:,:,:) = zTsfcn(:,:,:) - rt0                   ! Convert in Celsius
  
Index: /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 8178)
+++ /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 8179)
@@ -177,4 +177,7 @@
             IF ( ln_pnd_fw ) &
                wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj)
+
+            ! add the snow melt water to snow mass flux to the ocean
+            wfx_snw(:,:)      = wfx_snw(:,:) + wfx_snw_sum(:,:) 
 
             ! mass flux at the ocean/ice interface
Index: /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8178)
+++ /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8179)
@@ -463,4 +463,5 @@
          
          CALL tab_2d_1d( nbpb, wfx_snw_1d (1:nbpb), wfx_snw         , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, wfx_snw_sum_1d(1:nbpb), wfx_snw_sum  , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, wfx_sub_1d (1:nbpb), wfx_sub         , jpi, jpj, npb(1:nbpb) )
          
@@ -514,4 +515,5 @@
          
          CALL tab_1d_2d( nbpb, wfx_snw       , npb, wfx_snw_1d(1:nbpb)   , jpi, jpj )
+         CALL tab_1d_2d( nbpb, wfx_snw_sum   , npb, wfx_snw_sum_1d(1:nbpb),jpi, jpj )
          CALL tab_1d_2d( nbpb, wfx_sub       , npb, wfx_sub_1d(1:nbpb)   , jpi, jpj )
          
Index: /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8178)
+++ /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8179)
@@ -169,5 +169,5 @@
             hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_1d(ji,1) * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
             ! Contribution to mass flux
-            wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
+            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * ht_s_1d(ji) * a_i_1d(ji) * r1_rdtice
             ! updates
             ht_s_1d(ji)   = 0._wp
@@ -233,5 +233,5 @@
          hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_rdtice
          ! snow melting only = water into the ocean (then without snow precip), >0
-         wfx_snw_1d(ji) = wfx_snw_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
          ! updates available heat + precipitations after melting
          zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )      
@@ -255,5 +255,5 @@
             hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * q_s_1d(ji,jk) * r1_rdtice 
             ! snow melting only = water into the ocean (then without snow precip)
-            wfx_snw_1d(ji)   = wfx_snw_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
             ! updates available heat + thickness
             zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )
@@ -607,5 +607,5 @@
          hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * q_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0)
          ! Contribution to mass flux
-         wfx_snw_1d(ji)  =  wfx_snw_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
          !    
          ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
Index: /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8178)
+++ /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8179)
@@ -54,4 +54,5 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_1d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_snw_sum_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sub_1d
 
@@ -142,5 +143,5 @@
       ii = ii + 1
       ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_1d    (jpij) ,                     &
-         &      fhtur_1d   (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) ,                     &
+         &      fhtur_1d   (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) , wfx_snw_sum_1d(jpij) , &
          &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d (jpij) , wfx_bom_1d(jpij) ,  &
          &      wfx_sum_1d(jpij)  , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) , wfx_res_1d(jpij) ,  &
Index: /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90
===================================================================
--- /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 8178)
+++ /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 8179)
@@ -167,9 +167,4 @@
          ENDIF 
 
-         ! Note that this way of coding things means that the thickness
-         ! dependence of bare ice albedo is also altering the ponded ice albedo
-         ! I'm not sure that this is valid, but at least, this maintains
-         ! consistency with previous formulation (see below)
-     
          DO jl = 1, ijpl
             DO jj = 1, jpj
Index: /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8178)
+++ /branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8179)
@@ -666,5 +666,5 @@
 
       ! MV MP 2016
-      wfx_pnd(:,:) = 0._wp
+      wfx_pnd(:,:) = 0._wp   ;   wfx_snw_sum(:,:) = 0._wp
       ! END MV MP 2016
       
