Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8141)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8142)
@@ -46,5 +46,5 @@
 !------------------------------------------------------------------------------
                   ! -- limistate -- !
-   ln_limini      = .false.         !  activate ice initialization (T) or not (F)
+   ln_limini      = .true.          !  activate ice initialization (T) or not (F)
    ln_limini_file = .false.         !  netcdf file provided for initialization (T) or not (F)
    rn_thres_sst   =  0.5            !  maximum water temperature with initial ice (degC)
@@ -169,8 +169,10 @@
 !------------------------------------------------------------------------------
                  ! -- limmp      -- !
-   ln_pnd         = .true.          !  activate melt ponds
-   nn_pnd_scheme  = 1               !  type of melt pond scheme:   =0 prescribed ( Tsu>0 ), =1 empirical, =2 topographic
-   nn_pnd_cpl     = 3               !  type of melt pond coupling: =0 passive, =1 full, =2 radiation only, =3 freshwater only
-   rn_apnd        = 0.0             !  prescribed pond fraction (0<rn_apnd<1), only if nn_pnd_scheme = 0
+   ln_pnd         = .false.         !  active melt ponds
+   ln_pnd_rad     = .false.         !  active melt ponds radiative coupling
+   ln_pnd_fw      = .false.         !  active melt ponds freshwater coupling
+   nn_pnd_scheme  = 0               !  type of melt pond scheme  : =0 prescribed ( Tsu=0 ), =1 empirical, =2 topographic
+   rn_apnd        = 0.2             !  prescribed pond fraction, at Tsu=0  : (0<rn_apnd<1, nn_pnd_scheme = 0)
+   rn_hpnd        = 0.05            !  prescribed pond depth, at Tsu=0     : (0<rn_apnd<1, nn_pnd_scheme = 0)
 /
 !------------------------------------------------------------------------------
@@ -188,7 +190,9 @@
    rn_por_rdg     =   0.3           !     (ln_ridging = T) porosity of newly ridged ice (Lepparanta et al., 1995)
    rn_fsnowrdg    =   0.5           !     (ln_ridging = T) snow volume fraction that survives in ridging
+   rn_fpondrdg    =   1.0           !     (ln_ridging = T) pond fraction that survives in ridging (small a priori)
    ln_rafting     =   .true.        !  rafting activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
    rn_hraft       =   0.75          !     (ln_rafting = T) threshold thickness for rafting (m)
    rn_craft       =   5.0           !     (ln_rafting = T) squeezing coefficient used in the rafting function
    rn_fsnowrft    =   0.5           !     (ln_rafting = T) snow volume fraction that survives in rafting
+   rn_fpondrft    =   1.0           !     (ln_rafting = T) pond fraction that survives in rafting (0.5 a priori)
 /
Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8141)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8142)
@@ -25,5 +25,5 @@
 !-----------------------------------------------------------------------
    nn_no       =       0   !  job number (no more used...)
-   cn_exp      =  "ORCA2"  !  experience name
+   cn_exp      =  "ORCA2"   !  experience name
    nn_it000    =       1   !  first time step
    nn_itend    =    5475   !  last  time step (std 5475)
@@ -506,11 +506,14 @@
 &namsbc_alb    !   albedo parameters
 !-----------------------------------------------------------------------
