Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/CONFIG/SHARED/field_def.xml
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 5047)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/CONFIG/SHARED/field_def.xml	(revision 5048)
@@ -297,4 +297,8 @@
          <field id="vfxsub"       long_name="snw sublimation"                                              unit="m/day"   />
          <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="m/day"   />
+
+         <field id="afxtot"       long_name="area tendency (total)"                                        unit="s-1"     />
+         <field id="afxdyn"       long_name="area tendency (dynamics)"                                     unit="s-1"     />
+         <field id="afxthd"       long_name="area tendency (thermo)"                                       unit="s-1"     />
 
          <field id="hfxsum"   long_name="heat fluxes causing surface ice melt"            unit="W/m2"  />
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 5047)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 5048)
@@ -2,7 +2,7 @@
 !! NEMO/LIM3 :  1 - dynamics/advection/thermo          (namicerun)
 !! namelists    2 - ice intialisation                  (namiceini)
-!!              3 - ice dynamic                        (namicedyn)
+!!              3 - ice dynamics                       (namicedyn)
 !!              4 - ice advection                      (namicetrp)
-!!              5 - thermodynamic                      (namicethd)
+!!              5 - thermodynamics                     (namicethd)
 !!              6 - ice salinity                       (namicesal)
 !!              7 - mechanical redistribution of ice   (namiceitdme)
@@ -17,86 +17,83 @@
    cn_icerst_out = "restart_ice"   !  suffix of ice restart name (output)
    ln_limdyn     = .true.          !  ice dynamics (T) or thermodynamics only (F)
-   amax          = 0.999           !  maximum ice concentration
-   ln_nicep      = .false.         !  Ice points output for debug (yes or no)
-   ln_limdiahsb  = .false.          !  check the heat and salt budgets (T) or not (F)
+   amax          = 0.999           !  maximum tolerated ice concentration
+   ln_nicep      = .false.         !  ice points output for debug (yes or no)
+   ln_limdiahsb  = .false.         !  check the heat and salt budgets (T) or not (F)
    ln_limdiaout  = .true.          !  output the heat and salt budgets (T) or not (F)
 /
 !-----------------------------------------------------------------------
-&namiceini     !   ice initialisation
+&namiceini     !   ice initialization
 !-----------------------------------------------------------------------
    ln_limini   = .false.   !  activate ice initialization (T) or not (F)
    thres_sst   =  0.0      !  threshold water temperature for initial sea ice
-   hts_ini_n   =  0.3      !  initial snow thickness in the north
-   hts_ini_s   =  0.3      !        "            "          south
-   hti_ini_n   =  3.0      !  initial ice thickness in the north
-   hti_ini_s   =  1.0      !        "            "         south
-   ati_ini_n   =  0.9      !  initial ice concentration in the north
-   ati_ini_s   =  0.9      !        "            "             south
-   smi_ini_n   =  6.3      !  initial ice salinity in the north
-   smi_ini_s   =  6.3      !        "            "    south
-   tmi_ini_n   =  270.     !  initial ice/snw temp in the north
-   tmi_ini_s   =  270.     !  initial ice/snw temp in the south
+   hts_ini_n   =  0.3      !  initial snow thickness (North)
+   hts_ini_s   =  0.3      !        "            "   (South)
+   hti_ini_n   =  3.0      !  initial ice thickness  (North)
+   hti_ini_s   =  1.0      !        "            "   (South)
+   ati_ini_n   =  0.9      !  initial ice concentration (North)
+   ati_ini_s   =  0.9      !        "            "      (South)
+   smi_ini_n   =  6.3      !  initial ice salinity   (North)
+   smi_ini_s   =  6.3      !        "            "   (South)
+   tmi_ini_n   =  270.     !  initial ice/snw temp   (North)
+   tmi_ini_s   =  270.     !  initial ice/snw temp   (South)
 /
 !-----------------------------------------------------------------------
-&namicedyn     !   ice dynamic
+&namicedyn     !   ice dynamics
 !-----------------------------------------------------------------------
-   cw          =   5.0e-03 !  drag coefficient for oceanic stress
-   pstar       =   2.0e+04 !  1st bulk-rheology parameter
+   cw          =   5.0e-03 !  ice-ocean drag coefficient
+   pstar       =   2.0e+04 !  1st bulk-rheology parameter (N/m)
    c_rhg       =  20.0     !  2nd bulk-rhelogy parameter
-   creepl      =   1.0e-12 !  creep limit
+   creepl      =   1.0e-12 !  creep limit (s-1)
    ecc         =   2.0     !  eccentricity of the elliptical yield curve
-   ahi0        = 350.e0    !  horizontal eddy diffusivity coefficient for sea-ice [m2/s] ! depends on resolution
-   nevp        = 120       !  number of iterations for subcycling in EVP
-   relast      = 0.333     !  ratio of elastic timescale over ice time step (1/3 if nevp=120 ; 1/9 if nevp=300)
-   hminrhg     =   0.001   !  ice volume (a*h in m) below which ice velocity equal ocean velocity
+   ahi0        = 350.e0    !  horizontal eddy diffusivity coefficient for sea-ice (m2/s), depends on resolution
+   nevp        = 120       !  number of EVP subcycles
+   relast      =   0.333   !  ratio of elastic timescale to ice time step (1/3 if nevp=120, 1/9 if nevp=300)
+   hminrhg     =   0.001   !  ice volume (a*h in m) below which ice velocity equals ocean velocity
 /
 !-----------------------------------------------------------------------
-&namicethd     !   ice thermodynamic
+&namicethd     !   ice thermodynamics
 !-----------------------------------------------------------------------
-   hiccrit     = 0.1       !  ice thickness for lateral accretion 
-                           !         caution 1.0, 1.0 best value to be used!!! (gilles G.)  ????
-   fraz_swi    = 0         !  use of frazil ice collection thickness in function of wind (1.0) or not (0.0)
+   hiccrit     = 0.1       !  thickness for new ice formation in open water
+   fraz_swi    = 0         !  use frazil ice collection thickness as a function of wind (1.0) or not (0.0)
    maxfrazb    = 0.0       !  maximum portion of frazil ice collecting at the ice bottom
