Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml	(revision 8597)
@@ -132,5 +132,5 @@
          <field id="icevmp"       long_name="melt pond volume"                                             standard_name="sea_ice_meltpond_volume"                            unit="m"            /> 
 
-         <field id="iceconc_cat"  long_name="Sea-ice area fractions in thickness categories"               unit=""   grid_ref="grid_T_3D_ncatice" />
+         <field id="iceconc_cat"  long_name="Sea-ice concentration in thickness categories"                unit=""   grid_ref="grid_T_3D_ncatice" />
          <field id="icethic_cat"  long_name="Sea-ice thickness in thickness categories"                    unit="m"  grid_ref="grid_T_3D_ncatice" />
          <field id="snowthic_cat" long_name="Snow thickness in thickness categories"                       unit="m"  grid_ref="grid_T_3D_ncatice" />
@@ -140,6 +140,8 @@
          <field id="icetemp_cat"  long_name="Ice temperature for categories"                               unit="degC"   grid_ref="grid_T_3D_ncatice" />
          <field id="snwtemp_cat"  long_name="Snow temperature for categories"                              unit="degC"   grid_ref="grid_T_3D_ncatice" />
-         <field id="iceamp_cat"   long_name="Ice melt pond fraction for categories"                        unit="%"      grid_ref="grid_T_3D_ncatice" /> 
+         <field id="iceamp_cat"   long_name="Ice melt pond concentration for categories"                   unit="%"      grid_ref="grid_T_3D_ncatice" /> 
          <field id="icevmp_cat"   long_name="Ice melt pond volume for categories"                          unit="m"      grid_ref="grid_T_3D_ncatice" /> 
+         <field id="icehmp_cat"   long_name="Ice melt pond thickness for categories"                       unit="m"      grid_ref="grid_T_3D_ncatice" /> 
+         <field id="iceafp_cat"   long_name="Ice melt pond fraction for categories"                        unit="m"      grid_ref="grid_T_3D_ncatice" /> 
 
          <field id="micet"        long_name="Mean ice temperature"                                         unit="degC"     />
@@ -186,4 +188,5 @@
          <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="kg/m2/s"   />
          <field id="vfxthin"      long_name="thermo ice prod. for thin ice(20cm) + open water"             unit="kg/m2/s"   />
+         <field id="vfxpnd"       long_name="melt pond water flux to ocean"                                unit="kg/m2/s"   />
 
          <field id="afxtot"       long_name="area tendency (total)"                                        unit="s-1"   />
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8597)
@@ -125,5 +125,5 @@
    ln_icedA         = .true.          !  activate lateral melting param. (T) or not (F)
    ln_icedO         = .true.          !  activate ice growth in open-water (T) or not (F)
-   ln_icedS         = .true.          !  activate gravity drainage and flushing (T) or not (F)
+   ln_icedS         = .true.          !  activate brine drainage (T) or not (F)
 /
 !------------------------------------------------------------------------------
@@ -176,10 +176,10 @@
 &nammp      !   Melt ponds
 !------------------------------------------------------------------------------
-   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)
+   ln_pnd_H12       = .false.         !  activate evolutive melt ponds (from Holland et al 2012)
+      ln_pnd_fwb    = .false.         !     melt ponds store freshwater or not
+   ln_pnd_CST       = .false.         !  activate constant melt ponds
+      rn_apnd       =   0.2           !     prescribed pond fraction, at Tsu=0
+      rn_hpnd       =   0.05          !     prescribed pond depth, at Tsu=0
+   ln_pnd_alb       = .false.         !  melt ponds affect albedo or not
 /
 !------------------------------------------------------------------------------
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8597)
@@ -204,13 +204,12 @@
    REAL(wp), PUBLIC ::   rn_simin         !: minimum ice salinity [PSU]
 
-   ! MV MP 2016
-   !                                     !!** melt pond namelist (nammp)
-   LOGICAL , PUBLIC ::   ln_pnd           !: activate ponds or not
-   LOGICAL , PUBLIC ::   ln_pnd_rad       !: ponds radiatively active or not
-   LOGICAL , PUBLIC ::   ln_pnd_fw        !: ponds active wrt meltwater or not
-   INTEGER , PUBLIC ::   nn_pnd_scheme    !: type of melt pond scheme:   =0 prescribed, =1 empirical, =2 topographic
-   REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1), only if nn_pnd_scheme = 0
-   REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1), only if nn_pnd_scheme = 0
-   ! END MV MP 2016
+   !                                     !!** namelist namthd_pnd
+   LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al 2012
+   LOGICAL , PUBLIC ::   ln_pnd_fwb       !: melt ponds store freshwater
+   LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth
+   REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1)
+   REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1)
+   LOGICAL , PUBLIC ::   ln_pnd_alb       !: melt ponds affect albedo
+
    !                                     !!** ice-diagnostics namelist (namdia) **
    LOGICAL , PUBLIC ::   ln_icediachk     !: flag for ice diag (T) or not (F)
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90	(revision 8597)
@@ -49,5 +49,5 @@
 CONTAINS
 
-   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, pafrac_pnd, ph_pnd, palb_cs, palb_os )
+   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, palb_cs, palb_os )
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE ice_alb  ***
@@ -98,4 +98,5 @@
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice       !  sea-ice thickness
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw       !  snow depth
+      LOGICAL , INTENT(in   )                   ::   ld_pnd_alb   !  effect of melt ponds on albedo
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area)
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth
@@ -125,5 +126,9 @@
                IF( ph_snw(ji,jj,jl) == 0._wp ) THEN
                   zafrac_snw = 0._wp
-                  zafrac_pnd = pafrac_pnd(ji,jj,jl)
+                  IF( ld_pnd_alb ) THEN
+                     zafrac_pnd = pafrac_pnd(ji,jj,jl)
+                  ELSE
+                     zafrac_pnd = 0._wp
+                  ENDIF
                   zafrac_ice = 1._wp - zafrac_pnd
                ELSE
@@ -134,5 +139,5 @@
                !
                !                       !--- Bare ice albedo (for hi > 150cm)
-               IF( zafrac_pnd > 0._wp ) THEN
+               IF( ld_pnd_alb ) THEN
                   zalb_ice = rn_alb_idry
                ELSE
@@ -154,5 +159,5 @@
                ENDIF
                !                       !--- Ponded ice albedo
-               IF( zafrac_pnd > 0._wp ) THEN
+               IF( ld_pnd_alb ) THEN
                   zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd(ji,jj,jl) * z1_href_pnd ) 
                ELSE
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90	(revision 8597)
@@ -172,9 +172,7 @@
             END DO
          END DO
-      END DO
-            
-      IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i
-         WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
-      ENDIF
+      END DO 
+      !                                           !-- correct pond fraction to avoid a_ip > a_i
+      WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:)
       !
    END SUBROUTINE Hbig
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90	(revision 8597)
@@ -120,5 +120,5 @@
             z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content
          END DO
-         IF ( nn_pnd_scheme > 0 ) THEN
+         IF ( ln_pnd_H12 ) THEN
             z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction
             z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume
@@ -167,5 +167,5 @@
                      &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) )
                END DO
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction --
                      &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
@@ -220,5 +220,5 @@
                      &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) )
                END DO
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction ---
                      &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  )
@@ -249,5 +249,5 @@
          END DO
          ! MV MP 2016
-         IF ( nn_pnd_scheme > 0 ) THEN
+         IF ( ln_pnd_H12 ) THEN
             pa_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:)
             pv_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:)
@@ -755,5 +755,5 @@
                   sxyage(:,:,jl)= z2d(:,:)
                END DO
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   DO jl = 1, jpl 
                      WRITE(zchar,'(I2.2)') jl
@@ -833,5 +833,5 @@
                syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp
                sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
                   syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
@@ -854,5 +854,5 @@
             syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp
             sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp
-            IF ( nn_pnd_scheme > 0 ) THEN
+            IF ( ln_pnd_H12 ) THEN
                sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp 
                syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp 
@@ -989,5 +989,5 @@
             END DO
          END DO
-         IF ( nn_pnd_scheme > 0 ) THEN
+         IF ( ln_pnd_H12 ) THEN
             DO jl = 1, jpl 
                WRITE(zchar,'(I2.2)') jl
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_umx.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_umx.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_umx.F90	(revision 8597)
@@ -124,5 +124,5 @@
             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume
             CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content
-            IF ( nn_pnd_scheme > 0 ) THEN
+            IF ( ln_pnd_H12 ) THEN
                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction
                CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90	(revision 8597)
@@ -253,8 +253,6 @@
          CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
          CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
-         IF ( nn_pnd_scheme > 0 ) THEN
-            CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
-            CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
-         ENDIF
+         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
          DO jl = 1, jpl
             DO jk = 1, nlay_s
@@ -270,7 +268,5 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
          CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
-         IF ( nn_pnd_scheme > 0 ) THEN
-            CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
-         ENDIF
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
          
          !-----------------------------------------------------------------------------!
@@ -320,8 +316,6 @@
          CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) )
          CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) )
-         IF ( nn_pnd_scheme > 0 ) THEN
-            CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
-            CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
-         ENDIF
+         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
          DO jl = 1, jpl
             DO jk = 1, nlay_s
@@ -337,7 +331,5 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) )
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) )
-         IF ( nn_pnd_scheme > 0 ) THEN
-            CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
-         ENDIF
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) )
 
       ENDIF ! npti > 0
@@ -651,5 +643,5 @@
 
                !MV MP 2016
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   aprdg1     = a_ip_2d(ji,jl1) * afrdg
                   aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
@@ -679,5 +671,5 @@
                ! Put the melt pond water into the ocean
                !------------------------------------------            
-               IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN
+               IF ( ln_pnd_fwb ) THEN
                   wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean
                      &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice
@@ -702,5 +694,5 @@
                oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1
                ! MV MP 2016
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
                   v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
@@ -766,5 +758,5 @@
                      &                                    vsrft (ji) * rn_fsnwrft * zswitch(ji) )
                   ! MV MP 2016
-                  IF ( nn_pnd_scheme > 0 ) THEN
+                  IF ( ln_pnd_H12 ) THEN
                      v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   &
                         &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   )
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90	(revision 8597)
@@ -130,5 +130,5 @@
 
       ! --- cloud-sky and overcast-sky ice albedos --- !
-      CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, zalb_cs, zalb_os )
+      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os )
 
       ! albedo depends on cloud fraction because of non-linear spectral effects
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8597)
@@ -94,5 +94,5 @@
       REAL(wp) ::   ztmelts, zdh
       INTEGER  ::   i_hemis, i_fill, jl0
-      REAL(wp) ::   zarg, zV, zconv, zdv
+      REAL(wp) ::   zarg, zV, zconv, zdv, zfac
       INTEGER , DIMENSION(4)           ::   itest
       REAL(wp), DIMENSION(jpi,jpj)     ::   z2d
@@ -358,20 +358,13 @@
          tn_ice (:,:,:) = t_su (:,:,:)
 
-         ! MV MP 2016
          ! Melt pond volume and fraction
-         IF ( ln_pnd ) THEN
-            DO jl = 1, jpl
-               a_ip_frac(:,:,jl) = 0.2 *  zswitch(:,:)
-               h_ip     (:,:,jl) = 0.05 * zswitch(:,:)
-               a_ip(:,:,jl)      = a_ip_frac(:,:,jl) * a_i (:,:,jl) 
-               v_ip(:,:,jl)      = h_ip     (:,:,jl) * a_ip(:,:,jl)
-            END DO
-         ELSE
-            a_ip(:,:,:)      = 0._wp
-            v_ip(:,:,:)      = 0._wp
-            a_ip_frac(:,:,:) = 0._wp
-            h_ip     (:,:,:) = 0._wp
-         ENDIF
-         ! END MV MP 2016
+         IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN   ;   zfac = 1._wp
+         ELSE                                     ;   zfac = 0._wp   ;   ENDIF 
+         DO jl = 1, jpl
+            a_ip_frac(:,:,jl) = rn_apnd * zswitch(:,:) * zfac
+            h_ip     (:,:,jl) = rn_hpnd * zswitch(:,:) * zfac
+         END DO
+         a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:) 
+         v_ip(:,:,:) = h_ip     (:,:,:) * a_ip(:,:,:)
 
       ELSE ! if ln_iceini=false
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8597)
@@ -293,9 +293,5 @@
             IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN
                a_i_1d (ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 
-               ! MV MP 2016
-               IF ( nn_pnd_scheme > 0 ) THEN
-                  a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
-               ENDIF
-               ! END MV MP 2016
+               IF ( ln_pnd_H12 )    a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin
                h_i_1d(ji) = rn_himin
             ENDIF
@@ -457,6 +453,5 @@
                zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans
                !  
-               ! MV MP 2016 
-               IF ( nn_pnd_scheme > 0 ) THEN
+               IF ( ln_pnd_H12 ) THEN
                   !                                                  ! Pond fraction
                   ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work
@@ -468,5 +463,4 @@
                   v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
                ENDIF
-               ! END MV MP 2016
                !
             ENDIF   ! jl1 >0
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90	(revision 8597)
@@ -189,6 +189,7 @@
       !! ** purpose  :   read restart file
       !!----------------------------------------------------------------------
-      INTEGER  :: jk, jl
-      REAL(wp) :: zfice, ziter
+      INTEGER  ::   jk, jl
+      INTEGER  ::   id1            ! local integer
+      REAL(wp) ::   zfice, ziter
       REAL(wp), DIMENSION(jpi,jpj) ::   z2d
       CHARACTER(len=25) ::   znam
@@ -244,13 +245,20 @@
       END DO
 
-      DO jl = 1, jpl 
-         WRITE(zchar,'(I2.2)') jl
-         znam = 'a_ip'//'_htc'//zchar
-         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-         a_ip(:,:,jl) = z2d(:,:)
-         znam = 'v_ip'//'_htc'//zchar
-         CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
-         v_ip(:,:,jl) = z2d(:,:)
-      END DO
+      id1 = iom_varid( numrir, 'a_ip_htc01' , ldstop = .FALSE. )
+      IF( id1 > 0 ) THEN                       ! fields exist (melt ponds)
+         DO jl = 1, jpl 
+            WRITE(zchar,'(I2.2)') jl
+            znam = 'a_ip'//'_htc'//zchar
+            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+            a_ip(:,:,jl) = z2d(:,:)
+            znam = 'v_ip'//'_htc'//zchar
+            CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
+            v_ip(:,:,jl) = z2d(:,:)
+         END DO
+      ELSE                                     ! start from rest
+         IF(lwp) WRITE(numout,*) '   ==>>   previous run without melt ponds output then set it'
+         a_ip(:,:,:) = 0._wp
+         v_ip(:,:,:) = 0._wp
+      ENDIF
 
       DO jl = 1, jpl 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8597)
