Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 3517)
@@ -199,5 +199,5 @@
             !
             !                          !------------------------------------------!
-            !                          !      mass flux at the ocean surface      !
+            !                          !  mass and salt flux at the ocean surface !
             !                          !------------------------------------------!
             !
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 3517)
@@ -158,9 +158,10 @@
    !! * Share Module variables
    !!--------------------------------------------------------------------------
-   INTEGER , PUBLIC ::   nstart    !: iteration number of the begining of the run 
-   INTEGER , PUBLIC ::   nlast     !: iteration number of the end of the run 
-   INTEGER , PUBLIC ::   nitrun    !: number of iteration
-   INTEGER , PUBLIC ::   numit     !: iteration number
-   REAL(wp), PUBLIC ::   rdt_ice   !: ice time step
+   INTEGER , PUBLIC ::   nstart      !: iteration number of the begining of the run 
+   INTEGER , PUBLIC ::   nlast       !: iteration number of the end of the run 
+   INTEGER , PUBLIC ::   nitrun      !: number of iteration
+   INTEGER , PUBLIC ::   numit       !: iteration number
+   REAL(wp), PUBLIC ::   rdt_ice     !: ice time step
+   REAL(wp), PUBLIC ::   r1_rdtice   !: = 1. / rdt_ice
 
    !                                          !!** ice-dynamic namelist (namicedyn) **
@@ -201,8 +202,7 @@
    !                                              !!** ice-salinity namelist (namicesal) **
    INTEGER , PUBLIC ::   num_sal     = 1           !: salinity configuration used in the model
-   !                                               ! 1 - s constant in space and time
+   !                                               ! 1 - constant salinity in both space and time
    !                                               ! 2 - prognostic salinity (s(z,t))
    !                                               ! 3 - salinity profile, constant in time
-   !                                               ! 4 - salinity variations affect only ice thermodynamics
    INTEGER , PUBLIC ::   sal_prof    = 1           !: salinity profile or not 
    INTEGER , PUBLIC ::   thcon_i_swi = 1           !: thermal conductivity: =1 Untersteiner (1964) ; =2 Pringle et al (2007)
@@ -278,11 +278,11 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qfvbq       !: store energy in case of total lateral ablation (?)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dmgwi       !: Variation of the mass of snow ice
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsalt_res   !: Residual salt flux due to correction of ice thickness
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsbri       !: Salt flux due to brine rejection
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsalt_rpo   !: Salt flux associated with porous ridged ice formation
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_rpo   !: Heat flux associated with porous ridged ice formation
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsalt_res   !: Residual salt flux due to correction of ice thickness [???]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsbri       !: Salt flux due to brine rejection                      [???]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fsalt_rpo   !: Salt flux associated with porous ridged ice formation [???]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_rpo   !: Heat flux associated with porous ridged ice formation [???]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhbri       !: heat flux due to brine rejection
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fmmec       !: Mass flux due to snow loss during compression
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fseqv       !: Equivalent salt flux due to ice growth/melt
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fseqv       !: salt flux due to ice growth/melt   [units ???]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhmec       !: Heat flux due to snow loss during compression
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fheat_res   !: Residual heat flux due to correction of ice thickness
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 3517)
@@ -79,5 +79,6 @@
       CALL lim_thd_sal_init            ! set ice salinity parameters
       !
-      rdt_ice = nn_fsbc * rdttra(1)    ! sea-ice timestep
+      rdt_ice   = nn_fsbc * rdttra(1)  ! sea-ice timestep
+      r1_rdtice = 1._wp / rdt_ice      ! sea-ice timestep inverse
       !
       CALL lim_msh                     ! ice mesh initialization
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90	(revision 3517)
@@ -88,5 +88,5 @@
             zs2new  = MIN(  2.0 * zslpmax - 0.3334 * ABS( zs1new ),      &
                &            MAX( ABS( zs1new ) - zslpmax, psxx(ji,jj) )  )
-            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
+            zin0    = ( 1.0 - MAX( rzero, SIGN( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
 
             ps0 (ji,jj) = zslpmax  
@@ -273,5 +273,5 @@
             zs2new  = MIN(  ( 2.0 * zslpmax - 0.3334 * ABS( zs1new ) ),   &
                &             MAX( ABS( zs1new )-zslpmax, psyy(ji,jj) )  )
-            zin0    = ( 1.0 - MAX( rzero, sign ( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
+            zin0    = ( 1.0 - MAX( rzero, SIGN( rone, -zslpmax) ) ) * tms(ji,jj)   ! Case of empty boxes & Apply mask
             !
             ps0 (ji,jj) = zslpmax  
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90	(revision 3517)
@@ -180,5 +180,5 @@
                vinfor(43) = vinfor(43) + diag_dyn_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp 
                vinfor(45) = vinfor(45) + dv_dt_thd(ji,jj,5)*aire(ji,jj) * 1.e-12_wp
-               vinfor(47) = vinfor(47) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp / rdt_ice ! volume acc in OW
+               vinfor(47) = vinfor(47) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp * r1_rdtice   ! volume acc in OW
             ENDIF
          END DO
@@ -240,8 +240,6 @@
       vinfor(84) = 0.0
       DO ji = 134, 138
-         vinfor(83) = vinfor(83) - v_ice(ji,jj) * & 
-            e1t(ji,jj)*at_i(ji,jj)*rdt_ice * 1.e-12_wp
-         vinfor(84) = vinfor(84) - v_ice(ji,jj) * & 
-            e1t(ji,jj)*vt_i(ji,jj)*rdt_ice * 1.e-12_wp
+         vinfor(83) = vinfor(83) - v_ice(ji,jj) * e1t(ji,jj)*at_i(ji,jj)*rdt_ice * 1.e-12_wp
+         vinfor(84) = vinfor(84) - v_ice(ji,jj) * e1t(ji,jj)*vt_i(ji,jj)*rdt_ice * 1.e-12_wp
       END DO
 
@@ -331,5 +329,5 @@
                vinfor(44) = vinfor(44) + diag_dyn_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp 
                vinfor(46) = vinfor(46) + dv_dt_thd(ji,jj,5)*aire(ji,jj) * 1.e-12_wp
-               vinfor(48) = vinfor(48) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp / rdt_ice ! volume acc in OW
+               vinfor(48) = vinfor(48) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp * r1_rdtice   ! volume acc in OW
             ENDIF
          END DO
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90	(revision 3517)
@@ -38,7 +38,7 @@
    PUBLIC   lim_itd_me_alloc        ! called by iceini.F90
 
-   REAL(wp)  ::   epsi11 = 1.e-11_wp   ! constant values
-   REAL(wp)  ::   epsi10 = 1.e-10_wp   ! constant values
-   REAL(wp)  ::   epsi06 = 1.e-06_wp   ! constant values
+   REAL(wp) ::   epsi11 = 1.e-11_wp   ! constant values
+   REAL(wp) ::   epsi10 = 1.e-10_wp   ! constant values
+   REAL(wp) ::   epsi06 = 1.e-06_wp   ! constant values
 
    !-----------------------------------------------------------------------
@@ -47,8 +47,8 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   asum     ! sum of total ice and open water area
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aksum    ! ratio of area removed to area ridged
-
+   !
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   athorn   ! participation function; fraction of ridging/
    !                                                           !  closing associated w/ category n
-
+   !
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmin    ! minimum ridge thickness
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   hrmax    ! maximum ridge thickness
@@ -70,4 +70,5 @@
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   dvirdgdt   ! rate of ice volume ridged (m/s)
    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   opening    ! rate of opening due to divergence/shear (1/s)
+   !
    !!----------------------------------------------------------------------
    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
@@ -126,7 +127,7 @@
       INTEGER ::   ji, jj, jk, jl   ! dummy loop index
       INTEGER ::   niter, nitermax = 20   ! local integer 
-      LOGICAL  ::   asum_error              ! flag for asum .ne. 1
-      INTEGER  ::   iterate_ridging         ! if true, repeat the ridging
-      REAL(wp) ::   w1, tmpfac, dti         ! local scalar
+      LOGICAL  ::   asum_error            ! flag for asum .ne. 1
+      INTEGER  ::   iterate_ridging       ! if true, repeat the ridging
+      REAL(wp) ::   w1, tmpfac            ! local scalar
       CHARACTER (len = 15) ::   fieldid
       REAL(wp), POINTER, DIMENSION(:,:) ::   closing_net     ! net rate at which area is removed    (1/s)
@@ -152,12 +153,12 @@
       ! 1) Thickness categories boundaries, ice / o.w. concentrations, init_ons
       !-----------------------------------------------------------------------------!
-      ! Set hi_max(ncat) to a big value to ensure that all ridged ice 
-      ! is thinner than hi_max(ncat).
+      ! Set hi_max(ncat) to a big value to ensure that all ridged ice is thinner than hi_max(ncat).
 
       hi_max(jpl) = 999.99
 
-      Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0      ! proport const for PE
-      CALL lim_itd_me_ridgeprep ! prepare ridging
-
+      Cp = 0.5 * grav * (rau0-rhoic) * rhoic / rau0                ! proport const for PE
+      !
+      CALL lim_itd_me_ridgeprep                                    ! prepare ridging
+      !
       IF( con_i)   CALL lim_column_sum( jpl, v_i, vt_i_init )      ! conservation check
 
@@ -166,8 +167,8 @@
             msnow_mlt(ji,jj) = 0._wp
             esnow_mlt(ji,jj) = 0._wp
-            dardg1dt (ji,jj)  = 0._wp
-            dardg2dt (ji,jj)  = 0._wp
-            dvirdgdt (ji,jj)  = 0._wp
-            opening  (ji,jj)  = 0._wp
+            dardg1dt (ji,jj) = 0._wp
+            dardg2dt (ji,jj) = 0._wp
+            dvirdgdt (ji,jj) = 0._wp
+            opening  (ji,jj) = 0._wp
 
             !-----------------------------------------------------------------------------!
@@ -201,5 +202,5 @@
             ! to give asum = 1.0 after ridging.
 
-            divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) / rdt_ice  ! asum found in ridgeprep
+            divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep
 
             IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) )
@@ -207,6 +208,5 @@
             ! 2.3 opning
             !------------
-            ! Compute the (non-negative) opening rate that will give 
-            ! asum = 1.0 after ridging.
+            ! Compute the (non-negative) opening rate that will give asum = 1.0 after ridging.
             opning(ji,jj) = closing_net(ji,jj) + divu_adv(ji,jj)
          END DO
@@ -257,5 +257,5 @@
                   IF ( a_i(ji,jj,jl) > epsi11 .AND. athorn(ji,jj,jl) > 0._wp )THEN
                      w1 = athorn(ji,jj,jl) * closing_gross(ji,jj) * rdt_ice
-                     IF ( w1 > a_i(ji,jj,jl) ) THEN
+                     IF ( w1  >  a_i(ji,jj,jl) ) THEN
                         tmpfac = a_i(ji,jj,jl) / w1
                         closing_gross(ji,jj) = closing_gross(ji,jj) * tmpfac
@@ -291,5 +291,5 @@
                ELSE
                   iterate_ridging    = 1
-                  divu_adv   (ji,jj) = (1._wp - asum(ji,jj)) / rdt_ice
+                  divu_adv   (ji,jj) = (1._wp - asum(ji,jj)) * r1_rdtice
                   closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) )
                   opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) )
@@ -308,5 +308,5 @@
 
          IF( iterate_ridging == 1 ) THEN
-            IF( niter .GT. nitermax ) THEN
+            IF( niter  >  nitermax ) THEN
                WRITE(numout,*) ' ALERTE : non-converging ridging scheme '
                WRITE(numout,*) ' niter, iterate_ridging ', niter, iterate_ridging
@@ -323,6 +323,4 @@
       ! Update fresh water and heat fluxes due to snow melt.
 
-      dti = 1._wp / rdt_ice
-
       asum_error = .false. 
 
@@ -330,16 +328,16 @@
          DO ji = 1, jpi
 
-            IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) asum_error = .true.
-
-            dardg1dt(ji,jj) = dardg1dt(ji,jj) * dti
-            dardg2dt(ji,jj) = dardg2dt(ji,jj) * dti
-            dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * dti
-            opening (ji,jj) = opening (ji,jj) * dti
+            IF( ABS( asum(ji,jj) - 1.0 ) > epsi11 )   asum_error = .true.
+
+            dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice
+            dardg2dt(ji,jj) = dardg2dt(ji,jj) * r1_rdtice
+            dvirdgdt(ji,jj) = dvirdgdt(ji,jj) * r1_rdtice
+            opening (ji,jj) = opening (ji,jj) * r1_rdtice
 
             !-----------------------------------------------------------------------------!
             ! 5) Heat, salt and freshwater fluxes
             !-----------------------------------------------------------------------------!
-            fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * dti     ! fresh water source for ocean
-            fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * dti     ! heat sink for ocean
+            fmmec(ji,jj) = fmmec(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean
+            fhmec(ji,jj) = fhmec(ji,jj) + esnow_mlt(ji,jj) * r1_rdtice     ! heat sink for ocean
 
          END DO
@@ -349,5 +347,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF (ABS(asum(ji,jj) - 1.0) .GT. epsi11) THEN ! there is a bug
+            IF( ABS( asum(ji,jj) - 1._wp )  >  epsi11 ) THEN   ! there is a bug
                WRITE(numout,*) ' '
                WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj)
@@ -391,5 +389,5 @@
       d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:)
       d_smv_i_trp(:,:,:)   = 0._wp
-      IF(  num_sal == 2  .OR.  num_sal == 4  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:)
+      IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:)
 
       IF(ln_ctl) THEN     ! Control print
@@ -430,13 +428,13 @@
 
       ! update of fields will be made later in lim update
-      u_ice(:,:)    = old_u_ice(:,:)
-      v_ice(:,:)    = old_v_ice(:,:)
-      a_i(:,:,:)    = old_a_i(:,:,:)
-      v_s(:,:,:)    = old_v_s(:,:,:)
-      v_i(:,:,:)    = old_v_i(:,:,:)
-      e_s(:,:,:,:)  = old_e_s(:,:,:,:)
-      e_i(:,:,:,:)  = old_e_i(:,:,:,:)
-      oa_i(:,:,:)   = old_oa_i(:,:,:)
-      IF(  num_sal == 2  .OR.  num_sal == 4  )   smv_i(:,:,:)  = old_smv_i(:,:,:)
+      u_ice(:,:)   = old_u_ice(:,:)
+      v_ice(:,:)   = old_v_ice(:,:)
+      a_i(:,:,:)   = old_a_i(:,:,:)
+      v_s(:,:,:)   = old_v_s(:,:,:)
+      v_i(:,:,:)   = old_v_i(:,:,:)
+      e_s(:,:,:,:) = old_e_s(:,:,:,:)
+      e_i(:,:,:,:) = old_e_i(:,:,:,:)
+      oa_i(:,:,:)  = old_oa_i(:,:,:)
+      IF(  num_sal == 2  )   smv_i(:,:,:) = old_smv_i(:,:,:)
 
       !----------------------------------------------------!
