Index: branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4146)
+++ branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 4161)
@@ -46,5 +46,5 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)
+   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
    !! $Id$
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -78,11 +78,10 @@
       !!               update ht_s_b, ht_i_b and tbif_1d(:,:)      
       !!------------------------------------------------------------------------
-      INTEGER  ::   ji,jj,jk,jl,jm   ! dummy loop indices
-      INTEGER  ::   layer, nbpac     ! local integers 
-      INTEGER  ::   zji, zjj, iter   !   -       -
-      REAL(wp) ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde   ! local scalars
+      INTEGER ::   ji,jj,jk,jl,jm   ! dummy loop indices
+      INTEGER ::   layer, nbpac     ! local integers 
+      INTEGER ::   ii, ij, iter   !   -       -
+      REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde  ! local scalars
       REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      -
       REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      -
-      REAL(wp) ::   zcoef                                                       !   -      -
       LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness
       CHARACTER (len = 15) :: fieldid
@@ -160,5 +159,5 @@
                DO ji = 1, jpi
                   !Energy of melting q(S,T) [J.m-3]
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * nlay_i
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * REAL( nlay_i )
                   zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) )  )   !0 if no ice and 1 if yes
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb
@@ -342,7 +341,7 @@
          CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)]
             DO ji = 1, nbpac
-               zji =   MOD( npac(ji) - 1 , jpi ) + 1
-               zjj =      ( npac(ji) - 1 ) / jpi + 1
-               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(zji,zjj)  )
+               ii =   MOD( npac(ji) - 1 , jpi ) + 1
+               ij =      ( npac(ji) - 1 ) / jpi + 1
+               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij)  )
             END DO
          CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
@@ -389,24 +388,15 @@
          END DO
 
-         !---------------------------------
-         ! Salt flux due to new ice growth
-         !---------------------------------
-         ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine)
-         DO ji = 1, nbpac
-            sfx_thd_1d(ji) = sfx_thd_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice
-            rdm_ice_1d(ji) = rdm_ice_1d(ji) +                 rhoic * zv_newice(ji)
-         END DO ! ji
-
          !------------------------------------
          ! Diags for energy conservation test
          !------------------------------------
          DO ji = 1, nbpac
-            zji = MOD( npac(ji) - 1 , jpi ) + 1
-            zjj =    ( npac(ji) - 1 ) / jpi + 1
+            ii = MOD( npac(ji) - 1 , jpi ) + 1
+            ij =    ( npac(ji) - 1 ) / jpi + 1
             !
-            zde = ze_newice(ji) / unit_fac * area(zji,zjj) * zv_newice(ji)
+            zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji)
             !
-            vt_i_init(zji,zjj) = vt_i_init(zji,zjj) + zv_newice(ji)             ! volume
-            et_i_init(zji,zjj) = et_i_init(zji,zjj) + zde                       ! Energy
+            vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji)             ! volume
+            et_i_init(ii,ij) = et_i_init(ii,ij) + zde                       ! Energy
 
          END DO
@@ -419,8 +409,8 @@
          !-----------------
          DO ji = 1, nbpac
-            zji = MOD( npac(ji) - 1 , jpi ) + 1
-            zjj =    ( npac(ji) - 1 ) / jpi + 1
+            ii = MOD( npac(ji) - 1 , jpi ) + 1
+            ij =    ( npac(ji) - 1 ) / jpi + 1
             za_newice(ji) = zv_newice(ji) / zh_newice(ji)
-            diag_lat_gr(zji,zjj) = zv_newice(ji) * r1_rdtice
+            diag_lat_gr(ii,ij) = diag_lat_gr(ii,ij) + zv_newice(ji) * r1_rdtice ! clem
          END DO !ji
 
@@ -441,6 +431,6 @@
          ! we keep the excessive volume in memory and attribute it later to bottom accretion
          DO ji = 1, nbpac
-            IF ( za_newice(ji)  >  ( 1._wp - zat_i_ac(ji) ) ) THEN
-               zda_res(ji)   = za_newice(ji) - (1.0 - zat_i_ac(ji) )
+            IF ( za_newice(ji) >  ( amax - zat_i_ac(ji) ) ) THEN
+               zda_res(ji)   = za_newice(ji) - ( amax - zat_i_ac(ji) )
                zdv_res(ji)   = zda_res  (ji) * zh_newice(ji) 
                za_newice(ji) = za_newice(ji) - zda_res  (ji)
@@ -473,8 +463,8 @@
          DO ji = 1, nbpac
             jl = zcatac(ji)                                                           ! categroy in which new ice is put
-            zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) ) )             ! zindb=1 if ice =0 otherwise
+            zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , -za_old(ji,jl) + epsi10 ) )             ! zindb=1 if ice =0 otherwise
             zhice_old(ji,jl) = zv_old(ji,jl) / MAX( za_old(ji,jl) , epsi10 ) * zindb  ! old ice thickness
             zdhex    (ji) = MAX( 0._wp , zh_newice(ji) - zhice_old(ji,jl) )           ! difference in thickness
-            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi11 ) )   ! ice totally new in jl category
+            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) )   ! ice totally new in jl category
          END DO
 
@@ -482,11 +472,11 @@
             DO ji = 1, nbpac
                jl = zcatac(ji)
-               zqold   = ze_i_ac(ji,jk,jl)      ! [ J.m-3 ]
-               zalphai = MIN( zhice_old(ji,jl) *   jk       / nlay_i , zh_newice(ji) )   &
-                  &    - MIN( zhice_old(ji,jl) * ( jk - 1 ) / nlay_i , zh_newice(ji) )
+               zqold   = ze_i_ac(ji,jk,jl) ! [ J.m-3 ]
+               zalphai = MIN( zhice_old(ji,jl) * REAL( jk )     / REAL( nlay_i ), zh_newice(ji) )   &
+                  &    - MIN( zhice_old(ji,jl) * REAL( jk - 1 ) / REAL( nlay_i ), zh_newice(ji) )
                ze_i_ac(ji,jk,jl) = zswinew(ji) * ze_newice(ji)                                     &
