Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8426)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8486)
@@ -31,8 +31,19 @@
    !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
    !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
+   !!            3.5  ! 2012    (M. Vancoppenolle)  add ice_var_itd
+   !!            3.6  ! 2014-01 (C. Rousset) add ice_var_zapsmall, rewrite ice_var_itd
    !!----------------------------------------------------------------------
 #if defined key_lim3
    !!----------------------------------------------------------------------
    !!   'key_lim3'                                      LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   ice_var_agg       : integrate variables over layers and categories
+   !!   ice_var_glo2eqv   : transform from VGLO to VEQV
+   !!   ice_var_eqv2glo   : transform from VEQV to VGLO
+   !!   ice_var_salprof   : salinity profile in the ice
+   !!   ice_var_bv        : brine volume
+   !!   ice_var_salprof1d : salinity profile in the ice 1D
+   !!   ice_var_zapsmall  : remove very small area and volume
+   !!   ice_var_itd       : convert 1-cat to multiple cat
    !!----------------------------------------------------------------------
    USE par_oce        ! ocean parameters
@@ -59,5 +70,5 @@
 
    !!----------------------------------------------------------------------
-   !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011)
+   !! NEMO/ICE 4.0 , NEMO Consortium (2017)
    !! $Id: icevar.F90 8422 2017-08-08 13:58:05Z clem $
    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
@@ -69,84 +80,62 @@
       !!                ***  ROUTINE ice_var_agg  ***
       !!
-      !! ** Purpose :   aggregates ice-thickness-category variables to all-ice variables
-      !!              i.e. it turns VGLO into VAGG
-      !! ** Method  :
-      !!
-      !! ** Arguments : n = 1, at_i vt_i only
-      !!                n = 2 everything
-      !!
-      !! note : you could add an argument when you need only at_i, vt_i
-      !!        and when you need everything
-      !!------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kn     ! =1 at_i & vt only ; = what is needed
-      !
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      !!------------------------------------------------------------------
-
-      ! integrated values
-      vt_i (:,:) = SUM( v_i, dim=3 )
-      vt_s (:,:) = SUM( v_s, dim=3 )
-      at_i (:,:) = SUM( a_i, dim=3 )
+      !! ** Purpose :   aggregates ice-thickness-category variables to 
+      !!              all-ice variables, i.e. it turns VGLO into VAGG
+      !!------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kn     ! =1 state variables only
+      !                                 ! >1 state variables + others
+      !
+      INTEGER ::   ji, jj, jk, jl   ! dummy loop indices
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   z1_at_i, z1_vt_i
+      !!------------------------------------------------------------------
+      !
+      !                                      ! integrated values
+      vt_i(:,:) =       SUM( v_i(:,:,:)           , dim=3 )
+      vt_s(:,:) =       SUM( v_s(:,:,:)           , dim=3 )
+      at_i(:,:) =       SUM( a_i(:,:,:)           , dim=3 )
       et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 )
       et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 )
 
       ! MV MP 2016
-      IF ( ln_pnd ) THEN
-         at_ip(:,:) = SUM( a_ip, dim=3 )
-         vt_ip(:,:) = SUM( v_ip, dim=3 )
+      IF ( ln_pnd ) THEN                     ! Melt pond
+         at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 )
+         vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 )
       ENDIF
       ! END MP 2016
 
-      ! open water fraction
-      DO jj = 1, jpj
+      DO jj = 1, jpj                         ! open water fraction
          DO ji = 1, jpi
             ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )   
          END DO
       END DO
+!!gm I think this should do the work :
+!      ato_i(:,:) = MAX( 1._wp - at_i(:,:), 0._wp )  
+!!gm end
 
       IF( kn > 1 ) THEN