@@ -465,23 +463,23 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. &
-                  ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN
-                  old_v_i(ji,jj,jl)     = d_v_i_trp(ji,jj,jl)
-                  d_v_i_trp(ji,jj,jl)   = 0._wp
-                  old_a_i(ji,jj,jl)     = d_a_i_trp(ji,jj,jl)
-                  d_a_i_trp(ji,jj,jl)   = 0._wp
-                  old_v_s(ji,jj,jl)     = d_v_s_trp(ji,jj,jl)
-                  d_v_s_trp(ji,jj,jl)   = 0._wp
-                  old_e_s(ji,jj,1,jl)   = d_e_s_trp(ji,jj,1,jl)
-                  d_e_s_trp(ji,jj,1,jl) = 0._wp
-                  old_oa_i(ji,jj,jl)    = d_oa_i_trp(ji,jj,jl)
-                  d_oa_i_trp(ji,jj,jl)  = 0._wp
-                  IF(  num_sal == 2  .OR.  num_sal == 4  )   old_smv_i(ji,jj,jl)   = d_smv_i_trp(ji,jj,jl)
-                  d_smv_i_trp(ji,jj,jl) = 0._wp
+               IF(  old_v_i  (ji,jj,jl) < epsi06  .AND. &
+                    d_v_i_trp(ji,jj,jl) > epsi06    ) THEN
+                  old_v_i   (ji,jj,jl)   = d_v_i_trp(ji,jj,jl)
+                  d_v_i_trp (ji,jj,jl)   = 0._wp
+                  old_a_i   (ji,jj,jl)   = d_a_i_trp(ji,jj,jl)
+                  d_a_i_trp (ji,jj,jl)   = 0._wp
+                  old_v_s   (ji,jj,jl)   = d_v_s_trp(ji,jj,jl)
+                  d_v_s_trp (ji,jj,jl)   = 0._wp
+                  old_e_s   (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl)
+                  d_e_s_trp (ji,jj,1,jl) = 0._wp
+                  old_oa_i  (ji,jj,jl)   = d_oa_i_trp(ji,jj,jl)
+                  d_oa_i_trp(ji,jj,jl)   = 0._wp
+                  IF(  num_sal == 2  )   old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl)
+                  d_smv_i_trp(ji,jj,jl)  = 0._wp
                ENDIF
             END DO
          END DO
       END DO
-
+      !
       CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final )
       !
@@ -612,13 +610,13 @@
                   ! present
                   zworka(ji,jj) = 4.0 * strength(ji,jj)              &
-                     + strength(ji-1,jj) * tms(ji-1,jj) &  
-                     + strength(ji+1,jj) * tms(ji+1,jj) &  
-                     + strength(ji,jj-1) * tms(ji,jj-1) &  
-                     + strength(ji,jj+1) * tms(ji,jj+1)    
+                     &          + strength(ji-1,jj) * tms(ji-1,jj) &  
+                     &          + strength(ji+1,jj) * tms(ji+1,jj) &  
+                     &          + strength(ji,jj-1) * tms(ji,jj-1) &  
+                     &          + strength(ji,jj+1) * tms(ji,jj+1)    
 
                   zw1 = 4.0 + tms(ji-1,jj) + tms(ji+1,jj) + tms(ji,jj-1) + tms(ji,jj+1)
                   zworka(ji,jj) = zworka(ji,jj) / zw1
                ELSE
-                  zworka(ji,jj) = 0.0
+                  zworka(ji,jj) = 0._wp
                ENDIF
             END DO
@@ -1048,6 +1046,6 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF (aicen_init(ji,jj,jl1) .GT. epsi11 .AND. athorn(ji,jj,jl1) .GT. 0.0       &
-                  .AND. closing_gross(ji,jj) > 0.0) THEN
+               IF( aicen_init(ji,jj,jl1)  >  epsi11  .AND.  athorn(ji,jj,jl1)  >  0._wp       &
+                  .AND. closing_gross(ji,jj) > 0._wp ) THEN
                   icells = icells + 1
                   indxi(icells) = ji
@@ -1130,5 +1128,5 @@
             ! Salinity
             !-------------
-            smsw(ji,jj)  = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0       ! salt content of water frozen in voids
+            smsw(ji,jj)  = sss_m(ji,jj) * vsw(ji,jj) * rhoic / rau0       ! salt content of seawater frozen in voids
 
             zsrdg2       = srdg1(ji,jj) + smsw(ji,jj)                     ! salt content of new ridge
@@ -1137,6 +1135,8 @@
             
             !                                                             ! excess of salt is flushed into the ocean
-            fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic / rdt_ice
-
+            fsalt_rpo(ji,jj) = fsalt_rpo(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice
+
+            rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic / rau0   ! increase in ice volume du to seawater frozen in voids
+            
             !------------------------------------            
             ! 3.6 Increment ridging diagnostics
@@ -1148,6 +1148,6 @@
             dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj)
             dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj)
-            diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) / rdt_ice
-            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj)*rdt_ice
+            diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice
+            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice
 
             IF( con_i )   vice_init(ji,jj) = vice_init(ji,jj) + vrdg2(ji,jj) - vrdg1(ji,jj)
@@ -1156,5 +1156,4 @@
             ! 3.7 Put the snow somewhere in the ocean
             !------------------------------------------            
-
             !  Place part of the snow lost by ridging into the ocean. 
             !  Note that esnow_mlt < 0; the ocean must cool to melt snow.
@@ -1179,7 +1178,6 @@
             !           ij looping 1-icells
 
-            dhr(ji,jj)  = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)
+            dhr (ji,jj) = hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1)
             dhr2(ji,jj) = hrmax(ji,jj,jl1) * hrmax(ji,jj,jl1) - hrmin(ji,jj,jl1) * hrmin(ji,jj,jl1)
-
 
          END DO                 ! ij
@@ -1211,5 +1209,5 @@
 
                ! heat flux
-               fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) / rdt_ice
+               fheat_rpo(ji,jj) = fheat_rpo(ji,jj) + ( ersw(ji,jj,jk) - zdummy0 ) * r1_rdtice
 
                ! Correct dimensions to avoid big values
@@ -1275,23 +1273,23 @@
                ! Transfer area, volume, and energy accordingly.
 
-               IF (hrmin(ji,jj,jl1) .GE. hi_max(jl2) .OR.        &
-                  hrmax(ji,jj,jl1) .LE. hi_max(jl2-1)) THEN
-                  hL = 0.0
-                  hR = 0.0
+               IF( hrmin(ji,jj,jl1) >= hi_max(jl2)  .OR.        &
+                   hrmax(ji,jj,jl1) <= hi_max(jl2-1) ) THEN
+                  hL = 0._wp
+                  hR = 0._wp
                ELSE
-                  hL = MAX (hrmin(ji,jj,jl1), hi_max(jl2-1))
-                  hR = MIN (hrmax(ji,jj,jl1), hi_max(jl2))
+                  hL = MAX( hrmin(ji,jj,jl1), hi_max(jl2-1) )
+                  hR = MIN( hrmax(ji,jj,jl1), hi_max(jl2)   )
                ENDIF
 
                ! fraction of ridged ice area and volume going to n2
-               farea = (hR-hL) / dhr(ji,jj) 
-               fvol(ji,jj) = (hR*hR - hL*hL) / dhr2(ji,jj)
-
-               a_i  (ji,jj,jl2)   = a_i  (ji,jj,jl2)   + ardg2 (ji,jj) * farea
-               v_i  (ji,jj,jl2)   = v_i  (ji,jj,jl2)   + vrdg2 (ji,jj) * fvol(ji,jj)
-               v_s  (ji,jj,jl2)   = v_s  (ji,jj,jl2)   + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg
+               farea = ( hR - hL ) / dhr(ji,jj) 
+               fvol(ji,jj) = ( hR*hR - hL*hL ) / dhr2(ji,jj)
+
+               a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + ardg2 (ji,jj) * farea
+               v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + vrdg2 (ji,jj) * fvol(ji,jj)
+               v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrdg (ji,jj) * fvol(ji,jj) * fsnowrdg
                e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrdg (ji,jj) * fvol(ji,jj) * fsnowrdg
-               smv_i(ji,jj,jl2)   = smv_i(ji,jj,jl2)   + srdg2 (ji,jj) * fvol(ji,jj)
-               oa_i (ji,jj,jl2)   = oa_i (ji,jj,jl2)   + oirdg2(ji,jj) * farea
+               smv_i(ji,jj  ,jl2) = smv_i(ji,jj  ,jl2) + srdg2 (ji,jj) * fvol(ji,jj)
+               oa_i (ji,jj  ,jl2) = oa_i (ji,jj  ,jl2) + oirdg2(ji,jj) * farea
 
             END DO ! ij
@@ -1317,15 +1315,15 @@
                ! Compute the fraction of rafted ice area and volume going to 
                ! thickness category jl2, transfer area, volume, and energy accordingly.
-
-               IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND.        &
-                  hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN
-                  a_i(ji,jj,jl2) = a_i(ji,jj,jl2) + arft2(ji,jj)
-                  v_i(ji,jj,jl2) = v_i(ji,jj,jl2) + virft(ji,jj)
-                  v_s(ji,jj,jl2) = v_s(ji,jj,jl2) + vsrft(ji,jj)*fsnowrft
-                  e_s(ji,jj,1,jl2) = e_s(ji,jj,1,jl2) + esrft(ji,jj)*fsnowrft
-                  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)    
+               !
+               IF( hraft(ji,jj,jl1) <= hi_max(jl2)  .AND.        &
+                   hraft(ji,jj,jl1) >  hi_max(jl2-1) ) THEN
+                  a_i  (ji,jj  ,jl2) = a_i  (ji,jj  ,jl2) + arft2 (ji,jj)
+                  v_i  (ji,jj  ,jl2) = v_i  (ji,jj  ,jl2) + virft (ji,jj)
+                  v_s  (ji,jj  ,jl2) = v_s  (ji,jj  ,jl2) + vsrft (ji,jj) * fsnowrft
+                  e_s  (ji,jj,1,jl2) = e_s  (ji,jj,1,jl2) + esrft (ji,jj) * fsnowrft
+                  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
-
+               !
             END DO ! ij
 
@@ -1336,6 +1334,6 @@
                   ji = indxi(ij)
                   jj = indxj(ij)
-                  IF (hraft(ji,jj,jl1) .LE. hi_max(jl2) .AND.        &
-                     hraft(ji,jj,jl1) .GT. hi_max(jl2-1)) THEN
+                  IF(  hraft(ji,jj,jl1)  <=  hi_max(jl2)   .AND.        &
+                       hraft(ji,jj,jl1)  >   hi_max(jl2-1)  ) THEN
                      e_i(ji,jj,jk,jl2) = e_i(ji,jj,jk,jl2) + eirft(ji,jj,jk)
                   ENDIF
@@ -1504,5 +1502,5 @@
             DO jj = 1 , jpj
                DO ji = 1 , jpi
-!!gm                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) / rdt_ice
+!!gm                  xtmp = e_i(ji,jj,jk,jl) / area(ji,jj) * r1_rdtice
 !!gm                  xtmp = xtmp * unit_fac
                   ! fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
@@ -1524,8 +1522,8 @@
                ! fluxes are positive to the ocean
                ! here the flux has to be negative for the ocean
-!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice
+!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice
                !           fheat_res(ji,jj) = fheat_res(ji,jj) - xtmp
 
-!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) / rdt_ice !RB   ???????
+!!gm               xtmp = ( rhosn*cpic*( rtt-t_s(ji,jj,1,jl) ) + rhosn*lfus ) * r1_rdtice !RB   ???????
 
                t_s(ji,jj,1,jl) = rtt * zmask(ji,jj) + t_s(ji,jj,1,jl) * ( 1 - zmask(ji,jj) )
@@ -1540,8 +1538,8 @@
 
                !           fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj)                  ) * & 
-               !                               rhosn * v_s(ji,jj,jl) / rdt_ice
+               !                               rhosn * v_s(ji,jj,jl) * r1_rdtice
 
                !           fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * & 
-               !                               rhoic * v_i(ji,jj,jl) / rdt_ice
+               !                               rhoic * v_i(ji,jj,jl) * r1_rdtice
 
                !           sfx (i,j)      = sfx (i,j)      + xtmp
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90	(revision 3517)
@@ -101,12 +101,12 @@
 
       !- Trend terms
-      d_a_i_thd (:,:,:)  = a_i(:,:,:)   - old_a_i(:,:,:) 
-      d_v_s_thd (:,:,:)  = v_s(:,:,:)   - old_v_s(:,:,:)
-      d_v_i_thd (:,:,:)  = v_i(:,:,:)   - old_v_i(:,:,:)  
+      d_a_i_thd(:,:,:)   = a_i(:,:,:)   - old_a_i(:,:,:) 
+      d_v_s_thd(:,:,:)   = v_s(:,:,:)   - old_v_s(:,:,:)
+      d_v_i_thd(:,:,:)   = v_i(:,:,:)   - old_v_i(:,:,:)  
       d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:) 
       d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:)
 
       d_smv_i_thd(:,:,:) = 0._wp
-      IF( num_sal == 2 .OR. num_sal == 4 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
+      IF( num_sal == 2  )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
 
       IF(ln_ctl) THEN   ! Control print
@@ -143,11 +143,11 @@
 
       !- Recover Old values
-      a_i(:,:,:)   = old_a_i (:,:,:)
-      v_s(:,:,:)   = old_v_s (:,:,:)
-      v_i(:,:,:)   = old_v_i (:,:,:)
-      e_s(:,:,:,:) = old_e_s (:,:,:,:)
-      e_i(:,:,:,:) = old_e_i (:,:,:,:)
-      !
-      IF( num_sal == 2 .OR. num_sal == 4 )   smv_i(:,:,:)       = old_smv_i (:,:,:)
+      a_i(:,:,:)   = old_a_i(:,:,:)
+      v_s(:,:,:)   = old_v_s(:,:,:)
+      v_i(:,:,:)   = old_v_i(:,:,:)
+      e_s(:,:,:,:) = old_e_s(:,:,:,:)
+      e_i(:,:,:,:) = old_e_i(:,:,:,:)
+      !
+      IF( num_sal == 2 )   smv_i(:,:,:) = old_smv_i(:,:,:)
       !
    END SUBROUTINE lim_itd_th
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 3517)
@@ -402,7 +402,7 @@
                      zsmax = 4.5_wp
                      zsmin = 3.5_wp
-                     IF( sm_i(ji,jj,jl) .LT. zsmin ) THEN
+                     IF(     sm_i(ji,jj,jl) < zsmin ) THEN
                         zalpha = 1._wp
