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 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 5070)
@@ -147,6 +147,6 @@
    !! tm_i        |      -      |    Mean sea ice temperature     | K     |
    !! ot_i        !      -      !    Sea ice areal age content    | day   |
-   !! et_i        !      -      !    Total ice enthalpy           | 10^9 J| 
-   !! et_s        !      -      !    Total snow enthalpy          | 10^9 J| 
+   !! et_i        !      -      !    Total ice enthalpy           | J/m2  | 
+   !! et_s        !      -      !    Total snow enthalpy          | J/m2  | 
    !! bv_i        !      -      !    Mean relative brine volume   | ???   | 
    !!=====================================================================
@@ -172,7 +172,9 @@
 
    !                                     !!** ice-dynamics namelist (namicedyn) **
+   LOGICAL , PUBLIC ::   ln_icestr_bvf    !: use brine volume to diminish ice strength
    INTEGER , PUBLIC ::   nn_icestr        !: ice strength parameterization (0=Hibler79 1=Rothrock75)
-   INTEGER , PUBLIC ::   nn_icestr_bvf    !: use brine volume to diminish ice strength
    INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling
+   INTEGER , PUBLIC ::   nn_ahi0          !: sea-ice hor. eddy diffusivity coeff. (3 ways of calculation)
+   REAL(wp), PUBLIC ::   rn_pe_rdg        !: ridging work divided by pot. energy change in ridging, nn_icestr = 1
    REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress
    REAL(wp), PUBLIC ::   rn_pstar         !: determines ice strength (N/M), Hibler JPO79
@@ -202,5 +204,4 @@
    !                                     !!** ice-mechanical redistribution namelist (namiceitdme)
    REAL(wp), PUBLIC ::   rn_cs            !: fraction of shearing energy contributing to ridging            
-   REAL(wp), PUBLIC ::   rn_pe_rdg        !: ridging work divided by pot. energy change in ridging, nn_ice str = 1
    REAL(wp), PUBLIC ::   rn_fsnowrdg      !: fractional snow loss to the ocean during ridging
    REAL(wp), PUBLIC ::   rn_fsnowrft      !: fractional snow loss to the ocean during ridging
@@ -217,5 +218,5 @@
 
    !                                     !!** ice-mechanical redistribution namelist (namiceitdme)
-   INTEGER , PUBLIC ::   nn_rafting      !: rafting of ice or not                        
+   LOGICAL , PUBLIC ::   ln_rafting      !: rafting of ice or not                        
    INTEGER , PUBLIC ::   nn_partfun      !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007)
 
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 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90	(revision 5070)
@@ -173,45 +173,55 @@
       IF( icount == 0 ) THEN
 
-         zvi_b  = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tmask(:,:,1) )
+         zfs_b  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
+            &                  sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &
+            &                ) *  e12t(:,:) * tmask(:,:,1) )
+
+         zfw_b  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
+            &                  wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    &
+            &                ) *  e12t(:,:) * tmask(:,:,1) )
+
+         zft_b  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
+            &                - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
+            &                ) *  e12t(:,:) * tmask(:,:,1) * zconv )
+
+         zvi_b  = glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * e12t(:,:) * tmask(:,:,1) )
+
          zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) )
-         zei_b  = glob_sum( SUM( SUM(   e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv +  &
-            &               SUM( SUM(   e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv )
-         zfw_b  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
-            &                   wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    &
-            &             ) * e12t(:,:) * tmask(:,:,1) )
-         zfs_b  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                   sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  &
-            &                 ) * e12t(:,:) * tmask(:,:,1) )
-         zft_b  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
-            &                 - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
-            &                  ) * e12t(:,:) * tmask(:,:,1) * zconv )
+
+         zei_b  = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  &
+            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    &
+                            ) * e12t(:,:) * tmask(:,:,1) * zconv )
 
       ELSEIF( icount == 1 ) THEN
 
-         zfs  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
-            &                 sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
-            &                ) * e12t(:,:) * tmask(:,:,1) ) - zfs_b
-         zfw  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
-            &                 wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    &
-            &                ) * e12t(:,:) * tmask(:,:,1) ) - zfw_b
-         zft  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
-            &               - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
-            &                ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b
+         zfs  = glob_sum(  ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  &
+            &                sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:)                                  & 
+            &              ) * e12t(:,:) * tmask(:,:,1) ) - zfs_b
+
+         zfw  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  &
+            &                wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:)    &
+            &              ) * e12t(:,:) * tmask(:,:,1) ) - zfw_b
+
+         zft  = glob_sum(  ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  & 
+            &              - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:)   &
+            &              ) * e12t(:,:) * tmask(:,:,1) * zconv ) - zft_b
  
