Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5078)
@@ -165,7 +165,4 @@
 
    !                                     !!** ice-thickness distribution namelist (namiceitd) **
-   INTEGER , PUBLIC ::   jpl              !: number of ice  categories 
-   INTEGER , PUBLIC ::   nlay_i           !: number of ice  layers 
-   INTEGER , PUBLIC ::   nlay_s           !: number of snow layers 
    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)
@@ -223,8 +220,10 @@
    REAL(wp), PUBLIC ::   usecc2           !:  = 1.0 / ( rn_ecc * rn_ecc )
    REAL(wp), PUBLIC ::   rhoco            !: = rau0 * cio
-
+   REAL(wp), PUBLIC ::   r1_nlay_i        !: 1 / nlay_i
+   REAL(wp), PUBLIC ::   r1_nlay_s        !: 1 / nlay_s 
+   !
    !                                     !!** switch for presence of ice or not 
    REAL(wp), PUBLIC ::   rswitch
-
+   !
    !                                     !!** define some parameters 
    REAL(wp), PUBLIC, PARAMETER ::   epsi06   = 1.e-06_wp  !: small number 
@@ -368,9 +367,12 @@
    !!--------------------------------------------------------------------------
    !                                                  !!: ** Namelist namicerun read in sbc_lim_init **
-   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input)
-   CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output)
-   LOGICAL               , PUBLIC ::   ln_limdyn       !: flag for ice dynamics (T) or not (F)
-   LOGICAL               , PUBLIC ::   ln_nicep        !: flag for sea-ice points output (T) or not (F)
-   REAL(wp)              , PUBLIC ::   rn_amax         !: maximum ice concentration
+   INTEGER          , PUBLIC ::   jpl             !: number of ice  categories 
+   INTEGER          , PUBLIC ::   nlay_i          !: number of ice  layers 
+   INTEGER          , PUBLIC ::   nlay_s          !: number of snow layers 
+   CHARACTER(len=32), PUBLIC ::   cn_icerst_in    !: suffix of ice restart name (input)
+   CHARACTER(len=32), PUBLIC ::   cn_icerst_out   !: suffix of ice restart name (output)
+   LOGICAL          , PUBLIC ::   ln_limdyn       !: flag for ice dynamics (T) or not (F)
+   LOGICAL          , PUBLIC ::   ln_nicep        !: flag for sea-ice points output (T) or not (F)
+   REAL(wp)         , PUBLIC ::   rn_amax         !: maximum ice concentration
    !
    !!--------------------------------------------------------------------------
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 5078)
@@ -210,5 +210,5 @@
             &                    * e12t(:,:) * tmask(:,:,1) ) - zvi_b ) * r1_rdtice - zfw 
 
-         zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic )
+         zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) ) - zsmv_b ) * r1_rdtice + ( zfs * r1_rhoic )
 
          zei  =   glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  &
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 5078)
@@ -326,5 +326,5 @@
                ! recompute ht_i, ht_s avoiding out of bounds values
                ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
-               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn )
+               ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic * r1_rhosn )
 
                ! ice volume, salt content, age content
@@ -347,5 +347,5 @@
 
                    ! Mutliply by volume, and divide by number of layers to get heat content in J/m2
-                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) / nlay_s
+                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
                END DO ! ji
             END DO ! jj
@@ -368,5 +368,5 @@
 
                    ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
-                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) / nlay_i
+                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i
                END DO ! ji
             END DO ! jj
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 5078)
@@ -1124,5 +1124,5 @@
                ! clem: if sst>0, then ersw <0 (is that possible?)
                zsstK  = sst_m(ji,jj) + rt0
-               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) / REAL( nlay_i )
+               ersw(ji,jj,jk)   = - rhoic * vsw(ji,jj) * rcp * ( zsstK - rt0 ) * r1_nlay_i
 
                ! heat flux to the ocean
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5078)
@@ -337,5 +337,5 @@
       DO jl = 1, jpl
          DO jk = 1, nlay_i
-            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) / REAL( nlay_i )
+            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * a_i(:,:,jl) * ht_i(:,:,jl) * r1_nlay_i
          END DO
       END DO
@@ -347,5 +347,5 @@
       DO jl = 1, jpl
          DO jk = 1, nlay_s
-            e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) / REAL( nlay_s )
+            e_s(:,:,jk,jl) = e_s(:,:,jk,jl) * a_i(:,:,jl) * ht_s(:,:,jl) * r1_nlay_s
          END DO
       END DO
