Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icealb.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icealb.F90	(revision 12381)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icealb.F90	(revision 12382)
@@ -45,5 +45,5 @@
 CONTAINS
 
-   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, palb_cs, palb_os )
+   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, plh_pnd, palb_cs, palb_os )
       !!----------------------------------------------------------------------
       !!               ***  ROUTINE ice_alb  ***
@@ -97,6 +97,11 @@
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area)
       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   plh_pnd      !  melt pond lid thickness
       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_cs      !  albedo of ice under clear    sky
       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_os      !  albedo of ice under overcast sky
+      !
+      REAL(wp), PARAMETER :: pnd_lid_max = 0.015                  !  pond lid thickness above which the ponds disappear from the albedo calculation
+      REAL(wp), PARAMETER :: pnd_lid_min = 0.005                  !  pond lid thickness below which the full pond area is used in the albedo calculation
+                                                                  ! Note: these two variables are mirrored in sbccpl.F90 (maybe put them in one place...)
       !
       INTEGER  ::   ji, jj, jl                ! dummy loop indices
@@ -106,4 +111,5 @@
       REAL(wp) ::   zalb_ice, zafrac_ice      ! bare sea ice albedo & relative ice fraction
       REAL(wp) ::   zalb_snw, zafrac_snw      ! snow-covered sea ice albedo & relative snow fraction
+      REAL(wp) ::   lfrac_pnd                 ! The fraction of the meltpond exposed (not inder a frozen lid)
       !!---------------------------------------------------------------------
       !
@@ -123,5 +129,14 @@
                   zafrac_snw = 0._wp
                   IF( ld_pnd_alb ) THEN
-                     zafrac_pnd = pafrac_pnd(ji,jj,jl)
+                     IF ( plh_pnd(ji,jj,jl) > pnd_lid_max ) THEN
+                        lfrac_pnd = 0._wp
+                     ELSE
+                        IF ( plh_pnd(ji,jj,jl) < pnd_lid_min ) THEN
+                           lfrac_pnd = 1._wp
+                        ELSE
+                           lfrac_pnd = ( plh_pnd(ji,jj,jl) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min)
+                        END IF
+                     END IF
+                     zafrac_pnd = pafrac_pnd(ji,jj,jl) * lfrac_pnd
                   ELSE
                      zafrac_pnd = 0._wp
Index: /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icesbc.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icesbc.F90	(revision 12381)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icesbc.F90	(revision 12382)
@@ -128,5 +128,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_frac, h_ip, lh_ip, zalb_cs, zalb_os )
 
       ! albedo depends on cloud fraction because of non-linear spectral effects
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 12381)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/icethd_pnd.F90	(revision 12382)
@@ -129,11 +129,18 @@
       REAL(wp), PARAMETER ::   zpnd_aspect = 0.174_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) ::   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) ::   z1_Tp            ! inverse reference temperature
       REAL(wp) ::   z1_rhow          ! inverse freshwater density
       REAL(wp) ::   z1_zpnd_aspect   ! inverse pond aspect ratio
       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
@@ -142,4 +149,7 @@
       z1_zpnd_aspect = 1._wp / zpnd_aspect
       z1_Tp          = 1._wp / zTp 
+
+      ! Define time-independent field for use in refreezing
+      omega_dt = 2.0_wp * rcnd_i * rdtice / (rLfus * rhow)
 
       DO ji = 1, npti
@@ -163,25 +173,52 @@
             !
             ! 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)
+            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
             !
             !--- 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 )
-            !
-            ! melt pond mass flux (<0)
-            IF( zdv_mlt > 0._wp ) THEN
-               zfac = zfr_mlt * zdv_mlt * rhow * r1_rdtice
-               wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
-               !
-               ! adjust ice/snow melting flux to balance melt pond flux (>0)
-               zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )
-               wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
-               wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum)
-            ENDIF
+            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
+            !
             !
             !--- 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 ( t_su_1d(ji) < freezing_t .AND. v_ip_1d(ji) > 0._wp ) THEN