-   nn_ice_alb  =    0   !  parameterization of ice/snow albedo
-                        !     0: Shine & Henderson-Sellers (JGR 1985)
-                        !     1: "home made" based on Brandt et al. (J. Climate 2005)
-                        !                         and Grenfell & Perovich (JGR 2004)
-   rn_albice   =  0.53  !  albedo of bare puddled ice (values from 0.49 to 0.58)
-                        !     0.53 (default) => if nn_ice_alb=0
-                        !     0.50 (default) => if nn_ice_alb=1
+   nn_ice_alb   =    0   !  parameterization of ice/snow albedo
+                         !     0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo
+                         !     1: "home made" based on Brandt et al. (JClim 2005) and Grenfell & Perovich (JGR 2004),
+                         !        giving cloud-sky albedo
+                         !     2: "homme made", "1" rewritten, more naturally based on specific areas of ponds, snow & bare ice
+   rn_alb_sdry  =  0.85  !  dry snow albedo         : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1/2); obs 0.85-0.87 (cloud-sky)
+   rn_alb_smlt  =  0.75  !  melting snow albedo     : 0.65 ( '' )          ; 0.75 ( '' )            ; obs 0.72-0.82 ( '' )
+   rn_alb_idry  =  0.60  !  dry ice albedo          : 0.72 ( '' )          ; 0.60 ( '' )            ; obs 0.54-0.65 ( '' )
+   rn_alb_imlt  =  0.50  !  bare puddled ice albedo : 0.53 ( '' )          ; 0.50 ( '' )            ; obs 0.49-0.58 ( '' ) -> used only if no ponds
+   rn_alb_dpnd  =  0.27  !  bare puddled ice albedo : 0.25 ( '' )          ; 0.27 ( '' )            ; obs 0.10-0.30 ( '' )
 /
 !-----------------------------------------------------------------------
Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 8141)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 8142)
@@ -509,4 +509,6 @@
       ! MV MP 2016
       REAL(wp), POINTER, DIMENSION(:) ::   vprdg            ! pond volume of ridging ice
+      REAL(wp), POINTER, DIMENSION(:) ::   aprdg1           ! pond area of ridging ice
+      REAL(wp), POINTER, DIMENSION(:) ::   aprdg2           ! pond area of ridging ice
       ! END MV MP 2016
       REAL(wp), POINTER, DIMENSION(:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
@@ -525,4 +527,6 @@
       ! MV MP 2016
       REAL(wp), POINTER, DIMENSION(:) ::   vprft            ! pond volume of rafting ice
+      REAL(wp), POINTER, DIMENSION(:) ::   aprft1           ! pond area of rafted ice
+      REAL(wp), POINTER, DIMENSION(:) ::   aprft2           ! pond area of new rafted ice
       ! END MV MP 2016
       REAL(wp), POINTER, DIMENSION(:) ::   esrft , smrft    ! snow energy & salinity of rafting ice
@@ -539,10 +543,11 @@
       ! MV MP 2016
       !CALL wrk_alloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
-      CALL wrk_alloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, vprdg, dhr, dhr2 )
+      CALL wrk_alloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, vprdg, aprdg1, aprdg2, dhr, dhr2 )
       ! END MV MP 2016
       CALL wrk_alloc( jpij,        vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 )
       ! MV MP 2016
       !CALL wrk_alloc( jpij,        afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
-      CALL wrk_alloc( jpij,        afrft, arft1, arft2, virft, vsrft, esrft, vprft, smrft, oirft1, oirft2 )
+      CALL wrk_alloc(  jpij,        afrft, arft1, arft2, virft, vsrft, esrft, aprft1, aprft2)
+      CALL wrk_alloc ( jpij,        vprft, smrft, oirft1, oirft2 )
       ! END MV MP 2016
       CALL wrk_alloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw )
@@ -598,5 +603,5 @@
 
             !--------------------------------------------------------------------------
-            ! 3.4) Subtract area, volume, and energy from ridging 
+            ! 3.4) Substract area, volume, and energy from ridging 
             !     / rafting category n1.
             !--------------------------------------------------------------------------
@@ -609,5 +614,7 @@
             !MV MP 2016
             IF ( nn_pnd_scheme > 0 ) THEN
-               vprdg(ij) = v_ip(ji,jj, jl1) * afrac(ij)
+               aprdg1(ij) = a_ip(ji,jj, jl1) * afrac(ij)
+               aprdg2(ij) = a_ip(ji,jj, jl1) * afrac(ij) * krdg(ji,jj,jl1)
+               vprdg(ij)  = v_ip(ji,jj, jl1) * afrac(ij)
             ENDIF
             ! END MV MP 2016
@@ -621,5 +628,7 @@
             !MV MP 2016
             IF ( nn_pnd_scheme > 0 ) THEN