-   vfrazb      = 0.4166667 !  thresold drift speed for frazil ice collecting at the ice bottom
+   vfrazb      = 0.417     !  thresold drift speed for frazil ice collecting at the ice bottom
    Cfrazb      = 5.0       !  squeezing coefficient for frazil ice collecting at the ice bottom
    hiclim      = 0.10      !  minimum thickness of the first thickness category, must be smaller than hiccrit
    betas       = 0.6       !  exponent in lead-ice fractionation of snow precipitation 0.66
                            !        betas = 1 -> equipartition, betas < 1 -> more on leads
-   kappa_i     = 1.0       !  extinction radiation parameter in sea ice (1.0)
+   kappa_i     = 1.0       !  extinction radiation parameter in sea ice (m-1)
    nconv_i_thd = 50        !  maximal number of iterations for heat diffusion computation
-   maxer_i_thd = 0.0001    !  maximal error in temperature for heat diffusion computation
-   thcon_i_swi = 1         !  switch for computation of thermal conductivity in the ice
-                           !        (0) Untersteiner (1964), (1) Pringle et al. (2007)
+   maxer_i_thd = 0.0001    !  temperature error tolerance after heat diffusion 
+   thcon_i_swi = 1         !  sea ice thermal conductivity, Untersteiner_64 (0) or Pringle_et_al_07 (1)
+   nn_monocat  = 0         !  virtual ITD mono-category parameterizations (1) or not (0); 1 ok for jpl=1 only
+                           !  2=replace ridging by piling, 3=activate G(he) only, 4=activate lateral melting only) --- those are temporary test options
 /
 !-----------------------------------------------------------------------
 &namicesal     !   ice salinity
 !-----------------------------------------------------------------------
-   num_sal     =  2        !  salinity option: 1 -> S = bulk_sal
-                           !                   2 -> S = S(z,t) with a simple parameterization
-                           !                   3 -> S = S(z) profile of Scwharzacher [1959]
-                           !                   4 -> S = S(h) Cox and Weeks [1974]
+   num_sal     =  2        !  S=bulk_sal (1), S(z,t) - (2), S(z) Schwarzacher_59 (3) 
    bulk_sal    =  4.0      !  if 1 is used, it represents the ice salinity
-   sal_G       =  5.00     !  restoring salinity for GD
-   time_G      =  1.728e+6 !  restoring time for GD
-   sal_F       =  2.00     !  restoring salinity for flushing
-   time_F      =  8.640e+5 !  restoring time for flushing
-   s_i_max     = 20.0      !  Maximum salinity 
-   s_i_min     =  0.1      !  Minimum tolerated ice salinity
+   sal_G       =  5.00     !  restoring salinity, gravity drainage (g/kg)
+   time_G      =  1.728e+6 !  restoring time scale, gravity drainage  (s)
+   sal_F       =  2.00     !  restoring salinity, flushing (g/kg)
+   time_F      =  8.640e+5 !  restoring time, flushing (s)
+   s_i_max     = 20.0      !  maximum tolerated ice salinity (g/kg)
+   s_i_min     =  0.1      !  minimum tolerated ice salinity (g/kg)
 /
 !-----------------------------------------------------------------------
 &namiceitdme   !   parameters for mechanical redistribution of ice 
 !-----------------------------------------------------------------------
-   ridge_scheme_swi = 0      !  which ridging scheme using (1=Rothrock,else=Hibler79)
-   Cs               =   0.50 !  shearing energy contribution to ridging
-   Cf               =  17.0  !  ratio of ridging work to PE change in ridging
-   fsnowrdg         =   0.5  !  snow fraction that survives in ridging
-   fsnowrft         =   0.5  !  snow fraction that survives in rafting
-   Gstar            =   0.15 !  fractional area of thin ice being ridged
-   astar            =   0.05 !  equivalent of gstar (0.05 for TH75 and 0.03 for weaker ice)
-   Hstar            = 100.0  !  parameter determining the maximum thickness of ridged ice
-   raft_swi         =   1    !  rafting or not
-   hparmeter        =   0.75 !  threshold thickness for rafting or not
+   ridge_scheme_swi =   0    !  ice strength parameteriztaion, Hibler_79 (0), Rothrock_75 (1)
+   Cs               =   0.50 !  ratio of shearing energy contributing to ridging
+   Cf               =  17.0  !  ratio of ridging work to potential energy change in ridging
+   fsnowrdg         =   0.5  !  snow volume fraction that survives in ridging
+   fsnowrft         =   0.5  !  snow volume fraction that survives in rafting
+   Gstar            =   0.15 !  fractional area of thin ice being ridged (partfun_swi = 0)
+   astar            =   0.05 !  measures fractional area of thin ice being ridged (partfun_swi = 1)
+   Hstar            = 100.0  !  determines the maximum thickness of ridged ice (m)
+   raft_swi         =   1    !  rafting activated (1) or not (0)
+   hparmeter        =   0.75 !  threshold thickness for rafting (m)
    Craft            =   5.0  !  coefficient used in the rafting function
-   ridge_por        =   0.3  !  initial porosity of the ridged ice (typically 0.30)
-   partfun_swi      =   1    !  participation function linear, TH75 (0) or exponential Letal07 (1)
-   brinstren_swi    =   0    !  (1) use brine volume to diminish ice strength
+   ridge_por        =   0.3  !  porosity of newly ridged ice
+   partfun_swi      =   1    !  participation function: linear, Thorndike et al 75 (0) or exponential, Lipscomb 07 (1)
+   brinstren_swi    =   0    !  ice strength is a function of brine volume fraction (1) or not (0)
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5047)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5048)
@@ -165,5 +165,9 @@
    REAL(wp), PUBLIC ::   r1_rdtice        !: = 1. / rdt_ice
 
-   !                                     !!** ice-dynamic namelist (namicedyn) **
+   ! Clem je sais pas ou les mettre
+   INTEGER , PUBLIC ::   nn_hibnd         !: thickness category boundaries: tanh function (1), or h^(-alpha) function (2)
+   REAL(wp), PUBLIC ::   rn_hibnd         !: mean thickness of the domain (used to compute category boundaries, nn_hibnd = 2 only)
+
+   !                                     !!** ice-dynamics namelist (namicedyn) **
    INTEGER , PUBLIC ::   nevp             !: number of iterations for subcycling
    REAL(wp), PUBLIC ::   cw               !: drag coefficient for oceanic stress