-
-         ! mean ice/snow thickness
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
-               htm_i(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
-               htm_s(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch
-            ENDDO
-         ENDDO
-
-         ! mean temperature (K), salinity and age
-         smt_i(:,:) = 0._wp
-         tm_i(:,:)  = 0._wp
-         tm_su(:,:) = 0._wp
-         tm_si(:,:) = 0._wp
-         om_i (:,:) = 0._wp
-         DO jl = 1, jpl
-            
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )
-                  tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
-                  tm_si(ji,jj) = tm_si(ji,jj) + rswitch * ( t_si(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 )
-                  om_i (ji,jj) = om_i (ji,jj) + rswitch *   oa_i(ji,jj,jl)                         / MAX( at_i(ji,jj) , epsi06 )
-               END DO
-            END DO
-            
-            DO jk = 1, nlay_i
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
-                     tm_i(ji,jj)  = tm_i(ji,jj)  + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  &
-                        &            / MAX( vt_i(ji,jj) , epsi10 )
-                     smt_i(ji,jj) = smt_i(ji,jj) + r1_nlay_i * rswitch * s_i(ji,jj,jk,jl) * v_i(ji,jj,jl)  &
-                        &            / MAX( vt_i(ji,jj) , epsi10 )
-                  END DO
-               END DO
-            END DO
-         END DO
-         tm_i  = tm_i  + rt0
-         tm_su = tm_su + rt0
-         tm_si = tm_si + rt0
-         !
+         !
+         ALLOCATE( z1_at_i(jpi,jpj) , z1_vt_i(jpi,jpj) )
+         WHERE( at_i(:,:) > epsi10 )   ;   z1_at_i(:,:) = 1._wp / at_i(:,:)
+         ELSEWHERE                     ;   z1_at_i(:,:) = 0._wp
+         END WHERE
+         WHERE( vt_i(:,:) > epsi10 )   ;   z1_vt_i(:,:) = 1._wp / vt_i(:,:)
+         ELSEWHERE                     ;   z1_vt_i(:,:) = 0._wp
+         END WHERE
+         !
+         !                          ! mean ice/snow thickness
+         htm_i(:,:) = vt_i(:,:) * z1_at_i(:,:)
+         htm_s(:,:) = vt_s(:,:) * z1_at_i(:,:)
+         !         
+         !                          ! mean temperature (K), salinity and age
+         tm_su(:,:) = SUM( t_su(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(ji,jj)
+         tm_si(:,:) = SUM( t_si(:,:,:) * a_i(:,:,:) , dim=3 ) * z1_at_i(ji,jj)
+         om_i (:,:) = SUM( oa_i(:,:,:)              , dim=3 ) * z1_at_i(ji,jj)
+         !
+         tm_i (:,:) = r1_nlay_i * SUM( SUM( t_i(:,:,:,:) * v_i(:,:,:), dim=4 ) , dim=3 ) * z1_vt_i(:,:)
+         smt_i(:,:) = r1_nlay_i * SUM( SUM( s_i(:,:,:,:) * v_i(:,:,:), dim=4 ) , dim=3 ) * z1_vt_i(:,:)
+!
+!!gm  QUESTION 1 : why salinity is named smt_i  and not just sm_i ?   since the 4D field is named s_i. (NB for temp: tm_i, t_i)
+         !
+         DEALLOCATE( z1_at_i , z1_vt_i )
       ENDIF
       !
@@ -158,53 +147,48 @@
       !!                ***  ROUTINE ice_var_glo2eqv ***
       !!
-      !! ** Purpose :   computes equivalent variables as function of global variables 
-      !!              i.e. it turns VGLO into VEQV
+      !! ** Purpose :   computes equivalent variables as function of  
+      !!              global variables, i.e. it turns VGLO into VEQV
       !!------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      REAL(wp) ::   ze_i, zaaa, zbbb, zccc, zdiscrim     ! local scalars
-      REAL(wp) ::   ztmelts, ze_s, zfac1, zfac2   !   -      -
-      !!------------------------------------------------------------------
+      REAL(wp) ::   ze_i, z1_CpR, zdiscrim, zaaa, z1_2aaa             ! local scalars
+      REAL(wp) ::   ze_s, zL_Cp , ztmelts , zbbb, zccc                !   -      -
+      REAL(wp) ::   zhmax, z1_zhmax, zsm_i, zcpMcpic, zt_i   !   -      -
+      REAL(wp) ::   zlay_i, zlay_s   !   -      -
+      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i, z1_v_i
+      !!------------------------------------------------------------------
+
+!!gm Question 1:  here use epsi20 , in ice_var_agg it is epsi10 which is used...  why ??
+
+!!gm Question 2:  It is possible to define existence of sea-ice in a common way between 
+!!                ice area and ice volume ?
+!!                the idea is to be able to define one for all at the begining of this routine
+!!                a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 )
 
       !-------------------------------------------------------
       ! Ice thickness, snow thickness, ice salinity, ice age
       !-------------------------------------------------------
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes
-               ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
-            END DO
-         END DO
-      END DO
-      ! Force the upper limit of ht_i to always be < hi_max (99 m).
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            rswitch = MAX( 0._wp , SIGN( 1._wp, ht_i(ji,jj,jpl) - epsi20 ) )
-            ht_i(ji,jj,jpl) = MIN( ht_i(ji,jj,jpl) , hi_max(jpl) )
-            a_i (ji,jj,jpl) = v_i(ji,jj,jpl) / MAX( ht_i(ji,jj,jpl) , epsi20 ) * rswitch
-         END DO
-      END DO
-
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes
-               ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
-               o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch
-            END DO
-         END DO
-      END DO
+      !                                            !--- inverse of the ice area
+      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
+      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp
+      END WHERE
+      !
+      ht_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)    !--- ice thickness
+
+	  zhmax    =          hi_max(jpl)
+      z1_zhmax =  1._wp / hi_max(jpl)               
+      WHERE( ht_i(:,:,jpl) > zhmax )               !--- bound ht_i by hi_max (i.e. 99 m) with associated update of ice area
+         ht_i  (:,:,jpl) = zhmax
+         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 
+         z1_a_i(:,:,jpl) = zhmax / v_i(:,:,jpl)          ! NB: v_i always /=0 as ht_i > hi_max
+      END WHERE
+
+      ht_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)    !--- snow thickness
       
-      IF(  nn_icesal == 2  )THEN
-         DO jl = 1, jpl
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi20 ) )   !0 if no ice and 1 if yes
-                  sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * rswitch
-                  !                                      ! bounding salinity
-                  sm_i(ji,jj,jl) = MAX( sm_i(ji,jj,jl), rn_simin )
-               END DO
-            END DO
-         END DO
+      o_i(:,:,:)  = oa_i(:,:,:) * z1_a_i(:,:,:)    !--- ice age
+
+      IF( nn_icesal == 2 )THEN                     !--- salinity (with a minimum value imposed everywhere)
+         WHERE( v_i(:,:,:) > epsi20 )   ;   sm_i(:,:,:) = MAX( rn_simin , smv_i(:,:,:) / v_i(:,:,:) )
+         ELSEWHERE                      ;   sm_i(:,:,:) = rn_simin
+         END WHERE
       ENDIF
 
