Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 8313)
@@ -17,4 +17,8 @@
    nlay_i           =    2          !  number of ice  layers
    nlay_s           =    1          !  number of snow layers (only 1 is working)
+   nn_monocat       =    0          !  virtual ITD mono-category parameterizations (1-4 => jpl = 1 only) or not (0)
+                                    !     2: simple piling instead of ridging    --- temporary option
+                                    !     3: activate G(he) only                 --- temporary option
+                                    !     4: activate extra lateral melting only --- temporary option
    rn_amax_n        =   0.997       !  maximum tolerated ice concentration NH
    rn_amax_s        =   0.997       !  maximum tolerated ice concentration SH
@@ -70,8 +74,5 @@
 &namiceitd     !   Ice discretization
 !------------------------------------------------------------------------------
-   nn_catbnd      =    2           !  computation of ice category boundaries based on
-                                   !      1: tanh function
-                                   !      2: h^(-alpha), function of rn_himean
-   rn_himean      =    2.0         !  expected domain-average ice thickness (m), nn_catbnd = 2 only
+   rn_himean      =    2.0         !  expected domain-average ice thickness (m)
 /
 !------------------------------------------------------------------------------
@@ -91,4 +92,5 @@
                                     !
             ! -- limdyn & limrhg -- !
+   rn_ishlat      =    2.           !  free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2)
    rn_cio         =    5.0e-03      !  ice-ocean drag coefficient (-)
    rn_creepl      =    1.0e-12      !  creep limit (s-1)
@@ -110,16 +112,10 @@
                  ! -- limthd_dif -- !
    rn_kappa_i     = 1.0             !  radiation attenuation coefficient in sea ice (m-1)
-   nn_conv_dif    = 50              !  maximal number of iterations for heat diffusion computation
-   rn_terr_dif    = 1.0e-04         !  maximum temperature after heat diffusion (degC)
    nn_ice_thcon   = 1               !  sea ice thermal conductivity
                                     !     0: k = k0 + beta.S/T            (Untersteiner, 1964)
                                     !     1: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007)
-   ln_it_qnsice   = .true.          !  iterate the surface non-solar flux with surface temperature (T) or not (F)
-   nn_monocat     = 0               !  virtual ITD mono-category parameterizations (1, jpl = 1 only) or not (0)
-                                    !     2: simple piling instead of ridging    --- temporary option
-                                    !     3: activate G(he) only                 --- temporary option
-                                    !     4: activate extra lateral melting only --- temporary option
-   rn_cdsn     = 0.31              !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971)
-                                   !  Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
+   ln_dqnsice     = .true.          !  change the surface non-solar flux with surface temperature (T) or not (F)
+   rn_cdsn        = 0.31            !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971)
+                                    !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013)
                   ! -- limthd_dh -- !
    ln_limdH       = .true.          !  activate ice thickness change from growing/melting (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 8313)
@@ -201,7 +201,6 @@
                            !  =1 use observed ice-cover      ,
                            !  =3 to 4 :  ice-model used (LIM3 or CICE)                         ("key_lim3" or "key_cice")
-   nn_ice_embd = 1         !  =0 levitating ice (no mass exchange, concentration/dilution effect)
-                           !  =1 levitating ice with mass and salt exchange but no presure effect
-                           !  =2 embedded sea-ice (full salt and mass exchanges and pressure)
+   ln_ice_embd = .false.   !  =F levitating ice with mass and salt exchange but no presure effect
+                           !  =T embedded sea-ice (full salt and mass exchanges and pressure)
                      ! Misc. options of sbc : 
    ln_traqsr   = .true.    !  Light penetration in the ocean            (T => fill namtra_qsr)
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8313)
@@ -196,5 +196,4 @@
    
    !                                     !!** ice-thickness distribution namelist (namiceitd) **
-   INTEGER , PUBLIC ::   nn_catbnd        !: categories distribution following: tanh function (1), or h^(-alpha) function (2)
    REAL(wp), PUBLIC ::   rn_himean        !: mean thickness of the domain (used to compute the distribution, nn_itdshp = 2 only)
 
@@ -210,4 +209,5 @@
    LOGICAL , PUBLIC ::   ln_icestr_bvf    !: use brine volume to diminish ice strength
                                           ! -- limdyn & limrhg -- !
+   REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice
    REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress
    REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9
@@ -223,9 +223,7 @@
                                           ! -- limthd_dif -- !
    REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m]
-   REAL(wp), PUBLIC ::   nn_conv_dif      !: maximal number of iterations for heat diffusion
-   REAL(wp), PUBLIC ::   rn_terr_dif      !: maximal tolerated error (C) for heat diffusion
    INTEGER , PUBLIC ::   nn_ice_thcon     !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007)
-   LOGICAL , PUBLIC ::   ln_it_qnsice     !: iterate surface flux with changing surface temperature or not (F)
-   INTEGER , PUBLIC ::   nn_monocat       !: virtual ITD mono-category parameterizations (1) or not (0)
+   LOGICAL , PUBLIC ::   ln_dqnsice       !: change non-solar surface flux with changing surface temperature (T) or not (F)
+   INTEGER , PUBLIC ::   nn_monocat       !: virtual ITD mono-category parameterizations (1-4) or not (0)
    REAL(wp), PUBLIC ::   rn_cdsn          !: thermal conductivity of the snow [W/m/K]
                                           ! -- limthd_dh -- !
@@ -310,7 +308,4 @@
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]     
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frld        !: Leads fraction = 1 - ice fraction
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   pfrld       !: Leads fraction at previous time  
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   phicif      !: Old ice thickness
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qlead       !: heat balance of the lead (or of the open ocean)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhtur       !: net downward heat flux from the ice to the ocean
@@ -534,5 +529,5 @@
 
       ii = ii + 1
