Index: /NEMO/branches/2020/SI3-05_MeltPonds_topo/cfgs/ORCA2_SAS_ICE/EXPREF/namelist_cfg
===================================================================
--- /NEMO/branches/2020/SI3-05_MeltPonds_topo/cfgs/ORCA2_SAS_ICE/EXPREF/namelist_cfg	(revision 13849)
+++ /NEMO/branches/2020/SI3-05_MeltPonds_topo/cfgs/ORCA2_SAS_ICE/EXPREF/namelist_cfg	(revision 13850)
@@ -23,5 +23,5 @@
    cn_exp      =  "ORCA2_SAS"  !  experience name
    nn_it000    =       1   !  first time step
-   nn_itend    =     100   !  last  time step (std 5475)
+   nn_itend    =    5840   !  last  time step (std 5840)
 /
 !-----------------------------------------------------------------------
Index: /NEMO/branches/2020/SI3-05_MeltPonds_topo/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/2020/SI3-05_MeltPonds_topo/cfgs/SHARED/namelist_ice_ref	(revision 13849)
+++ /NEMO/branches/2020/SI3-05_MeltPonds_topo/cfgs/SHARED/namelist_ice_ref	(revision 13850)
@@ -195,5 +195,6 @@
 !------------------------------------------------------------------------------
    ln_pnd            = .true.         !  activate melt ponds or not
-      ln_pnd_LEV     = .true.         !  level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012)
+      ln_pnd_TOPO    = .false.        !  topographic melt ponds 
+      ln_pnd_LEV     = .true.         !  level ice melt ponds 
          rn_apnd_min =   0.15         !     minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ??
          rn_apnd_max =   0.85         !     maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ??
Index: /NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/ice.F90	(revision 13849)
+++ /NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/ice.F90	(revision 13850)
@@ -208,5 +208,6 @@
    !                                     !!** ice-ponds namelist (namthd_pnd)
    LOGICAL , PUBLIC ::   ln_pnd           !: Melt ponds (T) or not (F)
-   LOGICAL , PUBLIC ::   ln_pnd_LEV       !: Melt ponds scheme from Holland et al (2012), Flocco et al (2007, 2010)
+   LOGICAL , PUBLIC ::   ln_pnd_TOPO      !: Topographic Melt ponds scheme (Flocco et al 2007, 2010)
+   LOGICAL , PUBLIC ::   ln_pnd_LEV       !: Simple melt pond scheme
    REAL(wp), PUBLIC ::   rn_apnd_min      !: Minimum ice fraction that contributes to melt ponds
    REAL(wp), PUBLIC ::   rn_apnd_max      !: Maximum ice fraction that contributes to melt ponds
@@ -308,4 +309,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t1_ice          !: temperature of the first layer          (ln_cndflx=T) [K]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   cnd_ice         !: effective conductivity of the 1st layer (ln_cndflx=T) [W.m-2.K-1]
+
+   ! meltwater arrays to save for melt ponds
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_i_sum_2d   !: surface melt (2d arrays for ponds)       [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_s_mlt_2d   !: snow surf melt (2d arrays for ponds)     [m]
 
    !!----------------------------------------------------------------------
@@ -458,4 +463,5 @@
       ii = ii + 1
       ALLOCATE( qtr_ice_bot(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) ,  &
+         &      dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) ,                    &
          &      h_i        (jpi,jpj,jpl) , a_i    (jpi,jpj,jpl) , v_i   (jpi,jpj,jpl) ,  &
          &      v_s        (jpi,jpj,jpl) , h_s    (jpi,jpj,jpl) , t_su  (jpi,jpj,jpl) ,  &
Index: /NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icestp.F90
===================================================================
--- /NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icestp.F90	(revision 13849)
+++ /NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icestp.F90	(revision 13850)
@@ -414,4 +414,7 @@
             wfx_pnd(ji,jj) = 0._wp
 
+            dh_i_sum_2d(:,:,:) = 0._wp
+            dh_s_mlt_2d(:,:,:) = 0._wp
+
             hfx_thd(ji,jj) = 0._wp   ;
             hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp
Index: /NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd.F90	(revision 13849)
+++ /NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd.F90	(revision 13850)
@@ -247,5 +247,4 @@
             IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- !
                               CALL ice_thd_dh                           ! Ice-Snow thickness   
-                              CALL ice_thd_pnd                          ! Melt ponds formation
                               CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping
             ENDIF
@@ -268,4 +267,6 @@
       IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft)
       !                   
+      IF ( ln_pnd  )          CALL ice_thd_pnd                      ! --- Melt ponds
+
       IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- !
       !
@@ -536,4 +537,7 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), cnd_ice_1d(1:npti), cnd_ice(:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), t1_ice_1d (1:npti), t1_ice (:,:,kl) )
+         ! ponds
+         CALL tab_1d_2d( npti, nptidx(1:npti), dh_i_sum  (1:npti) , dh_i_sum_2d(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), dh_s_mlt  (1:npti) , dh_s_mlt_2d(:,:,kl) )
          ! SIMIP diagnostics         
          CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d       (1:npti), t_si       (:,:,kl) )
Index: /NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90
===================================================================
--- /NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90	(revision 13849)
+++ /NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90	(revision 13850)
@@ -1,4 +1,5 @@
 MODULE icethd_pnd 
    !!======================================================================
