Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/cfgs/SHARED/field_def_nemo-ice.xml	(revision 12448)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/cfgs/SHARED/field_def_nemo-ice.xml	(revision 12449)
@@ -115,4 +115,5 @@
           <field id="sfxres"       long_name="ice-ocean salt flux from undiagnosed processes"               unit="kg/m2/s" />
           <field id="sfxsub"       long_name="ice-ocean salt flux from ice sublimation"                     unit="kg/m2/s" />
+          <field id="sfxpnd"       long_name="ice-ocean salt flux from meltponds leaking into the ocean"    unit="kg/m2/s" />
 
 	  <!-- mass fluxes -->
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice.F90	(revision 12448)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice.F90	(revision 12449)
@@ -266,4 +266,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: salt flux due to correction on ice thick. (residual) [pss.kg.m-2.s-1 => g.m-2.s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sub     !: salt flux due to ice sublimation                     [pss.kg.m-2.s-1 => g.m-2.s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_pnd     !: salt flux due to melt ponds leaking into the ocean   [pss.kg.m-2.s-1 => g.m-2.s-1]
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth           [W.m-2]
@@ -418,4 +419,5 @@
          &      sfx_res    (jpi,jpj) , sfx_bri   (jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) ,  &
          &      sfx_bog    (jpi,jpj) , sfx_bom   (jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,  &
+         &      sfx_pnd    (jpi,jpj) ,                                                                 &
          &      hfx_res    (jpi,jpj) , hfx_snw   (jpi,jpj) , hfx_sub(jpi,jpj) ,                        & 
          &      qt_atm_oi  (jpi,jpj) , qt_oce_ai (jpi,jpj) , fhld   (jpi,jpj) ,                        &
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice1d.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice1d.F90	(revision 12448)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/ice1d.F90	(revision 12449)
@@ -93,4 +93,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_lam_1d
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_dyn_1d
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_pnd_1d 
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d
@@ -207,5 +208,5 @@
          &      sfx_bri_1d    (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) , sfx_sum_1d (jpij),  &
          &      sfx_sni_1d    (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , sfx_sub_1d (jpij),  &
-         &      sfx_lam_1d    (jpij) , sfx_dyn_1d(jpij)  , STAT=ierr(ii) )
+         &      sfx_lam_1d    (jpij) , sfx_dyn_1d(jpij)  , sfx_pnd_1d (jpij) , STAT=ierr(ii) )
       !
       ii = ii + 1
@@ -216,4 +217,5 @@
          &      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) ,                                                                        &
          &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) )
       !
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icectl.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icectl.F90	(revision 12448)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icectl.F90	(revision 12449)
@@ -87,6 +87,6 @@
          pdiag_fs = glob_sum( 'icectl',                                                                     &
             &                  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                    sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    &
-            &                  ) *  e1e2t(:,:) ) * zconv 
+            &                    sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) +  &
+            &                    sfx_pnd(:,:) ) *  e1e2t(:,:) ) * zconv 
          !
          !                          ! heat flux
@@ -116,6 +116,6 @@
          zfs = glob_sum( 'icectl',                                                                       &
             &              ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:)    & 
-            &              ) * e1e2t(:,:) ) * zconv - pdiag_fs
+            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) +  & 
+            &                sfx_pnd(:,:) ) * e1e2t(:,:) ) * zconv - pdiag_fs
 
          ! heat flux
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icestp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icestp.F90	(revision 12448)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icestp.F90	(revision 12449)
@@ -486,4 +486,5 @@
       sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp
       sfx_res(:,:) = 0._wp   ;   sfx_sub(:,:) = 0._wp
+      sfx_pnd(:,:) = 0._wp
       !
       wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd.F90	(revision 12448)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd.F90	(revision 12449)
@@ -455,4 +455,5 @@
          CALL tab_2d_1d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub          )
          CALL tab_2d_1d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam          )
+         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_pnd_1d (1:npti), sfx_pnd          )
          !
          CALL tab_2d_1d( npti, nptidx(1:npti), hfx_thd_1d    (1:npti), hfx_thd       )