-      ALLOCATE( t_bo   (jpi,jpj) , frld       (jpi,jpj) , pfrld      (jpi,jpj) , phicif     (jpi,jpj) ,  &
+      ALLOCATE( t_bo   (jpi,jpj) ,                                                                       &
          &      wfx_snw(jpi,jpj) , wfx_snw_dyn(jpi,jpj) , wfx_snw_sum(jpi,jpj) , wfx_snw_sub(jpi,jpj) ,  &
          &      wfx_ice(jpi,jpj) , wfx_sub    (jpi,jpj) , wfx_ice_sub(jpi,jpj) , wfx_lam    (jpi,jpj) ,  &
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90	(revision 8313)
@@ -385,6 +385,6 @@
                WRITE(numout,*) '   ~~~~~~~~~~~~~~~~ '
                WRITE(numout,*) ' - Heat fluxes in and out the ice ***'
-               WRITE(numout,*) ' qsr_ini       : ', pfrld(ji,jj) * qsr(ji,jj) + SUM( a_i_b(ji,jj,:) * qsr_ice(ji,jj,:) )
-               WRITE(numout,*) ' qns_ini       : ', pfrld(ji,jj) * qns(ji,jj) + SUM( a_i_b(ji,jj,:) * qns_ice(ji,jj,:) )
+               WRITE(numout,*) ' qsr_ini       : ', (1._wp-at_i_b(ji,jj)) * qsr(ji,jj) + SUM( a_i_b(ji,jj,:) * qsr_ice(ji,jj,:) )
+               WRITE(numout,*) ' qns_ini       : ', (1._wp-at_i_b(ji,jj)) * qns(ji,jj) + SUM( a_i_b(ji,jj,:) * qns_ice(ji,jj,:) )
                WRITE(numout,*)
                WRITE(numout,*) 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 8313)
@@ -109,6 +109,7 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicedyn/ nn_limadv, nn_limadv_ord,  &
-         &                nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, &
+      NAMELIST/namicedyn/ nn_limadv, nn_limadv_ord,                                &
+         &                nn_icestr, rn_pe_rdg, rn_pstar, rn_crhg, ln_icestr_bvf,  &
+         &                rn_ishlat, rn_cio, rn_creepl, rn_ecc,                    &
          &                nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr, rn_lfrelax
       !!-------------------------------------------------------------------
@@ -130,11 +131,13 @@
          WRITE(numout,*)'    choose the advection scheme (-1=Prather, 0=Ulimate-Macho)   nn_limadv     = ', nn_limadv 
          WRITE(numout,*)'    choose the order of the scheme (if ultimate)                nn_limadv_ord = ', nn_limadv_ord  
-         ! limrhg
+         ! limitd_me
          WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)         nn_icestr     = ', nn_icestr 
-         WRITE(numout,*)'    Including brine volume in ice strength comp.                ln_icestr_bvf = ', ln_icestr_bvf
          WRITE(numout,*)'    Ratio of ridging work to PotEner change in ridging          rn_pe_rdg     = ', rn_pe_rdg 
-         WRITE(numout,*) '   drag coefficient for oceanic stress                         rn_cio        = ', rn_cio
          WRITE(numout,*) '   first bulk-rheology parameter                               rn_pstar      = ', rn_pstar
          WRITE(numout,*) '   second bulk-rhelogy parameter                               rn_crhg       = ', rn_crhg
+         WRITE(numout,*)'    Including brine volume in ice strength comp.                ln_icestr_bvf = ', ln_icestr_bvf
+         ! limrhg
+         WRITE(numout,*) '   lateral boundary condition for sea ice dynamics             rn_ishlat     = ', rn_ishlat
+         WRITE(numout,*) '   drag coefficient for oceanic stress                         rn_cio        = ', rn_cio
          WRITE(numout,*) '   creep limit                                                 rn_creepl     = ', rn_creepl
          WRITE(numout,*) '   eccentricity of the elliptical yield curve                  rn_ecc        = ', rn_ecc
@@ -142,7 +145,13 @@
          WRITE(numout,*) '   ratio of elastic timescale over ice time step               rn_relast     = ', rn_relast
          WRITE(numout,*) '   Landfast: param (T or F)                                    ln_landfast   = ', ln_landfast
-         WRITE(numout,*) '   Landfast: fraction of ocean depth that ice must reach       rn_gamma      = ', rn_gamma
-         WRITE(numout,*) '   Landfast: maximum bottom stress per unit area of contact    rn_icebfr     = ', rn_icebfr
-         WRITE(numout,*) '   Landfast: relax time scale (s-1) to reach static friction   rn_lfrelax    = ', rn_lfrelax
+         WRITE(numout,*) '      T: fraction of ocean depth that ice must reach           rn_gamma      = ', rn_gamma
+         WRITE(numout,*) '      T: maximum bottom stress per unit area of contact        rn_icebfr     = ', rn_icebfr
+         WRITE(numout,*) '      T: relax time scale (s-1) to reach static friction       rn_lfrelax    = ', rn_lfrelax
+      ENDIF
+      !
+      IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  free-slip '
+      ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  no-slip '
+      ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  partial-slip '
+      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  strong-slip '
       ENDIF
       !
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 8313)
@@ -150,5 +150,4 @@
       REAL(wp), PARAMETER               ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
       REAL(wp), PARAMETER               ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
-      REAL(wp), PARAMETER               ::   zshlat = 2._wp                  ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip)
       !!-------------------------------------------------------------------
 
@@ -182,5 +181,5 @@
          DO ji = fs_2, fs_jpim1   ! vector opt.
             IF( zfmask(ji,jj) == 0._wp ) THEN
-               zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
+               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
             ENDIF
          END DO
@@ -188,16 +187,16 @@
       DO jj = 2, jpjm1
          IF( zfmask(1,jj) == 0._wp ) THEN
-            zfmask(1  ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
+            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
          ENDIF
          IF( zfmask(jpi,jj) == 0._wp ) THEN
-            zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
+            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
          ENDIF
       END DO
       DO ji = 2, jpim1
          IF( zfmask(ji,1) == 0._wp ) THEN
-            zfmask(ji,1  ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
+            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
          ENDIF
          IF( zfmask(ji,jpj) == 0._wp ) THEN
-            zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
+            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
          ENDIF
       END DO
@@ -251,5 +250,5 @@
       !------------------------------------------------------------------------------!
 
-      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
+      IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
          !                                            
          ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 8313)
@@ -112,6 +112,6 @@
       ! --- case we bypass ice thermodynamics --- !
       IF( .NOT. ln_limthd ) THEN   ! we suppose ice is impermeable => ocean is isolated from atmosphere
-         hfx_in   (:,:)   = pfrld(:,:) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
-         hfx_out  (:,:)   = pfrld(:,:) *   qns_oce(:,:)                  + qemp_oce(:,:)
+         hfx_in   (:,:)   = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
+         hfx_out  (:,:)   = ( 1._wp - at_i_b(:,:) ) *   qns_oce(:,:)                  + qemp_oce(:,:)
          ftr_ice  (:,:,:) = 0._wp
          emp_ice  (:,:)   = 0._wp
@@ -195,15 +195,13 @@
          &     + sfx_res(:,:) + sfx_dyn(:,:) + sfx_bri(:,:) + sfx_sub(:,:) + sfx_lam(:,:)
 
-      !-------------------------------------------------------------!
-      !   mass of snow and ice per unit area for embedded sea-ice   !
-      !-------------------------------------------------------------!
-      IF( nn_ice_embd /= 0 ) THEN
-         ! save mass from the previous ice time step
-         snwice_mass_b(:,:) = snwice_mass(:,:)                  
-         ! new mass per unit area
-         snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  ) 
-         ! time evolution of snow+ice mass
-         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice
-      ENDIF
+      !----------------------------------------!
+      !   mass of snow and ice per unit area   !
+      !----------------------------------------!
+      ! save mass from the previous ice time step
+      snwice_mass_b(:,:) = snwice_mass(:,:)                  
+      ! new mass per unit area
+      snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  ) 
+      ! time evolution of snow+ice mass
+      snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice
 
       !-----------------------------------------------!
