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 7306)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 7319)
@@ -272,8 +272,10 @@
    REAL(wp), PUBLIC ::   rn_por_rdg       !: initial porosity of ridges (0.3 regular value)
    REAL(wp), PUBLIC ::   rn_fsnowrdg      !: fractional snow loss to the ocean during ridging
+   REAL(wp), PUBLIC ::   rn_fpondrdg      !: fractional melt pond loss to the ocean during ridging
    LOGICAL , PUBLIC ::   ln_rafting       !: rafting of ice or not                        
    REAL(wp), PUBLIC ::   rn_hraft         !: threshold thickness (m) for rafting / ridging 
    REAL(wp), PUBLIC ::   rn_craft         !: coefficient for smoothness of the hyperbolic tangent in rafting
    REAL(wp), PUBLIC ::   rn_fsnowrft      !: fractional snow loss to the ocean during ridging
+   REAL(wp), PUBLIC ::   rn_fpondrft      !: fractional snow loss to the ocean during rafting
 
    ! MV MP 2016
@@ -316,4 +318,7 @@
 
    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_pnd     !: melt pond-ocean mass exchange   [kg.m-2.s-1]
+   ! END MV MP 2016
    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]
@@ -498,4 +503,7 @@
       ALLOCATE( t_bo   (jpi,jpj) , frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,     &
          &      wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) , wfx_lam(jpi,jpj) ,     &
+         ! MV MP 2016
+         &      wfx_pnd(jpi,jpj) ,                                                              &
+         ! END MV MP 2016
          &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     &
          &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,     &
Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 7306)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 7319)
@@ -418,6 +418,6 @@
       !
       IF ( ln_limMP ) THEN
-         a_ip(:,:,:) = 0._wp
-         v_ip(:,:,:) = 0._wp
+         a_ip(:,:,:) = 0.1_wp
+         v_ip(:,:,:) = 0.1_wp
          h_ip(:,:,:) = 0._wp
          a_ip_frac(:,:,:) = 0._wp
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 7306)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 7319)
@@ -507,4 +507,7 @@
       REAL(wp), POINTER, DIMENSION(:) ::   ardg1 , ardg2    ! area of ice ridged & new ridges
       REAL(wp), POINTER, DIMENSION(:) ::   vsrdg , esrdg    ! snow volume & energy of ridging ice
+      ! MV MP 2016
+      REAL(wp), POINTER, DIMENSION(:) ::   vprdg            ! pond volume of ridging ice
+      ! END MV MP 2016
       REAL(wp), POINTER, DIMENSION(:) ::   dhr   , dhr2     ! hrmax - hrmin  &  hrmax^2 - hrmin^2
 
@@ -520,4 +523,7 @@
       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
+      ! MV MP 2016
+      REAL(wp), POINTER, DIMENSION(:) ::   vprft            ! pond volume of rafting ice
+      ! END MV MP 2016
       REAL(wp), POINTER, DIMENSION(:) ::   esrft , smrft    ! snow energy & salinity of rafting ice
       REAL(wp), POINTER, DIMENSION(:) ::   oirft1, oirft2   ! ice age of ice rafted
@@ -531,7 +537,13 @@
       CALL wrk_alloc( jpij,        indxi, indxj )
       CALL wrk_alloc( jpij,        zswitch, fvol )
-      CALL wrk_alloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
+      ! 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 )
+      ! END MV MP 2016
       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 )
+      ! 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 )
+      ! END MV MP 2016
       CALL wrk_alloc( jpij,nlay_i, eirft, erdg1, erdg2, ersw )
 
@@ -595,4 +607,9 @@
             vsrdg (ij) = v_s  (ji,jj,  jl1) * afrac(ij)
             esrdg (ij) = e_s  (ji,jj,1,jl1) * afrac(ij)
+            !MV MP 2016
+            IF ( ln_limMP ) THEN
+               vprdg(ij) = v_ip(ji,jj, jl1) * afrac(ij)
+            ENDIF
+            ! END MV MP 2016
             srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij)
             oirdg1(ij) = oa_i (ji,jj,  jl1) * afrac(ij)
