Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/ice.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/ice.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/ice.F90	(revision 12636)
@@ -70,4 +70,5 @@
    !! a_ip        |      -      |    Ice pond concentration       |       |
    !! v_ip        |      -      |    Ice pond volume per unit area| m     |
+   !! lh_ip       !    lh_ip_1d !    Ice pond lid thickness       ! m     !
    !!                                                                     |
    !!-------------|-------------|---------------------------------|-------|
@@ -195,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) **
@@ -331,5 +339,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_ip       !: melt pond volume per grid cell area      [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_frac  !: melt pond fraction (a_ip/a_i)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_eff   !: melt pond effective fraction (not covered up by lid) (a_ip/a_i)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_ip       !: melt pond depth                          [m]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   lh_ip      !: melt pond lid thickness                  [m]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond concentration
@@ -448,5 +458,6 @@
 
       ii = ii + 1
-      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl) , STAT = ierr(ii) )
+      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl),  &
+         &      lh_ip(jpi,jpj,jpl), a_ip_eff(jpi,jpj,jpl) , STAT = ierr(ii) )
 
       ii = ii + 1
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/ice1d.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/ice1d.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/ice1d.F90	(revision 12636)
@@ -128,4 +128,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ip_1d       !:
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_frac_1d  !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_eff_1d   !:
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   lh_ip_1d      !: Ice pond lid thickness   [m]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s
@@ -157,4 +159,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_ip_2d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_ip_2d 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   lh_ip_2d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_su_2d 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_2d
@@ -208,6 +211,6 @@
          &      dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm  (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) ,  &    
          &      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) ,                   &
-         &      h_ip_1d (jpij) , a_ip_frac_1d(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) , a_ip_eff_1d(jpij) ,                               &
          &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) )
       !
@@ -226,5 +229,5 @@
       ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) ,  &
          &      v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) ,  &
-         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) ,                      &
+         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , lh_ip_2d(jpij,jpl),  &
          &      STAT=ierr(ii) )
 
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icealb.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icealb.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icealb.F90	(revision 12636)
@@ -96,4 +96,5 @@
       LOGICAL , INTENT(in   )                   ::   ld_pnd_alb   !  effect of melt ponds on albedo
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area)
+                                                                  !  This is the effective fraction not covered up by a pond lid
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth
       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_cs      !  albedo of ice under clear    sky
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn.F90	(revision 12636)
@@ -101,4 +101,5 @@
       ELSEWHERE
          h_ip(:,:,:) = 0._wp
+         lh_ip(:,:,:) = 0._wp
       END WHERE
       !
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv.F90	(revision 12636)
@@ -84,10 +84,10 @@
          !                             !-----------------------!
          CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, &
-            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
+            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, lh_ip, e_s, e_i )
          !                             !-----------------------!
       CASE( np_advPRA )                ! PRATHER scheme        !
          !                             !-----------------------!
          CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, &
-            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i )
+            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, lh_ip, e_s, e_i )
       END SELECT
 
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv_pra.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv_pra.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv_pra.F90	(revision 12636)
@@ -44,4 +44,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume
+   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxlhp, sylhp, sxxlhp, syylhp, sxylhp   ! melt pond lid thickness
 
    !! * Substitutions
@@ -55,5 +56,5 @@
 
    SUBROUTINE ice_dyn_adv_pra( kt, pu_ice, pv_ice,  &
-      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
       !!----------------------------------------------------------------------
       !!                **  routine ice_dyn_adv_pra  **
@@ -78,4 +79,5 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
@@ -88,5 +90,5 @@
       REAL(wp), DIMENSION(jpi,jpj,1)          ::   z0opw
       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi
-      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp
+      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp, z0lhp
       REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es
       REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei
@@ -129,4 +131,5 @@
             z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction
             z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume
+            z0lhp(:,:,jl)  = plh_ip(:,:,jl) * e1e2t(:,:)   ! Melt pond lid thickness
          ENDIF
       END DO
@@ -167,4 +170,6 @@
                CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
                CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
+               CALL adv_x( zdt , pu_ice , 1._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp)    !--- melt pond lid thickness
+               CALL adv_y( zdt , pv_ice , 0._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp)
             ENDIF
          END DO
@@ -202,4 +207,6 @@
                CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume
                CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )
+               CALL adv_y( zdt , pv_ice , 1._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp)    !--- melt pond lid thickness
+               CALL adv_x( zdt , pu_ice , 0._wp , zarea , z0lhp, sxlhp, sxxlhp, sylhp, syylhp, sxylhp) 
             ENDIF
          END DO
@@ -225,4 +232,5 @@
             pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
             pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
+            plh_ip(:,:,jl) = z0lhp (:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1)
          ENDIF
       END DO
@@ -231,5 +239,5 @@
       !     Remove negative values (conservation is ensured)
       !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
-      CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+      CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
       !
       ! --- Ensure snow load is not too big --- !
@@ -653,4 +661,5 @@
          &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   &
          &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   &
+         &      sxlhp(jpi,jpj,jpl) , sylhp (jpi,jpj,jpl), sxxlhp (jpi,jpj,jpl), syylhp (jpi,jpj,jpl), sxylhp (jpi,jpj,jpl),   &
          !
          &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , &
@@ -765,4 +774,10 @@
                CALL iom_get( numrir, jpdom_autoglo, 'syyvp', syyvp )
                CALL iom_get( numrir, jpdom_autoglo, 'sxyvp', sxyvp )
+               !                                                     ! melt pond lid thickness
+               CALL iom_get( numrir, jpdom_autoglo, 'sxlhp' , sxlhp  )
+               CALL iom_get( numrir, jpdom_autoglo, 'sylhp' , sylhp  )
+               CALL iom_get( numrir, jpdom_autoglo, 'sxxlhp', sxxlhp )
+               CALL iom_get( numrir, jpdom_autoglo, 'syylhp', syylhp )
+               CALL iom_get( numrir, jpdom_autoglo, 'sxylhp', sxylhp )
             ENDIF
             !
@@ -782,4 +797,5 @@
                sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction
                sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume
+               sxlhp  = 0._wp  ;   sylhp  = 0._wp  ;   sxxlhp  = 0._wp  ;   syylhp  = 0._wp  ;   sxylhp  = 0._wp  ! melt pond lid thickness
             ENDIF
          ENDIF
@@ -862,4 +878,10 @@
             CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp )
             CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp )
+            !                                                        ! melt pond lid thickness
+            CALL iom_rstput( iter, nitrst, numriw, 'sxlhp' , sxlhp  )
+            CALL iom_rstput( iter, nitrst, numriw, 'sylhp' , sylhp  )
+            CALL iom_rstput( iter, nitrst, numriw, 'sxxlhp', sxxlhp )
+            CALL iom_rstput( iter, nitrst, numriw, 'syylhp', syylhp )
+            CALL iom_rstput( iter, nitrst, numriw, 'sxylhp', sxylhp )
          ENDIF
          !
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv_umx.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv_umx.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_adv_umx.F90	(revision 12636)
@@ -60,5 +60,5 @@
 
    SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  &
-      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
       !!----------------------------------------------------------------------
       !!                  ***  ROUTINE ice_dyn_adv_umx  ***
@@ -85,4 +85,5 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond concentration
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
@@ -334,4 +335,10 @@
             CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
                &                                      zhvar, pv_ip, zua_ups, zva_ups )
+            ! lid thickness
+            zamsk = 0._wp
+            zhvar(:,:,:) = plh_ip(:,:,:) * z1_aip(:,:,:)
+            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, &
+               &                                      zhvar, plh_ip, zua_ups, zva_ups )
+            
          ENDIF
          !
@@ -350,5 +357,5 @@
          ! Remove negative values (conservation is ensured)
          !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20)
-         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
          !
          ! Make sure ice thickness is not too big
@@ -1519,5 +1526,5 @@
       !!
       !! ** Purpose : Thickness correction in case advection scheme creates