+   !! --- closest to CICE version
    !!                     ***  MODULE  icethd_pnd   ***
    !!   sea-ice: Melt ponds on top of sea ice
@@ -20,4 +21,5 @@
    USE ice1D          ! sea-ice: thermodynamics variables
    USE icetab         ! sea-ice: 1D <==> 2D transformation
+   USE sbc_ice        ! surface energy budget
    !
    USE in_out_manager ! I/O manager
@@ -34,7 +36,8 @@
    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 ice pond scheme
-   INTEGER, PARAMETER ::   np_pndLEV = 2   ! Level ice pond scheme
+   INTEGER, PARAMETER ::   np_pndNO   = 0   ! No pond scheme
+   INTEGER, PARAMETER ::   np_pndCST  = 1   ! Constant ice pond scheme
+   INTEGER, PARAMETER ::   np_pndLEV  = 2   ! Level ice pond scheme
+   INTEGER, PARAMETER ::   np_pndTOPO = 3   ! Level ice pond scheme
 
    !! * Substitutions
@@ -52,4 +55,17 @@
       !!               
       !! ** Purpose :   change melt pond fraction and thickness
+      !!
+      !! Note: Melt ponds affect only radiative transfer for now
+      !!
+      !!       No heat, no salt.
+      !!       The melt water they carry is collected but 
+      !!       not removed from fw budget or released to the ocean
+      !!
+      !!       A wfx_pnd has been coded for diagnostic purposes
+      !!       It is not fully consistent yet.
+      !!
+      !!       The current diagnostic lacks a contribution from drainage
+      !!
+      !!       Each time wfx_pnd is updated, wfx_sum / wfx_snw_sum must be updated
       !!                
       !!-------------------------------------------------------------------
@@ -57,7 +73,9 @@
       SELECT CASE ( nice_pnd )
       !
-      CASE (np_pndCST)   ;   CALL pnd_CST    !==  Constant melt ponds  ==!
+      CASE (np_pndCST)   ;   CALL pnd_CST                              !==  Constant melt ponds  ==!
          !
-      CASE (np_pndLEV)   ;   CALL pnd_LEV    !==  Level ice melt ponds  ==!
+      CASE (np_pndLEV)   ;   CALL pnd_LEV                              !==  Level ice melt ponds  ==!
+         !
+      CASE (np_pndTOPO)  ;   CALL pnd_TOPO                  !&         !==  Topographic melt ponds  ==!
          !
       END SELECT
@@ -75,5 +93,5 @@
       !!                to non-zero values when t_su = 0C
       !!
-      !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
+      !! ** Tunable parameters : Pond fraction (rn_apnd) & depth (rn_hpnd)
       !!                
       !! ** Note   : Coupling with such melt ponds is only radiative
@@ -147,10 +165,10 @@
       !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min
       !! 
-      !! ** Note       :   mostly stolen from CICE
-      !!
-      !! ** References :   Flocco and Feltham (JGR, 2007)
-      !!                   Flocco et al       (JGR, 2010)
-      !!                   Holland et al      (J. Clim, 2012)
+      !! ** Note       :   mostly stolen from CICE but not only
+      !!
+      !! ** References :   Holland et al      (J. Clim, 2012)
+      !!                   
       !!-------------------------------------------------------------------
+      
       REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array
       !!
@@ -169,137 +187,1039 @@
       REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse
       !!
-      INTEGER  ::   ji, jk                            ! loop indices
+      INTEGER  ::   ji, jj, jk, jl                    ! loop indices
+      
       !!-------------------------------------------------------------------
+      
       z1_rhow   = 1._wp / rhow 
       z1_aspect = 1._wp / zaspect
       z1_Tp     = 1._wp / zTp 
-
-      DO ji = 1, npti
-         !                                                            !----------------------------------------------------!
-         IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction !
-            !                                                         !----------------------------------------------------!
-            !--- Remove ponds on thin ice or tiny ice fractions
-            a_ip_1d(ji)      = 0._wp
-            h_ip_1d(ji)      = 0._wp
-            h_il_1d(ji)      = 0._wp
-            !                                                         !--------------------------------!
-         ELSE                                                         ! Case ice thickness >= rn_himin !
-            !                                                         !--------------------------------!
-            v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness
-            v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji)
-            !
-            !------------------!
-            ! case ice melting !
-            !------------------!
-            !
-            !--- available meltwater for melt ponding ---!
-            zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji)
-            zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) = fraction of melt water that is not flushed
-            zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 
-            !
-            !--- overflow ---!
-            ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
-            !    a_ip_max = zfr_mlt * a_i
-            !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 
-            zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect
-            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
-
-            ! If pond depth exceeds half the ice thickness then reduce the pond volume
-            !    h_ip_max = 0.5 * h_i
-            !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 
-            zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji)
-            zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
+      
+      !-----------------------------------------------------------------------------------------------
+      !  Identify grid cells with ice
+      !-----------------------------------------------------------------------------------------------
+      at_i(:,:) = SUM( a_i, dim=3 )
+      !
+      npti = 0   ;   nptidx(:) = 0
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF ( at_i(ji,jj) > epsi10 ) THEN
+               npti = npti + 1
+               nptidx( npti ) = (jj - 1) * jpi + ji
+            ENDIF
+         END DO
+      END DO
+      
+      !-----------------------------------------------------------------------------------------------
+      ! Prepare 1D arrays
+      !-----------------------------------------------------------------------------------------------
+
+      IF( npti > 0 ) THEN
+      
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum   )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd          )
+      
+         DO jl = 1, jpl
+         
+            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d      (1:npti), a_i      (:,:,jl) )
+            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d      (1:npti), h_i      (:,:,jl) )
+            CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d     (1:npti), t_su     (:,:,jl) )
+            CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,jl) )
+            CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,jl) )
+            CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,jl) )
+
+            CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum    (1:npti), dh_i_sum_2d     (:,:,jl) )
+            CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt    (1:npti), dh_s_mlt_2d     (:,:,jl) )
             