-                     ELSEIF( sm_i(ji,jj,jl) .LT.zsmax ) THEN
+                     ELSEIF( sm_i(ji,jj,jl) < zsmax ) THEN
                         zalpha = sm_i(ji,jj,jl) / ( zsmin - zsmax ) + zsmax / ( zsmax - zsmin )
                      ELSE
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 3517)
@@ -9,5 +9,6 @@
    !!            3.3  ! 2010-05 (G. Madec) decrease ocean & ice reference salinities in the Baltic sea
    !!                 !                  + simplification of the ice-ocean stress calculation
-   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
+   !!            3.5  ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -44,5 +45,4 @@
    PUBLIC   lim_sbc_tau    ! called by sbc_ice_lim
 
-   REAL(wp)  ::   r1_rdtice            ! = 1. / rdt_ice 
    REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values
    REAL(wp)  ::   rzero  = 0._wp    
@@ -88,5 +88,5 @@
       !!              - qns     : sea heat flux: non solar
       !!              - emp     : freshwater budget: volume flux 
-      !!              - sfx     : freshwater budget: concentration/dillution 
+      !!              - sfx     : salt flux 
       !!              - fr_i    : ice fraction
       !!              - tn_ice  : sea-ice surface temperature
@@ -102,7 +102,7 @@
       INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
       INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
-      REAL(wp) ::   zinda, zfons, zpme              ! local scalars
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes
-      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp   ! 2D/3D workspace
+      REAL(wp) ::   zinda, zemp, zemp_snow, zfmm    ! local scalars
+      REAL(wp), POINTER, DIMENSION(:,:)   ::   zfcm1 , zfcm2   ! solar/non solar heat fluxes
+      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace
       !!---------------------------------------------------------------------
       
@@ -153,10 +153,9 @@
                !              &           + iflt * ffltbif(ji,jj) !!! only if one category is used
                &           + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   &
-               &           + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice   &
-               &           + fhmec(ji,jj)     & ! new contribution due to snow melt in ridging!!
+               &           + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice   &
+               &           + fhmec(ji,jj)     & ! new contribution due to snow melt when ridging!!
                &           + fheat_rpo(ji,jj) & ! contribution from ridge formation
                &           + fheat_res(ji,jj)
-            ! fscmbq  Part of the solar radiation transmitted through the ice and going to the ocean
-            !         computed in limthd_zdf.F90
+            ! fscmbq  Part of the solar radiation transmitted through the ice and going to the ocean computed in limthd_zdf.F90
             ! ffltbif Total heat content of the ice (brine pockets+ice) / delta_t
             ! qcmif   Energy needed to bring the ocean surface layer until its freezing (ok)
@@ -167,9 +166,7 @@
             ! qdtcn Energy from the turbulent oceanic heat flux heat flux coming in the lead
 
-            IF ( num_sal == 2 ) zfcm2(ji,jj) = zfcm2(ji,jj) + &
-               fhbri(ji,jj) ! new contribution due to brine drainage 
-
-            ! bottom radiative component is sent to the computation of the
-            ! oceanic heat flux
+            IF( num_sal == 2 )   zfcm2(ji,jj) = zfcm2(ji,jj) + fhbri(ji,jj)    ! add contribution due to brine drainage 
+
+            ! bottom radiative component is sent to the computation of the oceanic heat flux
             fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     
 
@@ -179,5 +176,5 @@
             !                           ! fdtcn : turbulent oceanic heat flux
 
-            !!gm   this IF prevents the vertorisation of the whole loop
+!!gm   this IF prevents the vertorisation of the whole loop
             IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN
                WRITE(numout,*) ' lim_sbc : heat fluxes '
@@ -208,5 +205,5 @@
                WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)
             ENDIF
-            !!gm   end
+!!gm   end
          END DO
       END DO
@@ -229,43 +226,29 @@
 
             !  computing freshwater exchanges at the ice/ocean interface
-            zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction
-               &   + tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean
-               &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice
-               &   - rdm_snw(ji,jj) * r1_rdtice                       &   ! freshwaterflux due to snow melting 
-               &   + fmmec(ji,jj)                                         ! snow falling when ridging
-
-
-            !  computing salt exchanges at the ice/ocean interface
-            !  sice should be the same as computed with the ice model
-            zfons =  ( soce_0(ji,jj) - sice_0(ji,jj) ) * rdm_ice(ji,jj) * r1_rdtice 
-            ! SOCE
-            zfons =  ( sss_m (ji,jj) - sice_0(ji,jj) ) * rdm_ice(ji,jj) * r1_rdtice
-
-            !CT useless            !  salt flux for constant salinity
-            !CT useless            fsalt(ji,jj)      =  zfons / ( sss_m(ji,jj) + epsi16 ) + fsalt_res(ji,jj)
-            !  salt flux for variable salinity
-            zinda             = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
-            !  correcting brine and salt fluxes
-            fsbri(ji,jj)      =  zinda*fsbri(ji,jj)
-            !  converting the salt fluxes from ice to a freshwater flux from ocean
-            fsalt_res(ji,jj)  =  fsalt_res(ji,jj) / ( sss_m(ji,jj) + epsi16 )
-            fseqv(ji,jj)      =  fseqv(ji,jj)     / ( sss_m(ji,jj) + epsi16 )
-            fsbri(ji,jj)      =  fsbri(ji,jj)     / ( sss_m(ji,jj) + epsi16 )
-            fsalt_rpo(ji,jj)  =  fsalt_rpo(ji,jj) / ( sss_m(ji,jj) + epsi16 )
-
-            !  freshwater mass exchange (positive to the ice, negative for the ocean ?)
-            !  actually it's a salt flux (so it's minus freshwater flux)
-            !  if sea ice grows, zfons is positive, fsalt also
-            !  POSITIVE SALT FLUX FROM THE ICE TO THE OCEAN
-            !  POSITIVE FRESHWATER FLUX FROM THE OCEAN TO THE ICE [kg.m-2.s-1]
-
-            emp(ji,jj) = - zpme 
+            zemp =   emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction
+               &   - tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean
+               &   + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice
+               &   - fmmec(ji,jj)                                         ! snow falling when ridging
+
+            ! mass flux at the ocean/ice interface (sea ice fraction)
+            zemp_snw = rdm_snw(ji,jj) * r1_rdtice                         ! snow melting = pure water that enters the ocean
+            zfmm     = rdm_ice(ji,jj) * r1_rdtice                         ! Freezing minus mesting  
+
+            emp(ji,jj) = zemp + zemp_snw + zfmm  ! mass flux + F/M mass flux (always ice/ocean mass exchange)
+            
+            !  correcting brine salt fluxes   (zinda = 1  if pfrld=1 , =0 otherwise)
+            zinda        = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
+            fsbri(ji,jj) = zinda * fsbri(ji,jj)
          END DO
       END DO
 
+      !------------------------------------------!
+      !      salt flux at the ocean surface      !
+      !------------------------------------------!
+
       IF( num_sal == 2 ) THEN      ! variable ice salinity: brine drainage included in the salt flux
-         sfx (:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
+         sfx (:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + fsbri(:,:)
       ELSE                         ! constant ice salinity:
-         sfx (:,:) =              fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
+         sfx (:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:)
       ENDIF
       !-----------------------------------------------!
@@ -277,5 +260,5 @@
          snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  ) 
          !                                                      ! time evolution of snow+ice mass
-         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) / rdt_ice
+         snwice_fmass (:,:) = ( snwice_mass(:,:) - snwice_mass_b(:,:) ) * r1_rdtice
       ENDIF
 
@@ -403,6 +386,4 @@
       !                                      ! allocate lim_sbc array
       IF( lim_sbc_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'lim_sbc_init : unable to allocate standard arrays' )
-      !
-      r1_rdtice = 1. / rdt_ice
       !
       soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating sea-ice case
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 3517)
@@ -110,5 +110,5 @@
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * nlay_i
                   !0 if no ice and 1 if yes
-                  zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_i(ji,jj,jl) ) ) 
+                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) )  ) 
                   !convert units ! very important that this line is here
                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 
@@ -122,5 +122,5 @@
                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * nlay_s
                   !0 if no ice and 1 if yes
-                  zindb = 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s(ji,jj,jl) ) ) 
+                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) )  ) 
                   !convert units ! very important that this line is here
                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * unit_fac * zindb 
@@ -284,6 +284,6 @@
             CALL tab_2d_1d( nbpb, fbif_1d    (1:nbpb), fbif       , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, qldif_1d   (1:nbpb), qldif      , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, rdmicif_1d (1:nbpb), rdm_ice    , jpi, jpj, npb(1:nbpb) )
-            CALL tab_2d_1d( nbpb, rdmsnif_1d (1:nbpb), rdm_snw    , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, rdm_ice_1d (1:nbpb), rdm_ice    , jpi, jpj, npb(1:nbpb) )
+            CALL tab_2d_1d( nbpb, rdm_snw_1d (1:nbpb), rdm_snw    , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, dmgwi_1d   (1:nbpb), dmgwi      , jpi, jpj, npb(1:nbpb) )
             CALL tab_2d_1d( nbpb, qlbbq_1d   (1:nbpb), zqlbsbq    , jpi, jpj, npb(1:nbpb) )
@@ -352,6 +352,6 @@
             CALL tab_1d_2d( nbpb, qldif  , npb, qldif_1d  (1:nbpb), jpi, jpj )
             CALL tab_1d_2d( nbpb, qfvbq  , npb, qfvbq_1d  (1:nbpb), jpi, jpj )
-            CALL tab_1d_2d( nbpb, rdm_ice, npb, rdmicif_1d(1:nbpb), jpi, jpj )
-            CALL tab_1d_2d( nbpb, rdm_snw, npb, rdmsnif_1d(1:nbpb), jpi, jpj )
+            CALL tab_1d_2d( nbpb, rdm_ice, npb, rdm_ice_1d(1:nbpb), jpi, jpj )
+            CALL tab_1d_2d( nbpb, rdm_snw, npb, rdm_snw_1d(1:nbpb), jpi, jpj )
             CALL tab_1d_2d( nbpb, dmgwi  , npb, dmgwi_1d  (1:nbpb), jpi, jpj )
             CALL tab_1d_2d( nbpb, rdvosif, npb, dvsbq_1d  (1:nbpb), jpi, jpj )
@@ -419,5 +419,5 @@
       !--------------------------------------------
       d_v_i_thd(:,:,:) = v_i      (:,:,:) - old_v_i(:,:,:)    ! ice volumes 
-      dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) / rdt_ice * 86400.0
+      dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday
 
       IF( con_i )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:)
@@ -488,7 +488,7 @@
       !
       IF(lwp) WRITE(numout,*) ' lim_thd_glohec '
-      IF(lwp) WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) / rdt_ice
-      IF(lwp) WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) / rdt_ice
-      IF(lwp) WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) / rdt_ice
+      IF(lwp) WRITE(numout,*) ' qt_i_in : ', eti(jiindex_1d,jl) * r1_rdtice
+      IF(lwp) WRITE(numout,*) ' qt_s_in : ', ets(jiindex_1d,jl) * r1_rdtice
+      IF(lwp) WRITE(numout,*) ' qt_in   : ', ( eti(jiindex_1d,jl) + ets(jiindex_1d,jl) ) * r1_rdtice
       !
    END SUBROUTINE lim_thd_glohec
@@ -538,5 +538,5 @@
       !--------------------
       DO ji = kideb, kiut
-         cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
+         cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) )
       END DO
 
@@ -597,5 +597,5 @@
             WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
             WRITE(numout,*) ' surf_error : ', surf_error(ji,jl)
-            WRITE(numout,*) ' dq_i       : ', - dq_i(ji,jl) / rdt_ice
+            WRITE(numout,*) ' dq_i       : ', - dq_i(ji,jl) * r1_rdtice
             WRITE(numout,*) ' Fdt        : ', sum_fluxq(ji,jl)
             WRITE(numout,*)
@@ -631,9 +631,9 @@
             WRITE(numout,*)
             WRITE(numout,*) ' Layer by layer ... '
-            WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice
+            WRITE(numout,*) ' dq_snow : ', ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) * r1_rdtice
             WRITE(numout,*) ' dfc_snow  : ', fc_s(ji,1) - fc_s(ji,0)
             DO jk = 1, nlay_i
                WRITE(numout,*) ' layer  : ', jk
-               WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) / rdt_ice  
+               WRITE(numout,*) ' dq_ice : ', dq_i_layer(ji,jk) * r1_rdtice  
                WRITE(numout,*) ' radab  : ', radab(ji,jk)
                WRITE(numout,*) ' dfc_i  : ', fc_i(ji,jk) - fc_i(ji,jk-1)
@@ -681,5 +681,5 @@
          fatm      (ji,jl) = qnsr_ice_1d(ji) + qsr_ice_1d(ji)                       ! total heat flux
          sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(zji,zjj,jl) 
-         cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
+         cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) )
       END DO
 
@@ -688,5 +688,5 @@
       !--------------------
       DO ji = kideb, kiut
-         cons_error(ji,jl) = ABS( dq_i(ji,jl) / rdt_ice + sum_fluxq(ji,jl) )
+         cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) )
       END DO
 
@@ -722,7 +722,7 @@
             WRITE(numout,*) ' * '
             WRITE(numout,*) ' Ftotal     : ', sum_fluxq(ji,jl)
-            WRITE(numout,*) ' dq_t       : ', - dq_i(ji,jl) / rdt_ice
-            WRITE(numout,*) ' dq_i       : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) / rdt_ice
-            WRITE(numout,*) ' dq_s       : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / rdt_ice
+            WRITE(numout,*) ' dq_t       : ', - dq_i(ji,jl) * r1_rdtice
+            WRITE(numout,*) ' dq_i       : ', - ( qt_i_fin(ji,jl) - qt_i_in(ji,jl) ) * r1_rdtice
+            WRITE(numout,*) ' dq_s       : ', - ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) * r1_rdtice
             WRITE(numout,*) ' cons_error : ', cons_error(ji,jl)
             WRITE(numout,*) ' * '
@@ -734,10 +734,10 @@
             WRITE(numout,*) ' * '
             WRITE(numout,*) ' Heat contents --- : '
