Index: /NEMO/branches/UKMO/NEMO_4.0.4_fix_topo_ponds/src/ICE/icethd_pnd.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_fix_topo_ponds/src/ICE/icethd_pnd.F90	(revision 14668)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_fix_topo_ponds/src/ICE/icethd_pnd.F90	(revision 14669)
@@ -520,8 +520,8 @@
          zfsurf,  &      ! net heat flux, excluding conduction and transmitted radiation (W/m2)
          zTp,     &      ! pond freezing temperature (C)
+         lid_Ttop,&      ! lid top temperature (K)
          zrhoi_L, &      ! volumetric latent heat of sea ice (J/m^3)
          zfr_mlt, &      ! fraction and volume of available meltwater retained for melt ponding
          z1_rhow, &      ! inverse water density
-         zv_pnd  , &     ! volume of meltwater contributing to ponds
          zv_mlt          ! total amount of meltwater produced
 
@@ -530,5 +530,6 @@
                                         zvolp_res    !! remaining melt pond water available after drainage
                                         
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, &
+                                            zv_pnd       !! volume of meltwater contributing to ponds
 
       INTEGER  ::   ji, jj, jk, jl                    ! loop indices
@@ -594,5 +595,5 @@
                       ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area
                   zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj)                  ! fraction of surface meltwater going to ponds
-                  zv_pnd  = zfr_mlt * zv_mlt                                                           ! contributing meltwater volume for category jl
+                  zv_pnd(ji,jj,jl)  = zfr_mlt * zv_mlt                                                           ! contributing meltwater volume for category jl
 
                   diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_rdtice                     ! diags
@@ -608,5 +609,5 @@
                   
                   !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage
-                  v_ip(ji,jj,jl) = v_ip(ji,jj,jl) +  zv_pnd                                            ! use pond water to increase thickness
+                  v_ip(ji,jj,jl) = v_ip(ji,jj,jl) +  zv_pnd(ji,jj,jl)                                            ! use pond water to increase thickness
                   h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
                   
@@ -652,94 +653,63 @@
                DO jl = 1, jpl-1
             
-                  IF ( v_il(ji,jj,jl) > epsi10 ) THEN
                
-                     !----------------------------------------------------------------
-                     ! Lid melting: floating upper ice layer melts in whole or part
-                     !----------------------------------------------------------------
-                     ! Use Tsfc for each category
-                     IF ( t_su(ji,jj,jl) > zTp ) THEN
-
-                        zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) )
-                        
-                        IF ( zdvice > epsi10 ) THEN
-                        
-                           v_il (ji,jj,jl) = v_il (ji,jj,jl)   - zdvice
-                           v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)    + zdvice ! MV: not sure i understand dh_i_sum seems counted twice - 
-                                                                        ! as it is already counted in surface melt
-!                          zvolp(ji,jj)     = zvolp(ji,jj)     + zdvice ! pointless to calculate total volume (done in icevar)
-!                          zfpond(ji,jj)    = fpond(ji,jj)     + zdvice ! pointless to follow fw budget (ponds have no fw)
-                     
-                           IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN
-                           ! ice lid melted and category is pond covered
-                              v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)  + v_il(ji,jj,jl) 
-!                             zfpond(ji,jj)    = zfpond (ji,jj)    + v_il(ji,jj,jl) 
-                              v_il(ji,jj,jl)   = 0._wp
-                           ENDIF
-                           h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here
-                           
-                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice   ! diag
-                           
-                        ENDIF
-                        
-                     !----------------------------------------------------------------
-                     ! Freeze pre-existing lid 
-                     !----------------------------------------------------------------
-
-                     ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp
-
-                        ! differential growth of base of surface floating ice layer
-                        ! zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0  ! line valid for temp in C 
-                        zdTice = MAX( - ( t_su(ji,jj,jl) - zTd ) , 0._wp ) ! > 0   
-                        zomega = rcnd_i * zdTice / zrhoi_L
-                        zdHui  = SQRT( 2._wp * zomega * rdt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) &
-                               - v_il(ji,jj,jl) / a_i(ji,jj,jl)
-                        zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
+                  !----------------------------------------------------------------
+                  ! Lid melting: floating upper ice layer melts in whole or part
+                  !----------------------------------------------------------------
+
+                  zdvice = MIN( zv_pnd(ji,jj,jl), v_il(ji,jj,jl) )
+
+                  v_il (ji,jj,jl) = v_il (ji,jj,jl)   - zdvice
+
+                  IF ( v_il(ji,jj,jl) < epsi10 ) THEN
+                     ! ice lid melted
+                     v_il(ji,jj,jl)   = 0._wp
+                  ENDIF
+
+                  diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice   ! diag
+
+                  !----------------------------------------------------------------
+                  ! Freeze lid 
+                  !----------------------------------------------------------------
+
+                  ! The temperature of the lid's top surface is warmer than the ice top temperature
+                  ! by about 2 degrees (warmed by the melt pond under the lid)
+                  lid_Ttop = t_su(ji,jj,jl) + 2.0_wp
+
+                  ! differential growth of base of surface floating ice layer
+                  zdTice = MAX( - ( lid_Ttop - zTp ) , 0._wp ) ! > 0   
+                  zomega = rcnd_i * zdTice / zrhoi_L
+                  zdvice = SQRT( 2._wp * zomega * rdt_ice * a_ip(ji,jj,jl)**2 + v_il(ji,jj,jl)**2 ) &
+                         - v_il(ji,jj,jl)
+                  zdvice = min( zdvice , v_ip(ji,jj,jl) )
+
+                  IF ( zdvice > epsi10 ) THEN
+                     v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice
+                     v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice
+
+                     ! Possibly this check might not be needed to avoid divide by zero errors as at
+                     ! the end of ice_thd_pnd_area all ponds with an area less than epsi10 were zerod anyway
+                     IF ( a_ip(ji,jj,jl) > epsi10 ) THEN
+                        h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
+                     ENDIF
+
+                     diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice    ! diag
+
+                  ENDIF
                   
