Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8424)
@@ -63,12 +63,12 @@
    !!   meaningful and/or used in ice thermodynamics                      |
    !!                                                                     |
-   !! Routines in limvar.F90 perform conversions                          |
-   !!  - lim_var_glo2eqv  : from global to equivalent variables           |
-   !!  - lim_var_eqv2glo  : from equivalent to global variables           |
+   !! Routines in icevar.F90 perform conversions                          |
+   !!  - ice_var_glo2eqv  : from global to equivalent variables           |
+   !!  - ice_var_eqv2glo  : from equivalent to global variables           |
    !!                                                                     |
    !! For various purposes, the sea ice state variables have sometimes    |
    !! to be aggregated over all ice thickness categories. This operation  |
    !! is done in :                                                        |
-   !!  - lim_var_agg                                                      |
+   !!  - ice_var_agg                                                      |
    !!                                                                     |
    !! in icestp.F90, the routines that compute the changes in the ice     |
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceadv.F90	(revision 8424)
@@ -18,5 +18,5 @@
    USE sbc_oce , ONLY : nn_fsbc
    USE ice            ! ice variables
-   USE limvar         ! 
+   USE icevar         ! 
    USE iceadv_prather ! advection scheme (Prather)
    USE iceadv_umx     ! advection scheme (ultimate-macho)
@@ -96,5 +96,5 @@
       ENDIF
       
-      CALL lim_var_agg( 1 ) ! integrated values + ato_i
+      CALL ice_var_agg( 1 ) ! integrated values + ato_i
 
       !-------------------------------------!
@@ -414,5 +414,5 @@
 
          ! zap small areas
-         CALL lim_var_zapsmall
+         CALL ice_var_zapsmall
            
          !--- Thickness correction in case too high --- !
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceerr1.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceerr1.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceerr1.F90	(revision 8424)
@@ -18,5 +18,5 @@
    USE ice1D           ! LIM thermodynamic sea-ice variables
    USE iceitd
-   USE limvar
+   USE icevar
    USE icecons         ! conservation tests
    USE icectl          ! control prints
@@ -109,5 +109,5 @@
       ! zap small values
       !-----------------
-      CALL lim_var_zapsmall
+      CALL ice_var_zapsmall
 
       ! -------------------------------------------------
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceerr2.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceerr2.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceerr2.F90	(revision 8424)
@@ -18,5 +18,5 @@
    USE ice1D           ! LIM thermodynamic sea-ice variables
    USE iceitd
-   USE limvar
+   USE icevar
    USE icecons         ! conservation tests
    USE icectl
@@ -123,5 +123,5 @@
       ! zap small values
       !-----------------
-      CALL lim_var_zapsmall
+      CALL ice_var_zapsmall
 
       !------------------------------------------------------------------------------
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90	(revision 8424)
@@ -23,5 +23,5 @@
    USE ice            ! sea-ice variables
    USE par_oce        ! ocean parameters
-   USE limvar         ! lim_var_salprof
+   USE icevar         ! ice_var_salprof
    !
    USE in_out_manager ! I/O manager
@@ -324,5 +324,5 @@
          ! for constant salinity in time
          IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN
-            CALL lim_var_salprof
+            CALL ice_var_salprof
             smv_i = sm_i * v_i
          ENDIF
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90	(revision 8424)
@@ -23,5 +23,4 @@
    USE ice1D            ! LIM-3 thermodynamic variables
    USE ice              ! LIM-3 variables
-   USE limvar           ! LIM-3 variables
    USE icecons          ! conservation tests
    USE icetab
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerdgrft.F90	(revision 8424)
@@ -18,5 +18,5 @@
    USE ice1D            ! LIM thermodynamics
    USE ice              ! LIM variables
-   USE limvar           ! LIM
+   USE icevar           ! LIM
    USE icecons          ! conservation tests
    USE icectl           ! control prints
@@ -275,5 +275,5 @@
       END DO !! on the do while over iter
 
-      CALL lim_var_agg( 1 ) 
+      CALL ice_var_agg( 1 ) 
 
       !-----------------------------------------------------------------------------!
@@ -311,5 +311,5 @@
 
       ! Zero out categories with very small areas
-      CALL lim_var_zapsmall
+      CALL ice_var_zapsmall
 
       ! Ice thickness needed for rafting
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90	(revision 8424)
@@ -21,5 +21,5 @@
    USE icecons          ! conservation tests
    USE icectl           ! control prints
-   USE limvar
+   USE icevar
    !
    USE lbclnk           ! lateral boundary conditions - MPP exchanges
@@ -62,5 +62,5 @@
       IF( nn_timing == 1 )  CALL timing_start('icerhg')
 
-      CALL lim_var_agg(1)   ! aggregate ice categories
+      CALL ice_var_agg(1)   ! aggregate ice categories
       !
       ! conservation test
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8424)
@@ -42,5 +42,5 @@
    USE iceerr1         ! Ice corrections after dynamics
    USE iceerr2         ! Ice corrections after thermo
-   USE limvar          ! Ice variables switch
+   USE icevar          ! Ice variables switch
    USE icectl          !
    ! MV MP 2016
@@ -167,6 +167,6 @@
          IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice( kt )       ! -- bdy ice thermo 
          ! previous lead fraction and ice volume for flux calculations
-                                      CALL lim_var_glo2eqv     ! ht_i and ht_s for ice albedo calculation
-                                      CALL lim_var_agg(1)      ! at_i for coupling 
+                                      CALL ice_var_glo2eqv     ! ht_i and ht_s for ice albedo calculation
+                                      CALL ice_var_agg(1)      ! at_i for coupling 
                                       CALL ice_bef
 