@@ -602,4 +619,10 @@
             virft (ij) = v_i  (ji,jj,  jl1) * afrft(ij)
             vsrft (ij) = v_s  (ji,jj,  jl1) * afrft(ij)
+            !MV MP 2016
+            IF ( ln_limMP ) THEN
+               vprft(ij) = v_ip(ji,jj,jl1) * afrft(ij)
+            ENDIF
+            ! END MV MP 2016
+            srdg1 (ij) = smv_i(ji,jj,  jl1) * afrac(ij)
             esrft (ij) = e_s  (ji,jj,1,jl1) * afrft(ij)
             smrft (ij) = smv_i(ji,jj,  jl1) * afrft(ij) 
@@ -637,4 +660,15 @@
                &                                - esrft(ij) * ( 1._wp - rn_fsnowrft ) ) * r1_rdtice        ! heat sink for ocean (<0, W.m-2)
 
+            ! MV MP 2016
+            !------------------------------------------            
+            ! 3.X Put the melt pond water in the ocean
+            !------------------------------------------            
+            !  Place part of the melt pond volume into the ocean. 
+            IF ( ln_limMP .AND. ( ( nn_limMP .EQ. 3 ) .OR. ( nn_limMP .EQ. 1 ) ) )  THEN
+               wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( rhofw * vprdg(ij) * ( 1._wp - rn_fpondrdg )   & 
+               &                                 + rhofw * vprft(ij) * ( 1._wp - rn_fpondrft ) ) * r1_rdtice  ! fresh water source for ocean
+            ENDIF
+            ! END MV MP 2016
+
             !-----------------------------------------------------------------
             ! 3.8 Compute quantities used to apportion ice among categories
@@ -652,4 +686,10 @@
             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)
+
+            ! MV MP 2016
+            IF ( ln_limMP ) THEN
+               v_ip (ji,jj,jl1) = v_ip (ji,jj,jl1) - vprdg (ij) - vprft (ij)
+            ENDIF
+            ! END MV MP 2016
 
          END DO
@@ -716,4 +756,10 @@
                e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + ( esrdg (ij) * rn_fsnowrdg * fvol(ij)  +  &
                   &                                        esrft (ij) * rn_fsnowrft * zswitch(ij) )
+               ! MV MP 2016
+               IF ( ln_limMP ) THEN
+                  v_ip (ji,jj,jl2) = v_ip (ji,jj,jl2)  + ( vprdg (ij) * rn_fpondrdg * fvol(ij)  +  &
+                  &                                        vprft (ij) * rn_fpondrft * zswitch(ij) )
+               ENDIF
+               ! END MV MP 2016
 
             END DO
@@ -734,7 +780,12 @@
       CALL wrk_dealloc( jpij,        indxi, indxj )
       CALL wrk_dealloc( jpij,        zswitch, fvol )
-      CALL wrk_dealloc( jpij,        afrac, ardg1, ardg2, vsrdg, esrdg, dhr, dhr2 )
+      ! 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 )
+      ! END MV MP 2016
       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 )
+      ! 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,nlay_i, eirft, erdg1, erdg2, ersw )
       !
@@ -915,5 +966,6 @@
       INTEGER :: ios                 ! Local integer output status for namelist read
       NAMELIST/namiceitdme/ rn_cs, nn_partfun, rn_gstar, rn_astar,             & 
-        &                   ln_ridging, rn_hstar, rn_por_rdg, rn_fsnowrdg, ln_rafting, rn_hraft, rn_craft, rn_fsnowrft
+        &                   ln_ridging, rn_hstar, rn_por_rdg, rn_fsnowrdg, rn_fpondrdg, & 
+                            ln_rafting, rn_hraft, rn_craft,   rn_fsnowrft, rn_fpondrft
       !!-------------------------------------------------------------------
       !
@@ -939,8 +991,10 @@
          WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
          WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