@@ -243,5 +241,5 @@
       !!                      tmod_io = rhoco * | U_ice-U_oce |
       !!                - update the modulus of stress at ocean surface
-      !!                      taum = frld * taum + (1-frld) * tmod_io * | U_ice-U_oce |
+      !!                      taum = (1-a) * taum + a * tmod_io * | U_ice-U_oce |
       !!              * at each ocean time step (every kt): 
       !!                  compute linearized ice-ocean stresses as
@@ -336,13 +334,9 @@
       !
       IF( .NOT. ln_rstart ) THEN
-         !                                      ! embedded sea ice
-         IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
-            snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )
-            snwice_mass_b(:,:) = snwice_mass(:,:)
-         ELSE
-            snwice_mass  (:,:) = 0._wp          ! no mass exchanges
-            snwice_mass_b(:,:) = 0._wp          ! no mass exchanges
-         ENDIF
-         IF( nn_ice_embd == 2 ) THEN            ! full embedment (case 2) deplete the initial ssh below sea-ice area
+         !
+         snwice_mass  (:,:) = tmask(:,:,1) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )   ! snow+ice mass
+         snwice_mass_b(:,:) = snwice_mass(:,:)
+         !
+         IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
             sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
             sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 8313)
@@ -172,5 +172,6 @@
             ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- !
             zqld =  tmask(ji,jj,1) * rdt_ice *  &
-               &    ( pfrld(ji,jj) * qsr_oce(ji,jj) * frq_m(ji,jj) + pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
+               &    ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) +  &
+               &      ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) )
 
             ! --- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- !
@@ -210,5 +211,5 @@
       IF( .NOT. ln_limdO )  qlead(:,:) = 0._wp
       ! In case we bypass growing/melting from top and bottom: we suppose ice is impermeable => ocean is isolated from atmosphere
-      IF( .NOT. ln_limdH )  hfx_in(:,:) = pfrld(:,:) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
+      IF( .NOT. ln_limdH )  hfx_in(:,:) = ( 1._wp - at_i_b(:,:) ) * ( qns_oce(:,:) + qsr_oce(:,:) ) + qemp_oce(:,:)
       IF( .NOT. ln_limdH )  fhtur (:,:) = 0._wp  ;  fhld  (:,:) = 0._wp
 
@@ -221,9 +222,9 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            hfx_out(ji,jj) =   pfrld(ji,jj) * qns_oce(ji,jj) + qemp_oce(ji,jj)  &  ! Non solar heat flux received by the ocean               
-               &             - qlead(ji,jj) * r1_rdtice                         &  ! heat flux taken from the ocean where there is open water ice formation
-               &             - at_i(ji,jj) * fhtur(ji,jj)                       &  ! heat flux taken by turbulence
-               &             - at_i(ji,jj) *  fhld(ji,jj)                          ! heat flux taken during bottom growth/melt 
-                                                                                   !    (fhld should be 0 while bott growth)
+            hfx_out(ji,jj) = ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj)  &  ! Non solar heat flux received by the ocean               
+               &             - qlead(ji,jj) * r1_rdtice                                    &  ! heat flux taken from the ocean where there is open water ice formation
+               &             - at_i(ji,jj) * fhtur(ji,jj)                                  &  ! heat flux taken by turbulence
+               &             - at_i(ji,jj) *  fhld(ji,jj)                                     ! heat flux taken during bottom growth/melt 
+                                                                                              !    (fhld should be 0 while bott growth)
          END DO
       END DO
@@ -590,7 +591,7 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicethd/ rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, ln_it_qnsice, nn_monocat, rn_cdsn,  &
-         &                ln_limdH, rn_betas,                                                                     &
-         &                ln_limdA, rn_beta, rn_dmin,                                                             &
+      NAMELIST/namicethd/ rn_kappa_i, nn_ice_thcon, ln_dqnsice, rn_cdsn,                                  &
+         &                ln_limdH, rn_betas,                                                             &
+         &                ln_limdA, rn_beta, rn_dmin,                                                     &
          &                ln_limdO, rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb, rn_himin
       !!-------------------------------------------------------------------
@@ -605,9 +606,4 @@
       IF(lwm) WRITE ( numoni, namicethd )
       !
-      IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN 
-         nn_monocat = 0
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
-      ENDIF
       !
       IF(lwp) THEN                          ! control print
@@ -616,9 +612,6 @@
          WRITE(numout,*)'   -- limthd_dif --'
          WRITE(numout,*)'      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i
-         WRITE(numout,*)'      maximal n. of iter. for heat diffusion computation      nn_conv_dif  = ', nn_conv_dif
-         WRITE(numout,*)'      maximal err. on T for heat diffusion computation        rn_terr_dif  = ', rn_terr_dif
          WRITE(numout,*)'      switch for comp. of thermal conductivity in the ice     nn_ice_thcon = ', nn_ice_thcon
-         WRITE(numout,*)'      iterate the surface non-solar flux (T) or not (F)       ln_it_qnsice = ', ln_it_qnsice
-         WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat
+         WRITE(numout,*)'      change the surface non-solar flux with Tsu or not       ln_dqnsice   = ', ln_dqnsice
          WRITE(numout,*)'      thermal conductivity of the snow                        rn_cdsn      = ', rn_cdsn
          WRITE(numout,*)'   -- limthd_dh --'
@@ -641,4 +634,6 @@
       ENDIF
       !
+      IF ( rn_hnewice < rn_himin )   CALL ctl_stop( 'STOP', 'lim_thd_init : rn_hnewice should be >= rn_himin' )
+      !
    END SUBROUTINE lim_thd_init
 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 8313)