-         zvi  = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 )  &
+         zvi  = ( glob_sum( SUM( v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 )  &
             &                    * 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 )
-         zei  =   glob_sum( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv +   &
-            &               SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) * e12t(:,:) * tmask(:,:,1) * zconv )  &
-            &                         * r1_rdtice - zei_b * r1_rdtice + zft
-
-         zvmin = glob_min(v_i)
-         zamax = glob_max(SUM(a_i,dim=3))
-         zamin = glob_min(a_i)
+
+         zei  =   glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) +  &
+            &                 SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )    &
+            &                ) * e12t(:,:) * tmask(:,:,1) * zconv ) * r1_rdtice - zei_b * r1_rdtice + zft
+
+         zvmin = glob_min( v_i )
+         zamax = glob_max( SUM( a_i, dim=3 ) )
+         zamin = glob_min( a_i )
        
          IF(lwp) THEN
             IF ( ABS( zvi    ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (',cd_routine,') = ',(zvi * rday)
             IF ( ABS( zsmv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday)
-            IF ( ABS( zei    ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',(zei)
+            IF ( ABS( zei    ) >  1     ) WRITE(numout,*) 'violation enthalpy [GW]       (',cd_routine,') = ',(zei)
             IF ( zvmin <  -epsi10       ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin)
             IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > rn_amax+epsi10 ) THEN
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90	(revision 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90	(revision 5070)
@@ -180,5 +180,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF( ABS( sfx (ji,jj) ) .GT. 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth
+            IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth
                !CALL lim_prt( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' )
                !DO jl = 1, jpl
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90	(revision 5070)
@@ -113,6 +113,6 @@
 
       ! Heat budget
-      zbg_ihc      = glob_sum( et_i(:,:) * 1.e-20 ) ! ice heat content  [1.e-20 J]
-      zbg_shc      = glob_sum( et_s(:,:) * 1.e-20 ) ! snow heat content [1.e-20 J]
+      zbg_ihc      = glob_sum( et_i(:,:) * e12t(:,:) * 1.e-20 ) ! ice heat content  [1.e20 J]
+      zbg_shc      = glob_sum( et_s(:,:) * e12t(:,:) * 1.e-20 ) ! snow heat content [1.e20 J]
       zbg_hfx_dhc  = glob_sum( diag_heat_dhc(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
       zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W]
@@ -137,5 +137,6 @@
       z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) ! salt fluxes
       z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + &
-                          &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes
+                          &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + &
+                          &     wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes
       !
       frc_vol  = frc_vol  + z_frc_vol  * rdt_ice
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 5070)
@@ -172,5 +172,5 @@
       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
          !
-         zcoef = SQRT( 0.5_wp ) / rau0
+         zcoef = SQRT( 0.5_wp ) * r1_rau0
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
@@ -243,5 +243,6 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicedyn/ nn_icestr, nn_icestr_bvf, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, nn_nevp, rn_relast, rn_ahi0_ref
+      NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio,  rn_creepl, rn_ecc, &
+         &                nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref
       INTEGER  ::   ji, jj
       REAL(wp) ::   za00, zd_max
@@ -261,14 +262,16 @@
          WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
          WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)  nn_icestr =', nn_icestr 
-         WRITE(numout,*)'    Including brine volume in ice strength comp.     nn_icestr_bvf=', nn_icestr_bvf
-         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,*) '   creep limit                                      rn_creepl = ', rn_creepl
-         WRITE(numout,*) '   eccentricity of the elliptical yield curve       rn_ecc    = ', rn_ecc
-         WRITE(numout,*) '   number of iterations for subcycling              nn_nevp   = ', nn_nevp
-         WRITE(numout,*) '   ratio of elastic timescale over ice time step    rn_relast = ', rn_relast
-         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)       rn_ahi0_ref   = ', rn_ahi0_ref
+         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,*) '   creep limit                                          rn_creepl     = ', rn_creepl
+         WRITE(numout,*) '   eccentricity of the elliptical yield curve           rn_ecc        = ', rn_ecc
+         WRITE(numout,*) '   number of iterations for subcycling                  nn_nevp       = ', nn_nevp
+         WRITE(numout,*) '   ratio of elastic timescale over ice time step        rn_relast     = ', rn_relast
+         WRITE(numout,*) '   horizontal diffusivity calculation                   nn_ahi0       = ', nn_ahi0
+         WRITE(numout,*) '   horizontal diffusivity coeff. (orca2 grid)           rn_ahi0_ref   = ', rn_ahi0_ref
       ENDIF
       !
@@ -277,20 +280,46 @@
       !
       !  Diffusion coefficients
-      zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
-      IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
-           
-      za00 = rn_ahi0_ref / ( 1.e05_wp )              ! 1.e05 = 100km = max grid space at 60° latitude in orca2
-                                              !                    (60° = min latitude for ice cover)  
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)
-            ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)
+      SELECT CASE( nn_ahi0 )
+
+      CASE( 0 )
+         ahiu(:,:) = rn_ahi0_ref
+         ahiv(:,:) = rn_ahi0_ref
+
+         IF(lwp) WRITE(numout,*) ''
+         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim constant = rn_ahi0_ref'
+
+      CASE( 1 ) 
+
+         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+         IF( lk_mpp )   CALL mpp_max( zd_max )          ! max over the global domain
+         
+         ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp   ! 1.e05 = 100km = max grid space at 60° latitude in orca2
+                                                        !                    (60° = min latitude for ice cover)  
+         ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp
+
+         IF(lwp) WRITE(numout,*) ''
+         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')'
+         IF(lwp) WRITE(numout,*) '   value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp 
+         
+      CASE( 2 ) 
+
+         zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )
+         IF( lk_mpp )   CALL mpp_max( zd_max )   ! max over the global domain
+         
+         za00 = rn_ahi0_ref * 1.e-05_wp          ! 1.e05 = 100km = max grid space at 60° latitude in orca2
+                                                 !                    (60° = min latitude for ice cover)  
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)
+               ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)
+            END DO
          END DO
-      END DO
-      !
-      IF(lwp) WRITE(numout,*) ''
-      IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1'
-      IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max
- 
+         !
+         IF(lwp) WRITE(numout,*) ''
+         IF(lwp) WRITE(numout,*) '   laplacian operator: ahim proportional to e1'
+         IF(lwp) WRITE(numout,*) '   maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max
+         
+      END SELECT
+
    END SUBROUTINE lim_dyn_init
 
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90	(revision 5070)
@@ -72,5 +72,5 @@
          DO jj = 2, jpjm1  
             DO ji = fs_2 , fs_jpim1   ! vector opt.
-               efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) / ( e1t(ji,jj) * e2t(ji,jj) )
+               efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) * r1_e12t(ji,jj)
             END DO
          END DO
@@ -98,6 +98,6 @@
          DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
             DO ji = 1 , fs_jpim1   ! vector opt.
-               zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
-               zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
+               zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
+               zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
             END DO
          END DO
@@ -105,6 +105,5 @@
          DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
             DO ji = fs_2 , fs_jpim1   ! vector opt. 
-               zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   &
-                  &            + zflv(ji,jj) - zflv(ji  ,jj-1)  ) / ( e1t (ji,jj) * e2t (ji,jj) )
+               zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
             END DO
          END DO
@@ -116,6 +115,6 @@
                zrlxint = (   ztab0(ji,jj)    &
                   &       +  rdt_ice * (           zalfa   * ( zdiv(ji,jj) + efact(ji,jj) * ptab(ji,jj) )   &
-                  &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )  )                             & 
-                  &    / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )
+                  &                      + ( 1.0 - zalfa ) *   zdiv0(ji,jj) )                               & 
+                  &      ) / ( 1.0 + zalfa * rdt_ice * efact(ji,jj) )
                zrlx(ji,jj) = ptab(ji,jj) + zrelax * ( zrlxint - ptab(ji,jj) )
             END DO
@@ -139,6 +138,6 @@
       DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction
          DO ji = 1 , fs_jpim1   ! vector opt.
-            zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
-            zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
+            zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) * r1_e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) )
+            zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) * r1_e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) )
          END DO
       END DO
@@ -146,6 +145,5 @@
       DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes
          DO ji = fs_2 , fs_jpim1   ! vector opt. 
-            zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   &
-                 &            + zflv(ji,jj) - zflv(ji  ,jj-1)  ) / ( e1t (ji,jj) * e2t (ji,jj) )
+            zdiv(ji,jj) = ( zflu(ji,jj) - zflu(ji-1,jj) + zflv(ji,jj) - zflv(ji,jj-1) ) * r1_e12t(ji,jj)
             ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) )
          END DO
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 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 5070)
@@ -157,5 +157,5 @@
       ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
       !-----------------------------------------------------------------------------!