@@ -479,5 +479,5 @@
             ! Conversion q(S,T) -> T (second order equation)
             zaaa          =  cpic
-            zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) / rhoic - lfus
+            zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_1d(ji,jk) * r1_rhoic - lfus
             zccc          =  lfus * ( ztmelts - rtt )
             zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) )
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5078)
@@ -145,5 +145,5 @@
       DO jk = 1, nlay_i
          DO ji = kideb, kiut
-            h_i_old (ji,jk) = ht_i_1d(ji) / REAL( nlay_i )
+            h_i_old (ji,jk) = ht_i_1d(ji) * r1_nlay_i
             qh_i_old(ji,jk) = q_i_1d(ji,jk) * h_i_old(ji,jk)
          ENDDO
@@ -188,5 +188,5 @@
       !
       DO ji = kideb, kiut     
-         zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )
+         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s
       END DO
       !
@@ -199,5 +199,5 @@
       DO jk = 1, nlay_i
          DO ji = kideb, kiut
-            zh_i(ji,jk) = ht_i_1d(ji) / REAL( nlay_i )
+            zh_i(ji,jk) = ht_i_1d(ji) * r1_nlay_i
             zqh_i(ji)   = zqh_i(ji) + q_i_1d(ji,jk) * zh_i(ji,jk)
          END DO
@@ -228,5 +228,5 @@
          ! thickness change
          zcoeff = ( 1._wp - ( 1._wp - at_i_1d(ji) )**rn_betas ) / at_i_1d(ji) 
-         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
+         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice * r1_rhosn
          ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K)
          zqprec   (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )   
@@ -255,5 +255,5 @@
          zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) )      
          ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdh_s_mel(ji) )
-         zh_s  (ji) = ht_s_1d(ji) / REAL( nlay_s )
+         zh_s  (ji) = ht_s_1d(ji) * r1_nlay_s
 
          ENDIF
@@ -311,5 +311,5 @@
       DO ji = kideb, kiut
          dh_s_tot(ji)   = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
-         zh_s(ji)       = ht_s_1d(ji) / REAL( nlay_s )
+         zh_s(ji)       = ht_s_1d(ji) * r1_nlay_s
       END DO ! ji
 
@@ -335,5 +335,5 @@
       DO jk = 1, nlay_i
          DO ji = kideb, kiut 
-            zEi            = - q_i_1d(ji,jk) / rhoic                ! Specific enthalpy of layer k [J/kg, <0]
+            zEi            = - q_i_1d(ji,jk) * r1_rhoic             ! Specific enthalpy of layer k [J/kg, <0]
 
             ztmelts        = - tmut * s_i_1d(ji,jk) + rtt           ! Melting point of layer k [K]
@@ -345,5 +345,5 @@
             zfmdt          = - zq_su(ji) / zdE                     ! Mass flux to the ocean [kg/m2, >0]
 
-            zdeltah(ji,jk) = - zfmdt / rhoic                       ! Melt of layer jk [m, <0]
+            zdeltah(ji,jk) = - zfmdt * r1_rhoic                    ! Melt of layer jk [m, <0]
 
             zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) )    ! Melt of layer jk cannot exceed the layer thickness [m, <0]
@@ -506,5 +506,5 @@
                IF( t_i_1d(ji,jk) >= ztmelts ) THEN !!! Internal melting
 
-                  zEi               = - q_i_1d(ji,jk) / rhoic        ! Specific enthalpy of melting ice (J/kg, <0)
+                  zEi               = - q_i_1d(ji,jk) * r1_rhoic    ! Specific enthalpy of melting ice (J/kg, <0)
 
                   !!zEw               = rcp * ( t_i_1d(ji,jk) - rtt )  ! Specific enthalpy of meltwater at T = t_i_1d (J/kg, <0)
@@ -535,13 +535,13 @@
                ELSE                               !!! Basal melting
 