-               vprft(ij) = v_ip(ji,jj,jl1) * afrft(ij)
+               aprft1(ij) = a_ip (ji,jj,  jl1) * afrft(ij)
+               aprft2(ij) = a_ip (ji,jj,  jl1) * afrft(ij) * kraft
+               vprft(ij)  = v_ip(ji,jj,jl1)    * afrft(ij)
             ENDIF
             ! END MV MP 2016
@@ -690,4 +699,5 @@
             IF ( nn_pnd_scheme > 0 ) THEN
                v_ip (ji,jj,jl1) = v_ip (ji,jj,jl1) - vprdg (ij) - vprft (ij)
+               a_ip (ji,jj,jl1) = a_ip (ji,jj,jl1) - aprdg1(ij) - aprft1(ij)
             ENDIF
             ! END MV MP 2016
@@ -758,6 +768,8 @@
                ! MV MP 2016
                IF ( nn_pnd_scheme > 0 ) THEN
-                  v_ip (ji,jj,jl2) = v_ip (ji,jj,jl2)  + ( vprdg (ij) * rn_fpondrdg * fvol(ij)  +  &
-                  &                                        vprft (ij) * rn_fpondrft * zswitch(ij) )
+                  v_ip (ji,jj,jl2) = v_ip (ji,jj,jl2)  + ( vprdg (ij) * rn_fpondrdg * fvol(ij)  +   &
+                                                       &   vprft (ij) * rn_fpondrft * zswitch(ij) )
+                  a_ip (ji,jj,jl2) = a_ip(ji,jj,jl2)   + ( aprdg2(ij) * rn_fpondrdg * farea +       & 
+                                                       &   aprft2(ij) * rn_fpondrft * zswitch(ji) )
                ENDIF
                ! END MV MP 2016
@@ -782,10 +794,11 @@
       ! MV MP 2016
       !CALL wrk_dealloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
-      CALL wrk_dealloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, vprdg, dhr, dhr2 )
+      CALL wrk_dealloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, vprdg, aprdg1, aprdg2,  dhr, dhr2 )
       ! END MV MP 2016
       CALL wrk_dealloc( jpij,        vrdg1, vrdg2, vsw  , srdg1, srdg2, smsw, oirdg1, oirdg2 )
       ! MV MP 2016
       !CALL wrk_dealloc( jpij,        afrft, arft1, arft2, virft, vsrft, esrft, smrft, oirft1, oirft2 )
-      CALL wrk_dealloc( jpij,        afrft, arft1, arft2, virft, vsrft, esrft, vprft, smrft, oirft1, oirft2 )
+      CALL wrk_dealloc( jpij,        afrft, arft1, arft2, virft, vsrft, esrft, aprft1, aprft2, vprft )
+      CALL wrk_dealloc( jpij,        smrft, oirft1, oirft2 )
       CALL wrk_dealloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw )
       !
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 8141)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8142)
@@ -6,5 +6,7 @@
    !! history :       ! Original code by Daniela Flocco and Adrian Turner
    !!            1.0  ! 2012    (O. Lecomte) Adaptation for scientific tests (NEMO3.1)
-   !!            2.0  ! 2016    (O. Lecomte, C. Rousset, M. Vancoppenolle) Implementation in NEMO3.6
+   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6
+   !!                 ! NB: Only lim_mp_cstt and lim_mp_cesm work in this
+   !!                       version
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -137,5 +139,4 @@
    
 
-
    SUBROUTINE lim_mp( kt )
       !!-------------------------------------------------------------------
@@ -153,9 +154,5 @@
       INTEGER  ::   ji, jj, jl     ! dummy loop indices
 
-!     REAL(wp), POINTER, DIMENSION(:,:) ::   zfsurf  ! surface heat flux(obsolete, should be droped)
-      !
       !!-------------------------------------------------------------------
-
-      IF( nn_timing == 1 )  CALL timing_start('limthd')
 
       IF( nn_timing == 1 )  CALL timing_start('lim_mp')
@@ -180,9 +177,5 @@
       END SELECT
 
-      ! we should probably not aggregate here since we do it in lim_var_agg
-      ! before output, unless we need the total volume and faction else where
-
-      ! we should also make sure a_ip and v_ip are properly updated at the end
-      ! of the routine
+      IF( nn_timing == 1 )  CALL timing_stop('lim_mp')
 
    END SUBROUTINE lim_mp 
