Index: /NEMO/branches/2020/SI3-05_MP/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/2020/SI3-05_MP/cfgs/SHARED/namelist_ice_ref	(revision 13813)
+++ /NEMO/branches/2020/SI3-05_MP/cfgs/SHARED/namelist_ice_ref	(revision 13814)
@@ -195,4 +195,5 @@
 !------------------------------------------------------------------------------
    ln_pnd            = .true.         !  activate melt ponds or not
+      ln_pnd_TOPO    = .false.        !  constant  melt ponds
       ln_pnd_LEV     = .true.         !  level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012)
          rn_apnd_min =   0.15         !     minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ??
Index: /NEMO/branches/2020/SI3-05_MP/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/2020/SI3-05_MP/src/ICE/ice.F90	(revision 13813)
+++ /NEMO/branches/2020/SI3-05_MP/src/ICE/ice.F90	(revision 13814)
@@ -208,5 +208,5 @@
    !                                     !!** 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_LEV       !: Simple "level-ice" melt pond scheme, hybrid CESM (Holland et al 2012) + lids and drainage
    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
@@ -214,4 +214,5 @@
    REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1)
    REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1)
+   LOGICAL , PUBLIC ::   ln_pnd_TOPO      !: Topographic melt pond scheme (Flocco and Feltham, 2007)
    LOGICAL,  PUBLIC ::   ln_pnd_lids      !: Allow ponds to have frozen lids
    LOGICAL , PUBLIC ::   ln_pnd_alb       !: melt ponds affect albedo
Index: /NEMO/branches/2020/SI3-05_MP/src/ICE/icethd_pnd.F90
===================================================================
--- /NEMO/branches/2020/SI3-05_MP/src/ICE/icethd_pnd.F90	(revision 13813)
+++ /NEMO/branches/2020/SI3-05_MP/src/ICE/icethd_pnd.F90	(revision 13814)
@@ -69,9 +69,15 @@
       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  ==!
+      CASE (np_pndTOPO)  ;   CALL pnd_TOPO                  &          !==  Topographic melt ponds  ==!
+                          (   at_i,      a_i,               &
+                  &           vt_i,      v_i,               &
+                  &           v_s,                          &
+                  &           t_i,       s_i,               &
+                  &           a_ip_frac, h_ip,              &
+                  &                      t_su )   
          !
       END SELECT
@@ -419,5 +425,5 @@
        REAL (wp), PARAMETER :: &
           ! MV ouate de phoque!!! constants hard coded ????
-          ! 917 = rhoice
+          ! 917 = rhoi
           ! 0.334 =Lfus!!!!
           zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3)
@@ -440,9 +446,9 @@
       zTsfcn(:,:,:) = Tsfc(:,:,:) - rt0         ! Convert in Celsius
 
-      IF ( ln_pnd_fw ) THEN
-         v_ip_b(:,:,:) = v_ip(:,:,:)
-      ELSE
-         v_ip_b(:,:,:) = 0._wp
-      ENDIF
+!     IF ( ln_pnd_fw ) THEN
+!        v_ip_b(:,:,:) = v_ip(:,:,:)
+!     ELSE
+!        v_ip_b(:,:,:) = 0._wp
+!     ENDIF
 
       !------------------------------------------------------------------
@@ -490,6 +496,6 @@
       !!! zrfrac(:,:,:) = rn_pnd_fracmin + ( rn_pnd_fracmax - rn_pnd_fracmin ) * aicen(:,:,:)
       ! MV2020
-      zrfrac(:,:,:) = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(:) !  = ( 1 - r ) = fraction of melt water that is not flushed
-      zwfx_tmp(:,:) = 0._wp
+      !!!jzrfrac(:,:,:) = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(:) !  = ( 1 - r ) = fraction of melt water that is not flushed
+      !!!zwfx_tmp(:,:) = 0._wp
 
       !--- Add retained melt water to melt ponds
@@ -507,7 +513,7 @@
       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)
+!        zvolp(:,:) = zvolp(:,:) - zrfrac(:,:,jl)                                                  & 
+!                                   * ( dh_i_pnd(:,:,jl)*rhoi + dh_s_pnd(:,:,jl)*rhosn )          & 
+!                                   * z1_rhow * a_i(:,:,jl)
 ! MV ---> use expression from level ice melt ponds (dv_mlt)
                                             
@@ -530,9 +536,9 @@
 
 
-      ! Fresh water flux going into the ponds
-      wfx_pnd_in(:,:) = wfx_pnd_in(:,:) + rhow * ( zvolp(:,:) - vt_ip(:,:) ) * r1_rdtice
+!     ! 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
+!     IF ( ln_pnd_fw ) THEN
          !wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrfrac(:,:) )
          !wfx_sum(:,:)     = wfx_sum(:,:)     *  ( 1. - zrfrac(:,:) )
@@ -543,13 +549,13 @@
          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
+!        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
+!     ENDIF
       
       
@@ -676,11 +682,11 @@
        END DO ! ij
 
-       IF ( ln_pnd_fw ) THEN
+!      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
+!         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
@@ -728,5 +734,4 @@
           hsnon, &
           asnon, &
-          rhos,  &  ! OLI 07/2017 : for now this is useless, but will be useful with new snow scheme
           alfan, &
           betan, &
@@ -781,6 +786,4 @@
        !-------------------------------------------------------------------
 
-       rhos(:) = rhosn ! OLI 07/2017 : same has above
-
        DO n = 1, jpl
 
@@ -867,5 +870,5 @@
              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
+                   - 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)), z0)
@@ -897,5 +900,5 @@
        ! height and area corresponding to the remaining volume
 
-      CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
+      CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index)
 
        DO n=1, m_index
@@ -917,7 +920,7 @@
        msno = z0
        DO n=1,jpl
-         msno = msno + vsnon(n) * rhos(n)
+         msno = msno + vsnon(n) * rhos
        END DO
-       floe_weight = (msno + rhoic*vice + rau0*zvolp) / aice
+       floe_weight = (msno + rhoi*vice + rau0*zvolp) / aice
        hsl_rel = floe_weight / rau0 &
                - ((sum(betan(:)*aicen(:))/aice) + alfan(1))
@@ -950,5 +953,5 @@
           IF (permflag > 0) THEN
              ! recompute pond depth
-            CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
+            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)
@@ -970,5 +973,5 @@
        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))
+                   - (rhos/rhow) * asnon(n) * min(hsnon(n), zhpondn(n))
        END DO
 
@@ -1015,5 +1018,5 @@
 
 
-    SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, rhos, alfan, zvolp, cum_max_vol, hpond, m_index)
+    SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index)
        !!-------------------------------------------------------------------
        !!                ***  ROUTINE ice_thd_pnd_depth  ***
@@ -1026,5 +1029,4 @@
           asnon, &
           hsnon, &
-          rhos,  &
           alfan, &
           cum_max_vol
@@ -1161,5 +1163,5 @@
 
            area = sum(aicetl(:)) - &                 ! total area of sub-layer
-                (rhos(n)/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
+                (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