@@ -702,5 +702,5 @@
    !!--------------------------------------------------------------------------
    SUBROUTINE lim_thd_snwblow_2d( pin, pout )
-      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( pfrld or (1. - a_i_b) )
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) :: pin   ! previous fraction lead ( 1. - a_i_b )
       REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
       pout = ( 1._wp - ( pin )**rn_betas )
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 8313)
@@ -92,10 +92,11 @@
 
       !! * Local variables
-      INTEGER ::   ji          ! spatial loop index
-      INTEGER ::   ii, ij      ! temporary dummy loop index
-      INTEGER ::   numeq       ! current reference number of equation
-      INTEGER ::   jk       ! vertical dummy loop index 
-      INTEGER ::   nconv       ! number of iterations in iterative procedure
+      INTEGER ::   ji             ! spatial loop index
+      INTEGER ::   ii, ij         ! temporary dummy loop index
+      INTEGER ::   numeq          ! current reference number of equation
+      INTEGER ::   jk             ! vertical dummy loop index 
       INTEGER ::   minnumeqmin, maxnumeqmax
+      INTEGER ::   iconv          ! number of iterations in iterative procedure
+      INTEGER ::   iconv_max = 50 ! max number of iterations in iterative procedure
       
       INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation
@@ -109,7 +110,8 @@
       REAL(wp) ::   zkimin    =  0.10_wp      ! minimum ice thermal conductivity
       REAL(wp) ::   ztsu_err  =  1.e-5_wp     ! range around which t_su is considered as 0°C 
-      REAL(wp) ::   ztmelt_i    ! ice melting temperature
-      REAL(wp) ::   zerritmax   ! current maximal error on temperature 
+      REAL(wp) ::   ztmelt_i                  ! ice melting temperature
       REAL(wp) ::   zhsu
+      REAL(wp) ::   zdti_max                  ! current maximal error on temperature 
+      REAL(wp) ::   zdti_bnd = 1.e-4_wp       ! maximal authorized error on temperature 
       
       REAL(wp), POINTER, DIMENSION(:)     ::   isnow       ! switch for presence (1) or absence (0) of snow
@@ -122,5 +124,5 @@
       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(:)     ::   zdti        ! 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
@@ -168,5 +170,5 @@
       CALL wrk_alloc( jpij, numeqmin, numeqmax )
       CALL wrk_alloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw )
-      CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe )
+      CALL wrk_alloc( jpij, zf, dzf, zqns_ice_b, zdti, 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 )
@@ -290,5 +292,5 @@
          ztsubit(ji) =  t_su_1d(ji)                              ! temperature at the previous iter
          t_su_1d(ji) =  MIN( t_su_1d(ji), rt0 - ztsu_err )       ! necessary
-         zerrit (ji) =  1000._wp                                 ! initial value of error
+         zdti   (ji) =  1000._wp                                 ! initial value of error
       END DO
 
@@ -305,10 +307,10 @@
       END DO
 
-      nconv     =  0           ! number of iterations
-      zerritmax =  1000._wp    ! maximal value of error on all points
-
-      DO WHILE ( zerritmax > rn_terr_dif .AND. nconv < nn_conv_dif )
-         !
-         nconv = nconv + 1
+      iconv    =  0           ! number of iterations
+      zdti_max =  1000._wp    ! maximal value of error on all points
+
+      DO WHILE ( zdti_max > zdti_bnd .AND. iconv < iconv_max )
+         !
+         iconv = iconv + 1
          !
          !------------------------------------------------------------------------------|
@@ -452,5 +454,5 @@
          !------------------------------------------------------------------------------|
          !
-         IF ( ln_it_qnsice ) THEN 
+         IF ( ln_dqnsice ) THEN 
             DO ji = kideb , kiut
                ! update of the non solar flux according to the update in T_su
@@ -703,8 +705,8 @@
          !
          ! check that nowhere it has started to melt
-         ! zerrit(ji) is a measure of error, it has to be under terr_dif
+         ! zdti(ji) is a measure of error, it has to be under zdti_bnd
          DO ji = kideb , kiut
             t_su_1d(ji) =  MAX(  MIN( t_su_1d(ji) , rt0 ) , 190._wp  )
-            zerrit(ji)  =  ABS( t_su_1d(ji) - ztsubit(ji) )     
+            zdti   (ji) =  ABS( t_su_1d(ji) - ztsubit(ji) )     
          END DO
 
@@ -712,5 +714,5 @@
             DO ji = kideb , kiut
                t_s_1d(ji,jk) = MAX(  MIN( t_s_1d(ji,jk), rt0 ), 190._wp  )
-               zerrit(ji)    = MAX( zerrit(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) )
+               zdti  (ji)    = MAX( zdti(ji), ABS( t_s_1d(ji,jk) - ztstemp(ji,jk) ) )
             END DO
          END DO
@@ -720,5 +722,5 @@
                ztmelt_i      = -tmut * s_i_1d(ji,jk) + rt0 
                t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), 190._wp )
-               zerrit(ji)    =  MAX( zerrit(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) )
+               zdti  (ji)    =  MAX( zdti(ji), ABS( t_i_1d(ji,jk) - ztitemp(ji,jk) ) )
             END DO
          END DO
@@ -726,9 +728,9 @@
          ! Compute spatial maximum over all errors
          ! note that this could be optimized substantially by iterating only the non-converging points
-         zerritmax = 0._wp
+         zdti_max = 0._wp
          DO ji = kideb, kiut
-            zerritmax = MAX( zerritmax, zerrit(ji) )   
-         END DO
-         IF( lk_mpp ) CALL mpp_max( zerritmax, kcom=ncomm_ice )
+            zdti_max = MAX( zdti_max, zdti(ji) )   
+         END DO
+         IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice )
 
       END DO  ! End of the do while iterative procedure
@@ -745,6 +747,6 @@
 
       IF( ln_limctl .AND. lwp ) THEN
-         WRITE(numout,*) ' zerritmax : ', zerritmax
-         WRITE(numout,*) ' nconv     : ', nconv
+         WRITE(numout,*) ' zdti_max : ', zdti_max
+         WRITE(numout,*) ' iconv    : ', iconv
       ENDIF
 
@@ -772,5 +774,5 @@
 
       ! --- diagnose the change in non-solar flux due to surface temperature change --- !
-      IF ( ln_it_qnsice ) THEN
+      IF ( ln_dqnsice ) THEN
          DO ji = kideb, kiut
             hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji) 