@@ -191,4 +195,5 @@
    !                                         ! 3 - salinity profile, constant in time
    INTEGER , PUBLIC ::   thcon_i_swi         !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007)
+   INTEGER , PUBLIC ::   nn_monocat          !: virtual ITD mono-category parameterizations (1) or not (0)
 
    !                                     !!** ice-mechanical redistribution namelist (namiceitdme)
@@ -260,4 +265,8 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res    !: residual component of wfx_ice [kg/m2]
 
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_tot     !: ice concentration tendency (total) [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_thd     !: ice concentration tendency (thermodynamics) [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   afx_dyn     !: ice concentration tendency (dynamics) [s-1]
+
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bom     !: salt flux due to ice growth/melt                      [PSU/m2/s]
@@ -290,4 +299,7 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice   !: transmitted solar radiation under ice
+
+   ! lateral melting (mono-category)
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_i_melt !: net melting (m)
 
    ! temporary arrays for dummy version of the code
@@ -436,5 +448,6 @@
          &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     &
          &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,  qlead  (jpi,jpj) ,     &
-         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl) ,      &
+         &      afx_tot(jpi,jpj) , afx_thd(jpi,jpj),  afx_dyn(jpi,jpj) , &
+         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl) , dh_i_melt(jpi,jpj,jpl) , &
          &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,      &
          &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,   &
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 5047)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 5048)
@@ -172,6 +172,7 @@
       !!              limistate (only) and is changed to 99 m in ice_init
       !!------------------------------------------------------------------
-      INTEGER  ::   jl                   ! dummy loop index
-      REAL(wp) ::   zc1, zc2, zc3, zx1   ! local scalars
+      INTEGER  ::   jl                        ! dummy loop index
+      REAL(wp) ::   zc1, zc2, zc3, zx1        ! local scalars
+      REAL(wp) ::   zhmax, znum, zden, zalpha !
       !!------------------------------------------------------------------
 
@@ -189,14 +190,42 @@
       !- Thickness categories boundaries 
       !----------------------------------
+
+      !  Clem - je sais pas encore dans quelle namelist les mettre, ca depend des chgts liés à bdy
+      nn_hibnd  = 2                !  thickness category boundaries: tanh function (1) h^(-alpha) (2)
+      rn_himean = 2.5              !  mean thickness of the domain (used to compute category boundaries, nn_hibnd = 2 only)
+
       hi_max(:) = 0._wp
 
-      zc1 =  3._wp / REAL( jpl, wp )
-      zc2 = 10._wp * zc1
-      zc3 =  3._wp
-
-      DO jl = 1, jpl
-         zx1 = REAL( jl-1, wp ) / REAL( jpl, wp )
-         hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) )
-      END DO
+      SELECT CASE ( nn_hbnd  )       
+                                   !----------------------
+         CASE (1)                  ! tanh function (CICE)
+                                   !----------------------
+
+         zc1 =  3._wp / REAL( jpl, wp )
+         zc2 = 10._wp * zc1
+         zc3 =  3._wp
+
+         DO jl = 1, jpl
+            zx1 = REAL( jl-1, wp ) / REAL( jpl, wp )
+            hi_max(jl) = hi_max(jl-1) + zc1 + zc2 * (1._wp + TANH( zc3 * (zx1 - 1._wp ) ) )
+         END DO
+
+                                   !----------------------
+         CASE (2)                  ! h^(-alpha) function
+                                   !----------------------
+
+         zalpha = 0.05             ! exponent of the transform function
+
+         zhmax  = 3.*rn_himean
+
+         DO jl = 1, jpl 
+            znum = jpl * ( zhmax+1 )**zalpha
+            zden = ( jpl - jl ) * ( zhmax+1 )**zalpha + jl
+            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
+         END DO
+
+      END SELECT
+
+      ! mv- would be nice if we could put the instruction hi_max(jpl) = 99 here
 
       IF(lwp) WRITE(numout,*) ' Thickness category boundaries '
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 5047)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 5048)
@@ -93,4 +93,7 @@
       CALL lim_thd_lac
       CALL lim_var_glo2eqv    ! only for info
+
+      ! MV: Could put lateral melting here, would be better I think ???
+
      
       IF(ln_ctl) THEN   ! Control print
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
       !
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5047)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5048)
@@ -688,4 +688,15 @@
          END DO
       END DO
+      
+      !
+      !-------------------------------------------------------------
+      ! Net melting (input for thin ice melting, mono-category case)
+      !-------------------------------------------------------------
+      !
+      IF ( ( ( nn_monocat .EQ. 1 ) .OR. ( nn_monocat .EQ. 4 ) ) .AND. ( jpl == 1 ) ) THEN
+         DO ji = kideb, kiut
+            dh_i_melt_1d(ji) = MIN( dh_i_surf(ji) + dh_i_bott(ji) + dh_snowice(ji), 0._wp )
+         END DO
+      ENDIF
 
       CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_rema )
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5047)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5048)
@@ -100,7 +100,9 @@
       INTEGER ::   nconv       ! number of iterations in iterative procedure
       INTEGER ::   minnumeqmin, maxnumeqmax
+      
       INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation
       INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation
       INTEGER, POINTER, DIMENSION(:) ::   isnow      ! switch for presence (1) or absence (0) of snow
+      
       REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system
       REAL(wp) ::   zg1       =  2._wp        !
@@ -113,46 +115,64 @@
       REAL(wp) ::   zerritmax   ! current maximal error on temperature 
       REAL(wp) ::   zhsu