@@ -201,6 +201,6 @@
          IF( .NOT. Agrif_Root() )     CALL agrif_update_lim3( kt )
 # endif
-                                      CALL lim_var_glo2eqv      ! necessary calls (at least for coupling)
-                                      CALL lim_var_agg( 2 )     ! necessary calls (at least for coupling)
+                                      CALL ice_var_glo2eqv      ! necessary calls (at least for coupling)
+                                      CALL ice_var_agg( 2 )     ! necessary calls (at least for coupling)
                                       !
 !! clem: one should switch the calculation of the fluxes onto the parent grid but the following calls do not work
@@ -286,6 +286,6 @@
          CALL ice_rst_read
       ENDIF
-      CALL lim_var_agg(2)
-      CALL lim_var_glo2eqv
+      CALL ice_var_agg(2)
+      CALL ice_var_glo2eqv
       !
       CALL ice_update_init                 ! ice surface boundary condition
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90	(revision 8424)
@@ -35,5 +35,5 @@
    USE iceitd         ! remapping thickness distribution
    USE icetab         ! 1D <==> 2D transformation
-   USE limvar         !
+   USE icevar         !
    USE icecons        ! conservation tests
    USE icectl         ! control print
@@ -101,5 +101,5 @@
       IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'icethd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-      CALL lim_var_glo2eqv
+      CALL ice_var_glo2eqv
 
       !---------------------------------------------!