@@ -810,5 +812,5 @@
       CALL wrk_dealloc( jpij, numeqmin, numeqmax )
       CALL wrk_dealloc( jpij, isnow, ztsub, ztsubit, zh_i, zh_s, zfsw )
-      CALL wrk_dealloc( jpij, zf, dzf, zqns_ice_b, zerrit, zdifcase, zftrice, zihic, zghe )
+      CALL wrk_dealloc( jpij, zf, dzf, zqns_ice_b, zdti, 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 )
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 8313)
@@ -67,5 +67,5 @@
       !!             - Updating ice internal temperature
       !!             - Computation of variation of ice volume and mass
-      !!             - Computation of frldb after lateral accretion and 
+      !!             - Computation of a_i after lateral accretion and 
       !!               update ht_s_1d, ht_i_1d and tbif_1d(:,:)      
       !!------------------------------------------------------------------------
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 8313)
@@ -109,11 +109,10 @@
       !----------------------------------------
       ! fluxes 
-      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD)
-      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface
-      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface
+      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * ( 1._wp - at_i_b(:,:) ) )                      !     solar flux at ocean surface
+      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) )      ! non-solar flux at ocean surface
       IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface
       IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface
       IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice
-      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) )  
+      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * ( 1._wp - at_i_b(:,:) ) + qemp_oce(:,:) )
       IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   &
          &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) )
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 8313)
@@ -82,5 +82,4 @@
 
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d    !: <==> the 2D  sprecip
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   frld_1d       !: <==> the 2D  frld
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_1d        !: <==> the 2D  at_i
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhtur_1d      !: <==> the 2D  fhtur
@@ -153,5 +152,5 @@
       !
       ii = ii + 1