@@ -212,21 +196,29 @@
 
       !-------------------
-      ! Ice temperatures
+      ! Ice temperature   [K]   (with a minimum value (rt0 - 100.) imposed everywhere)
       !-------------------
+      zlay_i   = REAL( nlay_i , wp )    ! number of layers
+      zaaa     = cpic                   ! Conversion q(S,T) -> T (second order equation)
+      z1_2aaa  = 1._wp / 2._wp *zaaa
+      zcpMcpic = rcp - cpic
       DO jl = 1, jpl
          DO jk = 1, nlay_i
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  !                                                              ! Energy of melting e(S,T) [J.m-3]
-                  rswitch = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes
-                  ze_i    = rswitch * e_i(ji,jj,jk,jl) / MAX( v_i(ji,jj,jl) , epsi20 ) * REAL(nlay_i,wp) 
-                  ztmelts = -tmut * s_i(ji,jj,jk,jl) + rt0              ! Ice layer melt temperature
-                  !
-                  zaaa       =  cpic                  ! Conversion q(S,T) -> T (second order equation)
-                  zbbb       =  ( rcp - cpic ) * ( ztmelts - rt0 ) + ze_i * r1_rhoic - lfus
-                  zccc       =  lfus * (ztmelts-rt0)
-                  zdiscrim   =  SQRT( MAX(zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp) )
-                  t_i(ji,jj,jk,jl) = rt0 + rswitch *( - zbbb - zdiscrim ) / ( 2.0 *zaaa )
-                  t_i(ji,jj,jk,jl) = MIN( ztmelts, MAX( rt0 - 100._wp, t_i(ji,jj,jk,jl) ) )  ! -100 < t_i < ztmelts
+                  IF ( v_i(ji,jj,jl) > epsi20 ) THEN     !--- icy area 
+                     !
+                     ze_i    =   e_i(ji,jj,jk,jl) / v_i(ji,jj,jl) * zlay_i               ! Energy of melting e(S,T) [J.m-3]
+                     ztmelts = - s_i(ji,jj,jk,jl) * tmut                                 ! Ice layer melt temperature [C]
+                     !
+                     zbbb     =  zcpMcpic * ztmelts + ze_i * r1_rhoic - lfus
+                     zccc     =  lfus * ztmelts
+                     zdiscrim =  SQRT(  MAX( zbbb*zbbb - 4._wp*zaaa*zccc , 0._wp)  )
+                     zt_i     = - ( zbbb + zdiscrim ) * z1_2aaa                          ! [C]
+                     t_i(ji,jj,jk,jl) = MAX( -100._wp , MIN( zt_i , ztmelts )  ) + rt0   ! [K] with bounds: -100 < zt_i < ztmelts
+                     !
+                  ELSE                                !--- no ice
+                     t_i(ji,jj,jk,jl) =  rt0 - 100._wp                                   ! impose 173.15 K (i.e. -100 C)
+!!clem: I think we should impose rt0 instead
+                  ENDIF
                END DO
             END DO
