Index: /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/cfgs/SHARED/namelist_ice_ref	(revision 12499)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/cfgs/SHARED/namelist_ice_ref	(revision 12500)
@@ -177,4 +177,10 @@
    ln_pnd           = .false.         !  activate melt ponds or not
      ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012)
+       rn_pnd_min   = 0.15            !  Minimum ice fraction that contributes to melt ponds
+       rn_pnd_max   = 0.7             !  Maximum ice fraction that contributes to melt ponds
+       ln_pnd_totfrac = .false.       !  Use total ice fraction instead of category ice fraction in melt pond contribution fracton
+       ln_pnd_overflow = .false.      !  Allow ponds to overflow and have vertical flushing
+       ln_pnd_lids  = .false.         !  Melt ponds can have frozen lids
+       ln_use_pndmass = .true.        !  Use melt pond mass flux diagnostic, passing it to the ocean for emp
      ln_pnd_CST     = .false.         !  activate constant  melt ponds
        rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC
Index: /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/doc/namelists/namthd_pnd
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/doc/namelists/namthd_pnd	(revision 12499)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/doc/namelists/namthd_pnd	(revision 12500)
@@ -4,4 +4,10 @@
    ln_pnd           = .false.         !  activate melt ponds or not
      ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012)
+       rn_pnd_min   = 0.15            !  Minimum ice fraction that contributes to melt ponds
+       rn_pnd_max   = 0.7             !  Maximum ice fraction that contributes to melt ponds
+       ln_pnd_totfrac = .false.       !  Use total ice fraction instead of category ice fraction in melt pond contribution fracton
+       ln_pnd_overflow = .false.      !  Allow ponds to overflow and have vertical flushing
+       ln_pnd_lids  = .false.         !  Melt ponds can have frozen lids
+       ln_use_pndmass = .true.        !  Use melt pond mass flux diagnostic, passing it to the ocean for emp
      ln_pnd_CST     = .false.         !  activate constant  melt ponds
        rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC
Index: /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/ice.F90	(revision 12499)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/ice.F90	(revision 12500)
@@ -196,4 +196,11 @@
    REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1)
    LOGICAL , PUBLIC ::   ln_pnd_alb       !: melt ponds affect albedo
+   REAL(wp), PUBLIC ::   rn_pnd_min       !: Minimum ice fraction that contributes to melt ponds
+   REAL(wp), PUBLIC ::   rn_pnd_max       !: Maximum ice fraction that contributes to melt ponds
+   LOGICAL,  PUBLIC ::   ln_pnd_totfrac   !: Use total ice fraction instead of category ice fraction
+                                          !: when determining ice fraction that contributes to melt ponds
+   LOGICAL,  PUBLIC ::   ln_pnd_overflow  !: Allow ponds to overflow and have vertical flushing
+   LOGICAL,  PUBLIC ::   ln_pnd_lids      !: Melt ponds can have frozen lids
+   LOGICAL,  PUBLIC ::   ln_use_pndmass   !: Use melt pond mass flux diagnostic, passing it to the ocean for emp
 
    !                                     !!** ice-diagnostics namelist (namdia) **
Index: /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/ice1d.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/ice1d.F90	(revision 12499)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/ice1d.F90	(revision 12500)
@@ -111,6 +111,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_tot      !: Snow accretion/ablation        [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_sum      !: Ice surface ablation [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_pnd      !: Ice surface adlation into ponds [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_s_pnd      !: Snow surface adlation into ponds [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_itm      !: Ice internal ablation [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bom      !: Ice bottom ablation  [m]
@@ -131,5 +129,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_frac_1d  !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   lh_ip_1d      !: Ice pond lid thickness   [m]
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_pnd_avail_1d !: Fraction of sea ice available for melt ponding
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s
@@ -214,6 +211,5 @@
          &      dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d  (jpij) , s_i_new (jpij) ,  &
          &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) , lh_ip_1d(jpij) ,  &
-         &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) , dh_i_pnd(jpij) , a_pnd_avail_1d(jpij) ,           &
-         &      dh_s_pnd(jpij) ,                                                                        &
+         &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) ,                                                   &
          &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) )
       !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icerst.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icerst.F90	(revision 12499)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icerst.F90	(revision 12500)