@@ -259,5 +259,5 @@
       IF( ln_limdiachk ) CALL ice_cons_hsm(1, 'icethd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
       !
-                         CALL lim_var_zapsmall           ! --- remove very small ice concentration (<1e-10) --- !
+                         CALL ice_var_zapsmall           ! --- remove very small ice concentration (<1e-10) --- !
       !                                                  !     & make sure at_i=SUM(a_i) & ato_i=1 where at_i=0
       !                   
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_ent.F90	(revision 8424)
@@ -25,5 +25,4 @@
    USE ice            ! LIM variables
    USE ice1D          ! LIM thermodynamics
-   USE limvar         ! LIM variables
    !
    USE in_out_manager ! I/O manager
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_lac.F90	(revision 8424)
@@ -25,5 +25,5 @@
    USE icecons        ! LIM conservation
    USE icethd_ent
-   USE limvar
+   USE icevar
    !
    USE in_out_manager ! I/O manager
@@ -114,6 +114,6 @@
       IF( ln_limdiachk ) CALL ice_cons_hsm(0, 'icethd_lac', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
 
-      CALL lim_var_agg(1)
-      CALL lim_var_glo2eqv
+      CALL ice_var_agg(1)
+      CALL ice_var_glo2eqv
 
       !------------------------------------------------------------------------------!
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icethd_sal.F90	(revision 8424)
@@ -18,5 +18,5 @@
    USE ice            ! LIM variables
    USE ice1D          ! LIM thermodynamics
-   USE limvar         ! LIM variables
+   USE icevar         ! LIM variables
    !
    USE in_out_manager ! I/O manager
@@ -94,5 +94,5 @@
 
          ! Salinity profile
-         CALL lim_var_salprof1d
+         CALL ice_var_salprof1d
          !
       ENDIF 
@@ -101,5 +101,5 @@
       !  3) vertical profile of salinity, constant in time                           |
       !------------------------------------------------------------------------------|
-      IF(  nn_icesal == 3  )   CALL lim_var_salprof1d
+      IF(  nn_icesal == 3  )   CALL ice_var_salprof1d
 
       !
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 8424)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90	(revision 8424)
@@ -0,0 +1,784 @@
+MODULE icevar
+   !!======================================================================
+   !!                       ***  MODULE icevar ***
+   !!                 Different sets of ice model variables 
+   !!                   how to switch from one to another
+   !!
+   !!                 There are three sets of variables
+   !!                 VGLO : global variables of the model
+   !!                        - v_i (jpi,jpj,jpl)
+   !!                        - v_s (jpi,jpj,jpl)
+   !!                        - a_i (jpi,jpj,jpl)
+   !!                        - t_s (jpi,jpj,jpl)
+   !!                        - e_i (jpi,jpj,nlay_i,jpl)
+   !!                        - smv_i(jpi,jpj,jpl)
+   !!                        - oa_i (jpi,jpj,jpl)
+   !!                 VEQV : equivalent variables sometimes used in the model
+   !!                        - ht_i(jpi,jpj,jpl)
+   !!                        - ht_s(jpi,jpj,jpl)
+   !!                        - t_i (jpi,jpj,nlay_i,jpl)
+   !!                        ...
+   !!                 VAGG : aggregate variables, averaged/summed over all
+   !!                        thickness categories
+   !!                        - vt_i(jpi,jpj)
+   !!                        - vt_s(jpi,jpj)
+   !!                        - at_i(jpi,jpj)
+   !!                        - et_s(jpi,jpj)  !total snow heat content
+   !!                        - et_i(jpi,jpj)  !total ice thermal content 
+   !!                        - smt_i(jpi,jpj) !mean ice salinity
+   !!                        - tm_i (jpi,jpj) !mean ice temperature
+   !!======================================================================
+   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
+   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   USE par_oce        ! ocean parameters
+   USE phycst         ! physical constants (ocean directory) 
+   USE sbc_oce , ONLY : sss_m
+   USE ice            ! ice variables
+   USE ice1D          ! ice variables (thermodynamics)
+   !
+   USE in_out_manager ! I/O manager
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_var_agg          
+   PUBLIC   ice_var_glo2eqv      
+   PUBLIC   ice_var_eqv2glo      
+   PUBLIC   ice_var_salprof      
+   PUBLIC   ice_var_bv           
+   PUBLIC   ice_var_salprof1d    
+   PUBLIC   ice_var_zapsmall
+   PUBLIC   ice_var_itd
+
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011)
+   !! $Id: icevar.F90 8422 2017-08-08 13:58:05Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_var_agg( kn )
+      !!------------------------------------------------------------------
+      !!                ***  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 )
+      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 )
+      ENDIF
+      ! END MP 2016
+
+      ! open water fraction
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )   
+         END DO
+      END DO
+
+      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
+         !
+      ENDIF
+      !
+   END SUBROUTINE ice_var_agg
+
+
+   SUBROUTINE ice_var_glo2eqv
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE ice_var_glo2eqv ***
+      !!
+      !! ** 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   !   -      -
+      !!------------------------------------------------------------------
+
+      !-------------------------------------------------------
+      ! 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
+      
+      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
+      ENDIF
+
+      CALL ice_var_salprof      ! salinity profile
+
+      !-------------------
+      ! Ice temperatures
+      !-------------------
+      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
+               END DO
+            END DO
+         END DO
+      END DO
+
+      !--------------------
+      ! Snow temperatures
+      !--------------------
+      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
+      END DO
+
+      ! integrated values 
+      vt_i (:,:) = SUM( v_i, dim=3 )
+      vt_s (:,:) = SUM( v_s, dim=3 )
+      at_i (:,:) = SUM( a_i, dim=3 )
+
+      ! MV MP 2016
+      ! probably should resum for melt ponds ???
+
+      !
+   END SUBROUTINE ice_var_glo2eqv
+
+
+   SUBROUTINE ice_var_eqv2glo
+      !!------------------------------------------------------------------
+      !!                ***  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(:,:,:)
+      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
+      !
+   END SUBROUTINE ice_var_eqv2glo
+
+
+   SUBROUTINE ice_var_salprof
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE ice_var_salprof ***
+      !!
+      !! ** Purpose :   computes salinity profile in function of bulk salinity     
+      !!
+      !! ** Method  : If bulk salinity greater than zsi1, 
+      !!              the profile is assumed to be constant (S_inf)
+      !!              If bulk salinity lower than zsi0,
+      !!              the profile is linear with 0 at the surface (S_zero)
+      !!              If it is between zsi0 and zsi1, it is a
+      !!              alpha-weighted linear combination of s_inf and s_zero
+      !!
+      !! ** References : Vancoppenolle et al., 2007
+      !!------------------------------------------------------------------
+      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), PARAMETER :: zsi0 = 3.5_wp
+      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
+      !!------------------------------------------------------------------
+
+      !---------------------------------------
+      ! Vertically constant, constant in time
+      !---------------------------------------
+      IF(  nn_icesal == 1  )  THEN
+         s_i (:,:,:,:) = rn_icesal
+         sm_i(:,:,:)   = rn_icesal
+      ENDIF
+
+      !-----------------------------------
+      ! Salinity profile, varying in time
+      !-----------------------------------
+      IF(  nn_icesal == 2  ) THEN
+         !
+         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
+         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 )
+               END DO
+            END DO
+         END DO
+
+         ! Computation of the profile
+         DO jl = 1, jpl
+            DO jk = 1, nlay_i
+               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 ) )
+                  END DO
+               END DO
+            END DO
+         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)
+         !
+         sm_i(:,:,:) = 2.30_wp
+         !
+         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
+      !
+   END SUBROUTINE ice_var_salprof
+
+
+   SUBROUTINE ice_var_bv
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE ice_var_bv ***
+      !!
+      !! ** Purpose :   computes mean brine volume (%) in sea ice
+      !!
+      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
+      !!
+      !! References : Vancoppenolle et al., JGR, 2007
+      !!------------------------------------------------------------------
+      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
+      !!------------------------------------------------------------------
+      !
+      bvm_i(:,:)   = 0._wp
+      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
+         END DO
+      END DO
+      !
+   END SUBROUTINE ice_var_bv
+
+
+   SUBROUTINE ice_var_salprof1d
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_var_salprof1d  ***
+      !!
+      !! ** Purpose :   1d computation of the sea ice salinity profile
+      !!                Works with 1d vectors and is used by thermodynamic modules
+      !!-------------------------------------------------------------------
+      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), 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
+            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 )
+         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 = (  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)
+         !
+         sm_i_1d(:) = 2.30_wp
+         !
+         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 ) ) ) )
+            DO ji = 1, nidx
+               s_i_1d(ji,jk) = zsal
+            END DO
+         END DO
+         !
+      ENDIF
+      !
+   END SUBROUTINE ice_var_salprof1d
+
+   SUBROUTINE ice_var_zapsmall
+      !!-------------------------------------------------------------------
+      !!                   ***  ROUTINE ice_var_zapsmall ***
+      !!
+      !! ** 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
+
+         !-----------------------------------------------------------------
+         ! Zap ice energy and use ocean heat to melt ice
+         !-----------------------------------------------------------------
+         DO jk = 1, nlay_i
+            DO jj = 1 , jpj
+               DO ji = 1 , jpi
+                  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
+                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
+                     &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
+                  zei              = e_i(ji,jj,jk,jl)
+                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
+                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
+                  ! update exchanges with ocean
+                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
+               END DO
+            END DO
+         END DO
+
+         DO jj = 1 , jpj
+            DO ji = 1 , jpi
+               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
+               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
+                  &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
+               zsal = smv_i(ji,jj,  jl)
+               zvi  = v_i  (ji,jj,  jl)
+               zvs  = v_s  (ji,jj,  jl)
+               zes  = e_s  (ji,jj,1,jl)
+               IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw )   zvp  = v_ip (ji,jj  ,jl)
+               !-----------------------------------------------------------------
+               ! Zap snow energy 
+               !-----------------------------------------------------------------
+               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
+               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
+
+               !-----------------------------------------------------------------
+               ! zap ice and snow volume, add water and salt to ocean
+               !-----------------------------------------------------------------
+               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
+               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch
+               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch
+               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch
+               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
+               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
+               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
+
+               ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * rswitch
+               ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * rswitch
+
+               ! MV MP 2016
+               IF ( ln_pnd ) THEN 
+                  a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch
+                  v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch
+                  IF ( ln_pnd_fw )   wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_ip(ji,jj,jl)  - zvp  ) * rhofw * r1_rdtice
+               ENDIF
+               ! END MV MP 2016
+
+               ! update exchanges with ocean
+               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
+               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
+               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
+               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
+            END DO
+         END DO
+      END DO 
+
+      ! to be sure that at_i is the sum of a_i(jl)
+      at_i (:,:) = 0._wp
+      DO jl = 1, jpl
+         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
+      END DO
+
+      ! open water = 1 if at_i=0
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
+            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
+         END DO
+      END DO
+
+      !
+   END SUBROUTINE ice_var_zapsmall
+
+   SUBROUTINE ice_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
+      !!------------------------------------------------------------------
+      !!                ***  ROUTINE ice_var_itd   ***
+      !!
+      !! ** Purpose :  converting 1-cat ice to multiple ice categories
+      !!
+      !!                  ice thickness distribution follows a gaussian law
+      !!               around the concentration of the most likely ice thickness
+      !!                           (similar as iceistate.F90)
+      !!
+      !! ** Method:   Iterative procedure
+      !!                
+      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
+      !!
+      !!               2) Check whether the distribution conserves area and volume, positivity and
+      !!                  category boundaries
+      !!              
+      !!               3) If not (input ice is too thin), the last category is empty and
+      !!                  the number of categories is reduced (jpl-1)
+      !!
+      !!               4) Iterate until ok (SUM(itest(:) = 4)
+      !!
+      !! ** Arguments : zhti: 1-cat ice thickness
+      !!                zhts: 1-cat snow depth
+      !!                zai : 1-cat ice concentration
+      !!
+      !! ** Output    : jpl-cat 
+      !!
+      !!  (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  
+      REAL(wp) :: zarg, zV, zconv, zdh, zdv
+      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
+      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)
+      zht_i(1:ijpij,1:jpl) = 0._wp
+      zht_s(1:ijpij,1:jpl) = 0._wp
+      za_i (1:ijpij,1:jpl) = 0._wp
+
+      ! ----------------------------------------
+      ! distribution over the jpl ice categories
+      ! ----------------------------------------
+      DO ji = 1, ijpij
+         
+         IF( zhti(ji) > 0._wp ) THEN
+
+            ! find which category (jl0) the input ice thickness falls into
+            jl0 = jpl
+            DO jl = 1, jpl
+               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
+                  jl0 = jl
+                  CYCLE
+               ENDIF
+            END DO
+
+            ! initialisation of tests
+            itest(:)  = 0
+         
+            i_fill = jpl + 1                                             !====================================
+            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
+               ! iteration                                               !====================================
+               i_fill = i_fill - 1
+               
+               ! initialisation of ice variables for each try
+               zht_i(ji,1:jpl) = 0._wp
+               za_i (ji,1:jpl) = 0._wp
+               itest(:)        = 0      
+               
+               ! *** case very thin ice: fill only category 1
+               IF ( i_fill == 1 ) THEN
+                  zht_i(ji,1) = zhti(ji)
+                  za_i (ji,1) = zai (ji)
+                  
+               ! *** case ice is thicker: fill categories >1
+               ELSE
+
+                  ! Fill ice thicknesses in the (i_fill-1) cat by hmean 
+                  DO jl = 1, i_fill - 1
+                     zht_i(ji,jl) = hi_mean(jl)
+                  END DO
+                  
+                  ! Concentrations in the (i_fill-1) categories 
+                  za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
+                  DO jl = 1, i_fill - 1
+                     IF ( jl /= jl0 ) THEN
+                        zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
+                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
+                     ENDIF
+                  END DO
+                  
+                  ! Concentration in the last (i_fill) category
+                  za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
+                  
+                  ! Ice thickness in the last (i_fill) category
+                  zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
+                  zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
+                  
+                  ! clem: correction if concentration of upper cat is greater than lower cat
+                  !       (it should be a gaussian around jl0 but sometimes it is not)
+                  IF ( jl0 /= jpl ) THEN
+                     DO jl = jpl, jl0+1, -1
+                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
+                           zdv = zht_i(ji,jl) * za_i(ji,jl)
+                           zht_i(ji,jl    ) = 0._wp
+                           za_i (ji,jl    ) = 0._wp
+                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
+                        END IF
+                     ENDDO
+                  ENDIF
+               
+               ENDIF ! case ice is thick or thin
+            
+               !---------------------
+               ! Compatibility tests
+               !--------------------- 
+               ! Test 1: area conservation
+               zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
+               IF ( zconv < epsi06 ) itest(1) = 1
+            
+               ! Test 2: volume conservation
+               zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
+               IF ( zconv < epsi06 ) itest(2) = 1
+               
+               ! Test 3: thickness of the last category is in-bounds ?
+               IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
+               
+               ! Test 4: positivity of ice concentrations
+               itest(4) = 1
+               DO jl = 1, i_fill
+                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
+               END DO
+               !                                         !============================
+            END DO                                       ! end iteration on categories
+               !                                         !============================
+         ENDIF ! if zhti > 0
+      END DO ! i loop
+
+      ! ------------------------------------------------
+      ! Adding Snow in each category where za_i is not 0
+      ! ------------------------------------------------ 
+      DO jl = 1, jpl
+         DO ji = 1, ijpij
+            IF( za_i(ji,jl) > 0._wp ) THEN
+               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
+               ! In case snow load is in excess that would lead to transformation from snow to ice
+               ! Then, transfer the snow excess into the ice (different from icethd_dh)
+               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
+               ! recompute ht_i, ht_s avoiding out of bounds values
+               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
+               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
+            ENDIF
+         ENDDO
+      ENDDO
+      !
+    END SUBROUTINE ice_var_itd
+
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option         Dummy module          NO  LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE ice_var_itd
+   END SUBROUTINE ice_var_itd
+#endif
+
+   !!======================================================================
+END MODULE icevar
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90	(revision 8424)
@@ -18,5 +18,5 @@
    USE sbc_ice         ! Surface boundary condition: ice fields
    USE ice
-   USE limvar
+   USE icevar
    !
    USE in_out_manager
@@ -78,5 +78,5 @@
 
       ! brine volume
-      CALL lim_var_bv 
+      CALL ice_var_bv 
 
       ! tresholds for outputs
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 8422)
+++ 	(revision )
@@ -1,798 +1,0 @@
-MODULE limvar
-   !!======================================================================
-   !!                       ***  MODULE limvar ***
-   !!                 Different sets of ice model variables 
-   !!                   how to switch from one to another
-   !!
-   !!                 There are three sets of variables
-   !!                 VGLO : global variables of the model
-   !!                        - v_i (jpi,jpj,jpl)
-   !!                        - v_s (jpi,jpj,jpl)
-   !!                        - a_i (jpi,jpj,jpl)
-   !!                        - t_s (jpi,jpj,jpl)
-   !!                        - e_i (jpi,jpj,nlay_i,jpl)
-   !!                        - smv_i(jpi,jpj,jpl)
-   !!                        - oa_i (jpi,jpj,jpl)
-   !!                 VEQV : equivalent variables sometimes used in the model
-   !!                        - ht_i(jpi,jpj,jpl)
-   !!                        - ht_s(jpi,jpj,jpl)
-   !!                        - t_i (jpi,jpj,nlay_i,jpl)
-   !!                        ...
-   !!                 VAGG : aggregate variables, averaged/summed over all
-   !!                        thickness categories
-   !!                        - vt_i(jpi,jpj)
-   !!                        - vt_s(jpi,jpj)
-   !!                        - at_i(jpi,jpj)
-   !!                        - et_s(jpi,jpj)  !total snow heat content
-   !!                        - et_i(jpi,jpj)  !total ice thermal content 
-   !!                        - smt_i(jpi,jpj) !mean ice salinity
-   !!                        - tm_i (jpi,jpj) !mean ice temperature
-   !!======================================================================
-   !! History :   -   ! 2006-01 (M. Vancoppenolle) Original code
-   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                      LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   USE par_oce        ! ocean parameters
-   USE phycst         ! physical constants (ocean directory) 
-   USE sbc_oce , ONLY : sss_m
-   USE ice            ! ice variables
-   USE ice1D          ! ice variables (thermodynamics)
-   !
-   USE in_out_manager ! I/O manager
-   USE lib_mpp        ! MPP library
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   lim_var_agg          
-   PUBLIC   lim_var_glo2eqv      
-   PUBLIC   lim_var_eqv2glo      
-   PUBLIC   lim_var_salprof      
-   PUBLIC   lim_var_bv           
-   PUBLIC   lim_var_salprof1d    
-   PUBLIC   lim_var_zapsmall
-   PUBLIC   lim_var_itd
-
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 3.5 , UCL - NEMO Consortium (2011)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE lim_var_agg( kn )
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_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 )
-      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 )
-      ENDIF
-      ! END MP 2016
-
-      ! open water fraction
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            ato_i(ji,jj) = MAX( 1._wp - at_i(ji,jj), 0._wp )   
-         END DO
-      END DO
-
-      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
-         !
-      ENDIF
-      !
-   END SUBROUTINE lim_var_agg
-
-
-   SUBROUTINE lim_var_glo2eqv
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_var_glo2eqv ***
-      !!
-      !! ** 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   !   -      -
-      !!------------------------------------------------------------------
-
-      !-------------------------------------------------------
-      ! 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
-      
-      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
-      ENDIF
-
-      CALL lim_var_salprof      ! salinity profile
-
-      !-------------------
-      ! Ice temperatures
-      !-------------------
-      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
-               END DO
-            END DO
-         END DO
-      END DO
-
-      !--------------------
-      ! Snow temperatures
-      !--------------------
-      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
-      END DO
-
-      ! integrated values 
-      vt_i (:,:) = SUM( v_i, dim=3 )
-      vt_s (:,:) = SUM( v_s, dim=3 )
-      at_i (:,:) = SUM( a_i, dim=3 )
-
-      ! MV MP 2016
-      ! probably should resum for melt ponds ???
-
-      !
-   END SUBROUTINE lim_var_glo2eqv
-
-
-   SUBROUTINE lim_var_eqv2glo
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_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(:,:,:)
-      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
-      !
-   END SUBROUTINE lim_var_eqv2glo
-
-
-   SUBROUTINE lim_var_salprof
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_var_salprof ***
-      !!
-      !! ** Purpose :   computes salinity profile in function of bulk salinity     
-      !!
-      !! ** Method  : If bulk salinity greater than zsi1, 
-      !!              the profile is assumed to be constant (S_inf)
-      !!              If bulk salinity lower than zsi0,
-      !!              the profile is linear with 0 at the surface (S_zero)
-      !!              If it is between zsi0 and zsi1, it is a
-      !!              alpha-weighted linear combination of s_inf and s_zero
-      !!
-      !! ** References : Vancoppenolle et al., 2007
-      !!------------------------------------------------------------------
-      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), PARAMETER :: zsi0 = 3.5_wp
-      REAL(wp), PARAMETER :: zsi1 = 4.5_wp
-      !!------------------------------------------------------------------
-
-      !---------------------------------------
-      ! Vertically constant, constant in time
-      !---------------------------------------
-      IF(  nn_icesal == 1  )  THEN
-         s_i (:,:,:,:) = rn_icesal
-         sm_i(:,:,:)   = rn_icesal
-      ENDIF
-
-      !-----------------------------------
-      ! Salinity profile, varying in time
-      !-----------------------------------
-      IF(  nn_icesal == 2  ) THEN
-         !
-         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
-         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 )
-               END DO
-            END DO
-         END DO
-
-         ! Computation of the profile
-         DO jl = 1, jpl
-            DO jk = 1, nlay_i
-               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 ) )
-                  END DO
-               END DO
-            END DO
-         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)
-         !
-         sm_i(:,:,:) = 2.30_wp
-         !
-         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
-      !
-   END SUBROUTINE lim_var_salprof
-
-
-   SUBROUTINE lim_var_bv
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_var_bv ***
-      !!
-      !! ** Purpose :   computes mean brine volume (%) in sea ice
-      !!
-      !! ** Method  : e = - 0.054 * S (ppt) / T (C)
-      !!
-      !! References : Vancoppenolle et al., JGR, 2007
-      !!------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices
-      !!------------------------------------------------------------------
-      !
-      bvm_i(:,:)   = 0._wp
-      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
-         END DO
-      END DO
-      !
-   END SUBROUTINE lim_var_bv
-
-
-   SUBROUTINE lim_var_salprof1d
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE lim_var_salprof1d  ***
-      !!
-      !! ** Purpose :   1d computation of the sea ice salinity profile
-      !!                Works with 1d vectors and is used by thermodynamic modules
-      !!-------------------------------------------------------------------
-      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), 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
-            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 )
-         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 = (  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)
-         !
-         sm_i_1d(:) = 2.30_wp
-         !
-         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 ) ) ) )
-            DO ji = 1, nidx
-               s_i_1d(ji,jk) = zsal
-            END DO
-         END DO
-         !
-      ENDIF
-      !
-   END SUBROUTINE lim_var_salprof1d
-
-   SUBROUTINE lim_var_zapsmall
-      !!-------------------------------------------------------------------
-      !!                   ***  ROUTINE lim_var_zapsmall ***
-      !!
-      !! ** 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
-
-         !-----------------------------------------------------------------
-         ! Zap ice energy and use ocean heat to melt ice
-         !-----------------------------------------------------------------
-         DO jk = 1, nlay_i
-            DO jj = 1 , jpj
-               DO ji = 1 , jpi
-                  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
-                  rswitch          = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
-                     &                                       / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
-                  zei              = e_i(ji,jj,jk,jl)
-                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * rswitch
-                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * rswitch + rt0 * ( 1._wp - rswitch )
-                  ! update exchanges with ocean
-                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * r1_rdtice ! W.m-2 <0
-               END DO
-            END DO
-         END DO
-
-         DO jj = 1 , jpj
-            DO ji = 1 , jpi
-               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
-               rswitch = MAX( 0._wp , SIGN( 1._wp, v_i(ji,jj,jl) * rswitch  &
-                  &                              / MAX( a_i(ji,jj,jl), epsi10 ) - epsi10 ) ) * rswitch
-               zsal = smv_i(ji,jj,  jl)
-               zvi  = v_i  (ji,jj,  jl)
-               zvs  = v_s  (ji,jj,  jl)
-               zes  = e_s  (ji,jj,1,jl)
-               IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw )   zvp  = v_ip (ji,jj  ,jl)
-               !-----------------------------------------------------------------
-               ! Zap snow energy 
-               !-----------------------------------------------------------------
-               t_s(ji,jj,1,jl) = t_s(ji,jj,1,jl) * rswitch + rt0 * ( 1._wp - rswitch )
-               e_s(ji,jj,1,jl) = e_s(ji,jj,1,jl) * rswitch
-
-               !-----------------------------------------------------------------
-               ! zap ice and snow volume, add water and salt to ocean
-               !-----------------------------------------------------------------
-               ato_i(ji,jj)    = a_i  (ji,jj,jl) * ( 1._wp - rswitch ) + ato_i(ji,jj)
-               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * rswitch
-               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * rswitch
-               v_s  (ji,jj,jl) = v_s  (ji,jj,jl) * rswitch
-               t_su (ji,jj,jl) = t_su (ji,jj,jl) * rswitch + t_bo(ji,jj) * ( 1._wp - rswitch )
-               oa_i (ji,jj,jl) = oa_i (ji,jj,jl) * rswitch
-               smv_i(ji,jj,jl) = smv_i(ji,jj,jl) * rswitch
-
-               ht_i (ji,jj,jl) = ht_i (ji,jj,jl) * rswitch
-               ht_s (ji,jj,jl) = ht_s (ji,jj,jl) * rswitch
-
-               ! MV MP 2016
-               IF ( ln_pnd ) THEN 
-                  a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * rswitch
-                  v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * rswitch
-                  IF ( ln_pnd_fw )   wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_ip(ji,jj,jl)  - zvp  ) * rhofw * r1_rdtice
-               ENDIF
-               ! END MV MP 2016
-
-               ! update exchanges with ocean
-               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice
-               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice
-               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice
-               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * r1_rdtice ! W.m-2 <0
-            END DO
-         END DO
-      END DO 
-
-      ! to be sure that at_i is the sum of a_i(jl)
-      at_i (:,:) = 0._wp
-      DO jl = 1, jpl
-         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-      END DO
-
-      ! open water = 1 if at_i=0
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            rswitch      = MAX( 0._wp , SIGN( 1._wp, - at_i(ji,jj) ) )
-            ato_i(ji,jj) = rswitch + (1._wp - rswitch ) * ato_i(ji,jj)
-         END DO
-      END DO
-
-      !
-   END SUBROUTINE lim_var_zapsmall
-
-   SUBROUTINE lim_var_itd( zhti, zhts, zai, zht_i, zht_s, za_i )
-      !!------------------------------------------------------------------
-      !!                ***  ROUTINE lim_var_itd   ***
-      !!
-      !! ** Purpose :  converting 1-cat ice to multiple ice categories
-      !!
-      !!                  ice thickness distribution follows a gaussian law
-      !!               around the concentration of the most likely ice thickness
-      !!                           (similar as iceistate.F90)
-      !!
-      !! ** Method:   Iterative procedure
-      !!                
-      !!               1) Try to fill the jpl ice categories (bounds hi_max(0:jpl)) with a gaussian
-      !!
-      !!               2) Check whether the distribution conserves area and volume, positivity and
-      !!                  category boundaries
-      !!              
-      !!               3) If not (input ice is too thin), the last category is empty and
-      !!                  the number of categories is reduced (jpl-1)
-      !!
-      !!               4) Iterate until ok (SUM(itest(:) = 4)
-      !!
-      !! ** Arguments : zhti: 1-cat ice thickness
-      !!                zhts: 1-cat snow depth
-      !!                zai : 1-cat ice concentration
-      !!
-      !! ** Output    : jpl-cat 
-      !!
-      !!  (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  
-      REAL(wp) :: zarg, zV, zconv, zdh, zdv
-      REAL(wp), DIMENSION(:),   INTENT(in)    ::   zhti, zhts, zai    ! input ice/snow variables
-      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)
-      zht_i(1:ijpij,1:jpl) = 0._wp
-      zht_s(1:ijpij,1:jpl) = 0._wp
-      za_i (1:ijpij,1:jpl) = 0._wp
-
-      ! ----------------------------------------
-      ! distribution over the jpl ice categories
-      ! ----------------------------------------
-      DO ji = 1, ijpij
-         
-         IF( zhti(ji) > 0._wp ) THEN
-
-            ! find which category (jl0) the input ice thickness falls into
-            jl0 = jpl
-            DO jl = 1, jpl
-               IF ( ( zhti(ji) >= hi_max(jl-1) ) .AND. ( zhti(ji) < hi_max(jl) ) ) THEN
-                  jl0 = jl
-                  CYCLE
-               ENDIF
-            END DO
-
-            ! initialisation of tests
-            itest(:)  = 0
-         
-            i_fill = jpl + 1                                             !====================================
-            DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) )  ! iterative loop on i_fill categories
-               ! iteration                                               !====================================
-               i_fill = i_fill - 1
-               
-               ! initialisation of ice variables for each try
-               zht_i(ji,1:jpl) = 0._wp
-               za_i (ji,1:jpl) = 0._wp
-               itest(:)        = 0      
-               
-               ! *** case very thin ice: fill only category 1
-               IF ( i_fill == 1 ) THEN
-                  zht_i(ji,1) = zhti(ji)
-                  za_i (ji,1) = zai (ji)
-                  
-               ! *** case ice is thicker: fill categories >1
-               ELSE
-
-                  ! Fill ice thicknesses in the (i_fill-1) cat by hmean 
-                  DO jl = 1, i_fill - 1
-                     zht_i(ji,jl) = hi_mean(jl)
-                  END DO
-                  
-                  ! Concentrations in the (i_fill-1) categories 
-                  za_i(ji,jl0) = zai(ji) / SQRT(REAL(jpl))
-                  DO jl = 1, i_fill - 1
-                     IF ( jl /= jl0 ) THEN
-                        zarg        = ( zht_i(ji,jl) - zhti(ji) ) / ( zhti(ji) * 0.5_wp )
-                        za_i(ji,jl) =   za_i (ji,jl0) * EXP(-zarg**2)
-                     ENDIF
-                  END DO
-                  
-                  ! Concentration in the last (i_fill) category
-                  za_i(ji,i_fill) = zai(ji) - SUM( za_i(ji,1:i_fill-1) )
-                  
-                  ! Ice thickness in the last (i_fill) category
-                  zV = SUM( za_i(ji,1:i_fill-1) * zht_i(ji,1:i_fill-1) )
-                  zht_i(ji,i_fill) = ( zhti(ji) * zai(ji) - zV ) / MAX( za_i(ji,i_fill), epsi10 ) 
-                  
-                  ! clem: correction if concentration of upper cat is greater than lower cat
-                  !       (it should be a gaussian around jl0 but sometimes it is not)
-                  IF ( jl0 /= jpl ) THEN
-                     DO jl = jpl, jl0+1, -1
-                        IF ( za_i(ji,jl) > za_i(ji,jl-1) ) THEN
-                           zdv = zht_i(ji,jl) * za_i(ji,jl)
-                           zht_i(ji,jl    ) = 0._wp
-                           za_i (ji,jl    ) = 0._wp
-                           za_i (ji,1:jl-1) = za_i(ji,1:jl-1) + zdv / MAX( REAL(jl-1) * zhti(ji), epsi10 )
-                        END IF
-                     ENDDO
-                  ENDIF
-               
-               ENDIF ! case ice is thick or thin
-            
-               !---------------------
-               ! Compatibility tests
-               !--------------------- 
-               ! Test 1: area conservation
-               zconv = ABS( zai(ji) - SUM( za_i(ji,1:jpl) ) )
-               IF ( zconv < epsi06 ) itest(1) = 1
-            
-               ! Test 2: volume conservation
-               zconv = ABS( zhti(ji)*zai(ji) - SUM( za_i(ji,1:jpl)*zht_i(ji,1:jpl) ) )
-               IF ( zconv < epsi06 ) itest(2) = 1
-               
-               ! Test 3: thickness of the last category is in-bounds ?
-               IF ( zht_i(ji,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1
-               
-               ! Test 4: positivity of ice concentrations
-               itest(4) = 1
-               DO jl = 1, i_fill
-                  IF ( za_i(ji,jl) < 0._wp ) itest(4) = 0
-               END DO
-               !                                         !============================
-            END DO                                       ! end iteration on categories
-               !                                         !============================
-         ENDIF ! if zhti > 0
-      END DO ! i loop
-
-      ! ------------------------------------------------
-      ! Adding Snow in each category where za_i is not 0
-      ! ------------------------------------------------ 
-      DO jl = 1, jpl
-         DO ji = 1, ijpij
-            IF( za_i(ji,jl) > 0._wp ) THEN
-               zht_s(ji,jl) = zht_i(ji,jl) * ( zhts(ji) / zhti(ji) )
-               ! In case snow load is in excess that would lead to transformation from snow to ice
-               ! Then, transfer the snow excess into the ice (different from icethd_dh)
-               zdh = MAX( 0._wp, ( rhosn * zht_s(ji,jl) + ( rhoic - rau0 ) * zht_i(ji,jl) ) * r1_rau0 ) 
-               ! recompute ht_i, ht_s avoiding out of bounds values
-               zht_i(ji,jl) = MIN( hi_max(jl), zht_i(ji,jl) + zdh )
-               zht_s(ji,jl) = MAX( 0._wp, zht_s(ji,jl) - zdh * rhoic * r1_rhosn )
-            ENDIF
-         ENDDO
-      ENDDO
-      !
-    END SUBROUTINE lim_var_itd
-
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option         Dummy module          NO  LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_var_agg          ! Empty routines
-   END SUBROUTINE lim_var_agg
-   SUBROUTINE lim_var_glo2eqv      ! Empty routines
-   END SUBROUTINE lim_var_glo2eqv
-   SUBROUTINE lim_var_eqv2glo      ! Empty routines
-   END SUBROUTINE lim_var_eqv2glo
-   SUBROUTINE lim_var_salprof      ! Empty routines
-   END SUBROUTINE lim_var_salprof
-   SUBROUTINE lim_var_bv           ! Emtpy routines
-   END SUBROUTINE lim_var_bv
-   SUBROUTINE lim_var_salprof1d    ! Emtpy routines
-   END SUBROUTINE lim_var_salprof1d
-   SUBROUTINE lim_var_zapsmall
-   END SUBROUTINE lim_var_zapsmall
-   SUBROUTINE lim_var_itd
-   END SUBROUTINE lim_var_itd
-#endif
-
-   !!======================================================================
-END MODULE limvar
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdydta.F90	(revision 8424)
@@ -27,5 +27,5 @@
 #if defined key_lim3
    USE ice
-   USE limvar          ! redistribute ice input into categories
+   USE icevar          ! redistribute ice input into categories
 #endif
    USE sbcapr
@@ -344,5 +344,5 @@
 #if defined key_lim3
                IF( .NOT. ll_bdylim3 .AND. cn_ice_lim(ib_bdy) /= 'none' .AND. nn_ice_lim_dta(ib_bdy) == 1 ) THEN ! bdy ice input (case input is 1cat)
-     	          CALL lim_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
+     	          CALL ice_var_itd ( bf(jfld_hti)%fnow(:,1,1), bf(jfld_hts)%fnow(:,1,1), bf(jfld_ai)%fnow(:,1,1), &
                                   & dta_bdy(ib_bdy)%ht_i,     dta_bdy(ib_bdy)%ht_s,     dta_bdy(ib_bdy)%a_i     )
                ENDIF
Index: branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90
===================================================================
--- branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90	(revision 8422)
+++ branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90	(revision 8424)
@@ -20,5 +20,5 @@
    USE oce             ! ocean dynamics and tracers variables
    USE ice             ! LIM_3 ice variables
-   USE limvar
+   USE icevar
    USE icectl
    USE par_oce         ! ocean parameters
@@ -56,5 +56,5 @@
       !!----------------------------------------------------------------------
       !
-      CALL lim_var_glo2eqv
+      CALL ice_var_glo2eqv
       !
       DO ib_bdy=1, nb_bdy
@@ -71,6 +71,6 @@
       END DO
       !
-                        CALL lim_var_zapsmall
-                        CALL lim_var_agg(1)
+                        CALL ice_var_zapsmall
+                        CALL ice_var_agg(1)
       IF( ln_limctl )   CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )
       !