-      !!              abnormally tick ice or snow
+      !!              abnormally thick ice or snow
       !!
       !! ** Method  : 1- check whether ice thickness is larger than the surrounding 9-points
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_rdgrft.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_rdgrft.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icedyn_rdgrft.F90	(revision 12636)
@@ -503,6 +503,6 @@
       REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1
       REAL(wp)                  ::   airft1, oirft1, aprft1
-      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg  ! area etc of new ridges
-      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft  ! area etc of rafted ice
+      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, lhprdg  ! area etc of new ridges
+      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, lhprft  ! area etc of rafted ice
       !
       REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges
@@ -578,7 +578,9 @@
                   aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1)
                   vprdg (ji) = v_ip_2d(ji,jl1) * afrdg
+                  lhprdg(ji) = lh_ip_2d(ji,jl1) * afrdg
                   aprft1     = a_ip_2d(ji,jl1) * afrft
                   aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft
                   vprft (ji) = v_ip_2d(ji,jl1) * afrft
+                  lhprft(ji) = lh_ip_2d(ji,jl1) * afrft
                ENDIF
 
@@ -610,4 +612,5 @@
                   a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1
                   v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji)
+                  lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - lhprdg(ji) - lhprft(ji)
                ENDIF
             ENDIF
@@ -706,4 +709,6 @@
                      a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         & 
                         &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   )
+                     lh_ip_2d (ji,jl2) = lh_ip_2d(ji,jl2) + (   lhprdg (ji) * rn_fpndrdg * fvol   (ji)   &
+                        &                                   + lhprft (ji) * rn_fpndrft * zswitch(ji)   )
                   ENDIF
                   
@@ -736,5 +741,5 @@
       !----------------
       ! In case ridging/rafting lead to very small negative values (sometimes it happens)
-      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )
+      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, lh_ip_2d, ze_s_2d, ze_i_2d )
       !
    END SUBROUTINE rdgrft_shift
@@ -848,4 +853,5 @@
          CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
          CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
+         CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) )
          DO jl = 1, jpl
             DO jk = 1, nlay_s
@@ -874,4 +880,5 @@
          CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) )
          CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) )
+         CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip(:,:,:) )
          DO jl = 1, jpl
             DO jk = 1, nlay_s
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceistate.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceistate.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceistate.F90	(revision 12636)
@@ -155,5 +155,7 @@
       a_ip     (:,:,:) = 0._wp
       v_ip     (:,:,:) = 0._wp
+      lh_ip    (:,:,:) = 0._wp
       a_ip_frac(:,:,:) = 0._wp
+      a_ip_eff (:,:,:) = 0._wp
       h_ip     (:,:,:) = 0._wp
       !
@@ -362,4 +364,5 @@
             a_ip_frac(:,:,:) = 0._wp
          END WHERE
+         a_ip_eff(:,:,:) = a_ip_frac(:,:,:)
          v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:)
           
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceitd.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceitd.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceitd.F90	(revision 12636)
@@ -411,4 +411,5 @@
       CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
       CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
+      CALL tab_3d_2d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip )
       CALL tab_3d_2d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
       DO jl = 1, jpl
@@ -483,4 +484,8 @@
                   v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans
                   v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans
+                  !                                              
+                  ztrans          = lh_ip_2d(ji,jl1) * zworka(ji)     ! Pond lid thickness
+                  lh_ip_2d(ji,jl1) = lh_ip_2d(ji,jl1) - ztrans
+                  lh_ip_2d(ji,jl2) = lh_ip_2d(ji,jl2) + ztrans
                ENDIF
                !
@@ -527,5 +532,5 @@
       ! clem: The transfer between one category to another can lead to very small negative values (-1.e-20)
       !       because of truncation error ( i.e. 1. - 1. /= 0 )
-      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d )
+      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, lh_ip_2d, ze_s_2d, ze_i_2d )
 
       ! at_i must be <= rn_amax