-      ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_1d    (jpij) ,                     &
+      ALLOCATE( sprecip_1d (jpij) , at_i_1d    (jpij) ,                     &
          &      fhtur_1d   (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) , wfx_snw_sum_1d(jpij) , &
          &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d (jpij) , wfx_bom_1d(jpij) ,  &
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/ASM/asmbkg.F90	(revision 8313)
@@ -146,9 +146,9 @@
             CALL iom_rstput( kt, nitdin_r, inum, 'sshn'   , sshn              )
 #if defined key_lim3
-            IF( nn_ice == 2  .OR.  nn_ice == 3 ) THEN
-	            IF( ALLOCATED(frld) ) THEN
-                  CALL iom_rstput( kt, nitdin_r, inum, 'iceconc', 1._wp - frld(:,:)   )
+            IF( nn_ice == 3 ) THEN
+	            IF( ALLOCATED(at_i) ) THEN
+                  CALL iom_rstput( kt, nitdin_r, inum, 'iceconc', at_i(:,:)   )
                ELSE
-		            CALL ctl_warn('Ice concentration not written to background as ice variable frld not allocated on this timestep')
+		            CALL ctl_warn('Ice concentration not written to background as ice variable at_i not allocated on this timestep')
 	            ENDIF
             ENDIF
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/ASM/asminc.F90	(revision 8313)
@@ -36,5 +36,5 @@
    USE lib_mpp          ! MPP library
 #if defined key_lim3
-   USE ice,   ONLY: htm_i, at_i, pfrld, phicif
+   USE ice,   ONLY: htm_i, at_i, at_i_b
 #endif
    USE sbc_oce          ! Surface boundary condition variables.
@@ -837,6 +837,6 @@
             zohicif(:,:) = htm_i(:,:)
             !
-            at_i(:,:) = 1. - MIN( MAX( 1.-at_i(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
-            pfrld     =      MIN( MAX( pfrld(:,:)   - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
+            at_i  (:,:) = 1. - MIN( MAX( 1.-at_i  (:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
+            at_i_b(:,:) = 1. - MIN( MAX( 1.-at_i_b(:,:) - seaice_bkginc(:,:) * zincwgt, 0.0_wp), 1.0_wp)
             fr_i(:,:) = at_i(:,:)        ! adjust ice fraction
             !
@@ -852,5 +852,4 @@
             ! nudge ice depth
             htm_i (:,:) = htm_i (:,:) + zhicifinc(:,:)
-            phicif(:,:) = phicif(:,:) + zhicifinc(:,:)       
             !
             ! seaice salinity balancing (to add)
@@ -889,5 +888,5 @@
             ! Initialize the now fields the background + increment
             at_i(:,:) = 1. - MIN( MAX( 1.-at_i(:,:) - seaice_bkginc(:,:), 0.0_wp), 1.0_wp)
-            pfrld(:,:) = 1.-at_i(:,:) 
+            at_i_b(:,:) = at_i(:,:) 
             fr_i(:,:) = at_i(:,:)        ! adjust ice fraction
             !
@@ -903,5 +902,4 @@
             ! nudge ice depth
             htm_i (:,:) = htm_i (:,:) + zhicifinc(:,:)
-            phicif(:,:) = phicif(:,:)       
             !
             ! seaice salinity balancing (to add)
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/DYN/dynspg.F90	(revision 8313)
@@ -89,5 +89,5 @@
       IF(      ln_apr_dyn                                                &   ! atmos. pressure
          .OR.  ( .NOT.ln_dynspg_ts .AND. (ln_tide_pot .AND. ln_tide) )   &   ! tide potential (no time slitting)
-         .OR.  nn_ice_embd == 2  ) THEN                                      ! embedded sea-ice
+         .OR.  ln_ice_embd ) THEN                                            ! embedded sea-ice
          !
          DO jj = 2, jpjm1
@@ -123,5 +123,5 @@
          ENDIF
          !
-         IF( nn_ice_embd == 2 ) THEN          !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
+         IF( ln_ice_embd ) THEN          !== embedded sea ice: Pressure gradient due to snow-ice mass ==!
             CALL wrk_alloc( jpi,jpj,   zpice )
             !                                            
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/OBS/diaobs.F90	(revision 8313)
@@ -539,5 +539,5 @@
 #if defined  key_lim3
       USE ice, ONLY : &            ! LIM3 Ice model variables
-         & frld
+         & at_i
 #endif
       IMPLICIT NONE
@@ -564,5 +564,5 @@
          & zgphi2                  ! Model latitudes for prof variable 2
 #if ! defined key_lim3
-      REAL(wp), POINTER, DIMENSION(:,:) :: frld
+      REAL(wp), POINTER, DIMENSION(:,:) :: at_i
 #endif
       LOGICAL :: llnightav        ! Logical for calculating night-time average
@@ -579,5 +579,5 @@
       CALL wrk_alloc( jpi, jpj, zgphi2 )
 #if ! defined key_lim3
-      CALL wrk_alloc(jpi,jpj,frld) 
+      CALL wrk_alloc(jpi,jpj,at_i) 
 #endif
 
@@ -591,8 +591,8 @@
 
       !-----------------------------------------------------------------------
-      ! No LIM => frld == 0.0_wp
+      ! No LIM => at_i == 0.0_wp
       !-----------------------------------------------------------------------
 #if ! defined key_lim3
-      frld(:,:) = 0.0_wp
+      at_i(:,:) = 0.0_wp
 #endif
       !-----------------------------------------------------------------------
@@ -674,5 +674,5 @@
                   CYCLE
                ELSE
-                  zsurfvar(:,:) = 1._wp - frld(:,:)
+                  zsurfvar(:,:) = at_i(:,:)
                ENDIF
 
@@ -699,5 +699,5 @@
       CALL wrk_dealloc( jpi, jpj, zgphi2 )
 #if ! defined key_lim3
-      CALL wrk_dealloc(jpi,jpj,frld)
+      CALL wrk_dealloc(jpi,jpj,at_i)
 #endif
 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 8313)
@@ -47,8 +47,7 @@
    LOGICAL , PUBLIC ::   ln_apr_dyn     !: Atmospheric pressure forcing used on dynamics (ocean & ice)
    INTEGER , PUBLIC ::   nn_ice         !: flag for ice in the surface boundary condition (=0/1/2/3)
-   INTEGER , PUBLIC ::   nn_ice_embd    !: flag for levitating/embedding sea-ice in the ocean
-   !                                             !: =0 levitating ice (no mass exchange, concentration/dilution effect)
-   !                                             !: =1 levitating ice with mass and salt exchange but no presure effect
-   !                                             !: =2 embedded sea-ice (full salt and mass exchanges and pressure)
+   LOGICAL , PUBLIC ::   ln_ice_embd    !: flag for levitating/embedding sea-ice in the ocean
+   !                                             !: =F levitating ice with mass and salt exchange but no presure effect
+   !                                             !: =T embedded sea-ice (full salt and mass exchanges and pressure)
    INTEGER , PUBLIC ::   nn_components  !: flag for sbc module (including sea-ice) coupling mode (see component definition below) 
    INTEGER , PUBLIC ::   nn_limflx      !: LIM3 Multi-category heat flux formulation
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90	(revision 8313)
@@ -40,5 +40,5 @@
    USE lib_fortran    ! to use key_nosignedzero
 #if defined key_lim3
-   USE ice     , ONLY :   u_ice, v_ice, jpl, pfrld, a_i_b, at_i_b
+   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b
    USE limthd_dh      ! for CALL lim_thd_snwblow
 #endif
@@ -789,11 +789,11 @@
       ! --- evaporation minus precipitation --- !
       zsnw(:,:) = 0._wp
-      CALL lim_thd_snwblow( pfrld, zsnw )  ! snow distribution over ice after wind blowing
-      emp_oce(:,:) = pfrld(:,:) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
+      CALL lim_thd_snwblow( (1.-at_i_b(:,:)), zsnw )  ! snow distribution over ice after wind blowing
+      emp_oce(:,:) = ( 1._wp - at_i_b(:,:) ) * zevap(:,:) - ( tprecip(:,:) - sprecip(:,:) ) - sprecip(:,:) * (1._wp - zsnw )
       emp_ice(:,:) = SUM( a_i_b(:,:,:) * evap_ice(:,:,:), dim=3 ) - sprecip(:,:) * zsnw
       emp_tot(:,:) = emp_oce(:,:) + emp_ice(:,:)
 
       ! --- heat flux associated with emp --- !
-      qemp_oce(:,:) = - pfrld(:,:) * zevap(:,:) * sst_m(:,:) * rcp                               & ! evap at sst
+      qemp_oce(:,:) = - ( 1._wp - at_i_b(:,:) ) * zevap(:,:) * sst_m(:,:) * rcp                  & ! evap at sst
          &          + ( tprecip(:,:) - sprecip(:,:) ) * ( sf(jp_tair)%fnow(:,:,1) - rt0 ) * rcp  & ! liquid precip at Tair
          &          +   sprecip(:,:) * ( 1._wp - zsnw ) *                                        & ! solid precip at min(Tair,Tsnow)
@@ -803,6 +803,7 @@
 
       ! --- total solar and non solar fluxes --- !
-      qns_tot(:,:) = pfrld(:,:) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 ) + qemp_ice(:,:) + qemp_oce(:,:)
-      qsr_tot(:,:) = pfrld(:,:) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
+      qns_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qns_oce(:,:) + SUM( a_i_b(:,:,:) * qns_ice(:,:,:), dim=3 )  &
+         &           + qemp_ice(:,:) + qemp_oce(:,:)
+      qsr_tot(:,:) = ( 1._wp - at_i_b(:,:) ) * qsr_oce(:,:) + SUM( a_i_b(:,:,:) * qsr_ice(:,:,:), dim=3 )
 
       ! --- heat content of precip over ice in J/m3 (to be used in 1D-thermo) --- !
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 8313)
@@ -1557,8 +1557,8 @@
       !!
       !! ** Details
-      !!             qns_tot = pfrld * qns_oce + ( 1 - pfrld ) * qns_ice   => provided
+      !!             qns_tot = (1-a) * qns_oce + a * qns_ice               => provided
       !!                     + qemp_oce + qemp_ice                         => recalculated and added up to qns
       !!
-      !!             qsr_tot = pfrld * qsr_oce + ( 1 - pfrld ) * qsr_ice   => provided
+      !!             qsr_tot = (1-a) * qsr_oce + a * qsr_ice               => provided
       !!
       !!             emp_tot = emp_oce + emp_ice                           => calving is provided and added to emp_tot (and emp_oce).
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_cice.F90	(revision 8313)
@@ -230,16 +230,12 @@
       CALL lbc_lnk ( fr_iv , 'V', 1. )
 
-      !                                      ! embedded sea ice
-      IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
-         CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
-         CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
-         snwice_mass  (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:)  )
-         snwice_mass_b(:,:) = snwice_mass(:,:)
-      ELSE
-         snwice_mass  (:,:) = 0.0_wp         ! no mass exchanges
-         snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges
-      ENDIF
+      ! set the snow+ice mass
+      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
+      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
+      snwice_mass  (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:)  )
+      snwice_mass_b(:,:) = snwice_mass(:,:)
+
       IF( .NOT.ln_rstart ) THEN
-         IF( nn_ice_embd == 2 ) THEN            ! full embedment (case 2) deplete the initial ssh below sea-ice area
+         IF( ln_ice_embd ) THEN            ! embedded sea-ice: deplete the initial ssh below sea-ice area
             sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
             sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
@@ -473,5 +469,5 @@
       CALL nemo2cice(ztmp,vocn,'F', -1. )
 
-      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==!
+      IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
           !
           ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
@@ -676,12 +672,10 @@
       CALL lbc_lnk ( fr_iv , 'V', 1. )
 
-      !                                      ! embedded sea ice
-      IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
-         CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
-         CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
-         snwice_mass  (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:)  )
-         snwice_mass_b(:,:) = snwice_mass(:,:)
-         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
-      ENDIF
+      ! set the snow+ice mass
+      CALL cice2nemo(vsno(:,:,:),ztmp1,'T', 1. )
+      CALL cice2nemo(vice(:,:,:),ztmp2,'T', 1. )
+      snwice_mass  (:,:) = ( rhosn * ztmp1(:,:) + rhoic * ztmp2(:,:)  )
+      snwice_mass_b(:,:) = snwice_mass(:,:)
+      snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / dt
 
 ! Release work space
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 8313)
@@ -196,8 +196,5 @@
                                       CALL sbc_lim_bef                        
                                       CALL lim_var_glo2eqv     ! ht_i and ht_s for ice albedo calculation
-                                      CALL lim_var_agg(1)      ! at_i for coupling (via pfrld) 
-         !
-         pfrld(:,:)   = 1._wp - at_i(:,:)
-         phicif(:,:)  = vt_i(:,:)
+                                      CALL lim_var_agg(1)      ! at_i for coupling 
 
          !------------------------------------------------------!
@@ -221,10 +218,10 @@
                alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
                                       CALL blk_ice_flx( t_su, alb_ice )
-               IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
+               IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( p_frld=(1.-at_i_b), palbi=alb_ice, psst=sst_m, pist=t_su )
                IF( nn_limflx /= 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
             CASE ( jp_purecpl )
                ! albedo depends on cloud fraction because of non-linear spectral effects
                alb_ice(:,:,:) = ( 1. - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:)
-                                      CALL sbc_cpl_ice_flx( p_frld=pfrld, palbi=alb_ice, psst=sst_m, pist=t_su )
+                                      CALL sbc_cpl_ice_flx( p_frld=(1.-at_i_b), palbi=alb_ice, psst=sst_m, pist=t_su )
                IF( nn_limflx == 2 )   CALL ice_lim_flx( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_limflx )
          END SELECT
@@ -374,5 +371,5 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir,   &
+      NAMELIST/namicerun/ jpl, nlay_i, nlay_s, nn_monocat, rn_amax_n, rn_amax_s, cn_icerst_in, cn_icerst_indir,   &
          &                cn_icerst_out, cn_icerst_outdir, ln_limthd, ln_limdyn, nn_limdyn, rn_uice, rn_vice  
       NAMELIST/namicediag/ ln_limdiachk, ln_limdiahsb, ln_limctl, iiceprt, jiceprt  
@@ -401,17 +398,18 @@
          WRITE(numout,*) 'lim_run_init : ice share parameters for dynamics/advection/thermo of sea-ice'
          WRITE(numout,*) ' ~~~~~~'
-         WRITE(numout,*) '   number of ice  categories                               = ', jpl
-         WRITE(numout,*) '   number of ice  layers                                   = ', nlay_i
-         WRITE(numout,*) '   number of snow layers                                   = ', nlay_s
-         WRITE(numout,*) '   maximum ice concentration for NH                        = ', rn_amax_n 
-         WRITE(numout,*) '   maximum ice concentration for SH                        = ', rn_amax_s
-         WRITE(numout,*) '   Ice thermodynamics (T) or not (F)            ln_limthd  = ', ln_limthd
-         WRITE(numout,*) '   Ice dynamics       (T) or not (F)            ln_limdyn  = ', ln_limdyn
-         WRITE(numout,*) '     (ln_limdyn=T) Ice dynamics switch          nn_limdyn  = ', nn_limdyn
+         WRITE(numout,*) '   number of ice  categories                              jpl    = ', jpl
+         WRITE(numout,*) '   number of ice  layers                                  nlay_i = ', nlay_i
+         WRITE(numout,*) '   number of snow layers                                  nlay_s = ', nlay_s
+         WRITE(numout,*) '   virtual ITD mono-category param (1-4) or not (0)   nn_monocat = ', nn_monocat
+         WRITE(numout,*) '   maximum ice concentration for NH                              = ', rn_amax_n 
+         WRITE(numout,*) '   maximum ice concentration for SH                              = ', rn_amax_s
+         WRITE(numout,*) '   Ice thermodynamics (T) or not (F)                  ln_limthd  = ', ln_limthd
+         WRITE(numout,*) '   Ice dynamics       (T) or not (F)                  ln_limdyn  = ', ln_limdyn
+         WRITE(numout,*) '     (ln_limdyn=T) Ice dynamics switch                nn_limdyn  = ', nn_limdyn
          WRITE(numout,*) '       2: total'
          WRITE(numout,*) '       1: advection only (no diffusion, no ridging/rafting)'
          WRITE(numout,*) '       0: advection only (as 1 + prescribed velocity, bypass rheology)'
-         WRITE(numout,*) '     (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0)     = ', rn_uice
-         WRITE(numout,*) '     (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0)     = ', rn_vice
+         WRITE(numout,*) '     (ln_limdyn=T) prescribed u-vel (case nn_limdyn=0)           = ', rn_uice
+         WRITE(numout,*) '     (ln_limdyn=T) prescribed v-vel (case nn_limdyn=0)           = ', rn_vice
          WRITE(numout,*)
          WRITE(numout,*) '...and ice diagnostics'
@@ -424,4 +422,13 @@
       ENDIF
       !
+      IF ( ( jpl > 1 ) .AND. ( nn_monocat == 1 ) ) THEN
+         nn_monocat = 0
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen'
+      ENDIF
+      IF ( ( jpl == 1 ) .AND. ( nn_monocat == 0 ) ) THEN
+         CALL ctl_stop( 'STOP', 'lim_run_init : if jpl=1 then nn_monocat should be between 1 and 4' )
+      ENDIF
+      !
       ! sea-ice timestep and inverse
       rdt_ice   = REAL(nn_fsbc) * rdt  
@@ -449,5 +456,5 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namiceitd/ nn_catbnd, rn_himean
+      NAMELIST/namiceitd/ rn_himean
       !
       INTEGER  ::   jl                   ! dummy loop index
@@ -469,6 +476,5 @@
          WRITE(numout,*) 'lim_itd_init : Initialization of ice cat distribution '
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   shape of ice categories distribution                          nn_catbnd = ', nn_catbnd
-         WRITE(numout,*) '   mean ice thickness in the domain (only active if nn_catbnd=2) rn_himean = ', rn_himean
+         WRITE(numout,*) '   mean ice thickness in the domain            rn_himean = ', rn_himean
       ENDIF
       !
@@ -479,25 +485,13 @@
       hi_max(:) = 0._wp
       !
-      SELECT CASE ( nn_catbnd  )    ! type of ice categories distribution
-      !
-      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_wp
-         zhmax  = 3._wp * rn_himean
-         DO jl = 1, jpl
-            znum = jpl * ( zhmax+1 )**zalpha
-            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
-            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
-         END DO
-         !
-      END SELECT
+      !==  h^(-alpha) function  ==!
+      zalpha = 0.05_wp
+      zhmax  = 3._wp * rn_himean
+      DO jl = 1, jpl
+         znum = jpl * ( zhmax+1 )**zalpha
+         zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp )
+         hi_max(jl) = ( znum / zden )**(1./zalpha) - 1
+      END DO
+      !
       !
       DO jl = 1, jpl                ! mean thickness by category
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 8313)
@@ -90,5 +90,5 @@
          &             ln_usr   , ln_flx   , ln_blk       ,                          &
          &             ln_cpl   , ln_mixcpl, nn_components, nn_limflx,               &