-                        IF ( zdvice > epsi10 ) THEN
-                           v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice
-                           v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice
-!                          zvolp(ji,jj)    = zvolp(ji,jj)     - zdvice
-!                          zfpond(ji,jj)   = zfpond(ji,jj)    - zdvice
-                           h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
-                           
-                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice    ! diag
-                           
-                        ENDIF
-                  
-                     ENDIF ! Tsfcn(i,j,n)
-
                   !----------------------------------------------------------------
-                  ! Freeze new lids
+                  ! Remove ponds if lids are much larger than ponds
                   !----------------------------------------------------------------
-                  !  upper ice layer begins to form
-                  ! note: albedo does not change
-
-                  ELSE ! v_il < epsi10
-                    
-                     ! thickness of newly formed ice
-                     ! the surface temperature of a meltpond is the same as that
-                     ! of the ice underneath (0C), and the thermodynamic surface 
-                     ! flux is the same
-                     
-                     !!! we need net surface energy flux, excluding conduction
-                     !!! fsurf is summed over categories in CICE
-                     !!! we have the category-dependent flux, let us use it ?
-                     zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl)                     
-                     zdHui  = MAX ( -zfsurf * rdt_ice / zrhoi_L , 0._wp )
-                     zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) )
-
-                     IF ( zdvice > epsi10 ) THEN
-
-                        v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice
-                        v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice
-                        diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice      ! diag
-
-!                       zvolp(ji,jj)     = zvolp(ji,jj)     - zdvice
-!                       zfpond(ji,jj)    = zfpond(ji,jj)    - zdvice
-
-                        h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) ! MV - in principle, this is useless as h_ip is computed in icevar
+
+                  IF ( a_ip(ji,jj,jl) > epsi10) THEN
+                     h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl)
+                     IF ( h_il(ji,jj,jl) > h_ip(ji,jj,jl) * 10._wp ) THEN
+                        a_ip(ji,jj,jl)      = 0._wp
+                        h_ip(ji,jj,jl)      = 0._wp
+                        v_ip(ji,jj,jl)      = 0._wp
+                        h_il(ji,jj,jl)      = 0._wp
+                        v_il(ji,jj,jl)      = 0._wp
                      ENDIF
-               
-                  ENDIF  ! v_il
+                  ENDIF
             
                END DO ! jl