+               t_grad = freezing_t - t_su_1d(ji)
+               
+               ! The following equation is a rearranged form of:
+               ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdtice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow)
+               ! where: lid_thickness_start = lh_ip_1d(ji)
+               !        lid_thickness_end = lh_ip_end
+               !        omega_dt is a bunch of terms in the equation that do not change
+               ! note the use of rhow instead of rhoi as we are working with volumes and it is easier if the water and ice volumes (for the lid and the pond) are the same
+               ! (have the same density). The lid will eventually remelt anyway so it doesn't matter if we make this simplification.
+                             
+               lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2)
+               actual_frz = lh_ip_end - lh_ip_1d(ji)
+
+               ! Pond shrinking
+               v_ip_1d(ji) = v_ip_1d(ji) - actual_frz * a_ip_1d(ji)
+
+               ! Lid growing
+               lh_ip_1d(ji) = lh_ip_1d(ji) + actual_frz
+            ELSE
+               actual_frz = 0._wp
+            END IF
+
+            ! melt pond mass flux diagnostic (melt only)
+            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)
+            !
+            ! 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
@@ -190,4 +227,13 @@
             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 lid thickness is ten times greater than pond thickness then remove pond            
+            IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
+                a_ip_1d(ji)      = 0._wp
+                a_ip_frac_1d(ji) = 0._wp
+                h_ip_1d(ji)      = 0._wp
+                lh_ip_1d(ji)     = 0._wp
+                v_ip_1d(ji)      = 0._wp
+            END IF
             !
          ENDIF
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 12381)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/ICE/iceupdate.F90	(revision 12382)
@@ -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_frac, h_ip, lh_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_add_pond_lids_prints/src/OCE/SBC/sbccpl.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/SBC/sbccpl.F90	(revision 12381)
+++ /NEMO/branches/UKMO/NEMO_4.0_add_pond_lids_prints/src/OCE/SBC/sbccpl.F90	(revision 12382)
@@ -2122,4 +2122,8 @@
       INTEGER, INTENT(in) ::   kt
       !
+      REAL(wp), PARAMETER :: pnd_lid_max = 0.015                  !  pond lid thickness above which the ponds disappear from the albedo calculation
+      REAL(wp), PARAMETER :: pnd_lid_min = 0.005                  !  pond lid thickness below which the full pond area is used in the albedo calculation
+                                                                  ! Note: these two variables are mirrored in icealb.F90 (maybe put them in one place...)
+      !
       INTEGER ::   ji, jj, jl   ! dummy loop indices
       INTEGER ::   isec, info   ! local integer
@@ -2127,4 +2131,6 @@
       REAL(wp), DIMENSION(jpi,jpj)     ::   zfr_l, ztmp1, ztmp2, zotx1, zoty1, zotz1, zitx1, zity1, zitz1
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   ztmp3, ztmp4   
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   lfrac_pnd                 ! The fraction of the meltpond exposed (not inder a frozen lid)
+
       !!----------------------------------------------------------------------
       !
@@ -2341,5 +2347,15 @@
             SELECT CASE( sn_snd_mpnd%clcat )  
             CASE( 'yes' )  
-               ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl)
+
+               ! Calculate how much meltpond is exposed (not under a frozen lid)
+               lfrac_pnd(:,:,1:jpl) = 1.0
+               WHERE( lh_ip(:,:,1:jpl) > pnd_lid_max )
+                    lfrac_pnd(:,:,1:jpl) = 0.0
+               END WHERE
+               WHERE( lh_ip(:,:,1:jpl) > pnd_lid_min .AND. lh_ip(:,:,1:jpl) <= pnd_lid_max )
+                    lfrac_pnd(:,:,1:jpl) = ( lh_ip(:,:,1:jpl) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min)
+               END WHERE
+
+               ztmp3(:,:,1:jpl) =  a_ip_frac(:,:,1:jpl) * lfrac_pnd(:,:,1:jpl)
                ztmp4(:,:,1:jpl) =  h_ip(:,:,1:jpl)
             CASE( 'no' )  