-            WRITE(numout,*) ' qt_s_in    : ', qt_s_in(ji,jl) / rdt_ice
-            WRITE(numout,*) ' qt_i_in    : ', qt_i_in(ji,jl) / rdt_ice
-            WRITE(numout,*) ' qt_in      : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) / rdt_ice
-            WRITE(numout,*) ' qt_s_fin   : ', qt_s_fin(ji,jl) / rdt_ice
-            WRITE(numout,*) ' qt_i_fin   : ', qt_i_fin(ji,jl) / rdt_ice
-            WRITE(numout,*) ' qt_fin     : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) / rdt_ice
+            WRITE(numout,*) ' qt_s_in    : ', qt_s_in(ji,jl) * r1_rdtice
+            WRITE(numout,*) ' qt_i_in    : ', qt_i_in(ji,jl) * r1_rdtice
+            WRITE(numout,*) ' qt_in      : ', ( qt_i_in(ji,jl) + qt_s_in(ji,jl) ) * r1_rdtice
+            WRITE(numout,*) ' qt_s_fin   : ', qt_s_fin(ji,jl) * r1_rdtice
+            WRITE(numout,*) ' qt_i_fin   : ', qt_i_fin(ji,jl) * r1_rdtice
+            WRITE(numout,*) ' qt_fin     : ', ( qt_i_fin(ji,jl) + qt_s_fin(ji,jl) ) * r1_rdtice
             WRITE(numout,*) ' * '
             WRITE(numout,*) ' Ice variables --- : '
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 3517)
@@ -7,5 +7,6 @@
    !!                 ! 2005-06 (M. Vancoppenolle) 3D version 
    !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif & rdmicif
-   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation
+   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
+   !!            3.5  ! 2012-10 (G. Madec & co) salt flux + bug fixes 
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -71,5 +72,5 @@
       !! 
       INTEGER  ::   ji , jk        ! dummy loop indices
-      INTEGER  ::   zji, zjj       ! 2D corresponding indices to ji
+      INTEGER  ::   ii, ij       ! 2D corresponding indices to ji
       INTEGER  ::   isnow          ! switch for presence (1) or absence (0) of snow
       INTEGER  ::   isnowic        ! snow ice formation not
@@ -145,8 +146,8 @@
       !
       DO ji = kideb, kiut
-         isnow         = INT( 1.0 - MAX ( 0.0 , SIGN ( 1.0 , - ht_s_b(ji) ) ) )
-         ztfs(ji)      = isnow * rtt + ( 1.0 - isnow ) * rtt
-         z_f_surf(ji)  = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)
-         z_f_surf(ji)  = MAX( zzero , z_f_surf(ji) ) * MAX( zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) ) )
+         isnow         = INT(  1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s_b(ji) )  )  )
+         ztfs     (ji) = isnow * rtt + ( 1.0 - isnow ) * rtt
+         z_f_surf (ji) = qnsr_ice_1d(ji) + ( 1.0 - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)
+         z_f_surf (ji) = MAX(  zzero , z_f_surf(ji)  ) * MAX(  zzero , SIGN( zone , t_su_b(ji) - ztfs(ji) )  )
          zfdt_init(ji) = ( z_f_surf(ji) + MAX( fbif_1d(ji) + qlbbq_1d(ji) + fc_bo_i(ji),0.0 ) ) * rdt_ice
       END DO ! ji
@@ -240,10 +241,10 @@
          zhsnew         =  ht_s_b(ji) + dh_s_tot(ji)
          ! If snow is still present zhn = 1, else zhn = 0
-         zhn            =  1.0 - MAX( zzero , SIGN( zone , - zhsnew ) )
+         zhn            =  1.0 - MAX(  zzero , SIGN( zone , - zhsnew )  )
          ht_s_b(ji)     =  MAX( zzero , zhsnew )
          ! Volume and mass variations of snow
-         dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_mel(ji) )
+         dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) )
          dvsbq_1d  (ji) =  MIN( zzero, dvsbq_1d(ji) )
-         rdmsnif_1d(ji) =  rdmsnif_1d(ji) + rhosn * dvsbq_1d(ji)
+         rdm_snw_1d(ji) =  rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji)
       END DO ! ji
 
@@ -253,5 +254,5 @@
       DO ji = kideb, kiut 
          dh_i_surf(ji) =  0._wp
-         z_f_surf (ji) =  zqfont_su(ji) / rdt_ice ! heat conservation test
+         z_f_surf (ji) =  zqfont_su(ji) * r1_rdtice   ! heat conservation test
          zdq_i    (ji) =  0._wp
       END DO ! ji
@@ -262,17 +263,14 @@
             zdeltah  (ji,jk) = - zqfont_su(ji) / q_i_b(ji,jk)
             !                                                    ! recompute heat available
-            zqfont_su(ji)    = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 
+            zqfont_su(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk) 
             !                                                    ! melt of layer jk cannot be higher than its thickness
             zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_i(ji) )
             !                                                    ! update surface melt
-            dh_i_surf(ji)    = dh_i_surf(ji) + zdeltah(ji,jk) 
+            dh_i_surf(ji   ) = dh_i_surf(ji) + zdeltah(ji,jk) 
             !                                                    ! for energy conservation
-            zdq_i    (ji)    = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice
+            zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice
             !
-            ! contribution to ice-ocean salt flux 
-            zji = MOD( npb(ji) - 1 , jpi ) + 1
-            zjj =    ( npb(ji) - 1 ) / jpi + 1
-            zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji) ) * a_i_b(ji)    &
-               &                              * MIN( zdeltah(ji,jk) , 0.e0 ) * rhoic / rdt_ice 
+            !                                                    ! contribution to ice-ocean salt flux 
+            zfsalt_melt(ji)  = zfsalt_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 
          END DO
       END DO
@@ -290,5 +288,5 @@
             IF( z_f_surf(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN!
                WRITE(numout,*) ' ALERTE heat loss for surface melt '
-               WRITE(numout,*) ' zji, zjj, jl :', zji, zjj, jl
+               WRITE(numout,*) ' ii, ij, jl :', ii, ij, jl
                WRITE(numout,*) ' ht_i_b       : ', ht_i_b(ji)
                WRITE(numout,*) ' z_f_surf     : ', z_f_surf(ji)
@@ -299,5 +297,5 @@
                WRITE(numout,*) ' qlbbq_1d     : ', qlbbq_1d(ji)
                WRITE(numout,*) ' s_i_new      : ', s_i_new(ji)
-               WRITE(numout,*) ' sss_m        : ', sss_m(zji,zjj)
+               WRITE(numout,*) ' sss_m        : ', sss_m(ii,ij)
             ENDIF
          END DO
@@ -338,5 +336,5 @@
          DO ji = kideb, kiut
             ! In case of disparition of the snow, we have to update the snow temperatures
-            zhisn  =  MAX( zzero , SIGN( zone, - ht_s_b(ji) ) )
+            zhisn  =  MAX(  zzero , SIGN( zone, - ht_s_b(ji) )  )
             t_s_b(ji,jk) = ( 1.0 - zhisn ) * t_s_b(ji,jk) + zhisn * rtt
             q_s_b(ji,jk) = ( 1.0 - zhisn ) * q_s_b(ji,jk)
@@ -358,7 +356,7 @@
       ! 4.1 Basal growth - (a) salinity not varying in time 
       !-----------------------------------------------------
-      IF(  num_sal /= 2  .AND.  num_sal /= 4  ) THEN
-         DO ji = kideb, kiut
-            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0.0  ) THEN
+      IF(  num_sal /= 2  ) THEN   ! ice salinity constant in time
+         DO ji = kideb, kiut
+            IF(  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) < 0._wp  ) THEN
                s_i_new(ji)         =  sm_i_b(ji)
                ! Melting point in K
@@ -371,5 +369,5 @@
                   &                           - rcp  * ( ztmelts - rtt )                                 )
                ! Basal growth rate = - F*dt / q
-               dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
+               dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
             ENDIF
          END DO
@@ -379,9 +377,8 @@
       ! 4.1 Basal growth - (b) salinity varying in time 
       !-------------------------------------------------
-      IF(  num_sal == 2 .OR.  num_sal == 4  ) THEN
-         ! the growth rate (dh_i_bott) is function of the new ice
-         ! heat content (q_i_b(nlay_i+1)). q_i_b depends on the new ice 
-         ! salinity (snewice). snewice depends on dh_i_bott
-         ! it converges quickly, so, no problem
+      IF(  num_sal == 2  ) THEN
+         ! the growth rate (dh_i_bott) is function of the new ice heat content (q_i_b(nlay_i+1)). 
+         ! q_i_b depends on the new ice salinity (snewice). 
+         ! snewice depends on dh_i_bott ; it converges quickly, so, no problem
          ! See Vancoppenolle et al., OM08 for more info on this
 
@@ -394,6 +391,6 @@
             DO ji = kideb, kiut
                IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) < 0.e0  ) THEN
-                  zji = MOD( npb(ji) - 1, jpi ) + 1
-                  zjj = ( npb(ji) - 1 ) / jpi + 1
+                  ii = MOD( npb(ji) - 1, jpi ) + 1
+                  ij = ( npb(ji) - 1 ) / jpi + 1
                   ! Melting point in K
                   ztmelts             =   - tmut * s_i_new(ji) + rtt 
@@ -408,5 +405,5 @@
                   ! zswi12 (1) if dh_i_bott/rdt .LT. 3.6e-7 and .GT. 2.0e-8
                   ! zswi1  (1) if dh_i_bott/rdt .LT. 2.0e-8
-                  zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) / rdt_ice , epsi13 ) )
+                  zgrr   = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi13 ) )
                   zswi2  = MAX( zzero , SIGN( zone , zgrr - 3.6e-7 ) ) 
                   zswi12 = MAX( zzero , SIGN( zone , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
@@ -414,6 +411,6 @@
                   zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
                      &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 
-                  zds         = zfracs * sss_m(zji,zjj) - s_i_new(ji)
-                  s_i_new(ji) = zfracs * sss_m(zji,zjj)
+                  zds         = zfracs * sss_m(ii,ij) - s_i_new(ji)
+                  s_i_new(ji) = zfracs * sss_m(ii,ij)
                ENDIF ! fc_bo_i
             END DO ! ji
@@ -432,5 +429,5 @@
                   &                            - rcp * ( ztmelts - rtt )                                    )
                ! Basal growth rate = - F*dt / q
-               dh_i_bott(ji)       =  - rdt_ice*( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1) 
+               dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)
                ! Salinity update
                ! entrapment during bottom growth
@@ -453,5 +450,5 @@
             s_i_new(ji)   =  s_i_b(ji,nlay_i)
             zqfont_bo(ji) =  rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) )
-            zfbase(ji)    =  zqfont_bo(ji) / rdt_ice     ! heat conservation test
+            zfbase(ji)    =  zqfont_bo(ji) * r1_rdtice     ! heat conservation test
             zdq_i(ji)     =  0._wp
             dh_i_bott(ji) =  0._wp
@@ -461,22 +458,19 @@
       DO jk = nlay_i, 1, -1
          DO ji = kideb, kiut
-            IF (  ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .GE. 0.0  ) THEN
-               ztmelts            =   - tmut * s_i_b(ji,jk) + rtt 
-               IF( t_i_b(ji,jk) >= ztmelts ) THEN
-                  zdeltah(ji,jk)  = - zh_i(ji)
-                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
-                  zinnermelt(ji)   = 1._wp
-               ELSE  ! normal ablation
-                  zdeltah(ji,jk)  = - zqfont_bo(ji) / q_i_b(ji,jk)
-                  zqfont_bo(ji)   = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)
-                  zdeltah(ji,jk)  = MAX(zdeltah(ji,jk), - zh_i(ji) )
-                  dh_i_bott(ji)   = dh_i_bott(ji) + zdeltah(ji,jk)
-                  zdq_i(ji)       = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) / rdt_ice
-                  ! contribution to salt flux
-                  zji             = MOD( npb(ji) - 1, jpi ) + 1
-                  zjj             = ( npb(ji) - 1 ) / jpi + 1
-                  zfsalt_melt(ji) = zfsalt_melt(ji) + ( sss_m(zji,zjj) - sm_i_b(ji)   ) * a_i_b(ji)   &
-                     &                              * MIN( zdeltah(ji,jk) , 0.0 ) * rhoic / rdt_ice 
+            IF(  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  >=  0._wp  ) THEN
+               ztmelts = - tmut * s_i_b(ji,jk) + rtt 
+               IF( t_i_b(ji,jk) >= ztmelts ) THEN   !!gm : a comment is needed
+                  zdeltah   (ji,jk) = - zh_i(ji)
+                  dh_i_bott (ji   ) = dh_i_bott(ji) + zdeltah(ji,jk)
+                  zinnermelt(ji   ) = 1._wp
+               ELSE                                  ! normal ablation
+                  zdeltah  (ji,jk) = - zqfont_bo(ji) / q_i_b(ji,jk)
+                  zqfont_bo(ji   ) = MAX( 0.0 , - zh_i(ji) - zdeltah(ji,jk) ) * q_i_b(ji,jk)
+                  zdeltah  (ji,jk) = MAX(zdeltah(ji,jk), - zh_i(ji) )
+                  dh_i_bott(ji   ) = dh_i_bott(ji) + zdeltah(ji,jk)
+                  zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice
                ENDIF
+               ! contribution to salt flux
+               zfsalt_melt(ji) = zfsalt_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 
             ENDIF
          END DO ! ji
@@ -493,5 +487,5 @@
                ENDIF
                IF ( zfbase(ji) + zdq_i(ji) .GE. 1.0e-3  ) THEN
-                  WRITE(numout,*) ' ALERTE heat loss for basal melt : zji, zjj, jl :', zji, zjj, jl
+                  WRITE(numout,*) ' ALERTE heat loss for basal melt : ii, ij, jl :', ii, ij, jl
                   WRITE(numout,*) ' ht_i_b    : ', ht_i_b(ji)
                   WRITE(numout,*) ' zfbase    : ', zfbase(ji)
@@ -502,5 +496,5 @@
                   WRITE(numout,*) ' qlbbq_1d  : ', qlbbq_1d(ji)
                   WRITE(numout,*) ' s_i_new   : ', s_i_new(ji)
-                  WRITE(numout,*) ' sss_m     : ', sss_m(zji,zjj)
+                  WRITE(numout,*) ' sss_m     : ', sss_m(ii,ij)
                   WRITE(numout,*) ' dh_i_bott : ', dh_i_bott(ji)
                   WRITE(numout,*) ' innermelt : ', INT( zinnermelt(ji) )
@@ -531,5 +525,5 @@
          !                     ! excessive energy is sent to lateral ablation
          fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 )   &
-            &                          * ( zdhbf - dh_i_bott(ji) ) / rdt_ice
+            &                          * ( zdhbf - dh_i_bott(ji) ) * r1_rdtice
          dh_i_bott(ji)  = zdhbf
          !                     !since ice volume is only used for outputs, we keep it global for all categories
@@ -538,9 +532,9 @@
          zhgnew   (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji)
          !                     ! diagnostic ( bottom ice growth )
-         zji = MOD( npb(ji) - 1, jpi ) + 1
-         zjj = ( npb(ji) - 1 ) / jpi + 1
-         diag_bot_gr(zji,zjj) = diag_bot_gr(zji,zjj) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice
-         diag_sur_me(zji,zjj) = diag_sur_me(zji,zjj) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) / rdt_ice
-         diag_bot_me(zji,zjj) = diag_bot_me(zji,zjj) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) / rdt_ice
+         ii = MOD( npb(ji) - 1, jpi ) + 1
+         ij = ( npb(ji) - 1 ) / jpi + 1
+         diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice
+         diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice
+         diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice
       END DO
 
