Index: branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90
===================================================================
--- branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8125)
+++ branches/2016/dev_r6859_LIM3_meltponds/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8142)
@@ -6,5 +6,7 @@
    !! history :       ! Original code by Daniela Flocco and Adrian Turner
    !!            1.0  ! 2012    (O. Lecomte) Adaptation for scientific tests (NEMO3.1)
-   !!            2.0  ! 2016    (O. Lecomte, C. Rousset, M. Vancoppenolle) Implementation in NEMO3.6
+   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6
+   !!                 ! NB: Only lim_mp_cstt and lim_mp_cesm work in this
+   !!                       version
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -137,5 +139,4 @@
    
 
-
    SUBROUTINE lim_mp( kt )
       !!-------------------------------------------------------------------
@@ -153,9 +154,5 @@
       INTEGER  ::   ji, jj, jl     ! dummy loop indices
 
-!     REAL(wp), POINTER, DIMENSION(:,:) ::   zfsurf  ! surface heat flux(obsolete, should be droped)
-      !
       !!-------------------------------------------------------------------
-
-      IF( nn_timing == 1 )  CALL timing_start('limthd')
 
       IF( nn_timing == 1 )  CALL timing_start('lim_mp')
@@ -180,9 +177,5 @@
       END SELECT
 
-      ! we should probably not aggregate here since we do it in lim_var_agg
-      ! before output, unless we need the total volume and faction else where
-
-      ! we should also make sure a_ip and v_ip are properly updated at the end
-      ! of the routine
+      IF( nn_timing == 1 )  CALL timing_stop('lim_mp')
 
    END SUBROUTINE lim_mp 
@@ -237,5 +230,5 @@
        !!                 surface is freezing.
        !!
-       !! ** Tunable parameters : (no expertise yet)
+       !! ** Tunable parameters : (no real expertise yet, ideas?)
        !! 
        !! ** Note       : Stolen from CICE for quick test of the melt pond
@@ -272,108 +265,122 @@
        !!-------------------------------------------------------------------
 