@@ -170,5 +170,5 @@
          IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics      
          !
-         IF ( ln_pnd )                  CALL lim_mp( kt )             ! -- Melt ponds
+                                        CALL lim_mp( kt )             ! -- Melt ponds
          !
          IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90	(revision 8597)
@@ -166,5 +166,5 @@
                &           + 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)
+            IF ( ln_pnd_fwb )   wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj)
 
             ! add the snow melt water to snow mass flux to the ocean
@@ -199,5 +199,5 @@
       ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
       !------------------------------------------------------------------
-      CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
+      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
       !
       alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
@@ -246,5 +246,5 @@
       CALL iom_put( "vfxlam"     , wfx_lam              )        ! lateral melt 
       CALL iom_put( "vfxice"     , wfx_ice              )        ! total ice growth/melt 
-      IF ( ln_pnd )   CALL iom_put( "vfxpnd", wfx_pnd   )        ! melt pond water flux
+      IF ( ln_pnd_fwb )   CALL iom_put( "vfxpnd", wfx_pnd   )        ! melt pond water flux
 
       IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations  
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8597)
@@ -97,10 +97,6 @@
       et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
 
-      ! MV MP 2016
-      IF ( ln_pnd ) THEN                     ! Melt pond
-         at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 )
-         vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
-      ENDIF
-      ! END MP 2016
+      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds
+      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
 
       ato_i(:,:) = 1._wp - at_i(:,:)    ! open water fraction  
@@ -479,4 +475,7 @@
                wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice
                hfx_res(ji,jj)  = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0
+               IF( ln_pnd_fwb ) THEN 
+                  wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
+               ENDIF
                !-----------------------------------------------------------------
                ! Zap snow energy 
@@ -499,18 +498,10 @@
                h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj)
 
-            END DO
-         END DO
-
-         IF( ln_pnd ) THEN 
-            DO jj = 1 , jpj
-               DO ji = 1 , jpi
-                  IF( ln_pnd_fw )   &
-                     &   wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice
-                  a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
-                  v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
-               END DO
-            END DO
-         ENDIF
-         
+               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
+               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
+
+            END DO
+         END DO
+        
       END DO 
 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90	(revision 8597)
@@ -128,8 +128,6 @@
       CALL iom_put( "snowvol"     , vt_s  * zswi      )        ! snow volume
       
-      IF ( ln_pnd ) THEN
-         CALL iom_put( "iceamp"  , at_ip  * zswi        )   ! melt pond total fraction
-         CALL iom_put( "icevmp"  , vt_ip  * zswi        )   ! melt pond total volume per unit area
-      ENDIF
+      CALL iom_put( "iceamp"  , at_ip  * zswi        )   ! melt pond total fraction
+      CALL iom_put( "icevmp"  , vt_ip  * zswi        )   ! melt pond total volume per unit area
 
       !----------------------------------
@@ -145,10 +143,8 @@
       IF ( iom_use('brinevol_cat') )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 )          ! brine volume
 
-      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('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
+      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
 
       !--------------------------------
@@ -316,14 +312,8 @@
       CALL histdef( kid, "sidive", "Ice divergence"          , "10-8s-1",   &
       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
-
-      ! MV MP 2016
-      IF ( ln_pnd ) THEN
-         CALL histdef( kid, "si_amp", "Melt pond fraction"      , "%"      ,   &
-      &         jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-         CALL histdef( kid, "si_vmp", "Melt pond volume"        ,  "m"     ,   &
-      &         jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
-      ENDIF
-      ! END MV MP 2016
-
+      CALL histdef( kid, "si_amp", "Melt pond fraction"      , "%"      ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
+      CALL histdef( kid, "si_vmp", "Melt pond volume"        ,  "m"     ,   &
+      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
       CALL histdef( kid, "vfxbog", "Ice bottom production"   , "m/s"    ,   &
       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )
@@ -370,10 +360,6 @@
       CALL histwrite( kid, "sidive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) )
 
-      ! MV MP 2016
-      IF ( ln_pnd ) THEN
-         CALL histwrite( kid, "si_amp", kt, at_ip         , jpi*jpj, (/1/) )
-         CALL histwrite( kid, "si_vmp", kt, vt_ip         , jpi*jpj, (/1/) )
-      ENDIF
-      ! END MV MP 2016
+      CALL histwrite( kid, "si_amp", kt, at_ip         , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "si_vmp", kt, vt_ip         , jpi*jpj, (/1/) )
 
       CALL histwrite( kid, "vfxbog", kt, wfx_bog        , jpi*jpj, (/1/) )
@@ -384,6 +370,5 @@
       CALL histwrite( kid, "vfxbom", kt, wfx_bom        , jpi*jpj, (/1/) )
       CALL histwrite( kid, "vfxsum", kt, wfx_sum        , jpi*jpj, (/1/) )
-      IF ( ln_pnd ) &
-         CALL histwrite( kid, "vfxpnd", kt, wfx_pnd     , jpi*jpj, (/1/) )
+      CALL histwrite( kid, "vfxpnd", kt, wfx_pnd        , jpi*jpj, (/1/) )
 
       CALL histwrite( kid, "sithicat", kt, h_i         , jpi*jpj*jpl, (/1/) )    
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8592)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8597)
@@ -7,6 +7,4 @@
    !!            1.0  ! 2012    (O. Lecomte) Adaptation for scientific tests (NEMO3.1)
    !!            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
@@ -16,18 +14,8 @@
    !!   lim_mp_init      : some initialization and namelist read
    !!   lim_mp           : main calling routine
-   !!   lim_mp_topo      : main melt pond routine for the "topographic" formulation (FloccoFeltham)
-   !!   lim_mp_area      : ??? compute melt pond fraction per category
-   !!   lim_mp_perm      : computes permeability (should be a FUNCTION!)
-   !!   calc_hpond       : computes melt pond depth
-   !!   permeability_phy : computes permeability
-
    !!----------------------------------------------------------------------
    USE phycst           ! physical constants
    USE dom_oce          ! ocean space and time domain
-!  USE sbc_ice          ! Surface boundary condition: ice   fields
    USE ice              ! LIM-3 variables
-!  USE icecons          ! conservation tests
-!  USE limctl           ! control prints
-!  USE limvar
    !
    USE lbclnk           ! lateral boundary conditions - MPP exchanges
@@ -37,24 +25,4 @@
    USE timing           ! Timing
 
-!OLI_CODE    USE ice_oce, ONLY: rdt_ice, tatm_ice
-!OLI_CODE    USE phycst
-!OLI_CODE    USE dom_ice
-!OLI_CODE    USE dom_oce
-!OLI_CODE    USE sbc_oce
-!OLI_CODE    USE sbc_ice
-!OLI_CODE    USE par_ice
-!OLI_CODE    USE par_oce
-!OLI_CODE    USE ice
-!OLI_CODE    USE thd_ice
-!OLI_CODE    USE in_out_manager
-!OLI_CODE    USE lbclnk
-!OLI_CODE    USE lib_mpp
-!OLI_CODE 
-!OLI_CODE    IMPLICIT NONE
-!OLI_CODE    PRIVATE
-!OLI_CODE 
-!OLI_CODE    PUBLIC   lim_mp_init
-!OLI_CODE    PUBLIC   lim_mp
-
    IMPLICIT NONE
    PRIVATE
@@ -62,4 +30,10 @@
    PUBLIC   lim_mp_init    ! routine called by icestp.F90
    PUBLIC   lim_mp         ! routine called by icestp.F90
+
+   INTEGER ::              nice_pnd   ! choice of the type of pond scheme
+   !                                               ! associated indices:
+   INTEGER, PARAMETER ::   np_pndNO  = 0   ! No pond scheme
+   INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant pond scheme
+   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012)
 
    !! * Substitutions
@@ -71,79 +45,4 @@
    !!----------------------------------------------------------------------
 CONTAINS
-
-   SUBROUTINE lim_mp_init 
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE lim_mp_init   ***
-      !!
-      !! ** Purpose : Physical constants and parameters linked to melt ponds
-      !!      over sea ice
-      !!
-      !! ** Method  :  Read the nammp  namelist and check the melt pond  
-      !!       parameter values called at the first timestep (nit000)
-      !!
-      !! ** input   :   Namelist nammp  
-      !!-------------------------------------------------------------------
-      INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/nammp/  ln_pnd, ln_pnd_rad, ln_pnd_fw, nn_pnd_scheme, rn_apnd, rn_hpnd
-      !!-------------------------------------------------------------------
-
-      REWIND( numnam_ice_ref )              ! Namelist nammp  in reference namelist : Melt Ponds  
-      READ  ( numnam_ice_ref, nammp, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp  in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )              ! Namelist nammp  in configuration namelist : Melt Ponds
-      READ  ( numnam_ice_cfg, nammp, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, nammp )
-      
-      IF(lwp) THEN                        ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_mp_init : ice parameters for melt ponds'
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   Active melt ponds                                           ln_pnd        = ', ln_pnd
-         WRITE(numout,*) '   Active melt ponds radiative coupling                        ln_pnd_rad    = ', ln_pnd_rad
-         WRITE(numout,*) '   Active melt ponds freshwater coupling                       ln_pnd_fw     = ', ln_pnd_fw
-         WRITE(numout,*) '   Type of melt pond scheme =0 presc, =1 empirical = 2 topo    nn_pnd_scheme = ', nn_pnd_scheme
-         WRITE(numout,*) '   Prescribed pond fraction                                    rn_apnd       = ', rn_apnd
-         WRITE(numout,*) '   Prescribed pond depth                                       rn_hpnd       = ', rn_hpnd
-      ENDIF
-
-      IF ( .NOT. ln_pnd ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) ' Melt ponds are not activated '
-         IF(lwp) WRITE(numout,*) ' ln_pnd_rad and ln_pnd_fw set to .FALSE. '
-         IF(lwp) WRITE(numout,*) ' nn_pnd_scheme, rn_apnd, rn_hpnd set to zero '
-         ln_pnd_rad    = .FALSE.
-         ln_pnd_fw     = .FALSE.
-         nn_pnd_scheme = 0
-         rn_apnd       = 0._wp
-         rn_hpnd       = 0._wp
-
-         IF(lwp) THEN                     ! control print
-            WRITE(numout,*) '    Active melt ponds radiative coupling                        ln_pnd_rad    = ', ln_pnd_rad
-            WRITE(numout,*) '    Active melt ponds freshwater coupling                       ln_pnd_fw     = ', ln_pnd_fw
-            WRITE(numout,*) '    Type of melt pond scheme =0 presc, =1 empirical = 2 topo    nn_pnd_scheme = ', nn_pnd_scheme
-            WRITE(numout,*) '    Prescribed pond fraction                                    rn_apnd       = ', rn_apnd
-            WRITE(numout,*) '    Prescribed pond depth                                       rn_hpnd       = ', rn_hpnd
-         ENDIF
-      ENDIF
-
-      IF ( ln_pnd .AND. ( nn_pnd_scheme == 0 ) .AND. ( ln_pnd_fw ) ) THEN
-         IF(lwp) 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
-         IF(lwp) WRITE(numout,*) ' Topographic melt ponds are incompatible with jpl = 1 '
-         IF(lwp) WRITE(numout,*) ' Run aborted '
-         CALL ctl_stop( 'STOP', 'lim_mp_init: uncompatible options, reset namelist_ice_ref ' )
-      ENDIF
-
-      !
-   END SUBROUTINE lim_mp_init
-   
 
    SUBROUTINE lim_mp( kt )
@@ -158,29 +57,20 @@
       !!              - 
       !!------------------------------------------------------------------------------------
-
-      INTEGER, INTENT(in) :: kt    ! number of iteration
-      INTEGER  ::   ji, jj, jl     ! dummy loop indices
-
+      INTEGER, INTENT(in) ::   kt    ! number of iteration
+      INTEGER ::   ji, jj, jl     ! dummy loop indices
       !!-------------------------------------------------------------------
 
       IF( nn_timing == 1 )  CALL timing_start('lim_mp')
 
-      SELECT CASE ( nn_pnd_scheme )
-
-         CASE (0)
-
-            CALL lim_mp_cstt       ! staircase melt ponds
-
-         CASE (1)
-
-            CALL lim_mp_cesm       ! empirical melt ponds
-
-         CASE (2)
-
-            CALL lim_mp_topo   &   ! topographic melt ponds 
-                      &          (at_i, a_i,                                       &
-                      &          vt_i, v_i, v_s,            t_i, sz_i, a_ip_frac,   &
-                      &          h_ip,     t_su)
-
+      SELECT CASE ( nice_pnd )
+
+      CASE (np_pndCST)
+         !                             !------------------------------!
+         CALL lim_mp_CST               ! Constant melt ponds          !
+         !                             !------------------------------!
+      CASE (np_pndH12)
+         !                             !------------------------------!
+         CALL lim_mp_H12               ! Holland et al 2012 melt ponds!
+         !                             !------------------------------!
       END SELECT
 
@@ -189,7 +79,7 @@
    END SUBROUTINE lim_mp 
 
-   SUBROUTINE lim_mp_cstt 
-       !!-------------------------------------------------------------------
-       !!                ***  ROUTINE lim_mp_cstt  ***
+   SUBROUTINE lim_mp_CST 
+       !!-------------------------------------------------------------------
+       !!                ***  ROUTINE lim_mp_CST  ***
        !!
        !! ** Purpose    : Compute melt pond evolution
@@ -206,11 +96,7 @@
        !!    
        !!-------------------------------------------------------------------
-       INTEGER                             :: ji, jj, jl
-       REAL(wp)                            :: z1_jpl            ! 1/jpl
-       !!-------------------------------------------------------------------
-
-       z1_jpl     = 1. / FLOAT(jpl)
 
        WHERE ( ( a_i > epsi10 ) .AND. ( t_su >= rt0-epsi06 ) ) 
+!!       WHERE ( ( a_i > 0._wp ) .AND. ( t_su >= rt0 ) ) 
           a_ip_frac = rn_apnd
           h_ip      = rn_hpnd    
@@ -226,9 +112,9 @@
        wfx_pnd(:,:) = 0._wp
 
-   END SUBROUTINE lim_mp_cstt
-
-   SUBROUTINE lim_mp_cesm
-       !!-------------------------------------------------------------------
-       !!                ***  ROUTINE lim_mp_cesm  ***
+   END SUBROUTINE lim_mp_CST
+
+   SUBROUTINE lim_mp_H12
+       !!-------------------------------------------------------------------
+       !!                ***  ROUTINE lim_mp_H12  ***
        !!
        !! ** Purpose    : Compute melt pond evolution
@@ -249,25 +135,25 @@
        !!-------------------------------------------------------------------
 
-       INTEGER, DIMENSION(jpij)         :: indxi             ! compressed indices for cells with ice melting
-       INTEGER, DIMENSION(jpij)         :: indxj             !
-
-       REAL(wp), DIMENSION(jpi,jpj)     :: zwfx_mlw          ! available meltwater for melt ponding
-       REAL(wp), DIMENSION(jpi,jpj,jpl) :: 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) :: zdTs              ! dummy temperature difference
-       REAL(wp) :: z1_rhofw          ! inverse freshwater density
-       REAL(wp) :: z1_zpnd_aspect    ! inverse pond aspect ratio
-       REAL(wp) :: zvpold            ! dummy pond volume
-
-       INTEGER  :: ji, jj, jl, ij    ! loop indices
-       INTEGER  :: icells            ! size of dummy array
+       INTEGER, DIMENSION(jpij)         ::   indxi             ! compressed indices for cells with ice melting
+       INTEGER, DIMENSION(jpij)         ::   indxj             !
+
+       REAL(wp), DIMENSION(jpi,jpj)     ::   zwfx_mlw          ! available meltwater for melt ponding
+       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   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) ::   zdTs              ! dummy temperature difference
+       REAL(wp) ::   z1_rhofw          ! inverse freshwater density
+       REAL(wp) ::   z1_zpnd_aspect    ! inverse pond aspect ratio
+       REAL(wp) ::   zvpold            ! dummy pond volume
+
+       INTEGER  ::   ji, jj, jl, ij    ! loop indices
+       INTEGER  ::   icells            ! size of dummy array
        !!-------------------------------------------------------------------
         z1_rhofw       = 1. / rhofw 
@@ -301,9 +187,7 @@
            ! 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
@@ -313,6 +197,6 @@
                    indxj(icells) = jj
                 ENDIF
-             END DO                 ! ji
-          END DO                    ! jj
+             END DO
+          END DO
  
           DO ij = 1, icells
@@ -331,5 +215,5 @@
                 h_ip(ji,jj,jl)      = 0._wp
  
-                IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean
+                IF ( ln_pnd_fwb ) & !--- Give freshwater to the ocean
                    wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + v_ip(ji,jj,jl) 
  
@@ -352,5 +236,5 @@
                 !--- Dump meltwater due to refreezing ( of course this is wrong
                 !--- but this parameterization is too simple )
-                IF ( ln_pnd_fw ) &
+                IF ( ln_pnd_fwb ) &
                    wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice
  
@@ -361,11 +245,4 @@
  
                 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
@@ -377,1762 +254,72 @@
         !--- Remove retained meltwater from surface fluxes 
  
-        IF ( ln_pnd_fw ) THEN
- 
-            wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) ) 
-
-            wfx_sum(:,:)     = wfx_sum(:,:) * ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
+        IF ( ln_pnd_fwb ) THEN
+ 
+           wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) ) 
+           wfx_sum(:,:)     = wfx_sum(:,:)     *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) )
 
         ENDIF
 