+         WRITE(numout,*)'   Fraction of pond volume conserved during ridging        rn_fpondrdg = ', rn_fpondrdg 
          WRITE(numout,*)'   Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
          WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
          WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft  
          WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
+         WRITE(numout,*)'   Fraction of pond volume conserved during rafting        rn_fpondrft = ', rn_fpondrft 
       ENDIF
       !
Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 7306)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 7319)
@@ -357,4 +357,9 @@
          IF ( a_i(ii,ij,1) > epsi10 .AND. ht_i(ii,ij,1) < rn_himin ) THEN
             a_i (ii,ij,1) = a_i(ii,ij,1) * ht_i(ii,ij,1) / rn_himin 
+            ! MV MP 2016
+            IF ( ln_limMP ) THEN
+               a_ip(ii,ij,1) = a_ip(ii,ij,1) * ht_i(ii,ij,1) / rn_himin
+            ENDIF
+            ! END MV MP 2016
             ht_i(ii,ij,1) = rn_himin
          ENDIF
@@ -485,4 +490,8 @@
       REAL(wp) ::   zdo_aice           ! ice age times volume transferred
       REAL(wp) ::   zdaTsf             ! aicen*Tsfcn transferred
+      ! MV MP 2016 
+      REAL(wp) ::   zdapnd             ! pond fraction transferred
+      REAL(wp) ::   zdvpnd             ! pond volume transferred
+      ! END MV MP 2016 
 
       INTEGER, POINTER, DIMENSION(:) ::   nind_i, nind_j   ! compressed indices for i/j directions
@@ -576,4 +585,23 @@
             zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf
             zaTsfn(ii,ij,jl2)  = zaTsfn(ii,ij,jl2) + zdaTsf 
+
+            ! MV MP 2016 
+            IF ( ln_limMP ) THEN
+            !---------------------
+            ! Pond fraction
+            !---------------------
+               zdapnd             = a_ip(ii,ij,jl1) * zdaice(ii,ij,jl)
+               a_ip(ii,ij,jl1)    = a_ip(ii,ij,jl1) - zdapnd
+               a_ip(ii,ij,jl2)    = a_ip(ii,ij,jl2) + zdapnd
+
+            !---------------------
+            ! Pond volume
+            !---------------------
+               zdvpnd             = v_ip(ii,ij,jl1) * zdaice(ii,ij,jl)
+               v_ip(ii,ij,jl1)    = v_ip(ii,ij,jl1) - zdvpnd
+               v_ip(ii,ij,jl2)    = v_ip(ii,ij,jl2) + zdvpnd
+
+            ENDIF
+            ! END MV MP 2016 
 
          END DO
Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 7306)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 7319)
@@ -76,5 +76,12 @@
       ! --- diffusion --- !
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zhdfptab
-      INTEGER , PARAMETER                    ::   ihdf_vars  = 6 ! Number of variables in which we apply horizontal diffusion
+      ! MV MP 2016
+      ! With melt ponds, we have to diffuse them
+      ! We hard code the number of variables to diffuse
+      ! since we can't put an IF ( ln_limMP ) for a declaration
+      ! ideally, the ihdf_vars should probably be passed as an argument and
+      ! defined somewhere depending on ln_limMP
+      ! END MV MP 2016
+      INTEGER , PARAMETER                    ::   ihdf_vars  = 8 ! Number of variables in which we apply horizontal diffusion
                                                                  !  inside limtrp for each ice category , not counting the 
                                                                  !  variables corresponding to ice_layers 
@@ -86,4 +93,8 @@
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0opw
       REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ice, z0snw, z0ai, z0es , z0smi , z0oi
+      ! MV MP 2016
+      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   z0ap , z0vp
+      REAL(wp) ::   za_old
+      ! END MV MP 2016
       REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   z0ei
       !!
@@ -205,4 +216,10 @@
                CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_s(:,:,jl) )      ! Snow volume
                CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, e_s(:,:,1,jl) )    ! Snow heat content