-            !--- Pond growing ---!
-            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
-            !
-            !--- Lid melting ---!
-            IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0
-            !
-            !--- mass flux ---!
-            IF( zdv_mlt > 0._wp ) THEN
-               zfac = zdv_mlt * rhow * r1_rdtice                        ! melt pond mass flux < 0 [kg.m-2.s-1]
-               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
-               !
-               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux
-               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
-               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum)
-            ENDIF
-
-            !-------------------!
-            ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0)
-            !-------------------!
-            !
-            zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp )
-            !
-            !--- Pond contraction (due to refreezing) ---!
-            IF( ln_pnd_lids ) THEN
-               !
-               !--- Lid growing and subsequent pond shrinking ---! 
-               zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0
-                  &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors
-               
-               ! Lid growing
-               v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz )
-               
-               ! Pond shrinking
-               v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz )
-
-            ELSE
-               ! Pond shrinking
-               v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6)
-            ENDIF
-            !
-            !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
-            ! v_ip     = h_ip * a_ip
-            ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip)
-            a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i
-            h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji)
-
-            !---------------!            
-            ! Pond flushing !
-            !---------------!
-            ! height of top of the pond above sea-level
-            zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0
+            DO jk = 1, nlay_i
+               CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl)  )
+            END DO
             
-            ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010)
-            DO jk = 1, nlay_i
-               zsbr = - 1.2_wp                                  &
-                  &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    &
-                  &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 &
-                  &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3
-               ztmp(jk) = sz_i_1d(ji,jk) / zsbr
-            END DO
-            zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 )
-            
-            ! Do the drainage using Darcy's law
-            zdv_flush   = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji)
-            zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) )
-            v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush
-            
-            !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
-            a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i
-            h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji)
-
-            !--- Corrections and lid thickness ---!
-            IF( ln_pnd_lids ) THEN
-               !--- retrieve lid thickness from volume ---!
-               IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji)
-               ELSE                              ;   h_il_1d(ji) = 0._wp
-               ENDIF
-               !--- remove ponds if lids are much larger than ponds ---!
-               IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
+            !-----------------------------------------------------------------------------------------------
+            ! Go for ponds
+            !-----------------------------------------------------------------------------------------------
+
+
+            DO ji = 1, npti
+               !                                                            !----------------------------------------------------!
+               IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction !
+                  !                                                         !----------------------------------------------------!
+                  !--- Remove ponds on thin ice or tiny ice fractions
                   a_ip_1d(ji)      = 0._wp
                   h_ip_1d(ji)      = 0._wp
                   h_il_1d(ji)      = 0._wp
+                  !                                                         !--------------------------------!
+               ELSE                                                         ! Case ice thickness >= rn_himin !
+                  !                                                         !--------------------------------!
+                  v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness
+                  v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji)
+                  !
+                  !------------------!
+                  ! case ice melting !
+                  !------------------!
+                  !
+                  !--- available meltwater for melt ponding ---!
+                  zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji)
+                  zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) = fraction of melt water that is not flushed
+                  zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 
+                  !
+                  !--- overflow ---!
+                  ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
+                  !    a_ip_max = zfr_mlt * a_i
+                  !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 
+                  zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect
+                  zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
+
+                  ! If pond depth exceeds half the ice thickness then reduce the pond volume
+                  !    h_ip_max = 0.5 * h_i
+                  !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 
+                  zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) ! MV dimensions are wrong here or comment is unclear
+                  zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
+            
+                  !--- Pond growing ---!
+                  v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
+                  !
+                  !--- Lid melting ---!
+                  IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0
+                  !
+                  !--- mass flux ---!
+                  ! MV I would recommend to remove that
+                  ! Since melt ponds carry no freshwater there is no point in modifying water fluxes
+                  
+                  IF( zdv_mlt > 0._wp ) THEN
+                     zfac = zdv_mlt * rhow * r1_rdtice                        ! melt pond mass flux < 0 [kg.m-2.s-1]
+                     wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
+                     !
+                     zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux
+                     wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
+                     wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum)
+                  ENDIF
+
+                  !-------------------!
+                  ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0)
+                  !-------------------!
+                  !
+                  zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp )
+                  !
+                  !--- Pond contraction (due to refreezing) ---!
+                  IF( ln_pnd_lids ) THEN
+                     !
+                     !--- Lid growing and subsequent pond shrinking ---! 
+                     zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0
+                        &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors
+               
+                     ! Lid growing
+                     v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz )
+               
+                     ! Pond shrinking
+                     v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz )
+
+                  ELSE
+                     ! Pond shrinking
+                     v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6)
+                  ENDIF
+                  !
+                  !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
+                  ! v_ip     = h_ip * a_ip
+                  ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip)
+                  a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i
+                  h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji)
+
+                  !---------------!            
+                  ! Pond flushing !
+                  !---------------!
+                  ! height of top of the pond above sea-level
+                  zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0
+            
+                  ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010)
+                  DO jk = 1, nlay_i
+                     ! MV Assur is inconsistent with SI3
+                     zsbr = - 1.2_wp                                  &
+                        &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    &
+                        &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 &
+                        &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3
+                     ! MV linear expression more consistent & simpler              zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt
+                     ztmp(jk) = sz_i_1d(ji,jk) / zsbr
+                  END DO
+                  zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 )
+            
+                  ! Do the drainage using Darcy's law
+                  zdv_flush   = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji)
+                  zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) )
+                  zdv_flush   = 0._wp ! MV remove pond drainage for now
+                  v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush
+            
+                  ! MV --- why pond drainage does not give back water into freshwater flux ?
+            
+                  !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
+                  a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i
+                  h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji)
+
+                  !--- Corrections and lid thickness ---!
+                  IF( ln_pnd_lids ) THEN
+                     !--- retrieve lid thickness from volume ---!
+                     IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji)
+                     ELSE                              ;   h_il_1d(ji) = 0._wp
+                     ENDIF
+                     !--- remove ponds if lids are much larger than ponds ---!
+                     IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
+                        a_ip_1d(ji)      = 0._wp
+                        h_ip_1d(ji)      = 0._wp
+                        h_il_1d(ji)      = 0._wp
+                     ENDIF
+                  ENDIF
+                  !
                ENDIF