@@ -235,22 +227,27 @@
 
       !--------------------
-      ! Snow temperatures
+      ! Snow temperature   [K]   (with a minimum value (rt0 - 100.) imposed everywhere)
       !--------------------
-      zfac1 = 1._wp / ( rhosn * cpic )
-      zfac2 = lfus / cpic  
-      DO jl = 1, jpl
-         DO jk = 1, nlay_s
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  !Energy of melting q(S,T) [J.m-3]
-                  rswitch = MAX( 0._wp , SIGN( 1._wp , v_s(ji,jj,jl) - epsi20 ) )     ! rswitch = 0 if no ice and 1 if yes
-                  ze_s  = rswitch * e_s(ji,jj,jk,jl) / MAX( v_s(ji,jj,jl) , epsi20 ) * REAL(nlay_s,wp)
-                  !
-                  t_s(ji,jj,jk,jl) = rt0 + rswitch * ( - zfac1 * ze_s + zfac2 )
-                  t_s(ji,jj,jk,jl) = MIN( rt0, MAX( rt0 - 100._wp , t_s(ji,jj,jk,jl) ) )     ! -100 < t_s < rt0
-               END DO
-            END DO
-         END DO
+      z1_CpR = 1._wp / ( cpic * rhosn )
+      zL_Cp  = lfus  /   cpic
+      zlay_s = REAL( nlay_s , wp )
+      DO jk = 1, nlay_s
+         WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
+            t_s(:,:,jk,:) = MAX(  -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp )  ) + rt0       
+         ELSEWHERE                           !--- no ice
+            t_s(:,:,jk,:) =  rt0 - 100._wp         ! impose 173.15 K (i.e. -100 C)
+         END WHERE
       END DO
+!!gm perhaps better like this (?) :  (if this compile, yes! one test instead of nlay_s tests)
+!      WHERE( v_s(:,:,:) > epsi20 )        !--- icy area
+!         DO jk = 1, nlay_s
+!            t_s(:,:,jk,:) = MAX(  -100._wp , MIN( - z1_CpR * ( e_s(:,:,jk,:)/v_s(:,:,:)*zlay_s ) + zL_Cp , 0._wp )  ) + rt0       
+!         END DO
+!      ELSEWHERE                           !--- no ice
+!         DO jk = 1, nlay_s
+!            t_s(:,:,jk,:) =  rt0 - 100._wp      ! impose 173.15 K (i.e. -100 C)
+!         END DO
+!      END WHERE
+!!gm end better ?
 
       ! integrated values 
@@ -259,6 +256,6 @@
       at_i (:,:) = SUM( a_i, dim=3 )
 
-      ! MV MP 2016
-      ! probably should resum for melt ponds ???
+! MV MP 2016
+! probably should resum for melt ponds ???
 
       !
@@ -270,13 +267,10 @@
       !!                ***  ROUTINE ice_var_eqv2glo ***
       !!
-      !! ** Purpose :   computes global variables as function of equivalent variables
-      !!                i.e. it turns VEQV into VGLO
-      !! ** Method  :
-      !!
-      !! ** History :  (01-2006) Martin Vancoppenolle, UCL-ASTR
-      !!------------------------------------------------------------------
-      !
-      v_i(:,:,:)   = ht_i(:,:,:) * a_i(:,:,:)
-      v_s(:,:,:)   = ht_s(:,:,:) * a_i(:,:,:)
+      !! ** Purpose :   computes global variables as function of 
+      !!              equivalent variables,  i.e. it turns VEQV into VGLO
+      !!------------------------------------------------------------------
+      !
+      v_i  (:,:,:) = ht_i(:,:,:) * a_i(:,:,:)
+      v_s  (:,:,:) = ht_s(:,:,:) * a_i(:,:,:)
       smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
       !