-      Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0                ! proport const for PE
+      Cp = 0.5 * grav * (rau0-rhoic) * rhoic * r1_rau0                ! proport const for PE
       !
       CALL lim_itd_me_ridgeprep                                    ! prepare ridging
@@ -191,5 +191,5 @@
             !  (thick, newly ridged ice).
 
-            closing_net(ji,jj) = rn_cs * 0.5 * ( Delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
+            closing_net(ji,jj) = rn_cs * 0.5 * ( delta_i(ji,jj) - ABS( divu_i(ji,jj) ) ) - MIN( divu_i(ji,jj), 0._wp )
 
             ! 2.2 divu_adv
@@ -235,7 +235,7 @@
                ! Reduce the closing rate if more than 100% of the open water 
                ! would be removed.  Reduce the opening rate proportionately.
-               IF ( ato_i(ji,jj) .GT. epsi10 .AND. athorn(ji,jj,0) .GT. 0.0 ) THEN
+               IF ( ato_i(ji,jj) > epsi10 .AND. athorn(ji,jj,0) > 0.0 ) THEN
                   w1 = athorn(ji,jj,0) * closing_gross(ji,jj) * rdt_ice
-                  IF ( w1 .GT. ato_i(ji,jj)) THEN
+                  IF ( w1 > ato_i(ji,jj)) THEN
                      tmpfac = ato_i(ji,jj) / w1
                      closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac
@@ -286,5 +286,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF (ABS(asum(ji,jj) - kamax ) .LT. epsi10) THEN
+               IF (ABS(asum(ji,jj) - kamax ) < epsi10) THEN
                   closing_net(ji,jj) = 0._wp
                   opning     (ji,jj) = 0._wp
@@ -357,8 +357,7 @@
                   WRITE(numout,*) jl, a_i(ji,jj,jl), athorn(ji,jj,jl)
                END DO
-            ENDIF  ! asum
-
-         END DO !ji
-      END DO !jj
+            ENDIF
+         END DO
+      END DO
 
       ! Conservation check
@@ -482,9 +481,9 @@
                         * z1_3 * (hrmax(ji,jj,jl)**3 - hrmin(ji,jj,jl)**3) / ( hrmax(ji,jj,jl)-hrmin(ji,jj,jl) )   
 !!gm Optimization:  (a**3-b**3)/(a-b) = a*a+ab+b*b   ==> less costly operations even if a**3 is replaced by a*a*a...                    
-                  ENDIF            ! aicen > epsi10
+                  ENDIF
                   !
-               END DO ! ji
-            END DO !jj
-         END DO !jl
+               END DO
+            END DO
+         END DO
 
          zzc = rn_pe_rdg * Cp     ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation
@@ -510,9 +509,9 @@
       ! CAN BE REMOVED
       !
-      IF ( nn_icestr_bvf == 1 ) THEN
+      IF( ln_icestr_bvf ) THEN
 
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF ( bv_i(ji,jj) .GT. 0.0 ) THEN
+               IF ( bv_i(ji,jj) > 0.0 ) THEN
                   zdummy = MIN ( bv_i(ji,jj), 0.10 ) * MIN( bv_i(ji,jj), 0.10 )
                ELSE
@@ -520,6 +519,6 @@
                ENDIF
                strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0)))
-            END DO              ! j
-         END DO                 ! i
+            END DO
+         END DO
 
       ENDIF
@@ -539,6 +538,5 @@
          DO jj = 2, jpj - 1
             DO ji = 2, jpi - 1
-               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN ! ice is
-                  ! present
+               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN ! ice is present
                   zworka(ji,jj) = 4.0 * strength(ji,jj)              &
                      &          + strength(ji-1,jj) * tmask(ji-1,jj,1) &  
@@ -562,5 +560,5 @@
          CALL lbc_lnk( strength, 'T', 1. )
 
-      ENDIF ! ksmooth
+      ENDIF
 
       !--------------------
@@ -579,8 +577,8 @@
          DO jj = 1, jpj - 1
             DO ji = 1, jpi - 1
-               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) .GT. epsi10) THEN       ! ice is present
+               IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > epsi10) THEN       ! ice is present
                   numts_rm = 1 ! number of time steps for the running mean
-                  IF ( strp1(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
-                  IF ( strp2(ji,jj) .GT. 0.0 ) numts_rm = numts_rm + 1
+                  IF ( strp1(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
+                  IF ( strp2(ji,jj) > 0.0 ) numts_rm = numts_rm + 1
                   zp = ( strength(ji,jj) + strp1(ji,jj) + strp2(ji,jj) ) / numts_rm
                   strp2(ji,jj) = strp1(ji,jj)
@@ -661,6 +659,6 @@
          DO jj = 1, jpj 
             DO ji = 1, jpi
-               IF( a_i(ji,jj,jl) .GT. epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
-               ELSE                                   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1)
+               IF( a_i(ji,jj,jl) > epsi10 ) THEN   ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1) + a_i(ji,jj,jl)
+               ELSE                                ;   Gsum(ji,jj,jl) = Gsum(ji,jj,jl-1)
                ENDIF
             END DO
@@ -716,10 +714,10 @@
       ENDIF ! nn_partfun
 
-      IF( nn_rafting == 1 ) THEN      ! Ridging and rafting ice participation functions
+      IF( ln_rafting ) THEN      ! Ridging and rafting ice participation functions
          !
          DO jl = 1, jpl
             DO jj = 1, jpj 
                DO ji = 1, jpi
-                  IF ( athorn(ji,jj,jl) .GT. 0._wp ) THEN
+                  IF ( athorn(ji,jj,jl) > 0._wp ) THEN
 !!gm  TANH( -X ) = - TANH( X )  so can be computed only 1 time....
                      aridge(ji,jj,jl) = ( TANH (  rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) + 1.0 ) * 0.5 * athorn(ji,jj,jl)
@@ -727,10 +725,10 @@
                      IF ( araft(ji,jj,jl) < epsi06 )   araft(ji,jj,jl)  = 0._wp
                      aridge(ji,jj,jl) = MAX( athorn(ji,jj,jl) - araft(ji,jj,jl), 0.0 )
-                  ENDIF ! athorn
-               END DO ! ji
-            END DO ! jj
-         END DO ! jl
-
-      ELSE  ! nn_rafting = 0
+                  ENDIF
+               END DO
+            END DO
+         END DO
+
+      ELSE
          !
          DO jl = 1, jpl
@@ -740,11 +738,11 @@
       ENDIF
 
-      IF ( nn_rafting == 1 ) THEN
-
-         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) .GT. epsi10 ) THEN
+      IF( ln_rafting ) THEN
+
+         IF( MAXVAL(aridge + araft - athorn(:,:,1:jpl)) > epsi10 ) THEN
             DO jl = 1, jpl
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) .GT. epsi10 ) THEN
+                     IF ( aridge(ji,jj,jl) + araft(ji,jj,jl) - athorn(ji,jj,jl) > epsi10 ) THEN
                         WRITE(numout,*) ' ALERTE 96 : wrong participation function ... '
                         WRITE(numout,*) ' ji, jj, jl : ', ji, jj, jl