-         &             nn_ice   , nn_ice_embd,                                       &
+         &             nn_ice   , ln_ice_embd,                                       &
          &             ln_traqsr, ln_dm2dc ,                                         &
          &             ln_rnf   , nn_fwb   , ln_ssr   , ln_isf    , ln_apr_dyn ,     &
@@ -139,6 +139,6 @@
          WRITE(numout,*) '         Multicategory heat flux formulation (LIM3) nn_limflx     = ', nn_limflx
          WRITE(numout,*) '      Sea-ice : '
-         WRITE(numout,*) '         ice management in the sbc (=0/1/2/3)       nn_ice        = ', nn_ice
-         WRITE(numout,*) '         ice-ocean embedded/levitating (=0/1/2)     nn_ice_embd   = ', nn_ice_embd
+         WRITE(numout,*) '         ice management in the sbc (=0/1/3/4)       nn_ice        = ', nn_ice
+         WRITE(numout,*) '         ice embedded into ocean                    ln_ice_embd   = ', ln_ice_embd
          WRITE(numout,*) '      Misc. options of sbc : '
          WRITE(numout,*) '         Light penetration in temperature Eq.       ln_traqsr     = ', ln_traqsr
@@ -199,8 +199,6 @@
       CASE( 1 )                        !- Ice-cover climatology ("Ice-if" model)  
       CASE( 3 )                        !- LIM3 ice model
