Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8496)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8498)
@@ -157,8 +157,8 @@
       !!------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   ze_i, z1_CpR, zdiscrim, zaaa, z1_2aaa             ! local scalars
-      REAL(wp) ::   ze_s, zL_Cp , ztmelts , zbbb, zccc                !   -      -
-      REAL(wp) ::   zhmax, z1_zhmax, zsm_i, zcpMcpic, zt_i   !   -      -
-      REAL(wp) ::   zlay_i, zlay_s   !   -      -
+      REAL(wp) ::   ze_i, z1_cp, z1_2cp             ! local scalars
+      REAL(wp) ::   ze_s, ztmelts, zbbb, zccc       !   -      -
+      REAL(wp) ::   zhmax, z1_zhmax, zsm_i          !   -      -
+      REAL(wp) ::   zlay_i, zlay_s                  !   -      -
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
       !!------------------------------------------------------------------
@@ -181,5 +181,5 @@
       ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)    !--- ice thickness
 
-	  zhmax    =          hi_max(jpl)
+      zhmax    =          hi_max(jpl)
       z1_zhmax =  1._wp / hi_max(jpl)               
       WHERE( ht_i(:,:,jpl) > zhmax )               !--- bound ht_i by hi_max (i.e. 99 m) with associated update of ice area
@@ -205,7 +205,5 @@
       !-------------------
       zlay_i   = REAL( nlay_i , wp )    ! number of layers
-      zaaa     = cpic                   ! Conversion q(S,T) -> T (second order equation)
-      z1_2aaa  = 1._wp / ( 2._wp * zaaa )
-      zcpMcpic = rcp - cpic
+      z1_2cp  = 1._wp / ( 2._wp * cpic )
       DO jl = 1, jpl
          DO jk = 1, nlay_i
@@ -214,15 +212,13 @@
                   IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area 
                      !
-                     ze_i    =   e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3]
-                     ztmelts = - s_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
-                     !
-                     zbbb     =  zcpMcpic * ztmelts + ze_i * r1_rhoic - lfus
-                     zccc     =  lfus * ztmelts
-                     zdiscrim =  SQRT(  MAX( zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp)  )
-                     zt_i     = - ( zbbb + zdiscrim ) * z1_2aaa                          ! [C]
-                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( zt_i , ztmelts )  ) + rt0   ! [K] with bounds: -100 < zt_i < ztmelts
+                     ze_i             =   e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3]
+                     ztmelts          = - s_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
+                     ! Conversion q(S,T) -> T (second order equation)
+                     zbbb             = ( rcp - cpic ) * ztmelts + ze_i * r1_rhoic - lfus
+                     zccc             = SQRT( MAX( zbbb * zbbb - 4._wp * cpic * lfus * ztmelts , 0._wp) )
+                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( -( zbbb + zccc ) * z1_2cp , ztmelts ) ) + rt0   ! [K] with bounds: -100 < t_i < ztmelts
                      !
                   ELSE                                !--- no ice
-                     t_i(ji,jj,jk,jl) =  rt0 - 100._wp                                   ! impose 173.15 K (i.e. -100 C)
+                     t_i(ji,jj,jk,jl) = rt0 - 100._wp                                   ! impose 173.15 K (i.e. -100 C)
 !!clem: I think we should impose rt0 instead
                   ENDIF
@@ -235,25 +231,13 @@
       ! Snow temperature   [K]   (with a minimum value (rt0 - 100.) imposed everywhere)
       !--------------------
-      z1_CpR = 1._wp / ( cpic * rhosn )
-      zL_Cp  = lfus  /   cpic
       zlay_s = REAL( nlay_s , wp )
+      z1_cp  = 1._wp / cpic
       DO jk = 1, nlay_s
          WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
-            t_s(:,:,jk,:) = MAX(  -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp )  ) + rt0       
+            t_s(:,:,jk,:) = MAX( -100._wp , MIN( z1_cp * ( -r1_rhosn * (e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s) + lfus ) , 0._wp ) ) + rt0
          ELSEWHERE                           !--- no ice
             t_s(:,:,jk,:) =  rt0 - 100._wp         ! impose 173.15 K (i.e. -100 C)
          END WHERE
       END DO
-!!gm perhaps better like this (?) :  (if this compile, yes! one test instead of nlay_s tests)
-!      WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
-!         DO jk = 1, nlay_s
-!            t_s(:,:,jk,:) = MAX(  -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp )  ) + rt0       
-!         END DO
-!      ELSEWHERE                           !--- no ice
-!         DO jk = 1, nlay_s
-!            t_s(:,:,jk,:) =  rt0 - 100._wp      ! impose 173.15 K (i.e. -100 C)
-!         END DO
-!      END WHERE
-!!gm end better ?
 
       ! integrated values 