@@ -792,5 +790,5 @@
             DO ji = 1, jpi
 
-               IF (a_i(ji,jj,jl) .GT. epsi10 .AND. athorn(ji,jj,jl) .GT. 0.0 ) THEN
+               IF (a_i(ji,jj,jl) > epsi10 .AND. athorn(ji,jj,jl) > 0.0 ) THEN
                   hi = v_i(ji,jj,jl) / a_i(ji,jj,jl)
                   hrmean          = MAX(SQRT(rn_hstar*hi), hi*krdgmin)
@@ -988,5 +986,4 @@
          large_afrft = .false.
 
-!CDIR NODEP
          DO ij = 1, icells
             ji = indxi(ij)
@@ -1115,5 +1112,4 @@
          !--------------------------------------------------------------------
          DO jk = 1, nlay_i
-!CDIR NODEP
             DO ij = 1, icells
                ji = indxi(ij)
@@ -1142,5 +1138,4 @@
          IF( con_i ) THEN
             DO jk = 1, nlay_i
-!CDIR NODEP
                DO ij = 1, icells
                   ji = indxi(ij)
@@ -1152,5 +1147,4 @@
 
          IF( large_afrac ) THEN   ! there is a bug
-!CDIR NODEP
             DO ij = 1, icells
                ji = indxi(ij)
@@ -1164,5 +1158,4 @@
          ENDIF
          IF( large_afrft ) THEN  ! there is a bug
-!CDIR NODEP
             DO ij = 1, icells
                ji = indxi(ij)
@@ -1182,5 +1175,4 @@
          DO jl2  = 1, jpl 
             ! over categories to which ridged ice is transferred
-!CDIR NODEP
             DO ij = 1, icells
                ji = indxi(ij)
@@ -1215,5 +1207,4 @@
             ! Transfer ice energy to category jl2 by ridging
             DO jk = 1, nlay_i
-!CDIR NODEP
                DO ij = 1, icells
                   ji = indxi(ij)
@@ -1227,5 +1218,4 @@
          DO jl2 = 1, jpl 
 
-!CDIR NODEP
             DO ij = 1, icells
                ji = indxi(ij)
@@ -1242,11 +1232,10 @@
                   smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + smrft (ji,jj)    
                   oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirft2(ji,jj)    
-               ENDIF ! hraft
+               ENDIF
                !
-            END DO ! ij
+            END DO
 
             ! Transfer rafted ice energy to category jl2 
             DO jk = 1, nlay_i
-!CDIR NODEP
                DO ij = 1, icells
                   ji = indxi(ij)
@@ -1256,8 +1245,8 @@
                      e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)
                   ENDIF
-               END DO           ! ij
-            END DO !jk
-
-         END DO ! jl2
+               END DO
+            END DO
+
+         END DO
 
       END DO ! jl1 (deforming categories)
@@ -1331,6 +1320,6 @@
       !!-------------------------------------------------------------------
       INTEGER :: ios                 ! Local integer output status for namelist read
-      NAMELIST/namiceitdme/ rn_cs, rn_pe_rdg, rn_fsnowrdg, rn_fsnowrft,              & 
-        &                   rn_gstar, rn_astar, rn_hstar, nn_rafting, rn_hraft, rn_craft, rn_por_rdg, &
+      NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft,              & 
+        &                   rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, &
         &                   nn_partfun
       !!-------------------------------------------------------------------
@@ -1349,16 +1338,15 @@
          WRITE(numout,*)' lim_itd_me_init : ice parameters for mechanical ice redistribution '
          WRITE(numout,*)' ~~~~~~~~~~~~~~~'
-         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs              ', rn_cs 
-         WRITE(numout,*)'   Ratio of ridging work to PotEner change in ridging      rn_pe_rdg              ', rn_pe_rdg 
-         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg        ', rn_fsnowrdg 
-         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft        ', rn_fsnowrft 
-         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar           ', rn_gstar
-         WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar           ', rn_astar
-         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar           ', rn_hstar
-         WRITE(numout,*)'   Rafting of ice sheets or not                            nn_rafting        ', nn_rafting
-         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft       ', rn_hraft
-         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft           ', rn_craft  
-         WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg       ', rn_por_rdg
-         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun     ', nn_partfun
+         WRITE(numout,*)'   Fraction of shear energy contributing to ridging        rn_cs       = ', rn_cs 
+         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrdg = ', rn_fsnowrdg 
+         WRITE(numout,*)'   Fraction of snow volume conserved during ridging        rn_fsnowrft = ', rn_fsnowrft 
+         WRITE(numout,*)'   Fraction of total ice coverage contributing to ridging  rn_gstar    = ', rn_gstar
+         WRITE(numout,*)'   Equivalent to G* for an exponential part function       rn_astar    = ', rn_astar
+         WRITE(numout,*)'   Quantity playing a role in max ridged ice thickness     rn_hstar    = ', rn_hstar
+         WRITE(numout,*)'   Rafting of ice sheets or not                            ln_rafting  = ', ln_rafting
+         WRITE(numout,*)'   Parmeter thickness (threshold between ridge-raft)       rn_hraft    = ', rn_hraft
+         WRITE(numout,*)'   Rafting hyperbolic tangent coefficient                  rn_craft    = ', rn_craft  
+         WRITE(numout,*)'   Initial porosity of ridges                              rn_por_rdg  = ', rn_por_rdg
+         WRITE(numout,*)'   Switch for part. function (0) linear (1) exponential    nn_partfun  = ', nn_partfun
       ENDIF
       !
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 5070)
@@ -160,6 +160,6 @@
                zremap_flag(ji,jj) = 0
             ENDIF
-         END DO !ji
-      END DO !jj
+         END DO
+      END DO
 
       !-----------------------------------------------------------------------------------------------
@@ -201,5 +201,5 @@
          END DO
 
-      END DO !jl
+      END DO
 
       !-----------------------------------------------------------------------------------------------
@@ -244,6 +244,5 @@
       !-----------------------------------------------------------------------------------------------
       !- 7.1 g(h) for category 1 at start of time step
