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 8060)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8061)
@@ -83,5 +83,5 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicemp/  ln_limMP, nn_limMP
+      NAMELIST/namicemp/  ln_pnd, nn_pnd_scheme, nn_pnd_cpl, rn_apnd
       !!-------------------------------------------------------------------
 
@@ -99,9 +99,13 @@
          WRITE(numout,*) 'lim_mp_init : ice parameters for melt ponds'
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*)'    Activate melt ponds                                         ln_limMP      = ', ln_limMP  
-         WRITE(numout,*)'    Type of melt pond implementation (1=all,2=rad, 3=fw, 4=no)  nn_limMP      = ', nn_limMP  
+         WRITE(numout,*)'    Activate melt ponds                                         ln_pnd        = ', ln_pnd
+         WRITE(numout,*)'    Type of melt pond scheme =0 presc, =1 empirical = 2 topo    nn_pnd_scheme = ', nn_pnd_scheme
+         WRITE(numout,*)'    Type of melt pond coupling =0 pass., =1 full, =2 rad, 3=fw  nn_pnd_cpl    = ', nn_pnd_cp
+         WRITE(numout,*)'    Prescribed pond fraction                                    rn_apnd       = ', rn_apnd
       ENDIF
       !
    END SUBROUTINE lim_mp_init
+
+
 
    SUBROUTINE lim_mp( kt )
@@ -128,10 +132,17 @@
       IF( nn_timing == 1 )  CALL timing_start('lim_mp')
 
-      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)
-
+      SELECT CASE ( nn_pnd_scheme )
+
+          CASE (1)
+
+             CALL lim_mp_cesm ! empirical melt ponds
+
+          CASE (2)
+
+             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)
+
+      END SELECT
 
       ! we should probably not aggregate here since we do it in lim_var_agg
@@ -142,4 +153,6 @@
 
    END SUBROUTINE lim_mp 
+
+
 
    SUBROUTINE lim_mp_cesm
@@ -149,5 +162,7 @@
        !! ** Purpose    : Compute melt pond evolution
        !!
-       !! ** Method     : Accumulation of meltwater and exponential release
+       !! ** Method     : Empirical method. A fraction of meltwater is accumulated 
+       !!                 in pond volume. It is then released exponentially when
+       !!                 surface is freezingAccumulation of meltwater and exponential release
        !!
        !! ** Tunable parameters :
@@ -155,5 +170,5 @@
        !! 
        !! ** Note       : Stolen from CICE for quick test of the melt pond
-       !!                 interface
+       !!                 radiation and freshwater interfaces
        !!
        !! ** References : Holland, M. M. et al (J Clim 2012)
@@ -161,6 +176,6 @@
        !!-------------------------------------------------------------------
 
-       INTEGER, POINTER, DIMENSION(:,:)    :: indxi             ! compressed indices for cells with ice melting
-       INTEGER, POINTER, DIMENSION(:,:)    :: indxj             !
+       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
@@ -175,17 +190,20 @@
        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
 
+       INTEGER                             :: ji, jj, jl, ij    ! loop indices
+       INTEGER                             :: icells            ! size of dummy array
+
        !!-------------------------------------------------------------------
 
-       CALL wrk_alloc( jpi,jpj, zwfx_mlw )
        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 + rt0.
+       zTp            = -2. 
 
        !------------------------------------------------------------------
@@ -215,6 +233,6 @@
                IF ( a_i(ji,jj,jl) > epsi10 ) THEN
                   icells = icells + 1
-                  indxi(icells) = i
-                  indxj(icells) = j
+                  indxi(icells) = ji
+                  indxj(icells) = jj
                ENDIF
             END DO                 ! ji
@@ -231,21 +249,33 @@
             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
+               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 ( ( nn_pnd_cpl .EQ. 1 ) .OR. ( nn_pnd_cpl .EQ. 3 ) ) & !--- 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
-               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)
+               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) , 0. )
+               zdTs                = MAX ( zTp - t_su(ji,jj,jl) + rt0 , 0. )
                
-               v_ip(ji,jj,jl) = v_ip(ji,jj,jl) * EXP( rexp * zdTs / zTp )
-           
+               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 ( ( nn_pnd_cpl .EQ. 1 ) .OR. ( nn_pnd_cpl .EQ. 3 ) ) THEN
+
+                  wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice
+
+               ENDIF
+
                a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) )
 
@@ -263,7 +293,17 @@
             ENDIF
 
-         END DO
-
-      END DO ! jpl
+          END DO
+
+       END DO ! jpl
+
+       !--- Remove retained meltwater from surface fluxes 
+
+       IF ( ( nn_pnd_cpl .EQ. 1 ) .OR. ( nn_pnd_cpl .EQ. 3 ) ) THEN
+
+           wfx_snw(:,:) = wfx_snw(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
+
+           wfx_sum(:,:) = wfx_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
+
+       ENDIF
 
    END SUBROUTINE lim_mp_cesm