@@ -548,11 +542,10 @@
       ! 5.2 More than available ice melts
       !-----------------------------------
-      ! then heat applied minus heat content at previous time step
-      ! should equal heat remaining 
+      ! then heat applied minus heat content at previous time step should equal heat remaining 
       !
       DO ji = kideb, kiut
          ! Adapt the remaining energy if too much ice melts
          !--------------------------------------------------
-         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice
+         zihgnew =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice
          ! 0 if no more ice
          zhgnew    (ji) =         zihgnew   * zhgnew(ji)      ! ice thickness is put to 0
@@ -562,10 +555,10 @@
          ! If snow remains, energy is used to melt snow
          zhni =  ht_s_b(ji)      ! snow depth at previous time step
-         zihg =  MAX( zzero , SIGN ( zone , - ht_s_b(ji) ) ) ! 0 if snow 
+         zihg =  MAX(  zzero , SIGN ( zone , - ht_s_b(ji) )  )   ! =0 if snow 
 
          ! energy of melting of remaining snow
          zqt_s(ji) =    ( 1. - zihg ) * zqt_s(ji) / MAX( zhni, epsi13 )
          zdhnm     =  - ( 1. - zihg ) * ( 1. - zihgnew ) * zfdt_final(ji) / MAX( zqt_s(ji) , epsi13 )
-         zhnfi          =  zhni + zdhnm
+         zhnfi     =  zhni + zdhnm
          zfdt_final(ji) =  MAX( zfdt_final(ji) + zqt_s(ji) * zdhnm , 0.0 )
          ht_s_b(ji)     =  MAX( zzero , zhnfi )
@@ -581,14 +574,14 @@
          !
          !                                              ! mass variation cumulated over category
-         rdmsnif_1d(ji) = rdmsnif_1d(ji) + zzfmass_s                     ! snow 
-         rdmicif_1d(ji) = rdmicif_1d(ji) + zzfmass_i                     ! ice 
+         rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s                     ! snow 
+         rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i                     ! ice 
 
          ! Remaining heat to the ocean 
          !---------------------------------
-         focea(ji)  = - zfdt_final(ji) / rdt_ice         ! focea is in W.m-2 * dt
-
-      END DO
-
-      ftotal_fin (:) = zfdt_final(:)  / rdt_ice
+         focea(ji)  = - zfdt_final(ji) * r1_rdtice         ! focea is in W.m-2 * dt
+
+      END DO
+
+      ftotal_fin (:) = zfdt_final(:)  * r1_rdtice
 
       !---------------------------
@@ -596,26 +589,18 @@
       !---------------------------
       DO ji = kideb, kiut
-         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   !1 if ice
-
+         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice
+         !
          ! Salt flux
-         zji = MOD( npb(ji) - 1, jpi ) + 1
-         zjj = ( npb(ji) - 1 ) / jpi + 1
-         ! new lines
-         IF( num_sal == 4 ) THEN
-            fseqv_1d(ji) = fseqv_1d(ji) +        zihgnew  * zfsalt_melt(ji)                                &
-               &                        + (1.0 - zihgnew) * zfmass_i(ji) * ( sss_m(zji,zjj) - bulk_sal   ) / rdt_ice
-         ELSE
-            fseqv_1d(ji) = fseqv_1d(ji) +        zihgnew  * zfsalt_melt(ji)                                &
-               &                        + (1.0 - zihgnew) * zfmass_i(ji) * ( sss_m(zji,zjj) - sm_i_b(ji) ) / rdt_ice
-         ENDIF
+         fseqv_1d(ji) = fseqv_1d(ji) +        zihgnew  * zfsalt_melt(ji)                                &
+            &                        - (1.0 - zihgnew) * zfmass_i   (ji) * sm_i_b(ji) ) * r1_rdtice
+         !
          ! Heat flux
          ! excessive bottom ablation energy (fsup) - 0 except if jpl = 1
-         ! excessive total ablation energy (focea) sent to the ocean
+         ! excessive total  ablation energy (focea) sent to the ocean
          qfvbq_1d(ji)  = qfvbq_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji) * a_i_b(ji) * rdt_ice
 
-         zihic   = 1.0 - MAX( zzero , SIGN( zone , -ht_i_b(ji) ) )
-         ! equals 0 if ht_i = 0, 1 if ht_i gt 0
+         zihic   = 1.0 - MAX(  zzero , SIGN( zone , -ht_i_b(ji) )  )      ! equals 0 if ht_i = 0, 1 if ht_i gt 0
          fscbq_1d(ji) =  a_i_b(ji) * fstbif_1d(ji)
-         qldif_1d(ji)  = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea(ji)    * a_i_b(ji) * rdt_ice   &
+         qldif_1d(ji)  = qldif_1d(ji) + fsup(ji) + ( 1.0 - zihgnew ) * focea   (ji) * a_i_b(ji) * rdt_ice   &
             &                                    + ( 1.0 - zihic   ) * fscbq_1d(ji)             * rdt_ice
       END DO  ! ji
@@ -656,29 +641,24 @@
          dmgwi_1d  (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn
 
-         rdmicif_1d(ji) = rdmicif_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic 
-         rdmsnif_1d(ji) = rdmsnif_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn
+         rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic 
+         rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn
 
          !        Equivalent salt flux (1) Snow-ice formation component
          !        -----------------------------------------------------
-         zji = MOD( npb(ji) - 1, jpi ) + 1
-         zjj =    ( npb(ji) - 1 ) / jpi + 1
-
-         IF( num_sal /= 2 ) THEN   ;   zsm_snowice = sm_i_b(ji)
-         ELSE                      ;   zsm_snowice = ( rhoic - rhosn ) / rhoic * sss_m(zji,zjj) 
+         ii = MOD( npb(ji) - 1, jpi ) + 1
+         ij =    ( npb(ji) - 1 ) / jpi + 1
+
+         IF( num_sal == 2 ) THEN   ;   zsm_snowice = sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic
+         ELSE                      ;   zsm_snowice = sm_i_b(ji)   
          ENDIF
-         IF( num_sal == 4 ) THEN
-            fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - bulk_sal    ) * a_i_b(ji)   &
-               &                        * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
-         ELSE
-            fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - zsm_snowice ) * a_i_b(ji)   &
-               &                        * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic / rdt_ice
-         ENDIF
+         fseqv_1d(ji) = fseqv_1d(ji) - zsm_snowice  * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice
+         !
          ! entrapment during snow ice formation
          i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + 1.0e-6 ) )
          isnowic      = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji)      ) ) * i_ice_switch
-         IF(  num_sal == 2  .OR.  num_sal == 4  )   &
-            dsm_i_si_1d(ji) = ( zsm_snowice*dh_snowice(ji) &
-            &               + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13)   &
-            &               - sm_i_b(ji) ) * isnowic     
+         IF(  num_sal == 2  )   &
+            dsm_i_si_1d(ji) = (  zsm_snowice * dh_snowice(ji)    &
+            &                  + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13 )   &
+            &                  - sm_i_b(ji)  ) * isnowic     
 
          !  Actualize new snow and ice thickness.
@@ -690,7 +670,7 @@
 
          ! diagnostic ( snow ice growth )
-         zji = MOD( npb(ji) - 1, jpi ) + 1
-         zjj =    ( npb(ji) - 1 ) / jpi + 1
-         diag_sni_gr(zji,zjj)  = diag_sni_gr(zji,zjj) + dh_snowice(ji)*a_i_b(ji) / rdt_ice
+         ii = MOD( npb(ji) - 1, jpi ) + 1
+         ij =    ( npb(ji) - 1 ) / jpi + 1
+         diag_sni_gr(ii,ij)  = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice
          !
       END DO !ji
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 3517)
@@ -147,6 +147,5 @@
       REAL(wp), DIMENSION(kiut,jkmax+2) ::   zdiagbis
       REAL(wp), DIMENSION(kiut,jkmax+2,3) ::   ztrid   ! tridiagonal system terms
-      !!------------------------------------------------------------------
-      
+      !!------------------------------------------------------------------     
       ! 
       !------------------------------------------------------------------------------!
@@ -156,5 +155,5 @@
       DO ji = kideb , kiut
          ! is there snow or not
-         isnow(ji)= INT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) )  )
+         isnow(ji)= INT(  1._wp - MAX(  0._wp , SIGN(1._wp, - ht_s_b(ji) )  )  )
          ! surface temperature of fusion
 !!gm ???  ztfs(ji) = rtt !!!????
@@ -201,5 +200,5 @@
       DO ji = kideb , kiut
          ! switches
-         isnow(ji) = INT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )  ) 
+         isnow(ji) = INT(  1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_b(ji) )  )  ) 
          ! hs > 0, isnow = 1
          zhsu (ji) = hnzst  ! threshold for the computation of i0
@@ -262,6 +261,6 @@
       ! just to check energy conservation
       DO ji = kideb, kiut
-         ii                = MOD( npb(ji) - 1, jpi ) + 1
-         ij                = ( npb(ji) - 1 ) / jpi + 1
+         ii = MOD( npb(ji) - 1 , jpi ) + 1
+         ij =    ( npb(ji) - 1 ) / jpi + 1
          fstroc(ii,ij,jl) = zradtr_i(ji,nlay_i)
       END DO
@@ -273,5 +272,4 @@
          END DO
       END DO
-
 
       !
@@ -662,7 +660,7 @@
 
             ! surface temperature
-            isnow(ji)     = INT(1.0-max(0.0,sign(1.0,-ht_s_b(ji))))
+            isnow(ji)     = INT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) )  )  )
             ztsuoldit(ji) = t_su_b(ji)
-            IF (t_su_b(ji) .LT. ztfs(ji)) &
+            IF( t_su_b(ji) < ztfs(ji) )   &
                t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( isnow(ji)*t_s_b(ji,1)   &
                &          + (1.0-isnow(ji))*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))  
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90	(revision 3517)
@@ -408,13 +408,13 @@
       IF ( con_i ) THEN
          DO ji = kideb, kiut
-            IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN
+            IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN
                zji                 = MOD( npb(ji) - 1, jpi ) + 1
                zjj                 = ( npb(ji) - 1 ) / jpi + 1
                WRITE(numout,*) ' violation of heat conservation : ',             &
-                  ABS ( zqts_in(ji) - zqts_fin(ji) ) / rdt_ice
+                  ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice
                WRITE(numout,*) ' ji, jj   : ', zji, zjj
                WRITE(numout,*) ' ht_s_b   : ', ht_s_b(ji)
-               WRITE(numout,*) ' zqts_in  : ', zqts_in(ji) / rdt_ice
-               WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) / rdt_ice
+               WRITE(numout,*) ' zqts_in  : ', zqts_in (ji) * r1_rdtice
+               WRITE(numout,*) ' zqts_fin : ', zqts_fin(ji) * r1_rdtice
                WRITE(numout,*) ' dh_snowice : ', dh_snowice(ji)
                WRITE(numout,*) ' dh_s_tot : ', dh_s_tot(ji)
@@ -526,5 +526,5 @@
          ! bottom formation temperature
          ztform = t_i_b(ji,nlay_i)
-         IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) ztform = t_bo_b(ji)
+         IF(  num_sal == 2  )   ztform = t_bo_b(ji)
          qm0(ji,nbot0(ji)) = ( 1.0 - icboswi(ji) )*qm0(ji,nbot0(ji))             &   ! case of melting ice
             &              + icboswi(ji) * rhoic * ( cpic*(ztmelts-ztform)       &   ! case of forming ice
@@ -622,12 +622,12 @@
       !
       DO ji = kideb, kiut
-         IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice .GT. 1.0e-6 ) THEN
+         IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN
             zji                 = MOD( npb(ji) - 1, jpi ) + 1
             zjj                 = ( npb(ji) - 1 ) / jpi + 1
-            WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) / rdt_ice
+            WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice
             WRITE(numout,*) ' ji, jj   : ', zji, zjj
             WRITE(numout,*) ' ht_i_b   : ', ht_i_b(ji)
-            WRITE(numout,*) ' zqti_in  : ', zqti_in(ji) / rdt_ice
-            WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) / rdt_ice
+            WRITE(numout,*) ' zqti_in  : ', zqti_in (ji) * r1_rdtice
+            WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice
             WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji)
             WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji)
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90	(revision 3517)
@@ -143,9 +143,9 @@
       ! 1) Conservation check and changes in each ice category
       !------------------------------------------------------------------------------!
-      IF ( con_i ) THEN
-         CALL lim_column_sum (jpl, v_i, vt_i_init)
-         CALL lim_column_sum (jpl, v_s, vt_s_init)
-         CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init)
-         CALL lim_column_sum (jpl,   e_s(:,:,1,:) , et_s_init)
+      IF( con_i ) THEN
+         CALL lim_column_sum        ( jpl, v_i          , vt_i_init)
+         CALL lim_column_sum        ( jpl, v_s          , vt_s_init)
+         CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init)
+         CALL lim_column_sum        ( jpl, e_s(:,:,1,:) , et_s_init)
       ENDIF
 
@@ -158,9 +158,7 @@
                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
-                  zindb      = 1.0-MAX(0.0,SIGN(1.0,-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
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * 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
                END DO
             END DO
@@ -182,19 +180,15 @@
       ! 
 
-      zvrel(:,:) = 0.0
+      zvrel(:,:) = 0._wp
 
       ! Default new ice thickness 
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            hicol(ji,jj) = hiccrit(1)
-         END DO
-      END DO
-
-      IF (fraz_swi.eq.1.0) THEN
+      hicol(:,:) = hiccrit(1)
+
+      IF( fraz_swi == 1._wp ) THEN
 
          !--------------------
          ! Physical constants
          !--------------------
-         hicol(:,:) = 0.0
+         hicol(:,:) = 0._wp
 
          zhicrit = 0.04 ! frazil ice thickness
@@ -306,33 +300,21 @@
       IF ( nbpac > 0 ) THEN
 
-         CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         ,       &
-            jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) )
          DO jl = 1, jpl
-            CALL tab_2d_1d( nbpac, za_i_ac(1:nbpac,jl)  , a_i(:,:,jl)  ,       &
-               jpi, jpj, npac(1:nbpac) )
-            CALL tab_2d_1d( nbpac, zv_i_ac(1:nbpac,jl)  , v_i(:,:,jl)  ,       &
-               jpi, jpj, npac(1:nbpac) )
-            CALL tab_2d_1d( nbpac, zoa_i_ac(1:nbpac,jl) , oa_i(:,:,jl) ,       &
-               jpi, jpj, npac(1:nbpac) )
-            CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl),       &
-               jpi, jpj, npac(1:nbpac) )
+            CALL tab_2d_1d( nbpac, za_i_ac(1:nbpac,jl)  , a_i(:,:,jl)  , jpi, jpj, npac(1:nbpac) )
+            CALL tab_2d_1d( nbpac, zv_i_ac(1:nbpac,jl)  , v_i(:,:,jl)  , jpi, jpj, npac(1:nbpac) )
+            CALL tab_2d_1d( nbpac, zoa_i_ac(1:nbpac,jl) , oa_i(:,:,jl) , jpi, jpj, npac(1:nbpac) )
+            CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) )
             DO jk = 1, nlay_i
