Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4155)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4220)
@@ -399,5 +399,5 @@
    REAL(wp)              , PUBLIC ::   cao           = 1.00e-3            !: atmospheric drag over ocean
    REAL(wp)              , PUBLIC ::   amax          = 0.99               !: maximum ice concentration
-   !                                                                      !
+   !
    !!--------------------------------------------------------------------------
    !! * Ice diagnostics
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4155)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4220)
@@ -87,7 +87,4 @@
       CALL lim_itd_ini                 ! ice thickness distribution initialization
       !
-      CALL lim_sbc_init                ! ice surface boundary condition   
-
-
       !                                ! Initial sea-ice state
       IF( .NOT.ln_rstart ) THEN              ! start from rest
@@ -104,5 +101,8 @@
       ENDIF
       !
+      CALL lim_sbc_init                ! ice surface boundary condition   
+      !
       fr_i(:,:) = at_i(:,:)           ! initialisation of sea-ice fraction
+      tn_ice(:,:,:) = t_su(:,:,:)
       !
       nstart = numit  + nn_fsbc      
@@ -239,5 +239,4 @@
       END DO
       !
-      tn_ice(:,:,:) = t_su(:,:,:)
       !
    END SUBROUTINE lim_itd_ini
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4155)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4220)
@@ -123,5 +123,5 @@
 
       REAL(wp) ::   dtevp              ! time step for subcycling
-      REAL(wp) ::   dtotel, ecc2       ! square of yield ellipse eccenticity
+      REAL(wp) ::   dtotel, ecc2, ecci ! square of yield ellipse eccenticity
       REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars
       REAL(wp) ::   zu_ice2, zv_ice1   !
@@ -199,5 +199,5 @@
             zc1(ji,jj)    = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )
 #if defined key_lim3
-            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) * 0.5_wp
+            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj)
 #endif
 #if defined key_lim2
@@ -325,4 +325,5 @@
       !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
       ecc2 = ecc * ecc
+      ecci = 1. / ecc2
 
       !-Initialise stress tensor 
@@ -425,14 +426,17 @@
                !-Calculate stress tensor components zs1 and zs2 
                !-at centre of grid cells (see section 3.5 of CICE user's guide).
-               zs1(ji,jj) = ( zs1(ji,jj) &
-                  &          - dtotel*( ( 1.0 - alphaevp) * zs1(ji,jj) +    &
-                  &            ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) &
-                  * zpresh(ji,jj) ) )                          &       
-                  &        / ( 1.0 + alphaevp * dtotel )
-
-               zs2(ji,jj) = ( zs2(ji,jj)   &
-                  &          - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj) -  &
-                  zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) &
-                  &        / ( 1.0 + alphaevp*ecc2*dtotel )
+               zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   &
+                  &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &       
+                  &          / ( 1._wp + alphaevp * dtotel )
+
+               zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   &
+                             zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   &
+                  &          / ( 1._wp + alphaevp * ecc2 * dtotel )
+
+               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
+               !zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  &
+               !   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
+               !zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  &
+               !   &         / ( 1._wp + dtotel )
 
             END DO
@@ -468,8 +472,12 @@
 
                !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
-               zs12(ji,jj) = ( zs12(ji,jj)      &
-                  &        - dtotel*( (1.0-alphaevp)*ecc2*zs12(ji,jj) - zds(ji,jj) / &
-                  &          ( 2.0*deltac(ji,jj) ) * zpreshc(ji,jj))) &
-                  &         / ( 1.0 + alphaevp*ecc2*dtotel ) 
+               zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  &
+                  &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+                  &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
+
+               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
+               !zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
+               !   &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+               !   &          / ( 1.0 + dtotel ) 
 
             END DO ! ji
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 4155)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 4220)
@@ -17,4 +17,5 @@
    !!----------------------------------------------------------------------
    USE ice            ! sea-ice variables
+   USE oce     , ONLY :  snwice_mass, snwice_mass_b
    USE par_ice        ! sea-ice parameters
    USE dom_oce        ! ocean domain
@@ -159,12 +160,12 @@
       END DO
 
-      CALL iom_rstput( iter, nitrst, numriw, 'u_ice'     , u_ice      )
-      CALL iom_rstput( iter, nitrst, numriw, 'v_ice'     , v_ice      )
-      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq'     , fsbbq      )
-      CALL iom_rstput( iter, nitrst, numriw, 'iatte'     , iatte      ) ! clem modif
-      CALL iom_rstput( iter, nitrst, numriw, 'oatte'     , oatte      ) ! clem modif
-      CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
-      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
-      CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
+      CALL iom_rstput( iter, nitrst, numriw, 'u_ice'        , u_ice      )
+      CALL iom_rstput( iter, nitrst, numriw, 'v_ice'        , v_ice      )
+      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq'        , fsbbq      )
+      CALL iom_rstput( iter, nitrst, numriw, 'stress1_i'    , stress1_i  )
+      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i'    , stress2_i  )
+      CALL iom_rstput( iter, nitrst, numriw, 'stress12_i'   , stress12_i )
+      CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass )   !clem modif
+      CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b ) !clem modif
 
       DO jl = 1, jpl 