-      CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd),         &
-         &                  g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   &
+      CALL lim_itd_fitline( klbnd, zhb0, zhb1, zht_i_b(:,:,klbnd), g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd),   &
          &                  hR(:,:,klbnd), zremap_flag )
 
@@ -254,34 +253,32 @@
 
          !ji
-         IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN
+         IF( a_i(ii,ij,klbnd) > epsi10 ) THEN
             zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category
             ! ji, a_i > epsi10
-            IF (zdh0 .lt. 0.0) THEN !remove area from category 1
+            IF( zdh0 < 0.0 ) THEN !remove area from category 1
                ! ji, a_i > epsi10; zdh0 < 0
-               zdh0 = MIN(-zdh0,hi_max(klbnd))
+               zdh0 = MIN( -zdh0, hi_max(klbnd) )
 
                !Integrate g(1) from 0 to dh0 to estimate area melted
-               zetamax = MIN(zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd)
-               IF (zetamax.gt.0.0) THEN
+               zetamax = MIN( zdh0, hR(ii,ij,klbnd) ) - hL(ii,ij,klbnd)
+               IF( zetamax > 0.0 ) THEN
                   zx1  = zetamax
-                  zx2  = 0.5 * zetamax*zetamax 
+                  zx2  = 0.5 * zetamax * zetamax 
                   zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed
                   ! Constrain new thickness <= ht_i
-                  zdamax = a_i(ii,ij,klbnd) * & 
-                     (1.0 - ht_i(ii,ij,klbnd)/zht_i_b(ii,ij,klbnd)) ! zdamax > 0
+                  zdamax = a_i(ii,ij,klbnd) * (1.0 - ht_i(ii,ij,klbnd) / zht_i_b(ii,ij,klbnd) ) ! zdamax > 0
                   !ice area lost due to melting of thin ice
-                  zda0   = MIN(zda0, zdamax)
+                  zda0   = MIN( zda0, zdamax )
 
                   ! Remove area, conserving volume
-                  ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) & 
-                     * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 )
+                  ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 )
                   a_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) - zda0
-                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem-useless ?
-               ENDIF     ! zetamax > 0
+                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) * ht_i(ii,ij,klbnd) ! clem-useless ?
+               ENDIF
                ! ji, a_i > epsi10
 
             ELSE ! if ice accretion
                ! ji, a_i > epsi10; zdh0 > 0
-               zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd)) 
+               zhbnew(ii,ij,klbnd-1) = MIN( zdh0, hi_max(klbnd) ) 
                ! zhbnew was 0, and is shifted to the right to account for thin ice
                ! growth in openwater (F0 = f1)
@@ -295,6 +292,6 @@
       !- 7.3 g(h) for each thickness category  
       DO jl = klbnd, kubnd
-         CALL lim_itd_fitline(jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), &
-            g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag)
+         CALL lim_itd_fitline( jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl),  &
+            &                  g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), zremap_flag )
       END DO
 
@@ -316,9 +313,9 @@
             ij = nind_j(ji)
 
-            IF (zhbnew(ii,ij,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1
+            IF (zhbnew(ii,ij,jl) > hi_max(jl)) THEN ! transfer from jl to jl+1
 
                ! left and right integration limits in eta space
-               zvetamin(ji) = MAX(hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl)
-               zvetamax(ji) = MIN(zhbnew(ii,ij,jl), hR(ii,ij,jl)) - hL(ii,ij,jl)
+               zvetamin(ji) = MAX( hi_max(jl), hL(ii,ij,jl) ) - hL(ii,ij,jl)
+               zvetamax(ji) = MIN (zhbnew(ii,ij,jl), hR(ii,ij,jl) ) - hL(ii,ij,jl)
                zdonor(ii,ij,jl) = jl
 
@@ -327,24 +324,24 @@
                ! left and right integration limits in eta space
                zvetamin(ji) = 0.0
-               zvetamax(ji) = MIN(hi_max(jl), hR(ii,ij,jl+1)) - hL(ii,ij,jl+1)
+               zvetamax(ji) = MIN( hi_max(jl), hR(ii,ij,jl+1) ) - hL(ii,ij,jl+1)
                zdonor(ii,ij,jl) = jl + 1
 
             ENDIF  ! zhbnew(jl) > hi_max(jl)
 
-            zetamax = MAX(zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin
+            zetamax = MAX( zvetamax(ji), zvetamin(ji) ) ! no transfer if etamax < etamin
             zetamin = zvetamin(ji)
 
             zx1  = zetamax - zetamin
-            zwk1 = zetamin*zetamin
-            zwk2 = zetamax*zetamax
-            zx2  = 0.5 * (zwk2 - zwk1)
+            zwk1 = zetamin * zetamin
+            zwk2 = zetamax * zetamax
+            zx2  = 0.5 * ( zwk2 - zwk1 )
             zwk1 = zwk1 * zetamin
             zwk2 = zwk2 * zetamax
-            zx3  = 1.0/3.0 * (zwk2 - zwk1)
+            zx3  = 1.0 / 3.0 * ( zwk2 - zwk1 )
             nd   = zdonor(ii,ij,jl)
             zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1
             zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + zdaice(ii,ij,jl)*hL(ii,ij,nd)
 
-         END DO ! ji
+         END DO
       END DO ! jl klbnd -> kubnd - 1
 
@@ -365,5 +362,5 @@
             ht_i(ii,ij,1) = rn_himin
          ENDIF
-      END DO !ji
+      END DO
 
       !!----------------------------------------------------------------------------------------------
@@ -401,6 +398,5 @@
 
 
-   SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice,   &
-      &                        g0, g1, hL, hR, zremap_flag )
+   SUBROUTINE lim_itd_fitline( num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag )
       !!------------------------------------------------------------------
       !!                ***  ROUTINE lim_itd_fitline ***
@@ -442,6 +438,6 @@
                ! Change hL or hR if hice falls outside central third of range
 
-               zh13 = 1.0/3.0 * (2.0*hL(ji,jj) + hR(ji,jj))
-               zh23 = 1.0/3.0 * (hL(ji,jj) + 2.0*hR(ji,jj))
+               zh13 = 1.0 / 3.0 * ( 2.0 * hL(ji,jj) + hR(ji,jj) )
+               zh23 = 1.0 / 3.0 * ( hL(ji,jj) + 2.0 * hR(ji,jj) )
 
                IF    ( hice(ji,jj) < zh13 ) THEN   ;   hR(ji,jj) = 3._wp * hice(ji,jj) - 2._wp * hL(ji,jj)
@@ -454,6 +450,6 @@
                zwk1 = 6._wp * a_i(ji,jj,num_cat) * zdhr
                zwk2 = ( hice(ji,jj) - hL(ji,jj) ) * zdhr
-               g0(ji,jj) = zwk1 * ( 2._wp/3._wp - zwk2 )
-               g1(ji,jj) = 2._wp * zdhr * zwk1 * (zwk2 - 0.5)
+               g0(ji,jj) = zwk1 * ( 2._wp / 3._wp - zwk2 )
+               g1(ji,jj) = 2._wp * zdhr * zwk1 * ( zwk2 - 0.5 )
                !
             ELSE                   ! remap_flag = .false. or a_i < epsi10 
@@ -516,5 +512,5 @@
 
       DO jl = klbnd, kubnd
-         zaTsfn(:,:,jl) = a_i(:,:,jl)*t_su(:,:,jl)
+         zaTsfn(:,:,jl) = a_i(:,:,jl) * t_su(:,:,jl)
       END DO
 
@@ -539,13 +535,11 @@
             DO ji = 1, jpi
 
-               IF (zdonor(ji,jj,jl) .GT. 0) THEN
+               IF (zdonor(ji,jj,jl) > 0) THEN
                   jl1 = zdonor(ji,jj,jl)
 
-                  IF (zdaice(ji,jj,jl) .LT. 0.0) THEN
-                     IF (zdaice(ji,jj,jl) .GT. -epsi10) THEN
-                        IF ( ( jl1.EQ.jl   .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) )           &
-                           .OR.                                      &
-                           ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) )           &  
-                           ) THEN                                                             
+                  IF (zdaice(ji,jj,jl) < 0.0) THEN
+                     IF (zdaice(ji,jj,jl) > -epsi10) THEN
+                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  &
+                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN
                            zdaice(ji,jj,jl) = a_i(ji,jj,jl1)  ! shift entire category
                            zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