@@ -300,53 +294,45 @@
       !!------------------------------------------------------------------
       INTEGER  ::   ji, jj, jk, jl   ! dummy loop index
-      REAL(wp) ::   zfac0, zfac1, zsal
-      REAL(wp) ::   zswi0, zswi01, zargtemp , zs_zero   
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z_slope_s, zalpha
+      REAL(wp) ::   zsal, z1_dS
+      REAL(wp) ::   zargtemp , zs0, zs
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   z_slope_s, zalpha    ! case 2 only
       REAL(wp), PARAMETER :: zsi0 = 3.5_wp
       REAL(wp), PARAMETER :: zsi1 = 4.5_wp
       !!------------------------------------------------------------------
 
-      !---------------------------------------
-      ! Vertically constant, constant in time
-      !---------------------------------------
-      IF(  nn_icesal == 1  )  THEN
+!!gm  much much more secure to defined when reading nn_icesal in the namelist integers =1, 2, 3  with explicit names
+!!       for example np_Scst_noProfile = 1  ;  np_Svar_linProfile = 2   ;   np_Scst_fixProfile
+
+!!gm Question: Remove the option 3 ?  How many years since it last use ? 
+
+      SELECT CASE ( nn_icesal )
+      !
+      !              !---------------------------------------!
+      CASE( 1 )      !  constant salinity in time and space  !
+         !           !---------------------------------------!
          s_i (:,:,:,:) = rn_icesal
          sm_i(:,:,:)   = rn_icesal
-      ENDIF
-
-      !-----------------------------------
-      ! Salinity profile, varying in time
-      !-----------------------------------
-      IF(  nn_icesal == 2  ) THEN
+         !
+         !           !---------------------------------------------!
+      CASE( 2 )      !  time varying salinity with linear profile  !
+         !           !---------------------------------------------!
+         !
+         ALLOCATE( z_slope_s(jpi,jpj,jpl) , zalpha(jpi,jpj,jpl) )
          !
          DO jk = 1, nlay_i
             s_i(:,:,jk,:)  = sm_i(:,:,:)
          END DO
-         !
-         DO jl = 1, jpl                               ! Slope of the linear profile 
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i(ji,jj,jl) - epsi20 ) )
-                  z_slope_s(ji,jj,jl) = rswitch * 2._wp * sm_i(ji,jj,jl) / MAX( epsi20 , ht_i(ji,jj,jl) )
-               END DO
-            END DO
-         END DO
-         !
-         zfac0 = 1._wp / ( zsi0 - zsi1 )       ! Weighting factor between zs_zero and zs_inf
-         zfac1 = zsi1  / ( zsi1 - zsi0 )
-         !
-         zalpha(:,:,:) = 0._wp
+         !                                      ! Slope of the linear profile 
+         WHERE( ht_i(:,:,:) > epsi20 )   ;   z_slope_s(:,:,:) = 2._wp * sm_i(:,:,:) / ht_i(:,:,:)
+         ELSEWHERE                       ;   z_slope_s(:,:,:) = 0._wp
+         END WHERE
+         !
+         z1_dS = 1._wp / ( zsi1 - zsi0 )
          DO jl = 1, jpl
             DO jj = 1, jpj
                DO ji = 1, jpi
-                  ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
-                  zswi0  = MAX( 0._wp   , SIGN( 1._wp  , zsi0 - sm_i(ji,jj,jl) ) ) 
-                  ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 
-                  zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp   , SIGN( 1._wp  , zsi1 - sm_i(ji,jj,jl) ) ) 
-                  ! If 2.sm_i GE sss_m then rswitch = 1
-                  ! this is to force a constant salinity profile in the Baltic Sea
-                  rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i(ji,jj,jl) - sss_m(ji,jj) ) )
-                  zalpha(ji,jj,jl) = zswi0  + zswi01 * ( sm_i(ji,jj,jl) * zfac0 + zfac1 )
-                  zalpha(ji,jj,jl) = zalpha(ji,jj,jl) * ( 1._wp - rswitch )
+                  zalpha(ji,jj,jl) = MAX(  0._wp , MIN( ( zsi1 - sm_i(ji,jj,jl) ) * z1_dS , 1._wp )  )
+                  !                             ! force a constant profile when SSS too low (Baltic Sea)
+                  IF( 2._wp * sm_i(ji,jj,jl) >= sss_m(ji,jj) )   zalpha(ji,jj,jl) = 0._wp  
                END DO
             END DO