@@ -237,5 +230,5 @@
        !!                 surface is freezing.
        !!
-       !! ** Tunable parameters : (no expertise yet)
+       !! ** Tunable parameters : (no real expertise yet, ideas?)
        !! 
        !! ** Note       : Stolen from CICE for quick test of the melt pond
@@ -272,108 +265,122 @@
        !!-------------------------------------------------------------------
 
-       CALL wrk_alloc( jpi*jpj, indxi, indxj)
-       CALL wrk_alloc( jpi,jpj,     zwfx_mlw )
-       CALL wrk_alloc( jpi,jpj,jpl, zrfrac   )
-
-       z1_rhofw       = 1. / rhofw 
-       z1_zpnd_aspect = 1. / zpnd_aspect
-       zTp            = -2. 
-
-       !------------------------------------------------------------------
-       ! Available melt water for melt ponding and corresponding fraction
-       !------------------------------------------------------------------
-
-       zwfx_mlw(:,:) = wfx_sum(:,:) + wfx_snw(:,:)        ! available meltwater for melt ponding
-
-       zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:)  
-
-       DO jl = 1, jpl   
-
-          ! v_ip(:,:,jl) ! Initialize things
-          ! a_ip(:,:,jl)
-          ! volpn(:,:) = hpnd(:,:) * apnd(:,:) * aicen(:,:)
-
-          !------------------------------------------------------------------------------
-          ! 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                 ! ji
-         END DO                    ! jj
-
-         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_fw ) & !--- Give freshwater to the ocean
-                  wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + v_ip(ji,jj,jl) 
-
-
-            ELSE                        !--- Update pond characteristics
-
-               !--- Add retained melt water to melt ponds
-               v_ip(ji,jj,jl)      = v_ip(ji,jj,jl) + 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 / zTp )
-
-               !--- Dump meltwater due to refreezing ( of course this is wrong
-               !--- but this parameterization is too simple )
-               IF ( ln_pnd_fw ) &
-                  wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice
-
-               a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) )
-
-               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)
-
-            !-----------------------------------------------------------
-            ! Limit pond depth
-            !-----------------------------------------------------------
-            ! hpondn = min(hpondn, dpthhi*hi)
-
-            !--- Give freshwater to the ocean ?
-
-            ENDIF
-
-          END DO
-
-       END DO ! jpl
-
-       !--- Remove retained meltwater from surface fluxes 
-
-       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(:,:) )
-
-       ENDIF
+        CALL wrk_alloc( jpi*jpj, indxi, indxj)
+        CALL wrk_alloc( jpi,jpj,     zwfx_mlw )
+        CALL wrk_alloc( jpi,jpj,jpl, zrfrac   )
+
+        z1_rhofw       = 1. / rhofw 
+        z1_zpnd_aspect = 1. / zpnd_aspect
+        zTp            = -2. 
+
+        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(:,:), 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
+ 
+        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                 ! ji
+          END DO                    ! jj
+ 
+          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_fw ) & !--- Give freshwater to the ocean
+                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + v_ip(ji,jj,jl) 
+ 
+ 
+             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 / zTp )
+ 
+                !--- Dump meltwater due to refreezing ( of course this is wrong
+                !--- but this parameterization is too simple )
+                IF ( ln_pnd_fw ) &
+                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice
+ 
+                a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) )
+ 
+                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)
+ 
+             !-----------------------------------------------------------
+             ! Limit pond depth
+             !-----------------------------------------------------------
+             ! The original version has pond depth limitation, which I did not
+             ! keep here. Maybe we need it later on
+             !
+             
+             ENDIF
+ 
+           END DO
+ 
+        END DO ! jpl
+ 
+        !--- Remove retained meltwater from surface fluxes 
+ 
+        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(:,:) )
+ 
+        ENDIF
+
+       CALL wrk_dealloc( jpi*jpj, indxi, indxj)
+       CALL wrk_dealloc( jpi,jpj,     zwfx_mlw )
+       CALL wrk_dealloc( jpi,jpj,jpl, zrfrac   )
 
    END SUBROUTINE lim_mp_cesm
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 8141)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 8142)
@@ -124,8 +124,8 @@
       ELSEWHERE                       ;  zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b
       END WHERE
