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 7325)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8060)
@@ -128,7 +128,10 @@
       IF( nn_timing == 1 )  CALL timing_start('lim_mp')
 
-      CALL lim_mp_topo    (at_i, a_i,                                       &
-                &          vt_i, v_i, v_s,            t_i, s_i, a_ip_frac,  &
-                &          h_ip,     t_su)
+      CALL lim_mp_cesm ! test routine with hyper simple melt ponds
+
+      !CALL lim_mp_topo    (at_i, a_i,                                       &
+      !          &          vt_i, v_i, v_s,            t_i, s_i, a_ip_frac,  &
+      !          &          h_ip,     t_su)
+
 
       ! we should probably not aggregate here since we do it in lim_var_agg
@@ -139,4 +142,130 @@
 
    END SUBROUTINE lim_mp 
+
+   SUBROUTINE lim_mp_cesm
+       !!-------------------------------------------------------------------
+       !!                ***  ROUTINE lim_mp_cesm  ***
+       !!
+       !! ** Purpose    : Compute melt pond evolution
+       !!
+       !! ** Method     : Accumulation of meltwater and exponential release
+       !!
+       !! ** Tunable parameters :
+       !!                
+       !! 
+       !! ** Note       : Stolen from CICE for quick test of the melt pond
+       !!                 interface
+       !!
+       !! ** References : Holland, M. M. et al (J Clim 2012)
+       !!    
+       !!-------------------------------------------------------------------
+
+       INTEGER, POINTER, DIMENSION(:,:)    :: indxi             ! compressed indices for cells with ice melting
+       INTEGER, POINTER, DIMENSION(:,:)    :: indxj             !
+
+       REAL(wp), POINTER, DIMENSION(:,:)   :: zwfx_mlw          ! available meltwater for melt ponding
+       REAL(wp), POINTER, DIMENSION(:,:,:) :: 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               ! reference temperature
+       REAL(wp)                            :: zrexp             ! rate constant to refreeze melt ponds
+       REAL(wp)                            :: zdTs              ! dummy temperature difference
+       REAL(wp)                            :: z1_rhofw          ! inverse freshwater density
+       REAL(wp)                            :: z1_zpnd_aspect    ! inverse pond aspect ratio
+
+       !!-------------------------------------------------------------------
+
+       CALL wrk_alloc( jpi,jpj, zwfx_mlw )
+       CALL wrk_alloc( jpi*jpj, indxi, indxj)
+
+       z1_rhofw       = 1. / rhofw 
+       z1_zpnd_aspect = 1. / zpnd_aspect
+       zTp            = -2 + rt0.
+
+       !------------------------------------------------------------------
+       ! 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) = i
+                  indxj(icells) = j
+               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
+
+               wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + v_ip(ji,jj,jl) !--- Give freshwater to the ocean
+
+               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
+
+            ELSE                        !--- Update pond characteristics
+
+               !--- Add retained melt water
+               v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + rfrac(ji,jj,jl) * z1_rhofw * zwfx_mlw(ji,jj) * a_i(ji,jj,jl)
+
+               !--- Shrink pond due to refreezing
+               zdTs = MAX ( zTp - t_su(ji,jj,jl) , 0. )
+               
+               v_ip(ji,jj,jl) = v_ip(ji,jj,jl) * EXP( rexp * zdTs / zTp )
+           
+               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
+
+   END SUBROUTINE lim_mp_cesm
 
    SUBROUTINE lim_mp_topo    (aice,      aicen,     &
@@ -304,10 +433,11 @@
           wfx_pnd(i,j) = wfx_pnd(i,j) + zdvn ! update flux from ponds to ocean
  
-          ! mean surface temperature
-          zTavg = z0
-          DO n = 1, jpl
-             zTavg = zTavg + zTsfcn(i,j,n)*aicen(i,j,n)
-          END DO
-          zTavg = zTavg / aice(i,j)
+          ! mean surface temperature MV - why do we need that ? --> for the lid
+
+          ! zTavg = z0
+          ! DO n = 1, jpl
+          !   zTavg = zTavg + zTsfcn(i,j,n)*aicen(i,j,n)
+          ! END DO
+          ! zTavg = zTavg / aice(i,j)
  
        END DO ! ij
@@ -438,4 +568,5 @@
                                   * (-0.024_wp*hicen(n) + 0.832_wp) 
              asnon(n) = reduced_aicen(n)  ! effective snow fraction (empirical)
+             ! MV should check whether this makes sense to have the same effective snow fraction in here
           END IF
  
@@ -448,7 +579,13 @@
  ! alfan = 60% of the ice volume) in each category lies above the reference line, 
  ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
- 
-          alfan(n) = p6 * hicen(n)
-          betan(n) = p4 * hicen(n)
+
+ ! MV:  
+ ! Note that this choice is not in the original FF07 paper and has been adopted in CICE
+ ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe
+
+ ! Where does that choice come from
+ 
+          alfan(n) = 0.6 * hicen(n)
+          betan(n) = 0.4 * hicen(n)
         
           cum_max_vol(n)     = z0
@@ -1874,4 +2011,6 @@
    SUBROUTINE lim_mp_topo     ! Empty routine
    END SUBROUTINE lim_mp_topo
+   SUBROUTINE lim_mp_cesm     ! Empty routine
+   END SUBROUTINE lim_mp_cesm
    SUBROUTINE lim_mp_area     ! Empty routine
    END SUBROUTINE lim_mp_area
