Index: NEMO/branches/2020/SI3-05_MP/src/ICE/icethd_pnd.F90
===================================================================
--- NEMO/branches/2020/SI3-05_MP/src/ICE/icethd_pnd.F90	(revision 13808)
+++ NEMO/branches/2020/SI3-05_MP/src/ICE/icethd_pnd.F90	(revision 13809)
@@ -34,7 +34,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
@@ -60,4 +61,6 @@
          !
       CASE (np_pndLEV)   ;   CALL pnd_LEV    !==  Level ice melt ponds  ==!
+         !
+      CASE (np_pndTOPO)  ;   CALL pnd_TOPO   !==  Topographic melt ponds  ==!
          !
       END SELECT
@@ -153,4 +156,5 @@
       !!                   Holland et al      (J. Clim, 2012)
       !!-------------------------------------------------------------------
+      
       REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array
       !!
@@ -170,5 +174,7 @@
       !!
       INTEGER  ::   ji, jk                            ! loop indices
+      
       !!-------------------------------------------------------------------
+      
       z1_rhow   = 1._wp / rhow 
       z1_aspect = 1._wp / zaspect
@@ -196,5 +202,5 @@
             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? 
+            zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! >0, max for roundoff errors? 
             !
             !--- overflow ---!
@@ -208,5 +214,5 @@
             !    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)
+            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) ) )
             
@@ -216,10 +222,17 @@
             !--- 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
+               ! MV add comment on what that mass flux means
+               ! water removed from fw flux due to melt pond growth
                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
                !
+               ! MV 
+               ! why surface melt and snow fluxes must be adjusted is not clear
+               ! sounds like things are counted twice
+               ! 
                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)
@@ -264,4 +277,9 @@
             
             ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010)
+            ! MV the expression for zsbr has wrong sign!!!!
+            ! MV note - the sign-corrected expression is inconsistent 
+            ! with the rest of the SI3 code which assumes linear liquidus
+            ! best expression (most consistent for SI3 should be)
+            ! zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt
             DO jk = 1, nlay_i
                zsbr = - 1.2_wp                                  &
@@ -269,7 +287,8 @@
                   &   - 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
+               
+               ztmp(jk) = sz_i_1d(ji,jk) / zsbr ! MV -> this is brine fraction and as it reads is always negative
             END DO
-            zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 )
+            zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) ! MV at present since ztmp is negative, this is always zero
             
             ! Do the drainage using Darcy's law
@@ -301,5 +320,915 @@
       !
    END SUBROUTINE pnd_LEV
