Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8597)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8598)
@@ -35,5 +35,4 @@
    USE icedyn         ! sea-ice: dynamics
    USE icethd         ! sea-ice: thermodynamics
-   USE limmp          ! sea-ice: melt ponds
    USE icecor         ! sea-ice: corrections
    USE iceupdate      ! sea-ice: sea surface boundary condition update
@@ -170,5 +169,4 @@
          IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics      
          !
-                                        CALL lim_mp( kt )             ! -- Melt ponds
          !
          IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections
@@ -238,6 +236,7 @@
       CALL ice_itd_init                ! ice thickness distribution initialization
       !
-      CALL lim_mp_init                 ! set melt ponds parameters (clem: important to be located here)
-      !
+      IF( ln_icethd ) THEN
+         CALL ice_thd_init             ! set ice thermodynics parameters (clem: important to call it first for melt ponds)
+      ENDIF   
       !                                ! Initial sea-ice state
       IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst
@@ -255,8 +254,4 @@
          CALL ice_dyn_init             ! set ice dynamics parameters
       ENDIF
-      !
-      IF( ln_icethd ) THEN
-         CALL ice_thd_init             ! set ice thermodynics parameters
-      ENDIF   
       !
       CALL ice_update_init             ! ice surface boundary condition
@@ -398,7 +393,5 @@
       wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp
       wfx_snw_sni(:,:) = 0._wp 
-      ! MV MP 2016
       wfx_pnd(:,:) = 0._wp
-      ! END MV MP 2016
 
       hfx_thd(:,:) = 0._wp   ;
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8597)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8598)
@@ -34,4 +34,5 @@
    USE icethd_ent     ! sea-ice: enthalpy redistribution
    USE icethd_do      ! sea-ice: growth in open water
+   USE icethd_pnd     ! sea-ice: melt ponds
    USE iceitd         ! sea-ice: remapping thickness distribution
    USE icetab         ! sea-ice: 1D <==> 2D transformation
@@ -86,5 +87,5 @@
       !!             - call ice_thd_rem  for remapping thickness distribution
       !!             - call ice_thd_do   for ice growth in leads
-      !!---------------------------------------------------------------------
+      !!-------------------------------------------------------------------
       INTEGER, INTENT(in) :: kt    ! number of iteration
       !
@@ -272,4 +273,6 @@
       !
       IF( ln_icedO )       CALL ice_thd_do                 ! --- frazil ice growing in leads --- !
+      !
+                           CALL ice_thd_pnd( kt )          ! --- Melt ponds --- !
       !
       ! controls
@@ -530,5 +533,5 @@
 
    SUBROUTINE ice_thd_init
-      !!-----------------------------------------------------------------------
+      !!-------------------------------------------------------------------
       !!                   ***  ROUTINE ice_thd_init *** 
       !!                 
@@ -570,4 +573,5 @@
       IF( ln_icedO )   CALL ice_thd_do_init    ! set ice growth in open water parameters
                        CALL ice_thd_sal_init   ! set ice salinity parameters