-      IF( iom_use('alb_ice' ) )  CALL iom_put( "alb_ice"  , zalb(:,:) )           ! ice albedo output
+      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 - at_i_b )      
-      IF( iom_use('albedo'  ) )  CALL iom_put( "albedo"  , zalb(:,:) )           ! ice albedo output
+      IF( iom_use('albedo'  ) )        CALL iom_put( "albedo"  , zalb(:,:) )           ! surface albedo output
 
       CALL wrk_dealloc( jpi,jpj, zalb )    
@@ -174,4 +174,7 @@
             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_fw ) &
+               wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj)
 
             ! mass flux at the ocean/ice interface
Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 8141)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 8142)
@@ -566,5 +566,5 @@
                zvs  = v_s  (ji,jj,  jl)
                zes  = e_s  (ji,jj,1,jl)
-               IF ( ln_pnd ) zvp  = v_ip (ji,jj  ,jl)
+               IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw )   zvp  = v_ip (ji,jj  ,jl)
                !-----------------------------------------------------------------
                ! Zap snow energy 
@@ -588,5 +588,5 @@
                   a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch
                   v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch
-                  wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_ip(ji,jj,jl)  - zvp  ) * rhofw * r1_rdtice
+                  IF ( ln_pnd_fw )   wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_ip(ji,jj,jl)  - zvp  ) * rhofw * r1_rdtice
                ENDIF
                ! END MV MP 2016
Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8141)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8142)
@@ -236,6 +236,8 @@
       ! MV MP 2016
       IF ( ln_pnd ) THEN
-         IF ( iom_use( "iceamp_cat"  ) )  CALL iom_put( "iceamp_cat"     , a_ip   * zswi2   )       ! melt pond frac for categories
-         IF ( iom_use( "icevmp_cat"  ) )  CALL iom_put( "icevmp_cat"     , v_ip   * zswi2   )       ! melt pond frac for categories
+         IF ( iom_use( "iceamp_cat"  ) )  CALL iom_put( "iceamp_cat"     , a_ip       * zswi2   )       ! melt pond frac for categories
+         IF ( iom_use( "icevmp_cat"  ) )  CALL iom_put( "icevmp_cat"     , v_ip       * zswi2   )       ! melt pond frac for categories
+         IF ( iom_use( "icehmp_cat"  ) )  CALL iom_put( "icehmp_cat"     , h_ip       * zswi2   )       ! melt pond frac for categories
+         IF ( iom_use( "iceafp_cat"  ) )  CALL iom_put( "iceafp_cat"     , a_ip_frac  * zswi2   )       ! melt pond frac for categories
       ENDIF
       ! END MV MP 2016
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 8141)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/albedo.F90	(revision 8142)
@@ -110,8 +110,4 @@
       !!---------------------------------------------------------------------
 
-      WRITE(numout,*) ' albedo_ice '
-      WRITE(numout,*) ' ld_pnd     : ', ld_pnd
-      WRITE(numout,*) ' nn_ice_alb : ', nn_ice_alb
-
       ijpl = SIZE( pt_ice, 3 )                     ! number of ice categories
       
@@ -326,6 +322,10 @@
 
           ! Depth-dependent ponded ice albedo
-          zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015)
-          zalb_pnd  = ralb_dp - ( ralb_dp - zalb_ice ) * EXP( - ph_pnd / zhref_pnd ) 
+          IF ( ld_pnd ) THEN
+             zhref_pnd = 0.05        ! Characteristic length scale for thickness dependence of ponded ice albedo, Lecomte et al (2015)
+             zalb_pnd  = ralb_dp - ( ralb_dp - zalb_ice ) * EXP( - ph_pnd / zhref_pnd ) 
+          ELSE
+             zalb_pnd  = ralb_dp
+          ENDIF
 
           ! Surface albedo is weighted mean of snow, ponds and bare ice contributions