-
+   
+   SUBROUTINE pnd_TOPO    (aice,      aicen,     &
+                           vice,      vicen,     &
+                           vsnon,                &
+                           ticen,     salin,     &
+                           a_ip_frac, h_ip,      &
+                                      Tsfc )
+                                         
+      !!-------------------------------------------------------------------
+      !!                ***  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.
+      !!
+      !!-------------------------------------------------------------------
+
+       !js 190423: the lid on melt ponds appears only in the analog subroutine of CICE 5.1.2
+
+       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 temperature per category (K)
+          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)
+
+       REAL (wp), DIMENSION (jpi,jpj,jpl) :: &
+          zrfrac, & ! fraction of available meltwater retained for melt ponding
+          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)
+          zwfx_tmp    ! temporary array for melt water
+
+       REAL (wp) :: &
+          zhi,      & ! ice thickness (m)
+          zTavg,    & ! mean surface temperature across categories (C)
+          z1_rhow, & ! inverse freshwater density
+          zdTs,     & ! temperature difference for freeze-up (C)
+          zvpold,   & ! dummy pond volume
+          zdvn        ! change in melt pond volume for fresh water budget
+          
+       INTEGER, DIMENSION (jpi*jpj) :: &
+          indxi, indxj    ! compressed indices for cells with ice melting
+
+       INTEGER :: ij,icells,ji,jj,jl ! loop indices
+
+       REAL (wp), PARAMETER :: &
+          zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3)
+          zTp       = -0.15_wp, & ! pond freezing temperature (C)
+          zmin_volp = 1.e-4_wp, & ! minimum pond volume (m)
+          zrexp     = 0.01_wp,  & ! constant melt pond freeze-up rate
+          z01       = 0.01_wp,  &
+          z25       = 0.25_wp,  &
+          z5        = 0.5_wp
+
+      z1_rhow     = 1. / rhow
+
+      !---------------------------------------------------------------
+      ! Initialization
+      !---------------------------------------------------------------
+      zhpondn  (:,:,:) = 0._wp
+      zapondn  (:,:,:) = 0._wp
+      zvolpn   (:,:,:) = 0._wp
+
+      zTsfcn(:,:,:) = Tsfc(:,:,:) - rt0         ! Convert in Celsius
+
+      IF ( ln_pnd_fw ) THEN
+         v_ip_b(:,:,:) = v_ip(:,:,:)
+      ELSE
+         v_ip_b(:,:,:) = 0._wp
+      ENDIF
+
+      !------------------------------------------------------------------
+      ! Available melt water for melt ponding and corresponding fraction
+      !------------------------------------------------------------------
+      !js 03/05/19 unset restriction on sign of wfx_pnd_in;  mask values close to zero for future division
+      !wfx_pnd_in(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:), 0._wp )        ! available meltwater for melt ponding
+      !wfx_pnd_in(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:) , epsi10 ) * MAX( 0._wp, SIGN( 1._wp, wfx_sum(:,:) + wfx_snw_sum(:,:) - epsi10 ) )
+      !wfx_pnd_in(:,:) = (wfx_sum(:,:) + wfx_snw_sum(:,:)) * MAX( 0._wp, SIGN( 1._wp, ABS(wfx_sum(:,:) + wfx_snw_sum(:,:)) - epsi10 ) )
+
+      ! MV
+      ! NB: wfx_pnd_in can be slightly negative for very small values (why?)
+      ! This can in some occasions give negative
+      ! v_ip in the first category, which then gives crazy pond
+      ! fractions and crashes the code as soon as the melt-pond
+      ! radiative coupling is activated
+      ! if we understand and remove why wfx_sum or wfx_snw could be
+      ! negative, then, we can remove the MAX
+      ! NB: I now changed to wfx_snw_sum, this may fix the problem.
+      ! We should check
+
+
+      ! OLI 07/2017: when we (MV & OL) first started the inclusion of melt
+      ! ponds in the model, we removed the Holland et al. (2012, see
+      ! CESM scheme above) parameterization. I put it back here,
+      ! because I think it is needed. In summary, the sinks of FW for
+      ! ponds are: 1/ Runoff through cracks/leads => depends on the
+      !               total ice area only
+      !            2/ overflow, including Lüthje et al. (2006) limitation
+      !               (max a_ip fraction function of h_i). This is in
+      !               fact an other form of runoff that depends on the
+      !               ITD
+      !            3/ Flushing - losses by permeability
+      !            4/ Refreezing
+      !            5/ Removal of ponds on thin ice
+      ! I think 1 is needed because it is different from 2. However,
+      ! test runs should/could be done, to check the sensitivity and
+      ! the real usefulness of that stuff.
+      ! Note : the Holland et al. param was wrongly wired in NEMO3.1 (using
+      ! a_i instead of at_i), which might well explain why I had a too
+      ! weak melt pond cover in my simulations (compared to MODIS, in
+      ! situ obs. and CICE simulations.
+
+      !js 23/04/19: rewired back to a fraction with a_i
+      zrfrac(:,:,:) = rn_pnd_fracmin + ( rn_pnd_fracmax - rn_pnd_fracmin ) * aicen(:,:,:)
+      zwfx_tmp(:,:) = 0._wp
+
+! MV ---> use expression from level-ice ponds, clarify what is needed
+
+      !--- Add retained melt water to melt ponds
+      ! v_ip should never be negative, otherwise code crashes
+      ! Rq MV: as far as I saw, UM5 can create very small negative v_ip values
+      ! hence I added the max, which was not required with Prather (1 yr run)
+      ! OLI: Here I use vt_ip, so I don't know if the max is
+      ! required...
+      !zvolp(:,:) = MAX(vt_ip(:,:),0._wp) + zrfrac(:,:) * wfx_pnd_in(:,:) * z1_rhow * rdt_ice  ! Total available melt water, to be distributed as melt ponds
+      ! OLI: 07/2017 Bugfix above, removed " * aice(:,:)"
+      !js 19/04/18: change zrfrac to use aicen
+
+      zvolp(:,:) = vt_ip(:,:)
+
+      DO jl = 1, jpl
+         ! Melt water, to be distributed as melt ponds
+         zvolp(:,:) = zvolp(:,:) - zrfrac(:,:,jl)                                                  & 
+                                    * ( dh_i_pnd(:,:,jl)*rhoic + dh_s_pnd(:,:,jl)*rhosn )          & 
+                                    * z1_rhow * a_i(:,:,jl)
+      END DO
+
+! MV ---> use expression from level ice melt ponds (dv_mlt)
+
+      !js 03/05/19: we truncate negative values after calculating zvolp, in a
+      ! similar manner to the subroutine lim_mp_cesm. Variation dh_i_pnd and
+      ! dh_s_pnd are negative, indicating a loss of ice or snow. But we can expect them
+      ! to be negative for some reasons. We keep this behaviour as it is, for
+      ! fluxes conservation reasons. If some dh are positive, then we remove water
+      ! indirectly from the ponds.
+      zvolp(:,:) = MAX( zvolp(:,:) , 0._wp )
+
+
+      ! Fresh water flux going into the ponds
+      wfx_pnd_in(:,:) = wfx_pnd_in(:,:) + rhow * ( zvolp(:,:) - vt_ip(:,:) ) * r1_rdtice
+
+      !--- Remove retained meltwater from surface fluxes
+      IF ( ln_pnd_fw ) THEN
+         !wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrfrac(:,:) )
+         !wfx_sum(:,:)     = wfx_sum(:,:)     *  ( 1. - zrfrac(:,:) )
+
+         !js 190419: we change the code to use a_i in zrfrac. To be tested, but
+         !it should be conservative. zwfx_tmp is the flux accumulated in the
+         !ponds. wfx_pnd_in is the total surface melt fluxes.
+         zwfx_tmp(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:) , epsi10 )                           &
+            &            * MAX( 0._wp, SIGN( 1._wp, ABS(wfx_sum(:,:) + wfx_snw_sum(:,:)) - epsi10 ) )
+         WHERE ( ABS(zwfx_tmp(:,:)) > epsi10 )
+            zwfx_tmp(:,:) = wfx_pnd_in(:,:) / zwfx_tmp(:,:)
+         ELSEWHERE
+            zwfx_tmp(:,:) = 0._wp
+         ENDWHERE
+
+         wfx_sum(:,:)     = ( 1._wp - zwfx_tmp(:,:) ) * wfx_sum(:,:)
+         wfx_snw_sum(:,:) = ( 1._wp - zwfx_tmp(:,:) ) * wfx_snw_sum(:,:)
+      ENDIF
+
+       !-----------------------------------------------------------------
+       ! Identify grid cells with ponds
+       !-----------------------------------------------------------------
+
+       icells = 0
+       DO jj = 1, jpj
+         DO ji = 1, jpi
+            IF ( aice(ji,jj) > epsi10 ) THEN
+               zhi = vice(ji,jj) / aice(ji,jj)
+            ELSE
+               zhi = 0._wp
+            END IF
+
+            IF ( aice(ji,jj) > z01 .and. zhi > rn_himin .and. &
+                 zvolp(ji,jj) > zmin_volp*aice(ji,jj)) THEN
+               icells = icells + 1
+               indxi(icells) = ji
+               indxj(icells) = jj
+            ELSE  ! remove ponds on thin ice, or too small ponds
+               zvolpn(ji,jj,:) = 0._wp
+               zvolp (ji,jj) = 0._wp
+
+               a_ip(ji,jj,:)      = 0._wp
+               v_ip(ji,jj,:)      = 0._wp
+               a_ip_frac(ji,jj,:) = 0._wp
+               h_ip(ji,jj,:)      = 0._wp
+
+               vt_ip(ji,jj)       = 0._wp
+               at_ip(ji,jj)       = 0._wp
+
+!               IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean
+!                   wfx_pnd_out(ji,jj)   = wfx_pnd_out(ji,jj) + zvolp(ji,jj) * rhow * r1_rdtice
+            END IF
+         END DO                     ! ji
+      END DO                     ! jj
+
+
+       DO ij = 1, icells
+       
+          ji = indxi(ij)
+          jj = indxj(ij)
+
+          !--------------------------------------------------------------
+          ! Shrink pond due to refreezing
+          !--------------------------------------------------------------
+          ! OLI 07/2017: Done like for empirical melt pond scheme (CESM).
+          ! Therefore, I chose to put this part of the code before the main
+          ! routines lim_mp_area/depth (contrary to the original code),
+          ! seeing the freeze-up as a global sink of
+          ! freshwater for melt ponds in the whole grid cell. If this was done
+          ! after, I would need to make an additional assumption on the shape of
+          ! melt ponds, which I don't want to do (for the CESM scheme, this
+          ! assumption was on the aspect ratio). So I remove some water due to
+          ! refreezing first (using zTavg instead of zTsfcn in each category) and
+          ! then let the FF07 routines do their job for the fractional areas and
+          ! depths of melt ponds.
+          ! The whole ice lid related stuff from FF07 was thus removed and replaced
+          ! by this. As mentionned below, this should be improved, but is much
+          ! easier to conserve heat and freshwater this way.
+
+          ! Average surface temperature is needed to compute freeze-up at the cell
+          ! scale
+          zTavg = 0._wp
+          DO jl = 1, jpl
+             zTavg = zTavg + zTsfcn(ji,jj,jl)*aicen(ji,jj,jl)
+          END DO
+          zTavg = zTavg / aice(ji,jj)
+
+          ! The freezing temperature for meltponds is assumed slightly below 0C,
+          ! as if meltponds had a little salt in them (hence the use of zTp).
+          ! The salt budget is not
+          ! altered for meltponds, but if it were then an actual pond freezing
+          ! temperature could be computed.
+
+          zdTs        = MAX ( zTp - zTavg, 0. )
+
+          zvpold      = zvolp(ji,jj)
+
+          zvolp(ji,jj)  = zvolp(ji,jj) * EXP( zrexp * zdTs / zTp )
+
+          !--- Dump meltwater due to refreezing ( of course this is wrong
+          !--- but this parameterization is too simple )
+!          IF ( ln_pnd_fw ) &
+!             wfx_pnd_out(ji,jj)   = wfx_pnd_out(ji,jj) + rhow * ( zvpold - zvolp(ji,jj) ) * r1_rdtice
+!             ! OLI 07/2017 : Bugfix above, zvpold - zvolp instead of the
+!             ! opposite, otherwise negative contribution
+
+          !--------------------------------------------------------------
+          ! calculate pond area and depth
+          !--------------------------------------------------------------
+          zdvn = 0._wp
+          CALL lim_mp_area(aice(ji,jj),vice(ji,jj), &
+                    aicen(ji,jj,:), vicen(ji,jj,:), vsnon(ji,jj,:), &
+                    ticen(ji,jj,:,:), salin(ji,jj,:,:), &
+                    zvolpn(ji,jj,:), zvolp(ji,jj), &
+                    zapondn(ji,jj,:),zhpondn(ji,jj,:), zdvn)
+          ! outputs are
+          ! - zdvn
+          ! - zvolpn
+          ! - zvolp
+          ! - zapondn
+          ! - zhpondn
+
+!          IF ( ln_pnd_fw ) &
+!             wfx_pnd_out(ji,jj) = wfx_pnd_out(ji,jj) + zdvn * rhow * r1_rdtice ! update flux from ponds to ocean
+
+          !---------------------------------------------------------------
+          ! Update pond volume and fraction
+          !---------------------------------------------------------------
+          DO jl = 1, jpl
+             a_ip(ji,jj,jl) = zapondn(ji,jj,jl)
+             v_ip(ji,jj,jl) = zvolpn(ji,jj,jl)
+             a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / MAX(aicen(ji,jj,jl), epsi10) &
+                                    * MAX(0._wp, SIGN(1._wp, aicen(ji,jj,jl) - epsi10))
+             h_ip     (ji,jj,jl) = zhpondn(ji,jj,jl)
+          END DO
+       END DO ! ij
+
+       IF ( ln_pnd_fw ) THEN
+          !js 15/05/19: water going out of the ponds give a positive freshwater
+          ! flux. 
+          wfx_pnd_out(:,:) = SUM(MAX(0._wp, v_ip_b(:,:,:) - v_ip(:,:,:)), DIM=3) * rhow * r1_rdtice
+       ELSE
+         wfx_pnd_out(:,:) = 0._wp
+       ENDIF
+
+   END SUBROUTINE pnd_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, &
+          rhos,  &  ! OLI 07/2017 : for now this is useless, but will be useful with new snow scheme
+          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
+
+      !-----------|
+      !           |
+      !           |-----------|
+      !___________|___________|______________________________________sea-level
+      !           |           |
+      !           |           |---^--------|
+      !           |           |   |        |
+      !           |           |   |        |-----------|              |-------
+      !           |           |   |alfan(n)|           |              |
+      !           |           |   |        |           |--------------|
+      !           |           |   |        |           |              |
+      !---------------------------v-------------------------------------------
+      !           |           |   ^        |           |              |
+      !           |           |   |        |           |--------------|
+      !           |           |   |betan(n)|           |              |
+      !           |           |   |        |-----------|              |-------
+      !           |           |   |        |
+      !           |           |---v------- |
+      !           |           |
+      !           |-----------|
+      !           |
+      !-----------|
+
+       !-------------------------------------------------------------------
+       ! initialize
+       !-------------------------------------------------------------------
+
+       rhos(:) = rhosn ! OLI 07/2017 : same has above
+
+       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
+             asnon(n) = z0          !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) = c1 ! 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)     = 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_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) > 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) &
+                   - rhos(ns)/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)), 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 lim_mp_depth(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 ?
+          zhpondn(n) = max((hpond - alfan(n) + alfan(1)), z0)      !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 = z0
+       DO n=1,jpl
+         msno = msno + vsnon(n) * rhos(n)
+       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
+             !IF (hicen(n) > z0) THEN           !js: from CICE 5.1.2
+                perm = 0._wp ! MV ugly dummy patch
+                CALL lim_mp_perm(ticen(:,n), salin(:,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 lim_mp_depth(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
+             !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(n)/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) = 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
+
+
+    SUBROUTINE lim_mp_depth(aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
+       !!-------------------------------------------------------------------
+       !!                ***  ROUTINE lim_mp_depth  ***
+       !!
+       !! ** Purpose :   Compute melt pond depth
+       !!-------------------------------------------------------------------
+
+       REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
+          aicen, &
+          asnon, &
+          hsnon, &
+          rhos,  &
+          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(n)/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 lim_mp_depth
+
+
+    SUBROUTINE lim_mp_perm(ticen, salin, 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,  &     ! 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
+       !-----------------------------------------------------------------
+
+       IF (maxval(Tin) <= -2.0_wp) THEN
+
+          ! Assur 1958
+          DO k = 1,nlay_i
+             Sbr = - 1.2_wp                 &
+                   - 21.8_wp    * Tin(k)    &
+                   - 0.919_wp   * Tin(k)**2 &
+                   - 0.01878_wp * Tin(k)**3
+             phi(k) = salin(k) / Sbr ! liquid fraction
+
+             !js: Sbr must not be zero
+         END DO ! k
+
+       ELSE
+
+          ! Notz 2005 thesis eq. 3.2
+          !js: "interstitial brine Sbr (in ppt) as a function of temperature T (in degC)".
+          DO k = 1,nlay_i
+             Sbr = - 17.6_wp   * Tin(k)    &
+                   - 0.389_wp  * Tin(k)**2 &
+                   - 0.00362_wp* Tin(k)**3
+             phi(k) = salin(k) / Sbr ! liquid fraction
+
+             !js: Sbr must not be zero
+          END DO
+
+       END IF
+
+       !-----------------------------------------------------------------
+       ! permeability
+       !-----------------------------------------------------------------
+
+       perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007)
+
+   END SUBROUTINE lim_mp_perm
+   
+   
+!----------------------------------------------------------------------------------------------------------------------
 
    SUBROUTINE ice_thd_pnd_init 
@@ -319,4 +1248,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 +1266,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