@@ -371,52 +372,5 @@
          CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
          t_su(:,:,jl) = z2d(:,:)
-         tn_ice (:,:,:) = t_su (:,:,:)
-      END DO
-
-      DO jl = 1, jpl 
-         CALL lbc_lnk( smv_i(:,:,jl) , 'T' ,  1. )
-         CALL lbc_lnk( v_i  (:,:,jl) , 'T' ,  1. )
-         CALL lbc_lnk( a_i  (:,:,jl) , 'T' ,  1. )
-      END DO
-
-      ! we first with bulk ice salinity
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zindb          = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - 1.0e-4 ) ) 
-               sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX(v_i(ji,jj,jl),1.0e-6) * zindb
-               ht_i(ji,jj,jl) = v_i(ji,jj,jl)   / MAX(a_i(ji,jj,jl),1.0e-6) * zindb
-            END DO
-         END DO
-      END DO
-
-      DO jk = 1, nlay_i
-         s_i(:,:,jk,:) = sm_i(:,:,:)
-      END DO
-
-      IF( num_sal == 2 ) THEN      ! Salinity profile
-         DO jl = 1, jpl
-            DO jk = 1, nlay_i
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zs_inf        = sm_i(ji,jj,jl)
-                     z_slope_s     = 2._wp * sm_i(ji,jj,jl) / MAX( 0.01_wp , ht_i(ji,jj,jl) )
-                     !- slope of the salinity profile
-                     zs_zero(jk)   = z_slope_s * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) / REAL(nlay_i,wp)
-                     zsmax = 4.5_wp
-                     zsmin = 3.5_wp
-                     IF(     sm_i(ji,jj,jl) < zsmin ) THEN
-                        zalpha = 1._wp
-                     ELSEIF( sm_i(ji,jj,jl) < zsmax ) THEN
-                        zalpha = sm_i(ji,jj,jl) / ( zsmin - zsmax ) + zsmax / ( zsmax - zsmin )
-                     ELSE
-                        zalpha = 0._wp
-                     ENDIF
-                     s_i(ji,jj,jk,jl) = zalpha * zs_zero(jk) + ( 1._wp - zalpha ) * zs_inf
-                  END DO
-               END DO
-            END DO
-         END DO
-      ENDIF
+      END DO
 
       DO jl = 1, jpl 
@@ -440,9 +394,9 @@
       CALL iom_get( numrir, jpdom_autoglo, 'v_ice'     , v_ice      )
       CALL iom_get( numrir, jpdom_autoglo, 'fsbbq'     , fsbbq      )
-      CALL iom_get( numrir, jpdom_autoglo, 'iatte'     , iatte      ) ! clem modif
-      CALL iom_get( numrir, jpdom_autoglo, 'oatte'     , oatte      ) ! clem modif
       CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
       CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
       CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
+      CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass )   !clem modif
+      CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b ) !clem modif
 
       DO jl = 1, jpl 
@@ -568,5 +522,5 @@
       END DO
       !
-      !clem CALL iom_close( numrir )
+      !CALL iom_close( numrir ) !clem: closed in sbcice_lim.F90
       !
       CALL wrk_dealloc( nlay_i, zs_zero )
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4155)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4220)
@@ -38,5 +38,5 @@
    USE cpl_oasis3, ONLY : lk_cpl
    USE traqsr           ! clem: add penetration of solar flux into the calculation of heat budget
-   USE oce,        ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n
+   USE oce,        ONLY : iatte, oatte, sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n
    USE dom_ice,    ONLY : tms
    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
@@ -249,5 +249,5 @@
             ! mass flux at the ocean/ice interface (sea ice fraction)
             zemp_snw = rdm_snw(ji,jj) * r1_rdtice                         ! snow melting = pure water that enters the ocean
-            zfmm     = rdm_ice(ji,jj) * r1_rdtice                         ! Freezing minus mesting  
+            zfmm     = rdm_ice(ji,jj) * r1_rdtice                         ! Freezing minus melting  
 
             fmmflx(ji,jj) = zfmm                                     ! F/M mass flux save at least for biogeochemical model
@@ -416,51 +416,55 @@
       ENDIF
       ! clem modif