-               CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , &
-                  jpi, jpj, npac(1:nbpac) )
+               CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) )
             END DO ! jk
          END DO ! jl
 
-         CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif ,              &
-            jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif ,              &
-            jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo  ,              &
-            jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, fseqv_1d  (1:nbpac)     , fseqv ,              &
-            jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol ,              &
-            jpi, jpj, npac(1:nbpac) )
-         CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel ,              &
-            jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, qldif_1d  (1:nbpac)     , qldif , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, qcmif_1d  (1:nbpac)     , qcmif , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, t_bo_b    (1:nbpac)     , t_bo  , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, fseqv_1d  (1:nbpac)     , fseqv , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol , jpi, jpj, npac(1:nbpac) )
+         CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel , jpi, jpj, npac(1:nbpac) )
 
          !------------------------------------------------------------------------------!
@@ -344,7 +326,7 @@
          !----------------------
          DO ji = 1, nbpac
-            zh_newice(ji)     = hiccrit(1)
-         END DO
-         IF ( fraz_swi .EQ. 1.0 ) zh_newice(:) = hicol_b(:)
+            zh_newice(ji) = hiccrit(1)
+         END DO
+         IF( fraz_swi == 1.0 )   zh_newice(:) = hicol_b(:)
 
          !----------------------
@@ -352,22 +334,17 @@
          !----------------------
 
-         IF ( num_sal .EQ. 1 ) THEN
-            zs_newice(:)      =   bulk_sal
-         ENDIF ! num_sal
-
-         IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN
-
-            DO ji = 1, nbpac
-               zs_newice(ji)  =   MIN( 4.606 + 0.91 / zh_newice(ji) , s_i_max )
-               zji            =   MOD( npac(ji) - 1, jpi ) + 1
-               zjj            =   ( npac(ji) - 1 ) / jpi + 1
-               zs_newice(ji)  =   MIN( 0.5*sss_m(zji,zjj) , zs_newice(ji) )
-            END DO ! jl
-
-         ENDIF ! num_sal
-
-         IF ( num_sal .EQ. 3 ) THEN
-            zs_newice(:)      =   2.3
-         ENDIF ! num_sal
+         SELECT CASE ( num_sal )
+         CASE ( 1 )                    ! Sice = constant 
+            zs_newice(:) = bulk_sal
+         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)  )
+            END DO
+         CASE ( 3 )                    ! Sice = F(z) [multiyear ice]
+            zs_newice(:) =   2.3
+         END SELECT
+
 
          !-------------------------
@@ -376,12 +353,10 @@
          ! We assume that new ice is formed at the seawater freezing point
          DO ji = 1, nbpac
-            ztmelts           = - tmut * zs_newice(ji) + rtt ! Melting point (K)
-            ze_newice(ji)     =   rhoic * ( cpic * ( ztmelts - t_bo_b(ji) )    &
-               + lfus * ( 1.0 - ( ztmelts - rtt )   &
-               / ( t_bo_b(ji) - rtt ) )           &
-               - rcp * ( ztmelts-rtt ) )
-            ze_newice(ji)     =   MAX( ze_newice(ji) , 0.0 ) +                 &
-               MAX( 0.0 , SIGN( 1.0 , - ze_newice(ji) ) )   & 
-               * rhoic * lfus
+            ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K)
+            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             &
+               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   &
+               &                       - rcp  *         ( ztmelts - rtt )  )
+            ze_newice(ji) =   MAX( ze_newice(ji) , 0._wp )    &
+               &          +   MAX(  0.0 , SIGN( 1.0 , - ze_newice(ji) )  ) * rhoic * lfus
          END DO ! ji
          !----------------
@@ -389,5 +364,5 @@
          !----------------
          DO ji = 1, nbpac
-            zo_newice(ji)     = 0.0
+            zo_newice(ji) = 0._wp
          END DO ! ji
 
@@ -396,5 +371,5 @@
          !--------------------------
          DO ji = 1, nbpac
-            zqbgow(ji)        = qldif_1d(ji) - qcmif_1d(ji) !<0
+            zqbgow(ji) = qldif_1d(ji) - qcmif_1d(ji)     !<0
          END DO ! ji
 
@@ -403,10 +378,9 @@
          !-------------------
          DO ji = 1, nbpac
-            zv_newice(ji)     = - zqbgow(ji) / ze_newice(ji)
+            zv_newice(ji) = - zqbgow(ji) / ze_newice(ji)
 
             ! A fraction zfrazb of frazil ice is accreted at the ice bottom
-            zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) )     & 
-               + 1.0 ) / 2.0 * maxfrazb
-            zdh_frazb(ji) = zfrazb*zv_newice(ji)
+            zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb
+            zdh_frazb(ji) =         zfrazb   * zv_newice(ji)
             zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji)
          END DO
@@ -415,21 +389,8 @@
          ! Salt flux due to new ice growth
          !---------------------------------
-         IF ( ( num_sal .EQ. 4 ) ) THEN 
-            DO ji = 1, nbpac
-               zji            = MOD( npac(ji) - 1, jpi ) + 1
-               zjj            = ( npac(ji) - 1 ) / jpi + 1
-               fseqv_1d(ji)   = fseqv_1d(ji) +                                     &
-                  ( sss_m(zji,zjj) - bulk_sal      ) * rhoic *       &
-                  zv_newice(ji) / rdt_ice
-            END DO
-         ELSE
-            DO ji = 1, nbpac
-               zji            = MOD( npac(ji) - 1, jpi ) + 1
-               zjj            = ( npac(ji) - 1 ) / jpi + 1
-               fseqv_1d(ji)   = fseqv_1d(ji) +                                     &
-                  ( sss_m(zji,zjj) - zs_newice(ji) ) * rhoic *       &
-                  zv_newice(ji) / rdt_ice
-            END DO ! ji
-         ENDIF
+         ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine)
+         DO ji = 1, nbpac
+            fseqv_1d(ji) = fseqv_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice
+         END DO ! ji
 
          !------------------------------------
@@ -448,6 +409,5 @@
 
          ! keep new ice volume in memory
-         CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , &
-            jpi, jpj )
+         CALL tab_1d_2d( nbpac, v_newice , npac(1:nbpac), zv_newice(1:nbpac) , jpi, jpj )
 
          !-----------------
@@ -459,5 +419,5 @@
             zji                  = MOD( npac(ji) - 1, jpi ) + 1
             zjj                  = ( npac(ji) - 1 ) / jpi + 1
-            diag_lat_gr(zji,zjj) = zv_newice(ji) / rdt_ice
+            diag_lat_gr(zji,zjj) = zv_newice(ji) * r1_rdtice
          END DO !ji
 
@@ -628,8 +588,8 @@
          ! Update salinity
          !-----------------
-         IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN
+         IF(  num_sal == 2  ) THEN      ! Sice = F(z,t)
             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) ) )   ! 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
@@ -645,5 +605,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  .OR.  num_sal == 4  )   &
+            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
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90	(revision 3517)
@@ -44,14 +44,12 @@
       !! ** Purpose :   computes new salinities in the ice
       !!
-      !! ** Method  :  4 possibilities
-      !!               -> num_sal = 1 -> constant salinity for z,t
-      !!               -> num_sal = 2 -> S = S(z,t) [simple Vancoppenolle et al 2005]
-      !!               -> num_sal = 3 -> S = S(z)   [multiyear ice]
-      !!               -> num_sal = 4 -> S = S(h)   [Cox and Weeks 74]
+      !! ** Method  :  3 possibilities
+      !!               -> num_sal = 1 -> Sice = cst    [ice salinity constant in both time & space] 
+      !!               -> num_sal = 2 -> Sice = S(z,t) [Vancoppenolle et al. 2005]
+      !!               -> num_sal = 3 -> Sice = S(z)   [multiyear ice]
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::  kideb, kiut   ! thickness category index
+      INTEGER, INTENT(in) ::   kideb, kiut   ! thickness category index
       !
       INTEGER  ::   ji, jk     ! dummy loop indices 
-      INTEGER  ::   zji, zjj   ! local integers
       REAL(wp) ::   zsold, iflush, iaccrbo, igravdr, isnowic, i_ice_switch,  ztmelts   ! local scalars
       REAL(wp) ::   zaaa, zbbb, zccc, zdiscrim   ! local scalars
@@ -64,17 +62,10 @@
       ! 1) Constant salinity, constant in time                                       |
       !------------------------------------------------------------------------------|
-!!gm comment: if num_sal = 1 s_i_b and sm_i_b can be set to bulk_sal one for all in the initialisation phase !!
-      IF( num_sal == 1 ) THEN
-         !
-         DO jk = 1, nlay_i
-            DO ji = kideb, kiut
-               s_i_b(ji,jk) =  bulk_sal
-            END DO ! ji
-         END DO ! jk
-         !
-         DO ji = kideb, kiut
-            sm_i_b(ji)      =  bulk_sal 
-         END DO ! ji
-         !
+!!gm comment: if num_sal = 1 s_i_new, s_i_b and sm_i_b can be set to bulk_sal one for all in the initialisation phase !!
+!!gm           ===>>>   simplification of almost all test on num_sal value
+      IF(  num_sal == 1  ) THEN
+            s_i_b  (kideb:kiut,1:nlay_i) =  bulk_sal
+            sm_i_b (kideb:kiut)          =  bulk_sal 
+            s_i_new(kideb:kiut)          =  bulk_sal
       ENDIF
 
@@ -83,5 +74,5 @@
       !------------------------------------------------------------------------------|
 
-      IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN
+      IF(  num_sal == 2  ) THEN
 
          !---------------------------------
@@ -118,5 +109,5 @@
             dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice 
             !                                   ! drainage by flushing  
-            dsm_i_fl_1d(ji) = - iflush * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice
+            dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice
 
             !-----------------
@@ -133,6 +124,6 @@
          END DO ! ji
 
-         ! Salinity profile
-         CALL lim_var_salprof1d( kideb, kiut )
+         CALL lim_var_salprof1d( kideb, kiut )         ! Salinity profile
+
 
          !----------------------------
@@ -143,9 +134,9 @@
 !!gm useless
             ! iflush  : 1 if summer 
-            iflush  =  MAX( 0._wp , SIGN ( 1._wp , t_su_b(ji) - rtt ) ) 
+            iflush  =  MAX( 0._wp , SIGN( 1._wp , t_su_b(ji) - rtt        )  ) 
             ! igravdr : 1 if t_su lt t_bo
-            igravdr =  MAX( 0._wp , SIGN ( 1._wp , t_bo_b(ji) - t_su_b(ji) ) ) 
+            igravdr =  MAX( 0._wp , SIGN( 1._wp , t_bo_b(ji) - t_su_b(ji) )  ) 
             ! iaccrbo : 1 if bottom accretion
-            iaccrbo =  MAX( 0._wp , SIGN ( 1._wp , dh_i_bott(ji) ) )
+            iaccrbo =  MAX( 0._wp , SIGN( 1._wp , dh_i_bott(ji)           )  )
 !!gm end useless
             !
@@ -157,8 +148,7 @@
          !----------------------------
          DO ji = kideb, kiut
-            i_ice_switch = 1._wp - MAX ( 0._wp, SIGN( 1._wp , - ht_i_b(ji) ) )
+            i_ice_switch = 1._wp - MAX(  0._wp, SIGN( 1._wp , - ht_i_b(ji) )  )
             fsbri_1d(ji) = fsbri_1d(ji) - i_ice_switch * rhoic * a_i_b(ji) * ht_i_b(ji)         &
-               &         * ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), sm_i_b(ji) - zsiold(ji) ) ) / rdt_ice
-            IF( num_sal == 4 ) fsbri_1d(ji) = 0._wp
+               &         * ( MAX(dsm_i_gd_1d(ji) + dsm_i_fl_1d(ji), sm_i_b(ji) - zsiold(ji) ) ) * r1_rdtice
          END DO ! ji
 
@@ -179,5 +169,5 @@
          END DO
          !
-      ENDIF ! num_sal .EQ. 2
+      ENDIF 
 
       !------------------------------------------------------------------------------|
@@ -185,49 +175,14 @@
       !------------------------------------------------------------------------------|
 
-      IF( num_sal == 3 )   CALL lim_var_salprof1d( kideb, kiut )
-
-      !------------------------------------------------------------------------------|
-      !  Module 4 : Constant salinity varying in time                                |
-      !------------------------------------------------------------------------------|
-
-      IF( num_sal == 5 ) THEN      ! Cox and Weeks, 1974
-         !
-         DO ji = kideb, kiut
-            zsold = sm_i_b(ji)
-            IF( ht_i_b(ji) < 0.4 ) THEN
-               sm_i_b(ji) = 14.24 - 19.39 * ht_i_b(ji) 
-            ELSE
-               sm_i_b(ji) =  7.88 - 1.59 * ht_i_b(ji)
-               sm_i_b(ji) = MIN( sm_i_b(ji) , zsold )  
-            ENDIF
-            IF( ht_i_b(ji) > 3.06918239 ) THEN 
-               sm_i_b(ji) = 3._wp
-            ENDIF
-            DO jk = 1, nlay_i
-               s_i_b(ji,jk)   = sm_i_b(ji)
-            END DO
-         END DO
-         !
-      ENDIF ! num_sal
+      IF(  num_sal == 3  )   CALL lim_var_salprof1d( kideb, kiut )
+
 
       !------------------------------------------------------------------------------|
       ! 5) Computation of salt flux due to Bottom growth
       !------------------------------------------------------------------------------|
-
-      IF ( num_sal == 4 ) THEN
-         DO ji = kideb, kiut
-            zji = MOD( npb(ji) - 1 , jpi ) + 1
-            zjj =    ( npb(ji) - 1 ) / jpi + 1
-            fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - bulk_sal    )               &
-               &                        * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice
-         END DO
-      ELSE
-         DO ji = kideb, kiut
-            zji = MOD( npb(ji) - 1 , jpi ) + 1
-            zjj =    ( npb(ji) - 1 ) / jpi + 1
-            fseqv_1d(ji) = fseqv_1d(ji) + ( sss_m(zji,zjj) - s_i_new(ji) )               &
-               &                        * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0.0 ) / rdt_ice
-         END DO
-      ENDIF
+      ! note: s_i_new = bulk_sal in constant salinity case
+      DO ji = kideb, kiut
+         fseqv_1d(ji) = fseqv_1d(ji) - s_i_new(ji) * rhoic * a_i_b(ji) * MAX( dh_i_bott(ji) , 0._wp ) * r1_rdtice
+      END DO
       !
       CALL wrk_dealloc( jpij, ze_init, zhiold, zsiold )
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 3517)
@@ -392,12 +392,10 @@
 
                   ! Ice salinity and age