-      REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! ice melting point
-      REAL(wp), POINTER, DIMENSION(:) ::   ztsub       ! old surface temperature (before the iterative procedure )
-      REAL(wp), POINTER, DIMENSION(:) ::   ztsubit     ! surface temperature at previous iteration
-      REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness
-      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness
-      REAL(wp), POINTER, DIMENSION(:) ::   zfsw        ! solar radiation absorbed at the surface
-      REAL(wp), POINTER, DIMENSION(:) ::   zf          ! surface flux function
-      REAL(wp), POINTER, DIMENSION(:) ::   dzf         ! derivative of the surface flux function
-      REAL(wp), POINTER, DIMENSION(:) ::   zerrit      ! current error on temperature
-      REAL(wp), POINTER, DIMENSION(:) ::   zdifcase    ! case of the equation resolution (1->4)
-      REAL(wp), POINTER, DIMENSION(:) ::   zftrice     ! solar radiation transmitted through the ice
-      REAL(wp), POINTER, DIMENSION(:) ::   zihic
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztcond_i    ! Ice thermal conductivity
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_i    ! Radiation transmitted through the ice
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_i    ! Radiation absorbed in the ice
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_i    ! Kappa factor in the ice
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztib        ! Old temperature in the ice
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_i      ! Eta factor in the ice
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztitemp     ! Temporary temperature in the ice to check the convergence
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zspeche_i   ! Ice specific heat
-      REAL(wp), POINTER, DIMENSION(:,:) ::   z_i         ! Vertical cotes of the layers in the ice
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_s    ! Radiation transmited through the snow
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_s    ! Radiation absorbed in the snow
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_s    ! Kappa factor in the snow
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_s      ! Eta factor in the snow
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztstemp     ! Temporary temperature in the snow to check the convergence
-      REAL(wp), POINTER, DIMENSION(:,:) ::   ztsb        ! Temporary temperature in the snow
-      REAL(wp), POINTER, DIMENSION(:,:) ::   z_s         ! Vertical cotes of the layers in the snow
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zswiterm    ! Independent term
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitbis    ! temporary independent term
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiagbis
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid     ! tridiagonal system terms
+      
+      REAL(wp), POINTER, DIMENSION(:)     ::   ztfs        ! ice melting point
+      REAL(wp), POINTER, DIMENSION(:)     ::   ztsub       ! old surface temperature (before the iterative procedure )
+      REAL(wp), POINTER, DIMENSION(:)     ::   ztsubit     ! surface temperature at previous iteration
+      REAL(wp), POINTER, DIMENSION(:)     ::   zh_i        ! ice layer thickness
+      REAL(wp), POINTER, DIMENSION(:)     ::   zh_s        ! snow layer thickness
+      REAL(wp), POINTER, DIMENSION(:)     ::   zfsw        ! solar radiation absorbed at the surface
+      REAL(wp), POINTER, DIMENSION(:)     ::   zf          ! surface flux function
+      REAL(wp), POINTER, DIMENSION(:)     ::   dzf         ! derivative of the surface flux function
+      REAL(wp), POINTER, DIMENSION(:)     ::   zerrit      ! current error on temperature
+      REAL(wp), POINTER, DIMENSION(:)     ::   zdifcase    ! case of the equation resolution (1->4)
+      REAL(wp), POINTER, DIMENSION(:)     ::   zftrice     ! solar radiation transmitted through the ice
+      REAL(wp), POINTER, DIMENSION(:)     ::   zihic
+      
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztcond_i    ! Ice thermal conductivity
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_i    ! Radiation transmitted through the ice
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_i    ! Radiation absorbed in the ice
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_i    ! Kappa factor in the ice
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztib        ! Old temperature in the ice
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_i      ! Eta factor in the ice
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztitemp     ! Temporary temperature in the ice to check the convergence
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zspeche_i   ! Ice specific heat
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_i         ! Vertical cotes of the layers in the ice
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradtr_s    ! Radiation transmited through the snow
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zradab_s    ! Radiation absorbed in the snow
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zkappa_s    ! Kappa factor in the snow
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeta_s      ! Eta factor in the snow
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztstemp     ! Temporary temperature in the snow to check the convergence
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   ztsb        ! Temporary temperature in the snow
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   z_s         ! Vertical cotes of the layers in the snow
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindterm    ! 'Ind'ependent term
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zindtbis    ! Temporary 'ind'ependent term
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zdiagbis    ! Temporary 'dia'gonal term
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid       ! Tridiagonal system terms
+      
       ! diag errors on heat
-      REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx_err
+      REAL(wp), POINTER, DIMENSION(:)     :: zdq, zq_ini, zhfx_err
+      
+      ! Mono-category
+      REAL(wp)                            :: zepsilon      ! determines thres. above which computation of G(h) is done
+      REAL(wp)                            :: zratio_s      ! dummy factor
+      REAL(wp)                            :: zratio_i      ! dummy factor
+      REAL(wp)                            :: zh_thres      ! thickness thres. for G(h) computation
+      REAL(wp)                            :: zhe           ! dummy factor
+      REAL(wp)                            :: zswitch       ! dummy switch
+      REAL(wp)                            :: zkimean       ! mean sea ice thermal conductivity
+      REAL(wp)                            :: zfac          ! dummy factor
+      REAL(wp)                            :: zihe          ! dummy factor
+      REAL(wp)                            :: zheshth       ! dummy factor
+      
+      REAL(wp), POINTER, DIMENSION(:)     :: zghe          ! G(he), th. conduct enhancement factor, mono-cat
+      
       !!------------------------------------------------------------------     
       ! 
       CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow )
       CALL wrk_alloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw )
-      CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic )
+      CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe )
       CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0)
       CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart=0)
-      CALL wrk_alloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis  )
+      CALL wrk_alloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis  )
       CALL wrk_alloc( jpij, nlay_i+3, 3, ztrid )
 
@@ -200,5 +220,5 @@
       !
       !------------------------------------------------------------------------------|
-      ! 2) Radiations                                                                |
+      ! 2) Radiation																 |
       !------------------------------------------------------------------------------|
       !
@@ -213,7 +233,6 @@
       ! zftrice = io.qsr_ice       is below the surface 
       ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice 
+      ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover
       zhsu = 0.1_wp ! threshold for the computation of i0
-         !fr1_i0_1d = i0 for a thin ice surface
-         !fr1_i0_2d = i0 for a thick ice surface
       DO ji = kideb , kiut
          ! switches
@@ -339,43 +358,87 @@
             END DO
          ENDIF