@@ -555,4 +560,5 @@
       CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip )
       CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip )
+      CALL tab_2d_3d( npti, nptidx(1:npti), lh_ip_2d(1:npti,1:jpl), lh_ip )
       CALL tab_2d_3d( npti, nptidx(1:npti), t_su_2d(1:npti,1:jpl), t_su )
       DO jl = 1, jpl
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icerst.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icerst.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icerst.F90	(revision 12636)
@@ -132,4 +132,5 @@
       CALL iom_rstput( iter, nitrst, numriw, 'a_ip' , a_ip  )
       CALL iom_rstput( iter, nitrst, numriw, 'v_ip' , v_ip  )
+      CALL iom_rstput( iter, nitrst, numriw, 'lh_ip', lh_ip )
       ! Snow enthalpy
       DO jk = 1, nlay_s 
@@ -171,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
@@ -250,4 +251,12 @@
             v_ip(:,:,:) = 0._wp
          ENDIF
+         a_ip_eff(:,:,:) = a_ip(:,:,:)
+         ! 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
          ! fields needed for Met Office (Jules) coupling
          IF( ln_cpl ) THEN
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icesbc.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icesbc.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icesbc.F90	(revision 12636)
@@ -132,5 +132,5 @@
 
       ! --- cloud-sky and overcast-sky ice albedos --- !
-      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os )
+      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, zalb_cs, zalb_os )
 
       ! albedo depends on cloud fraction because of non-linear spectral effects
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd.F90	(revision 12636)
@@ -355,4 +355,6 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
          CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) )
+         CALL tab_2d_1d( npti, nptidx(1:npti), lh_ip_1d     (1:npti), lh_ip     (:,:,kl) )
          !
          CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d  (1:npti), qprec_ice            )
@@ -461,4 +463,6 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) )
          CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_eff_1d (1:npti), a_ip_eff (:,:,kl) )
+         CALL tab_1d_2d( npti, nptidx(1:npti), lh_ip_1d    (1:npti), lh_ip    (:,:,kl) )
          !
          CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni )
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd_pnd.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd_pnd.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd_pnd.F90	(revision 12636)
@@ -38,4 +38,6 @@
    INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012)
 
+   REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
@@ -89,8 +91,10 @@
          IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
             a_ip_frac_1d(ji) = rn_apnd
+            a_ip_eff_1d(ji)  = rn_apnd
             h_ip_1d(ji)      = rn_hpnd    
             a_ip_1d(ji)      = a_ip_frac_1d(ji) * a_i_1d(ji)
          ELSE
             a_ip_frac_1d(ji) = 0._wp
+            a_ip_eff_1d(ji)  = 0._wp
             h_ip_1d(ji)      = 0._wp    
             a_ip_1d(ji)      = 0._wp
@@ -129,43 +133,81 @@
       REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio
       REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature
-      !
+      REAL(wp), PARAMETER ::   max_h_diff_s = -1.0E-6 ! Maximum meltpond depth change due to leaking or overflow (m s-1)
+      REAL(wp), PARAMETER ::   pnd_lid_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation
+      REAL(wp), PARAMETER ::   pnd_lid_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation
+      !
+      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
+      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding (equivalent volume change)
       REAL(wp) ::   z1_Tp            ! inverse reference temperature
       REAL(wp) ::   z1_rhow          ! inverse freshwater density
       REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio
+      REAL(wp) ::   z1_rhoi          ! inverse ice density
       REAL(wp) ::   zfac, zdum