@@ -358,10 +344,8 @@
                DO jj = 1, jpj
                   DO ji = 1, jpi
-                     !                                      ! linear profile with 0 at the surface
-                     zs_zero = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
-                     !                                      ! weighting the profile
-                     s_i(ji,jj,jk,jl) = zalpha(ji,jj,jl) * zs_zero + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)
-                     !                                      ! bounding salinity
-                     s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( s_i(ji,jj,jk,jl), rn_simin ) )
+                     !                          ! linear profile with 0 surface value
+                     zs0 = z_slope_s(ji,jj,jl) * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) * r1_nlay_i
+                     zs  = zalpha(ji,jj,jl) * zs0 + ( 1._wp - zalpha(ji,jj,jl) ) * sm_i(ji,jj,jl)     ! weighting the profile
+                     s_i(ji,jj,jk,jl) = MIN( rn_simax, MAX( zs, rn_simin ) )
                   END DO
                END DO
@@ -369,23 +353,26 @@
          END DO
          !
-      ENDIF ! nn_icesal
-
-      !-------------------------------------------------------
-      ! Vertically varying salinity profile, constant in time
-      !-------------------------------------------------------
-
-      IF(  nn_icesal == 3  ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
+         DEALLOCATE( z_slope_s , zalpha )
+         !
+         !           !-------------------------------------------!
+      CASE( 3 )      ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
+         !           !-------------------------------------------!                                   (mean = 2.30)
          !
          sm_i(:,:,:) = 2.30_wp
+!!gm Remark: if we keep the case 3, then compute an store one for all time-step
+!!           a array  S_prof(1:nlay_i) containing the calculation and just do:
+!         DO jk = 1, nlay_i
+!            s_i(:,:,jk,:) = S_prof(jk)
+!         END DO
+!!gm end
          !
          DO jl = 1, jpl
             DO jk = 1, nlay_i
                zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
-               zsal =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
-               s_i(:,:,jk,jl) =  zsal
-            END DO
-         END DO
-         !
-      ENDIF ! nn_icesal
+               s_i(:,:,jk,jl) =  1.6_wp * (  1._wp - COS( rpi * zargtemp**(0.407_wp/(0.573_wp+zargtemp)) )  )
+            END DO
+         END DO
+         !
+      END SELECT
       !
    END SUBROUTINE ice_var_salprof
@@ -405,25 +392,18 @@
       !!------------------------------------------------------------------
       !
-      bvm_i(:,:)   = 0._wp
+!!gm I prefere to use WHERE / ELSEWHERE  to set it to zero only where needed   <<<=== to be done
+!!   instead of setting everything to zero as just below
       bv_i (:,:,:) = 0._wp
       DO jl = 1, jpl
          DO jk = 1, nlay_i
-            DO jj = 1, jpj
-               DO ji = 1, jpi
-                  rswitch        = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  )
-                  bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i  &
-                     &                            / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )
-               END DO
-            END DO
-         END DO
-         
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               rswitch      = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) )
-               bvm_i(ji,jj) = bvm_i(ji,jj) + rswitch * bv_i(ji,jj,jl) * v_i(ji,jj,jl) / MAX( vt_i(ji,jj), epsi10 )
-            END DO
+            WHERE( t_i(:,:,jk,jl) < rt0 - epsi10 )   
+               bv_i(:,:,jl) = bv_i(:,:,jl) - tmut * s_i(:,:,jk,jl) * r1_nlay_i / ( t_i(:,:,jk,jl) - rt0 )
+            END WHERE
          END DO
       END DO
-      !
+      WHERE( vt_i(:,:) > epsi20 )   bvm_i(:,:) = SUM( bv_i(:,:,:) * v_i(:,:,:) , dim=3 ) / vt_i(:,:)
+      ELSEWHERE                     bvm_i(:,:) = 0._wp
+	  END WHERE
+	  !
    END SUBROUTINE ice_var_bv
 
@@ -437,62 +417,51 @@
       !!-------------------------------------------------------------------
       INTEGER  ::   ji, jk    ! dummy loop indices