-         !
-         !------------------------------------------------------------------------------|
-         !  5) kappa factors                                                            |
-         !------------------------------------------------------------------------------|
-         !
+         
+         !
+         !------------------------------------------------------------------------------|
+         !  6) G(he) - enhancement of thermal conductivity in mono-category case        |
+         !------------------------------------------------------------------------------|
+         !
+         ! Computation of effective thermal conductivity G(h)
+         ! Used in mono-category case only to simulate an ITD implicitly
+         ! Fichefet and Morales Maqueda, JGR 1997
+
+         zghe(:) = 0._wp
+
+         SELECT CASE ( nn_monocat )
+
+         CASE (0,2,4)
+
+            zghe(kideb:kiut) = 1._wp
+
+         CASE (1,3) ! LIM3
+
+            zepsilon = 0.1
+            zh_thres = EXP( 1._wp ) * zepsilon / 2.
+
+            DO ji = kideb, kiut
+   
+               ! Mean sea ice thermal conductivity
+               zkimean  = SUM( ztcond_i(ji,0:nlay_i) ) / REAL(nlay_i+1,wp)
+
+               ! Effective thickness he (zhe)
+               zfac     = 1._wp / ( rcdsn + zkimean )
+               zratio_s = rcdsn   * zfac
+               zratio_i = zkimean * zfac
+               zhe      = zratio_s * ht_i_1d(ji) + zratio_i * ht_s_1d(ji)
+
+               ! G(he)
+               zswitch  = MAX( 0._wp , SIGN( 1._wp , zhe - zh_thres ) )  ! =0 if zhe < zh_thres, if >
+               zghe(ji) = ( 1.0 - zswitch ) + zswitch * ( 0.5 + 0.5 * LOG( 2.*zhe / zepsilon ) )
+   
+               ! Impose G(he) < 2.
+               zghe(ji) = MIN( zghe(ji), 2.0 )
+
+            END DO
+
+         END SELECT
+
+         !
+         !------------------------------------------------------------------------------|
+         !  7) kappa factors                                                            |
+         !------------------------------------------------------------------------------|
+         !
+         !--- Snow
          DO ji = kideb, kiut
-
-            !-- Snow kappa factors
-            zkappa_s(ji,0)         = rcdsn / MAX(epsi10,zh_s(ji))
-            zkappa_s(ji,nlay_s)    = rcdsn / MAX(epsi10,zh_s(ji))
+            zfac                  =  1. / MAX( epsi10 , zh_s(ji) )
+            zkappa_s(ji,0)        = zghe(ji) * rcdsn * zfac
+            zkappa_s(ji,nlay_s)   = zghe(ji) * rcdsn * zfac
          END DO
 
          DO jk = 1, nlay_s-1
             DO ji = kideb , kiut
-               zkappa_s(ji,jk)  = 2.0 * rcdsn / &
-                  MAX(epsi10,2.0*zh_s(ji))
-            END DO
-         END DO
-
+               zkappa_s(ji,jk)    = zghe(ji) * 2.0 * rcdsn / MAX( epsi10, 2.0*zh_s(ji) )
+            END DO
+         END DO
+
+         !--- Ice
          DO jk = 1, nlay_i-1
             DO ji = kideb , kiut
-               !-- Ice kappa factors
-               zkappa_i(ji,jk)  = 2.0*ztcond_i(ji,jk)/ &
-                  MAX(epsi10,2.0*zh_i(ji)) 
-            END DO
-         END DO
-
-         DO ji = kideb , kiut
-            zkappa_i(ji,0)        = ztcond_i(ji,0)/MAX(epsi10,zh_i(ji))
-            zkappa_i(ji,nlay_i)   = ztcond_i(ji,nlay_i) / MAX(epsi10,zh_i(ji))
-            !-- Interface
-            zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(epsi10, &
-               (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji)))
-            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) &
-               + zkappa_i(ji,0)*REAL( 1 - isnow(ji) )
-         END DO
-         !
-         !------------------------------------------------------------------------------|
-         ! 6) Sea ice specific heat, eta factors                                        |
+               zkappa_i(ji,jk)    = zghe(ji) * 2.0 * ztcond_i(ji,jk) / MAX( epsi10 , 2.0*zh_i(ji) )
+            END DO
+         END DO
+
+         !--- Snow-ice interface
+         DO ji = kideb , kiut
+            zfac                  = 1./ MAX( epsi10 , zh_i(ji) )
+            zkappa_i(ji,0)        = zghe(ji) * ztcond_i(ji,0) * zfac
+            zkappa_i(ji,nlay_i)   = zghe(ji) * ztcond_i(ji,nlay_i) * zfac
+            zkappa_s(ji,nlay_s)   = zghe(ji) * zghe(ji) * 2.0 * rcdsn * ztcond_i(ji,0) / & 
+           &                        MAX( epsi10, ( zghe(ji) * ztcond_i(ji,0) * zh_s(ji) + zghe(ji) * rcdsn * zh_i(ji) ) )
+            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji), wp ) + zkappa_i(ji,0)*REAL( 1 - isnow(ji), wp )
+         END DO
+
+         !
+         !------------------------------------------------------------------------------|
+         ! 8) Sea ice specific heat, eta factors                                        |
          !------------------------------------------------------------------------------|
          !
@@ -383,8 +446,6 @@
             DO ji = kideb , kiut
                ztitemp(ji,jk)   = t_i_1d(ji,jk)
-               zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ &
-                  MAX((t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt),epsi10)
-               zeta_i(ji,jk)    = rdt_ice / MAX(rhoic*zspeche_i(ji,jk)*zh_i(ji), &
-                  epsi10)
+               zspeche_i(ji,jk) = cpic + zgamma*s_i_1d(ji,jk)/ MAX( (t_i_1d(ji,jk)-rtt)*(ztib(ji,jk)-rtt) , epsi10 )
+               zeta_i(ji,jk)    = rdt_ice / MAX( rhoic*zspeche_i(ji,jk)*zh_i(ji), epsi10 )
             END DO
          END DO
@@ -396,7 +457,8 @@
             END DO
          END DO
-         !
-         !------------------------------------------------------------------------------|
-         ! 7) surface flux computation                                                  |
+
+         !
+         !------------------------------------------------------------------------------|
+         ! 9) surface flux computation                                                  |
          !------------------------------------------------------------------------------|
          !
@@ -418,5 +480,5 @@
          !
          !------------------------------------------------------------------------------|
-         ! 8) tridiagonal system terms                                                  |
+         ! 10) tridiagonal system terms                                                 |
          !------------------------------------------------------------------------------|
          !
@@ -433,6 +495,6 @@
                ztrid(ji,numeq,2) = 0.
                ztrid(ji,numeq,3) = 0.