+               ! MV MP 2016
+               IF ( ln_limMP ) THEN
+                  CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, a_ip(:,:,jl) )  ! Melt pond fraction
+                  CALL lim_adv_umx( kt, zdt, zudy, zvdx, zcu_box, zcv_box, v_ip(:,:,jl) )  ! Melt pond volume
+               ENDIF
+               ! END MV MP 2016
             END DO
          END DO
@@ -222,4 +239,5 @@
          CALL wrk_alloc( jpi,jpj,1,          z0opw )
          CALL wrk_alloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
+         CALL wrk_alloc( jpi,jpj,jpl,        z0ap , z0vp )
          CALL wrk_alloc( jpi,jpj,nlay_i,jpl, z0ei )
          
@@ -246,4 +264,11 @@
                z0ei  (:,:,jk,jl) = e_i  (:,:,jk,jl) * e12t(:,:) ! Ice  heat content
             END DO
+            ! MV MP 2016
+            IF ( ln_limMP ) THEN
+               z0ap  (:,:,jl)  = a_ip (:,:,jl) * e12t(:,:) ! Melt pond fraction
+               z0vp  (:,:,jl)  = v_ip (:,:,jl) * e12t(:,:) ! Melt pond volume
+            ENDIF
+            ! END MV MP 2016
+
          END DO
          
@@ -288,4 +313,14 @@
                         &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
                   END DO