@@ -559,10 +553,8 @@
                   ENDIF
 
-                  IF (zdvice(ji,jj,jl) .LT. 0.0) THEN
-                     IF (zdvice(ji,jj,jl) .GT. -epsi10 ) THEN
-                        IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) )     &
-                           .OR.                                     &
-                           ( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) &
-                           ) THEN
+                  IF (zdvice(ji,jj,jl) < 0.0) THEN
+                     IF (zdvice(ji,jj,jl) > -epsi10 ) THEN
+                        IF ( ( jl1 == jl   .AND. ht_i(ji,jj,jl1) >  hi_max(jl) ) .OR.  &
+                             ( jl1 == jl+1 .AND. ht_i(ji,jj,jl1) <= hi_max(jl) ) ) THEN
                            zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category
                            zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
@@ -577,6 +569,6 @@
 
                   ! If daice is close to aicen, set daice = aicen.
-                  IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - epsi10 ) THEN
-                     IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+epsi10) THEN
+                  IF (zdaice(ji,jj,jl) > a_i(ji,jj,jl1) - epsi10 ) THEN
+                     IF (zdaice(ji,jj,jl) < a_i(ji,jj,jl1)+epsi10) THEN
                         zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
                         zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
@@ -586,6 +578,6 @@
                   ENDIF
 
-                  IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-epsi10) THEN
-                     IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+epsi10) THEN
+                  IF (zdvice(ji,jj,jl) > v_i(ji,jj,jl1)-epsi10) THEN
+                     IF (zdvice(ji,jj,jl) < v_i(ji,jj,jl1)+epsi10) THEN
                         zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
                         zdvice(ji,jj,jl) = v_i(ji,jj,jl1) 
@@ -596,8 +588,8 @@
 
                ENDIF               ! donor > 0
-            END DO                   ! i
-         END DO                 ! j
-
-      END DO !jl
+            END DO
+         END DO
+
+      END DO
 
       !-------------------------------------------------------------------------------
@@ -609,9 +601,9 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF (zdaice(ji,jj,jl) .GT. 0.0 ) THEN ! daice(n) can be < puny
+               IF (zdaice(ji,jj,jl) > 0.0 ) THEN ! daice(n) can be < puny
                   nbrem = nbrem + 1
                   nind_i(nbrem) = ji
                   nind_j(nbrem) = jj
-               ENDIF ! tmask
+               ENDIF
             END DO
          END DO
@@ -622,8 +614,8 @@
 
             jl1 = zdonor(ii,ij,jl)
-            rswitch             = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) )
-            zworka(ii,ij)   = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * rswitch
+            rswitch       = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) )
+            zworka(ii,ij) = zdvice(ii,ij,jl) / MAX( v_i(ii,ij,jl1), epsi10 ) * rswitch
             IF( jl1 == jl) THEN   ;   jl2 = jl1+1
-            ELSE                    ;   jl2 = jl 
+            ELSE                  ;   jl2 = jl 
             ENDIF
 
@@ -682,5 +674,5 @@
             zaTsfn(ii,ij,jl2)  = zaTsfn(ii,ij,jl2) + zdaTsf 
 
-         END DO                 ! ji
+         END DO
 
          !------------------
@@ -689,5 +681,4 @@
 
          DO jk = 1, nlay_i
-!CDIR NODEP
             DO ji = 1, nbrem
                ii = nind_i(ji)
@@ -695,5 +686,5 @@
 
                jl1 = zdonor(ii,ij,jl)
-               IF (jl1 .EQ. jl) THEN
+               IF (jl1 == jl) THEN
                   jl2 = jl+1
                ELSE             ! n1 = n+1
@@ -704,6 +695,6 @@
                e_i(ii,ij,jk,jl1) =  e_i(ii,ij,jk,jl1) - zdeice
                e_i(ii,ij,jk,jl2) =  e_i(ii,ij,jk,jl2) + zdeice 
-            END DO              ! ji
-         END DO                 ! jk
+            END DO
+         END DO
 
       END DO                   ! boundaries, 1 to ncat-1
@@ -719,12 +710,12 @@
                   ht_i(ji,jj,jl)  =  v_i   (ji,jj,jl) / a_i(ji,jj,jl) 
                   t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl) 
-                  rswitch         =  1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes
+                  rswitch         =  1.0 - MAX( 0.0, SIGN( 1.0, -v_s(ji,jj,jl) + epsi10 ) ) !0 if no ice and 1 if yes
                ELSE
                   ht_i(ji,jj,jl)  = 0._wp
                   t_su(ji,jj,jl)  = rtt
                ENDIF
-            END DO                 ! ji
-         END DO                 ! jj
-      END DO                    ! jl
+            END DO
+         END DO
+      END DO
       !
       CALL wrk_dealloc( jpi,jpj,jpl, zaTsfn )
@@ -836,6 +827,6 @@
                   zdvice(ji,jj,jl)  = v_i(ji,jj,jl) - ( a_i(ji,jj,jl) - zdaice(ji,jj,jl) ) * ( hi_max(jl) - epsi10 )
                ENDIF