-                  zEi               = - q_i_1d(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
-
-                  zEw               = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0)
-
-                  zdE               = zEi - zEw              ! Specific enthalpy difference   (J/kg, <0)
-
-                  zfmdt             = - zq_bo(ji) / zdE  ! Mass flux x time step (kg/m2, >0)
-
-                  zdeltah(ji,jk)    = - zfmdt / rhoic        ! Gross thickness change
+                  zEi               = - q_i_1d(ji,jk) * r1_rhoic ! Specific enthalpy of melting ice (J/kg, <0)
+
+                  zEw               = rcp * ( ztmelts - rtt )    ! Specific enthalpy of meltwater (J/kg, <0)
+
+                  zdE               = zEi - zEw                  ! Specific enthalpy difference   (J/kg, <0)
+
+                  zfmdt             = - zq_bo(ji) / zdE          ! Mass flux x time step (kg/m2, >0)
+
+                  zdeltah(ji,jk)    = - zfmdt * r1_rhoic         ! Gross thickness change
 
                   zdeltah(ji,jk)    = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
@@ -627,5 +627,5 @@
          ! Salinity of snow ice
          ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
-         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji)
+         zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) * r1_rhoic + ( 1. - zswitch_sal ) * sm_i_1d(ji)
 
          ! entrapment during snow ice formation
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5078)
@@ -179,6 +179,6 @@
       zdq(:) = 0._wp ; zq_ini(:) = 0._wp      
       DO ji = kideb, kiut
-         zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  &
-            &           SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) ) 
+         zq_ini(ji) = ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
+            &           SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s ) 
       END DO
 
@@ -189,6 +189,6 @@
          isnow(ji)= 1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_1d(ji) ) )  ! is there snow or not
          ! layer thickness
-         zh_i(ji) = ht_i_1d(ji) / REAL( nlay_i )
-         zh_s(ji) = ht_s_1d(ji) / REAL( nlay_s )
+         zh_i(ji) = ht_i_1d(ji) * r1_nlay_i
+         zh_s(ji) = ht_s_1d(ji) * r1_nlay_s
       END DO
 
@@ -202,5 +202,5 @@
       DO jk = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer
          DO ji = kideb , kiut
-            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) / REAL( nlay_s )
+            z_s(ji,jk) = z_s(ji,jk-1) + ht_s_1d(ji) * r1_nlay_s
          END DO
       END DO
@@ -208,5 +208,5 @@
       DO jk = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer
          DO ji = kideb , kiut
-            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) / REAL( nlay_i )
+            z_i(ji,jk) = z_i(ji,jk-1) + ht_i_1d(ji) * r1_nlay_i
          END DO
       END DO
@@ -576,5 +576,5 @@
                !------------------------------------------------------------------------------|
                !
-               IF (t_su_1d(ji) < rtt) THEN
+               IF ( t_su_1d(ji) < rtt ) THEN
                   !
                   !------------------------------------------------------------------------------|
@@ -600,5 +600,5 @@
                   !!case of only one layer in the ice (surface & ice equations are altered)
 
-                  IF (nlay_i.eq.1) THEN
+                  IF ( nlay_i == 1 ) THEN
                      ztrid(ji,numeqmin(ji),1)    =  0.0
                      ztrid(ji,numeqmin(ji),2)    =  dzf(ji) - zkappa_i(ji,0) * 2.0
@@ -631,5 +631,5 @@
 
                   !!case of only one layer in the ice (surface & ice equations are altered)
-                  IF (nlay_i.eq.1) THEN
+                  IF ( nlay_i == 1 ) THEN
                      ztrid(ji,numeqmin(ji),1)  =  0.0
                      ztrid(ji,numeqmin(ji),2)  =  1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) )
@@ -760,6 +760,6 @@
       ! --- diag conservation imbalance on heat diffusion - PART 2 --- !
       DO ji = kideb, kiut
-         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) / REAL( nlay_i ) +  &
-            &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) / REAL( nlay_s ) )
+         zdq(ji)        = - zq_ini(ji) + ( SUM( q_i_1d(ji,1:nlay_i) ) * ht_i_1d(ji) * r1_nlay_i +  &
+            &                              SUM( q_s_1d(ji,1:nlay_s) ) * ht_s_1d(ji) * r1_nlay_s )
          IF( t_su_1d(ji) < rtt ) THEN  ! case T_su < 0degC
             zhfx_err(ji) = qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq(ji) * r1_rdtice 
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 5078)
@@ -102,5 +102,5 @@
       ! new layer thickesses
       DO ji = kideb, kiut
-         zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) / REAL( nlay_i )  
+         zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) * r1_nlay_i  
       ENDDO
 
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 5078)
@@ -362,7 +362,7 @@
          DO ji = 1, nbpac
 