+               
+            END DO ! ji
+
+            !-----------------------------------------------------------------------------------------------
+            ! Retrieve 2D arrays
+            !-----------------------------------------------------------------------------------------------
+            
+            v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
+            v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti)
+            CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,jl) )
+            CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,jl) )
+            CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,jl) )
+            CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d     (1:npti), v_ip     (:,:,jl) )
+            CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d     (1:npti), v_il     (:,:,jl) )
+            DO jk = 1, nlay_i
+               CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl)  )
+            END DO
+   
+         END DO ! ji
+                  
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        )
+                  
+      !
+      ENDIF
+      
+   END SUBROUTINE pnd_LEV
+   
+   SUBROUTINE pnd_TOPO    
+                                         
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE pnd_TOPO  ***
+      !!
+      !! ** Purpose : Compute melt pond evolution
+      !!
+      !! ** 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.
+      !!
+      !!-------------------------------------------------------------------
+
+      ! local variables
+      REAL(wp) :: &
+         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)
+         zfsurf,  &   ! net heat flux, excluding conduction and transmitted radiation (W/m2)
+         zTp,     &   ! pond freezing temperature (C)
+         zrhoi_L, &   ! volumetric latent heat of sea ice (J/m^3)
+         zdvn   , &   ! change in melt pond volume for fresh water budget (not used for now)
+         zfr_mlt, &   ! fraction and volume of available meltwater retained for melt ponding
+         z1_rhow, &   ! inverse water density
+         zpond  , &   ! dummy variable
+         zdum         ! dummy variable
+
+      REAL(wp), PARAMETER :: &
+         zhi_min = 0.1_wp        , & ! minimum ice thickness with ponds (m) 
+         zTd     = 0.15_wp       , & ! temperature difference for freeze-up (C)
+         zvp_min = 1.e-4_wp        ! minimum pond volume (m)
+
+      REAL(wp), DIMENSION(jpi,jpj) ::   zvolp !! total melt pond water available before redistribution and drainage
+
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i
+
+      INTEGER  ::   ji, jj, jk, jl                    ! loop indices
+
+      INTEGER  ::   i_test
+
+      ! Note
+      ! equivalent for CICE translation
+      ! a_ip      -> apond
+      ! a_ip_frac -> apnd
+      
+      !---------------------------------------------------------------
+      ! Initialize
+      !---------------------------------------------------------------
+
+      ! Parameters & constants (move to parameters)
+      zrhoi_L   = rhoi * rLfus      ! volumetric latent heat (J/m^3)
+      zTp       = rt0 - 0.15_wp          ! pond freezing point, slightly below 0C (ponds are bid saline)
+      z1_rhow   = 1._wp / rhow 
+
+      ! Set required ice variables (hard-coded here for now)
+!      zfpond(:,:) = 0._wp          ! contributing freshwater flux (?) 
+      
+      at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction
+      vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area
+      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area
+      vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area
+      
+      ! thickness
+      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
+      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp
+      END WHERE
+      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
+
+      !-----------------------------------------------------------------
+      ! Start pond loop
+      !-----------------------------------------------------------------      
+
+      zvolp(:,:) = 0._wp
+
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+         
+            !--------------------------------------------------------------
+            ! Collect meltwater on ponds
+            !--------------------------------------------------------------
+            DO jl = 1, jpl
+            
+               IF ( a_i(ji,jj,jl) > epsi10 ) THEN
+            
+                  !--- Available meltwater for melt ponding ---!
+                  zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * a_i(ji,jj,jl) !  = ( 1 - r ) = fraction of melt water that is not flushed
+                  zdum    = -( dh_i_sum_2d(ji,jj,jl)*rhoi + dh_s_mlt_2d(ji,jj,jl)*rhos ) * z1_rhow * a_i(ji,jj,jl)
+                  zpond   = zfr_mlt * zdum ! check
+
+                  !--- Create possible new ponds
+                  ! if pond does not exist, create new pond over full ice area
+                  IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN
+                     h_ip(ji,jj,jl)       = 0._wp
+                     a_ip_frac(ji,jj,jl) = 1.0_wp    ! pond fraction of sea ice (apnd for CICE)
+                  ENDIF
+                  
+                  !--- Deepen existing ponds before redistribution and drainage(later on)
+                  h_ip(ji,jj,jl) = ( zpond + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) ) / a_ip_frac(ji,jj,jl)
+                  zvolp(ji,jj)   = zvolp(ji,jj) + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) 
+!                 zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) !useless
+                   
+               ENDIF ! a_i
+
+            END DO
+
+            h_ip(ji,jj,:)  = 0._wp ! pond thickness
+            a_ip(ji,jj,:)  = 0._wp ! pond fraction per grid area      
+
+            !-----------------------------------------------------------------
+            ! Identify grid cells with ice
+            !-----------------------------------------------------------------
+
+            IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. vt_ip(ji,jj) > zvp_min *at_i(ji,jj) ) THEN
+      
+               !----------------------------------------------------------------------------
+               ! Calculate pond area and depth (redistribution of melt water on topography)
+               !----------------------------------------------------------------------------
+!              CALL ice_thd_pnd_area( at_i(ji,jj) , vt_i(ji,jj), &
+!                                     a_i (ji,jj,:), v_i(ji,jj,:), v_s(ji,jj,:), &
+!                                     t_i(ji,jj,:,:), sz_i(ji,jj,:,:), &
+!                                     v_ip(ji,jj,:) , zvolp(ji,jj), &
+!                                     a_ip(ji,jj,:), h_ip(ji,jj,:), zdvn)
+                                   
+!               zfpond(ji,jj) = zfpond(ji,jj) - zdvn ! ---> zdvn is drainage (useless for now since we dont change fw budget, no ?)
+            
+               !--------------------------
+               ! Pond lid growth and melt
+               !--------------------------
+               ! Mean surface temperature
+               zTavg = 0._wp
+               DO jl = 1, jpl
+                  zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl)
+               END DO
+               zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here
+         
+               DO jl = 1, jpl-1
+            
+                  IF ( v_il(ji,jj,jl) > epsi10 ) THEN
+               
+                     !----------------------------------------------------------------
+                     ! Lid melting: floating upper ice layer melts in whole or part
+                     !----------------------------------------------------------------
+                     ! Use Tsfc for each category
+                     IF ( t_su(ji,jj,jl) > zTp ) THEN
+
+                        zdvice = min( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl))
+                        IF ( zdvice > epsi10 ) then
+                           v_il (ji,jj,jl) = v_il (ji,jj,jl)   - zdvice
+                           v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)    + zdvice
+!                          zvolp(ji,jj)     = zvolp(ji,jj)     + zdvice ! pointless to calculate total volume (done in icevar)
+!                          zfpond(ji,jj)    = fpond(ji,jj)     + zdvice ! pointless to follow fw budget (ponds have no fw)
+                     
+                           IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN
+                           ! ice lid melted and category is pond covered
+                              v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)  + v_il(ji,jj,jl) 
+!                             zfpond(ji,jj)    = zfpond (ji,jj)    + v_il(ji,jj,jl) 
+                              v_il(ji,jj,jl)   = 0._wp
+                           ENDIF
+                           h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here
+                        ENDIF
+                  
+                     !----------------------------------------------------------------
+                     ! Freeze pre-existing lid 
+                     !----------------------------------------------------------------
+
+                     ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp
+
+                        ! differential growth of base of surface floating ice layer
+                        zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0   
+                        zomega = rcnd_i * zdTice / zrhoi_L
+                        zdHui  = SQRT( 2._wp * zomega * rdt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) &
+                               - v_il(ji,jj,jl) / a_i(ji,jj,jl)
+                        zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
+                  
+                        IF ( zdvice > epsi10 ) THEN
+                           v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice
+                           v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice
+!                          zvolp(ji,jj)    = zvolp(ji,jj)     - zdvice
+!                          zfpond(ji,jj)   = zfpond(ji,jj)    - zdvice
+                           h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
+                        ENDIF
+                  
+                     ENDIF ! Tsfcn(i,j,n)
+
+                  !----------------------------------------------------------------
+                  ! Freeze new lids
+                  !----------------------------------------------------------------
+                  !  upper ice layer begins to form
+                  ! note: albedo does not change
+
+                  ELSE ! v_il < epsi10
+                    
+                     ! thickness of newly formed ice
+                     ! the surface temperature of a meltpond is the same as that
+                     ! of the ice underneath (0C), and the thermodynamic surface 
+                     ! flux is the same
+                     
+                     !!!!! FSURF !!!!!
+                     !!! we need net surface energy flux, excluding conduction
+                     !!! fsurf is summed over categories in CICE
+                     !!! we have the category-dependent flux, let us use it ?
+                     zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl)                     
+                     zdHui  = MAX ( -zfsurf * rdt_ice/zrhoi_L , 0._wp )
+                     zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
+                     IF ( zdvice > epsi10 ) THEN
+                        v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice
+                        v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice
+!                       zvolp(ji,jj)     = zvolp(ji,jj)     - zdvice
+!                       zfpond(ji,jj)    = zfpond(ji,jj)    - zdvice
+                        h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) ! in principle, this is useless as h_ip is computed in icevar
+                     ENDIF
+               
+                  ENDIF  ! v_il
+            
+               END DO ! jl
+
+            ELSE  ! remove ponds on thin ice
+
+               v_ip(ji,jj,:) = 0._wp
+               v_il (ji,jj,:)  = 0._wp
+!              zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj)
+!		            zvolp(ji,jj)    = 0._wp         
+
             ENDIF