@@ -172,5 +172,5 @@
       INTEGER           ::   jk
       LOGICAL           ::   llok
-      INTEGER           ::   id0, id1, id2, id3, id4   ! local integer
+      INTEGER           ::   id0, id1, id2, id3, id4, id5   ! local integer
       CHARACTER(len=25) ::   znam
       CHARACTER(len=2)  ::   zchar, zchar1
@@ -246,9 +246,14 @@
             CALL iom_get( numrir, jpdom_autoglo, 'a_ip' , a_ip )
             CALL iom_get( numrir, jpdom_autoglo, 'v_ip' , v_ip )
-            CALL iom_get( numrir, jpdom_autoglo, 'lh_ip', lh_ip)
          ELSE                                     ! start from rest
             IF(lwp) WRITE(numout,*) '   ==>>   previous run without melt ponds output then set it to zero'
             a_ip(:,:,:) = 0._wp
             v_ip(:,:,:) = 0._wp
+         ENDIF
+         ! melt pond lids
+         id5 = iom_varid( numrir, 'lh_ip' , ldstop = .FALSE. )
+         IF( id5 > 0 ) THEN
+            CALL iom_get( numrir, jpdom_autoglo, 'lh_ip', lh_ip)
+         ELSE
             lh_ip(:,:,:) = 0._wp
          ENDIF
Index: /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd.F90	(revision 12499)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd.F90	(revision 12500)
@@ -220,6 +220,4 @@
             dh_i_sub  (1:npti) = 0._wp ; dh_i_bog(1:npti) = 0._wp
             dh_snowice(1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp
-            dh_i_pnd  (1:npti) = 0._wp
-            dh_s_pnd  (1:npti) = 0._wp
             !                                      
                               CALL ice_thd_zdf                      ! --- Ice-Snow temperature --- !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_dh.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_dh.F90	(revision 12499)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_dh.F90	(revision 12500)
@@ -24,5 +24,4 @@
    USE lib_mpp        ! MPP library
    USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
-   USE icethd_pnd, ONLY : a_pnd_avail
    
    IMPLICIT NONE
@@ -99,5 +98,4 @@
       REAL(wp), DIMENSION(jpij,nlay_i) ::   zh_i      ! ice layer thickness
       REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah
-      REAL(wp), DIMENSION(jpij,nlay_i) ::   zdeltah_pnd  ! Change in thickness due to melting into ponds
       INTEGER , DIMENSION(jpij,nlay_i) ::   icount    ! number of layers vanished by melting 
 
@@ -114,9 +112,4 @@
          CASE( 2 )       ;   zswitch_sal = 1._wp   ! varying salinity profile
       END SELECT
-
-      ! Initialise fraction of sea ice available for melt ponding
-      DO ji = 1, npti
-         a_pnd_avail_1d(ji) = a_pnd_avail * at_i_1d(ji)
-      END DO
 
       ! initialize layer thicknesses and enthalpies
@@ -249,5 +242,4 @@
                zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) )                   ! bound melting
                zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)
-               zdeltah_pnd(ji,jk) = zdeltah  (ji,jk) * a_pnd_avail_1d(ji)
                
                hfx_snw_1d(ji)     = hfx_snw_1d(ji)     - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_rdtice   ! heat used to melt snow(W.m-2, >0)
@@ -256,5 +248,4 @@
                ! updates available heat + thickness
                dh_s_mlt(ji)    = dh_s_mlt(ji) + zdeltah(ji,jk)
-               dh_s_pnd(ji)    = dh_s_pnd(ji) + zdeltah_pnd(ji,jk)                ! Cumulate surface melt on areas contributing to melt ponds
                zq_top  (ji)    = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
                h_s_1d  (ji)    = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
@@ -352,5 +343,4 @@
                
                zdeltah(ji,jk) = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0]
-               zdeltah_pnd(ji,jk) = - zfmdt * r1_rhoi * a_pnd_avail_1d(ji)                          ! Melt of layer jk [m, <0] into meltponds
                
                zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0]