-            zEi           = - ze_newice(ji) / rhoic                ! specific enthalpy of forming ice [J/kg]
-
-            zEw           = rcp * ( t_bo_1d(ji) - rt0 )             ! specific enthalpy of seawater at t_bo_1d [J/kg]
+            zEi           = - ze_newice(ji) * r1_rhoic             ! specific enthalpy of forming ice [J/kg]
+
+            zEw           = rcp * ( t_bo_1d(ji) - rt0 )            ! specific enthalpy of seawater at t_bo_1d [J/kg]
                                                                    ! clem: we suppose we are already at the freezing point (condition qlead<0 is satisfyied) 
                                                                    
@@ -371,5 +371,5 @@
             zfmdt         = - qlead_1d(ji) / zdE                   ! Fm.dt [kg/m2] (<0) 
                                                                    ! clem: we use qlead instead of zqld (limthd) because we suppose we are at the freezing point   
-            zv_newice(ji) = - zfmdt / rhoic
+            zv_newice(ji) = - zfmdt * r1_rhoic
 
             zQm           = zfmdt * zEw                            ! heat to the ocean >0 associated with mass flux  
@@ -458,5 +458,5 @@
             DO jk = 1, nlay_i
                DO ji = 1, nbpac
-                  h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i )
+                  h_i_old (ji,jk) = zv_i_1d(ji,jl) * r1_nlay_i
                   qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk)
                END DO
@@ -525,5 +525,5 @@
                DO ji = 1, jpi
                   ! heat content in J/m2
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) / REAL( nlay_i ,wp ) 
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * v_i(ji,jj,jl) * r1_nlay_i 
                END DO
             END DO
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 5078)
@@ -196,5 +196,5 @@
                   !
                   zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation)
-                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rtt ) + zq_i / rhoic - lfus
+                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rtt ) + zq_i * r1_rhoic - lfus
                   zccc       =  lfus * (ztmelts-rtt)
                   zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) )
@@ -280,6 +280,6 @@
       !!------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
-      REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac, zsal
-      REAL(wp) ::   zswi0, zswi01, zswibal, zargtemp , zs_zero   
+      REAL(wp) ::   zfac0, zfac1, zsal
+      REAL(wp) ::   zswi0, zswi01, zargtemp , zs_zero   
       REAL(wp), POINTER, DIMENSION(:,:,:) ::   z_slope_s, zalpha
       REAL(wp), PARAMETER :: zsi0 = 3.5_wp
@@ -311,6 +311,6 @@
          END DO
          !
-         dummy_fac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf
-         dummy_fac1 = zsi1 / ( zsi1 - zsi0 )
+         zfac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf
+         zfac1 = zsi1  / ( zsi1 - zsi0 )
          !
          zalpha(:,:,:) = 0._wp
@@ -322,14 +322,14 @@
                   ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 
                   zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp   , SIGN( 1._wp  , zsi1 - sm_i(ji,jj,jl) ) ) 
-                  ! If 2.sm_i GE sss_m then zswibal = 1
+                  ! If 2.sm_i GE sss_m then rswitch = 1
                   ! this is to force a constant salinity profile in the Baltic Sea
-                  zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
-                  zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * dummy_fac0 + dummy_fac1 )
-                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - zswibal )
-               END DO
-            END DO
-         END DO
-
-         dummy_fac = 1._wp / REAL( nlay_i )                   ! Computation of the profile
+                  rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
+                  zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 )
+                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch )
+               END DO
+            END DO
+         END DO
+
+         ! Computation of the profile
          DO jl = 1, jpl
             DO jk = 1, nlay_i
@@ -337,5 +337,5 @@
                   DO ji = 1, jpi
                      !                                      ! linear profile with 0 at the surface
-                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * dummy_fac
+                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
                      !                                      ! weighting the profile
                      s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
@@ -357,5 +357,5 @@
          DO jl = 1, jpl
             DO jk = 1, nlay_i
-               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)
+               zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
                zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
                s_i(:,:,jk,jl) =  zsal
@@ -387,5 +387,5 @@
                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
                   tm_i(ji,jj) = tm_i(ji,jj) + rswitch * t_i(ji,jj,jk,jl) * v_i(ji,jj,jl)   &
-                     &                      / (  REAL(nlay_i,wp) * MAX( vt_i(ji,jj) , epsi10 )  )
+                     &                      * r1_nlay_i / MAX( vt_i(ji,jj) , epsi10 )
                END DO
             END DO
