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 12483)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/ice1d.F90	(revision 12484)
@@ -111,4 +111,6 @@
    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]
@@ -129,4 +131,5 @@
    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
@@ -211,5 +214,6 @@
          &      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) ,                                                   &
+         &      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.1_fix_cpl/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd.F90	(revision 12483)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd.F90	(revision 12484)
@@ -220,4 +220,6 @@
             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 12483)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_dh.F90	(revision 12484)
@@ -24,4 +24,5 @@
    USE lib_mpp        ! MPP library
    USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
+   USE icethd_pnd, ONLY : a_pnd_avail
    
    IMPLICIT NONE
@@ -98,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 
 
@@ -112,4 +114,9 @@
          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
@@ -242,4 +249,5 @@
                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)
@@ -248,4 +256,5 @@
                ! 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) )
@@ -343,4 +352,5 @@
                
                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]
@@ -349,4 +359,5 @@
                
                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 12483)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_pnd.F90	(revision 12484)
@@ -37,4 +37,10 @@
    INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant pond scheme
    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
+
+   REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp
 
    !! * Substitutions
@@ -129,24 +135,37 @@
       REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp   ! pond aspect ratio
       REAL(wp), PARAMETER ::   zTp         = -2._wp   ! reference temperature
-      REAL(wp), PARAMETER ::   freezing_t  = 273.0_wp ! temperature below which refreezing occurs
+      REAL(wp), PARAMETER ::   max_h_diff_s = -1.0E-6  ! Maximum meltpond depth change due to leaking or overflow (m s-1)
       !
       REAL(wp) ::   zfr_mlt          ! fraction of available meltwater retained for melt ponding
-      REAL(wp) ::   zdv_mlt          ! available meltwater for melt ponding
-      REAL(wp) ::   actual_mlt       ! actual melt water used to make melt ponds (per m2).
-                                     ! Need to multiply this by ice area to work on volumes.
+      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
       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) ::   actual_frz       ! Amount of melt pond that freezes
-      !
-      INTEGER  ::   ji   ! loop indices
+      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) ::   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
+      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)
+      
+      !
+      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
@@ -154,8 +173,9 @@
 
       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
@@ -163,27 +183,29 @@
             lh_ip_1d(ji)     = 0._wp
 
-            actual_mlt = 0._wp
-            actual_frz = 0._wp
+            zdv_mlt = 0._wp
+            zdh_frz = 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
-            zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji)  ! from CICE doc
-            !zfr_mlt = zrmin + zrmax * a_i_1d(ji)             ! from Holland paper 
-            actual_mlt = zfr_mlt * zdv_mlt
+
+            ! 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 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) + actual_mlt * a_i_1d(ji)
-            !
-            !--- Lid shrinking ---!
-            lh_ip_1d(ji) = lh_ip_1d(ji) - actual_mlt
+            v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
+            !
+            !--- Lid shrinking. ---!
+            lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip
             !
             !
             !--- Pond contraction (due to refreezing) ---!
-            IF ( t_su_1d(ji) < freezing_t .AND. v_ip_1d(ji) > 0._wp ) THEN
-               t_grad = freezing_t - t_su_1d(ji)
+            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:
@@ -196,25 +218,15 @@
                              
                lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2)
-               actual_frz = lh_ip_end - lh_ip_1d(ji)
+               zdh_frz = lh_ip_end - lh_ip_1d(ji)
 
                ! Pond shrinking
-               v_ip_1d(ji) = v_ip_1d(ji) - actual_frz * a_ip_1d(ji)
+               v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji)
 
                ! Lid growing
-               lh_ip_1d(ji) = lh_ip_1d(ji) + actual_frz
+               lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz
             ELSE
-               actual_frz = 0._wp
+               zdh_frz = 0._wp
             END IF
 
-            ! melt pond mass flux (<0)
-            IF( zdv_mlt > 0._wp ) THEN
-               zfac = actual_mlt * a_i_1d(ji) * 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
             !
             ! Make sure pond volume or lid thickness has not gone negative
@@ -227,4 +239,60 @@
             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) = h_ip_1d(ji) + dh_ip_over
+                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)
+
+            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) = 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) = 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            
@@ -236,6 +304,15 @@
                 v_ip_1d(ji)      = 0._wp
             END IF
+
+            ! 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
+
+      END DO
       END DO
       !