-            END DO                 ! ji
-         END DO                 ! jj
+            END DO
+         END DO
          IF(lk_mpp)   CALL mpp_max( zshiftflag )
 
@@ -870,6 +861,6 @@
                   zdvice(ji,jj,jl) = v_i(ji,jj,jl+1)
                ENDIF
-            END DO                 ! ji
-         END DO                 ! jj
+            END DO
+         END DO
 
          IF(lk_mpp)   CALL mpp_max( zshiftflag )
@@ -890,9 +881,9 @@
 !                  a_i (ji,jj,jl+1) = v_i(ji,jj,jl+1) / ht_i(ji,jj,jl+1) 
 !               ENDIF
-!            END DO                 ! ji
-!         END DO                 ! jj
+!            END DO
+!         END DO
          ! clem-change end
 
-      END DO                    ! jl
+      END DO
 
       !------------------------------------------------------------------------------
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 5070)
@@ -291,6 +291,6 @@
             ! include it later
 
-            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) / e1u(ji,jj)
-            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) / e2v(ji,jj)
+            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
+            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
 
             za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
@@ -357,11 +357,11 @@
                   &            ) * r1_e12t(ji,jj)
 
-               zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
-                  &         - ( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+                  &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
                   &         ) * r1_e12t(ji,jj)
 
                !
-               zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
-                  &         + ( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
                   &         ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) )   &
                   &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
@@ -391,10 +391,10 @@
 
                !- Calculate Delta on corners
-               zddc  = (  ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
-                  &     + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
+               zddc  = (  ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
+                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
                   &    ) * r1_e12f(ji,jj)
 
-               zdtc  = (- ( v_ice1(ji,jj+1) / e1u(ji,jj+1) - v_ice1(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
-                  &     + ( u_ice2(ji+1,jj) / e2v(ji+1,jj) - u_ice2(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
+               zdtc  = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
+                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
                   &    ) * r1_e12f(ji,jj)
 
@@ -420,11 +420,11 @@
                !- contribution of zs1, zs2 and zs12 to zf1
                zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  &
-                  &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) / e2u(ji,jj)          &
-                  &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) / e1u(ji,jj)  &
+                  &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj)          &
+                  &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj)  &
                   &                ) * r1_e12u(ji,jj)
                ! contribution of zs1, zs2 and zs12 to zf2
                zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  &
-                  &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) / e1v(ji,jj)          &
-                  &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) / e2v(ji,jj)  &
+                  &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj)          &
+                  &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj)  &
                   &               )  * r1_e12v(ji,jj)
             END DO
@@ -609,11 +609,11 @@
                   &            ) * r1_e12t(ji,jj)
 
-               zdt(ji,jj) = ( ( u_ice(ji,jj) / e2u(ji,jj) - u_ice(ji-1,jj) / e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  &
-                  &          -( v_ice(ji,jj) / e1v(ji,jj) - v_ice(ji,jj-1) / e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  &
+               zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  &
+                  &          -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  &
                   &         ) * r1_e12t(ji,jj)
                !
                ! SB modif because ocean has no slip boundary condition 
-               zds(ji,jj) = ( ( u_ice(ji,jj+1) / e1u(ji,jj+1) - u_ice(ji,jj) / e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
-                  &          +( v_ice(ji+1,jj) / e2v(ji+1,jj) - v_ice(ji,jj) / e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
+               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
+                  &          +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
                   &         ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) )                                     &
                   &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
@@ -677,5 +677,5 @@
             DO jj = k_j1+1, k_jpj-1
                DO ji = 2, jpim1
-                  IF (zpresh(ji,jj) .GT. 1.0) THEN
+                  IF (zpresh(ji,jj) > 1.0) THEN
                      sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
                      sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )
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 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 5070)
@@ -631,5 +631,5 @@
             rswitch     = 1._wp - MAX( 0._wp , SIGN( 1._wp , - a_i_1d(ji) + epsi20 ) )
             ht_i_1d(ji) = rswitch * zv / MAX( a_i_1d(ji), epsi20 )
-            ! adjust thickness
+            ! retrieve total concentration
             at_i_1d(ji) = a_i_1d(ji)
          END IF
@@ -651,7 +651,6 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicethd/ rn_hnewice, nn_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,    &
-         &                rn_himin, parsub, rn_betas,                          & 
-         &                rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &
+      NAMELIST/namicethd/ rn_hnewice, ln_frazil, rn_maxfrazb, rn_vfrazb, rn_Cfrazb,                       &
+         &                rn_himin, parsub, rn_betas, rn_kappa_i, nn_conv_dif, rn_terr_dif, nn_ice_thcon, &
          &                nn_monocat
       !!-------------------------------------------------------------------
@@ -682,17 +681,17 @@
          WRITE(numout,*)
          WRITE(numout,*)'   Namelist of ice parameters for ice thermodynamic computation '
-         WRITE(numout,*)'      ice thick. for lateral accretion                        rn_hnewice      = ', rn_hnewice
-         WRITE(numout,*)'      Frazil ice thickness as a function of wind or not       nn_frazil     = ', nn_frazil
-         WRITE(numout,*)'      Maximum proportion of frazil ice collecting at bottom   rn_maxfrazb     = ', rn_maxfrazb
-         WRITE(numout,*)'      Thresold relative drift speed for collection of frazil  rn_vfrazb       = ', rn_vfrazb
-         WRITE(numout,*)'      Squeezing coefficient for collection of frazil          rn_Cfrazb       = ', rn_Cfrazb
-         WRITE(numout,*)'      minimum ice thickness                                   rn_himin       = ', rn_himin 
+         WRITE(numout,*)'      ice thick. for lateral accretion                        rn_hnewice   = ', rn_hnewice
+         WRITE(numout,*)'      Frazil ice thickness as a function of wind or not       ln_frazil    = ', ln_frazil
+         WRITE(numout,*)'      Maximum proportion of frazil ice collecting at bottom   rn_maxfrazb  = ', rn_maxfrazb
+         WRITE(numout,*)'      Thresold relative drift speed for collection of frazil  rn_vfrazb    = ', rn_vfrazb
+         WRITE(numout,*)'      Squeezing coefficient for collection of frazil          rn_Cfrazb    = ', rn_Cfrazb
+         WRITE(numout,*)'      minimum ice thickness                                   rn_himin     = ', rn_himin 
          WRITE(numout,*)'      numerical carac. of the scheme for diffusion in ice '
          WRITE(numout,*)'      switch for snow sublimation  (=1) or not (=0)           parsub       = ', parsub  
-         WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          rn_betas        = ', rn_betas
-         WRITE(numout,*)'      extinction radiation parameter in sea ice               rn_kappa_i      = ', rn_kappa_i
+         WRITE(numout,*)'      coefficient for ice-lead partition of snowfall          rn_betas     = ', rn_betas
+         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,*)'      switch for comp. of thermal conductivity in the ice     nn_ice_thcon = ', nn_ice_thcon
          WRITE(numout,*)'      check heat conservation in the ice/snow                 con_i        = ', con_i
          WRITE(numout,*)'      virtual ITD mono-category parameterizations (1) or not  nn_monocat   = ', nn_monocat
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 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 5070)
@@ -276,5 +276,5 @@
 
             ! updates available heat + thickness