-            !
+
+      !---------------------------------------------------------------
+      ! remove ice lid if there is no liquid pond
+      ! v_il may be nonzero on category ncat due to dynamics
+      !---------------------------------------------------------------
+
+      DO jl = 1, jpl
+      
+         IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 &
+                             .AND. v_il (ji,jj,jl) > epsi10) THEN
+            v_il(ji,jj,jl) = 0._wp
+         ENDIF
+
+         ! reload tracers
+         IF ( a_ip(ji,jj,jl) > epsi10) THEN
+            h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar
+         ELSE
+            v_il(ji,jj,jl) = 0._wp
+            h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as a_ip_frac computed in icevar
+         ENDIF
+         IF ( a_ip(ji,jj,jl) > epsi10 ) THEN
+            a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar
+            !h_ip(ji,jj,jl) = zhpondn(ji,jj,jl)
+         ELSE
+            a_ip_frac(ji,jj,jl) = 0._wp
+            h_ip(ji,jj,jl)      = 0._wp ! MV in principle, useless as computed in icevar
+            h_il(ji,jj,jl)      = 0._wp ! MV in principle, useless as omputed in icevar
          ENDIF
          
-      END DO
-      !
-   END SUBROUTINE pnd_LEV
-
+      END DO       ! jl
+
+      END DO   ! jj
+      END DO   ! ji
+
+   END SUBROUTINE pnd_TOPO
+
+   
+       SUBROUTINE ice_thd_pnd_area(aice,  vice, &
+                           aicen, vicen, vsnon, ticen, &
+                           salin, zvolpn, zvolp,         &
+                           zapondn,zhpondn,dvolp)
+
+       !!-------------------------------------------------------------------
+       !!                ***  ROUTINE ice_thd_pnd_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
+
+      !-----------|
+      !           |
+      !           |-----------|
+      !___________|___________|______________________________________sea-level
+      !           |           |
+      !           |           |---^--------|
+      !           |           |   |        |
+      !           |           |   |        |-----------|              |-------
+      !           |           |   |alfan(n)|           |              |
+      !           |           |   |        |           |--------------|
+      !           |           |   |        |           |              |
+      !---------------------------v-------------------------------------------
+      !           |           |   ^        |           |              |
+      !           |           |   |        |           |--------------|
+      !           |           |   |betan(n)|           |              |
+      !           |           |   |        |-----------|              |-------
+      !           |           |   |        |
+      !           |           |---v------- |
+      !           |           |
+      !           |-----------|
+      !           |
+      !-----------|
+
+       !-------------------------------------------------------------------
+       ! initialize
+       !-------------------------------------------------------------------
+
+       DO n = 1, jpl
+
+          zapondn(n) = 0._wp
+          zhpondn(n) = 0._wp
+
+          !----------------------------------------
+          ! X) compute the effective snow fraction
+          !----------------------------------------
+          IF (aicen(n) < epsi10)  THEN
+             hicen(n) =  0._wp
+             hsnon(n) =  0._wp
+             reduced_aicen(n) = 0._wp
+             asnon(n) = 0._wp         !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables
+          ELSE
+             hicen(n) = vicen(n) / aicen(n)
+             hsnon(n) = vsnon(n) / aicen(n)
+             reduced_aicen(n) = 1._wp ! n=jpl
+
+             !js: initial code in NEMO_DEV
+             !IF (n < jpl) reduced_aicen(n) = aicen(n) &
+             !                     * (-0.024_wp*hicen(n) + 0.832_wp)
+
+             !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large
+             IF (n < jpl) reduced_aicen(n) = aicen(n) &
+                                  * max(0.2_wp,(-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
+             ! OLI: it probably doesn't
+          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 ? => OLI : Coz' Chuck Norris said so...
+
+          alfan(n) = 0.6 * hicen(n)
+          betan(n) = 0.4 * hicen(n)
+
+          cum_max_vol(n)     = 0._wp
+          cum_max_vol_tmp(n) = 0._wp
+
+       END DO ! jpl
+
+       cum_max_vol_tmp(0) = 0._wp
+       drain = 0._wp
+       dvolp = 0._wp
+
+       !----------------------------------------------------------
+       ! 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_out
+       ! END MV
+       !js 18/04/19: XXX do something about this flux thing
+
+       DO n = 1, jpl-1 ! last category can not hold any volume
+
+          IF (alfan(n+1) >= alfan(n) .AND. alfan(n+1) > 0._wp ) 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) &
+                   - rhos/rhow * &   ! 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)), 0._wp)
+             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 = 0._wp
+          END IF
+       END IF
+
+       ! height and area corresponding to the remaining volume
+
+      CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, 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 ?
+          zhpondn(n) = max((hpond - alfan(n) + alfan(1)), 0._wp)      !js: from CICE 5.1.2
+          zapondn(n) = reduced_aicen(n)
+          ! in practise, pond fraction depends on the empirical snow fraction
+          ! so in turn on ice thickness
+       END DO
+       !zapond = sum(zapondn(1:m_index))     !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
+
+       !------------------------------------------------------------------------
+       ! Drainage through brine network (permeability)
+       !------------------------------------------------------------------------
+       !!! drainage due to ice permeability - Darcy's law
+
+       ! sea water level
+       msno =0._wp 
+       DO n=1,jpl
+         msno = msno + vsnon(n) * rhos
+       END DO
+       floe_weight = (msno + rhoi*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, 0._wp)
+
+       ! drain if ice is permeable
+       permflag = 0
+       IF (pressure_head > 0._wp) THEN
+          DO n = 1, jpl-1
+             IF (hicen(n) /= 0._wp) THEN
+             !IF (hicen(n) > 0._wp) THEN           !js: from CICE 5.1.2
+                perm = 0._wp ! MV ugly dummy patch
+                CALL ice_thd_pnd_perm(ticen(:,n), salin(:,n), perm)
+                IF (perm > 0._wp) permflag = 1
+
+                drain = perm*zapondn(n)*pressure_head*rdt_ice / &
+                                         (viscosity*hicen(n))
+                dvolp = dvolp + min(drain, zvolp)
+                zvolp = max(zvolp - drain, 0._wp)
+                IF (zvolp < epsi10) THEN
+                   dvolp = dvolp + zvolp
+                   zvolp = 0._wp
+                END IF
+            END IF
+         END DO
+
+          ! adjust melt pond dimensions
+          IF (permflag > 0) THEN
+             ! recompute pond depth
+            CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, 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
+             !zapond = sum(zapondn(1:m_index))       !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
+          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
+                   - (rhos/rhow) * 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) = 0._wp 
+           zhpondn(m_index) = 0._wp
+           zapondn(m_index) = 0._wp
+           ! 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) = 0._wp 
+             zvolpn(n)  = 0._wp
+             zapondn(n) = 0._wp
+          END IF
+       END DO
+       DO n = m_index+1, jpl
+          zhpondn(n) = 0._wp 
+          zapondn(n) = 0._wp 
+          zvolpn (n) = 0._wp 
+       END DO
+
+    END SUBROUTINE ice_thd_pnd_area
+
+
+    SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index)
+       !!-------------------------------------------------------------------
+       !!                ***  ROUTINE ice_thd_pnd_depth  ***
+       !!
+       !! ** Purpose :   Compute melt pond depth
+       !!-------------------------------------------------------------------
+
+       REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
+          aicen, &
+          asnon, &
+          hsnon, &
+          alfan, &
+          cum_max_vol
+
+       REAL (wp), INTENT(IN) :: &
+          zvolp
+
+       REAL (wp), INTENT(OUT) :: &
+          hpond
+
+       INTEGER, INTENT(OUT) :: &
+          m_index
+
+       INTEGER :: n, ns
+
+       REAL (wp), DIMENSION(0:jpl+1) :: &
+          hitl, &
+          aicetl
+
+       REAL (wp) :: &
+          rem_vol, &
+          area, &
+          vol, &
+          tmp, &
+          z0   = 0.0_wp
+
+       !----------------------------------------------------------------
+       ! hpond is zero if zvolp is zero - have we fully drained?
+       !----------------------------------------------------------------
+
+       IF (zvolp < epsi10) THEN
+        hpond = z0
+        m_index = 0
+       ELSE
+
+        !----------------------------------------------------------------
+        ! Calculate the category where water fills up to
+        !----------------------------------------------------------------
+
+        !----------|
+        !          |
+        !          |
+        !          |----------|                                     -- --
+        !__________|__________|_________________________________________ ^
+        !          |          |             rem_vol     ^                | Semi-filled
+        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer
+        !          |          |          |              |
+        !          |          |          |              |hpond
+        !          |          |          |----------|   |     |-------
+        !          |          |          |          |   |     |
+        !          |          |          |          |---v-----|
+        !          |          | m_index  |          |         |
+        !-------------------------------------------------------------
+
+        m_index = 0  ! 1:m_index categories have water in them
+        DO n = 1, jpl
+           IF (zvolp <= cum_max_vol(n)) THEN
+              m_index = n
+              IF (n == 1) THEN
+                 rem_vol = zvolp
+              ELSE
+                 rem_vol = zvolp - cum_max_vol(n-1)
+              END IF
+              exit ! to break out of the loop
+           END IF
+        END DO
+        m_index = min(jpl-1, m_index)
+
+        !----------------------------------------------------------------
+        ! semi-filled layer may have m_index different snow in it
+        !----------------------------------------------------------------
+
+        !-----------------------------------------------------------  ^
+        !                                                             |  alfan(m_index+1)
+        !                                                             |
+        !hitl(3)-->                             |----------|          |
+        !hitl(2)-->                |------------| * * * * *|          |
+        !hitl(1)-->     |----------|* * * * * * |* * * * * |          |
+        !hitl(0)-->-------------------------------------------------  |  ^
+        !                various snow from lower categories          |  |alfa(m_index)
+
+        ! hitl - heights of the snow layers from thinner and current categories
+        ! aicetl - area of each snow depth in this layer
+
+        hitl(:) = z0
+        aicetl(:) = z0
+        DO n = 1, m_index
+           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
+                                  alfan(m_index+1) - alfan(m_index)), z0)
+           aicetl(n) = asnon(n)
+
+           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
+        END DO
+
+        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
+        aicetl(m_index+1) = z0
+
+        !----------------------------------------------------------------
+        ! reorder array according to hitl
+        ! snow heights not necessarily in height order
+        !----------------------------------------------------------------
+
+        DO ns = 1, m_index+1
+           DO n = 0, m_index - ns + 1
+              IF (hitl(n) > hitl(n+1)) THEN ! swap order
+                 tmp = hitl(n)
+                 hitl(n) = hitl(n+1)
+                 hitl(n+1) = tmp
+                 tmp = aicetl(n)
+                 aicetl(n) = aicetl(n+1)
+                 aicetl(n+1) = tmp
+              END IF
+           END DO
+        END DO
+
+        !----------------------------------------------------------------
+        ! divide semi-filled layer into set of sublayers each vertically homogenous
+        !----------------------------------------------------------------
+
+        !hitl(3)----------------------------------------------------------------
+        !                                                       | * * * * * * * *
+        !                                                       |* * * * * * * * *
+        !hitl(2)----------------------------------------------------------------
+        !                                    | * * * * * * * *  | * * * * * * * *
+        !                                    |* * * * * * * * * |* * * * * * * * *
+        !hitl(1)----------------------------------------------------------------
+        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * *
+        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * *
+        !hitl(0)----------------------------------------------------------------
+        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3)
+
+        ! move up over layers incrementing volume
+        DO n = 1, m_index+1
+
+           area = sum(aicetl(:)) - &                 ! total area of sub-layer
+                (rhos/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
+
+           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area
+
+           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within
+              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1)
+
+              exit
+           ELSE  ! still in sub-layer below the sub-layer with the depth
+              rem_vol = rem_vol - vol
+           END IF
+
+        END DO
+
+       END IF
+
+    END SUBROUTINE ice_thd_pnd_depth
+
+
+    SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm)
+       !!-------------------------------------------------------------------
+       !!                ***  ROUTINE ice_thd_pnd_perm ***
+       !!
+       !! ** Purpose :   Determine the liquid fraction of brine in the ice
+       !!                and its permeability
+       !!-------------------------------------------------------------------
+
+       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
+          ticen,  &     ! internal ice temperature (K)
+          salin         ! salinity (ppt)     !js: ppt according to cice
+
+       REAL (wp), INTENT(OUT) :: &
+          perm      ! permeability
+
+       REAL (wp) ::   &
+          Sbr       ! brine salinity
+
+       REAL (wp), DIMENSION(nlay_i) ::   &
+          Tin, &    ! ice temperature
+          phi       ! liquid fraction
+
+       INTEGER :: k
+
+       !-----------------------------------------------------------------
+       ! Compute ice temperatures from enthalpies using quadratic formula
+       !-----------------------------------------------------------------
+
+       DO k = 1,nlay_i
+          Tin(k) = ticen(k) - rt0   !js: from K to degC
+       END DO
+
+       !-----------------------------------------------------------------
+       ! brine salinity and liquid fraction
+       !-----------------------------------------------------------------
+
+       DO k = 1, nlay_i
+       
+          Sbr    = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus)
+          ! Best expression to date is that one
+          Sbr  = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3
+          phi(k) = salin(k) / Sbr
+          
+       END DO
+
+       !-----------------------------------------------------------------
+       ! permeability
+       !-----------------------------------------------------------------
+
+       perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007)
+
+   END SUBROUTINE ice_thd_pnd_perm
+   
+   
+!----------------------------------------------------------------------------------------------------------------------
 
    SUBROUTINE ice_thd_pnd_init 
@@ -319,4 +1239,5 @@
       NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, &
          &                          ln_pnd_CST , rn_apnd, rn_hpnd,         &
+         &                          ln_pnd_TOPO ,                          &
          &                          ln_pnd_lids, ln_pnd_alb
       !!-------------------------------------------------------------------
@@ -336,4 +1257,5 @@
          WRITE(numout,*) '   Namelist namicethd_pnd:'
          WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd
+         WRITE(numout,*) '         Topographic melt pond scheme                             ln_pnd_TOPO  = ', ln_pnd_TOPO
          WRITE(numout,*) '         Level ice melt pond scheme                               ln_pnd_LEV   = ', ln_pnd_LEV
          WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min
@@ -351,4 +1273,5 @@
       IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
       IF( ln_pnd_LEV  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndLEV    ;   ENDIF
+      IF( ln_pnd_TOPO ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndTOPO   ;   ENDIF
       IF( ioptio /= 1 )   &
          & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV or ln_pnd_CST)' )