-      !
-      INTEGER  ::   ji   ! loop indices
+      REAL(wp) ::   t_grad           ! Temperature deficit for refreezing
+      REAL(wp) ::   omega_dt         ! Time independent accumulated variables used for freezing
+      REAL(wp) ::   lh_ip_end        ! Lid thickness at end of timestep (temporary variable)
+      REAL(wp) ::   zdh_frz          ! Amount of melt pond that freezes (m)
+      REAL(wp) ::   v_ip_old         ! Pond volume before leaking back to the ocean
+      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) ::   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
+      REAL(wp) ::   Sbr              ! Brine salinity
+      REAL(wp), DIMENSION(nlay_i) ::   phi     ! liquid fraction
+      REAL(wp) ::   perm             ! Permeability of the sea ice
+      REAL(wp) ::   za_ip            ! Temporary pond fraction
+      REAL(wp) ::   max_h_diff_ts    ! Maximum meltpond depth change due to leaking or overflow (m per ts)
+      REAL(wp) ::   lfrac_pnd        ! The fraction of the meltpond exposed (not inder a frozen lid)
+      
+      !
+      INTEGER  ::   ji, jk   ! loop indices
       !!-------------------------------------------------------------------
       z1_rhow        = 1._wp / rhow 
       z1_zpnd_aspect = 1._wp / zpnd_aspect
       z1_Tp          = 1._wp / zTp 
+      z1_rhoi        = 1._wp / rhoi
+      max_h_diff_ts  = max_h_diff_s * rdt_ice
+
+      ! Define time-independent field for use in refreezing
+      omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow)
 
       DO ji = 1, npti
-         !                                                        !--------------------------------!
-         IF( h_i_1d(ji) < rn_himin) THEN                          ! Case ice thickness < rn_himin  !
-            !                                                     !--------------------------------!
-            !--- Remove ponds on thin ice
+
+         !                                                            !----------------------------------------------------!
+         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
             a_ip_frac_1d(ji) = 0._wp
             h_ip_1d(ji)      = 0._wp
+            lh_ip_1d(ji)     = 0._wp
             !                                                     !--------------------------------!
          ELSE                                                     ! Case ice thickness >= rn_himin !
             !                                                     !--------------------------------!
             v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
-            !
-            ! available meltwater for melt ponding [m, >0] and fraction
-            zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji)
-            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc
-            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper 
+
+            ! To avoid divide by zero errors in some calculations we will use a temporary pond fraction variable that has a minimum value of epsi06
+            za_ip = a_ip_1d(ji)
+            IF ( za_ip < epsi06 ) za_ip = epsi06
+            !
+            ! 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 ---!
-            ! v_ip should never be negative, otherwise code crashes
-            v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt )
+            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
+            !
+            !--- Lid shrinking. ---!
+            IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip
             !
             ! melt pond mass flux (<0)
-            IF( zdv_mlt > 0._wp ) THEN
-               zfac = zfr_mlt * zdv_mlt * rhow * r1_rdtice
+            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
                !
@@ -177,5 +219,39 @@
             !
             !--- Pond contraction (due to refreezing) ---!
-            v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp )
+            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
+                  IF ( ln_pnd_lids )  lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz
+               ELSE
+                  zdh_frz = 0._wp
+               END IF
+
+            ELSE
+               ! 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
+
+            !
+            ! Make sure pond volume or lid thickness has not gone negative
+            IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp 
+            IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp
             !
             ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac
@@ -184,6 +260,96 @@
             a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
             h_ip_1d(ji)      = zpnd_aspect * a_ip_frac_1d(ji)