-         IF( nn_ice_embd == 0            )   CALL ctl_stop( 'sbc_init : LIM3 sea-ice models require nn_ice_embd = 1 or 2' )
       CASE( 4 )                        !- CICE ice model
          IF( .NOT.( ln_blk .OR. ln_cpl ) )   CALL ctl_stop( 'sbc_init : CICE sea-ice model requires ln_blk or ln_cpl = T' )
-         IF( nn_ice_embd == 0            )   CALL ctl_stop( 'sbc_init : CICE sea-ice models require nn_ice_embd = 1 or 2' )
          IF( lk_agrif                    )   CALL ctl_stop( 'sbc_init : CICE sea-ice model not currently available with AGRIF' ) 
       CASE DEFAULT                     !- not supported
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 8313)
@@ -48,5 +48,5 @@
    LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
-   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem)
+   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model
    INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
    REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands)
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/TOP_SRC/TRP/trcsbc.F90	(revision 8313)
@@ -63,5 +63,5 @@
       INTEGER  ::   ji, jj, jn                                     ! dummy loop indices
       REAL(wp) ::   zse3t, zrtrn, zratio, zfact                    ! temporary scalars
-      REAL(wp) ::   zswitch, zftra, zcd, zdtra, ztfx, ztra         ! temporary scalars
+      REAL(wp) ::   zftra, zcd, zdtra, ztfx, ztra                  ! temporary scalars
       CHARACTER (len=22) :: charout
       REAL(wp), POINTER, DIMENSION(:,:  ) :: zsfx
@@ -77,10 +77,4 @@
       !
       zrtrn = 1.e-15_wp
-
-      SELECT CASE( nn_ice_embd )         ! levitating or embedded sea-ice option
-         CASE( 0    )   ;   zswitch = 1  ! (0) standard levitating sea-ice : salt exchange only
-         CASE( 1, 2 )   ;   zswitch = 0  ! (1) levitating sea-ice: salt and volume exchange but no pressure effect                                
-      !                                  ! (2) embedded sea-ice : salt and volume fluxes and pressure
-      END SELECT
 
       IF( kt == nittrc000 ) THEN
@@ -145,5 +139,5 @@
                   ! tracer flux only       : add concentration dilution term in net tracer flux, no F-M in volume flux
                   ! tracer and mass fluxes : no concentration dilution term in net tracer flux, F-M term in volume flux
-                  ztfx  = zftra + zswitch * zcd                ! net tracer flux (+C/D if no ice/ocean mass exchange)
+                  ztfx  = zftra                             ! net tracer flux
    
                   zdtra = r1_rau0 * ( ztfx + zsfx(ji,jj) * trn(ji,jj,1,jn) ) 
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 8312)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/TOP_SRC/oce_trc.F90	(revision 8313)
@@ -63,5 +63,4 @@
    USE sbc_oce , ONLY :   ln_rnf     =>    ln_rnf     !: runoffs / runoff mouths
    USE sbc_oce , ONLY :   fr_i       =>    fr_i       !: ice fraction (between 0 to 1)
-   USE sbc_oce , ONLY :   nn_ice_embd => nn_ice_embd  !: flag for  levitating/embedding sea-ice in the ocean
    USE sbc_oce , ONLY :   atm_co2    =>    atm_co2    !  atmospheric pCO2
    USE traqsr  , ONLY :   rn_abs     =>    rn_abs     !: fraction absorbed in the very near surface