-   END SUBROUTINE lim_mp_cesm
-
-   SUBROUTINE lim_mp_topo    (aice,      aicen,     &
-                              vice,      vicen,     &
-                              vsnon,                &
-                              ticen,     salin,     &
-                              a_ip_frac, h_ip,      &
-                                         Tsfc )
-       !!-------------------------------------------------------------------
-       !!                ***  ROUTINE lim_mp_topo  ***
-       !!
-       !! ** Purpose :   Compute melt pond evolution based on the ice 
-       !!                topography as inferred from the ice thickness 
-       !!                distribution.  
-       !!
-       !! ** Method  :   This code is initially based on Flocco and Feltham 
-       !!                (2007) and Flocco et al. (2010). More to come...
-       !!
-       !! ** Tunable parameters :
-       !! 
-       !! ** Note :
-       !!
-       !! ** References
-       !!    Flocco, D. and D. L. Feltham, 2007.  A continuum model of melt pond 
-       !!    evolution on Arctic sea ice.  J. Geophys. Res. 112, C08016, doi: 
-       !!    10.1029/2006JC003836.
-       !!    Flocco, D., D. L. Feltham and A. K. Turner, 2010.  Incorporation of
-       !!    a physically based melt pond scheme into the sea ice component of a
-       !!    climate model.  J. Geophys. Res. 115, C08012, 
-       !!    doi: 10.1029/2009JC005568.
-       !!    
-       !!-------------------------------------------------------------------
- 
-       REAL (wp), DIMENSION (jpi,jpj), &
-          INTENT(IN) :: &
-          aice, &    ! total ice area fraction
-          vice       ! total ice volume (m)
- 
-       REAL (wp), DIMENSION (jpi,jpj,jpl), &
-          INTENT(IN) :: &
-          aicen, &   ! ice area fraction, per category
-          vsnon, &   ! snow volume, per category (m)
-          vicen      ! ice volume, per category (m)
- 
-       REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), &
-          INTENT(IN) :: &
-          ticen, &   ! ice enthalpy, per category
-          salin
- 
-       REAL (wp), DIMENSION (jpi,jpj,jpl), &
-          INTENT(INOUT) :: &
-          a_ip_frac , &   ! pond area fraction of ice, per ice category
-          h_ip       ! pond depth, per ice category (m)
- 
-       REAL (wp), DIMENSION (jpi,jpj,jpl), &
-          INTENT(IN) :: &
-          Tsfc       ! snow/sea ice surface temperature
- 
-       ! local variables
-       REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
-          zTsfcn, & ! ice/snow surface temperature (C)
-          zvolpn, & ! pond volume per unit area, per category (m)
-          zvuin     ! water-equivalent volume of ice lid on melt pond ('upper ice', m)
- 
-       REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
-          zapondn,& ! pond area fraction, per category
-          zhpondn   ! pond depth, per category (m)
- 
-       REAL (wp), DIMENSION (jpi,jpj) :: &
-          zvolp       ! total volume of pond, per unit area of pond (m)
- 
-       REAL (wp) :: &
-          zhi,    & ! ice thickness (m)
-          zdHui,  & ! change in thickness of ice lid (m)
-          zomega, & ! conduction
-          zdTice, & ! temperature difference across ice lid (C)
-          zdvice, & ! change in ice volume (m)
-          zTavg,  & ! mean surface temperature across categories (C)
-          zTp,    & ! pond freezing temperature (C)
-          zdvn      ! change in melt pond volume for fresh water budget
-       INTEGER, DIMENSION (jpi*jpj) :: &
-          indxi, indxj    ! compressed indices for cells with ice melting
- 
-       INTEGER :: n,k,i,j,ij,icells,indxij ! loop indices
- 
-       INTEGER, DIMENSION (jpl) :: &
-          kcells          ! cells where ice lid combines with vice
- 
-       INTEGER, DIMENSION (jpi*jpj,jpl) :: &
-          indxii, indxjj  ! i,j indices for kcells loop
- 
-       REAL (wp), parameter :: &
-          zhicemin  = 0.1_wp , & ! minimum ice thickness with ponds (m) 
-          zTd       = 0.15_wp, & ! temperature difference for freeze-up (C)
-          zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3)
-          zmin_volp = 1.e-4_wp, & ! minimum pond volume (m)
-          z0       = 0._wp,    & 
-          zTimelt   = 0._wp,    &
-          z01      = 0.01_wp,  &
-          z25      = 0.25_wp,  &
-          z5       = 0.5_wp
-
-       !---------------------------------------------------------------
-       ! Initialization
-       !---------------------------------------------------------------
- 
-       zhpondn(:,:,:) = 0._wp
-       zapondn(:,:,:) = 0._wp
-       indxii(:,:) = 0
-       indxjj(:,:) = 0
-       kcells(:)   = 0
-
-       zvolp(:,:) = wfx_sum(:,:) + wfx_snw_sum(:,:) + vt_ip(:,:) ! Total available melt water, to be distributed as melt ponds
-       zTsfcn(:,:,:) = zTsfcn(:,:,:) - rt0                   ! Convert in Celsius
- 
-       ! The freezing temperature for meltponds is assumed slightly below 0C,
-       ! as if meltponds had a little salt in them.  The salt budget is not
-       ! altered for meltponds, but if it were then an actual pond freezing 
-       ! temperature could be computed.
- 
-       ! zTp = zTimelt - zTd  ---> for lids
- 
-       !-----------------------------------------------------------------
-       ! Identify grid cells with ponds
-       !-----------------------------------------------------------------
- 
-       icells = 0
-       DO j = 1, jpj
-       DO i = 1, jpi
-          zhi = z0
-          IF (aice(i,j) > epsi10 ) zhi = vice(i,j)/aice(i,j)
-          IF ( aice(i,j) > z01 .and. zhi > zhicemin .and. &
-             zvolp(i,j) > zmin_volp*aice(i,j)) THEN
-             icells = icells + 1
-             indxi(icells) = i
-             indxj(icells) = j
-          ELSE  ! remove ponds on thin ice
-             !fpond(i,j) = fpond(i,j) - zvolp(i,j)
-             zvolpn(i,j,:) = z0
-             zvuin (i,j,:) = z0
-             zvolp (i,j) = z0
-          END IF
-       END DO                     ! i
-       END DO                     ! j
- 
-       DO ij = 1, icells
-          i = indxi(ij)
-          j = indxj(ij)
- 
-          !--------------------------------------------------------------
-          ! calculate pond area and depth
-          !--------------------------------------------------------------
-          CALL lim_mp_area(aice(i,j),vice(i,j), &
-                    aicen(i,j,:), vicen(i,j,:), vsnon(i,j,:), &
-                    ticen(i,j,:,:), salin(i,j,:,:), &
-                    zvolpn(i,j,:), zvolp(i,j), &
-                    zapondn(i,j,:),zhpondn(i,j,:), zdvn)
-         ! outputs are
-         ! - zdvn
-         ! - zvolpn
-         ! - zvolp
-         ! - zapondn
-         ! - zhpondn
-
-          wfx_pnd(i,j) = wfx_pnd(i,j) + zdvn ! update flux from ponds to ocean
- 
-          ! 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
- 
-       !---------------------------------------------------------------
-       ! Update pond volume and fraction
-       !---------------------------------------------------------------
- 
-       a_ip(:,:,:) = zapondn(:,:,:)
-       v_ip(:,:,:) = zapondn(:,:,:) * zhpondn(:,:,:)
-       a_ip_frac(:,:,:) = 0._wp
-       h_ip     (:,:,:) = 0._wp
-
-    END SUBROUTINE lim_mp_topo
-
-    SUBROUTINE lim_mp_area(aice,vice, &
-                         aicen, vicen, vsnon, ticen, &
-                         salin, zvolpn, zvolp,         &
-                         zapondn,zhpondn,dvolp)
-
-       !!-------------------------------------------------------------------
-       !!                ***  ROUTINE lim_mp_area ***
-       !!
-       !! ** Purpose : Given the total volume of meltwater, update
-       !!              pond fraction (a_ip) and depth (should be volume)
-       !!
-       !! **
-       !!              
-       !!------------------------------------------------------------------
-
-       REAL (wp), INTENT(IN) :: &
-          aice,vice
- 
-       REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
-          aicen, vicen, vsnon
- 
-       REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: &
-          ticen, salin
- 
-       REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: &
-          zvolpn
- 
-       REAL (wp), INTENT(INOUT) :: &
-          zvolp, dvolp
- 
-       REAL (wp), DIMENSION(jpl), INTENT(OUT) :: &
-          zapondn, zhpondn
- 
-       INTEGER :: &
-          n, ns,   &
-          m_index, &
-          permflag
- 
-       REAL (wp), DIMENSION(jpl) :: &
-          hicen, &
-          hsnon, &
-          asnon, &
-          alfan, &
-          betan, &
-          cum_max_vol, &
-          reduced_aicen        
- 
-       REAL (wp), DIMENSION(0:jpl) :: &
-          cum_max_vol_tmp
- 
-       REAL (wp) :: &
-          hpond, &
-          drain, &
-          floe_weight, &
-          pressure_head, &
-          hsl_rel, &
-          deltah, &
-          perm, &
-          msno
- 
-       REAL (wp), parameter :: & 
-          viscosity = 1.79e-3_wp, &  ! kinematic water viscosity in kg/m/s
-          z0        = 0.0_wp    , &
-          c1        = 1.0_wp    , &
-          p4        = 0.4_wp    , &
-          p6        = 0.6_wp    , &
-          epsi10      = 1.0e-11_wp
-          
-      !-----------|
-      !           |
-      !           |-----------|
-      !___________|___________|______________________________________sea-level
-      !           |           |
-      !           |           |---^--------|
-      !           |           |   |        |
-      !           |           |   |        |-----------|              |-------
-      !           |           |   |alfan(n)|           |              |
-      !           |           |   |        |           |--------------|
-      !           |           |   |        |           |              |
-      !---------------------------v-------------------------------------------
-      !           |           |   ^        |           |              |
-      !           |           |   |        |           |--------------|
-      !           |           |   |betan(n)|           |              |
-      !           |           |   |        |-----------|              |-------
-      !           |           |   |        |
-      !           |           |---v------- |
-      !           |           |
-      !           |-----------|
-      !           |
-      !-----------|
-     
-       !-------------------------------------------------------------------
-       ! initialize
-       !-------------------------------------------------------------------
- 
-       DO n = 1, jpl
- 
-          zapondn(n) = z0
-          zhpondn(n) = z0
- 
-          !----------------------------------------
-          ! X) compute the effective snow fraction
-          !----------------------------------------
-          IF (aicen(n) < epsi10)  THEN
-             hicen(n) =  z0 
-             hsnon(n) = z0
-             reduced_aicen(n) = z0
-          ELSE
-             hicen(n) = vicen(n) / aicen(n)
-             hsnon(n) = vsnon(n) / aicen(n)
-             reduced_aicen(n) = c1 ! n=jpl
-             IF (n < jpl) reduced_aicen(n) = aicen(n) &
-                                  * (-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
- 
- ! This choice for alfa and beta ignores hydrostatic equilibium of categories.
- ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
- ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
- ! categories.  alfa and beta partition the ITD - they are areas not thicknesses!
- ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
- ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 
- ! 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.
-
- ! 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
-          cum_max_vol_tmp(n) = z0
-     
-       END DO ! jpl
- 
-       cum_max_vol_tmp(0) = z0
-       drain = z0
-       dvolp = z0
-
-       !----------------------------------------------------------
-       ! x) Drain overflow water, update pond fraction and volume
-       !----------------------------------------------------------
-       
-       !--------------------------------------------------------------------------
-       ! the maximum amount of water that can be contained up to each ice category
-       !--------------------------------------------------------------------------
-
-       ! MV
-       ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW)
-       ! Then the excess volume cum_max_vol(jl) drains out of the system
-       ! It should be added to wfx_pnd
-       ! END MV
-     
-       DO n = 1, jpl-1 ! last category can not hold any volume
- 
-          IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN
- 
-             ! total volume in level including snow
-             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + &
-                (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n)) 
- 
-             ! subtract snow solid volumes from lower categories in current level
-             DO ns = 1, n 
-                cum_max_vol_tmp(n) = cum_max_vol_tmp(n) &
-                   - rhosn/rhofw * &   ! free air fraction that can be filled by water
-                     asnon(ns)  * &    ! effective areal fraction of snow in that category
-                     max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)- &
-                                   alfan(n)), z0)  
-             END DO
-
-          ELSE ! assume higher categories unoccupied
-             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1)
-          END IF
-          !IF (cum_max_vol_tmp(n) < z0) THEN
-          !   call abort_ice('negative melt pond volume')
-          !END IF
-       END DO
-       cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume
-       cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl)
-     
-       !----------------------------------------------------------------
-       ! is there more meltwater than can be held in the floe?
-       !----------------------------------------------------------------
-       IF (zvolp >= cum_max_vol(jpl)) THEN
-          drain = zvolp - cum_max_vol(jpl) + epsi10
-          zvolp = zvolp - drain ! update meltwater volume available
-          dvolp = drain         ! this is the drained water
-          IF (zvolp < epsi10) THEN
-             dvolp = dvolp + zvolp
-             zvolp = z0
-          END IF
-       END IF
-     
-       ! height and area corresponding to the remaining volume
- 
-!      call calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, &
-!           zvolp, cum_max_vol, hpond, m_index)
-     
-       DO n=1, m_index
-          zhpondn(n) = hpond - alfan(n) + alfan(1) ! here oui choulde update
-                                                   !  volume instead, no ?
-          zapondn(n) = reduced_aicen(n) 
-          ! in practise, pond fraction depends on the empirical snow fraction
-          ! so in turn on ice thickness
-       END DO
-     
-       !------------------------------------------------------------------------
-       ! Drainage through brine network (permeability)
-       !------------------------------------------------------------------------
-       !!! drainage due to ice permeability - Darcy's law
-     
-       ! sea water level 
-       msno = z0
-       DO n=1,jpl
-         msno = msno + vsnon(n) * rhosn
-       END DO
-       floe_weight = (msno + rhoic*vice + rau0*zvolp) / aice
-       hsl_rel = floe_weight / rau0 &
-               - ((sum(betan(:)*aicen(:))/aice) + alfan(1))
-     
-       deltah = hpond - hsl_rel
-       pressure_head = grav * rau0 * max(deltah, z0)
- 
-       ! drain IF ice is permeable    
-       permflag = 0
-       IF (pressure_head > z0) THEN
-       DO n = 1, jpl-1
-          IF (hicen(n) /= z0) THEN
-             perm = 0. ! MV ugly dummy patch
-             CALL lim_mp_perm(ticen(:,n), salin(:,n), vicen(n), perm)
-             IF (perm > z0) permflag = 1
-
-             drain = perm*zapondn(n)*pressure_head*rdt_ice / &
-                                      (viscosity*hicen(n))
-             dvolp = dvolp + min(drain, zvolp)
-             zvolp = max(zvolp - drain, z0)
-             IF (zvolp < epsi10) THEN
-                dvolp = dvolp + zvolp
-                zvolp = z0
-             END IF
-          END IF
-       END DO
-  
-       ! adjust melt pond DIMENSIONs
-       IF (permflag > 0) THEN
-          ! recompute pond depth    
-!         CALL calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, &
-!                         zvolp, cum_max_vol, hpond, m_index)
-          DO n=1, m_index
-             zhpondn(n) = hpond - alfan(n) + alfan(1)
-             zapondn(n) = reduced_aicen(n) 
-          END DO
-       END IF
-       END IF ! pressure_head
- 
-       !-------------------------------
-       ! X) remove water from the snow
-       !-------------------------------
-       !------------------------------------------------------------------------
-       ! total melt pond volume in category DOes not include snow volume
-       ! snow in melt ponds is not melted
-       !------------------------------------------------------------------------
- 
-       ! Calculate pond volume for lower categories
-       DO n=1,m_index-1
-          zvolpn(n) = zapondn(n) * zhpondn(n) & ! what is not in the snow
-                   - (rhosn/rhofw) * asnon(n) * min(hsnon(n), zhpondn(n))
-       END DO
- 
-       ! Calculate pond volume for highest category = remaining pond volume
-
-       ! The following is completely unclear to Martin at least
-       ! Could we redefine properly and recode in a more readable way ?
-
-       ! m_index = last category with melt pond
-
-       IF (m_index == 1) zvolpn(m_index) = zvolp ! volume of mw in 1st category is the total volume of melt water
-
-       IF (m_index > 1) THEN
-         IF (zvolp > sum(zvolpn(1:m_index-1))) THEN
-           zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1)) ! 
-         ELSE
-           zvolpn(m_index) = z0
-           zhpondn(m_index) = z0
-           zapondn(m_index) = z0
-           ! If remaining pond volume is negative reduce pond volume of 
-           ! lower category
-           IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) & 
-             zvolpn(m_index-1) = zvolpn(m_index-1)-sum(zvolpn(1:m_index-1))&
-                                + zvolp
-         END IF
-       END IF
- 
-       DO n=1,m_index
-          IF (zapondn(n) > epsi10) THEN
-              zhpondn(n) = zvolpn(n) / zapondn(n)
-          ELSE
-             dvolp = dvolp + zvolpn(n)
-             zhpondn(n) = z0
-             zvolpn(n) = z0
-             zapondn(n) = z0
-          end IF
-       END DO
-       DO n = m_index+1, jpl
-          zhpondn(n) = z0
-          zapondn(n) = z0
-          zvolpn (n) = z0
-       END DO
- 
-    END SUBROUTINE lim_mp_area
-
-!OLI_CODE   
-!OLI_CODE 
-!OLI_CODE    SUBROUTINE calc_hpond(aicen, asnon, hsnon, rhos, alfan, &
-!OLI_CODE                          zvolp, cum_max_vol,                &
-!OLI_CODE                          hpond, m_index)
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE       !!                ***  ROUTINE calc_hpond  ***
-!OLI_CODE       !!
-!OLI_CODE       !! ** Purpose :   Compute melt pond depth
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE     
-!OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
-!OLI_CODE          aicen, &
-!OLI_CODE          asnon, &
-!OLI_CODE          hsnon, &
-!OLI_CODE          rhos,  &
-!OLI_CODE          alfan, &
-!OLI_CODE          cum_max_vol
-!OLI_CODE     
-!OLI_CODE       REAL (wp), INTENT(IN) :: &
-!OLI_CODE          zvolp
-!OLI_CODE     
-!OLI_CODE       REAL (wp), INTENT(OUT) :: &
-!OLI_CODE          hpond
-!OLI_CODE     
-!OLI_CODE       INTEGER, INTENT(OUT) :: &
-!OLI_CODE          m_index
-!OLI_CODE     
-!OLI_CODE       INTEGER :: n, ns
-!OLI_CODE     
-!OLI_CODE       REAL (wp), DIMENSION(0:jpl+1) :: &
-!OLI_CODE          hitl, &
-!OLI_CODE          aicetl
-!OLI_CODE     
-!OLI_CODE       REAL (wp) :: &
-!OLI_CODE          rem_vol, &
-!OLI_CODE          area, &
-!OLI_CODE          vol, &
-!OLI_CODE          tmp, &
-!OLI_CODE          z0   = 0.0_wp,    &   
-!OLI_CODE          epsi10 = 1.0e-11_wp
-!OLI_CODE     
-!OLI_CODE       !----------------------------------------------------------------
-!OLI_CODE       ! hpond is zero if zvolp is zero - have we fully drained? 
-!OLI_CODE       !----------------------------------------------------------------
-!OLI_CODE     
-!OLI_CODE       IF (zvolp < epsi10) THEN
-!OLI_CODE        hpond = z0
-!OLI_CODE        m_index = 0
-!OLI_CODE       ELSE
-!OLI_CODE        
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        ! Calculate the category where water fills up to 
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        
-!OLI_CODE        !----------|
-!OLI_CODE        !          |
-!OLI_CODE        !          |
-!OLI_CODE        !          |----------|                                     -- --
-!OLI_CODE        !__________|__________|_________________________________________ ^
-!OLI_CODE        !          |          |             rem_vol     ^                | Semi-filled
-!OLI_CODE        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer
-!OLI_CODE        !          |          |          |              |
-!OLI_CODE        !          |          |          |              |hpond
-!OLI_CODE        !          |          |          |----------|   |     |-------
-!OLI_CODE        !          |          |          |          |   |     |
-!OLI_CODE        !          |          |          |          |---v-----|
-!OLI_CODE        !          |          | m_index  |          |         |
-!OLI_CODE        !-------------------------------------------------------------
-!OLI_CODE        
-!OLI_CODE        m_index = 0  ! 1:m_index categories have water in them
-!OLI_CODE        DO n = 1, jpl
-!OLI_CODE           IF (zvolp <= cum_max_vol(n)) THEN
-!OLI_CODE              m_index = n
-!OLI_CODE              IF (n == 1) THEN
-!OLI_CODE                 rem_vol = zvolp
-!OLI_CODE              ELSE
-!OLI_CODE                 rem_vol = zvolp - cum_max_vol(n-1)
-!OLI_CODE              END IF
-!OLI_CODE              exit ! to break out of the loop
-!OLI_CODE           END IF
-!OLI_CODE        END DO
-!OLI_CODE        m_index = min(jpl-1, m_index)
-!OLI_CODE        
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        ! semi-filled layer may have m_index different snow in it
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        
-!OLI_CODE        !-----------------------------------------------------------  ^
-!OLI_CODE        !                                                             |  alfan(m_index+1)
-!OLI_CODE        !                                                             |
-!OLI_CODE        !hitl(3)-->                             |----------|          |
-!OLI_CODE        !hitl(2)-->                |------------| * * * * *|          |
-!OLI_CODE        !hitl(1)-->     |----------|* * * * * * |* * * * * |          |
-!OLI_CODE        !hitl(0)-->-------------------------------------------------  |  ^
-!OLI_CODE        !                various snow from lower categories          |  |alfa(m_index)
-!OLI_CODE        
-!OLI_CODE        ! hitl - heights of the snow layers from thinner and current categories
-!OLI_CODE        ! aicetl - area of each snow depth in this layer
-!OLI_CODE        
-!OLI_CODE        hitl(:) = z0
-!OLI_CODE        aicetl(:) = z0
-!OLI_CODE        DO n = 1, m_index
-!OLI_CODE           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
-!OLI_CODE                                  alfan(m_index+1) - alfan(m_index)), z0)
-!OLI_CODE           aicetl(n) = asnon(n)
-!OLI_CODE           
-!OLI_CODE           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
-!OLI_CODE        END DO
-!OLI_CODE        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
-!OLI_CODE        aicetl(m_index+1) = z0
-!OLI_CODE        
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        ! reorder array according to hitl 
-!OLI_CODE        ! snow heights not necessarily in height order
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        
-!OLI_CODE        DO ns = 1, m_index+1
-!OLI_CODE           DO n = 0, m_index - ns + 1
-!OLI_CODE              IF (hitl(n) > hitl(n+1)) THEN ! swap order
-!OLI_CODE                 tmp = hitl(n)
-!OLI_CODE                 hitl(n) = hitl(n+1)
-!OLI_CODE                 hitl(n+1) = tmp
-!OLI_CODE                 tmp = aicetl(n)
-!OLI_CODE                 aicetl(n) = aicetl(n+1)
-!OLI_CODE                 aicetl(n+1) = tmp
-!OLI_CODE              END IF
-!OLI_CODE           END DO
-!OLI_CODE        END DO
-!OLI_CODE        
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        ! divide semi-filled layer into set of sublayers each vertically homogenous
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        
-!OLI_CODE        !hitl(3)----------------------------------------------------------------
-!OLI_CODE        !                                                       | * * * * * * * *  
-!OLI_CODE        !                                                       |* * * * * * * * * 
-!OLI_CODE        !hitl(2)----------------------------------------------------------------
-!OLI_CODE        !                                    | * * * * * * * *  | * * * * * * * *  
-!OLI_CODE        !                                    |* * * * * * * * * |* * * * * * * * * 
-!OLI_CODE        !hitl(1)----------------------------------------------------------------
-!OLI_CODE        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * *  
-!OLI_CODE        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 
-!OLI_CODE        !hitl(0)----------------------------------------------------------------
-!OLI_CODE        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3)            
-!OLI_CODE        
-!OLI_CODE        ! move up over layers incrementing volume
-!OLI_CODE        DO n = 1, m_index+1
-!OLI_CODE           
-!OLI_CODE           area = sum(aicetl(:)) - &                 ! total area of sub-layer
-!OLI_CODE                (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
-!OLI_CODE           
-!OLI_CODE           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area
-!OLI_CODE           
-!OLI_CODE           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within
-!OLI_CODE              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - &
-!OLI_CODE                           alfan(1)
-!OLI_CODE              exit
-!OLI_CODE           ELSE  ! still in sub-layer below the sub-layer with the depth
-!OLI_CODE              rem_vol = rem_vol - vol
-!OLI_CODE           END IF
-!OLI_CODE           
-!OLI_CODE        END DO
-!OLI_CODE        
-!OLI_CODE       END IF
-!OLI_CODE     
-!OLI_CODE    END SUBROUTINE calc_hpond
-!OLI_CODE   
-!OLI_CODE 
-    SUBROUTINE lim_mp_perm(ticen, salin, vicen, perm)
-       !!-------------------------------------------------------------------
-       !!                ***  ROUTINE lim_mp_perm ***
-       !!
-       !! ** Purpose :   Determine the liquid fraction of brine in the ice
-       !!                and its permeability
-       !!-------------------------------------------------------------------
-       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
-          ticen,  & ! energy of melting for each ice layer (J/m2)
-          salin
- 
-       REAL (wp), INTENT(IN) :: &
-          vicen     ! ice volume
-     
-       REAL (wp), INTENT(OUT) :: &
-          perm      ! permeability
- 
-       REAL (wp) ::   &
-          Sbr        ! brine salinity
- 
-       REAL (wp), DIMENSION(nlay_i) ::   &
-          Tin, &    ! ice temperature
-          phi       ! liquid fraction
- 
-       INTEGER :: k
-     
-       REAL (wp) :: &
-          c2    = 2.0_wp
-  
-       !-----------------------------------------------------------------
-       ! Compute ice temperatures from enthalpies using quadratic formula
-       !-----------------------------------------------------------------
- 
-       DO k = 1,nlay_i
-          Tin(k) = ticen(k) 
-       END DO
- 
-       !-----------------------------------------------------------------
-       ! brine salinity and liquid fraction
-       !-----------------------------------------------------------------
- 
-       IF (maxval(Tin-rtt) <= -c2) THEN
- 
-          DO k = 1,nlay_i
-             Sbr = - 1.2_wp                 &
-                   -21.8_wp     * (Tin(k)-rtt)    &
-                   - 0.919_wp   * (Tin(k)-rtt)**2 &
-                   - 0.01878_wp * (Tin(k)-rtt)**3
-             phi(k) = salin(k)/Sbr ! liquid fraction
-          END DO ! k
-        
-       ELSE
- 
-          DO k = 1,nlay_i
-             Sbr = -17.6_wp    * (Tin(k)-rtt)    &
-                   - 0.389_wp  * (Tin(k)-rtt)**2 &
-                   - 0.00362_wp* (Tin(k)-rtt)**3
-             phi(k) = salin(k)/Sbr ! liquid fraction
-          END DO
- 
-       END IF
-     
-       !-----------------------------------------------------------------
-       ! permeability
-       !-----------------------------------------------------------------
- 
-       perm = 3.0e-08_wp * (minval(phi))**3 ! REFERENCE PLEASE (this fucking
-                                            ! bastard of Golden)
-     
-    END SUBROUTINE lim_mp_perm
-!OLI_CODE   
-!OLI_CODE #else
-!OLI_CODE    !!----------------------------------------------------------------------
-!OLI_CODE    !!   Default option         Dummy Module          No LIM-3 sea-ice model
-!OLI_CODE    !!----------------------------------------------------------------------
-!OLI_CODE CONTAINS
-!OLI_CODE    SUBROUTINE lim_mp_init          ! Empty routine
-!OLI_CODE    END SUBROUTINE lim_mp_init   
-!OLI_CODE    SUBROUTINE lim_mp               ! Empty routine
-!OLI_CODE    END SUBROUTINE lim_mp
-!OLI_CODE    SUBROUTINE compute_mp_topo      ! Empty routine
-!OLI_CODE    END SUBROUTINE compute_mp_topo        
-!OLI_CODE    SUBROUTINE pond_area            ! Empty routine
-!OLI_CODE    END SUBROUTINE pond_area   
-!OLI_CODE    SUBROUTINE calc_hpond           ! Empty routine
-!OLI_CODE    END SUBROUTINE calc_hpond   
-!OLI_CODE    SUBROUTINE permeability_phy     ! Empty routine
-!OLI_CODE    END SUBROUTINE permeability_phy   
-!OLI_CODE #endif
-!OLI_CODE    !!======================================================================
-!OLI_CODE END MODULE limmp_topo
-
-
-!OLI_CODE MODULE limmp_topo
-!OLI_CODE    !!======================================================================
-!OLI_CODE    !!                       ***  MODULE limmp_topo ***
-!OLI_CODE    !! LIM-3 sea-ice :  computation of melt ponds' properties
-!OLI_CODE    !!======================================================================
-!OLI_CODE    !! History :  Original code by Daniela Flocco and Adrian Turner
-!OLI_CODE    !!            ! 2012-09 (O. Lecomte) Adaptation for routine inclusion in 
-!OLI_CODE    !!                      NEMO-LIM3.1
-!OLI_CODE    !!            ! 2016-11 (O. Lecomte, C. Rousset, M. Vancoppenolle) 
-!OLI_CODE    !!                      Adaptation for merge with NEMO-LIM3.6
-!OLI_CODE    !!----------------------------------------------------------------------
-!OLI_CODE #if defined key_lim3
-!OLI_CODE    !!----------------------------------------------------------------------
-!OLI_CODE    !!   'key_lim3'                                      LIM-3 sea-ice model
-!OLI_CODE    !!----------------------------------------------------------------------
-!OLI_CODE    !!   lim_mp_init     : melt pond properties initialization
-!OLI_CODE    !!   lim_mp          : melt pond routine caller
-!OLI_CODE    !!   compute_mp_topo : Actual melt pond routine
-!OLI_CODE    !!----------------------------------------------------------------------
-!OLI_CODE    USE ice_oce, ONLY: rdt_ice, tatm_ice
-!OLI_CODE    USE phycst
-!OLI_CODE    USE dom_ice
-!OLI_CODE    USE dom_oce
-!OLI_CODE    USE sbc_oce
-!OLI_CODE    USE sbc_ice
-!OLI_CODE    USE par_ice
-!OLI_CODE    USE par_oce
-!OLI_CODE    USE ice
-!OLI_CODE    USE thd_ice
-!OLI_CODE    USE in_out_manager
-!OLI_CODE    USE lbclnk
-!OLI_CODE    USE lib_mpp
-!OLI_CODE 
-!OLI_CODE    IMPLICIT NONE
-!OLI_CODE    PRIVATE
-!OLI_CODE 
-!OLI_CODE    PUBLIC   lim_mp_init
-!OLI_CODE    PUBLIC   lim_mp
-!OLI_CODE 
-!OLI_CODE CONTAINS
-!OLI_CODE 
-!OLI_CODE    SUBROUTINE lim_mp_init
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE       !!                ***  ROUTINE lim_mp_init  ***
-!OLI_CODE       !!
-!OLI_CODE       !! ** Purpose :   Initialize melt ponds
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE       a_ip_frac(:,:,:)    = 0._wp
-!OLI_CODE       a_ip(:,:,:)         = 0._wp
-!OLI_CODE       h_ip(:,:,:)         = 0._wp
-!OLI_CODE       v_ip(:,:,:)         = 0._wp
-!OLI_CODE       h_il(:,:,:)         = 0._wp
-!OLI_CODE       v_il(:,:,:)         = 0._wp
-!OLI_CODE         
-!OLI_CODE    END SUBROUTINE lim_mp_init
-!OLI_CODE 
-!OLI_CODE 
-!OLI_CODE    SUBROUTINE lim_mp
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE       !!                ***  ROUTINE lim_mp  ***
-!OLI_CODE       !!
-!OLI_CODE       !! ** Purpose :   Compute surface heat flux and call main melt pond
-!OLI_CODE       !!                routine
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE 
-!OLI_CODE       INTEGER  ::   ji, jj, jl   ! dummy loop indices
-!OLI_CODE 
-!OLI_CODE       fsurf(:,:) = 0.e0
-!OLI_CODE       DO jl = 1, jpl
-!OLI_CODE          DO jj = 1, jpj
-!OLI_CODE             DO ji = 1, jpi
-!OLI_CODE                fsurf(ji,jj) = fsurf(ji,jj) + a_i(ji,jj,jl) * &
-!OLI_CODE                      (qns_ice(ji,jj,jl) + (1.0 - izero(ji,jj,jl)) &
-!OLI_CODE                        * qsr_ice(ji,jj,jl))
-!OLI_CODE             END DO
-!OLI_CODE          END DO
-!OLI_CODE       END DO
-!OLI_CODE 
-!OLI_CODE       CALL compute_mp_topo(at_i, a_i,                               &
-!OLI_CODE                          vt_i, v_i, v_s, rhosn_glo, t_i, sz_i, a_ip_frac,  &
-!OLI_CODE                       h_ip, h_il, t_su, tatm_ice, diag_sur_me*rdt_ice, &
-!OLI_CODE                          fsurf, fwoc)
-!OLI_CODE 
-!OLI_CODE       at_ip(:,:) = 0.0
-!OLI_CODE       vt_ip(:,:) = 0.0
-!OLI_CODE       vt_il(:,:) = 0.0
-!OLI_CODE       DO jl = 1, jpl
-!OLI_CODE         DO jj = 1, jpj
-!OLI_CODE            DO ji = 1, jpi
-!OLI_CODE               a_ip(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * &
-!OLI_CODE                                                   a_i(ji,jj,jl))
-!OLI_CODE               v_ip(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * &
-!OLI_CODE                                   a_i(ji,jj,jl) * h_ip(ji,jj,jl))
-!OLI_CODE               v_il(ji,jj,jl) = MAX(0.0_wp, a_ip_frac(ji,jj,jl) * &
-!OLI_CODE                                   a_i(ji,jj,jl) * h_il(ji,jj,jl))
-!OLI_CODE               at_ip(ji,jj)   = at_ip(ji,jj) + a_ip(ji,jj,jl)
-!OLI_CODE               vt_ip(ji,jj)   = vt_ip(ji,jj) + v_ip(ji,jj,jl)
-!OLI_CODE               vt_il(ji,jj)   = vt_il(ji,jj) + v_il(ji,jj,jl)
-!OLI_CODE            END DO
-!OLI_CODE         END DO
-!OLI_CODE       END DO
-!OLI_CODE 
-!OLI_CODE    END SUBROUTINE lim_mp
-!OLI_CODE 
-!OLI_CODE 
-!OLI_CODE    SUBROUTINE compute_mp_topo(aice,      aicen,     &
-!OLI_CODE                               vice,      vicen,     &
-!OLI_CODE                               vsnon,     rhos,      &
-!OLI_CODE                               ticen,     salin,     &
-!OLI_CODE                               a_ip_frac, h_ip,      &
-!OLI_CODE                               h_il,      Tsfc,      &
-!OLI_CODE                               potT,      meltt,     &
-!OLI_CODE                               fsurf,     fwoc)
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE       !!                ***  ROUTINE compute_mp_topo  ***
-!OLI_CODE       !!
-!OLI_CODE       !! ** Purpose :   Compute melt pond evolution based on the ice 
-!OLI_CODE       !!                topography as inferred from the ice thickness 
-!OLI_CODE       !!                distribution.  
-!OLI_CODE       !!
-!OLI_CODE       !! ** Method  :   This code is initially based on Flocco and Feltham 
-!OLI_CODE       !!                (2007) and Flocco et al. (2010). More to come...
-!OLI_CODE       !!
-!OLI_CODE       !! ** Tunable parameters :
-!OLI_CODE       !! 
-!OLI_CODE       !! ** Note :
-!OLI_CODE       !!
-!OLI_CODE       !! ** References
-!OLI_CODE       !!    Flocco, D. and D. L. Feltham, 2007.  A continuum model of melt pond 
-!OLI_CODE       !!    evolution on Arctic sea ice.  J. Geophys. Res. 112, C08016, doi: 
-!OLI_CODE       !!    10.1029/2006JC003836.
-!OLI_CODE       !!    Flocco, D., D. L. Feltham and A. K. Turner, 2010.  Incorporation of
-!OLI_CODE       !!    a physically based melt pond scheme into the sea ice component of a
-!OLI_CODE       !!    climate model.  J. Geophys. Res. 115, C08012, 
-!OLI_CODE       !!    doi: 10.1029/2009JC005568.
-!OLI_CODE       !!    
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj), &
-!OLI_CODE          INTENT(IN) :: &
-!OLI_CODE          aice, &    ! total ice area fraction
-!OLI_CODE          vice       ! total ice volume (m)
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl), &
-!OLI_CODE          INTENT(IN) :: &
-!OLI_CODE          aicen, &   ! ice area fraction, per category
-!OLI_CODE          vsnon, &   ! snow volume, per category (m)
-!OLI_CODE          rhos,  &   ! equivalent snow density, per category (kg/m^3)
-!OLI_CODE          vicen      ! ice volume, per category (m)
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), &
-!OLI_CODE          INTENT(IN) :: &
-!OLI_CODE          ticen, &   ! ice enthalpy, per category
-!OLI_CODE          salin
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl), &
-!OLI_CODE          INTENT(INOUT) :: &
-!OLI_CODE          a_ip_frac , &   ! pond area fraction of ice, per ice category
-!OLI_CODE          h_ip , &   ! pond depth, per ice category (m)
-!OLI_CODE          h_il       ! Refrozen ice lid thickness, per ice category (m)
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj), &
-!OLI_CODE          INTENT(IN) :: &
-!OLI_CODE          potT,  &   ! air potential temperature
-!OLI_CODE          meltt, &   ! total surface meltwater flux
-!OLI_CODE          fsurf      ! thermodynamic heat flux at ice/snow surface (W/m^2)
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj), &
-!OLI_CODE          INTENT(INOUT) :: &
-!OLI_CODE          fwoc      ! fresh water flux to the ocean (from draining and other pond volume adjustments)
-!OLI_CODE                    ! (m)
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl), &
-!OLI_CODE          INTENT(IN) :: &
-!OLI_CODE          Tsfc       ! snow/sea ice surface temperature
-!OLI_CODE 
-!OLI_CODE       ! local variables
-!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
-!OLI_CODE          zTsfcn, & ! ice/snow surface temperature (C)
-!OLI_CODE          zvolpn, & ! pond volume per unit area, per category (m)
-!OLI_CODE          zvuin     ! water-equivalent volume of ice lid on melt pond ('upper ice', m)
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
-!OLI_CODE          zapondn,& ! pond area fraction, per category
-!OLI_CODE          zhpondn   ! pond depth, per category (m)
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION (jpi,jpj) :: &
-!OLI_CODE          zvolp       ! total volume of pond, per unit area of pond (m)
-!OLI_CODE 
-!OLI_CODE       REAL (wp) :: &
-!OLI_CODE          zhi,    & ! ice thickness (m)
-!OLI_CODE          zdHui,  & ! change in thickness of ice lid (m)
-!OLI_CODE          zomega, & ! conduction
-!OLI_CODE          zdTice, & ! temperature difference across ice lid (C)
-!OLI_CODE          zdvice, & ! change in ice volume (m)
-!OLI_CODE          zTavg,  & ! mean surface temperature across categories (C)
-!OLI_CODE          zTp,    & ! pond freezing temperature (C)
-!OLI_CODE          zdvn      ! change in melt pond volume for fresh water budget
-!OLI_CODE       INTEGER, DIMENSION (jpi*jpj) :: &
-!OLI_CODE          indxi, indxj    ! compressed indices for cells with ice melting
-!OLI_CODE 
-!OLI_CODE       INTEGER :: n,k,i,j,ij,icells,indxij ! loop indices
-!OLI_CODE 
-!OLI_CODE       INTEGER, DIMENSION (jpl) :: &
-!OLI_CODE          kcells          ! cells where ice lid combines with vice
-!OLI_CODE 
-!OLI_CODE       INTEGER, DIMENSION (jpi*jpj,jpl) :: &
-!OLI_CODE          indxii, indxjj  ! i,j indices for kcells loop
-!OLI_CODE 
-!OLI_CODE       REAL (wp), parameter :: &
-!OLI_CODE          zhicemin  = 0.1_wp , & ! minimum ice thickness with ponds (m) 
-!OLI_CODE          zTd       = 0.15_wp, & ! temperature difference for freeze-up (C)
-!OLI_CODE          zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3)
-!OLI_CODE          zmin_volp = 1.e-4_wp, & ! minimum pond volume (m)
-!OLI_CODE          z0       = 0._wp,    & 
-!OLI_CODE          zTimelt   = 0._wp,    &
-!OLI_CODE          z01      = 0.01_wp,  &
-!OLI_CODE          z25      = 0.25_wp,  &
-!OLI_CODE          z5       = 0.5_wp,   &
-!OLI_CODE          epsi10     = 1.0e-11_wp
-!OLI_CODE       !---------------------------------------------------------------
-!OLI_CODE       ! initialize
-!OLI_CODE       !---------------------------------------------------------------
-!OLI_CODE 
-!OLI_CODE       DO j = 1, jpj
-!OLI_CODE          DO i = 1, jpi
-!OLI_CODE             zvolp(i,j) = z0
-!OLI_CODE          END DO
-!OLI_CODE       END DO
-!OLI_CODE       DO n = 1, jpl
-!OLI_CODE          DO j = 1, jpj
-!OLI_CODE             DO i = 1, jpi
-!OLI_CODE                ! load tracers
-!OLI_CODE                zvolp(i,j) = zvolp(i,j) + h_ip(i,j,n) &
-!OLI_CODE                                      * a_ip_frac(i,j,n) * aicen(i,j,n)
-!OLI_CODE                zTsfcn(i,j,n) = Tsfc(i,j,n) - rtt ! convert in Celsius - Oli
-!OLI_CODE                zvuin (i,j,n) = h_il(i,j,n) &
-!OLI_CODE                             * a_ip_frac(i,j,n) * aicen(i,j,n)
-!OLI_CODE 
-!OLI_CODE                zhpondn(i,j,n) = z0     ! pond depth, per category
-!OLI_CODE                zapondn(i,j,n) = z0     ! pond area,  per category
-!OLI_CODE             END DO
-!OLI_CODE          END DO
-!OLI_CODE          indxii(:,n) = 0
-!OLI_CODE          indxjj(:,n) = 0
-!OLI_CODE          kcells  (n) = 0
-!OLI_CODE       END DO
-!OLI_CODE 
-!OLI_CODE       ! The freezing temperature for meltponds is assumed slightly below 0C,
-!OLI_CODE       ! as if meltponds had a little salt in them.  The salt budget is not
-!OLI_CODE       ! altered for meltponds, but if it were then an actual pond freezing 
-!OLI_CODE       ! temperature could be computed.
-!OLI_CODE 
-!OLI_CODE       zTp = zTimelt - zTd
-!OLI_CODE 
-!OLI_CODE       !-----------------------------------------------------------------
-!OLI_CODE       ! Identify grid cells with ponds
-!OLI_CODE       !-----------------------------------------------------------------
-!OLI_CODE 
-!OLI_CODE       icells = 0
-!OLI_CODE       DO j = 1, jpj
-!OLI_CODE       DO i = 1, jpi
-!OLI_CODE          zhi = z0
-!OLI_CODE          IF (aice(i,j) > epsi10) zhi = vice(i,j)/aice(i,j)
-!OLI_CODE          IF ( aice(i,j) > z01 .and. zhi > zhicemin .and. &
-!OLI_CODE             zvolp(i,j) > zmin_volp*aice(i,j)) THEN
-!OLI_CODE             icells = icells + 1
-!OLI_CODE             indxi(icells) = i
-!OLI_CODE             indxj(icells) = j
-!OLI_CODE          ELSE  ! remove ponds on thin ice
-!OLI_CODE             !fpond(i,j) = fpond(i,j) - zvolp(i,j)
-!OLI_CODE             zvolpn(i,j,:) = z0
-!OLI_CODE             zvuin (i,j,:) = z0
-!OLI_CODE             zvolp (i,j) = z0
-!OLI_CODE          END IF
-!OLI_CODE       END DO                     ! i
-!OLI_CODE       END DO                     ! j
-!OLI_CODE 
-!OLI_CODE       DO ij = 1, icells
-!OLI_CODE          i = indxi(ij)
-!OLI_CODE          j = indxj(ij)
-!OLI_CODE 
-!OLI_CODE          !--------------------------------------------------------------
-!OLI_CODE          ! calculate pond area and depth
-!OLI_CODE          !--------------------------------------------------------------
-!OLI_CODE          CALL pond_area(aice(i,j),vice(i,j),rhos(i,j,:), &
-!OLI_CODE                    aicen(i,j,:), vicen(i,j,:), vsnon(i,j,:), &
-!OLI_CODE                    ticen(i,j,:,:), salin(i,j,:,:), &
-!OLI_CODE                    zvolpn(i,j,:), zvolp(i,j), &
-!OLI_CODE                    zapondn(i,j,:),zhpondn(i,j,:), zdvn)
-!OLI_CODE 
-!OLI_CODE          fwoc(i,j) = fwoc(i,j) + zdvn ! -> Goes to fresh water budget
-!OLI_CODE 
-!OLI_CODE          ! mean surface temperature
-!OLI_CODE          zTavg = z0
-!OLI_CODE          DO n = 1, jpl
-!OLI_CODE             zTavg = zTavg + zTsfcn(i,j,n)*aicen(i,j,n)
-!OLI_CODE          END DO
-!OLI_CODE          zTavg = zTavg / aice(i,j)
-!OLI_CODE 
-!OLI_CODE          DO n = 1, jpl-1
-!OLI_CODE                        
-!OLI_CODE             IF (zvuin(i,j,n) > epsi10) THEN
-!OLI_CODE 
-!OLI_CODE          !----------------------------------------------------------------
-!OLI_CODE          ! melting: floating upper ice layer melts in whole or part
-!OLI_CODE          !----------------------------------------------------------------
-!OLI_CODE !               IF (zTsfcn(i,j,n) > zTp) THEN
-!OLI_CODE                IF (zTavg > zTp) THEN
-!OLI_CODE 
-!OLI_CODE                   zdvice = min(meltt(i,j)*zapondn(i,j,n), zvuin(i,j,n))
-!OLI_CODE                   IF (zdvice > epsi10) THEN
-!OLI_CODE                      zvuin (i,j,n) = zvuin (i,j,n) - zdvice
-!OLI_CODE                      zvolpn(i,j,n) = zvolpn(i,j,n) + zdvice
-!OLI_CODE                      zvolp (i,j)   = zvolp (i,j)   + zdvice
-!OLI_CODE                      !fwoc(i,j)   = fwoc(i,j)   + zdvice
-!OLI_CODE                        
-!OLI_CODE                      IF (zvuin(i,j,n) < epsi10 .and. zvolpn(i,j,n) > puny) THEN
-!OLI_CODE                         ! ice lid melted and category is pond covered
-!OLI_CODE                         zvolpn(i,j,n) = zvolpn(i,j,n) + zvuin(i,j,n)
-!OLI_CODE                         !fwoc(i,j)   = fwoc(i,j)   + zvuin(i,j,n)
-!OLI_CODE                         zvuin(i,j,n)  = z0
-!OLI_CODE                      END IF
-!OLI_CODE                      zhpondn(i,j,n) = zvolpn(i,j,n) / zapondn(i,j,n)
-!OLI_CODE                   END IF
-!OLI_CODE 
-!OLI_CODE          !----------------------------------------------------------------
-!OLI_CODE          ! freezing: existing upper ice layer grows
-!OLI_CODE          !----------------------------------------------------------------
-!OLI_CODE                ELSE IF (zvolpn(i,j,n) > epsi10) THEN ! zTavg <= zTp
-!OLI_CODE 
-!OLI_CODE                 ! dIFferential growth of base of surface floating ice layer
-!OLI_CODE                   zdTice = max(-zTavg, z0) ! > 0
-!OLI_CODE                   zomega = rcdic*zdTice * zr1_rlfus 
-!OLI_CODE                   zdHui = sqrt(zomega*rdt_ice + z25*(zvuin(i,j,n)/  &
-!OLI_CODE                         aicen(i,j,n))**2)- z5 * zvuin(i,j,n)/aicen(i,j,n)
-!OLI_CODE 
-!OLI_CODE                   zdvice = min(zdHui*zapondn(i,j,n), zvolpn(i,j,n))   
-!OLI_CODE                   IF (zdvice > epsi10) THEN
-!OLI_CODE                      zvuin (i,j,n) = zvuin (i,j,n) + zdvice
-!OLI_CODE                      zvolpn(i,j,n) = zvolpn(i,j,n) - zdvice
-!OLI_CODE                      zvolp (i,j)   = zvolp (i,j)   - zdvice
-!OLI_CODE                      !fwoc(i,j)   = fwoc(i,j)   - zdvice
-!OLI_CODE                      zhpondn(i,j,n) = zvolpn(i,j,n) / zapondn(i,j,n)
-!OLI_CODE                   END IF
-!OLI_CODE 
-!OLI_CODE                END IF ! zTavg
-!OLI_CODE 
-!OLI_CODE          !----------------------------------------------------------------
-!OLI_CODE          ! freezing: upper ice layer begins to form
-!OLI_CODE          ! note: albedo does not change
-!OLI_CODE          !----------------------------------------------------------------
-!OLI_CODE             ELSE ! zvuin < epsi10
-!OLI_CODE                     
-!OLI_CODE                ! thickness of newly formed ice
-!OLI_CODE                ! the surface temperature of a meltpond is the same as that
-!OLI_CODE                ! of the ice underneath (0C), and the thermodynamic surface 
-!OLI_CODE                ! flux is the same
-!OLI_CODE                zdHui = max(-fsurf(i,j)*rdt_ice*zr1_rlfus, z0)
-!OLI_CODE                zdvice = min(zdHui*zapondn(i,j,n), zvolpn(i,j,n))  
-!OLI_CODE                IF (zdvice > epsi10) THEN
-!OLI_CODE                   zvuin (i,j,n) = zdvice
-!OLI_CODE                   zvolpn(i,j,n) = zvolpn(i,j,n) - zdvice
-!OLI_CODE                   zvolp (i,j)   = zvolp (i,j)   - zdvice
-!OLI_CODE                   !fwoc(i,j)   = fwoc(i,j)   - zdvice
-!OLI_CODE                   zhpondn(i,j,n)= zvolpn(i,j,n) / zapondn(i,j,n)
-!OLI_CODE                END IF
-!OLI_CODE                     
-!OLI_CODE             END IF  ! zvuin
-!OLI_CODE 
-!OLI_CODE          END DO ! jpl
-!OLI_CODE 
-!OLI_CODE       END DO ! ij
-!OLI_CODE 
-!OLI_CODE       !---------------------------------------------------------------
-!OLI_CODE       ! remove ice lid if there is no liquid pond
-!OLI_CODE       ! zvuin may be nonzero on category jpl due to dynamics
-!OLI_CODE       !---------------------------------------------------------------
-!OLI_CODE 
-!OLI_CODE       DO j = 1, jpj
-!OLI_CODE       DO i = 1, jpi
-!OLI_CODE          DO n = 1, jpl
-!OLI_CODE             IF (aicen(i,j,n) > epsi10 .and. zvolpn(i,j,n) < puny &
-!OLI_CODE                                     .and. zvuin (i,j,n) > epsi10) THEN
-!OLI_CODE                kcells(n) = kcells(n) + 1
-!OLI_CODE                indxij    = kcells(n)
-!OLI_CODE                indxii(indxij,n) = i
-!OLI_CODE                indxjj(indxij,n) = j
-!OLI_CODE             END IF
-!OLI_CODE          END DO
-!OLI_CODE       END DO                     ! i
-!OLI_CODE       END DO                     ! j
-!OLI_CODE 
-!OLI_CODE       DO n = 1, jpl
-!OLI_CODE 
-!OLI_CODE          IF (kcells(n) > 0) THEN
-!OLI_CODE          DO ij = 1, kcells(n)
-!OLI_CODE             i = indxii(ij,n)
-!OLI_CODE             j = indxjj(ij,n)
-!OLI_CODE             fwoc(i,j) = fwoc(i,j) + rhoic/rauw * zvuin(i,j,n) ! Completely refrozen lid goes into ocean (to be changed)
-!OLI_CODE             zvuin(i,j,n) = z0 
-!OLI_CODE          END DO    ! ij
-!OLI_CODE          END IF
-!OLI_CODE 
-!OLI_CODE          ! reload tracers
-!OLI_CODE          DO j = 1, jpj
-!OLI_CODE             DO i = 1, jpi
-!OLI_CODE                IF (zapondn(i,j,n) > epsi10) THEN
-!OLI_CODE                   h_il(i,j,n) = zvuin(i,j,n) / zapondn(i,j,n)
-!OLI_CODE                ELSE
-!OLI_CODE                   zvuin(i,j,n) = z0
-!OLI_CODE                   h_il(i,j,n) = z0
-!OLI_CODE                END IF
-!OLI_CODE                IF (aicen(i,j,n) > epsi10) THEN
-!OLI_CODE                   a_ip_frac(i,j,n) = zapondn(i,j,n) / aicen(i,j,n) * &
-!OLI_CODE                         (1.0_wp - MAX(z0, SIGN(1.0_wp, -zvolpn(i,j,n))))
-!OLI_CODE                   h_ip(i,j,n) = zhpondn(i,j,n)
-!OLI_CODE                ELSE
-!OLI_CODE                   a_ip_frac(i,j,n) = z0
-!OLI_CODE                   h_ip(i,j,n) = z0
-!OLI_CODE                   h_il(i,j,n) = z0
-!OLI_CODE                END IF
-!OLI_CODE             END DO ! i
-!OLI_CODE          END DO    ! j
-!OLI_CODE 
-!OLI_CODE       END DO       ! n
-!OLI_CODE 
-!OLI_CODE    END SUBROUTINE compute_mp_topo
-!OLI_CODE 
-!OLI_CODE 
-!OLI_CODE    SUBROUTINE pond_area(aice,vice,rhos,             &
-!OLI_CODE                         aicen, vicen, vsnon, ticen, &
-!OLI_CODE                         salin, zvolpn, zvolp,         &
-!OLI_CODE                         zapondn,zhpondn,dvolp)
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE       !!                ***  ROUTINE pond_area  ***
-!OLI_CODE       !!
-!OLI_CODE       !! ** Purpose :   Compute melt pond area, depth and melting rates
-!OLI_CODE       !!------------------------------------------------------------------
-!OLI_CODE       REAL (wp), INTENT(IN) :: &
-!OLI_CODE          aice,vice
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
-!OLI_CODE          aicen, vicen, vsnon, rhos
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: &
-!OLI_CODE          ticen, salin
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: &
-!OLI_CODE          zvolpn
-!OLI_CODE 
-!OLI_CODE       REAL (wp), INTENT(INOUT) :: &
-!OLI_CODE          zvolp, dvolp
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(OUT) :: &
-!OLI_CODE          zapondn, zhpondn
-!OLI_CODE 
-!OLI_CODE       INTEGER :: &
-!OLI_CODE          n, ns,   &
-!OLI_CODE          m_index, &
-!OLI_CODE          permflag
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION(jpl) :: &
-!OLI_CODE          hicen, &
-!OLI_CODE          hsnon, &
-!OLI_CODE          asnon, &
-!OLI_CODE          alfan, &
-!OLI_CODE          betan, &
-!OLI_CODE          cum_max_vol, &
-!OLI_CODE          reduced_aicen        
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION(0:jpl) :: &
-!OLI_CODE          cum_max_vol_tmp
-!OLI_CODE 
-!OLI_CODE       REAL (wp) :: &
-!OLI_CODE          hpond, &
-!OLI_CODE          drain, &
-!OLI_CODE          floe_weight, &
-!OLI_CODE          pressure_head, &
-!OLI_CODE          hsl_rel, &
-!OLI_CODE          deltah, &
-!OLI_CODE          perm, &
-!OLI_CODE          apond, &
-!OLI_CODE          msno
-!OLI_CODE 
-!OLI_CODE       REAL (wp), parameter :: & 
-!OLI_CODE          viscosity = 1.79e-3_wp, &  ! kinematic water viscosity in kg/m/s
-!OLI_CODE          z0        = 0.0_wp    , &
-!OLI_CODE          c1        = 1.0_wp    , &
-!OLI_CODE          p4        = 0.4_wp    , &
-!OLI_CODE          p6        = 0.6_wp    , &
-!OLI_CODE          epsi10      = 1.0e-11_wp
-!OLI_CODE          
-!OLI_CODE      !-----------|
-!OLI_CODE      !           |
-!OLI_CODE      !           |-----------|
-!OLI_CODE      !___________|___________|______________________________________sea-level
-!OLI_CODE      !           |           |
-!OLI_CODE      !           |           |---^--------|
-!OLI_CODE      !           |           |   |        |
-!OLI_CODE      !           |           |   |        |-----------|              |-------
-!OLI_CODE      !           |           |   |alfan(n)|           |              |
-!OLI_CODE      !           |           |   |        |           |--------------|
-!OLI_CODE      !           |           |   |        |           |              |
-!OLI_CODE      !---------------------------v-------------------------------------------
-!OLI_CODE      !           |           |   ^        |           |              |
-!OLI_CODE      !           |           |   |        |           |--------------|
-!OLI_CODE      !           |           |   |betan(n)|           |              |
-!OLI_CODE      !           |           |   |        |-----------|              |-------
-!OLI_CODE      !           |           |   |        |
-!OLI_CODE      !           |           |---v------- |
-!OLI_CODE      !           |           |
-!OLI_CODE      !           |-----------|
-!OLI_CODE      !           |
-!OLI_CODE      !-----------|
-!OLI_CODE     
-!OLI_CODE       !-------------------------------------------------------------------
-!OLI_CODE       ! initialize
-!OLI_CODE       !-------------------------------------------------------------------
-!OLI_CODE 
-!OLI_CODE       DO n = 1, jpl
-!OLI_CODE 
-!OLI_CODE          zapondn(n) = z0
-!OLI_CODE          zhpondn(n) = z0
-!OLI_CODE 
-!OLI_CODE          IF (aicen(n) < epsi10)  THEN
-!OLI_CODE             hicen(n) =  z0 
-!OLI_CODE             hsnon(n) = z0
-!OLI_CODE             reduced_aicen(n) = z0
-!OLI_CODE          ELSE
-!OLI_CODE             hicen(n) = vicen(n) / aicen(n)
-!OLI_CODE             hsnon(n) = vsnon(n) / aicen(n)
-!OLI_CODE             reduced_aicen(n) = c1 ! n=jpl
-!OLI_CODE             IF (n < jpl) reduced_aicen(n) = aicen(n) &
-!OLI_CODE                                  * (-0.024_wp*hicen(n) + 0.832_wp) 
-!OLI_CODE             asnon(n) = reduced_aicen(n) 
-!OLI_CODE          END IF
-!OLI_CODE 
-!OLI_CODE ! This choice for alfa and beta ignores hydrostatic equilibium of categories.
-!OLI_CODE ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
-!OLI_CODE ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
-!OLI_CODE ! categories.  alfa and beta partition the ITD - they are areas not thicknesses!
-!OLI_CODE ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
-!OLI_CODE ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 
-!OLI_CODE ! alfan = 60% of the ice volume) in each category lies above the reference line, 
-!OLI_CODE ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
-!OLI_CODE 
-!OLI_CODE          alfan(n) = p6 * hicen(n)
-!OLI_CODE          betan(n) = p4 * hicen(n)
-!OLI_CODE        
-!OLI_CODE          cum_max_vol(n)     = z0
-!OLI_CODE          cum_max_vol_tmp(n) = z0
-!OLI_CODE     
-!OLI_CODE       END DO ! jpl
-!OLI_CODE 
-!OLI_CODE       cum_max_vol_tmp(0) = z0
-!OLI_CODE       drain = z0
-!OLI_CODE       dvolp = z0
-!OLI_CODE     
-!OLI_CODE       !--------------------------------------------------------------------------
-!OLI_CODE       ! the maximum amount of water that can be contained up to each ice category
-!OLI_CODE       !--------------------------------------------------------------------------
-!OLI_CODE     
-!OLI_CODE       DO n = 1, jpl-1 ! last category can not hold any volume
-!OLI_CODE 
-!OLI_CODE          IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN
-!OLI_CODE 
-!OLI_CODE             ! total volume in level including snow
-!OLI_CODE             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + &
-!OLI_CODE                (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n)) 
-!OLI_CODE 
-!OLI_CODE 
-!OLI_CODE             ! subtract snow solid volumes from lower categories in current level
-!OLI_CODE             DO ns = 1, n 
-!OLI_CODE                cum_max_vol_tmp(n) = cum_max_vol_tmp(n) &
-!OLI_CODE                   - rhos(ns)/rauw  * &    ! fraction of snow that is occupied by solid ??rauw
-!OLI_CODE                     asnon(ns)  * &    ! area of snow from that category
-!OLI_CODE                     max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)- &
-!OLI_CODE                                   alfan(n)), z0)  
-!OLI_CODE                                       ! thickness of snow from ns layer in n layer
-!OLI_CODE             END DO
-!OLI_CODE 
-!OLI_CODE          ELSE ! assume higher categories unoccupied
-!OLI_CODE             cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1)
-!OLI_CODE          END IF
-!OLI_CODE          !IF (cum_max_vol_tmp(n) < z0) THEN
-!OLI_CODE          !   call abort_ice('negative melt pond volume')
-!OLI_CODE          !END IF
-!OLI_CODE       END DO
-!OLI_CODE       cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume
-!OLI_CODE       cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl)
-!OLI_CODE     
-!OLI_CODE       !----------------------------------------------------------------
-!OLI_CODE       ! is there more meltwater than can be held in the floe?
-!OLI_CODE       !----------------------------------------------------------------
-!OLI_CODE       IF (zvolp >= cum_max_vol(jpl)) THEN
-!OLI_CODE          drain = zvolp - cum_max_vol(jpl) + epsi10
-!OLI_CODE          zvolp = zvolp - drain
-!OLI_CODE          dvolp = drain
-!OLI_CODE          IF (zvolp < epsi10) THEN
-!OLI_CODE             dvolp = dvolp + zvolp
-!OLI_CODE             zvolp = z0
-!OLI_CODE          END IF
-!OLI_CODE       END IF
-!OLI_CODE     
-!OLI_CODE       ! height and area corresponding to the remaining volume
-!OLI_CODE 
-!OLI_CODE       call calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, &
-!OLI_CODE            zvolp, cum_max_vol, hpond, m_index)
-!OLI_CODE     
-!OLI_CODE       DO n=1, m_index
-!OLI_CODE          zhpondn(n) = hpond - alfan(n) + alfan(1)
-!OLI_CODE          zapondn(n) = reduced_aicen(n) 
-!OLI_CODE       END DO
-!OLI_CODE       apond = sum(zapondn(1:m_index))
-!OLI_CODE     
-!OLI_CODE       !------------------------------------------------------------------------
-!OLI_CODE       ! drainage due to ice permeability - Darcy's law
-!OLI_CODE       !------------------------------------------------------------------------
-!OLI_CODE     
-!OLI_CODE       ! sea water level 
-!OLI_CODE       msno = z0
-!OLI_CODE       DO n=1,jpl
-!OLI_CODE         msno = msno + vsnon(n) * rhos(n)
-!OLI_CODE       END DO
-!OLI_CODE       floe_weight = (msno + rhoic*vice + rau0*zvolp) / aice
-!OLI_CODE       hsl_rel = floe_weight / rau0 &
-!OLI_CODE               - ((sum(betan(:)*aicen(:))/aice) + alfan(1))
-!OLI_CODE     
-!OLI_CODE       deltah = hpond - hsl_rel
-!OLI_CODE       pressure_head = grav * rau0 * max(deltah, z0)
-!OLI_CODE 
-!OLI_CODE       ! drain IF ice is permeable    
-!OLI_CODE       permflag = 0
-!OLI_CODE       IF (pressure_head > z0) THEN
-!OLI_CODE       DO n = 1, jpl-1
-!OLI_CODE          IF (hicen(n) /= z0) THEN
-!OLI_CODE             CALL permeability_phi(ticen(:,n), salin(:,n), vicen(n), perm)
-!OLI_CODE             IF (perm > z0) permflag = 1
-!OLI_CODE             drain = perm*zapondn(n)*pressure_head*rdt_ice / &
-!OLI_CODE                                      (viscosity*hicen(n))
-!OLI_CODE             dvolp = dvolp + min(drain, zvolp)
-!OLI_CODE             zvolp = max(zvolp - drain, z0)
-!OLI_CODE             IF (zvolp < epsi10) THEN
-!OLI_CODE                dvolp = dvolp + zvolp
-!OLI_CODE                zvolp = z0
-!OLI_CODE             END IF
-!OLI_CODE          END IF
-!OLI_CODE       END DO
-!OLI_CODE  
-!OLI_CODE       ! adjust melt pond DIMENSIONs
-!OLI_CODE       IF (permflag > 0) THEN
-!OLI_CODE          ! recompute pond depth    
-!OLI_CODE          CALL calc_hpond(reduced_aicen, asnon, hsnon, rhos, alfan, &
-!OLI_CODE                          zvolp, cum_max_vol, hpond, m_index)
-!OLI_CODE          DO n=1, m_index
-!OLI_CODE             zhpondn(n) = hpond - alfan(n) + alfan(1)
-!OLI_CODE             zapondn(n) = reduced_aicen(n) 
-!OLI_CODE          END DO
-!OLI_CODE          apond = sum(zapondn(1:m_index))
-!OLI_CODE       END IF
-!OLI_CODE       END IF ! pressure_head
-!OLI_CODE 
-!OLI_CODE       !------------------------------------------------------------------------
-!OLI_CODE       ! total melt pond volume in category DOes not include snow volume
-!OLI_CODE       ! snow in melt ponds is not melted
-!OLI_CODE       !------------------------------------------------------------------------
-!OLI_CODE 
-!OLI_CODE       ! Calculate pond volume for lower categories
-!OLI_CODE       DO n=1,m_index-1
-!OLI_CODE          zvolpn(n) = zapondn(n) * zhpondn(n) &
-!OLI_CODE                   - (rhos(n)/rauw) * asnon(n) * min(hsnon(n), zhpondn(n))! 
-!OLI_CODE       END DO
-!OLI_CODE 
-!OLI_CODE       ! Calculate pond volume for highest category = remaining pond volume
-!OLI_CODE       IF (m_index == 1) zvolpn(m_index) = zvolp
-!OLI_CODE       IF (m_index > 1) THEN
-!OLI_CODE         IF (zvolp > sum(zvolpn(1:m_index-1))) THEN
-!OLI_CODE           zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1))
-!OLI_CODE         ELSE
-!OLI_CODE           zvolpn(m_index) = z0
-!OLI_CODE           zhpondn(m_index) = z0
-!OLI_CODE           zapondn(m_index) = z0
-!OLI_CODE           ! If remaining pond volume is negative reduce pond volume of 
-!OLI_CODE           ! lower category
-!OLI_CODE           IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) & 
-!OLI_CODE             zvolpn(m_index-1) = zvolpn(m_index-1)-sum(zvolpn(1:m_index-1))&
-!OLI_CODE                                + zvolp
-!OLI_CODE         END IF
-!OLI_CODE       END IF
-!OLI_CODE 
-!OLI_CODE       DO n=1,m_index
-!OLI_CODE          IF (zapondn(n) > epsi10) THEN
-!OLI_CODE              zhpondn(n) = zvolpn(n) / zapondn(n)
-!OLI_CODE          ELSE
-!OLI_CODE             dvolp = dvolp + zvolpn(n)
-!OLI_CODE             zhpondn(n) = z0
-!OLI_CODE             zvolpn(n) = z0
-!OLI_CODE             zapondn(n) = z0
-!OLI_CODE          end IF
-!OLI_CODE       END DO
-!OLI_CODE       DO n = m_index+1, jpl
-!OLI_CODE          zhpondn(n) = z0
-!OLI_CODE          zapondn(n) = z0
-!OLI_CODE          zvolpn (n) = z0
-!OLI_CODE       END DO
-!OLI_CODE 
-!OLI_CODE    END SUBROUTINE pond_area
-!OLI_CODE   
-!OLI_CODE 
-!OLI_CODE    SUBROUTINE calc_hpond(aicen, asnon, hsnon, rhos, alfan, &
-!OLI_CODE                          zvolp, cum_max_vol,                &
-!OLI_CODE                          hpond, m_index)
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE       !!                ***  ROUTINE calc_hpond  ***
-!OLI_CODE       !!
-!OLI_CODE       !! ** Purpose :   Compute melt pond depth
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE     
-!OLI_CODE       REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
-!OLI_CODE          aicen, &
-!OLI_CODE          asnon, &
-!OLI_CODE          hsnon, &
-!OLI_CODE          rhos,  &
-!OLI_CODE          alfan, &
-!OLI_CODE          cum_max_vol
-!OLI_CODE     
-!OLI_CODE       REAL (wp), INTENT(IN) :: &
-!OLI_CODE          zvolp
-!OLI_CODE     
-!OLI_CODE       REAL (wp), INTENT(OUT) :: &
-!OLI_CODE          hpond
-!OLI_CODE     
-!OLI_CODE       INTEGER, INTENT(OUT) :: &
-!OLI_CODE          m_index
-!OLI_CODE     
-!OLI_CODE       INTEGER :: n, ns
-!OLI_CODE     
-!OLI_CODE       REAL (wp), DIMENSION(0:jpl+1) :: &
-!OLI_CODE          hitl, &
-!OLI_CODE          aicetl
-!OLI_CODE     
-!OLI_CODE       REAL (wp) :: &
-!OLI_CODE          rem_vol, &
-!OLI_CODE          area, &
-!OLI_CODE          vol, &
-!OLI_CODE          tmp, &
-!OLI_CODE          z0   = 0.0_wp,    &   
-!OLI_CODE          epsi10 = 1.0e-11_wp
-!OLI_CODE     
-!OLI_CODE       !----------------------------------------------------------------
-!OLI_CODE       ! hpond is zero if zvolp is zero - have we fully drained? 
-!OLI_CODE       !----------------------------------------------------------------
-!OLI_CODE     
-!OLI_CODE       IF (zvolp < epsi10) THEN
-!OLI_CODE        hpond = z0
-!OLI_CODE        m_index = 0
-!OLI_CODE       ELSE
-!OLI_CODE        
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        ! Calculate the category where water fills up to 
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        
-!OLI_CODE        !----------|
-!OLI_CODE        !          |
-!OLI_CODE        !          |
-!OLI_CODE        !          |----------|                                     -- --
-!OLI_CODE        !__________|__________|_________________________________________ ^
-!OLI_CODE        !          |          |             rem_vol     ^                | Semi-filled
-!OLI_CODE        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer
-!OLI_CODE        !          |          |          |              |
-!OLI_CODE        !          |          |          |              |hpond
-!OLI_CODE        !          |          |          |----------|   |     |-------
-!OLI_CODE        !          |          |          |          |   |     |
-!OLI_CODE        !          |          |          |          |---v-----|
-!OLI_CODE        !          |          | m_index  |          |         |
-!OLI_CODE        !-------------------------------------------------------------
-!OLI_CODE        
-!OLI_CODE        m_index = 0  ! 1:m_index categories have water in them
-!OLI_CODE        DO n = 1, jpl
-!OLI_CODE           IF (zvolp <= cum_max_vol(n)) THEN
-!OLI_CODE              m_index = n
-!OLI_CODE              IF (n == 1) THEN
-!OLI_CODE                 rem_vol = zvolp
-!OLI_CODE              ELSE
-!OLI_CODE                 rem_vol = zvolp - cum_max_vol(n-1)
-!OLI_CODE              END IF
-!OLI_CODE              exit ! to break out of the loop
-!OLI_CODE           END IF
-!OLI_CODE        END DO
-!OLI_CODE        m_index = min(jpl-1, m_index)
-!OLI_CODE        
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        ! semi-filled layer may have m_index different snow in it
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        
-!OLI_CODE        !-----------------------------------------------------------  ^
-!OLI_CODE        !                                                             |  alfan(m_index+1)
-!OLI_CODE        !                                                             |
-!OLI_CODE        !hitl(3)-->                             |----------|          |
-!OLI_CODE        !hitl(2)-->                |------------| * * * * *|          |
-!OLI_CODE        !hitl(1)-->     |----------|* * * * * * |* * * * * |          |
-!OLI_CODE        !hitl(0)-->-------------------------------------------------  |  ^
-!OLI_CODE        !                various snow from lower categories          |  |alfa(m_index)
-!OLI_CODE        
-!OLI_CODE        ! hitl - heights of the snow layers from thinner and current categories
-!OLI_CODE        ! aicetl - area of each snow depth in this layer
-!OLI_CODE        
-!OLI_CODE        hitl(:) = z0
-!OLI_CODE        aicetl(:) = z0
-!OLI_CODE        DO n = 1, m_index
-!OLI_CODE           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
-!OLI_CODE                                  alfan(m_index+1) - alfan(m_index)), z0)
-!OLI_CODE           aicetl(n) = asnon(n)
-!OLI_CODE           
-!OLI_CODE           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
-!OLI_CODE        END DO
-!OLI_CODE        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
-!OLI_CODE        aicetl(m_index+1) = z0
-!OLI_CODE        
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        ! reorder array according to hitl 
-!OLI_CODE        ! snow heights not necessarily in height order
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        
-!OLI_CODE        DO ns = 1, m_index+1
-!OLI_CODE           DO n = 0, m_index - ns + 1
-!OLI_CODE              IF (hitl(n) > hitl(n+1)) THEN ! swap order
-!OLI_CODE                 tmp = hitl(n)
-!OLI_CODE                 hitl(n) = hitl(n+1)
-!OLI_CODE                 hitl(n+1) = tmp
-!OLI_CODE                 tmp = aicetl(n)
-!OLI_CODE                 aicetl(n) = aicetl(n+1)
-!OLI_CODE                 aicetl(n+1) = tmp
-!OLI_CODE              END IF
-!OLI_CODE           END DO
-!OLI_CODE        END DO
-!OLI_CODE        
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        ! divide semi-filled layer into set of sublayers each vertically homogenous
-!OLI_CODE        !----------------------------------------------------------------
-!OLI_CODE        
-!OLI_CODE        !hitl(3)----------------------------------------------------------------
-!OLI_CODE        !                                                       | * * * * * * * *  
-!OLI_CODE        !                                                       |* * * * * * * * * 
-!OLI_CODE        !hitl(2)----------------------------------------------------------------
-!OLI_CODE        !                                    | * * * * * * * *  | * * * * * * * *  
-!OLI_CODE        !                                    |* * * * * * * * * |* * * * * * * * * 
-!OLI_CODE        !hitl(1)----------------------------------------------------------------
-!OLI_CODE        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * *  
-!OLI_CODE        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 
-!OLI_CODE        !hitl(0)----------------------------------------------------------------
-!OLI_CODE        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3)            
-!OLI_CODE        
-!OLI_CODE        ! move up over layers incrementing volume
-!OLI_CODE        DO n = 1, m_index+1
-!OLI_CODE           
-!OLI_CODE           area = sum(aicetl(:)) - &                 ! total area of sub-layer
-!OLI_CODE                (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
-!OLI_CODE           
-!OLI_CODE           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area
-!OLI_CODE           
-!OLI_CODE           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within
-!OLI_CODE              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - &
-!OLI_CODE                           alfan(1)
-!OLI_CODE              exit
-!OLI_CODE           ELSE  ! still in sub-layer below the sub-layer with the depth
-!OLI_CODE              rem_vol = rem_vol - vol
-!OLI_CODE           END IF
-!OLI_CODE           
-!OLI_CODE        END DO
-!OLI_CODE        
-!OLI_CODE       END IF
-!OLI_CODE     
-!OLI_CODE    END SUBROUTINE calc_hpond
-!OLI_CODE   
-!OLI_CODE 
-!OLI_CODE    SUBROUTINE permeability_phi(ticen, salin, vicen, perm)
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE       !!                ***  ROUTINE permeability_phi  ***
-!OLI_CODE       !!
-!OLI_CODE       !! ** Purpose :   Determine the liquid fraction of brine in the ice
-!OLI_CODE       !!                and its permeability
-!OLI_CODE       !!-------------------------------------------------------------------
-!OLI_CODE       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
-!OLI_CODE          ticen,  & ! energy of melting for each ice layer (J/m2)
-!OLI_CODE          salin
-!OLI_CODE 
-!OLI_CODE       REAL (wp), INTENT(IN) :: &
-!OLI_CODE          vicen     ! ice volume
-!OLI_CODE     
-!OLI_CODE       REAL (wp), INTENT(OUT) :: &
-!OLI_CODE          perm      ! permeability
-!OLI_CODE 
-!OLI_CODE       REAL (wp) ::   &
-!OLI_CODE          Sbr        ! brine salinity
-!OLI_CODE 
-!OLI_CODE       REAL (wp), DIMENSION(nlay_i) ::   &
-!OLI_CODE          Tin, &    ! ice temperature
-!OLI_CODE          phi       ! liquid fraction
-!OLI_CODE 
-!OLI_CODE       INTEGER :: k
-!OLI_CODE     
-!OLI_CODE       REAL (wp) :: &
-!OLI_CODE          c2    = 2.0_wp
-!OLI_CODE  
-!OLI_CODE       !-----------------------------------------------------------------
-!OLI_CODE       ! Compute ice temperatures from enthalpies using quadratic formula
-!OLI_CODE       !-----------------------------------------------------------------
-!OLI_CODE 
-!OLI_CODE       DO k = 1,nlay_i
-!OLI_CODE          Tin(k) = ticen(k) 
-!OLI_CODE       END DO
-!OLI_CODE 
-!OLI_CODE       !-----------------------------------------------------------------
-!OLI_CODE       ! brine salinity and liquid fraction
-!OLI_CODE       !-----------------------------------------------------------------
-!OLI_CODE 
-!OLI_CODE       IF (maxval(Tin-rtt) <= -c2) THEN
-!OLI_CODE 
-!OLI_CODE          DO k = 1,nlay_i
-!OLI_CODE             Sbr = - 1.2_wp                 &
-!OLI_CODE                   -21.8_wp     * (Tin(k)-rtt)    &
-!OLI_CODE                   - 0.919_wp   * (Tin(k)-rtt)**2 &
-!OLI_CODE                   - 0.01878_wp * (Tin(k)-rtt)**3
-!OLI_CODE             phi(k) = salin(k)/Sbr ! liquid fraction
-!OLI_CODE          END DO ! k
-!OLI_CODE        
-!OLI_CODE       ELSE
-!OLI_CODE 
-!OLI_CODE          DO k = 1,nlay_i
-!OLI_CODE             Sbr = -17.6_wp    * (Tin(k)-rtt)    &
-!OLI_CODE                   - 0.389_wp  * (Tin(k)-rtt)**2 &
-!OLI_CODE                   - 0.00362_wp* (Tin(k)-rtt)**3
-!OLI_CODE             phi(k) = salin(k)/Sbr ! liquid fraction
-!OLI_CODE          END DO
-!OLI_CODE 
-!OLI_CODE       END IF
-!OLI_CODE     
-!OLI_CODE       !-----------------------------------------------------------------
-!OLI_CODE       ! permeability
-!OLI_CODE       !-----------------------------------------------------------------
-!OLI_CODE 
-!OLI_CODE       perm = 3.0e-08_wp * (minval(phi))**3
-!OLI_CODE     
-!OLI_CODE    END SUBROUTINE permeability_phi
-!OLI_CODE   
-!OLI_CODE #else
-!OLI_CODE    !!----------------------------------------------------------------------
-!OLI_CODE    !!   Default option         Dummy Module          No LIM-3 sea-ice model
-!OLI_CODE    !!----------------------------------------------------------------------
-!OLI_CODE CONTAINS
-!OLI_CODE    SUBROUTINE lim_mp_init          ! Empty routine
-!OLI_CODE    END SUBROUTINE lim_mp_init   
-!OLI_CODE    SUBROUTINE lim_mp               ! Empty routine
-!OLI_CODE    END SUBROUTINE lim_mp
-!OLI_CODE    SUBROUTINE compute_mp_topo      ! Empty routine
-!OLI_CODE    END SUBROUTINE compute_mp_topo        
-!OLI_CODE    SUBROUTINE pond_area            ! Empty routine
-!OLI_CODE    END SUBROUTINE pond_area   
-!OLI_CODE    SUBROUTINE calc_hpond           ! Empty routine
-!OLI_CODE    END SUBROUTINE calc_hpond   
-!OLI_CODE    SUBROUTINE permeability_phy     ! Empty routine
-!OLI_CODE    END SUBROUTINE permeability_phy   
-!OLI_CODE #endif
-!OLI_CODE    !!======================================================================
-!OLI_CODE END MODULE limmp_topo
-!OLI_CODE 
+   END SUBROUTINE lim_mp_H12
+
+   SUBROUTINE lim_mp_init 
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE lim_mp_init   ***
+      !!
+      !! ** Purpose : Physical constants and parameters linked to melt ponds
+      !!      over sea ice
+      !!
+      !! ** Method  :  Read the nammp  namelist and check the melt pond  
+      !!       parameter values called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namelist nammp  
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ios, ioptio                 ! Local integer output status for namelist read
+      NAMELIST/nammp/  ln_pnd_H12, ln_pnd_fwb, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
+      !!-------------------------------------------------------------------
+
+      REWIND( numnam_ice_ref )              ! Namelist nammp  in reference namelist : Melt Ponds  
+      READ  ( numnam_ice_ref, nammp, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp  in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist nammp  in configuration namelist : Melt Ponds
+      READ  ( numnam_ice_cfg, nammp, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, nammp )
+      
+      IF(lwp) THEN                        ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
+         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namicethd_pnd:'
+         WRITE(numout,*) '      Evolutive melt pond fraction and depth (Holland et al 2012)  ln_pnd_H12 = ', ln_pnd_H12
+         WRITE(numout,*) '         Melt ponds store fresh water or not                       ln_pnd_fwb = ', ln_pnd_fwb
+         WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_Cst = ', ln_pnd_CST
+         WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
+         WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
+         WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
+      ENDIF
+      !
+      !                             !== set the choice of ice pond scheme ==!
+      ioptio = 0
+                                                            nice_pnd = np_pndNO
+      IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
+      IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
+      IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' )
+
+      SELECT CASE( nice_pnd )
+      CASE( np_pndNO )         
+         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when no ponds' ) ; ENDIF
+         IF(ln_pnd_alb) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
+      CASE( np_pndCST)
+         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when ln_pnd_CST=true' ) ; ENDIF
+      END SELECT
+      !
+   END SUBROUTINE lim_mp_init
+   
 #else
    !!----------------------------------------------------------------------
-   !!   Default option          Empty module           NO LIM sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_mp_init     ! Empty routine
-   END SUBROUTINE lim_mp_init
-   SUBROUTINE lim_mp          ! Empty routine
-   END SUBROUTINE lim_mp     
-   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
-   SUBROUTINE lim_mp_perm     ! Empty routine
-   END SUBROUTINE lim_mp_perm
+   !!   Default option          Empty module          NO ESIM sea-ice model
+   !!----------------------------------------------------------------------
 #endif 
 
