Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5047)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5048)
@@ -83,5 +83,6 @@
       !!
       INTEGER  :: ji, jj, jk, jl   ! dummy loop indices
-      INTEGER  :: nbpb             ! nb of icy pts for thermo. cal.
+      INTEGER  :: nbpb             ! nb of icy pts for vertical thermo calculations
+      INTEGER  :: nbplm            ! nb of icy pts for lateral melting calculations (mono-cat)
       INTEGER  :: ii, ij           ! temporary dummy loop index
       REAL(wp) :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04)
@@ -434,4 +435,8 @@
               CALL tab_1d_2d( nbpb, hfx_res       , npb, hfx_res_1d(1:nbpb)   , jpi, jpj )
               CALL tab_1d_2d( nbpb, hfx_err_rem   , npb, hfx_err_rem_1d(1:nbpb)   , jpi, jpj )
+
+            IF ( ( ( nn_monocat .EQ. 1 ) .OR. ( nn_monocat .EQ. 4 ) ) .AND. ( jpl == 1 ) ) THEN
+              CALL tab_1d_2d( nbpb, dh_i_melt(:,:,jl) , npb, dh_i_melt_1d(1:nbpb) , jpi, jpj )
+            ENDIF
             !
               CALL tab_1d_2d( nbpb, qns_ice(:,:,jl), npb, qns_ice_1d(1:nbpb) , jpi, jpj)
@@ -471,4 +476,52 @@
       !----------------------------------
       CALL lim_var_eqv2glo
+
+      !----------------------------------
+      ! 5.X) Lateral melting
+      !----------------------------------
+      !!! declare dh_i_melt (ok), dh_i_melt_1d (ok), nbplm (ok), nplm (ok), zda(ok)
+
+      IF ( ( ( nn_monocat .EQ. 1 ) .OR. ( nn_monocat .EQ. 4 ) ) .AND. ( jpl == 1 ) ) THEN
+
+         WRITE(numout,*) ' Lateral melting ON '
+
+         ! select points where lateral melting occurs
+         jl = 1
+
+         nbplm = 0
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               IF ( ( dh_i_melt(ji,jj,jl)                  .LT.-epsi10 ) .AND.    &
+     &              ( ht_i(ji,jj,jl) - dh_i_melt(ji,jj,jl) .GT. epsi10 ) .AND.    & 
+     &              ( ht_i(ji,jj,jl)                       .GT. epsi10 ) ) THEN     
+                  nbplm   = nbplm  + 1
+                  nplm(nbplm) = (jj - 1) * jpi + ji
+               ENDIF
+            END DO
+         END DO
+
+         IF( nbplm > 0 ) THEN  ! If there is no net melting, do nothing
+
+            ! Move to 1D arrays
+            !-------------------------
+
+            CALL tab_2d_1d( nbplm, a_i_1d      (1:nbplm), a_i(:,:,jl)       , jpi, jpj, nplm(1:nbplm) )
+            CALL tab_2d_1d( nbplm, ht_i_1d     (1:nbplm), ht_i(:,:,jl)      , jpi, jpj, nplm(1:nbplm) )
+            CALL tab_2d_1d( nbplm, dh_i_melt_1d(1:nbplm), dh_i_melt(:,:,jl) , jpi, jpj, nplm(1:nbplm) )
+
+            ! Compute lateral melting (dA = A/2h dh )
+            DO ji = 1, nbplm
+               zda        = a_i_1d(ji) * dh_i_melt_1d(ji) / ( 2._wp * ht_i_1d(ji) )
+               a_i_1d(ji) = a_i_1d(ji) + zda
+            END DO
+
+            ! Move back to 2D arrays
+            !-------------------------
+            CALL tab_1d_2d( nbplm, a_i (:,:,jl)  , nplm, a_i_1d     (1:nbplm)   , jpi, jpj )
+            at_i(:,:) = a_i(:,:,jl)
+
+         ENDIF
+
+      ENDIF
 
       !--------------------------------------------
@@ -563,7 +616,8 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicethd/ hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb,   &
+      NAMELIST/namicethd/ hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb,    &
          &                hiclim, parsub, betas,                          & 
-         &                kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi
+         &                kappa_i, nconv_i_thd, maxer_i_thd, thcon_i_swi, &
+         &                nn_monocat
       !!-------------------------------------------------------------------
       !
@@ -582,4 +636,9 @@
 902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicethd in configuration namelist', lwp )
       IF(lwm) WRITE ( numoni, namicethd )
+      !
+      IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN 
+         nn_monocat = 0
+         WRITE(numout, *) ' nn_monocat must be 0 in multi-category case '
+      ENDIF
 
       IF( lk_cpl .AND. parsub /= 0.0 )   CALL ctl_stop( 'In coupled mode, use parsub = 0. or send dqla' )
@@ -597,9 +656,10 @@
          WRITE(numout,*)'      switch for snow sublimation  (=1) or not (=0)           parsub       = ', parsub  
          WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          betas        = ', betas
-         WRITE(numout,*)'      extinction radiation parameter in sea ice (1.0)         kappa_i      = ', kappa_i
+         WRITE(numout,*)'      extinction radiation parameter in sea ice               kappa_i      = ', kappa_i
          WRITE(numout,*)'      maximal n. of iter. for heat diffusion computation      nconv_i_thd  = ', nconv_i_thd
          WRITE(numout,*)'      maximal err. on T for heat diffusion computation        maxer_i_thd  = ', maxer_i_thd
          WRITE(numout,*)'      switch for comp. of thermal conductivity in the ice     thcon_i_swi  = ', thcon_i_swi
          WRITE(numout,*)'      check heat conservation in the ice/snow                 con_i        = ', con_i
+         WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat
       ENDIF
       !