-                  zsal = MAX( MIN( (rhoic-rhosn)/rhoic*sss_m(ji,jj)  , &
-                     zusvoic * zs0sm(ji,jj,jl) ), s_i_min ) * v_i(ji,jj,jl)
-                  IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) & 
-                     smv_i(ji,jj,jl) = zindic*zsal + (1.0-zindic)*0.0
-
-                  zage = MAX( MIN( zbigval, zs0oi(ji,jj,jl) / & 
-                     MAX( a_i(ji,jj,jl), epsi16 )  ), 0.0 ) * a_i(ji,jj,jl)
-                  oa_i (ji,jj,jl)  = zindic*zage 
+                  zsal = MAX( MIN(  (rhoic-rhosn)/rhoic*sss_m(ji,jj) ,   &
+                     &              zusvoic * zs0sm(ji,jj,jl)         ) , s_i_min ) * v_i(ji,jj,jl)
+                  IF(  num_sal == 2  )   smv_i(ji,jj,jl) = zindic * zsal + (1.0-zindic) * 0._wp
+
+                  zage = MAX(  MIN( zbigval, zs0oi(ji,jj,jl) / MAX( a_i(ji,jj,jl), epsi16 ) ), 0._wp  ) * a_i(ji,jj,jl)
+                  oa_i (ji,jj,jl)  = zindic * zage 
 
                   ! Snow heat content
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90	(revision 3517)
@@ -190,24 +190,22 @@
 
                   ! is there any ice left ?
-                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 
+                  zindic = MAX( rzero , SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 
                   !=1 if hi > 1e-3 and 0 if not
-                  zdvres        = MAX(0.0,-v_i(ji,jj,jl)) !residual volume if too much ice was molten
+                  zdvres = MAX( 0.0 , -v_i(ji,jj,jl) ) !residual volume if too much ice was molten
                   !this quantity is positive
-                  v_i(ji,jj,jl) = zindic*v_i(ji,jj,jl)    !ice volume cannot be negative
+                  v_i(ji,jj,jl) = zindic * v_i(ji,jj,jl)    !ice volume cannot be negative
                   !correct thermodynamic ablation
-                  d_v_i_thd(ji,jj,jl)  = zindic  *  d_v_i_thd(ji,jj,jl) + (1.0-zindic) * (-zviold - d_v_i_trp(ji,jj,jl)) 
+                  d_v_i_thd(ji,jj,jl) = zindic  *  d_v_i_thd(ji,jj,jl) + (1.0-zindic) * (-zviold - d_v_i_trp(ji,jj,jl)) 
                   ! THIS IS NEW
-                  d_a_i_thd(ji,jj,jl)  = zindic  *  d_a_i_thd(ji,jj,jl) + & 
-                     (1.0-zindic) * (-old_a_i(ji,jj,jl)) 
+                  d_a_i_thd(ji,jj,jl) = zindic  *  d_a_i_thd(ji,jj,jl) + (1.0-zindic) * (-old_a_i(ji,jj,jl)) 
 
                   !residual salt flux if ice is over-molten
-                  fsalt_res(ji,jj)  = fsalt_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) ) * & 
-                     ( rhoic * zdvres / rdt_ice )
-                  !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhoic * lfus * zdvres / rdt_ice
+                  fsalt_res(ji,jj)  = fsalt_res(ji,jj) - sm_i(ji,jj,jl) * ( rhoic * zdvres * r1_rdtice )
+                  !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhoic * lfus * zdvres * r1_rdtice
 
                   ! is there any snow left ?
-                  zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) ) 
-                  zvsold        = v_s(ji,jj,jl)
-                  zdvres        = MAX(0.0,-v_s(ji,jj,jl)) !residual volume if too much ice was molten
+                  zindsn = MAX(  rzero, SIGN( rone , v_s(ji,jj,jl) - epsi10 )  ) 
+                  zvsold = v_s(ji,jj,jl)
+                  zdvres = MAX( 0.0 , -v_s(ji,jj,jl) )   !residual volume if too much ice was molten
                   !this quantity is positive
                   v_s(ji,jj,jl) = zindsn*v_s(ji,jj,jl)    !snow volume cannot be negative
@@ -215,10 +213,7 @@
                   d_v_s_thd(ji,jj,jl)  = zindsn  *  d_v_s_thd(ji,jj,jl) + & 
                      (1.0-zindsn) * (-zvsold - d_v_s_trp(ji,jj,jl)) 
-                  !unsure correction on salt flux.... maybe future will tell it was not that right
-
-                  !residual salt flux if snow is over-molten
-                  fsalt_res(ji,jj)  = fsalt_res(ji,jj) + sss_m(ji,jj) * ( rhosn * zdvres / rdt_ice )
-                  !this flux will be positive if snow was over-molten
-                  !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhosn * lfus * zdvres / rdt_ice
+
+                  !no salt flux when snow is over-molten
+                  !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhosn * lfus * zdvres * r1_rdtice
                ENDIF
             END DO !ji
@@ -306,5 +301,5 @@
       !--------------
 
-      IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN      ! general case
+      IF(  num_sal == 2  ) THEN      ! Prognostic salinity [Sice=F(z,t)]
          !
          IF( ln_nicep ) THEN  
@@ -317,5 +312,5 @@
             WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
          ENDIF
-
+         !
          smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_thd(:,:,:) + d_smv_i_trp(:,:,:)
          !
@@ -616,11 +611,10 @@
             DO ji = 1, jpi
                IF ( internal_melt(ji,jj,jl) == 1 ) THEN
-                  v_s(ji,jj,jl)   = 0.0
-                  e_s(ji,jj,1,jl) = 0.0
+                  v_s(ji,jj,jl)   = 0._wp
+                  e_s(ji,jj,1,jl) = 0._wp
                   !   ! release heat
-                  fheat_res(ji,jj) = fheat_res(ji,jj)  &
-                     + ze_s * v_s(ji,jj,jl) / rdt_ice
+                  fheat_res(ji,jj) = fheat_res(ji,jj) + ze_s  * v_s(ji,jj,jl) * r1_rdtice
                   ! release mass
-                  rdm_snw(ji,jj) =  rdm_snw(ji,jj) + rhosn * v_s(ji,jj,jl)
+                  rdm_snw  (ji,jj) = rdm_snw  (ji,jj) + rhosn * v_s(ji,jj,jl)
                ENDIF
             END DO
@@ -648,5 +642,5 @@
                !                ENDIF
                IF ((oa_i(ji,jj,jl)-1.0)*86400.0.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN
-                  oa_i(ji,jj,jl) = rdt_ice*numit/86400.0*a_i(ji,jj,jl)
+                  oa_i(ji,jj,jl) = rdt_ice * numit / 86400.0 * a_i(ji,jj,jl)
                ENDIF
                oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl)
@@ -657,5 +651,5 @@
                !             v_s(ji,jj,jl)  = MAX( zindb * v_s(ji,jj,jl), 0.0) 
                ! snow thickness cannot be smaller than 1e-6
-               v_s(ji,jj,jl)  = zindsn*v_s(ji,jj,jl)*zindb
+               v_s(ji,jj,jl)  = zindsn * v_s(ji,jj,jl) * zindb
                v_s(ji,jj,jl)  = v_s(ji,jj,jl) *  MAX( 0.0 , SIGN( 1.0 , v_s(ji,jj,jl) - epsi06 ) )
 
@@ -737,21 +731,18 @@
       !---------------------
 
-      IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) THEN ! general case
-
+      IF(  num_sal == 2  ) THEN      ! Prognostic salinity [Sice=F(z,t)]
+         !
          DO jl = 1, jpl
             DO jk = 1, nlay_i
                DO jj = 1, jpj 
-                  DO ji = 1, jpi
-                     ! salinity stays in bounds
-                     smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)), &
-                        0.1 * v_i(ji,jj,jl) )
+                  DO ji = 1, jpi           ! salinity stays in bounds
+                     smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)), 0.1 * v_i(ji,jj,jl) )
                      i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl)))
-                     smv_i(ji,jj,jl)  = i_ice_switch*smv_i(ji,jj,jl) + &
-                        0.1*(1.0-i_ice_switch)*v_i(ji,jj,jl)
+                     smv_i(ji,jj,jl)  = i_ice_switch*smv_i(ji,jj,jl) + 0.1*(1.0-i_ice_switch)*v_i(ji,jj,jl)
                   END DO ! ji
                END DO ! jj
             END DO !jk
          END DO !jl
-
+         !
       ENDIF
 
@@ -819,4 +810,5 @@
             ! 1) Count the number of existing categories
             DO jl  = 1, jpl
+!!cr : comment the second line of zindb definition, and use epsi04 in the 1st one
                zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi03 ) ) 
                zindb   =  MAX( rzero, SIGN( rone, v_i(ji,jj,jl) ) ) 
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 3517)
@@ -182,5 +182,5 @@
       END DO
 
-      IF(  num_sal == 2  .OR.  num_sal == 4  )THEN
+      IF(  num_sal == 2  )THEN
          DO jl = 1, jpl
             DO jj = 1, jpj
@@ -309,11 +309,10 @@
       ! Vertically constant, constant in time
       !---------------------------------------
-      IF( num_sal == 1 )   s_i(:,:,:,:) = bulk_sal
+      IF(  num_sal == 1  )   s_i(:,:,:,:) = bulk_sal
 
       !-----------------------------------
       ! Salinity profile, varying in time
       !-----------------------------------
-
-      IF(   num_sal == 2  .OR.   num_sal == 4   ) THEN
+      IF(  num_sal == 2  ) THEN
          !
          DO jk = 1, nlay_i
@@ -331,5 +330,5 @@
          dummy_fac0 = 1._wp / ( s_i_0 - s_i_1 )       ! Weighting factor between zs_zero and zs_inf
          dummy_fac1 = s_i_1 / ( s_i_1 - s_i_0 )
-
+         !
          zalpha(:,:,:) = 0._wp
          DO jl = 1, jpl
@@ -347,5 +346,5 @@
             END DO
          END DO
-
+         !
          dummy_fac = 1._wp / nlay_i                   ! Computation of the profile
          DO jl = 1, jpl
@@ -361,5 +360,5 @@
             END DO ! jk
          END DO ! jl
-
+         !
       ENDIF ! num_sal
 
@@ -368,5 +367,5 @@
       !-------------------------------------------------------
 
-      IF( num_sal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
+      IF(  num_sal == 3  ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
          !
          sm_i(:,:,:) = 2.30_wp
@@ -380,5 +379,5 @@
             END DO
          END DO
-
+         !
       ENDIF ! num_sal
       !
@@ -447,5 +446,5 @@
       !------------------------------------------------------
 
-      IF(  num_sal == 2  .OR.  num_sal == 4  ) THEN
+      IF(  num_sal == 2  ) THEN
          !
          DO ji = kideb, kiut          ! Slope of the linear profile zs_zero
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/limwri_dimg.h90	(revision 3517)
@@ -111,5 +111,5 @@
          zcmo(ji,jj,13) = qns(ji,jj)
          ! See thersf for the coefficient
-         zcmo(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce
+         zcmo(ji,jj,14) = - sfx (ji,jj) * rday      ! converted in Kg/m2/day = mm/day
          zcmo(ji,jj,15) = utau_ice(ji,jj)
          zcmo(ji,jj,16) = vtau_ice(ji,jj)
@@ -154,5 +154,5 @@
                rcmoy(ji,jj,13) = qns(ji,jj)
                ! See thersf for the coefficient
-               rcmoy(ji,jj,14) = - sfx (ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce
+               rcmoy(ji,jj,14) = - sfx (ji,jj) * rday      ! converted in mm/day
                rcmoy(ji,jj,15) = utau_ice(ji,jj)
                rcmoy(ji,jj,16) = vtau_ice(ji,jj)
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 3517)
@@ -66,6 +66,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_b        !: <==> the 2D  frld
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fbif_1d       !: <==> the 2D  fbif
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdmicif_1d    !: <==> the 2D  rdmicif
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdmsnif_1d    !: <==> the 2D  rdmsnif
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdm_ice_1d    !: <==> the 2D  rdm_ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   rdm_snw_1d    !: <==> the 2D  rdm_snw
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qlbbq_1d      !: <==> the 2D  qlbsbq
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   dmgwi_1d      !: <==> the 2D  dmgwi
@@ -160,5 +160,5 @@
       !
       ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_b     (jpij) ,     &
-         &      fbif_1d    (jpij) , rdmicif_1d (jpij) , rdmsnif_1d (jpij) ,     &
+         &      fbif_1d    (jpij) , rdm_ice_1d (jpij) , rdm_snw_1d (jpij) ,     &
          &      qlbbq_1d   (jpij) , dmgwi_1d   (jpij) , dvsbq_1d   (jpij) ,     &
          &      dvbbq_1d   (jpij) , dvlbq_1d   (jpij) , dvnbq_1d   (jpij) ,     &
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 3517)
@@ -175,5 +175,5 @@
          IF( ierr3 > 0 )   CALL ctl_stop( 'STOP', 'sbc_blk_clio: unable to allocate arrays' )
          !
-         sfx(:,:) = 0.0_wp                      ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
+         sfx(:,:) = 0._wp                       ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
          !
       ENDIF
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 3517)
@@ -179,5 +179,5 @@
          CALL fld_fill( sf, slf_i, cn_dir, 'sbc_blk_core', 'flux formulation for ocean surface boundary condition', 'namsbc_core' )
          !
-         sfx(:,:) = 0.0_wp                         ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
+         sfx(:,:) = 0._wp                          ! salt flux; zero unless ice is present (computed in limsbc(_2).F90)
          !
       ENDIF
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 3517)
@@ -779,5 +779,5 @@
          ! Stress module can be negative when received (interpolation problem)
          IF( llnewtau ) THEN 
-            frcv(jpr_taum)%z3(:,:,1) = MAX( 0.0e0, frcv(jpr_taum)%z3(:,:,1) )
+            frcv(jpr_taum)%z3(:,:,1) = MAX( 0._wp, frcv(jpr_taum)%z3(:,:,1) )
          ENDIF
       ENDIF
@@ -823,5 +823,5 @@
          !                                                   ! ========================= !
          !
-         !                                                       ! total freshwater fluxes over the ocean (emp, sfx )
+         !                                                       ! total freshwater fluxes over the ocean (emp)
          SELECT CASE( TRIM( sn_rcv_emp%cldes ) )                                    ! evaporation - precipitation
          CASE( 'conservative' )
@@ -1261,5 +1261,5 @@
       !                                                      ! ========================= !
       CASE( 'oce only' )