+                       CALL ice_thd_pnd_init   ! set melt ponds parameters
       !
       IF( ln_icedS .AND. nn_icesal == 1 ) THEN
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90	(revision 8598)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_pnd.F90	(revision 8598)
@@ -0,0 +1,324 @@
+MODULE icethd_pnd 
+   !!======================================================================
+   !!                     ***  MODULE  icethd_pnd   ***
+   !!   Melt ponds
+   !!======================================================================
+   !! history :       ! Original code by Daniela Flocco and Adrian Turner
+   !!            1.0  ! 2012    (O. Lecomte) Adaptation for scientific tests (NEMO3.1)
+   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3' :                                 LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   ice_thd_pnd_init      : some initialization and namelist read
+   !!   ice_thd_pnd           : main calling routine
+   !!----------------------------------------------------------------------
+   USE phycst           ! physical constants
+   USE dom_oce          ! ocean space and time domain
+   USE ice              ! LIM-3 variables
+   !
+   USE lbclnk           ! lateral boundary conditions - MPP exchanges
+   USE lib_mpp          ! MPP library
+   USE in_out_manager   ! I/O manager
+   USE lib_fortran      ! glob_sum
+   USE timing           ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_thd_pnd_init    ! routine called by icestp.F90
+   PUBLIC   ice_thd_pnd         ! routine called by icestp.F90
+
+   INTEGER ::              nice_pnd   ! choice of the type of pond scheme
+   !                                               ! associated indices:
+   INTEGER, PARAMETER ::   np_pndNO  = 0   ! No pond scheme
+   INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant pond scheme
+   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012)
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
+   !! $Id: limdyn.F90 6994 2016-10-05 13:07:10Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_thd_pnd( kt )
+      !!-------------------------------------------------------------------
+      !!               ***  ROUTINE ice_thd_pnd   ***
+      !!               
+      !! ** Purpose :   change melt pond fraction
+      !!                
+      !! ** Method  :   brutal force
+      !!
+      !! ** Action  : - 
+      !!              - 
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt    ! number of iteration
+      INTEGER ::   ji, jj, jl     ! dummy loop indices
+      !!-------------------------------------------------------------------
+
+      IF( nn_timing == 1 )  CALL timing_start('ice_thd_pnd')
+
+      SELECT CASE ( nice_pnd )
+
+      CASE (np_pndCST)
+         !                             !-------------------------------!
+         CALL pnd_CST                  ! Constant melt ponds           !
+         !                             !-------------------------------!
+      CASE (np_pndH12)
+         !                             !-------------------------------!
+         CALL pnd_H12                  ! Holland et al 2012 melt ponds !
+         !                             !-------------------------------!
+      END SELECT
+
+      IF( nn_timing == 1 )  CALL timing_stop('ice_thd_pnd')
+
+   END SUBROUTINE ice_thd_pnd 
+
+   SUBROUTINE pnd_CST 
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE pnd_CST  ***
+      !!
+      !! ** Purpose    : Compute melt pond evolution
+      !!
+      !! ** Method     : Melt pond fraction and thickness are prescribed 
+      !!                 to non-zero values when t_su = 0C
+      !!
+      !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
+      !!                
+      !! ** Note       : Coupling with such melt ponds is only radiative
+      !!                 Advection, ridging, rafting... are bypassed
+      !!
+      !! ** References : Bush, G.W., and Trump, D.J. (2017)
+      !!    
+      !!-------------------------------------------------------------------
+      
+      WHERE ( a_i(:,:,:) > 0._wp .AND. t_su(:,:,:) >= rt0 ) 
+         a_ip_frac = rn_apnd
+         h_ip      = rn_hpnd    
+         v_ip      = a_ip_frac * a_i * h_ip 
+         a_ip      = a_ip_frac * a_i
+      ELSEWHERE
+         a_ip      = 0._wp
+         h_ip      = 0._wp
+         v_ip      = 0._wp
+         a_ip_frac = 0._wp
+      END WHERE
+      
+   END SUBROUTINE pnd_CST
+
+   SUBROUTINE pnd_H12
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE pnd_H12  ***
+      !!
+      !! ** Purpose    : Compute melt pond evolution
+      !!
+      !! ** Method     : Empirical method. A fraction of meltwater is accumulated 
+      !!                 in pond volume. It is then released exponentially when
+      !!                 surface is freezing.
+      !!
+      !! ** Tunable parameters : (no real expertise yet, ideas?)
+      !! 
+      !! ** Note       : Stolen from CICE for quick test of the melt pond
+      !!                 radiation and freshwater interfaces
+      !!                 Coupling can be radiative AND freshwater
+      !!                 Advection, ridging, rafting are called
+      !!
+      !! ** References : Holland, M. M. et al (J Clim 2012)
+      !!    
+      !!-------------------------------------------------------------------
+      
+      INTEGER, DIMENSION(jpij)         ::   indxi             ! compressed indices for cells with ice melting
+      INTEGER, DIMENSION(jpij)         ::   indxj             !
+
+      REAL(wp), DIMENSION(jpi,jpj)     ::   zwfx_mlw          ! available meltwater for melt ponding
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zrfrac            ! fraction of available meltwater retained for melt ponding
+
+      REAL(wp), PARAMETER ::   zrmin  = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
+      REAL(wp), PARAMETER ::   zrmax  = 0.70_wp  ! maximum   ''           ''       ''        ''            ''
+      REAL(wp), PARAMETER ::   zrexp  = 0.01_wp  ! rate constant to refreeze melt ponds
+      REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp ! pond aspect ratio
+
+      REAL(wp) ::   zhi               ! dummy ice thickness
+      REAL(wp) ::   zhs               ! dummy snow depth
+      REAL(wp) ::   zTp, z1_Tp        ! reference temperature
+      REAL(wp) ::   zdTs              ! dummy temperature difference
+      REAL(wp) ::   z1_rhofw          ! inverse freshwater density
+      REAL(wp) ::   z1_zpnd_aspect    ! inverse pond aspect ratio
+      REAL(wp) ::   zvpold            ! dummy pond volume
+
+      INTEGER  ::   ji, jj, jl, ij    ! loop indices
+      INTEGER  ::   icells            ! size of dummy array
+      !!-------------------------------------------------------------------
+      z1_rhofw       = 1. / rhofw 
+      z1_zpnd_aspect = 1. / zpnd_aspect
+      zTp            = -2. 
+      z1_Tp          = 1._wp / zTp 
+
+      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_sum(:,:), 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
+      ! NB: I now changed to wfx_snw_sum, this may fix the problem. 
+      ! We should check
+
+      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
+         END DO
+
+         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_fwb )   wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + v_ip(ji,jj,jl)  !--- Give freshwater to the ocean
+
+            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 * z1_Tp )
+
+               !--- Dump meltwater due to refreezing ( of course this is wrong
+               !--- but this parameterization is too simple )
+               IF ( ln_pnd_fwb ) THEN
+                  wfx_pnd(ji,jj)   = wfx_pnd(ji,jj) + rhofw * ( v_ip(ji,jj,jl) - zvpold ) * r1_rdtice
+               ENDIF
+               
+               a_ip_frac(ji,jj,jl) = MIN( 1._wp , SQRT( v_ip(ji,jj,jl) * z1_zpnd_aspect / a_i(ji,jj,jl) ) )
+               !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative
+
+               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)
+
+            ENDIF
+
+         END DO
+
+      END DO
+
+      !--- Remove retained meltwater from surface fluxes 
+
+      IF ( ln_pnd_fwb ) THEN
+
+         wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) ) 
+         wfx_sum(:,:)     = wfx_sum(:,:)     *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) )
+
+      ENDIF
+
+   END SUBROUTINE pnd_H12
+
+   SUBROUTINE ice_thd_pnd_init 
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_thd_pnd_init   ***
+      !!
+      !! ** Purpose : Physical constants and parameters linked to melt ponds
+      !!              over sea ice
+      !!
+      !! ** Method  :  Read the namthd_pnd  namelist and check the melt pond  
+      !!               parameter values called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namelist namthd_pnd  
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ios, ioptio                 ! Local integer output status for namelist read
+      NAMELIST/namthd_pnd/  ln_pnd_H12, ln_pnd_fwb, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
+      !!-------------------------------------------------------------------
+
+      REWIND( numnam_ice_ref )              ! Namelist namthd_pnd  in reference namelist : Melt Ponds  
+      READ  ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd  in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namthd_pnd  in configuration namelist : Melt Ponds
+      READ  ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namthd_pnd )
+      
+      IF(lwp) THEN                        ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
+         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '   Namelist namicethd_pnd:'
+         WRITE(numout,*) '      Evolutive melt pond fraction and depth (Holland et al 2012)  ln_pnd_H12 = ', ln_pnd_H12
+         WRITE(numout,*) '         Melt ponds store fresh water or not                       ln_pnd_fwb = ', ln_pnd_fwb
+         WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_Cst = ', ln_pnd_CST
+         WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
+         WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
+         WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
+      ENDIF
+      !
+      !                             !== set the choice of ice pond scheme ==!
+      ioptio = 0
+                                                            nice_pnd = np_pndNO
+      IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
+      IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
+      IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' )
+
+      SELECT CASE( nice_pnd )
+      CASE( np_pndNO )         
+         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when no ponds' ) ; ENDIF
+         IF(ln_pnd_alb) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
+      CASE( np_pndCST)
+         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when ln_pnd_CST=true' ) ; ENDIF
+      END SELECT
+      !
+   END SUBROUTINE ice_thd_pnd_init
+   
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option          Empty module          NO ESIM sea-ice model
+   !!----------------------------------------------------------------------
+#endif 
+
+   !!======================================================================
+END MODULE icethd_pnd 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limmp.F90	(revision 8597)
+++ 	(revision )
@@ -1,327 +1,0 @@
-MODULE limmp 
-   !!======================================================================
-   !!                     ***  MODULE  limmp   ***
-   !!   Melt ponds
-   !!======================================================================
-   !! history :       ! Original code by Daniela Flocco and Adrian Turner
-   !!            1.0  ! 2012    (O. Lecomte) Adaptation for scientific tests (NEMO3.1)
-   !!            2.0  ! 2017    (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO3.6
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                 LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   lim_mp_init      : some initialization and namelist read
-   !!   lim_mp           : main calling routine
-   !!----------------------------------------------------------------------
-   USE phycst           ! physical constants
-   USE dom_oce          ! ocean space and time domain
-   USE ice              ! LIM-3 variables
-   !
-   USE lbclnk           ! lateral boundary conditions - MPP exchanges
-   USE lib_mpp          ! MPP library
-   USE in_out_manager   ! I/O manager
-   USE lib_fortran      ! glob_sum
-   USE timing           ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   lim_mp_init    ! routine called by icestp.F90
-   PUBLIC   lim_mp         ! routine called by icestp.F90
-
-   INTEGER ::              nice_pnd   ! choice of the type of pond scheme
-   !                                               ! associated indices:
-   INTEGER, PARAMETER ::   np_pndNO  = 0   ! No pond scheme
-   INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant pond scheme
-   INTEGER, PARAMETER ::   np_pndH12 = 2   ! Evolutive pond scheme (Holland et al. 2012)
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
-   !! $Id: limdyn.F90 6994 2016-10-05 13:07:10Z clem $
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE lim_mp( kt )
-      !!-------------------------------------------------------------------
-      !!               ***  ROUTINE lim_mp   ***
-      !!               
-      !! ** Purpose :   change melt pond fraction
-      !!                
-      !! ** Method  :   brutal force
-      !!
-      !! ** Action  : - 
-      !!              - 
-      !!------------------------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt    ! number of iteration
-      INTEGER ::   ji, jj, jl     ! dummy loop indices
-      !!-------------------------------------------------------------------
-
-      IF( nn_timing == 1 )  CALL timing_start('lim_mp')
-
-      SELECT CASE ( nice_pnd )
-
-      CASE (np_pndCST)
-         !                             !------------------------------!
-         CALL lim_mp_CST               ! Constant melt ponds          !
-         !                             !------------------------------!
-      CASE (np_pndH12)
-         !                             !------------------------------!
-         CALL lim_mp_H12               ! Holland et al 2012 melt ponds!
-         !                             !------------------------------!
-      END SELECT
-
-      IF( nn_timing == 1 )  CALL timing_stop('lim_mp')
-
-   END SUBROUTINE lim_mp 
-
-   SUBROUTINE lim_mp_CST 
-       !!-------------------------------------------------------------------
-       !!                ***  ROUTINE lim_mp_CST  ***
-       !!
-       !! ** Purpose    : Compute melt pond evolution
-       !!
-       !! ** Method     : Melt pond fraction and thickness are prescribed 
-       !!                 to non-zero values when t_su = 0C
-       !!
-       !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd)
-       !!                
-       !! ** Note       : Coupling with such melt ponds is only radiative
-       !!                 Advection, ridging, rafting... are bypassed
-       !!
-       !! ** References : Bush, G.W., and Trump, D.J. (2017)
-       !!    
-       !!-------------------------------------------------------------------
-
-       WHERE ( ( a_i > epsi10 ) .AND. ( t_su >= rt0-epsi06 ) ) 
-!!       WHERE ( ( a_i > 0._wp ) .AND. ( t_su >= rt0 ) ) 
-          a_ip_frac = rn_apnd
-          h_ip      = rn_hpnd    
-          v_ip      = a_ip_frac * a_i * h_ip 
-          a_ip      = a_ip_frac * a_i
-       ELSE WHERE
-          a_ip      = 0._wp
-          h_ip      = 0._wp
-          v_ip      = 0._wp
-          a_ip_frac = 0._wp
-       END WHERE
-
-       wfx_pnd(:,:) = 0._wp
-
-   END SUBROUTINE lim_mp_CST
-
-   SUBROUTINE lim_mp_H12
-       !!-------------------------------------------------------------------
-       !!                ***  ROUTINE lim_mp_H12  ***
-       !!
-       !! ** Purpose    : Compute melt pond evolution
-       !!
-       !! ** Method     : Empirical method. A fraction of meltwater is accumulated 
-       !!                 in pond volume. It is then released exponentially when
-       !!                 surface is freezing.
-       !!
-       !! ** Tunable parameters : (no real expertise yet, ideas?)
-       !! 
-       !! ** Note       : Stolen from CICE for quick test of the melt pond
-       !!                 radiation and freshwater interfaces
-       !!                 Coupling can be radiative AND freshwater
-       !!                 Advection, ridging, rafting are called
-       !!
-       !! ** References : Holland, M. M. et al (J Clim 2012)
-       !!    
-       !!-------------------------------------------------------------------
-
-       INTEGER, DIMENSION(jpij)         ::   indxi             ! compressed indices for cells with ice melting
-       INTEGER, DIMENSION(jpij)         ::   indxj             !
-
-       REAL(wp), DIMENSION(jpi,jpj)     ::   zwfx_mlw          ! available meltwater for melt ponding
-       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zrfrac            ! fraction of available meltwater retained for melt ponding
-
-       REAL(wp), PARAMETER ::   zrmin  = 0.15_wp  ! minimum fraction of available meltwater retained for melt ponding
-       REAL(wp), PARAMETER ::   zrmax  = 0.70_wp  ! maximum   ''           ''       ''        ''            ''
-       REAL(wp), PARAMETER ::   zrexp  = 0.01_wp  ! rate constant to refreeze melt ponds
-       REAL(wp), PARAMETER ::   zpnd_aspect = 0.8_wp ! pond aspect ratio
-
-       REAL(wp) ::   zhi               ! dummy ice thickness
-       REAL(wp) ::   zhs               ! dummy snow depth
-       REAL(wp) ::   zTp               ! reference temperature
-       REAL(wp) ::   zdTs              ! dummy temperature difference
-       REAL(wp) ::   z1_rhofw          ! inverse freshwater density
-       REAL(wp) ::   z1_zpnd_aspect    ! inverse pond aspect ratio
-       REAL(wp) ::   zvpold            ! dummy pond volume
-
-       INTEGER  ::   ji, jj, jl, ij    ! loop indices
-       INTEGER  ::   icells            ! size of dummy array
-       !!-------------------------------------------------------------------
-        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_sum(:,:), 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
-                ! NB: I now changed to wfx_snw_sum, this may fix the problem. 
-                ! We should check
- 
-        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
-          END DO
- 
-          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_fwb ) & !--- 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_fwb ) &
-                   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) ) )
-                   !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative
- 
-                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)
-             
-             ENDIF
- 
-           END DO
- 
-        END DO ! jpl
- 
-        !--- Remove retained meltwater from surface fluxes 
- 
-        IF ( ln_pnd_fwb ) THEN
- 
-           wfx_snw_sum(:,:) = wfx_snw_sum(:,:) *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) ) 
-           wfx_sum(:,:)     = wfx_sum(:,:)     *  ( 1. - zrmin - ( zrmax - zrmin ) * SUM( a_i(:,:,:), dim=3 ) )
-
-        ENDIF
-
-   END SUBROUTINE lim_mp_H12
-
-   SUBROUTINE lim_mp_init 
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE lim_mp_init   ***
-      !!
-      !! ** Purpose : Physical constants and parameters linked to melt ponds
-      !!      over sea ice
-      !!
-      !! ** Method  :  Read the nammp  namelist and check the melt pond  
-      !!       parameter values called at the first timestep (nit000)
-      !!
-      !! ** input   :   Namelist nammp  
-      !!-------------------------------------------------------------------
-      INTEGER  ::   ios, ioptio                 ! Local integer output status for namelist read
-      NAMELIST/nammp/  ln_pnd_H12, ln_pnd_fwb, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb
-      !!-------------------------------------------------------------------
-
-      REWIND( numnam_ice_ref )              ! Namelist nammp  in reference namelist : Melt Ponds  
-      READ  ( numnam_ice_ref, nammp, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp  in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )              ! Namelist nammp  in configuration namelist : Melt Ponds
-      READ  ( numnam_ice_cfg, nammp, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'nammp in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, nammp )
-      
-      IF(lwp) THEN                        ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
-         WRITE(numout,*) '~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) '   Namelist namicethd_pnd:'
-         WRITE(numout,*) '      Evolutive melt pond fraction and depth (Holland et al 2012)  ln_pnd_H12 = ', ln_pnd_H12
-         WRITE(numout,*) '         Melt ponds store fresh water or not                       ln_pnd_fwb = ', ln_pnd_fwb
-         WRITE(numout,*) '      Prescribed melt pond fraction and depth                      ln_pnd_Cst = ', ln_pnd_CST
-         WRITE(numout,*) '         Prescribed pond fraction                                  rn_apnd    = ', rn_apnd
-         WRITE(numout,*) '         Prescribed pond depth                                     rn_hpnd    = ', rn_hpnd
-         WRITE(numout,*) '      Melt ponds affect albedo or not                              ln_pnd_alb = ', ln_pnd_alb
-      ENDIF
-      !
-      !                             !== set the choice of ice pond scheme ==!
-      ioptio = 0
-                                                            nice_pnd = np_pndNO
-      IF( ln_pnd_CST ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF
-      IF( ln_pnd_H12 ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndH12    ;   ENDIF
-      IF( ioptio > 1 )   CALL ctl_stop( 'ice_thd_pnd_init: choose one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' )
-
-      SELECT CASE( nice_pnd )
-      CASE( np_pndNO )         
-         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when no ponds' ) ; ENDIF
-         IF(ln_pnd_alb) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
-      CASE( np_pndCST)
-         IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when ln_pnd_CST=true' ) ; ENDIF
-      END SELECT
-      !
-   END SUBROUTINE lim_mp_init
-   
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option          Empty module          NO ESIM sea-ice model
-   !!----------------------------------------------------------------------
-#endif 
-
-   !!======================================================================
-END MODULE limmp 