-      REAL(wp) ::   zfac0, zfac1, zargtemp, zsal   ! local scalars
-      REAL(wp) ::   zalpha, zswi0, zswi01, zs_zero              !   -      -
-      !
-      REAL(wp), DIMENSION(jpij) ::   z_slope_s
+      REAL(wp) ::   zargtemp, zsal, z1_dS   ! local scalars
+      REAL(wp) ::   zalpha, zs, zs0              !   -      -
+      !
+      REAL(wp), ALLOCATABLE, DIMENSION(:) ::   z_slope_s   !
       REAL(wp), PARAMETER :: zsi0 = 3.5_wp
       REAL(wp), PARAMETER :: zsi1 = 4.5_wp
       !!---------------------------------------------------------------------
-
-      !---------------------------------------
-      ! Vertically constant, constant in time
-      !---------------------------------------
-      IF( nn_icesal == 1 )   s_i_1d(:,:) = rn_icesal
-
-      !------------------------------------------------------
-      ! Vertically varying salinity profile, varying in time
-      !------------------------------------------------------
-
-      IF(  nn_icesal == 2  ) THEN
-         !
-         DO ji = 1, nidx          ! Slope of the linear profile zs_zero
+      !
+      SELECT CASE ( nn_icesal )
+      !
+      !              !---------------------------------------!
+      CASE( 1 )      !  constant salinity in time and space  !
+         !           !---------------------------------------!
+         s_i_1d(:,:) = rn_icesal
+         !
+         !           !---------------------------------------------!
+      CASE( 2 )      !  time varying salinity with linear profile  !
+         !           !---------------------------------------------!
+         !
+         ALLOCATE( z_slope_s(jpij) )
+         !
+         DO ji = 1, nidx          ! Slope of the linear profile
             rswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_1d(ji) - epsi20 ) )
             z_slope_s(ji) = rswitch * 2._wp * sm_i_1d(ji) / MAX( epsi20 , ht_i_1d(ji) )
          END DO
 
-         ! Weighting factor between zs_zero and zs_inf
-         !---------------------------------------------
-         zfac0 = 1._wp / ( zsi0 - zsi1 )
-         zfac1 = zsi1 / ( zsi1 - zsi0 )
+         z1_dS = 1._wp / ( zsi1 - zsi0 )
          DO jk = 1, nlay_i
             DO ji = 1, nidx
-               ! zswi0 = 1 if sm_i le zsi0 and 0 otherwise
-               zswi0  = MAX( 0._wp , SIGN( 1._wp  , zsi0 - sm_i_1d(ji) ) ) 
-               ! zswi01 = 1 if sm_i is between zsi0 and zsi1 and 0 othws 
-               zswi01 = ( 1._wp - zswi0 ) * MAX( 0._wp , SIGN( 1._wp , zsi1 - sm_i_1d(ji) ) ) 
-               ! if 2.sm_i GE sss_m then rswitch = 1
-               ! this is to force a constant salinity profile in the Baltic Sea
-               rswitch = MAX( 0._wp , SIGN( 1._wp , 2._wp * sm_i_1d(ji) - sss_1d(ji) ) )
+               zalpha = MAX(  0._wp , MIN(  ( zsi1 - sm_i_1d(ji) ) * z1_dS , 1._wp  )  )
+               IF( 2._wp * sm_i_1d(ji) >= sss_1d(ji) )   zalpha = 0._wp               ! cst profile when SSS too low (Baltic Sea)
                !