-            zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )
+            zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_1d(ji,jk) )
             ht_s_1d(ji) = MAX( 0._wp , ht_s_1d(ji) + zdeltah(ji,jk) )
 
@@ -573,6 +573,6 @@
            
             ENDIF
-         END DO ! ji
-      END DO ! jk
+         END DO
+      END DO
 
       !-------------------------------------------
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 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 5070)
@@ -517,5 +517,5 @@
 
          DO ji = kideb , kiut
-            IF ( ht_s_1d(ji).gt.0.0 ) THEN
+            IF ( ht_s_1d(ji) > 0.0 ) THEN
                !
                !------------------------------------------------------------------------------|
@@ -541,5 +541,5 @@
                ENDIF
 
-               IF ( t_su_1d(ji) .LT. rtt ) THEN
+               IF ( t_su_1d(ji) < rtt ) THEN
 
                   !------------------------------------------------------------------------------|
@@ -587,5 +587,5 @@
                !------------------------------------------------------------------------------|
                !
-               IF (t_su_1d(ji) .LT. rtt) THEN
+               IF (t_su_1d(ji) < rtt) THEN
                   !
                   !------------------------------------------------------------------------------|
@@ -705,5 +705,5 @@
          DO ji = kideb , kiut
             ! snow temperatures      
-            IF (ht_s_1d(ji).GT.0._wp) &
+            IF (ht_s_1d(ji) > 0._wp) &
                t_s_1d(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
                *  t_i_1d(ji,1))/zdiagbis(ji,nlay_s+1) &
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 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 5070)
@@ -156,5 +156,5 @@
       hicol(:,:) = rn_hnewice
 
-      IF( nn_frazil == 1 ) THEN
+      IF( ln_frazil ) THEN
 
          !--------------------
@@ -221,5 +221,5 @@
                   iterate_frazil = .true.
 
-                  DO WHILE ( iter .LT. 100 .AND. iterate_frazil ) 
+                  DO WHILE ( iter < 100 .AND. iterate_frazil ) 
                      zf = ( hicol(ji,jj) - zhicrit ) * ( hicol(ji,jj)**2 - zhicrit**2 ) &
                         - hicol(ji,jj) * zhicrit * ztwogp * zvrel2
@@ -321,5 +321,5 @@
             zh_newice(ji) = rn_hnewice
          END DO
-         IF( nn_frazil == 1 ) zh_newice(1:nbpac) = hicol_1d(1:nbpac)
+         IF( ln_frazil ) zh_newice(1:nbpac) = hicol_1d(1:nbpac)
 
          !----------------------
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 5070)
@@ -150,12 +150,12 @@
          WRITE(numout,*) 'lim_thd_sal_init : Ice parameters for salinity '
          WRITE(numout,*) '~~~~~~~~~~~~~~~~'
-         WRITE(numout,*) ' switch for salinity nn_icesal        : ', nn_icesal
-         WRITE(numout,*) ' bulk salinity value if nn_icesal = 1 : ', rn_icesal
-         WRITE(numout,*) ' restoring salinity for GD          : ', rn_sal_gd
-         WRITE(numout,*) ' restoring time for GD              : ', rn_time_gd
-         WRITE(numout,*) ' restoring salinity for flushing    : ', rn_sal_fl
-         WRITE(numout,*) ' restoring time for flushing        : ', rn_time_fl
-         WRITE(numout,*) ' Maximum tolerated ice salinity     : ', rn_simax
-         WRITE(numout,*) ' Minimum tolerated ice salinity     : ', rn_simin
+         WRITE(numout,*) '   switch for salinity nn_icesal        = ', nn_icesal
+         WRITE(numout,*) '   bulk salinity value if nn_icesal = 1 = ', rn_icesal
+         WRITE(numout,*) '   restoring salinity for GD            = ', rn_sal_gd
+         WRITE(numout,*) '   restoring time for GD                = ', rn_time_gd
+         WRITE(numout,*) '   restoring salinity for flushing      = ', rn_sal_fl
+         WRITE(numout,*) '   restoring time for flushing          = ', rn_time_fl
+         WRITE(numout,*) '   Maximum tolerated ice salinity       = ', rn_simax
+         WRITE(numout,*) '   Minimum tolerated ice salinity       = ', rn_simin
       ENDIF
       !
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 5070)
@@ -106,4 +106,5 @@
          ! mass and salt flux init
          zviold(:,:,:) = v_i(:,:,:)
+         zvsold(:,:,:) = v_s(:,:,:)
          zeiold(:,:)   = SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) 
          zesold(:,:)   = SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) 
@@ -131,6 +132,6 @@
 
          ! If ice drift field is too fast, use an appropriate time step for advection.         
-         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice / e1u(:,:) )         ! CFL test for stability
-         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice / e2v(:,:) ) )
+         zcfl  =            MAXVAL( ABS( u_ice(:,:) ) * rdt_ice * r1_e1u(:,:) )         ! CFL test for stability
+         zcfl  = MAX( zcfl, MAXVAL( ABS( v_ice(:,:) ) * rdt_ice * r1_e2v(:,:) ) )
          IF(lk_mpp )   CALL mpp_max( zcfl )
 
Index: branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 5070)
@@ -26,5 +26,5 @@
    REAL(wp), PUBLIC ::   rn_hnewice  !: thickness for new ice formation (m)
 
-   INTEGER , PUBLIC ::   nn_frazil   !: use of frazil ice collection as function of wind (1) or not (0)
+   LOGICAL , PUBLIC ::   ln_frazil   !: use of frazil ice collection as function of wind (T) or not (F)
 
    !!-----------------------------
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 5068)
+++ branches/2015/dev_r5044_CNRS_LIM3CLEAN/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 5070)
@@ -369,7 +369,7 @@
          WRITE(numout,*) 'ice_run : 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,*) '   number of ice  categories                               = ', jpl
+         WRITE(numout,*) '   number of ice  layers                                   = ', nlay_i
+         WRITE(numout,*) '   number of snow layers                                   = ', nlay_s
          WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn
          WRITE(numout,*) '   maximum ice concentration                               = ', rn_amax 