@@ -359,5 +349,4 @@
                
                dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)         ! Cumulate surface melt
-               dh_i_pnd(ji)   = dh_i_pnd(ji) + zdeltah_pnd(ji,jk)      ! Cumulate surface melt on areas contributing to melt ponds
                
                zfmdt          = - rhoi * zdeltah(ji,jk)               ! Recompute mass flux [kg/m2, >0]
Index: /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_pnd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_pnd.F90	(revision 12499)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_pnd.F90	(revision 12500)
@@ -38,5 +38,4 @@
    INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012)
 
-   REAL(wp), PUBLIC, PARAMETER :: a_pnd_avail = 0.7_wp   ! Fraction of sea ice available for melt ponding
    REAL(wp), PUBLIC, PARAMETER :: pnd_lid_max = 0.015_wp !  pond lid thickness above which the ponds disappear from the albedo calculation
    REAL(wp), PUBLIC, PARAMETER :: pnd_lid_min = 0.005_wp !  pond lid thickness below which the full pond area is used in the albedo calculation
@@ -137,4 +136,5 @@
       REAL(wp), PARAMETER ::   max_h_diff_s = -1.0E-6  ! Maximum meltpond depth change due to leaking or overflow (m s-1)
       !
+      REAL(wp) ::   tot_mlt          ! Total ice and snow surface melt (some goes into ponds, some into the ocean)
       REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding
       REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding (equivalent volume change)
@@ -151,5 +151,4 @@
       REAL(wp) ::   dh_ip_over       ! Pond thickness change due to overflow or leaking
       REAL(wp) ::   dh_i_pndleak     ! Grid box mean change in water depth due to leaking back to the ocean
-      REAL(wp) ::   weighted_lid_snow ! Lid to go on ponds under snow if snow thickness exceeds snow_lid_min
       REAL(wp) ::   h_gravity_head   ! Height above sea level of the top of the melt pond
       REAL(wp) ::   h_percolation    ! Distance between the base of the melt pond and the base of the sea ice
@@ -182,7 +181,4 @@
             h_ip_1d(ji)      = 0._wp
             lh_ip_1d(ji)     = 0._wp
-
-            zdv_mlt = 0._wp
-            zdh_frz = 0._wp
             !                                                     !--------------------------------!
          ELSE                                                     ! Case ice thickness >= rn_himin !
@@ -194,7 +190,12 @@
             IF ( za_ip < epsi06 ) za_ip = epsi06
             !
-            ! available meltwater for melt ponding
-            ! This is the change in ice volume due to melt
-            zdv_mlt = -(dh_i_pnd(ji)*rhoi + dh_s_pnd(ji)*rhos) * z1_rhow * a_i_1d(ji)
+            ! available meltwater for melt ponding [m, >0] and fraction 
+            tot_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow 
+            IF ( ln_pnd_totfrac ) THEN
+               zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * at_i_1d(ji) ! Use total ice fraction
+            ELSE
+               zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * a_i_1d(ji)  ! Use category ice fraction 
+            ENDIF 
+            zdv_mlt = zfr_mlt * tot_mlt 
             !
             !--- Pond gowth ---!
@@ -202,30 +203,48 @@
             !
             !--- Lid shrinking. ---!
-            lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip
-            !
+            IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip
+            !
+            ! melt pond mass flux (<0)
+            IF( ln_use_pndmass .AND. zdv_mlt > 0._wp ) THEN
+               zfac = zdv_mlt * rhow * r1_rdtice
+               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
+               !
+               ! adjust ice/snow melting flux to balance melt pond flux (>0)
+               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )
+               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
             !
             !--- Pond contraction (due to refreezing) ---!