-               zalpha = (  zswi0 + zswi01 * ( sm_i_1d(ji) * zfac0 + zfac1 )  ) * ( 1._wp - rswitch )
-               !
-               zs_zero = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i
-               ! weighting the profile
-               s_i_1d(ji,jk) = zalpha * zs_zero + ( 1._wp - zalpha ) * sm_i_1d(ji)
-               ! bounding salinity
-               s_i_1d(ji,jk) = MIN( rn_simax, MAX( s_i_1d(ji,jk), rn_simin ) )
-            END DO 
-         END DO 
-
-      ENDIF 
-
-      !-------------------------------------------------------
-      ! Vertically varying salinity profile, constant in time
-      !-------------------------------------------------------
-
-      IF( nn_icesal == 3 ) THEN      ! Schwarzacher (1959) multiyear salinity profile (mean = 2.30)
+               zs0 = z_slope_s(ji) * ( REAL(jk,wp) - 0.5_wp ) * ht_i_1d(ji) * r1_nlay_i   ! linear profile with 0 surface value
+               zs  = zalpha * zs0 + ( 1._wp - zalpha ) * sm_i_1d(ji)                      ! weighting the profile
+               s_i_1d(ji,jk) = MIN( rn_simax , MAX( zs , rn_simin ) )                     ! bounding salinity
+            END DO
+         END DO
+         !
+         DEALLOCATE( z_slope_s )
+
+         !           !-------------------------------------------!
+      CASE( 3 )      ! constant salinity with a fix profile      ! (Schwarzacher (1959) multiyear salinity profile
+         !           !-------------------------------------------!                                   (mean = 2.30)
          !
          sm_i_1d(:) = 2.30_wp
          !
+!!gm cf remark in ice_var_salprof routine, CASE( 3 )
          DO jk = 1, nlay_i
             zargtemp  = ( REAL(jk,wp) - 0.5_wp ) * r1_nlay_i
@@ -503,7 +472,8 @@
          END DO
          !
-      ENDIF
+      END SELECT
       !
    END SUBROUTINE ice_var_salprof1d
+
 
    SUBROUTINE ice_var_zapsmall
@@ -512,12 +482,11 @@
       !!
       !! ** Purpose :   Remove too small sea ice areas and correct fluxes
-      !!
-      !! history : LIM3.5 - 01-2014 (C. Rousset) original code
       !!-------------------------------------------------------------------
       INTEGER  ::   ji, jj, jl, jk   ! dummy loop indices
       REAL(wp) ::   zsal, zvi, zvs, zei, zes, zvp
       !!-------------------------------------------------------------------
-      DO jl = 1, jpl
-
+      !
+      DO jl = 1, jpl       !==  loop over the categories  ==!
+         !
          !-----------------------------------------------------------------
          ! Zap ice energy and use ocean heat to melt ice
@@ -526,4 +495,5 @@
             DO jj = 1 , jpj
                DO ji = 1 , jpi
+!!gm I think we can do better/faster  :  to be discussed
                   rswitch          = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi10 ) )
                   rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) - epsi10 ) ) * rswitch
@@ -588,10 +558,12 @@
 
       ! to be sure that at_i is the sum of a_i(jl)
-      at_i (:,:) = 0._wp
-      DO jl = 1, jpl
+      at_i (:,:) = a_i(:,:,1)
+	  vt_i (:,:) = v_i(:,:,1)
+      DO jl = 2, jpl
          at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
+	     vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)
       END DO
 
-      ! open water = 1 if at_i=0
+      ! open water = 1 if at_i=0 (no re-calculation of ato_i here)
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -600,7 +572,7 @@
          END DO
       END DO
-
       !
    END SUBROUTINE ice_var_zapsmall
+
 
    SUBROUTINE ice_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
@@ -633,10 +605,5 @@
       !!
       !!  (Example of application: BDY forcings when input are cell averaged)  
-      !!
       !!-------------------------------------------------------------------
-      !! History : LIM3.5 - 2012    (M. Vancoppenolle)  Original code
-      !!                    2014    (C. Rousset)        Rewriting
-      !!-------------------------------------------------------------------
-      !! Local variables
       INTEGER  :: ji, jk, jl             ! dummy loop indices
       INTEGER  :: ijpij, i_fill, jl0  
@@ -645,9 +612,10 @@
       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   zht_i, zht_s, za_i ! output ice/snow variables
       INTEGER , DIMENSION(4)                  ::   itest
- 
+      !!-------------------------------------------------------------------
+
       !--------------------------------------------------------------------
       ! initialisation of variables
       !--------------------------------------------------------------------
-      ijpij = SIZE(zhti,1)
+      ijpij = SIZE( zhti , 1 )
       zht_i(1:ijpij,1:jpl) = 0._wp
       zht_s(1:ijpij,1:jpl) = 0._wp
@@ -766,9 +734,8 @@
                zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
             ENDIF
-         ENDDO
-      ENDDO
+         END DO
+      END DO
       !
     END SUBROUTINE ice_var_itd
-
 
 #else
@@ -776,7 +743,4 @@
    !!   Default option         Dummy module          NO  LIM3 sea-ice model
    !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE ice_var_itd
-   END SUBROUTINE ice_var_itd
 #endif
 