-               zswiterm(ji,numeq)= 0.
-               zswitbis(ji,numeq)= 0.
+               zindterm(ji,numeq)= 0.
+               zindtbis(ji,numeq)= 0.
                zdiagbis(ji,numeq)= 0.
             ENDDO
@@ -446,5 +508,5 @@
                   zkappa_i(ji,jk))
                ztrid(ji,numeq,3)  =  - zeta_i(ji,jk)*zkappa_i(ji,jk)
-               zswiterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk)* &
+               zindterm(ji,numeq) =  ztib(ji,jk) + zeta_i(ji,jk)* &
                   zradab_i(ji,jk)
             END DO
@@ -458,5 +520,5 @@
                +  zkappa_i(ji,nlay_i-1) )
             ztrid(ji,numeq,3)  =  0.0
-            zswiterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* &
+            zindterm(ji,numeq) =  ztib(ji,nlay_i) + zeta_i(ji,nlay_i)* &
                ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i)*zg1 &
                *  t_bo_1d(ji) ) 
@@ -478,5 +540,5 @@
                      zkappa_s(ji,jk) )
                   ztrid(ji,numeq,3)   =  - zeta_s(ji,jk)*zkappa_s(ji,jk)
-                  zswiterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk)* &
+                  zindterm(ji,numeq)  =  ztsb(ji,jk) + zeta_s(ji,jk)* &
                      zradab_s(ji,jk)
                END DO
@@ -485,5 +547,5 @@
                IF ( nlay_i.eq.1 ) THEN
                   ztrid(ji,nlay_s+2,3)    =  0.0
-                  zswiterm(ji,nlay_s+2)   =  zswiterm(ji,nlay_s+2) + zkappa_i(ji,1)* &
+                  zindterm(ji,nlay_s+2)   =  zindterm(ji,nlay_s+2) + zkappa_i(ji,1)* &
                      t_bo_1d(ji) 
                ENDIF
@@ -502,5 +564,5 @@
                   ztrid(ji,1,2) = dzf(ji) - zg1s*zkappa_s(ji,0)
                   ztrid(ji,1,3) = zg1s*zkappa_s(ji,0)
-                  zswiterm(ji,1) = dzf(ji)*t_su_1d(ji)   - zf(ji)
+                  zindterm(ji,1) = dzf(ji)*t_su_1d(ji)   - zf(ji)
 
                   !!first layer of snow equation
@@ -508,5 +570,5 @@
                   ztrid(ji,2,2)  =  1.0 + zeta_s(ji,1)*(zkappa_s(ji,1) + zkappa_s(ji,0)*zg1s)
                   ztrid(ji,2,3)  =  - zeta_s(ji,1)* zkappa_s(ji,1)
-                  zswiterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
+                  zindterm(ji,2) =  ztsb(ji,1) + zeta_s(ji,1)*zradab_s(ji,1)
 
                ELSE 
@@ -525,5 +587,5 @@
                      zkappa_s(ji,0) * zg1s )
                   ztrid(ji,2,3)  =  - zeta_s(ji,1)*zkappa_s(ji,1) 
-                  zswiterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *            &
+                  zindterm(ji,2) = ztsb(ji,1) + zeta_s(ji,1) *            &
                      ( zradab_s(ji,1) +                         &
                      zkappa_s(ji,0) * zg1s * t_su_1d(ji) ) 
@@ -549,5 +611,5 @@
                   ztrid(ji,numeqmin(ji),2)   =  dzf(ji) - zkappa_i(ji,0)*zg1    
                   ztrid(ji,numeqmin(ji),3)   =  zkappa_i(ji,0)*zg1
-                  zswiterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji)
+                  zindterm(ji,numeqmin(ji))  =  dzf(ji)*t_su_1d(ji) - zf(ji)
 
                   !!first layer of ice equation
@@ -556,5 +618,5 @@
                      + zkappa_i(ji,0) * zg1 )
                   ztrid(ji,numeqmin(ji)+1,3) =  - zeta_i(ji,1)*zkappa_i(ji,1)  
-                  zswiterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)  
+                  zindterm(ji,numeqmin(ji)+1)=  ztib(ji,1) + zeta_i(ji,1)*zradab_i(ji,1)  
 
                   !!case of only one layer in the ice (surface & ice equations are altered)
@@ -569,5 +631,5 @@
                      ztrid(ji,numeqmin(ji)+1,3)  =  0.0
 
-                     zswiterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1)* &
+                     zindterm(ji,numeqmin(ji)+1) =  ztib(ji,1) + zeta_i(ji,1)* &
                         ( zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji) )
                   ENDIF
@@ -589,5 +651,5 @@
                      zg1)  
                   ztrid(ji,numeqmin(ji),3)      =  - zeta_i(ji,1) * zkappa_i(ji,1)
-                  zswiterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
+                  zindterm(ji,numeqmin(ji))     =  ztib(ji,1) + zeta_i(ji,1)*( zradab_i(ji,1) + &
                      zkappa_i(ji,0) * zg1 * t_su_1d(ji) ) 
 
@@ -598,5 +660,5 @@
                         zkappa_i(ji,1))
                      ztrid(ji,numeqmin(ji),3)  =  0.0
-                     zswiterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1)* &
+                     zindterm(ji,numeqmin(ji)) =  ztib(ji,1) + zeta_i(ji,1)* &
                         (zradab_i(ji,1) + zkappa_i(ji,1)*t_bo_1d(ji)) &
                         + t_su_1d(ji)*zeta_i(ji,1)*zkappa_i(ji,0)*2.0
@@ -610,5 +672,5 @@
          !
          !------------------------------------------------------------------------------|
-         ! 9) tridiagonal system solving                                                |
+         ! 11) tridiagonal system solving                                               |
          !------------------------------------------------------------------------------|
          !
@@ -622,5 +684,5 @@
 
          DO ji = kideb , kiut
-            zswitbis(ji,numeqmin(ji)) =  zswiterm(ji,numeqmin(ji))
+            zindtbis(ji,numeqmin(ji)) =  zindterm(ji,numeqmin(ji))
             zdiagbis(ji,numeqmin(ji)) =  ztrid(ji,numeqmin(ji),2)
             minnumeqmin               =  MIN(numeqmin(ji),minnumeqmin)