-            IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN
-               t_grad = (zTp+rt0) - t_su_1d(ji)
-               
-               ! The following equation is a rearranged form of:
-               ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow)
-               ! where: lid_thickness_start = lh_ip_1d(ji)
-               !        lid_thickness_end = lh_ip_end
-               !        omega_dt is a bunch of terms in the equation that do not change
-               ! note the use of rhow instead of rhoi as we are working with volumes and it is easier if the water and ice volumes (for the lid and the pond) are the same
-               ! (have the same density). The lid will eventually remelt anyway so it doesn't matter if we make this simplification.
-                             
-               lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2)
-               zdh_frz = lh_ip_end - lh_ip_1d(ji)
-
-               ! Pond shrinking
-               v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji)
-
-               ! Lid growing
-               lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz
+            IF ( ln_pnd_lids ) THEN
+
+               ! Code to use if using melt pond lids
+               IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN
+                  t_grad = (zTp+rt0) - t_su_1d(ji)
+
+                  ! The following equation is a rearranged form of:
+                  ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow)
+                  ! where: lid_thickness_start = lh_ip_1d(ji)
+                  !        lid_thickness_end = lh_ip_end
+                  !        omega_dt is a bunch of terms in the equation that do not change
+                  ! Note the use of rhow instead of rhoi as we are working with volumes and it is mathematically easier
+                  ! if the water and ice specific volumes (for the lid and the pond) are the same (have the same density). 
+
+                  lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2)
+                  zdh_frz = lh_ip_end - lh_ip_1d(ji)
+
+                  ! Pond shrinking
+                  v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji)
+
+                  ! Lid growing
+                  F ( ln_pnd_lids )  lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz
+               ELSE
+                  zdh_frz = 0._wp
+               END IF
+
             ELSE
-               zdh_frz = 0._wp
-            END IF
+               ! Code to use if not using melt pond lids
+               v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp )
+            ENDIF
 
             !
@@ -239,69 +258,72 @@
             a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
             h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
-
-            ! If pond area exceeds a_pnd_avail_1d(ji) * a_i_1d(ji) then reduce the pond volume
-            IF ( a_ip_1d(ji) > a_pnd_avail_1d(ji) * a_i_1d(ji) ) THEN
-                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
-                dh_ip_over = zpnd_aspect * a_pnd_avail_1d(ji) - h_ip_1d(ji)   ! This will be a negative number
-                dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit
-                h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
-                a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
-                a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
-                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
-
+            
+            !--- Pond overflow and vertical flushing ---!
+            IF ( ln_pnd_overflow ) THEN
+
+               ! If pond area exceeds a_pnd_avail_1d(ji) * a_i_1d(ji) then reduce the pond volume
+               IF ( a_ip_1d(ji) > a_pnd_avail_1d(ji) * a_i_1d(ji) ) THEN
+                   v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
+                   dh_ip_over = zpnd_aspect * a_pnd_avail_1d(ji) - h_ip_1d(ji)   ! This will be a negative number
+                   dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit
+                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
+                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
+                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
+                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
+               ENDIF
+
+               ! If pond depth exceeds half the ice thickness then reduce the pond volume
+               IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN
+                   v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
+                   dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji)                ! This will be a negative number
+                   dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit
+                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
+                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
+                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
+                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
+               ENDIF
+
+               !-- Vertical flushing of pond water --!
+               ! The height above sea level of the top of the melt pond is the ratios of density of ice and water times the ice depth.
+               ! This assumes the pond is sitting on top of the ice.
+               h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow
+
+               ! The depth through which water percolates is the distance under the melt pond
+               h_percolation = h_i_1d(ji) - h_ip_1d(ji)
+
+               ! Calculate the permeability of the ice (Assur 1958)
+               DO jk = 1, nlay_i
+                   Sbr = - 1.2_wp                         &
+                         - 21.8_wp     * (t_i_1d(ji,jk)-rt0)    &
+                         - 0.919_wp    * (t_i_1d(ji,jk)-rt0)**2 &
+                         - 0.01878_wp  * (t_i_1d(ji,jk)-rt0)**3
+                   phi(jk) = sz_i_1d(ji,jk)/Sbr
+               END DO
+               perm = 3.0e-08_wp * (minval(phi))**3
+
+               ! Do the drainage using Darcy's law
+               IF ( perm > 0._wp ) THEN
+                   dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation)  ! This should be a negative number
+                   dh_ip_over = MIN(dh_ip_over, 0._wp)      ! Make sure it is negative
+
+                   v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
+                   dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit
+                   h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
+                   a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
+                   a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
+                   v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
+               ENDIF
             ENDIF
 