@@ -548,4 +549,5 @@
          CALL tab_1d_2d( npti, nptidx(1:npti), sfx_sub_1d (1:npti), sfx_sub        )
          CALL tab_1d_2d( npti, nptidx(1:npti), sfx_lam_1d (1:npti), sfx_lam        )
+         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_pnd_1d (1:npti), sfx_pnd        )
          !
          CALL tab_1d_2d( npti, nptidx(1:npti), hfx_thd_1d    (1:npti), hfx_thd     )
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_dh.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_dh.F90	(revision 12448)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_dh.F90	(revision 12449)
@@ -99,4 +99,5 @@
       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 
 
@@ -281,15 +282,20 @@
                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)
-               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip)
+               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos           * a_i_1d(ji) * (zdeltah(ji,jk) - zdeltah_pnd(ji,jk)) * r1_rdtice   ! snow melting only = water into the ocean (then without snow precip)
                
                ! updates available heat + thickness
                dh_s_mlt(ji)    = dh_s_mlt(ji) + zdeltah(ji,jk)
-               dh_s_pnd(ji)    = dh_s_pnd(ji) + zdeltah(ji,jk) * a_pnd_avail_1d(ji)                   ! Cumulate surface melt on areas contributing to melt ponds
+               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) )
-               h_i_1d  (ji)    = MAX( 0._wp , h_i_1d(ji) - zdeltah(ji,jk) * a_pnd_avail_1d(ji) * rhos * r1_rhoi )   ! Snow melt trapped in ponds contrbutes to ice thickness
                zh_s    (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) )
+
+               ! Snow melt trapped in ponds contrbutes to ice thickness and snow to ice conversion
+               h_i_1d  (ji)    = MAX( 0._wp , h_i_1d(ji) - zdeltah_pnd(ji,jk) * rhos * r1_rhoi ) 
+               wfx_sni_1d(ji)     = wfx_sni_1d(ji)     - a_i_1d(ji) * zdeltah_pnd(ji,jk) * rhos * r1_rdtice
+               wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * zdeltah_pnd(ji,jk) * rhos * r1_rdtice
                !
             ENDIF
@@ -400,7 +406,8 @@
                zdE            =    zEi - zEw                          ! Specific enthalpy difference < 0
                
-               zfmdt          = - zq_top(ji) / zdE                    ! Mass flux to the ocean [kg/m2, >0]
-               
-               zdeltah(ji,jk) = - zfmdt * r1_rhoi                     ! Melt of layer jk [m, <0]
+               zfmdt          = - zq_top(ji) / zdE                    ! Mass flux to the ocean or meltponds [kg/m2, >0]
+               
+               zdeltah(ji,jk) = - zfmdt * r1_rhoi * (1._wp - a_pnd_avail_1d(ji))                    ! Melt of layer jk [m, <0] into the ocean
+               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]
@@ -408,6 +415,6 @@
                zq_top(ji)      = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat
                
-               dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk) * (1._wp - a_pnd_avail_1d(ji))         ! Cumulate surface melt on areas not contributing to melt ponds
-               dh_i_pnd(ji)   = dh_i_pnd(ji) + zdeltah(ji,jk) * a_pnd_avail_1d(ji)                   ! Cumulate surface melt on areas contributing to melt ponds
+               dh_i_sum(ji)   = dh_i_sum(ji) + zdeltah(ji,jk)          ! Cumulate surface melt on areas not contributing to melt ponds
+               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_add_pond_lids_prints/src/ICE/icethd_pnd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90	(revision 12448)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90	(revision 12449)
@@ -39,9 +39,8 @@
 
    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    !  pond lid thickness above which the ponds disappear from the albedo calculation
-   REAL(wp), PUBLIC, PARAMETER :: pnd_lid_min = 0.005    !  pond lid thickness below which the full pond area is used in the albedo calculation
-
-   REAL(wp), PARAMETER :: snow_lid_min = 0.15  ! Snow thickness above which form a lid of size pnd_lid_min on the melt ponds
-   REAL(wp), PARAMETER :: snow_lid_max = 0.25  ! Snow thickness above which form a lid of size pnd_lid_max on the melt ponds
+   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
+
+   REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp
 
    !! * Substitutions