+            
+            !--- Pond overflow and vertical flushing ---!
+            IF ( ln_pnd_overflow ) THEN
+
+               ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume
+               IF ( a_ip_1d(ji) > zfr_mlt * 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 * zfr_mlt - 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 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
+
+            ! If any of the previous changes has removed all the ice thickness then remove ice area.
+            IF ( h_i_1d(ji) == 0._wp ) THEN
+                a_i_1d(ji) = 0._wp
+                h_s_1d(ji) = 0._wp
+            ENDIF
+
             !
          ENDIF
+
+         ! Calculate how much meltpond is exposed to the atmosphere (not under a frozen lid)        
+         IF ( lh_ip_1d(ji) < pnd_lid_min ) THEN       ! Pond lid is very thin. Expose all the pond.
+            lfrac_pnd = 1.0
+         ELSE
+            IF ( lh_ip_1d(ji) > pnd_lid_max ) THEN    ! Pond lid is very thick. Cover all the pond up with ice and nosw.
+              lfrac_pnd = 0.0
+            ELSE                                      ! Pond lid is of intermediate thickness. Expose part of the melt pond.
+              lfrac_pnd = ( lh_ip_1d(ji) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min)
+            ENDIF
+         ENDIF
+         
+         ! Calculate the melt pond effective area
+         a_ip_eff_1d(ji) = a_ip_frac_1d(ji) * lfrac_pnd
+
       END DO
       !
@@ -205,5 +371,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
       !!-------------------------------------------------------------------
       !
@@ -223,4 +391,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
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceupdate.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceupdate.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/iceupdate.F90	(revision 12636)
@@ -185,5 +185,5 @@
       ! Snow/ice albedo (only if sent to coupler, useless in forced mode)
       !------------------------------------------------------------------
-      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
+      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos
       !
       alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icevar.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icevar.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icevar.F90	(revision 12636)
@@ -533,4 +533,5 @@
                a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj)
                v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj)
+               lh_ip (ji,jj,jl) = lh_ip (ji,jj,jl) * zswitch(ji,jj)
                !
             END DO
@@ -555,5 +556,5 @@
 
 
-   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i )
+   SUBROUTINE ice_var_zapneg( pdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_var_zapneg ***
@@ -570,4 +571,5 @@
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
@@ -637,9 +639,10 @@
       WHERE( pv_ip (:,:,:) < 0._wp )   pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)
       !                                                        but it does not change conservation, so keep it this way is ok
+      WHERE( plh_ip (:,:,:) < 0._wp )   plh_ip (:,:,:) = 0._wp
       !
    END SUBROUTINE ice_var_zapneg
 
 
-   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, pe_s, pe_i )
+   SUBROUTINE ice_var_roundoff( pa_i, pv_i, pv_s, psv_i, poa_i, pa_ip, pv_ip, plh_ip, pe_s, pe_i )
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_var_roundoff ***
@@ -654,4 +657,5 @@
       REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction
       REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume
+      REAL(wp), DIMENSION(:,:)  , INTENT(inout) ::   plh_ip     ! melt pond lid thickness
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_s       ! snw heat content
       REAL(wp), DIMENSION(:,:,:), INTENT(inout) ::   pe_i       ! ice heat content
@@ -668,4 +672,5 @@
          WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0
          WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0
+         WHERE( plh_ip(1:npti,:) < 0._wp .AND. plh_ip(1:npti,:) > -epsi10 )    plh_ip(1:npti,:)   = 0._wp   ! lh_ip must be >= 0
       ENDIF
       !
Index: NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icewri.F90
===================================================================
--- NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icewri.F90	(revision 12635)
+++ NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icewri.F90	(revision 12636)
@@ -149,4 +149,5 @@
 
       ! --- category-dependent fields --- !
+      IF( iom_use('icelhpnd_cat') )   CALL iom_put( 'icelhpnd_cat', lh_ip          * zmsk00l                                   ) ! melt pond lid thickness for categories
       IF( iom_use('icemask_cat' ) )   CALL iom_put( 'icemask_cat' ,                  zmsk00l                                   ) ! ice mask 0%
       IF( iom_use('iceconc_cat' ) )   CALL iom_put( 'iceconc_cat' , a_i            * zmsk00l                                   ) ! area for categories
@@ -165,4 +166,7 @@
       IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac for categories
       IF( iom_use('icealb_cat'  ) )   CALL iom_put( 'icealb_cat'  ,   alb_ice      * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice albedo for categories
+
+      IF( iom_use('icevol_cat'  ) )   CALL iom_put( "icevol_cat"  , v_i            * zmsk00l                                   ) ! volume for categories
+      IF( iom_use('iceaepnd_cat') )   CALL iom_put( 'iceaepnd_cat',   a_ip_eff     * zmsk00l                                   ) ! melt pond effective frac for categories
 
       !------------------