-            ! If pond depth exceeds half the ice thickness then reduce the pond volume
-            IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN
-                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
-                dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji)                ! This will be a negative number
-                dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit
-                h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
-                a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
-                a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
-                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
-
+            ! If lid thickness is ten times greater than pond thickness then remove pond 
+            IF ( ln_pnd_lids ) THEN           
+               IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
+                   a_ip_1d(ji)      = 0._wp
+                   a_ip_frac_1d(ji) = 0._wp
+                   h_ip_1d(ji)      = 0._wp
+                   lh_ip_1d(ji)     = 0._wp
+                   v_ip_1d(ji)      = 0._wp
+               ENDIF
             ENDIF
-
-            !-- Vertical flushing of pond water --!
-            ! The height above sea level of the top of the melt pond is the ratios of density of ice and water times the ice depth.
-            ! This assumes the pond is sitting on top of the ice.
-            h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow
-
-            ! The depth through which water percolates is the distance under the melt pond
-            h_percolation = h_i_1d(ji) - h_ip_1d(ji)
-
-            ! Calculate the permeability of the ice (Assur 1958)
-            DO jk = 1, nlay_i
-                Sbr = - 1.2_wp                         &
-                      - 21.8_wp     * (t_i_1d(ji,jk)-rt0)    &
-                      - 0.919_wp    * (t_i_1d(ji,jk)-rt0)**2 &
-                      - 0.01878_wp  * (t_i_1d(ji,jk)-rt0)**3
-                phi(jk) = sz_i_1d(ji,jk)/Sbr
-            END DO
-            perm = 3.0e-08_wp * (minval(phi))**3
-
-            ! Do the drainage using Darcy's law
-            IF ( perm > 0._wp ) THEN
-                dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation)  ! This should be a negative number
-                dh_ip_over = MIN(dh_ip_over, 0._wp)      ! Make sure it is negative
-
-                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
-                dh_ip_over = MAX(dh_ip_over,max_h_diff_ts)                       ! Apply a limit
-                h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)
-                a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect
-                a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)
-                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
-
-            ENDIF
-
-            ! If lid thickness is ten times greater than pond thickness then remove pond            
-            IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
-                a_ip_1d(ji)      = 0._wp
-                a_ip_frac_1d(ji) = 0._wp
-                h_ip_1d(ji)      = 0._wp
-                lh_ip_1d(ji)     = 0._wp
-                v_ip_1d(ji)      = 0._wp
-            END IF
 
             ! If any of the previous changes has removed all the ice thickness then remove ice area.
@@ -333,5 +355,7 @@
       INTEGER  ::   ios, ioptio   ! Local integer
       !!
-      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
+      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb,          &
+                            rn_pnd_min, rn_pnd_max, ln_pnd_overflow, ln_pnd_lids, ln_pnd_totfrac,  &
+                            ln_use_pndmass
       !!-------------------------------------------------------------------
       !
@@ -351,4 +375,10 @@
          WRITE(numout,*) '      Melt ponds activated or not                                     ln_pnd     = ', ln_pnd
          WRITE(numout,*) '         Evolutive  melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12
+         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds       rn_pnd_min = ', rn_pnd_min
+         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds       rn_pnd_max = ', rn_pnd_max
+         WRITE(numout,*) '            Use total ice fraction instead of category ice fraction   ln_pnd_totfrac  = ',ln_pnd_totfrac
+         WRITE(numout,*) '            Allow ponds to overflow and have vertical flushing        ln_pnd_overflow = ',ln_pnd_overflow
+         WRITE(numout,*) '            Melt ponds can have frozen lids                           ln_pnd_lids     = ',ln_pnd_lids
+         WRITE(numout,*) '            Use melt pond mass flux diagnostic, passing to ocean      ln_use_pndmass  = ',ln_use_pndmass
          WRITE(numout,*) '         Prescribed melt pond fraction and depth                      ln_pnd_CST = ', ln_pnd_CST
          WRITE(numout,*) '            Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