@@ -138,5 +137,5 @@
       !
       REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding
-      REAL(wp) ::   zdh_mlt          ! available meltwater for melt ponding (equivalent thickness change)
+      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
@@ -151,6 +150,14 @@
       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
-      !
-      INTEGER  ::   ji   ! loop indices
+      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) ::   drain            ! Amount of melt pond draining the sea ice per m2
+      REAL(wp) ::   za_ip            ! Temporary pond fraction
+ 
+      !
+      INTEGER  ::   ji, jk   ! loop indices
       !!-------------------------------------------------------------------
       z1_rhow        = 1._wp / rhow 
@@ -168,8 +175,8 @@
          END IF
 
-         !                                                        !--------------------------------!
-         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
@@ -177,5 +184,5 @@
             lh_ip_1d(ji)     = 0._wp
 
-            zdh_mlt = 0._wp
+            zdv_mlt = 0._wp
             zdh_frz = 0._wp
             !                                                     !--------------------------------!
@@ -183,15 +190,18 @@
             !                                                     !--------------------------------!
             v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! record pond volume at previous time step
+
+            ! 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
-            ! This is the change in ice thickness due to melt scaled up by the realive areas of the meltpond
-            ! and the area of sea ice feeding the melt ponds.
-            zdh_mlt = -(dh_i_pnd(ji)*rhoi + dh_s_pnd(ji)*rhos) * z1_rhow * a_i_1d(ji) / a_ip_1d(ji)
+            ! 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)
             !
             !--- Pond gowth ---!
-            v_ip_1d(ji) = v_ip_1d(ji) + zdh_mlt * a_ip_1d(ji)
+            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
             !
             !--- Lid shrinking. ---!
-            lh_ip_1d(ji) = lh_ip_1d(ji) - zdh_mlt
+            lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip
             !
             !
@@ -231,4 +241,90 @@
             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
+                h_ip_1d(ji) = zpnd_aspect * a_pnd_avail_1d(ji)
+                a_ip_1d(ji) = a_pnd_avail_1d(ji) * a_i_1d(ji)
+                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
+                a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
+
+                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean.
+                ! Scale this up to make water depth meaned over sea ice.
+                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number
+
+                ! Output this water loss as a mass flux diagnostic
+                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice
+
+                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference
+                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi )
+
+                ! Pass this as a salinity flux to the ocean
+                sfx_pnd_1d(ji) = sfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * s_i_1d(ji) * r1_rdtice
+            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
+                h_ip_1d(ji) = 0.5_wp * h_i_1d(ji)
+                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)
+
+                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean.
+                ! Scale this up to make water depth meaned over sea ice.
+                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number
+
+                ! Output this water loss as a mass flux diagnostic
+                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice
+
+                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference
+                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi )
+
+                ! Pass this as a salinity flux to the ocean
+                sfx_pnd_1d(ji) = sfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * s_i_1d(ji) * r1_rdtice
+            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
+                drain = perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation)
+
+                v_ip_old = v_ip_1d(ji)   ! Save original volume before leak for future use
+                h_ip_1d(ji) = h_ip_1d(ji) - drain
+                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)
+
+                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean.
+                ! Scale this up to make water depth meaned over sea ice.
+                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number
+
+                ! Output this water loss as a mass flux diagnostic
+                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice
+
+                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference
+                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi )
+
+                ! Pass this as a salinity flux to the ocean
+                sfx_pnd_1d(ji) = sfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * s_i_1d(ji) * r1_rdtice
+            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
@@ -240,43 +336,5 @@
             END IF
 