-       CALL wrk_alloc( jpi*jpj, indxi, indxj)
-       CALL wrk_alloc( jpi,jpj,     zwfx_mlw )
-       CALL wrk_alloc( jpi,jpj,jpl, zrfrac   )
-
-       z1_rhofw       = 1. / rhofw 
-       z1_zpnd_aspect = 1. / zpnd_aspect
-       zTp            = -2. 
-
-       !------------------------------------------------------------------
-       ! Available melt water for melt ponding and corresponding fraction
-       !------------------------------------------------------------------
-
-       zwfx_mlw(:,:) = wfx_sum(:,:) + wfx_snw(:,:)        ! available meltwater for melt ponding
-
-       zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:)  
-
-       DO jl = 1, jpl   
-
-          ! v_ip(:,:,jl) ! Initialize things
-          ! a_ip(:,:,jl)
-          ! volpn(:,:) = hpnd(:,:) * apnd(:,:) * aicen(:,:)
-
-          !------------------------------------------------------------------------------
-          ! Identify grid cells where ponds should be updated (can probably be improved)
-          !------------------------------------------------------------------------------
-
-          indxi(:) = 0
-          indxj(:) = 0
-          icells   = 0
-
-          DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF ( a_i(ji,jj,jl) > epsi10 ) THEN
-                  icells = icells + 1
-                  indxi(icells) = ji
-                  indxj(icells) = jj
-               ENDIF
-            END DO                 ! ji
-         END DO                    ! jj
-
-         DO ij = 1, icells
-
-            ji = indxi(ij)
-            jj = indxj(ij)
-
-            zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
-            zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)
-
-            IF ( zhi < rn_himin) THEN   !--- Remove ponds on thin ice if ice is too thin
-
-               a_ip(ji,jj,jl)      = 0._wp                               !--- Dump ponds
-               v_ip(ji,jj,jl)      = 0._wp
-               a_ip_frac(ji,jj,jl) = 0._wp
-               h_ip(ji,jj,jl)      = 0._wp
-
-               IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean
-                  wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + v_ip(ji,jj,jl) 
-
-
-            ELSE                        !--- Update pond characteristics
-
-               !--- Add retained melt water to melt ponds
-               v_ip(ji,jj,jl)      = v_ip(ji,jj,jl) + zrfrac(ji,jj,jl) * z1_rhofw * zwfx_mlw(ji,jj) * a_i(ji,jj,jl) * rdt_ice
-
-               !--- Shrink pond due to refreezing
-               zdTs                = MAX ( zTp - t_su(ji,jj,jl) + rt0 , 0. )
-               
-               zvpold              = v_ip(ji,jj,jl)
-
-               v_ip(ji,jj,jl)      = v_ip(ji,jj,jl) * EXP( zrexp * zdTs / zTp )
-
-               !--- Dump meltwater due to refreezing ( of course this is wrong
-               !--- but this parameterization is too simple )
-               IF ( ln_pnd_fw ) &
-                  wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice
-
-               a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) )
-
-               h_ip(ji,jj,jl)      = zpnd_aspect * a_ip_frac(ji,jj,jl)
-
-               a_ip(ji,jj,jl)      = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl)
-
-            !-----------------------------------------------------------
-            ! Limit pond depth
-            !-----------------------------------------------------------
-            ! hpondn = min(hpondn, dpthhi*hi)
-
-            !--- Give freshwater to the ocean ?
-
-            ENDIF
-
-          END DO
-
-       END DO ! jpl
-
-       !--- Remove retained meltwater from surface fluxes 
-
-       IF ( ln_pnd_fw ) THEN
-
-           wfx_snw(:,:) = wfx_snw(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
-
-           wfx_sum(:,:) = wfx_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
-
-       ENDIF
+        CALL wrk_alloc( jpi*jpj, indxi, indxj)
+        CALL wrk_alloc( jpi,jpj,     zwfx_mlw )
+        CALL wrk_alloc( jpi,jpj,jpl, zrfrac   )
+
+        z1_rhofw       = 1. / rhofw 
+        z1_zpnd_aspect = 1. / zpnd_aspect
+        zTp            = -2. 
+
+        a_ip_frac(:,:,:) = 0._wp
+        h_ip     (:,:,:) = 0._wp
+ 
+        !------------------------------------------------------------------
+        ! Available melt water for melt ponding and corresponding fraction
+        !------------------------------------------------------------------
+ 
+        zwfx_mlw(:,:) = MAX( wfx_sum(:,:) + wfx_snw(:,:), 0._wp )        ! available meltwater for melt ponding
+
+                ! NB: zwfx_mlw can be slightly negative for very small values (why?)
+                ! This can in some occasions give negative
+                ! v_ip in the first category, which then gives crazy pond
+                ! fractions and crashes the code as soon as the melt-pond
+                ! radiative coupling is activated
+                ! if we understand and remove why wfx_sum or wfx_snw could be
+                ! negative, then, we can remove the MAX
+ 
+        zrfrac(:,:,:) = zrmin + ( zrmax - zrmin ) * a_i(:,:,:)  
+ 
+        DO jl = 1, jpl   
+
+           !------------------------------------------------------------------------------
+           ! Identify grid cells where ponds should be updated (can probably be improved)
+           !------------------------------------------------------------------------------
+ 
+           indxi(:) = 0
+           indxj(:) = 0
+           icells   = 0
+ 
+           DO jj = 1, jpj
+             DO ji = 1, jpi
+                IF ( a_i(ji,jj,jl) > epsi10 ) THEN
+                   icells = icells + 1
+                   indxi(icells) = ji
+                   indxj(icells) = jj
+                ENDIF
+             END DO                 ! ji
+          END DO                    ! jj
+ 
+          DO ij = 1, icells
+ 
+             ji = indxi(ij)
+             jj = indxj(ij)
+ 
+             zhi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
+             zhs = v_s(ji,jj,jl) / a_i(ji,jj,jl)
+ 
+             IF ( zhi < rn_himin) THEN   !--- Remove ponds on thin ice if ice is too thin
+ 
+                a_ip(ji,jj,jl)      = 0._wp                               !--- Dump ponds
+                v_ip(ji,jj,jl)      = 0._wp
+                a_ip_frac(ji,jj,jl) = 0._wp
+                h_ip(ji,jj,jl)      = 0._wp
+ 
+                IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean
+                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + v_ip(ji,jj,jl) 
+ 
+ 
+             ELSE                        !--- Update pond characteristics
+ 
+                !--- Add retained melt water to melt ponds
+                ! v_ip should never be positive, otherwise code crashes
+                ! MV: as far as I saw, UM5 can create very small negative v_ip values
+                ! hence I added the max, which was not required with Prather (1 yr run)
+                v_ip(ji,jj,jl)      = MAX( v_ip(ji,jj,jl), 0._wp ) + zrfrac(ji,jj,jl) * z1_rhofw * zwfx_mlw(ji,jj) * a_i(ji,jj,jl) * rdt_ice
+ 
+                !--- Shrink pond due to refreezing
+                zdTs                = MAX ( zTp - t_su(ji,jj,jl) + rt0 , 0. )
+                
+                zvpold              = v_ip(ji,jj,jl)
+ 
+                v_ip(ji,jj,jl)      = v_ip(ji,jj,jl) * EXP( zrexp * zdTs / zTp )
+ 
+                !--- Dump meltwater due to refreezing ( of course this is wrong
+                !--- but this parameterization is too simple )
+                IF ( ln_pnd_fw ) &
+                   wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice
+ 
+                a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) )
+ 
+                h_ip(ji,jj,jl)      = zpnd_aspect * a_ip_frac(ji,jj,jl)
+ 
+                a_ip(ji,jj,jl)      = a_ip_frac(ji,jj,jl) * a_i(ji,jj,jl)
+ 
+             !-----------------------------------------------------------
+             ! Limit pond depth
+             !-----------------------------------------------------------
+             ! The original version has pond depth limitation, which I did not
+             ! keep here. Maybe we need it later on
+             !
+             
+             ENDIF
+ 
+           END DO
+ 
+        END DO ! jpl
+ 
+        !--- Remove retained meltwater from surface fluxes 
+ 
+        IF ( ln_pnd_fw ) THEN
+ 
+            wfx_snw(:,:) = wfx_snw(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
+ 
+            wfx_sum(:,:) = wfx_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * at_i(:,:) )
+ 
+        ENDIF
+
+       CALL wrk_dealloc( jpi*jpj, indxi, indxj)
+       CALL wrk_dealloc( jpi,jpj,     zwfx_mlw )
+       CALL wrk_dealloc( jpi,jpj,jpl, zrfrac   )
 
    END SUBROUTINE lim_mp_cesm