@@ -417,5 +417,5 @@
                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rtt) + epsi10 ) )  )
                   zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rtt, - epsi10 )   &
-                     &                   * v_i(ji,jj,jl)    / REAL(nlay_i,wp)
+                     &                   * v_i(ji,jj,jl) * r1_nlay_i
                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi10 ) )  )
                   bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi10 )
@@ -439,6 +439,6 @@
       INTEGER  ::   ji, jk    ! dummy loop indices
       INTEGER  ::   ii, ij    ! local integers
-      REAL(wp) ::   dummy_fac0, dummy_fac1, dummy_fac2, zargtemp, zsal   ! local scalars
-      REAL(wp) ::   zalpha, zswi0, zswi01, zswibal, zs_zero              !   -      -
+      REAL(wp) ::   zfac0, zfac1, zargtemp, zsal   ! local scalars
+      REAL(wp) ::   zalpha, zswi0, zswi01, zs_zero              !   -      -
       !
       REAL(wp), POINTER, DIMENSION(:) ::   z_slope_s
@@ -466,8 +466,6 @@
          ! Weighting factor between zs_zero and zs_inf
          !---------------------------------------------
-         dummy_fac0 = 1._wp / ( zsi0 - zsi1 )
-         dummy_fac1 = zsi1 / ( zsi1 - zsi0 )
-         dummy_fac2 = 1._wp / REAL(nlay_i,wp)
-
+         zfac0 = 1._wp / ( zsi0 - zsi1 )
+         zfac1 = zsi1 / ( zsi1 - zsi0 )
          DO jk = 1, nlay_i
             DO ji = kideb, kiut
@@ -478,11 +476,11 @@
                ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 
                zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 
-               ! if 2.sm_i GE sss_m then zswibal = 1
+               ! if 2.sm_i GE sss_m then rswitch = 1
                ! this is to force a constant salinity profile in the Baltic Sea
-               zswibal = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
+               rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_m(ii,ij) ) )
                !
-               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * dummy_fac0 + dummy_fac1 )  ) * ( 1.0 - zswibal )
+               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 )  ) * ( 1.0 - rswitch )
                !
-               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * dummy_fac2
+               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
                ! weighting the profile
                s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
@@ -501,5 +499,5 @@
          !
          DO jk = 1, nlay_i
-            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) / REAL(nlay_i,wp)
+            zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
             zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
             DO ji = kideb, kiut
@@ -734,5 +732,5 @@
                ! recompute ht_i, ht_s avoiding out of bounds values
                zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
-               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic / rhosn )
+               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
             ENDIF
          ENDDO
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90	(revision 5078)
@@ -266,5 +266,5 @@
                      zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0* &
                         ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rtt ), - epsi06 ) ) * &
-                        rswitch / nlay_i
+                        rswitch * r1_nlay_i
                   END DO
                END DO
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/DOM/phycst.F90	(revision 5078)
@@ -82,4 +82,8 @@
    REAL(wp), PUBLIC ::   xlic     =  300.33e+6_wp    !: volumetric latent heat fusion of ice                  [J/m3]
    REAL(wp), PUBLIC ::   xsn      =    2.8e+6_wp     !: volumetric latent heat of sublimation of snow         [J/m3]
+#endif
+#if defined key_lim3
+   REAL(wp), PUBLIC ::   r1_rhoic                    !: 1 / rhoic
+   REAL(wp), PUBLIC ::   r1_rhosn                    !: 1 / rhosn
 #endif
    !!----------------------------------------------------------------------
@@ -166,5 +170,8 @@
       lfus = xlsn / rhosn        ! latent heat of fusion of fresh ice
 #endif
-
+#if defined key_lim3
+      r1_rhoic = 1._wp / rhoic
+      r1_rhosn = 1._wp / rhosn
+#endif
       IF(lwp) THEN
          WRITE(numout,*)
Index: /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5077)
+++ /branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5078)
@@ -391,4 +391,8 @@
       rdt_ice   = nn_fsbc * rdttra(1)  
       r1_rdtice = 1._wp / rdt_ice 
+
+      ! inverse of nlay_i and nlay_s
+      r1_nlay_i = 1._wp / REAL( nlay_i, wp )
+      r1_nlay_s = 1._wp / REAL( nlay_s, wp )
       !
    END SUBROUTINE ice_run
