Index: /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/ice.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/ice.F90	(revision 12401)
@@ -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     !
    !!                                                                     |
    !!-------------|-------------|---------------------------------|-------|
@@ -332,4 +333,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_frac  !: melt pond fraction (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 +450,5 @@
 
       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) , STAT = ierr(ii) )
 
       ii = ii + 1
Index: /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/ice1d.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/ice1d.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/ice1d.F90	(revision 12401)
@@ -128,4 +128,5 @@
    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(:) ::   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 +158,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,5 +210,5 @@
          &      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) ,                   &
+         &      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) ,                                                   &
          &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) )
@@ -226,5 +228,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_add_pond_lids/src/ICE/icealb.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icealb.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icealb.F90	(revision 12401)
@@ -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.1_add_pond_lids/src/ICE/icedyn.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn.F90	(revision 12401)
@@ -101,4 +101,5 @@
       ELSEWHERE
          h_ip(:,:,:) = 0._wp
+         lh_ip(:,:,:) = 0._wp
       END WHERE
       !
Index: /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_adv.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_adv.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_adv.F90	(revision 12401)
@@ -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_add_pond_lids/src/ICE/icedyn_adv_pra.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_adv_pra.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_adv_pra.F90	(revision 12401)
@@ -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_add_pond_lids/src/ICE/icedyn_adv_umx.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_adv_umx.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_adv_umx.F90	(revision 12401)
@@ -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_add_pond_lids/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_rdgrft.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icedyn_rdgrft.F90	(revision 12401)
@@ -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_add_pond_lids/src/ICE/iceistate.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/iceistate.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/iceistate.F90	(revision 12401)
@@ -155,4 +155,5 @@
       a_ip     (:,:,:) = 0._wp
       v_ip     (:,:,:) = 0._wp
+      lh_ip    (:,:,:) = 0._wp
       a_ip_frac(:,:,:) = 0._wp
       h_ip     (:,:,:) = 0._wp
Index: /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/iceitd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/iceitd.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/iceitd.F90	(revision 12401)
@@ -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_add_pond_lids/src/ICE/icerst.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icerst.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icerst.F90	(revision 12401)
@@ -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 
@@ -245,8 +246,10 @@
             CALL iom_get( numrir, jpdom_autoglo, 'a_ip' , a_ip )
             CALL iom_get( numrir, jpdom_autoglo, 'v_ip' , v_ip )
+            CALL iom_get( numrir, jpdom_autoglo, 'lh_ip', lh_ip)
          ELSE                                     ! start from rest
             IF(lwp) WRITE(numout,*) '   ==>>   previous run without melt ponds output then set it to zero'
             a_ip(:,:,:) = 0._wp
             v_ip(:,:,:) = 0._wp
+            lh_ip(:,:,:) = 0._wp
          ENDIF
          ! fields needed for Met Office (Jules) coupling
Index: /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icesbc.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icesbc.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icesbc.F90	(revision 12401)
@@ -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_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.1_add_pond_lids/src/ICE/icethd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icethd.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icethd.F90	(revision 12401)
@@ -355,4 +355,5 @@
          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), 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 +462,5 @@
          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), 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_add_pond_lids/src/ICE/icethd_pnd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icethd_pnd.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icethd_pnd.F90	(revision 12401)
@@ -129,11 +129,18 @@
       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) ::   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 * rdt_ice / (rLfus * rhow)
 
       DO ji = 1, npti
@@ -151,4 +161,8 @@
             a_ip_frac_1d(ji) = 0._wp
             h_ip_1d(ji)      = 0._wp
+            lh_ip_1d(ji)     = 0._wp
+
+            actual_mlt = 0._wp
+            actual_frz = 0._wp
             !                                                     !--------------------------------!
          ELSE                                                     ! Case ice thickness >= rn_himin !
@@ -157,15 +171,43 @@
             !
             ! 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 )