-                  + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl)  * zqold * zhice_old(ji,jl) / nlay_i   &
+                  + ( 1.0 - zswinew(ji) ) * ( za_old(ji,jl)  * zqold * zhice_old(ji,jl) / REAL( nlay_i )  &
                   + za_newice(ji)  * ze_newice(ji) * zalphai                                       &
-                  + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / nlay_i ) / ( zv_i_ac(ji,jl) / nlay_i )
+                  + za_newice(ji)  * ze_newice(ji) * zdhex(ji) / REAL( nlay_i ) ) / ( ( zv_i_ac(ji,jl) ) / REAL( nlay_i ) )
             END DO
          END DO
@@ -513,5 +503,6 @@
             DO ji = 1, nbpac
                zindb = MAX( 0._wp, SIGN( 1._wp , zdv_res(ji) ) )
-               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 )
+               zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_lev(ji) - epsi06 ) )  ! clem
+               zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zindb * zinda * zdv_res(ji) * za_i_ac(ji,jl) / MAX( zat_i_lev(ji) , epsi06 )
             END DO
          END DO
@@ -524,5 +515,5 @@
          DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
             DO ji = 1, nbpac
-               zindb =  1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) ) )       ! zindb=1 if ice =0 otherwise
+               zindb =  1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl ) + epsi10 ) )       ! zindb=1 if ice =0 otherwise
                zhice_old(ji,jl) = zv_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb
                zdhicbot (ji,jl) = zdv_res(ji)    / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb    &
@@ -536,5 +527,5 @@
             DO jk = 1, nlay_i
                DO ji = 1, nbpac
-                  zthick0(ji,jk,jl) =  zhice_old(ji,jl) / nlay_i
+                  zthick0(ji,jk,jl) =  zhice_old(ji,jl) / REAL( nlay_i )
                   zqm0   (ji,jk,jl) =  ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl)
                END DO
@@ -555,9 +546,9 @@
                DO layer = 1, nlay_i + 1
                   DO ji = 1, nbpac
-                     zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) ) 
+                     zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) ) 
                      ! Redistributing energy on the new grid
-                     zweight = MAX (  MIN( zhice_old(ji,jl) * layer , zdummy(ji,jl) * jk )   &
-                        &    - MAX( zhice_old(ji,jl) * ( layer - 1 ) , zdummy(ji,jl) * ( jk - 1 ) ) , 0._wp )   &
-                        &    /( MAX(nlay_i * zthick0(ji,layer,jl),epsi10) ) * zindb
+                     zweight = MAX (  MIN( zhice_old(ji,jl) * REAL( layer ), zdummy(ji,jl) * REAL( jk ) )   &
+                        &    - MAX( zhice_old(ji,jl) * REAL( layer - 1 ) , zdummy(ji,jl) * REAL( jk - 1 ) ) , 0._wp )   &
+                        &    /( MAX(REAL(nlay_i) * zthick0(ji,layer,jl),epsi10) ) * zindb
                      ze_i_ac(ji,jk,jl) =  ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl)  
                   END DO ! ji
@@ -569,7 +560,7 @@
             DO jk = 1, nlay_i
                DO ji = 1, nbpac
-                  zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) ) 
+                  zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) ) 
                   ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl)   &
-                     &              / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * nlay_i * zindb
+                     &              / MAX( zv_i_ac(ji,jl) , epsi10) * za_i_ac(ji,jl) * REAL( nlay_i ) * zindb
                END DO
             END DO
@@ -581,5 +572,5 @@
          DO jl = 1, jpl
             DO ji = 1, nbpac
-               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) ) )  ! 0 if no ice and 1 if yes
+               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
                zoa_i_ac(ji,jl)  = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb   
             END DO 
@@ -589,13 +580,25 @@
          ! Update salinity
          !-----------------
-         IF(  num_sal == 2  ) THEN      ! Sice = F(z,t)
+         !clem IF(  num_sal == 2  ) THEN
             DO jl = 1, jpl
                DO ji = 1, nbpac
-                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) ) )   ! 0 if no ice and 1 if yes
+                  zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
                   zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
-                  zsmv_i_ac(ji,jl) = ( zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) ) * zindb
+                  zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) * zindb ! clem modif
                END DO
             END DO   
-         ENDIF
+         !clem ENDIF
+
+         !--------------------------------
+         ! Update mass/salt fluxes (clem)
+         !--------------------------------
+         DO jl = 1, jpl
+            DO ji = 1, nbpac
+               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes
+               zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl)
+               rdm_ice_1d(ji) = rdm_ice_1d(ji) + zdv * rhoic * zindb
+               sfx_thd_1d(ji)   =   sfx_thd_1d(ji) - zdv * rhoic * zs_newice(ji) * r1_rdtice * zindb
+           END DO
+         END DO
 
          !------------------------------------------------------------------------------!
@@ -606,5 +609,5 @@
             CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj )
             CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj )
-            IF (  num_sal == 2  )   &
+            !clem IF (  num_sal == 2  )   &
                CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj )
             DO jk = 1, nlay_i
@@ -622,5 +625,5 @@
       DO jl = 1, jpl
          DO jk = 1, nlay_i          ! heat content in 10^9 Joules
-            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / nlay_i / unit_fac 
+            e_i(:,:,jk,jl) = e_i(:,:,jk,jl) * area(:,:) * v_i(:,:,jl) / REAL( nlay_i )  / unit_fac 
          END DO
       END DO