@@ -633,6 +695,6 @@
                zdiagbis(ji,numeq)  =  ztrid(ji,numeq,2) - ztrid(ji,numeq,1)* &
                   ztrid(ji,numeq-1,3)/zdiagbis(ji,numeq-1)
-               zswitbis(ji,numeq)  =  zswiterm(ji,numeq) - ztrid(ji,numeq,1)* &
-                  zswitbis(ji,numeq-1)/zdiagbis(ji,numeq-1)
+               zindtbis(ji,numeq)  =  zindterm(ji,numeq) - ztrid(ji,numeq,1)* &
+                  zindtbis(ji,numeq-1)/zdiagbis(ji,numeq-1)
             END DO
          END DO
@@ -640,5 +702,5 @@
          DO ji = kideb , kiut
             ! ice temperatures
-            t_i_1d(ji,nlay_i)    =  zswitbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
+            t_i_1d(ji,nlay_i)    =  zindtbis(ji,numeqmax(ji))/zdiagbis(ji,numeqmax(ji))
          END DO
 
@@ -646,5 +708,5 @@
             DO ji = kideb , kiut
                jk    =  numeq - nlay_s - 1
-               t_i_1d(ji,jk)  =  (zswitbis(ji,numeq) - ztrid(ji,numeq,3)* &
+               t_i_1d(ji,jk)  =  (zindtbis(ji,numeq) - ztrid(ji,numeq,3)* &
                   t_i_1d(ji,jk+1))/zdiagbis(ji,numeq)
             END DO
@@ -654,5 +716,5 @@
             ! snow temperatures      
             IF (ht_s_1d(ji).GT.0._wp) &
-               t_s_1d(ji,nlay_s)     =  (zswitbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
+               t_s_1d(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
                *  t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) &
                *        MAX(0.0,SIGN(1.0,ht_s_1d(ji))) 
@@ -662,10 +724,10 @@
             ztsubit(ji) = t_su_1d(ji)
             IF( t_su_1d(ji) < ztfs(ji) ) &
-               t_su_1d(ji) = ( zswitbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1)   &
+               t_su_1d(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_1d(ji,1)   &
                &          + REAL( 1 - isnow(ji) )*t_i_1d(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))  
          END DO
          !
          !--------------------------------------------------------------------------
-         !  10) Has the scheme converged ?, end of the iterative procedure         |
+         !  12) Has the scheme converged ?, end of the iterative procedure         |
          !--------------------------------------------------------------------------
          !
@@ -709,5 +771,5 @@
       !
       !-------------------------------------------------------------------------!
-      !   11) Fluxes at the interfaces                                          !
+      !   13) Fluxes at the interfaces                                          !
       !-------------------------------------------------------------------------!
       DO ji = kideb, kiut
@@ -771,9 +833,9 @@
       CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow )
       CALL wrk_dealloc( jpij, ztfs, ztsub, ztsubit, zh_i, zh_s, zfsw )
-      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic )
+      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zghe )
       CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i,   &
          &              ztib, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 )
       CALL wrk_dealloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsb, zeta_s, ztstemp, z_s, kjstart = 0 )
-      CALL wrk_dealloc( jpij, nlay_i+3, zswiterm, zswitbis, zdiagbis )
+      CALL wrk_dealloc( jpij, nlay_i+3, zindterm, zindtbis, zdiagbis )
       CALL wrk_dealloc( jpij, nlay_i+3, 3, ztrid )
       CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx_err )
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5047)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5048)
@@ -420,4 +420,16 @@
          END DO
          ! -------------------------------------------------
+         
+         !--------------------------------------
+         ! Impose a_i < amax in mono-category
+         !--------------------------------------
+         !
+         IF ( ( nn_monocat == 2 ) .AND. ( jpl == 1 ) ) THEN ! simple conservative piling, comparable with LIM2
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  a_i(ji,jj,1)  = MIN( a_i(ji,jj,1), amax )
+               END DO
+            END DO
+         ENDIF
 
          ! --- diags ---
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 5047)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 5048)
@@ -200,15 +200,19 @@
 
       ztmp = rday / rhoic
-      CALL iom_put( "vfxres"     , wfx_res * ztmp  )             ! daily prod./melting due to limupdate 
-      CALL iom_put( "vfxopw"     , wfx_opw * ztmp  )             ! daily lateral thermodynamic ice production
-      CALL iom_put( "vfxsni"     , wfx_sni * ztmp  )             ! daily snowice ice production
-      CALL iom_put( "vfxbog"     , wfx_bog * ztmp  )             ! daily bottom thermodynamic ice production
-      CALL iom_put( "vfxdyn"     , wfx_dyn * ztmp  )             ! daily dynamic ice production (rid/raft)
-      CALL iom_put( "vfxsum"     , wfx_sum * ztmp  )             ! surface melt 
-      CALL iom_put( "vfxbom"     , wfx_bom * ztmp  )             ! bottom melt 
-      CALL iom_put( "vfxice"     , wfx_ice * ztmp  )             ! total ice growth/melt 
-      CALL iom_put( "vfxsnw"     , wfx_snw * ztmp  )             ! total snw growth/melt 
-      CALL iom_put( "vfxsub"     , wfx_sub * ztmp  )             ! sublimation (snow) 
-      CALL iom_put( "vfxspr"     , wfx_spr * ztmp  )             ! precip (snow) 
+      CALL iom_put( "vfxres"     , wfx_res * ztmp       )        ! daily prod./melting due to limupdate 
+      CALL iom_put( "vfxopw"     , wfx_opw * ztmp       )        ! daily lateral thermodynamic ice production
+      CALL iom_put( "vfxsni"     , wfx_sni * ztmp       )        ! daily snowice ice production
+      CALL iom_put( "vfxbog"     , wfx_bog * ztmp       )        ! daily bottom thermodynamic ice production
+      CALL iom_put( "vfxdyn"     , wfx_dyn * ztmp       )        ! daily dynamic ice production (rid/raft)
+      CALL iom_put( "vfxsum"     , wfx_sum * ztmp       )        ! surface melt 
+      CALL iom_put( "vfxbom"     , wfx_bom * ztmp       )        ! bottom melt 
+      CALL iom_put( "vfxice"     , wfx_ice * ztmp       )        ! total ice growth/melt 
+      CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt 
+      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow) 
+      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow)
+      
+      CALL iom_put( "afxtot"     , afx_tot              )        ! concentration tendency (total)
+      CALL iom_put( "afxdyn"     , afx_dyn              )        ! concentration tendency (dynamics)
+      CALL iom_put( "afxthd"     , afx_thd              )        ! concentration tendency (thermo)
 
       CALL iom_put ('hfxthd', hfx_thd(:,:) )   !  