-            ! 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
-                h_ip_1d(ji) = zpnd_aspect * a_pnd_avail_1d(ji)
-                a_ip_1d(ji) = a_pnd_avail_1d(ji) * a_i_1d(ji)
-                v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)
-                a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji)
-
-                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean.
-                ! Scale this up to make water depth meaned over sea ice.
-                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number
-
-                ! Output this water loss as a mass flux diagnostic
-                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice
-
-                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference
-                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi )
-            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
-                h_ip_1d(ji) = 0.5_wp * h_i_1d(ji)
-                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)
-
-                ! A change in volume of a melt pond is equivalent to a change in water depth in the grid box mean.
-                ! Scale this up to make water depth meaned over sea ice.
-                dh_i_pndleak = (v_ip_1d(ji) - v_ip_old) / a_i_1d(ji)   ! This will be a negative number
-
-                ! Output this water loss as a mass flux diagnostic
-                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - rhow * a_i_1d(ji) * dh_i_pndleak * r1_rdtice
-
-                ! Reduce the ice thickness using densities to convert from a water depth difference to a sea ice thickness difference
-                h_i_1d(ji) =  MAX( 0._wp , h_i_1d(ji) + dh_i_pndleak * rhow * z1_rhoi )
-            ENDIF
-
-            ! If any of the previous two IF tests has removed all the ice thickness then remove ice area.
+            ! 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
@@ -284,14 +342,4 @@
             ENDIF
 
-            ! If snow thickness exceeds snow_lid_min then form a very thin lid (so snow can go over the top)
-            IF ( h_s_1d(ji) > snow_lid_min .AND. h_s_1d(ji) < snow_lid_max ) THEN
-               weighted_lid_snow = (snow_lid_max - h_s_1d(ji))/(snow_lid_max-snow_lid_min) * pnd_lid_min +             &
-                                   (h_s_1d(ji) - snow_lid_min)/(snow_lid_max-snow_lid_min) * pnd_lid_max
-               lh_ip_1d(ji) = MAX(lh_ip_1d(ji), weighted_lid_snow)
-            ENDIF
-            IF ( h_s_1d(ji) >= snow_lid_max ) THEN
-               lh_ip_1d(ji) = MAX(lh_ip_1d(ji), pnd_lid_max)
-            ENDIF
-
             !
          ENDIF
@@ -299,8 +347,10 @@
          IF (to_print(ji) == 10) THEN
             write(numout,*)'icethd_pnd: h_ip_1d(ji), zpnd_aspect, a_ip_frac_1d(ji), a_ip_1d(ji) = ',h_ip_1d(ji), zpnd_aspect, a_ip_frac_1d(ji), a_ip_1d(ji)
-            write(numout,*)'icethd_pnd: a_i_1d(ji), v_ip_1d(ji), t_su_1d(ji), zfr_mlt, zdh_mlt = ',a_i_1d(ji), ' ', v_ip_1d(ji), ' ', t_su_1d(ji), ' ', zfr_mlt, ' ', zdh_mlt
+            write(numout,*)'icethd_pnd: a_i_1d(ji), v_ip_1d(ji), t_su_1d(ji), zfr_mlt, zdv_mlt = ',a_i_1d(ji), ' ', v_ip_1d(ji), ' ', t_su_1d(ji), ' ', zfr_mlt, ' ', zdv_mlt
             write(numout,*)'icethd_pnd: meltt = ', -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) / rhoi
-            write(numout,*)'icethd_pnd: lh_ip_1d(ji), zdh_mlt, zdh_frz, t_su_1d(ji) = ',lh_ip_1d(ji), ' ', zdh_mlt, ' ', zdh_frz, ' ', t_su_1d(ji)
+            write(numout,*)'icethd_pnd: lh_ip_1d(ji), sfx_pnd_1d(ji), zdh_frz, t_su_1d(ji) = ',lh_ip_1d(ji), ' ', sfx_pnd_1d(ji), ' ', zdh_frz, ' ', t_su_1d(ji)
             write(numout,*)'icethd_pnd: a_pnd_avail_1d(ji), at_i_1d(ji), wfx_pnd_1d(ji), h_i_1d(ji) = ', a_pnd_avail_1d(ji), ' ', at_i_1d(ji), ' ', wfx_pnd_1d(ji), ' ', h_i_1d(ji)