-      iatte(:,:) = 1._wp
-      oatte(:,:) = 1._wp
-      !
-      !                                      ! 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  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )
-         snwice_mass_b(:,:) = snwice_mass(:,:)
-      ELSE
-         snwice_mass  (:,:) = 0.0_wp         ! no mass exchanges
-         snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges
-      ENDIF
-      IF( nn_ice_embd == 2  .AND.         &  ! full embedment (case 2) & no restart
-         &  .NOT. ln_rstart ) THEN           ! deplete the initial ssh below sea-ice area
-         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
-         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
-         !
-         ! Note: Changed the initial values of sshb and sshn=>  need to recompute ssh[u,v,f]_[b,n] 
-         !       which were previously set in domvvl
-         IF ( lk_vvl ) THEN            ! Is this necessary? embd 2 should be restricted to vvl only???
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
-                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
-                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
-                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
-                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
-                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
-                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
-                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
-                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
-                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
-                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
-                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
+      IF( .NOT. ln_rstart ) THEN
+         iatte(:,:) = 1._wp
+         oatte(:,:) = 1._wp
+      ENDIF
+      !
+      ! clem: snwice_mass in the restart file now
+      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  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )
+            snwice_mass_b(:,:) = snwice_mass(:,:)
+         ELSE
+            snwice_mass  (:,:) = 0.0_wp         ! no mass exchanges
+            snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges
+         ENDIF
+         IF( nn_ice_embd == 2 ) THEN            ! full embedment (case 2) deplete the initial ssh below sea-ice area
+            sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
+            sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
+            !
+            ! Note: Changed the initial values of sshb and sshn=>  need to recompute ssh[u,v,f]_[b,n] 
+            !       which were previously set in domvvl
+            IF ( lk_vvl ) THEN            ! Is this necessary? embd 2 should be restricted to vvl only???
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
+                     zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
+                     zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
+                     zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
+                     sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
+                        &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
+                     sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
+                        &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
+                     sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
+                        &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
+                     sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
+                        &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
+                  END DO
                END DO
-            END DO
-            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
-            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
-                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
-                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
-                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
+               CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
+               CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1      ! NO Vector Opt.
+                     sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
+                          &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
+                          &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
+                          &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
+                  END DO
                END DO
-            END DO
-            CALL lbc_lnk( sshf_n, 'F', 1. )
-          ENDIF
-      ENDIF
+               CALL lbc_lnk( sshf_n, 'F', 1. )
+            ENDIF
+         ENDIF
+      ENDIF ! .NOT. ln_rstart
       !
 !!?      IF( .NOT. ln_rstart ) THEN           ! delete the initial ssh below sea-ice area
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4155)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4220)
@@ -22,4 +22,5 @@
    USE phycst         ! physical constants
    USE dom_oce        ! ocean space and time domain variables
+   USE oce     , ONLY :  iatte, oatte
    USE ice            ! LIM: sea-ice variables
    USE par_ice        ! LIM: sea-ice parameters
@@ -220,5 +221,5 @@
             !
             ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx)
-            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) )
+            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) )
             !
             ! oceanic heat flux (limthd_dh)
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4155)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4220)
@@ -656,5 +656,5 @@
          DO ji = kideb , kiut
             ! snow temperatures      
-            IF (ht_s_b(ji).GT.0) &
+            IF (ht_s_b(ji).GT.0._wp) &
                t_s_b(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
                *  t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) &
@@ -692,5 +692,5 @@
             DO ji = kideb , kiut
                ztmelt_i        = -tmut * s_i_b(ji,layer) + rtt 
-               t_i_b(ji,layer) =  MAX(MIN(t_i_b(ji,layer),ztmelt_i),190.0)
+               t_i_b(ji,layer) =  MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp)
                zerrit(ji)      =  MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer)))
             END DO
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 4155)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 4220)
@@ -74,5 +74,5 @@
       REAL(wp) ::   zindb  , zindsn , zindic, zindh, zinda      ! local scalar
       REAL(wp) ::   zusvosn, zusvoic, zbigval     !   -      -
-      REAL(wp) ::   zcfl , zusnit , zrtt          !   -      -
+      REAL(wp) ::   zcfl , zusnit                 !   -      -
       REAL(wp) ::   ze   , zsal   , zage          !   -      -
       !
@@ -450,5 +450,4 @@
                   zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi16 )
                   zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi16 )
-                  zrtt            = 173.15 * rone 
                   zindsn          = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) )
                   zindic          = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 4155)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 4220)
@@ -180,5 +180,4 @@
                ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb
                o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb
-               a_i(ji,jj,jl) = a_i (ji,jj,jl) * zindb ! clem correction
             END DO
          END DO