@@ -216,5 +220,5 @@
       CALL iom_put ('hfxres', hfx_res(:,:) )   !  
       CALL iom_put ('hfxout', hfx_out(:,:) )   !  
-      CALL iom_put ('hfxin' , hfx_in(:,:) )   !  
+      CALL iom_put ('hfxin' , hfx_in(:,:) )    !  
       CALL iom_put ('hfxsnw', hfx_snw(:,:) )   !  
       CALL iom_put ('hfxsub', hfx_sub(:,:) )   !  
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90	(revision 5047)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/par_ice.F90	(revision 5048)
@@ -12,5 +12,5 @@
 
    !                                             !!! ice thermodynamics
-   INTEGER, PUBLIC, PARAMETER ::   nlay_i   = 5   !: number of ice layers
+   INTEGER, PUBLIC, PARAMETER ::   nlay_i   = 2   !: number of ice layers
    INTEGER, PUBLIC, PARAMETER ::   nlay_s   = 1   !: number of snow layers
 
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 5047)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 5048)
@@ -35,6 +35,7 @@
    !: are the variables corresponding to 2d vectors
 
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npb    !: number of points where computations has to be done
-   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npac   !: correspondance between points (lateral accretion)
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npb    !: address vector for 1d vertical thermo computations
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   nplm   !: address vector for mono-category lateral melting
+   INTEGER , PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   npac   !: address vector for new ice formation
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qlead_1d     
@@ -109,4 +110,5 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_bott     !: Ice bottom accretion/ablation  [m]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_snowice    !: Snow ice formation             [m of ice]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dh_i_melt_1d  !: Net melting [m], for mono-cat lateral melting
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sm_i_1d       !: Ice bulk salinity [ppt]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   s_i_new       !: Salinity of new ice at the bottom
@@ -138,5 +140,5 @@
       !!---------------------------------------------------------------------!
 
-      ALLOCATE( npb      (jpij) , npac     (jpij),                          &
+      ALLOCATE( npb      (jpij) , nplm        (jpij) , npac     (jpij),   &
          !                                                                  !
          &      qlead_1d (jpij) , ftr_ice_1d  (jpij) ,     &
@@ -165,5 +167,6 @@
          &      ht_s_1d    (jpij) , fc_su    (jpij) , fc_bo_i  (jpij) ,    &    
          &      dh_s_tot  (jpij) , dh_i_surf(jpij) , dh_i_bott(jpij) ,    &    
-         &      dh_snowice(jpij) , sm_i_1d   (jpij) , s_i_new  (jpij) ,    &
+         &      dh_snowice(jpij) , dh_i_melt_1d(jpij) , &
+         &      sm_i_1d   (jpij) , s_i_new  (jpij) ,    &
          &      t_s_1d(jpij,nlay_s),                                       &
          &      t_i_1d(jpij,nlay_i+1), s_i_1d(jpij,nlay_i+1)                ,     &            
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5047)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5048)
@@ -191,15 +191,4 @@
             CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )
 
-            ! MV -> seb 
-!           CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=t_su    )
-
-!           IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, zalb_ice, qns_ice, qsr_ice ,   &
-!              &                                           dqns_ice, qla_ice, dqla_ice, nn_limflx )
-!           ! Latent heat flux is forced to 0 in coupled :
-!           !  it is included in qns (non-solar heat flux)
-!           qla_ice  (:,:,:) = 0._wp
-!           dqla_ice (:,:,:) = 0._wp
-            ! END MV -> seb
-            !
          END SELECT
          
@@ -236,4 +225,7 @@
          wfx_spr(:,:) = 0._wp   ;   
 
+         afx_tot(:,:) = at_i(:,:) ;  afx_dyn(:,:) = 0._wp
+         afx_thd(:,:) = 0._wp
+
          hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp
          hfx_thd(:,:) = 0._wp   ;   
@@ -252,9 +244,12 @@
          ! ----------------------------------------------
          IF( .NOT. lk_c1d ) THEN
+
+         IF ( ln_limdyn ) afx_dyn(:,:)     = at_i(:,:)
+
                           CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics )
                           CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion )
                           CALL lim_var_glo2eqv            ! equivalent variables, requested for rafting
          IF( ln_nicep )   CALL lim_prt_state( kt, jiindx, jjindx,-1, ' - ice dyn & trp - ' )   ! control print
-                          CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting)
+         IF( nn_monocat /= 2 ) CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting)
                           CALL lim_var_agg( 1 ) 
 #if defined key_bdy
@@ -278,7 +273,10 @@
                           oa_i_b (:,:,:)   = oa_i (:,:,:)
                           smv_i_b(:,:,:)   = smv_i(:,:,:)
+
+         IF ( ln_limdyn ) afx_dyn(:,:)     = ( at_i(:,:) - afx_dyn(:,:) ) * r1_rdtice
+                          afx_thd(:,:)     = at_i(:,:)
  
          ! ----------------------------------------------
-         ! ice thermodynamic
+         ! ice thermodynamics
          ! ----------------------------------------------
                           CALL lim_var_glo2eqv            ! equivalent variables
@@ -288,5 +286,4 @@
                           phicif(:,:)  = vt_i(:,:)
 
-                          ! MV -> seb
                           SELECT CASE( kblk )
                              CASE ( jp_cpl )
@@ -299,5 +296,4 @@
                              dqla_ice (:,:,:) = 0._wp
                           END SELECT
-                          ! END MV -> seb
                           !
                           CALL lim_var_bv                 ! bulk brine volume (diag)
@@ -320,4 +316,7 @@
          !                                           ! Diagnostics and outputs 
          IF (ln_limdiaout) CALL lim_diahsb
+
+                          afx_thd(:,:)     = ( at_i(:,:) - afx_thd(:,:) ) * r1_rdtice
+                          afx_tot(:,:)     = ( at_i(:,:) - afx_tot(:,:) ) * r1_rdtice
 
                           CALL lim_wri( 1  )              ! Ice outputs 