-         qsr_tot(:,:  ) = MAX(0.0,frcv(jpr_qsroce)%z3(:,:,1))
+         qsr_tot(:,:  ) = MAX( 0._wp , frcv(jpr_qsroce)%z3(:,:,1) )
       CASE( 'conservative' )
          qsr_tot(:,:  ) = frcv(jpr_qsrmix)%z3(:,:,1)
@@ -1364,5 +1364,5 @@
             ztmp3(:,:,1:jpl) = tn_ice(:,:,1:jpl) * a_i(:,:,1:jpl)
          CASE( 'no' )
-            ztmp3(:,:,:) = 0.0
+            ztmp3(:,:,:) = 0._wp
             DO jl=1,jpl
                ztmp3(:,:,1) = ztmp3(:,:,1) + tn_ice(:,:,jl) * a_i(:,:,jl)
@@ -1416,5 +1416,5 @@
             ztmp4(:,:,1:jpl) =  ht_s(:,:,1:jpl) * a_i(:,:,1:jpl)
          CASE( 'no' )
-            ztmp3(:,:,:) = 0.0   ;  ztmp4(:,:,:) = 0.0
+            ztmp3(:,:,:) = 0._wp   ;  ztmp4(:,:,:) = 0._wp
             DO jl=1,jpl
                ztmp3(:,:,1) = ztmp3(:,:,1) + ht_i(:,:,jl) * a_i(:,:,jl)
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 3517)
@@ -10,5 +10,5 @@
    !!             -   ! 2008-04  (G. Madec)  sltyle and lim_ctl routine
    !!            3.3  ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
-   !!            4.0  ! 2011-01  (A Porter)  dynamical allocation
+   !!            3.4  ! 2011-01  (A Porter)  dynamical allocation
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -170,25 +170,25 @@
 
          !                                           ! intialisation to zero    !!gm is it truly necessary ???
-         d_a_i_thd  (:,:,:)   = 0.e0   ;   d_a_i_trp  (:,:,:)   = 0.e0
-         d_v_i_thd  (:,:,:)   = 0.e0   ;   d_v_i_trp  (:,:,:)   = 0.e0
-         d_e_i_thd  (:,:,:,:) = 0.e0   ;   d_e_i_trp  (:,:,:,:) = 0.e0
-         d_v_s_thd  (:,:,:)   = 0.e0   ;   d_v_s_trp  (:,:,:)   = 0.e0
-         d_e_s_thd  (:,:,:,:) = 0.e0   ;   d_e_s_trp  (:,:,:,:) = 0.e0
-         d_smv_i_thd(:,:,:)   = 0.e0   ;   d_smv_i_trp(:,:,:)   = 0.e0
-         d_oa_i_thd (:,:,:)   = 0.e0   ;   d_oa_i_trp (:,:,:)   = 0.e0
-         !
-         fseqv    (:,:) = 0.e0
-         fsbri    (:,:) = 0.e0     ;   fsalt_res(:,:) = 0.e0
-         fsalt_rpo(:,:) = 0.e0
-         fhmec    (:,:) = 0.e0     ;   fhbri    (:,:) = 0.e0
-         fmmec    (:,:) = 0.e0     ;   fheat_res(:,:) = 0.e0
-         fheat_rpo(:,:) = 0.e0     ;   focea2D  (:,:) = 0.e0
-         fsup2D   (:,:) = 0.e0
+         d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp
+         d_v_i_thd  (:,:,:)   = 0._wp   ;   d_v_i_trp  (:,:,:)   = 0._wp
+         d_e_i_thd  (:,:,:,:) = 0._wp   ;   d_e_i_trp  (:,:,:,:) = 0._wp
+         d_v_s_thd  (:,:,:)   = 0._wp   ;   d_v_s_trp  (:,:,:)   = 0._wp
+         d_e_s_thd  (:,:,:,:) = 0._wp   ;   d_e_s_trp  (:,:,:,:) = 0._wp
+         d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp
+         d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp
+         !
+         fseqv    (:,:) = 0._wp
+         fsbri    (:,:) = 0._wp     ;   fsalt_res(:,:) = 0._wp
+         fsalt_rpo(:,:) = 0._wp
+         fhmec    (:,:) = 0._wp     ;   fhbri    (:,:) = 0._wp
+         fmmec    (:,:) = 0._wp     ;   fheat_res(:,:) = 0._wp
+         fheat_rpo(:,:) = 0._wp     ;   focea2D  (:,:) = 0._wp
+         fsup2D   (:,:) = 0._wp
          ! 
-         diag_sni_gr(:,:) = 0.e0   ;   diag_lat_gr(:,:) = 0.e0
-         diag_bot_gr(:,:) = 0.e0   ;   diag_dyn_gr(:,:) = 0.e0
-         diag_bot_me(:,:) = 0.e0   ;   diag_sur_me(:,:) = 0.e0
+         diag_sni_gr(:,:) = 0._wp   ;   diag_lat_gr(:,:) = 0._wp
+         diag_bot_gr(:,:) = 0._wp   ;   diag_dyn_gr(:,:) = 0._wp
+         diag_bot_me(:,:) = 0._wp   ;   diag_sur_me(:,:) = 0._wp
          ! dynamical invariants
-         delta_i(:,:) = 0.e0       ;   divu_i(:,:) = 0.e0       ;   shear_i(:,:) = 0.e0
+         delta_i(:,:) = 0._wp       ;   divu_i(:,:) = 0._wp       ;   shear_i(:,:) = 0._wp
 
                           CALL lim_rst_opn( kt )     ! Open Ice restart file
@@ -196,6 +196,5 @@
          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 1, ' - Beginning the time step - ' )   ! control print
          !
-         IF( .NOT. lk_c1d ) THEN
-                                                     ! Ice dynamics & transport (not in 1D case)
+         IF( .NOT. lk_c1d ) THEN                     ! Ice dynamics & transport (except in 1D case)
                           CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics )
                           CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion )
@@ -210,5 +209,5 @@
                           CALL lim_var_bv                 ! bulk brine volume (diag)
                           CALL lim_thd( kt )              ! Ice thermodynamics 
-                          zcoef = rdt_ice / 86400.e0      !  Ice natural aging
+                          zcoef = rdt_ice /rday           !  Ice natural aging
                           oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * zcoef
                           CALL lim_var_glo2eqv            ! this CALL is maybe not necessary (Martin)
@@ -268,5 +267,5 @@
 
       inb_altests = 10
-      inb_alp(:)  = 0
+      inb_alp(:)  =  0
 
       ! Alert if incompatible volume and concentration
@@ -277,5 +276,5 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF(  v_i(ji,jj,jl) /= 0.e0   .AND.   a_i(ji,jj,jl) == 0.e0   ) THEN
+               IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN
                   WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration '
                   WRITE(numout,*) ' at_i     ', at_i(ji,jj)
@@ -297,5 +296,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF(   ht_i(ji,jj,jl) .GT. 50.0   ) THEN
+            IF(   ht_i(ji,jj,jl)  >  50._wp   ) THEN
                CALL lim_prt_state( ji, jj, 2, ' ALERTE 3 :   Very thick ice ' )
                inb_alp(ialert_id) = inb_alp(ialert_id) + 1
@@ -309,6 +308,6 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) .GT. 0.5  .AND.  &
-               &  at_i(ji,jj) .GT. 0.e0   ) THEN
+            IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 0.5  .AND.  &
+               &  at_i(ji,jj) > 0._wp   ) THEN
                CALL lim_prt_state( ji, jj, 1, ' ALERTE 4 :   Very fast ice ' )
                WRITE(numout,*) ' ice strength             : ', strength(ji,jj)
@@ -332,5 +331,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF(   tms(ji,jj) .LE. 0.0   .AND.   at_i(ji,jj) .GT. 0.e0   ) THEN 
+            IF(   tms(ji,jj) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN 
                CALL lim_prt_state( ji, jj, 1, ' ALERTE 6 :   Ice on continents ' )
                WRITE(numout,*) ' masks s, u, v        : ', tms(ji,jj), tmu(ji,jj), tmv(ji,jj) 
@@ -356,7 +355,7 @@
             DO ji = 1, jpi
 !!gm  test twice sm_i ...  ????  bug?
-               IF( ( ( ABS( sm_i(ji,jj,jl) ) .LT. 0.50) .OR. &
-                     ( ABS( sm_i(ji,jj,jl) ) .LT. 0.50) ) .AND. &
-                             ( a_i(ji,jj,jl) .GT. 0.e0 ) ) THEN
+               IF( ( ( ABS( sm_i(ji,jj,jl) ) < 0.5 )   .OR.  &
+                     ( ABS( sm_i(ji,jj,jl) ) < 0.5 ) ) .AND. &
+                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN
 !                 CALL lim_prt_state(ji,jj,1, ' ALERTE 7 :   Very fresh ice ' )
 !                 WRITE(numout,*) ' sst                  : ', sst_m(ji,jj)
@@ -377,7 +376,7 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               IF ( ( ( ABS( o_i(ji,jj,jl) ) .GT. rdt_ice ) .OR. &
-                      ( ABS( o_i(ji,jj,jl) ) .LT. 0.00) ) .AND. &
-                             ( a_i(ji,jj,jl) .GT. 0.0 ) ) THEN
+               IF ( ( ( ABS( o_i(ji,jj,jl) ) > rdt_ice ) .OR. &
+                      ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. &
+                             ( a_i(ji,jj,jl) > 0._wp ) ) THEN
                   CALL lim_prt_state( ji, jj, 1, ' ALERTE 9 :   Wrong ice age ')
                   inb_alp(ialert_id) = inb_alp(ialert_id) + 1
@@ -412,5 +411,5 @@
       DO jj = 1, jpj
          DO ji = 1, jpi
-            IF(   ABS( qns(ji,jj) ) .GT. 1500.0   .AND.  ( at_i(ji,jj) .GT. 0.0 ) )  THEN
+            IF(   ABS( qns(ji,jj) ) > 1500._wp  .AND.  at_i(ji,jj) > 0._wp  )  THEN
                !
                WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux'
@@ -450,6 +449,6 @@
                DO ji = 1, jpi
                   ztmelts    =  -tmut * s_i(ji,jj,jk,jl) + rtt
-                  IF( t_i(ji,jj,jk,jl) .GE. ztmelts  .AND.  v_i(ji,jj,jl) .GT. 1.e-6   &
-                     &                               .AND.  a_i(ji,jj,jl) .GT. 0.e0    ) THEN
+                  IF( t_i(ji,jj,jk,jl) >= ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-6   &
+                     &                             .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN
                      WRITE(numout,*) ' ALERTE 10 :   Very warm ice'
                      WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl
Index: /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 3516)
+++ /branches/2012/dev_r3385_NOCS04_HAMF/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 3517)
@@ -82,4 +82,5 @@
    END FUNCTION sbc_rnf_alloc
 
+
    SUBROUTINE sbc_rnf( kt )
       !!----------------------------------------------------------------------
@@ -95,5 +96,5 @@
       !!----------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt          ! ocean time step
-      !!
+      !
       INTEGER  ::   ji, jj   ! dummy loop indices
       !!----------------------------------------------------------------------
@@ -126,5 +127,6 @@
          !
          IF( MOD( kt - 1, nn_fsbc ) == 0 ) THEN
-            rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )  
+            !
+            rnf(:,:) = rn_rfact * ( sf_rnf(1)%fnow(:,:,1) )       ! updated runoff value at time step kt
             !
             !                                                     ! set temperature & salinity content of runoffs
@@ -248,10 +250,10 @@
       INTEGER           ::   ji, jj, jk    ! dummy loop indices
       INTEGER           ::   ierror, inum  ! temporary integer
-      !! 
+      !
       NAMELIST/namsbc_rnf/ cn_dir, ln_rnf_emp, ln_rnf_depth, ln_rnf_tem, ln_rnf_sal,   &
          &                 sn_rnf, sn_cnf    , sn_s_rnf    , sn_t_rnf  , sn_dep_rnf,   &  
          &                 ln_rnf_mouth      , rn_hrnf     , rn_avt_rnf, rn_rfact  
       !!----------------------------------------------------------------------
-
+      !
       !                                   ! ============
       !                                   !   Namelist
@@ -269,5 +271,5 @@
       REWIND ( numnam )                         ! Read Namelist namsbc_rnf
       READ   ( numnam, namsbc_rnf )
-
+      !
       !                                         ! Control print
       IF(lwp) THEN
@@ -282,5 +284,5 @@
          WRITE(numout,*) '      multiplicative factor for runoff           rn_rfact     = ', rn_rfact    
       ENDIF
-
+      !
       !                                   ! ==================
       !                                   !   Type of runoff
@@ -391,6 +393,5 @@
             nkrnf = 2
             DO WHILE( nkrnf /= jpkm1 .AND. gdepw_0(nkrnf+1) < rn_hrnf )   ;   nkrnf = nkrnf + 1   ;   END DO
-            IF( ln_sco )   &
-               CALL ctl_warn( 'sbc_rnf: number of levels over which Kz is increased is computed for zco...' )
+            IF( ln_sco )   CALL ctl_warn( 'sbc_rnf: number of levels over which Kz is increased is computed for zco...' )
          ENDIF
          IF(lwp) WRITE(numout,*)
@@ -410,5 +411,5 @@
          nkrnf = 0
       ENDIF
-
+      !
    END SUBROUTINE sbc_rnf_init
 
@@ -434,5 +435,4 @@
       !!                rnfmsk_z vertical structure
       !!----------------------------------------------------------------------
-      !
       INTEGER            ::   inum        ! temporary integers
       CHARACTER(len=140) ::   cl_rnfile   ! runoff file name
@@ -442,17 +442,17 @@
       IF(lwp) WRITE(numout,*) 'rnf_mouth : river mouth mask'
       IF(lwp) WRITE(numout,*) '~~~~~~~~~ '
-
+      !
       cl_rnfile = TRIM( cn_dir )//TRIM( sn_cnf%clname )
       IF( .NOT. sn_cnf%ln_clim ) THEN   ;   WRITE(cl_rnfile, '(a,"_y",i4)' ) TRIM( cl_rnfile ), nyear    ! add year
          IF( sn_cnf%cltype == 'monthly' )   WRITE(cl_rnfile, '(a,"m",i2)'  ) TRIM( cl_rnfile ), nmonth   ! add month
       ENDIF
-  
+      !
       ! horizontal mask (read in NetCDF file)
       CALL iom_open ( cl_rnfile, inum )                           ! open file
       CALL iom_get  ( inum, jpdom_data, sn_cnf%clvar, rnfmsk )    ! read the river mouth array
       CALL iom_close( inum )                                      ! close file
-      
+      !
       IF( nclosea == 1 )    CALL clo_rnf( rnfmsk )                ! closed sea inflow set as ruver mouth
-
+      !
       rnfmsk_z(:)   = 0._wp                                        ! vertical structure 
       rnfmsk_z(1)   = 1.0