+                  ! MV MP 2016
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction --
+                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   & 
+                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
+                  CALL lim_adv_x( zusnit, u_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &    !--- melt pond volume   --
+                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
+                  CALL lim_adv_y( zusnit, v_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   & 
+                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
+                  ! END MV MP 2016
                END DO
             END DO
@@ -329,4 +364,16 @@
                         &                                         syye(:,:,jk,jl), sxye(:,:,jk,jl) )
                   END DO
+                  ! MV MP 2016
+                  IF ( ln_limMP ) THEN
+                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction ---
+                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
+                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &
+                     &                                         sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
+                     CALL lim_adv_y( zusnit, v_ice, 1._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &   !--- melt pond volume   ---
+                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
+                     CALL lim_adv_x( zusnit, u_ice, 0._wp, zarea, z0vp  (:,:,jl), sxvp (:,:,jl),   &
+                     &                                         sxxvp (:,:,jl), syvp (:,:,jl), syyvp (:,:,jl), sxyvp (:,:,jl)  )
+                  ENDIF
+                  ! END MV MP 2016
                END DO
             END DO
@@ -347,4 +394,10 @@
                e_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e12t(:,:)
             END DO
+            ! MV MP 2016
+            IF ( ln_limMP ) THEN
+               a_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e12t(:,:)
+               v_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e12t(:,:)
+            ENDIF
+            ! END MV MP 2016
          END DO
          
@@ -357,4 +410,7 @@
          CALL wrk_dealloc( jpi,jpj,1,          z0opw )
          CALL wrk_dealloc( jpi,jpj,jpl,        z0ice, z0snw, z0ai, z0es , z0smi , z0oi )
+         ! MV MP 2016
+         CALL wrk_dealloc( jpi,jpj,jpl,        z0ap , z0vp )
+         ! END MV MP 2016
          CALL wrk_dealloc( jpi,jpj,nlay_i,jpl, z0ei )
 
@@ -385,4 +441,10 @@
             zhdfptab(:,:,jm)= oa_i (:,:,  jl); jm = jm + 1
             zhdfptab(:,:,jm)= e_s  (:,:,1,jl); jm = jm + 1
+            ! MV MP 2016
+            IF ( ln_limMP ) THEN
+               zhdfptab(:,:,jm)= a_ip  (:,:,  jl); jm = jm + 1
+               zhdfptab(:,:,jm)= v_ip  (:,:,  jl); jm = jm + 1
+            ENDIF
+            ! END MV MP 2016
             ! Sample of adding more variables to apply lim_hdf (ihdf_vars must be increased)
             !   zhdfptab(:,:,jm) = variable_1 (:,:,1,jl); jm = jm + 1  
@@ -418,4 +480,9 @@
             oa_i (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
             e_s  (:,:,1,jl) = zhdfptab(:,:,jm); jm = jm + 1
+            ! MV MP 2016
+            IF ( ln_limMP ) THEN
+               a_ip (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
+               v_ip (:,:,  jl) = zhdfptab(:,:,jm); jm = jm + 1
+            ENDIF
             ! Sample of adding more variables to apply lim_hdf
             !   variable_1  (:,:,1,jl) = zhdfptab(:,:, jm  ) ; jm + 1 
@@ -470,4 +537,11 @@
                         e_s(ji,jj,1,jl)        = rswitch * e_s(ji,jj,1,jl)
                         e_i(ji,jj,1:nlay_i,jl) = rswitch * e_i(ji,jj,1:nlay_i,jl)
+
+                        ! MV MP 2016
+                        IF ( ln_limMP ) THEN
+                           a_ip (ji,jj,jl)        = rswitch * a_ip (ji,jj,jl)
+                           v_ip (ji,jj,jl)        = rswitch * v_ip (ji,jj,jl)
+                        ENDIF
+                        ! END MV MP 2016
                                                 
                      ENDIF
@@ -482,9 +556,19 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
+               ! MV MP 2016
+               za_old = a_i(ji,jj,jpl)
+               ! END MV MP 2016
                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
+               ! MV MP 2016
+               IF ( ln_limMP ) THEN
+                  ! correct pond fraction to avoid a_ip > a_i
+                  a_ip(ji,jj,jpl) = a_ip(ji,jj,jpl) * a_i(ji,jj,jpl) / MAX( za_old , epsi20 ) * rswitch
+               ENDIF
+               ! END MP 2016
+            END DO
+         END DO
+
 
       ENDIF
@@ -512,4 +596,11 @@
       vt_s(:,:) = SUM( v_s(:,:,:), dim=3 )
       at_i(:,:) = SUM( a_i(:,:,:), dim=3 )
+
+      ! MV MP 2016 (remove once we get rid of a_i_frac and ht_i)
+      IF ( ln_limMP ) THEN
+          at_ip(:,:) = SUM( a_ip(:,:,:), dim = 3 )
+          vt_ip(:,:) = SUM( v_ip(:,:,:), dim = 3 )
+      ENDIF
+      ! END MP 2016
       
       ! --- open water = 1 if at_i=0 --------------------------------
Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 7306)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 7319)
@@ -55,4 +55,7 @@
 
    REAL(wp), PUBLIC ::   rhosn    =  330._wp         !: volumic mass of snow          [kg/m3]
+   ! MV MP 2016
+   REAL(wp), PUBLIC ::   rhofw    = 1000._wp         !: volumic mass of freshwater in melt ponds [kg/m3]
+   ! END MV MP 2016
    REAL(wp), PUBLIC ::   emic     =    0.97_wp       !: emissivity of snow or ice
    REAL(wp), PUBLIC ::   sice     =    6.0_wp        !: salinity of ice               [psu]
@@ -192,4 +195,5 @@
          WRITE(numout,*) '          density of sea ice                        = ', rhoic   , ' kg/m^3'
          WRITE(numout,*) '          density of snow                           = ', rhosn   , ' kg/m^3'
+         WRITE(numout,*) '          density of freshwater (in melt ponds)     = ', rhofw   , ' kg/m^3'
          WRITE(numout,*) '          emissivity of snow or ice                 = ', emic  
          WRITE(numout,*) '          salinity of ice                           = ', sice    , ' psu'
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 7306)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 7319)
@@ -658,4 +658,8 @@
       wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp
       wfx_spr(:,:) = 0._wp   ;   wfx_lam(:,:) = 0._wp  
+
+      ! MV MP 2016
+      wfx_pnd(:,:) = 0._wp
+      ! END MV MP 2016
       
       hfx_thd(:,:) = 0._wp   ;   