-            !
+            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) ---!
+            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 * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow)
+               ! where: lid_thickness_start = lh_ip_1d(ji)
+               !        lid_thickness_end = lh_ip_end
+               !        omega_dt is a bunch of terms in the equation that do not change
+               ! note the use of rhow instead of rhoi as we are working with volumes and it is easier if the water and ice volumes (for the lid and the pond) are the same
+               ! (have the same density). The lid will eventually remelt anyway so it doesn't matter if we make this simplification.
+                             
+               lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2)
+               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 (<0)
             IF( zdv_mlt > 0._wp ) THEN
-               zfac = zfr_mlt * zdv_mlt * rhow * r1_rdtice
+               zfac = actual_mlt * a_i_1d(ji) * rhow * r1_rdtice
                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
                !
@@ -176,6 +218,7 @@
             ENDIF
             !
-            !--- 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 )
+            ! 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,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.1_add_pond_lids/src/ICE/iceupdate.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/iceupdate.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/iceupdate.F90	(revision 12401)
@@ -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.1_add_pond_lids/src/ICE/icevar.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icevar.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/ICE/icevar.F90	(revision 12401)
@@ -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_add_pond_lids/src/OCE/BDY/bdyice.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/OCE/BDY/bdyice.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/OCE/BDY/bdyice.F90	(revision 12401)
@@ -95,5 +95,5 @@
             ! exchange 3d arrays
             CALL lbc_lnk_multi( 'bdyice', a_i , 'T', 1., h_i , 'T', 1., h_s , 'T', 1., oa_i, 'T', 1. &
-                 &                      , a_ip, 'T', 1., v_ip, 'T', 1., s_i , 'T', 1., t_su, 'T', 1. &
+                 &                      , a_ip, 'T', 1., v_ip, 'T', 1., lh_ip,'T', 1., s_i , 'T', 1., t_su, 'T', 1. &
                  &                      , v_i , 'T', 1., v_s , 'T', 1., sv_i, 'T', 1.                &
                  &                      , kfillmode=jpfillnothing ,lsend=llsend1, lrecv=llrecv1      )
@@ -163,4 +163,5 @@
             a_ip(ji,jj,  jl) = ( a_ip(ji,jj,  jl) * zwgt1 + dta%aip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond concentration
             h_ip(ji,jj,  jl) = ( h_ip(ji,jj,  jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond depth
+            lh_ip(ji,jj, jl) = ( lh_ip(ji,jj, jl) * zwgt1 + dta%hip(i_bdy,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice  pond lid thickness
             !
             sz_i(ji,jj,:,jl) = s_i(ji,jj,jl)
@@ -170,4 +171,5 @@
                a_ip(ji,jj,jl) = 0._wp
                h_ip(ji,jj,jl) = 0._wp
+               lh_ip(ji,jj,jl) = 0._wp
             ENDIF
             !
@@ -231,4 +233,5 @@
                a_ip(ji,jj,  jl) = a_ip(ib,jb,  jl)
                h_ip(ji,jj,  jl) = h_ip(ib,jb,  jl)
+               lh_ip(ji,jj, jl) = lh_ip(ib,jb, jl)
                !
                sz_i(ji,jj,:,jl) = sz_i(ib,jb,:,jl)
@@ -280,4 +283,5 @@
                a_ip(ji,jj,  jl) = 0._wp
                v_ip(ji,jj,  jl) = 0._wp
+               lh_ip(ji,jj, jl) = 0._wp
                t_su(ji,jj,  jl) = rt0
                t_s (ji,jj,:,jl) = rt0
Index: /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/OCE/SBC/sbccpl.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/OCE/SBC/sbccpl.F90	(revision 12400)
+++ /NEMO/branches/UKMO/NEMO_4.0.1_add_pond_lids/src/OCE/SBC/sbccpl.F90	(revision 12401)
@@ -2047,4 +2047,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
@@ -2052,4 +2056,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)
+
       !!----------------------------------------------------------------------
       !
@@ -2256,6 +2262,17 @@
             SELECT CASE( sn_snd_mpnd%clcat )  
             CASE( 'yes' )  
-               ztmp3(:,:,1:jpl) =  a_ip(:,:,1:jpl)
-               ztmp4(:,:,1:jpl) =  v_ip(:,:,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' )  
                ztmp3(:,:,:) = 0.0  