+            write(numout,*)'icethd_pnd: drain, perm, h_percolation, h_gravity_head = ',drain, ' ', perm, ' ', h_percolation, ' ', h_gravity_head
+            write(numout,*)'icethd_pnd: t_i_1d(ji,1), sz_i_1d(ji,1) = ',t_i_1d(ji,1), ' ', sz_i_1d(ji,1)
          END IF
 
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_sal.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_sal.F90	(revision 12448)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_sal.F90	(revision 12449)
@@ -59,5 +59,5 @@
       INTEGER  ::   ji, jk                       ! dummy loop indices 
       REAL(wp) ::   iflush, igravdr              ! local scalars
-      REAL(wp) ::   zs_sni, zs_i_gd, zs_i_fl, zs_i_si, zs_i_bg   ! local scalars
+      REAL(wp) ::   zs_sni, zs_i_gd, zs_i_fl, zs_i_si, zs_i_bg, zs_i_sp   ! local scalars
       REAL(wp) ::   z1_time_gd, z1_time_fl
       !!---------------------------------------------------------------------
@@ -79,7 +79,8 @@
                zs_sni  = sss_1d(ji) * ( rhoi - rhos ) * r1_rhoi                     ! Salinity of snow ice
                zs_i_si = ( zs_sni      - s_i_1d(ji) ) * dh_snowice(ji) / h_i_1d(ji) ! snow-ice    
+               zs_i_sp = ( 0._wp       - s_i_1d(ji) ) * dh_s_pnd(ji)   / h_i_1d(ji) ! snow-pond
                zs_i_bg = ( s_i_new(ji) - s_i_1d(ji) ) * dh_i_bog  (ji) / h_i_1d(ji) ! bottom growth
                ! Update salinity (nb: salt flux already included in icethd_dh)
-               s_i_1d(ji) = s_i_1d(ji) + zs_i_bg + zs_i_si
+               s_i_1d(ji) = s_i_1d(ji) + zs_i_bg + zs_i_si + zs_i_sp
             ENDIF
             !
@@ -100,4 +101,5 @@
                sfx_bri_1d(ji) = sfx_bri_1d(ji) - rhoi * a_i_1d(ji) * h_i_1d(ji) * ( zs_i_fl + zs_i_gd ) * r1_rdtice
             ENDIF
+            !
          END DO
          !
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/iceupdate.F90	(revision 12448)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/iceupdate.F90	(revision 12449)
@@ -165,5 +165,5 @@
             !------------------------------------------
             sfx(ji,jj) = sfx_bog(ji,jj) + sfx_bom(ji,jj) + sfx_sum(ji,jj) + sfx_sni(ji,jj) + sfx_opw(ji,jj)   &
-               &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj)
+               &       + sfx_res(ji,jj) + sfx_dyn(ji,jj) + sfx_bri(ji,jj) + sfx_sub(ji,jj) + sfx_lam(ji,jj) + sfx_pnd(ji,jj)
             
             ! Mass of snow and ice per unit area   
@@ -197,5 +197,5 @@
       !
       ! --- salt fluxes [kg/m2/s] --- !
-      !                           ! sfxice =  sfxbog + sfxbom + sfxsum + sfxsni + sfxopw + sfxres + sfxdyn + sfxbri + sfxsub + sfxlam
+      !                           ! sfxice =  sfxbog + sfxbom + sfxsum + sfxsni + sfxopw + sfxres + sfxdyn + sfxbri + sfxsub + sfxlam + sfxpnd
       IF( iom_use('sfxice'  ) )   CALL iom_put( "sfxice", sfx     * 1.e-03 )   ! salt flux from total ice growth/melt
       IF( iom_use('sfxbog'  ) )   CALL iom_put( "sfxbog", sfx_bog * 1.e-03 )   ! salt flux from bottom growth
@@ -209,4 +209,5 @@
       IF( iom_use('sfxres'  ) )   CALL iom_put( "sfxres", sfx_res * 1.e-03 )   ! salt flux from undiagnosed processes
       IF( iom_use('sfxsub'  ) )   CALL iom_put( "sfxsub", sfx_sub * 1.e-03 )   ! salt flux from sublimation
+      IF( iom_use('sfxpnd'  ) )   CALL iom_put( "sfxpnd", sfx_pnd * 1.e-03 )   ! salt flux from melt ponds leaking into the ocean
 
       ! --- mass fluxes [kg/m2/s] --- !
