Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/cpp_ORCA2_LIM3.fcm
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/cpp_ORCA2_LIM3.fcm	(revision 4099)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/cpp_ORCA2_LIM3.fcm	(revision 4099)
@@ -0,0 +1,1 @@
+ bld::tool::fppkeys key_trabbl key_orca_r2 key_lim3 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_zdftke key_zdfddm key_zdftmx key_iomput key_mpp_mpi
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/cfg.txt
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/cfg.txt	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/cfg.txt	(revision 4099)
@@ -1,10 +1,10 @@
-AMM12 OPA_SRC
 GYRE OPA_SRC
 GYRE_BFM OPA_SRC TOP_SRC
 GYRE_PISCES OPA_SRC TOP_SRC
-ORCA2_LIM3 OPA_SRC LIM_SRC_3
 ORCA2_SAS_LIM OPA_SRC SAS_SRC LIM_SRC_2 NST_SRC
 ORCA2_LIM_CFC_C14b OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC
-ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC
 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC
 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC
+AMM12 OPA_SRC
+ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC
+ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4099)
@@ -395,5 +395,5 @@
    CHARACTER(len=32)     , PUBLIC ::   cn_icerst_out = "restart_ice"      !: suffix of ice restart name (output)
    LOGICAL               , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F)
-   LOGICAL               , PUBLIC ::   ln_nicep      = .TRUE.             !: flag for sea-ice points output (T) or not (F)
+   LOGICAL               , PUBLIC ::   ln_nicep      = .FALSE.             !: flag for sea-ice points output (T) or not (F)
    REAL(wp)              , PUBLIC ::   cai           = 1.40e-3            !: atmospheric drag over sea ice
    REAL(wp)              , PUBLIC ::   cao           = 1.00e-3            !: atmospheric drag over ocean
@@ -404,5 +404,6 @@
    !!--------------------------------------------------------------------------
    !! Check if everything down here is necessary
-   LOGICAL , PUBLIC                                      ::   ln_limdiahsb  = .TRUE. !: flag for ice diag (T) or not (F)
+   LOGICAL , PUBLIC                                      ::   ln_limdiahsb  = .FALSE. !: flag for ice diag (T) or not (F)
+   LOGICAL , PUBLIC                                      ::   ln_limdiaout  = .FALSE. !: flag for ice diag (T) or not (F)
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   v_newice   !: volume of ice formed in the leads
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dv_dt_thd  !: thermodynamic growth rates 
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4099)
@@ -126,5 +126,5 @@
       !! ** input   :   Namelist namicerun
       !!-------------------------------------------------------------------
-      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb
+      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb, ln_limdiaout
       !!-------------------------------------------------------------------
       !                    
@@ -147,4 +147,5 @@
          WRITE(numout,*) '   Several ice points in the ice or not in ocean.output    = ', ln_nicep
          WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb
+         WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout
       ENDIF
       !
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_1D.F90	(revision 4099)
@@ -154,14 +154,14 @@
                                                            !============================
          ! Check if tests have passed (i.e. volume conservation...)
-         IF ( ztests /= 4 ) THEN
-            WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
-            WRITE(numout,*) ' !! ALERT categories distribution !!'
-            WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
-            WRITE(numout,*) ' *** ztests is not equal to 4 '
-            WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4
-            WRITE(numout,*) 'i_fill=',i_fill
-            WRITE(numout,*) 'zai(ji)=',zai(ji)
-            WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:)
-         ENDIF
+         !IF ( ztests /= 4 ) THEN
+         !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
+         !   WRITE(numout,*) ' !! ALERT categories distribution !!'
+         !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
+         !   WRITE(numout,*) ' *** ztests is not equal to 4 '
+         !   WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4
+         !   WRITE(numout,*) 'i_fill=',i_fill
+         !   WRITE(numout,*) 'zai(ji)=',zai(ji)
+         !   WRITE(numout,*) 'za_i(ji,jpl)=',za_i(ji,:)
+         !ENDIF
             
       END DO ! i loop
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_2D.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_2D.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcat_2D.F90	(revision 4099)
@@ -155,11 +155,11 @@
                                                               !============================
             ! Check if tests have passed (i.e. volume conservation...)
-            IF ( ztests .NE. 4 ) THEN
-               WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
-               WRITE(numout,*) ' !! ALERT categories distribution !!'
-               WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
-               WRITE(numout,*) ' *** ztests is not equal to 4 '
-               WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4
-            ENDIF
+            !IF ( ztests .NE. 4 ) THEN
+            !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
+            !   WRITE(numout,*) ' !! ALERT categories distribution !!'
+            !   WRITE(numout,*) ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
+            !   WRITE(numout,*) ' *** ztests is not equal to 4 '
+            !   WRITE(numout,*) ' *** ztest (1:4) = ', ztest_1, ztest_2, ztest_3, ztest_4
+            !ENDIF
 
          END DO ! i loop
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90	(revision 4098)
+++ 	(revision )
@@ -1,592 +1,0 @@
-MODULE limdia
-   !!======================================================================
-   !!                       ***  MODULE limdia   ***
-   !!  LIM-3 sea ice model :   diagnostics of ice model 
-   !!======================================================================
-   !! History :  3.2  ! 2007-01  (M. Vancoppenolle)  Code adapted from LIM-2
-   !!             -   ! 2008-03  (M. Vancoppenolle)  add lim_dia_init
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                       LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   lim_dia       : computation and output of the time evolution of keys variables
-   !!   lim_dia_init  : initialization and namelist read
-   !!----------------------------------------------------------------------
-   USE ice            ! LIM-3: sea-ice variable
-   USE par_ice        ! LIM-3: ice parameters
-   USE dom_ice        ! LIM-3: sea-ice domain
-   USE dom_oce        ! ocean domain
-   USE sbc_oce        ! surface boundary condition: ocean fields
-   USE daymod         ! model calendar
-   USE phycst         ! physical constant
-   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_dia       ! called by ice_step
-
-   INTEGER, PUBLIC  ::   ntmoy   = 1   !: instantaneous values of ice evolution or averaging ntmoy
-   INTEGER, PUBLIC  ::   ninfo   = 1   !: frequency of ouputs on file ice_evolu in case of averaging
-
-   !                                              !!! Parameters for outputs to files "evolu"
-   INTEGER, PARAMETER ::   jpinfmx = 100           ! maximum number of key variables
-   INTEGER, PARAMETER ::   jpchinf = 5             ! ???
-   INTEGER, PARAMETER ::   jpchsep = jpchinf + 2   ! ???
-
-   INTEGER  ::   nfrinf  = 4         ! number of variables written in one line 
-   INTEGER  ::   nferme              ! last time step at which the var. are written on file
-   INTEGER  ::   nvinfo              ! number of total variables 
-   INTEGER  ::   nbvt                ! number of time variables
-   INTEGER  ::   naveg               ! number of step for accumulation before averaging
-   REAL(wp) ::   epsi06 = 1.e-6_wp   ! small number
-
-   CHARACTER(len= 8) ::   fmtinf = '1PE13.5 '   ! format of the output values  
-   CHARACTER(len=30) ::   fmtw                  ! formats
-   CHARACTER(len=30) ::   fmtr                  ! ???
-   CHARACTER(len=30) ::   fmtitr                ! ???
-
-   CHARACTER(len=jpchsep), DIMENSION(jpinfmx) ::   titvar   ! title of key variables
-
-   REAL(wp), DIMENSION(jpinfmx) ::   vinfom     ! temporary working space
-   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) ::   aire       ! masked grid cell area
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE lim_dia
-      !!--------------------------------------------------------------------
-      !!                  ***  ROUTINE lim_dia  ***
-      !!   
-      !! ** Purpose :   Computation and outputs on file ice.evolu 
-      !!              the temporal evolution of some key variables
-      !!-------------------------------------------------------------------
-      INTEGER  ::   jv, ji, jj, jl       ! dummy loop indices
-      INTEGER  ::   ii0, ii1, ij0, ij1   ! temporary integer
-      REAL(wp) ::   zshift_date          ! date from the minimum ice extent
-      REAL(wp) ::   zday, zday_min       ! current day, day of minimum extent
-      REAL(wp) ::   zafy, zamy           ! temporary area of fy and my ice
-      REAL(wp) ::   zindb
-      REAL(wp), DIMENSION(jpinfmx) ::   vinfor   ! 1D workspace 
-      !!-------------------------------------------------------------------
-
-      ! 0) date from the minimum of ice extent
-      !---------------------------------------
-      zday_min = 273._wp        ! zday_min = date of minimum extent, here September 30th
-      zday = REAL(numit-nit000,wp) * rdt_ice / ( 86400._wp * REAL(nn_fsbc,wp) )
-      !
-      IF( zday > zday_min ) THEN   ;   zshift_date  =  zday - zday_min
-      ELSE                         ;   zshift_date  =  zday - (365.0 - zday_min)
-      ENDIF
-
-      IF( numit == nstart )   CALL lim_dia_init   ! initialisation of ice_evolu file      
-
-      vinfor(1) = REAL(numit)       ! time diagnostics 
-      vinfor(2) = nyear
-
-      DO jv = nbvt + 1, nvinfo      ! put everything to zero
-         vinfor(jv) = 0._wp
-      END DO
-
-      !!-------------------------------------------------------------------
-      !! 1) Northern hemisphere
-      !!-------------------------------------------------------------------
-      !! 1.1) Diagnostics independent on age
-      !!------------------------------------
-      DO jj = njeq, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            IF( tms(ji,jj) == 1 ) THEN
-               vinfor(3)  = vinfor(3)  + at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !ice area
-               IF ( at_i(ji,jj) > 0.15 )   vinfor(5) = vinfor(5) + aire(ji,jj) * 1.e-12_wp !ice extent
-               vinfor(7)  = vinfor(7)  + vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !ice volume
-               vinfor(9)  = vinfor(9)  + vt_s(ji,jj)*aire(ji,jj) * 1.e-12_wp !snow volume
-               vinfor(15) = vinfor(15) + ot_i(ji,jj) *vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !mean age
-               vinfor(29) = vinfor(29) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !mean salinity
-               ! the computation of this diagnostic is not reliable
-               vinfor(31) = vinfor(31) + vt_i(ji,jj) * (  u_ice(ji,jj)*u_ice(ji,jj)   & 
-                  &                                     + v_ice(ji,jj)*v_ice(ji,jj) ) * aire(ji,jj) * 1.e-12 
-               vinfor(53) = vinfor(53) + sfx (ji,jj)*aire(ji,jj) * 1.e-12_wp !salt flux
-               vinfor(55) = vinfor(55) + sfx_bri(ji,jj)*aire(ji,jj) * 1.e-12_wp !brine drainage flux
-               vinfor(57) = vinfor(57) + sfx_thd(ji,jj)*aire(ji,jj) * 1.e-12_wp !equivalent salt flux
-               vinfor(59) = vinfor(59) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp  !SST
-               vinfor(61) = vinfor(61) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp  !SSS
-               vinfor(65) = vinfor(65) + et_s(ji,jj)/1.0e9*aire(ji,jj) * 1.e-12_wp  ! snow temperature
-               vinfor(67) = vinfor(67) + et_i(ji,jj)/1.0e9*aire(ji,jj) * 1.e-12_wp       ! ice heat content
-               vinfor(69) = vinfor(69) + v_i(ji,jj,1)*aire(ji,jj) * 1.e-12_wp !ice volume
-               vinfor(71) = vinfor(71) + v_i(ji,jj,2)*aire(ji,jj) * 1.e-12_wp !ice volume
-               vinfor(73) = vinfor(73) + v_i(ji,jj,3)*aire(ji,jj) * 1.e-12_wp !ice volume
-               vinfor(75) = vinfor(75) + v_i(ji,jj,4)*aire(ji,jj) * 1.e-12_wp !ice volume
-               vinfor(77) = vinfor(77) + v_i(ji,jj,5)*aire(ji,jj) * 1.e-12_wp !ice volume
-               vinfor(79) = 0.0
-               vinfor(81) = vinfor(81) + emp(ji,jj)*aire(ji,jj) * 1.e-12_wp ! mass flux
-            ENDIF
-         END DO
-      END DO
-
-      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
-         DO jj = njeq, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               IF( tms(ji,jj) == 1 ) THEN
-                  vinfor(11) = vinfor(11) + v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !undef def ice volume
-               ENDIF
-            END DO
-         END DO
-      END DO
-
-      vinfor(13) = 0._wp
-
-      vinfor(15) = vinfor(15) / MAX(vinfor(7),epsi06) ! these have to be divided by total ice volume to have the
-      vinfor(29) = vinfor(29) / MAX(vinfor(7),epsi06) ! right value
-      vinfor(31) = SQRT( vinfor(31) / MAX( vinfor(7) , epsi06 ) )
-      vinfor(67) = vinfor(67) / MAX(vinfor(7),epsi06)
-
-      vinfor(53) = vinfor(53) / MAX(vinfor(5),epsi06) ! these have to be divided by total ice extent to have the
-      vinfor(55) = vinfor(55) / MAX(vinfor(5),epsi06) ! right value 
-      vinfor(57) = vinfor(57) / MAX(vinfor(5),epsi06) ! 
-      vinfor(79) = vinfor(79) / MAX(vinfor(5),epsi06) !
-
-      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(3))) !
-      vinfor(59) = zindb*vinfor(59) / MAX(vinfor(3),epsi06) ! divide by ice area
-      vinfor(61) = zindb*vinfor(61) / MAX(vinfor(3),epsi06) !
-
-      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(9))) !
-      vinfor(65) = zindb*vinfor(65) / MAX(vinfor(9),epsi06) ! divide it by snow volume
-
-
-      DO jl = 1, jpl
-         DO jj = njeq, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               IF( tms(ji,jj) == 1 ) THEN
-                  vinfor(33) = vinfor(33) + d_v_i_trp(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !ice volume
-                  vinfor(35) = vinfor(35) + d_v_i_thd(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !ice volume
-               ENDIF
-            END DO
-         END DO
-      END DO
-
-      DO jj = njeq, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            IF( tms(ji,jj) == 1 ) THEN
-               vinfor(37) = vinfor(37) + diag_sni_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp !th growth rates
-               vinfor(39) = vinfor(39) + diag_lat_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp 
-               vinfor(41) = vinfor(41) + diag_bot_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp
-               vinfor(43) = vinfor(43) + diag_dyn_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp 
-               vinfor(45) = vinfor(45) + dv_dt_thd(ji,jj,5)*aire(ji,jj) * 1.e-12_wp
-               vinfor(47) = vinfor(47) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp * r1_rdtice   ! volume acc in OW
-            ENDIF
-         END DO
-      END DO
-
-      DO jl = 1, jpl
-         DO jj = njeq, jpjm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               IF( tms(ji,jj) == 1 ) THEN
-                  vinfor(63) = vinfor(63) + t_su(ji,jj,jl)*a_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp
-               ENDIF
-            END DO
-         END DO
-      END DO
-      vinfor(63) = vinfor(63) / MAX(vinfor(3),epsi06) ! these have to be divided by total ice area
-
-      !! 1.2) Diagnostics dependent on age
-      !!------------------------------------
-      DO jj = njeq, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            IF( tms(ji,jj) == 1 ) THEN
-               zafy = 0.0
-               zamy = 0.0
-               DO jl = 1, jpl
-                  IF ((o_i(ji,jj,jl) - zshift_date).LT.0.0) THEN
-                     vinfor(17) = vinfor(17) + a_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp ! FY ice area
-                     vinfor(25) = vinfor(25) + v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp ! FY ice volume
-                     vinfor(49) = vinfor(49) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !FY ice salinity
-                     zafy = zafy + a_i(ji,jj,jl)
-                  ENDIF
-                  IF ((o_i(ji,jj,jl) - zshift_date).GT.0.0) THEN
-                     vinfor(19) = vinfor(19) + a_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp    ! MY ice area
-                     vinfor(27) = vinfor(27) + v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp ! MY ice volume
-                     vinfor(51) = vinfor(51) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !MY ice salinity
-                     zamy = zamy + a_i(ji,jj,jl)
-                  ENDIF
-               END DO
-               IF ((at_i(ji,jj).GT.0.15).AND.(zafy.GT.zamy)) THEN
-                  vinfor(21) = vinfor(21) + aire(ji,jj) * 1.e-12_wp ! Seasonal ice extent
-               ENDIF
-               IF ((at_i(ji,jj).GT.0.15).AND.(zafy.LE.zamy)) THEN
-                  vinfor(23) = vinfor(23) + aire(ji,jj) * 1.e-12_wp ! Perennial ice extent
-               ENDIF
-            ENDIF
-         END DO
-      END DO
-      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(25))) !=0 if no multiyear ice 1 if yes
-      vinfor(49) = zindb*vinfor(49) / MAX(vinfor(25),epsi06)
-      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(27))) !=0 if no multiyear ice 1 if yes
-      vinfor(51) = zindb*vinfor(51) / MAX(vinfor(27),epsi06)
-
-      IF( cp_cfg == "orca" ) THEN   !* ORCA configuration : Fram Strait Export
-         SELECT CASE ( jp_cfg )
-         CASE ( 2 )                          ! ORCA_R2
-            ij0 = 136   ;   ij1 = 136              ! Fram strait : 83 = area export
-            ii0 = 134   ;   ii1 = 138              !               84 = volume export
-            DO jj = mj0(ij0),mj1(ij1)
-               DO ji = mi0(ii0),mi1(ii1)
-                  vinfor(83) = vinfor(83) - v_ice(ji,jj) * e1t(ji,jj)*at_i(ji,jj)*rdt_ice * 1.e-12_wp
-                  vinfor(84) = vinfor(84) - v_ice(ji,jj) * e1t(ji,jj)*vt_i(ji,jj)*rdt_ice * 1.e-12_wp
-               END DO
-            END DO
-         END SELECT
-!!gm   just above, this is NOT the correct way of evaluating the transport !
-!!gm        mass of snow is missing and v_ice should be the mean between jj and jj+1
-!!gm   Other ORCA configurations should be added
-      ENDIF
-
-      !!-------------------------------------------------------------------
-      !! 2) Southern hemisphere
-      !!-------------------------------------------------------------------
-      !! 2.1) Diagnostics independent on age
-      !!------------------------------------
-      DO jj = 2, njeqm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            IF( tms(ji,jj) == 1 ) THEN
-               vinfor(4)  = vinfor(4)  + at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !ice area
-               IF (at_i(ji,jj).GT.0.15) vinfor(6) = vinfor(6) + aire(ji,jj) * 1.e-12_wp !ice extent
-               vinfor(8)  = vinfor(8)  + vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !ice volume
-               vinfor(10) = vinfor(10) + vt_s(ji,jj)*aire(ji,jj) * 1.e-12_wp !snow volume
-               vinfor(16) = vinfor(16) + ot_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !mean age
-               vinfor(30) = vinfor(30) + smt_i(ji,jj)*vt_i(ji,jj)*aire(ji,jj) * 1.e-12_wp !mean salinity
-               ! this diagnostic is not well computed (weighted by vol instead
-               ! of area)
-               vinfor(32) = vinfor(32) + vt_i(ji,jj)*( u_ice(ji,jj)*u_ice(ji,jj) + & 
-                  v_ice(ji,jj)*v_ice(ji,jj) )*aire(ji,jj)/1.0e12 !ice vel
-               vinfor(54) = vinfor(54) + sfx (ji,jj)*aire(ji,jj) * 1.e-12_wp ! Total salt flux
-               vinfor(56) = vinfor(56) + sfx_bri(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Brine drainage salt flux
-               vinfor(58) = vinfor(58) + sfx_thd(ji,jj)*aire(ji,jj) * 1.e-12_wp ! Equivalent salt flux
-               vinfor(60) = vinfor(60) +(sst_m(ji,jj)+rt0)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp  !SST
-               vinfor(62) = vinfor(62) + sss_m(ji,jj)*at_i(ji,jj)*aire(ji,jj) * 1.e-12_wp  !SSS
-               vinfor(66) = vinfor(66) + et_s(ji,jj)/1.0e9*aire(ji,jj) * 1.e-12_wp ! snow temperature
-               vinfor(68) = vinfor(68) + et_i(ji,jj)/1.0e9*aire(ji,jj) * 1.e-12_wp ! ice enthalpy
-               vinfor(70) = vinfor(70) + v_i(ji,jj,1)*aire(ji,jj) * 1.e-12_wp !ice volume
-               vinfor(72) = vinfor(72) + v_i(ji,jj,2)*aire(ji,jj) * 1.e-12_wp !ice volume
-               vinfor(74) = vinfor(74) + v_i(ji,jj,3)*aire(ji,jj) * 1.e-12_wp !ice volume
-               vinfor(76) = vinfor(76) + v_i(ji,jj,4)*aire(ji,jj) * 1.e-12_wp !ice volume
-               vinfor(78) = vinfor(78) + v_i(ji,jj,5)*aire(ji,jj) * 1.e-12_wp !ice volume
-               vinfor(80) = 0.0
-               vinfor(82) = vinfor(82) + emp(ji,jj)*aire(ji,jj) * 1.e-12_wp ! mass flux
-            ENDIF
-         END DO
-      END DO
-
-      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
-         DO jj = 2, njeqm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               vinfor(12) = vinfor(12) + v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !undef def ice volume
-            END DO
-         END DO
-      END DO
-
-      vinfor(14) = 0.0
-
-      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(8))) 
-      vinfor(16) = zindb * vinfor(16) / MAX(vinfor(8),epsi06) ! these have to be divided by ice vol
-      vinfor(30) = zindb * vinfor(30) / MAX(vinfor(8),epsi06) ! 
-      vinfor(32) = zindb * SQRT( vinfor(32) / MAX( vinfor(8) , epsi06 ) )
-      vinfor(68) = zindb * vinfor(68) / MAX(vinfor(8),epsi06) ! 
-
-      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(6))) 
-      vinfor(54) = zindb * vinfor(54) / MAX(vinfor(6),epsi06) ! these have to be divided by ice extt
-      vinfor(56) = zindb * vinfor(56) / MAX(vinfor(6),epsi06) ! 
-      vinfor(58) = zindb * vinfor(58) / MAX(vinfor(6),epsi06) ! 
-      vinfor(80) = zindb * vinfor(80) / MAX(vinfor(6),epsi06) !
-      !      vinfor(84) = vinfor(84) / vinfor(6) !
-
-      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(4))) !
-      vinfor(60) = zindb*vinfor(60) / ( MAX(vinfor(4), epsi06) ) ! divide by ice area
-      vinfor(62) = zindb*vinfor(62) / ( MAX(vinfor(4), epsi06) ) !
-
-      zindb      = 1.0 - MAX(0.0,SIGN(1.0,-vinfor(10))) !
-      vinfor(66) = zindb*vinfor(66) / MAX(vinfor(10),epsi06) ! divide it by snow volume
-
-      DO jl = 1, jpl
-         DO jj = 2, njeqm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               IF( tms(ji,jj) == 1 ) THEN
-                  vinfor(34) = vinfor(34) + d_v_i_trp(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !ice volume
-                  vinfor(36) = vinfor(36) + d_v_i_thd(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !ice volume
-               ENDIF
-            END DO
-         END DO
-      END DO
-
-      DO jj = 2, njeqm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            IF( tms(ji,jj) == 1 ) THEN
-               vinfor(38) = vinfor(38) + diag_sni_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp !th growth rates
-               vinfor(40) = vinfor(40) + diag_lat_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp 
-               vinfor(42) = vinfor(42) + diag_bot_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp
-               vinfor(44) = vinfor(44) + diag_dyn_gr(ji,jj)*aire(ji,jj) * 1.e-12_wp 
-               vinfor(46) = vinfor(46) + dv_dt_thd(ji,jj,5)*aire(ji,jj) * 1.e-12_wp
-               vinfor(48) = vinfor(48) + v_newice(ji,jj) *aire(ji,jj) * 1.e-12_wp * r1_rdtice   ! volume acc in OW
-            ENDIF
-         END DO
-      END DO
-
-      DO jl = 1, jpl
-         DO jj = 2, njeqm1
-            DO ji = fs_2, fs_jpim1   ! vector opt.
-               IF( tms(ji,jj) == 1 ) THEN
-                  vinfor(64) = vinfor(64) + t_su(ji,jj,jl)*a_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp
-               ENDIF
-            END DO
-         END DO
-      END DO
-      zindb      = 1._wp - MAX(  0._wp , SIGN( 1._wp , -vinfor(4) )  )   !
-      vinfor(64) = zindb * vinfor(64) / MAX( vinfor(4) , epsi06 )   ! divide by ice extt
-      !! 2.2) Diagnostics dependent on age
-      !!------------------------------------
-      DO jj = 2, njeqm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            IF( tms(ji,jj) == 1 ) THEN
-               zafy = 0._wp
-               zamy = 0._wp
-               DO jl = 1, jpl
-                  IF( (o_i(ji,jj,jl) - zshift_date) < 0._wp ) THEN
-                     vinfor(18) = vinfor(18) + a_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp ! FY ice area
-                     vinfor(26) = vinfor(26) + v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp ! FY ice volume
-                     zafy = zafy + a_i(ji,jj,jl)
-                     vinfor(50) = vinfor(50) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !FY ice salinity
-                  ENDIF
-                  IF( (o_i(ji,jj,jl) - zshift_date) > 0._wp ) THEN
-                     vinfor(20) = vinfor(20) + a_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp    ! MY ice area
-                     vinfor(28) = vinfor(28) + v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp
-                     vinfor(52) = vinfor(52) + sm_i(ji,jj,jl)*v_i(ji,jj,jl)*aire(ji,jj) * 1.e-12_wp !FY ice salinity
-                     zamy = zamy + a_i(ji,jj,jl)
-                  ENDIF
-               END DO ! jl
-               IF ( at_i(ji,jj)  >  0.15  .AND. zafy  >  zamy ) THEN
-                  vinfor(22) = vinfor(22) + aire(ji,jj) * 1.e-12_wp ! Seasonal ice extent
-               ENDIF
-               IF ( at_i(ji,jj)  >  0.15  .AND. zafy <= zamy ) THEN
-                  vinfor(24) = vinfor(24) + aire(ji,jj) * 1.e-12_wp ! Perennial ice extent
-               ENDIF
-            ENDIF ! tms
-         END DO ! jj
-      END DO ! ji
-      zindb      = 1.0 - MAX(  0.0,SIGN( 1._wp , -vinfor(26) )  )     !=0 if no multiyear ice 1 if yes
-      vinfor(50) = zindb * vinfor(50) / MAX( vinfor(26) , epsi06 )
-      zindb      = 1.0 - MAX(  0._wp , SIGN( 1._wp , -vinfor(28) )  ) !=0 if no multiyear ice 1 if yes
-      vinfor(52) = zindb * vinfor(52) / MAX( vinfor(28) , epsi06 )
-
-      !  Accumulation before averaging 
-      DO jv = 1, nvinfo
-         vinfom(jv) = vinfom(jv) + vinfor(jv)
-      END DO
-      naveg = naveg + 1  
-
-      ! oututs on file ice_evolu    
-      !MV      IF( MOD( numit , ninfo ) == 0 ) THEN
-      WRITE(numevo_ice,fmtw) ( titvar(jv), vinfom(jv)/naveg, jv = 1, nvinfo )
-      naveg = 0
-      DO jv = 1, nvinfo
-         vinfom(jv) = 0._wp
-      END DO
-      !MV      ENDIF
-      !
-   END SUBROUTINE lim_dia
-
-
-   SUBROUTINE lim_dia_init
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE lim_dia_init  ***
-      !!             
-      !! ** Purpose : Preparation of the file ice_evolu for the output of
-      !!      the temporal evolution of key variables
-      !!
-      !! ** input   : Namelist namicedia
-      !!-------------------------------------------------------------------
-      INTEGER  ::   jv    ! dummy loop indice
-      INTEGER  ::   ierr, ntot , ndeb , irecl   ! local integers
-      REAL(wp) ::   zxx0, zxx1    ! local scalars
-      CHARACTER(len=jpchinf) ::   titinf
-      CHARACTER(len=50)      ::   clname
-      !!
-      NAMELIST/namicedia/fmtinf, nfrinf, ninfo, ntmoy
-      !!-------------------------------------------------------------------
-      !
-      REWIND( numnam_ice )             ! read namicedia namelist
-      READ  ( numnam_ice, namicedia )
-      !
-      IF(lwp) THEN                     ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_dia_init : ice parameters for ice diagnostics '
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         WRITE(numout,*) '   format of the output values                                 fmtinf = ', fmtinf
-         WRITE(numout,*) '   number of variables written in one line                     nfrinf = ', nfrinf 
-         WRITE(numout,*) '   Instantaneous values of ice evolution or averaging          ntmoy  = ', ntmoy
-         WRITE(numout,*) '   frequency of ouputs on file ice_evolu in case of averaging  ninfo  = ', ninfo
-      ENDIF
-
-      ALLOCATE( aire(jpi,jpj) , STAT=ierr )      ! masked grid cell area (interior domain only)
-      IF( lk_mpp    )   CALL mpp_sum( ierr )
-      IF( ierr /= 0 )   CALL ctl_stop( 'STOP', 'lim_dia_init_2 : unable to allocate arrays' )
-      aire(:,:) = area(:,:) * tms(:,:) * tmask_i(:,:)
-
-      ! Titles of ice key variables :
-      titvar(1) = 'NoIt'  ! iteration number
-      titvar(2) = 'T yr'  ! time step in years
-      nbvt      = 2       ! number of time variables
-
-      titvar(3) = 'AI_N'  ! sea ice area in the northern Hemisp.(10^12 km2)
-      titvar(4) = 'AI_S'  ! sea ice area in the southern Hemisp.(10^12 km2)
-      titvar(5) = 'EI_N'  ! sea ice extent (15%) in the northern Hemisp.(10^12 km2)
-      titvar(6) = 'EI_S'  ! sea ice extent (15%) in the southern Hemisp.(10^12 km2)
-      titvar(7) = 'VI_N'  ! sea ice volume in the northern Hemisp.(10^3 km3)
-      titvar(8) = 'VI_S'  ! sea ice volume in the southern Hemisp.(10^3 km3)
-      titvar(9) = 'VS_N'  ! snow volume over sea ice in the northern Hemisp.(10^3 km3)
-      titvar(10)= 'VS_S'  ! snow volume over sea ice in the northern Hemisp.(10^3 km3)
-      titvar(11)= 'VuIN'  ! undeformed sea ice volume in the northern Hemisp.(10^3 km3)
-      titvar(12)= 'VuIS'  ! undeformed sea ice volume in the southern Hemisp.(10^3 km3)
-      titvar(13)= 'VdIN'  ! deformed sea ice volume in the northern Hemisp.(10^3 km3)
-      titvar(14)= 'VdIS'  ! deformed sea ice volume in the southern Hemisp.(10^3 km3)
-      titvar(15)= 'OI_N'  ! sea ice mean age in the northern Hemisp.(years)
-      titvar(16)= 'OI_S'  ! sea ice mean age in the southern Hemisp.(years)
-      titvar(17)= 'AFYN'  ! total FY ice area northern Hemisp.(10^12 km2)
-      titvar(18)= 'AFYS'  ! total FY ice area southern Hemisp.(10^12 km2)
-      titvar(19)= 'AMYN'  ! total MY ice area northern Hemisp.(10^12 km2)
-      titvar(20)= 'AMYS'  ! total MY ice area southern Hemisp.(10^12 km2)
-      titvar(21)= 'EFYN'  ! total FY ice extent northern Hemisp.(10^12 km2) (with more 50% FY ice)
-      titvar(22)= 'EFYS'  ! total FY ice extent southern Hemisp.(10^12 km2) (with more 50% FY ice)
-      titvar(23)= 'EMYN'  ! total MY ice extent northern Hemisp.(10^12 km2) (with more 50% MY ice)
-      titvar(24)= 'EMYS'  ! total MY ice extent southern Hemisp.(10^12 km2) (with more 50% MY ice)
-      titvar(25)= 'VFYN'  ! total undeformed FY ice volume northern Hemisp.(10^3 km3) 
-      titvar(26)= 'VFYS'  ! total undeformed FY ice volume southern Hemisp.(10^3 km3)
-      titvar(27)= 'VMYN'  ! total undeformed MY ice volume northern Hemisp.(10^3 km3) 
-      titvar(28)= 'VMYS'  ! total undeformed MY ice volume southern Hemisp.(10^3 km3) 
-      titvar(29)= 'IS_N'  ! sea ice mean salinity in the northern hemisphere (ppt)  
-      titvar(30)= 'IS_S'  ! sea ice mean salinity in the southern hemisphere (ppt)  
-      titvar(31)= 'IVeN'  ! sea ice mean velocity in the northern hemisphere (m/s) 
-      titvar(32)= 'IVeS'  ! sea ice mean velocity in the southern hemisphere (m/s) 
-      titvar(33)= 'DVDN'  ! variation of sea ice volume due to dynamics in the northern hemisphere
-      titvar(34)= 'DVDS'  ! variation of sea ice volume due to dynamics in the southern hemisphere
-      titvar(35)= 'DVTN'  ! variation of sea ice volume due to thermo in the   northern hemisphere
-      titvar(36)= 'DVTS'  ! variation of sea ice volume due to thermo in the   southern hemisphere
-      titvar(37)= 'TG1N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 1  
-      titvar(38)= 'TG1S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 1  
-      titvar(39)= 'TG2N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 2  
-      titvar(40)= 'TG2S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 2  
-      titvar(41)= 'TG3N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 3  
-      titvar(42)= 'TG3S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 3  
-      titvar(43)= 'TG4N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 4  
-      titvar(44)= 'TG4S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 4  
-      titvar(45)= 'TG5N'  ! thermodynamic vertical growth rate in the northern hemisphere, cat 5  
-      titvar(46)= 'TG5S'  ! thermodynamic vertical growth rate in the souhtern hemisphere, cat 5  
-      titvar(47)= 'LA_N'  ! lateral accretion growth rate, northern hemisphere
-      titvar(48)= 'LA_S'  ! lateral accretion growth rate, southern hemisphere 
-      titvar(49)= 'SF_N'  ! Salinity FY, NH 
-      titvar(50)= 'SF_S'  ! Salinity FY, SH 
-      titvar(51)= 'SF_N'  ! Salinity MY, NH 
-      titvar(52)= 'SF_S'  ! Salinity MY, SH 
-      titvar(53)= 'Fs_N'  ! Total salt flux NH
-      titvar(54)= 'Fs_S'  ! Total salt flux SH
-      titvar(55)= 'FsbN'  ! Salt - brine drainage flux NH
-      titvar(56)= 'FsbS'  ! Salt - brine drainage flux SH
-      titvar(57)= 'FseN'  ! Salt - Equivalent salt flux NH
-      titvar(58)= 'FseS'  ! Salt - Equivalent salt flux SH
-      titvar(59)= 'SSTN'  ! SST, NH
-      titvar(60)= 'SSTS'  ! SST, SH
-      titvar(61)= 'SSSN'  ! SSS, NH
-      titvar(62)= 'SSSS'  ! SSS, SH
-      titvar(63)= 'TsuN'  ! Tsu, NH
-      titvar(64)= 'TsuS'  ! Tsu, SH
-      titvar(65)= 'TsnN'  ! Tsn, NH
-      titvar(66)= 'TsnS'  ! Tsn, SH
-      titvar(67)= 'ei_N'  ! ei, NH
-      titvar(68)= 'ei_S'  ! ei, SH
-      titvar(69)= 'vi1N'  ! vi1, NH
-      titvar(70)= 'vi1S'  ! vi1, SH
-      titvar(71)= 'vi2N'  ! vi2, NH
-      titvar(72)= 'vi2S'  ! vi2, SH
-      titvar(73)= 'vi3N'  ! vi3, NH
-      titvar(74)= 'vi3S'  ! vi3, SH
-      titvar(75)= 'vi4N'  ! vi4, NH
-      titvar(76)= 'vi4S'  ! vi4, SH
-      titvar(77)= 'vi5N'  ! vi5, NH
-      titvar(78)= 'vi5S'  ! vi5, SH
-      titvar(79)= 'vi6N'  ! vi6, NH
-      titvar(80)= 'vi6S'  ! vi6, SH
-      titvar(81)= 'fmaN'  ! mass flux in the ocean, NH
-      titvar(82)= 'fmaS'  ! mass flux in the ocean, SH
-      titvar(83)= 'AFSE'  ! Fram Strait Area export
-      titvar(84)= 'VFSE'  ! Fram Strait Volume export
-      nvinfo = 84
-
-      ! Definition et Ecriture de l'entete : nombre d'enregistrements 
-      ndeb   = ( nstart - 1 ) / ninfo
-      IF( nstart == 1 ) ndeb = -1
-
-      nferme = ( nstart - 1 + nitrun) / ninfo
-      ntot   = nferme - ndeb
-      ndeb   = ninfo * ( 1 + ndeb )
-      nferme = ninfo * nferme
-
-      ! definition of formats 
-      WRITE( fmtw  , '(A,I3,A2,I1,A)' )  '(', nfrinf, '(A', jpchsep, ','//fmtinf//'))'
-      WRITE( fmtr  , '(A,I3,A,I1,A)'  )  '(', nfrinf, '(', jpchsep, 'X,'//fmtinf//'))'
-      WRITE( fmtitr, '(A,I3,A,I1,A)'  )  '(', nvinfo, 'A', jpchinf, ')'
-
-      ! opening  "ice_evolu" file
-      IF( lk_mpp ) THEN   ;   WRITE(clname,FMT="('ice.evolu_',I4.4)") narea-1
-      ELSE                ;   clname = 'ice.evolu'
-      END IF
-      irecl = ( jpchinf + 1 ) * nvinfo 
-      CALL ctl_opn( numevo_ice, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL',    &
-         &         irecl, numout, lwp, narea )
-
-      !- ecriture de 2 lignes d''entete :
-      WRITE(numevo_ice,1000) fmtr, fmtw, fmtitr, nvinfo, ntot, 0, nfrinf
-      zxx0 = 0.001 * REAL(ninfo)
-      zxx1 = 0.001 * REAL(ndeb)
-      WRITE(numevo_ice,1111) REAL(jpchinf), 0., zxx1, zxx0, 0., 0., 0
-
-      !- ecriture de 2 lignes de titre :
-      WRITE(numevo_ice,'(A,I8,A,I8,A,I5)')                 &
-         'Evolution chronologique - Experience '//cexper   &
-         //'   de', ndeb, ' a', nferme, ' pas', ninfo
-      WRITE(numevo_ice,fmtitr) ( titvar(jv), jv = 1, nvinfo )
-
-      !--preparation de "titvar" pour l''ecriture parmi les valeurs numeriques :
-      DO jv = 2 , nvinfo
-         titinf     = titvar(jv)(:jpchinf)
-         titvar(jv) = '  '//titinf
-      END DO
-
-      !--Initialisation of the arrays for the accumulation
-      DO jv = 1, nvinfo
-         vinfom(jv) = 0._wp
-      END DO
-      naveg = 0
-
-1000  FORMAT( 3(A20),4(1x,I6) )
-1111  FORMAT( 3(F7.1,1X,F7.3,1X),I3,A )  
-      !
-   END SUBROUTINE lim_dia_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :                             NO LIM-3 sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_dia         ! Empty routine
-   END SUBROUTINE lim_dia
-#endif
-
-   !!======================================================================
-END MODULE limdia
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4099)
@@ -19,4 +19,5 @@
    USE dom_oce          ! ocean domain
    USE sbc_oce          ! Surface boundary condition: ocean fields
+   USE sbc_ice          ! Surface boundary condition: ice fields
    USE eosbn2           ! equation of state
    USE ice              ! sea-ice variables
@@ -27,6 +28,6 @@
    USE lbclnk           ! lateral boundary condition - MPP exchanges
    USE lib_mpp          ! MPP library
+   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
    USE wrk_nemo         ! work arrays
-   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
 
    IMPLICIT NONE
@@ -35,24 +36,22 @@
    PUBLIC   lim_istate      ! routine called by lim_init.F90
 
-   !                                  !!** init namelist (namiceini) **
-   REAL(wp) ::   ttest    = 2.0_wp     ! threshold water temperature for initial sea ice
-   REAL(wp) ::   hninn    = 0.5_wp     ! initial snow thickness in the north
-   REAL(wp) ::   hginn_u  = 2.5_wp     ! initial ice thickness in the north
-   REAL(wp) ::   aginn_u  = 0.7_wp     ! initial leads area in the north
-   REAL(wp) ::   hginn_d  = 5.0_wp     ! initial ice thickness in the north
-   REAL(wp) ::   aginn_d  = 0.25_wp    ! initial leads area in the north
-   REAL(wp) ::   hnins    = 0.1_wp     ! initial snow thickness in the south
-   REAL(wp) ::   hgins_u  = 1.0_wp     ! initial ice thickness in the south
-   REAL(wp) ::   agins_u  = 0.7_wp     ! initial leads area in the south
-   REAL(wp) ::   hgins_d  = 2.0_wp     ! initial ice thickness in the south
-   REAL(wp) ::   agins_d  = 0.2_wp     ! initial leads area in the south
-   REAL(wp) ::   sinn     = 6.301_wp   ! initial salinity 
-   REAL(wp) ::   sins     = 6.301_wp   !
-
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)
+   !! * Module variables
+   REAL(wp) ::             & !!! ** init namelist (namiceini) **
+      ttest    = 2.0  ,    &  ! threshold water temperature for initial sea ice
+      hninn    = 0.5  ,    &  ! initial snow thickness in the north
+      hnins    = 0.1  ,    &  ! initial snow thickness in the south
+      hginn    = 2.5  ,    &  ! initial ice thickness in the north
+      hgins    = 1.0  ,    &  ! initial ice thickness in the south
+      aginn    = 0.7  ,    &  ! initial leads area in the north
+      agins    = 0.7  ,    &  ! initial leads area in the south
+      sinn     = 6.301 ,   &  ! initial salinity 
+      sins     = 6.301
+
+   !!----------------------------------------------------------------------
+   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008)
    !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
 CONTAINS
 
@@ -63,22 +62,49 @@
       !! ** Purpose :   defined the sea-ice initial state
       !!
-      !! ** Method  :   restart from a state defined in a binary file
-      !!                or from arbitrary sea-ice conditions
-      !!-------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices
-      REAL(wp) ::   zeps6, zeps, ztmelts, epsi06   ! local scalars
-      REAL(wp) ::   zvol, zare, zh, zh1, zh2, zh3, zan, zbn, zas, zbs 
-      REAL(wp), POINTER, DIMENSION(:)   ::   zgfactorn, zhin 
-      REAL(wp), POINTER, DIMENSION(:)   ::   zgfactors, zhis
-      REAL(wp), POINTER, DIMENSION(:,:) ::   zidto      ! ice indicator
-      !--------------------------------------------------------------------
-
-      CALL wrk_alloc( jpm, zgfactorn, zgfactors, zhin, zhis )
+      !! ** Method  :   
+      !!                This routine will put some ice where ocean
+      !!                is at the freezing point, then fill in ice 
+      !!                state variables using prescribed initial 
+      !!                values in the namelist            
+      !!
+      !! ** Steps   :   
+      !!                1) Read namelist
+      !!                2) Basal temperature; ice and hemisphere masks
+      !!                3) Fill in the ice thickness distribution using gaussian
+      !!                4) Fill in space-dependent arrays for state variables
+      !!                5) Diagnostic arrays
+      !!                6) Lateral boundary conditions
+      !!
+      !! History :
+      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code
+      !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats
+      !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats
+      !!--------------------------------------------------------------------
+
+      !! * Local variables
+      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices
+      REAL(wp)   :: epsi06, epsi20, ztmelts
+      INTEGER    :: i_hemis, i_fill, jl0  
+      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv
+      REAL(wp), POINTER, DIMENSION(:)     :: zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zht_i_ini, za_i_ini, zv_i_ini
+      REAL(wp), POINTER, DIMENSION(:,:)   :: zidto    ! ice indicator
+      INTEGER,  POINTER, DIMENSION(:,:)   :: zhemis   ! hemispheric index
+      !--------------------------------------------------------------------
+
       CALL wrk_alloc( jpi, jpj, zidto )
-
-      !--------------------------------------------------------------------
-      ! 1) Preliminary things 
-      !--------------------------------------------------------------------
-      epsi06 = 1.e-6_wp
+      CALL wrk_alloc( jpi, jpj, zhemis )
+      CALL wrk_alloc( jpl,   2, zht_i_ini,  za_i_ini,  zv_i_ini )
+      CALL wrk_alloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini )
+
+      epsi06   = 1.0e-6
+      epsi20   = 1.0e-20
+      IF(lwp) WRITE(numout,*)
+      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
+      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
+
+      !--------------------------------------------------------------------
+      ! 1) Read namelist
+      !--------------------------------------------------------------------
 
       CALL lim_istate_init     !  reading the initials parameters of the ice
@@ -89,408 +115,308 @@
 
       !--------------------------------------------------------------------
-      ! 2) Ice initialization (hi,hs,frld,t_su,sm_i,t_i,t_s)              | 
-      !--------------------------------------------------------------------
-
-      IF(lwp) WRITE(numout,*)
-      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization '
-      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ '
-
+      ! 2) Basal temperature, ice mask and hemispheric index
+      !--------------------------------------------------------------------
+
+      ! Basal temperature is set to the freezing point of seawater in Celsius
       t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius]
 
       DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest
          DO ji = 1, jpi
-            IF( tsn(ji,jj,1,jp_tem)  - t_bo(ji,jj) >= ttest ) THEN   ;   zidto(ji,jj) = 0.e0      ! no ice
-            ELSE                                                     ;   zidto(ji,jj) = 1.e0      !    ice
+            IF( tsn(ji,jj,1,jp_tem)  - t_bo(ji,jj) >= ttest ) THEN   ;   zidto(ji,jj) = 0._wp      ! no ice
+            ELSE                                                     ;   zidto(ji,jj) = 1._wp      !    ice
             ENDIF
          END DO
       END DO
 
-      t_bo(:,:) = t_bo(:,:) + rt0                          ! t_bo converted from Celsius to Kelvin (rt0 over land)
-
-      ! constants for heat contents
-      zeps   = 1.e-20_wp
-      zeps6  = 1.e-06_wp
-
-      ! zgfactor for initial ice distribution
-      zgfactorn(:) = 0._wp
-      zgfactors(:) = 0._wp
-
-      ! first ice type
-      DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2)
-         zhin (1)     = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
-         zgfactorn(1) = zgfactorn(1) + exp(-(zhin(1)-hginn_u)*(zhin(1)-hginn_u) * 0.5_wp )
-         zhis (1)     = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
-         zgfactors(1) = zgfactors(1) + exp(-(zhis(1)-hgins_u)*(zhis(1)-hgins_u) * 0.5_wp )
-      END DO ! jl
-      zgfactorn(1) = aginn_u / zgfactorn(1)
-      zgfactors(1) = agins_u / zgfactors(1)
-
-      ! -------------
-      ! new distribution, polynom of second order, conserving area and volume
-      zh1 = 0._wp
-      zh2 = 0._wp
-      zh3 = 0._wp
-      DO jl = 1, jpl
-         zh = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp
-         zh1 = zh1 + zh
-         zh2 = zh2 + zh * zh
-         zh3 = zh3 + zh * zh * zh
-      END DO
-      IF(lwp) WRITE(numout,*) ' zh1 : ', zh1
-      IF(lwp) WRITE(numout,*) ' zh2 : ', zh2
-      IF(lwp) WRITE(numout,*) ' zh3 : ', zh3
-
-      zvol = aginn_u * hginn_u
-      zare = aginn_u
-      IF( jpl >= 2 ) THEN
-         zbn = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3)
-         zan = ( zare - zbn*zh1 ) / zh2
-      ENDIF
-
-      IF(lwp) WRITE(numout,*) ' zvol: ', zvol
-      IF(lwp) WRITE(numout,*) ' zare: ', zare
-      IF(lwp) WRITE(numout,*) ' zbn : ', zbn 
-      IF(lwp) WRITE(numout,*) ' zan : ', zan 
-
-      zvol = agins_u * hgins_u
-      zare = agins_u
-      IF( jpl >= 2 ) THEN
-         zbs = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3)
-         zas = ( zare - zbs*zh1 ) / zh2
-      ENDIF
-
-      IF(lwp) WRITE(numout,*) ' zvol: ', zvol
-      IF(lwp) WRITE(numout,*) ' zare: ', zare
-      IF(lwp) WRITE(numout,*) ' zbn : ', zbn 
-      IF(lwp) WRITE(numout,*) ' zan : ', zan 
-
-      !end of new lines
-      ! -------------
-!!!
-      ! retour a LIMA_MEC
-      !     ! second ice type
-      !     zdummy  = hi_max(ice_cat_bounds(2,1)-1)
-      !     hi_max(ice_cat_bounds(2,1)-1) = 0.0
-
-      !     ! here to change !!!!
-      !     jm = 2
-      !     DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2)
-      !        zhin (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
-      !        zhin (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + &
-      !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0
-      !        zgfactorn(2) = zgfactorn(2) + exp(-(zhin(2)-hginn_d)*(zhin(2)-hginn_d)/2.0)
-      !        zhis (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
-      !        zhis (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + &
-      !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0
-      !        zgfactors(2) = zgfactors(2) + exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0)
-      !     END DO ! jl
-      !     zgfactorn(2) = aginn_d / zgfactorn(2)
-      !     zgfactors(2) = agins_d / zgfactors(2)
-      !     hi_max(ice_cat_bounds(2,1)-1) = zdummy
-      ! END retour a LIMA_MEC
-!!!
-
-!!gm  optimisation :  loop over the ice categories inside the ji, jj loop !!!
-
+      t_bo(:,:) = t_bo(:,:) + rt0                          ! conversion to Kelvin
+
+      ! Hemispheric index
+      ! MV 2011 new initialization
       DO jj = 1, jpj
          DO ji = 1, jpi
-
-            !--- Northern hemisphere
-            !----------------------------------------------------------------
             IF( fcor(ji,jj) >= 0._wp ) THEN    
-
-               !-----------------------
-               ! Ice area / thickness
-               !-----------------------
-
-               IF ( jpl .EQ. 1) THEN ! one category
-
-                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
-                     a_i(ji,jj,jl)    = zidto(ji,jj) * aginn_u
-                     ht_i(ji,jj,jl)   = zidto(ji,jj) * hginn_u
-                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
-                  END DO
-
-               ELSE ! several categories
-
-                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
-                     zhin(1)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
-                     a_i(ji,jj,jl)    = zidto(ji,jj) * MAX( zgfactorn(1) * exp(-(zhin(1)-hginn_u)* & 
-                        (zhin(1)-hginn_u)/2.0) , epsi06)
-                     ! new line
-                     a_i(ji,jj,jl)    = zidto(ji,jj) * ( zan * zhin(1) * zhin(1) + zbn * zhin(1) )
-                     ht_i(ji,jj,jl)   = zidto(ji,jj) * zhin(1) 
-                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
-                  END DO
-
-               ENDIF
-
-
-!!!
-               ! retour a LIMA_MEC
-               !              !ridged ice
-               !              zdummy  = hi_max(ice_cat_bounds(2,1)-1)
-               !              hi_max(ice_cat_bounds(2,1)-1) = 0.0
-               !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) ! loop over ice thickness categories
-               !                 zhin(2)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
-               !                 a_i(ji,jj,jl)    = zidto(ji,jj) * MAX( zgfactorn(2) * exp(-(zhin(2)-hginn_d)* &
-               !                                    (zhin(2)-hginn_d)/2.0) , epsi06)
-               !                 ht_i(ji,jj,jl)   = zidto(ji,jj) * zhin(2) 
-               !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
-               !              END DO
-               !              hi_max(ice_cat_bounds(2,1)-1) = zdummy
-
-               !              !rafted ice
-               !              jl = 6
-               !              a_i(ji,jj,jl)       = 0.0
-               !              ht_i(ji,jj,jl)      = 0.0
-               !              v_i(ji,jj,jl)       = 0.0
-               ! END retour a LIMA_MEC
-!!!
-
-               DO jl = 1, jpl
-
-                  !-------------
-                  ! Snow depth
-                  !-------------
-                  ht_s(ji,jj,jl)   = zidto(ji,jj) * hninn
-                  v_s(ji,jj,jl)    = ht_s(ji,jj,jl)*a_i(ji,jj,jl)
-
-                  !---------------
-                  ! Ice salinity
-                  !---------------
-                  sm_i(ji,jj,jl)   = zidto(ji,jj) * sinn  + ( 1.0 - zidto(ji,jj) ) * 0.1
-                  smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
-
-                  !----------
-                  ! Ice age
-                  !----------
-                  o_i(ji,jj,jl)    = zidto(ji,jj) * 1.0   + ( 1.0 - zidto(ji,jj) )
-                  oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl)
-
-                  !------------------------------
-                  ! Sea ice surface temperature
-                  !------------------------------
-
-                  t_su(ji,jj,jl)   = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj)
-
-                  !------------------------------------
-                  ! Snow temperature and heat content
-                  !------------------------------------
-
-                  DO jk = 1, nlay_s
-                     t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt
-                     ! Snow energy of melting
-                     e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
-                     ! Change dimensions
-                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
-                     ! Multiply by volume, so that heat content in 10^9 Joules
-                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
-                        v_s(ji,jj,jl)  / nlay_s
-                  END DO !jk
-
-                  !-----------------------------------------------
-                  ! Ice salinities, temperature and heat content 
-                  !-----------------------------------------------
-
-                  DO jk = 1, nlay_i
-                     t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 
-                     s_i(ji,jj,jk,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1
-                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
-
-                     ! heat content per unit volume
-                     e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * &
-                        (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
-                        +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) &
-                        - rcp      * ( ztmelts - rtt ) &
-                        )
-
-                     ! Correct dimensions to avoid big values
-                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
-
-                     ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
-                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 
-                        area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / &
-                        nlay_i
-                  END DO ! jk
-
-               END DO ! jl 
-
-            ELSE ! on fcor 
-
-               !--- Southern hemisphere
-               !----------------------------------------------------------------
-
-               !-----------------------
-               ! Ice area / thickness
-               !-----------------------
-
-               IF ( jpl .EQ. 1) THEN ! one category
-
-                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories
-                     a_i(ji,jj,jl)    = zidto(ji,jj) * agins_u
-                     ht_i(ji,jj,jl)   = zidto(ji,jj) * hgins_u
-                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
-                  END DO
-
-               ELSE ! several categories
-
-                  !level ice
-                  DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) !over thickness categories
-
-                     zhis(1)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
-                     a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactors(1) * exp(-(zhis(1)-hgins_u) * & 
-                        (zhis(1)-hgins_u)/2.0) , epsi06 )
-                     ! new line square distribution volume conserving
-                     a_i(ji,jj,jl)    = zidto(ji,jj) * ( zas * zhis(1) * zhis(1) + zbs * zhis(1) )
-                     ht_i(ji,jj,jl)   = zidto(ji,jj) * zhis(1) 
-                     v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
-
-                  END DO ! jl
-
-               ENDIF
-
-!!!
-               ! retour a LIMA_MEC
-               !              !ridged ice
-               !              zdummy  = hi_max(ice_cat_bounds(2,1)-1)
-               !              hi_max(ice_cat_bounds(2,1)-1) = 0.0
-               !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) !over thickness categories
-               !                 zhis(2)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0
-               !                 a_i(ji,jj,jl) = zidto(ji,jj)*MAX( zgfactors(2)   &
-               !                    &          * exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0), epsi06 )
-               !                 ht_i(ji,jj,jl)   = zidto(ji,jj) * zhis(2) 
-               !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl)
-               !              END DO
-               !              hi_max(ice_cat_bounds(2,1)-1) = zdummy
-
-               !              !rafted ice
-               !              jl = 6
-               !              a_i(ji,jj,jl)       = 0.0
-               !              ht_i(ji,jj,jl)      = 0.0
-               !              v_i(ji,jj,jl)       = 0.0
-               ! END retour a LIMA_MEC
-!!!
-
-               DO jl = 1, jpl !over thickness categories
-
-                  !---------------
-                  ! Snow depth
-                  !---------------
-
-                  ht_s(ji,jj,jl)   = zidto(ji,jj) * hnins
-                  v_s(ji,jj,jl)    = ht_s(ji,jj,jl)*a_i(ji,jj,jl)
-
-                  !---------------
-                  ! Ice salinity
-                  !---------------
-
-                  sm_i(ji,jj,jl)   = zidto(ji,jj) * sins  + ( 1.0 - zidto(ji,jj) ) * 0.1
-                  smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
-
-                  !----------
-                  ! Ice age
-                  !----------
-
-                  o_i(ji,jj,jl)    = zidto(ji,jj) * 1.0   + ( 1.0 - zidto(ji,jj) )
-                  oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl)
-
-                  !------------------------------
-                  ! Sea ice surface temperature
-                  !------------------------------
-
-                  t_su(ji,jj,jl)   = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj)
-
-                  !----------------------------------
-                  ! Snow temperature / heat content
-                  !----------------------------------
-
-                  DO jk = 1, nlay_s
-                     t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt
-                     ! Snow energy of melting
-                     e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
-                     ! Change dimensions
-                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
-                     ! Multiply by volume, so that heat content in 10^9 Joules
-                     e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * &
-                        v_s(ji,jj,jl)  / nlay_s
-                  END DO
-
-                  !---------------------------------------------
-                  ! Ice temperature, salinity and heat content
-                  !---------------------------------------------
-
-                  DO jk = 1, nlay_i
-                     t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 
-                     s_i(ji,jj,jk,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1
-                     ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
-
-                     ! heat content per unit volume
-                     e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * &
-                        (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
-                        +   lfus  * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) &
-                        - rcp      * ( ztmelts - rtt ) &
-                        )
-
-                     ! Correct dimensions to avoid big values
-                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
-
-                     ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J
-                     e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * & 
-                        area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / &
-                        nlay_i
-                  END DO !jk
-
-               END DO ! jl 
-
-            ENDIF ! on fcor
-
+               zhemis(ji,jj) = 1 ! Northern hemisphere
+            ELSE
+               zhemis(ji,jj) = 2 ! Southern hemisphere
+            ENDIF
          END DO
       END DO
-
-      !--------------------------------------------------------------------
-      ! 3) Global ice variables for output diagnostics                    | 
-      !--------------------------------------------------------------------
-
-      fsbbq (:,:)     = 0.e0
-      u_ice (:,:)     = 0.e0
-      v_ice (:,:)     = 0.e0
-      stress1_i(:,:)  = 0.0
-      stress2_i(:,:)  = 0.0
-      stress12_i(:,:) = 0.0
-
-      !--------------------------------------------------------------------
-      ! 4) Moments for advection
-      !--------------------------------------------------------------------
-
-      sxopw (:,:) = 0.e0 
-      syopw (:,:) = 0.e0 
-      sxxopw(:,:) = 0.e0 
-      syyopw(:,:) = 0.e0 
-      sxyopw(:,:) = 0.e0
-
-      sxice (:,:,:)  = 0.e0   ;   sxsn (:,:,:)  = 0.e0   ;   sxa  (:,:,:)  = 0.e0
-      syice (:,:,:)  = 0.e0   ;   sysn (:,:,:)  = 0.e0   ;   sya  (:,:,:)  = 0.e0
-      sxxice(:,:,:)  = 0.e0   ;   sxxsn(:,:,:)  = 0.e0   ;   sxxa (:,:,:)  = 0.e0
-      syyice(:,:,:)  = 0.e0   ;   syysn(:,:,:)  = 0.e0   ;   syya (:,:,:)  = 0.e0
-      sxyice(:,:,:)  = 0.e0   ;   sxysn(:,:,:)  = 0.e0   ;   sxya (:,:,:)  = 0.e0
-
-      sxc0  (:,:,:)  = 0.e0   ;   sxe  (:,:,:,:)= 0.e0   
-      syc0  (:,:,:)  = 0.e0   ;   sye  (:,:,:,:)= 0.e0   
-      sxxc0 (:,:,:)  = 0.e0   ;   sxxe (:,:,:,:)= 0.e0   
-      syyc0 (:,:,:)  = 0.e0   ;   syye (:,:,:,:)= 0.e0   
-      sxyc0 (:,:,:)  = 0.e0   ;   sxye (:,:,:,:)= 0.e0   
-
-      sxsal  (:,:,:)  = 0.e0
-      sysal  (:,:,:)  = 0.e0
-      sxxsal (:,:,:)  = 0.e0
-      syysal (:,:,:)  = 0.e0
-      sxysal (:,:,:)  = 0.e0
-
-      !--------------------------------------------------------------------
-      ! 5) Lateral boundary conditions                                    | 
+      ! END MV 2011 new initialization
+
+      !--------------------------------------------------------------------
+      ! 3) Initialization of sea ice state variables
+      !--------------------------------------------------------------------
+
+      !-----------------------------
+      ! 3.1) Hemisphere-dependent arrays
+      !-----------------------------
+      ! assign initial thickness, concentration, snow depth and salinity to
+      ! an hemisphere-dependent array
+      zhm_i_ini(1) = hginn ; zhm_i_ini(2) = hgins  ! ice thickness
+      zat_i_ini(1) = aginn ; zat_i_ini(2) = agins  ! ice concentration
+      zvt_i_ini(:) = zhm_i_ini(:) * zat_i_ini(:)   ! ice volume
+      zhm_s_ini(1) = hninn ; zhm_s_ini(2) = hnins  ! snow depth
+      zsm_i_ini(1) = sinn  ; zsm_i_ini(2) = sins   ! bulk ice salinity
+
+      !---------------------------------------------------------------------
+      ! 3.2) Distribute ice concentration and thickness into the categories
+      !---------------------------------------------------------------------
+      ! a gaussian distribution for ice concentration is used
+      ! then we check whether the distribution fullfills
+      ! volume and area conservation, positivity and ice categories bounds
+      DO i_hemis = 1, 2
+
+      ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0
+
+      ! note for the great nemo engineers: 
+      ! only very few of the WRITE statements are necessary for the reference version
+      ! they were one day useful, but now i personally doubt of their
+      ! potential for bringing anything useful
+
+      DO i_fill = jpl, 1, -1
+         IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN
+            !----------------------------
+            ! fill the i_fill categories
+            !----------------------------
+            ! *** 1 category to fill
+            IF ( i_fill .EQ. 1 ) THEN
+               zht_i_ini(1,i_hemis)       = zhm_i_ini(i_hemis)
+               za_i_ini(1,i_hemis)        = zat_i_ini(i_hemis)
+               zht_i_ini(2:jpl,i_hemis)   = 0._wp
+               za_i_ini(2:jpl,i_hemis)    = 0._wp
+            ELSE
+
+            ! *** >1 categores to fill
+            !--- Ice thicknesses in the i_fill - 1 first categories
+               DO jl = 1, i_fill - 1
+                  zht_i_ini(jl,i_hemis)    = 0.5 * ( hi_max(jl) + hi_max(jl-1) )
+               END DO
+
+            !--- jl0: most likely index where cc will be maximum
+               DO jl = 1, jpl
+                  IF ( ( zhm_i_ini(i_hemis) .GT. hi_max(jl-1) ) .AND. &
+                       ( zhm_i_ini(i_hemis) .LE. hi_max(jl)   ) ) THEN
+                     jl0 = jl
+                  ENDIF
+               END DO
+               jl0 = MIN(jl0, i_fill)
+
+            !--- Concentrations
+               za_i_ini(jl0,i_hemis)      = zat_i_ini(i_hemis) / SQRT(REAL(jpl))
+               DO jl = 1, i_fill - 1
+                  IF ( jl .NE. jl0 ) THEN
+                     zsigma               = 0.5 * zhm_i_ini(i_hemis)
+                     zarg                 = ( zht_i_ini(jl,i_hemis) - zhm_i_ini(i_hemis) ) / zsigma
+                     za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2)
+                  ENDIF
+               END DO 
+
+               zA = 0. ! sum of the areas in the jpl categories 
+               DO jl = 1, i_fill - 1
+                 zA = zA + za_i_ini(jl,i_hemis)
+               END DO
+               za_i_ini(i_fill,i_hemis)   = zat_i_ini(i_hemis) - zA ! ice conc in the last category
+               IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
+         
+            !--- Ice thickness in the last category
+               zV = 0. ! sum of the volumes of the N-1 categories
+               DO jl = 1, i_fill - 1
+                  zV = zV + za_i_ini(jl,i_hemis)*zht_i_ini(jl,i_hemis)
+               END DO
+               zht_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis) 
+               IF ( i_fill .LT. jpl ) zht_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
+
+            !--- volumes
+               zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zht_i_ini(:,i_hemis)
+               IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp
+
+            ENDIF ! i_fill
+
+            !---------------------
+            ! Compatibility tests
+            !---------------------
+            ! Test 1: area conservation
+            zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons )
+            IF ( zconv .LT. 1.0e-6 ) THEN
+               ztest_1 = 1
+            ELSE 
+              ! this write is useful
+              IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis) 
+               ztest_1 = 0
+            ENDIF
+
+            ! Test 2: volume conservation
+            zV_cons = SUM(zv_i_ini(:,i_hemis))
+            zconv = ABS(zvt_i_ini(i_hemis) - zV_cons)
+
+            IF ( zconv .LT. 1.0e-6 ) THEN
+               ztest_2 = 1
+            ELSE
+              ! this write is useful
+              IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, &
+                            ' zvt_i_ini = ', zvt_i_ini(i_hemis)
+               ztest_2 = 0
+            ENDIF
+
+            ! Test 3: thickness of the last category is in-bounds ?
+            IF ( zht_i_ini(i_fill, i_hemis) .GT. hi_max(i_fill-1) ) THEN
+               ztest_3 = 1
+            ELSE
+               ! this write is useful
+               IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zht_i_ini(i_fill,i_hemis) = ', &
+               zht_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1)
+               ztest_3 = 0
+            ENDIF
+
+            ! Test 4: positivity of ice concentrations
+            ztest_4 = 1
+            DO jl = 1, jpl
+               IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN 
+                  ! this write is useful
+                  IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis)
+                  ztest_4 = 0
+               ENDIF
+            END DO
+
+         ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4
+ 
+         ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4
+
+      END DO ! i_fill
+
+      IF(lwp) THEN 
+         WRITE(numout,*), ' ztests : ', ztests
+         IF ( ztests .NE. 4 ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
+            WRITE(numout,*), ' !!!! RED ALERT                  !!! '
+            WRITE(numout,*), ' !!!! BIIIIP BIIIP BIIIIP BIIIIP !!!'
+            WRITE(numout,*), ' !!!! Something is wrong in the LIM3 initialization procedure '
+            WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! '
+            WRITE(numout,*)
+            WRITE(numout,*), ' *** ztests is not equal to 4 '
+            WRITE(numout,*), ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4
+            WRITE(numout,*), ' zat_i_ini : ', zat_i_ini(i_hemis)
+            WRITE(numout,*), ' zhm_i_ini : ', zhm_i_ini(i_hemis)
+         ENDIF ! ztests .NE. 4
+      ENDIF
+      
+      END DO ! i_hemis
+
+      !---------------------------------------------------------------------
+      ! 3.3) Space-dependent arrays for ice state variables
+      !---------------------------------------------------------------------
+
+      ! Ice concentration, thickness and volume, snow depth, ice
+      ! salinity, ice age, surface temperature
+      DO jl = 1, jpl ! loop over categories
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               a_i(ji,jj,jl)   = zidto(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration
+               ht_i(ji,jj,jl)  = zidto(ji,jj) * zht_i_ini(jl,zhemis(ji,jj))  ! ice thickness
+               ht_s(ji,jj,jl)  = zidto(ji,jj) * zhm_s_ini(zhemis(ji,jj))     ! snow depth
+               sm_i(ji,jj,jl)  = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity
+               o_i(ji,jj,jl)   = zidto(ji,jj) * 1._wp + ( 1._wp - zidto(ji,jj) ) ! age
+               t_su(ji,jj,jl)  = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * t_bo(ji,jj) ! surf temp
+ 
+               ! ice volume, snow volume, salt content, age content
+               v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume
+               v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume
+               smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content
+               oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content
+            END DO ! ji
+         END DO ! jj
+      END DO ! jl
+
+      ! Snow temperature and heat content
+      DO jk = 1, nlay_s
+         DO jl = 1, jpl ! loop over categories
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                   t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * rtt
+                   ! Snow energy of melting
+                   e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus )
+                   ! Change dimensions
+                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac
+                   ! Multiply by volume, so that heat content in 10^9 Joules
+                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s
+               END DO ! ji
+            END DO ! jj
+         END DO ! jl
+      END DO ! jk
+
+      ! Ice salinity, temperature and heat content
+      DO jk = 1, nlay_i
+         DO jl = 1, jpl ! loop over categories
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                   t_i(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1._wp - zidto(ji,jj) ) * rtt 
+                   s_i(ji,jj,jk,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min
+                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K
+
+                   ! heat content per unit volume
+                   e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
+                      +   lfus    * ( 1._wp - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) &
+                      -   rcp     * ( ztmelts - rtt ) )
+
+                   ! Correct dimensions to avoid big values
+                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac 
+
+                   ! Mutliply by ice volume, and divide by number of layers 
+                   ! to get heat content in 10^9 J
+                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i
+               END DO ! ji
+            END DO ! jj
+         END DO ! jl
+      END DO ! jk
+
+      !--------------------------------------------------------------------
+      ! 4) Global ice variables for output diagnostics                    | 
+      !--------------------------------------------------------------------
+      fsbbq (:,:)     = 0._wp
+      u_ice (:,:)     = 0._wp
+      v_ice (:,:)     = 0._wp
+      stress1_i(:,:)  = 0._wp
+      stress2_i(:,:)  = 0._wp
+      stress12_i(:,:) = 0._wp
+
+# if defined key_coupled
+      albege(:,:)   = 0.8 * tms(:,:)
+# endif
+
+      !--------------------------------------------------------------------
+      ! 5) Moments for advection
+      !--------------------------------------------------------------------
+
+      sxopw (:,:) = 0._wp 
+      syopw (:,:) = 0._wp 
+      sxxopw(:,:) = 0._wp 
+      syyopw(:,:) = 0._wp 
+      sxyopw(:,:) = 0._wp
+
+      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp
+      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp
+      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp
+      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp
+      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp
+
+      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp   
+      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp   
+      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp   
+      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp   
+      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp   
+
+      sxsal  (:,:,:)  = 0._wp
+      sysal  (:,:,:)  = 0._wp
+      sxxsal (:,:,:)  = 0._wp
+      syysal (:,:,:)  = 0._wp
+      sxysal (:,:,:)  = 0._wp
+
+      !--------------------------------------------------------------------
+      ! 6) Lateral boundary conditions                                    | 
       !--------------------------------------------------------------------
 
       DO jl = 1, jpl
+
          CALL lbc_lnk( a_i(:,:,jl)  , 'T', 1. )
          CALL lbc_lnk( v_i(:,:,jl)  , 'T', 1. )
@@ -498,5 +424,5 @@
          CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. )
          CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. )
-         !
+
          CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. )
          CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. )
@@ -515,4 +441,9 @@
          a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl)
       END DO
+      
+      at_i (:,:) = 0.0_wp
+      DO jl = 1, jpl
+         at_i (:,:) = at_i (:,:) + a_i (:,:,jl)
+      END DO
 
       CALL lbc_lnk( at_i , 'T', 1. )
@@ -521,9 +452,15 @@
       CALL lbc_lnk( fsbbq  , 'T', 1. )
       !
-      CALL wrk_dealloc( jpm, zgfactorn, zgfactors, zhin, zhis )
+      !--------------------------------------------------------------------
+      ! 6) ????                                                           | 
+      !--------------------------------------------------------------------
+      tn_ice (:,:,:) = t_su (:,:,:)
+
       CALL wrk_dealloc( jpi, jpj, zidto )
-      !
+      CALL wrk_dealloc( jpi, jpj, zhemis )
+      CALL wrk_dealloc( jpl,   2, zht_i_ini,  za_i_ini,  zv_i_ini )
+      CALL wrk_dealloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini )
+
    END SUBROUTINE lim_istate
-
 
    SUBROUTINE lim_istate_init
@@ -533,17 +470,24 @@
       !! ** Purpose : Definition of initial state of the ice 
       !!
-      !! ** Method :   Read the namiceini namelist and check the parameter 
-      !!             values called at the first timestep (nit000)
-      !!
-      !! ** input  :   namelist namiceini
+      !! ** Method : Read the namiceini namelist and check the parameter 
+      !!       values called at the first timestep (nit000)
+      !!
+      !! ** input : 
+      !!        Namelist namiceini
+      !!
+      !! history :
+      !!  8.5  ! 03-08 (C. Ethe) original code 
+      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization
       !!-----------------------------------------------------------------------------
-      NAMELIST/namiceini/ ttest, hninn, hginn_u, aginn_u, hginn_d, aginn_d, hnins,   &
-         &                hgins_u, agins_u, hgins_d, agins_d, sinn, sins
+      NAMELIST/namiceini/ ttest, hninn, hnins, hginn, hgins, aginn, agins, sinn, sins
       !!-----------------------------------------------------------------------------
-      !
-      REWIND ( numnam_ice )               ! Read Namelist namiceini 
+
+      ! Define the initial parameters
+      ! -------------------------
+
+      ! Read Namelist namiceini 
+      REWIND ( numnam_ice )
       READ   ( numnam_ice , namiceini )
-      !
-      IF(lwp) THEN                        ! control print
+      IF(lwp) THEN
          WRITE(numout,*)
          WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
@@ -551,17 +495,13 @@
          WRITE(numout,*) '   threshold water temp. for initial sea-ice    ttest      = ', ttest
          WRITE(numout,*) '   initial snow thickness in the north          hninn      = ', hninn
-         WRITE(numout,*) '   initial undef ice thickness in the north     hginn_u    = ', hginn_u
-         WRITE(numout,*) '   initial undef ice concentr. in the north     aginn_u    = ', aginn_u          
-         WRITE(numout,*) '   initial  def  ice thickness in the north     hginn_d    = ', hginn_d
-         WRITE(numout,*) '   initial  def  ice concentr. in the north     aginn_d    = ', aginn_d          
          WRITE(numout,*) '   initial snow thickness in the south          hnins      = ', hnins 
-         WRITE(numout,*) '   initial undef ice thickness in the north     hgins_u    = ', hgins_u
-         WRITE(numout,*) '   initial undef ice concentr. in the north     agins_u    = ', agins_u          
-         WRITE(numout,*) '   initial  def  ice thickness in the north     hgins_d    = ', hgins_d
-         WRITE(numout,*) '   initial  def  ice concentr. in the north     agins_d    = ', agins_d          
-         WRITE(numout,*) '   initial  ice salinity       in the north     sinn       = ', sinn
-         WRITE(numout,*) '   initial  ice salinity       in the south     sins       = ', sins
+         WRITE(numout,*) '   initial ice thickness  in the north          hginn      = ', hginn
+         WRITE(numout,*) '   initial ice thickness  in the south          hgins      = ', hgins
+         WRITE(numout,*) '   initial ice concentr.  in the north          aginn      = ', aginn
+         WRITE(numout,*) '   initial ice concentr.  in the north          agins      = ', agins
+         WRITE(numout,*) '   initial  ice salinity  in the north          sinn       = ', sinn
+         WRITE(numout,*) '   initial  ice salinity  in the south          sins       = ', sins
       ENDIF
-      !
+
    END SUBROUTINE lim_istate_init
 
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4099)
@@ -416,9 +416,11 @@
 
                delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zzdst*zzdst ) * usecc2 )  
-               deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl )
-!!gm faster to replace the line above with simply:
-!!                deltat(ji,jj) = MAX( delta, creepl )
-!!gm end 
-
+               ! MV rewriting
+               ! deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl )
+               !!gm faster to replace the line above with simply:
+               !!                deltat(ji,jj) = MAX( delta, creepl )
+               !!gm end  
+               deltat(ji,jj) = delta + creepl
+               ! END MV
                !-Calculate stress tensor components zs1 and zs2 
                !-at centre of grid cells (see section 3.5 of CICE user's guide).
@@ -741,8 +743,12 @@
                   &           - e1v( ji  , jj-1 ) * u_ice2(ji  ,jj-1)  ) / area(ji,jj)
 
-               deltat(ji,jj) = SQRT(    zdd(ji,jj)*zdd(ji,jj)   & 
-                  &                 + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 & 
-                  &                          ) + creepl
-
+!              deltat(ji,jj) = SQRT(    zdd(ji,jj)*zdd(ji,jj)   & 
+!                  &                 + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 & 
+!                  &                          ) + creepl
+               ! MV rewriting
+               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 )  
+               deltat(ji,jj) = delta + creepl
+               ! END MV
+            
             ENDIF ! zdummy
 
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 4099)
@@ -371,4 +371,5 @@
          CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
          t_su(:,:,jl) = z2d(:,:)
+         tn_ice (:,:,:) = t_su (:,:,:)
       END DO
 
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4099)
@@ -102,5 +102,5 @@
       INTEGER, INTENT(in) ::   kt    ! number of iteration
       !
-      INTEGER  ::   ji, jj           ! dummy loop indices
+      INTEGER  ::   ji, jj, jl           ! dummy loop indices
       INTEGER  ::   ierr, ifvt, i1mfr, idfr           ! local integer
       INTEGER  ::   iflt, ial , iadv , ifral, ifrdv   !   -      -
@@ -145,6 +145,16 @@
 
             !   computation the solar flux at ocean surface
-            zfcm1   = pfrld(ji,jj) * qsr(ji,jj)  + &
-                 &    ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) )
+            IF (lk_cpl) THEN ! be carfeful: not being tested yet
+               ! original line
+               !zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj)
+               ! new line to include solar penetration (not tested)
+               zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) )
+               DO jl = 1, jpl
+                  zfcm1 = zfcm1 - qsr_ice(ji,jj,jl) * a_i(ji,jj,jl)
+               END DO
+            ELSE
+               zfcm1   = pfrld(ji,jj) * qsr(ji,jj)  + &
+                    &    ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) )
+            ENDIF
             ! fstric     Solar flux transmitted trough the ice
             ! qsr        Net short wave heat flux on free ocean
@@ -227,8 +237,13 @@
 
             !  computing freshwater exchanges at the ice/ocean interface
-            zemp =   emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction
-               &   - tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean
-               &   + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice
-               &   - fmmec(ji,jj)                                         ! snow falling when ridging
+            IF (lk_cpl) THEN 
+               zemp = - emp_tot(ji,jj) + emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !
+                  &   - rdm_snw(ji,jj) / rdt_ice
+            ELSE
+               zemp =   emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction
+                  &   - tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean
+                  &   + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice
+                  &   - fmmec(ji,jj)                                         ! snow falling when ridging
+            ENDIF
 
             ! mass flux at the ocean/ice interface (sea ice fraction)
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90	(revision 4099)
@@ -84,5 +84,4 @@
       REAL(wp) ::   zdhnm, zhnnew, zhisn, zihic, zzc       !
       REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment
-      REAL(wp) ::   zds          ! increment of bottom ice salinity
       REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads
       REAL(wp) ::   zsm_snowice  ! snow-ice salinity
@@ -429,5 +428,5 @@
                   zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   &
                      &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) 
-                  zds         = zfracs * sss_m(ii,ij) - s_i_new(ji)
+                  zfracs = MIN( 0.5 , zfracs )
                   s_i_new(ji) = zfracs * sss_m(ii,ij)
                ENDIF ! fc_bo_i
@@ -438,5 +437,5 @@
          DO ji = kideb, kiut
             IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN
-               ! New ice salinity must not exceed 15 psu
+               ! New ice salinity must not exceed 20 psu
                s_i_new(ji) = MIN( s_i_new(ji), s_i_max )
                ! Metling point in K
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate.F90	(revision 4098)
+++ 	(revision )
@@ -1,1007 +1,0 @@
-MODULE limupdate
-   !!======================================================================
-   !!                     ***  MODULE  limupdate  ***
-   !!   LIM-3 : Update of sea-ice global variables at the end of the time step
-   !!======================================================================
-   !! History :  3.0  !  2006-04  (M. Vancoppenolle) Original code
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                      LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!    lim_update   : computes update of sea-ice global variables from trend terms
-   !!----------------------------------------------------------------------
-   USE dom_oce        ! ocean domain
-   USE oce            ! dynamics and tracers variables
-   USE sbc_oce        ! Surface boundary condition: ocean fields
-   USE sbc_ice        ! Surface boundary condition: ice fields
-   USE phycst         ! physical constants
-   USE ice            ! ice variables 
-   USE par_ice        ! ice parameters
-   USE thd_ice        ! ice thermodynamic variables
-   USE dom_ice        ! ice domain
-   USE limrhg         ! ice rheology
-   USE limdyn         ! ice dynamics
-   USE limtrp         ! ice transport
-   USE limthd         ! ice thermodynamics
-   USE limsbc         ! ice-oce surface boundary conditions
-   USE limdia         ! ice diagnostics
-   USE limwri         ! ice outputs
-   USE limrst         ! ice restart
-   USE limitd_th      ! ice thickness distribution (thermodynamics)
-   USE limvar         ! ice variables
-   USE prtctl         ! Print control
-   USE in_out_manager ! I/O manager
-   USE lbclnk         ! lateral boundary condition - MPP exchanges
-   USE wrk_nemo       ! work arrays
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   lim_update   ! routine called by ice_step
-
-      REAL(wp)  ::   epsi06 = 1.e-06_wp   ! module constants
-      REAL(wp)  ::   epsi04 = 1.e-04_wp   !    -       -
-      REAL(wp)  ::   epsi03 = 1.e-03_wp   !    -       -
-      REAL(wp)  ::   epsi10 = 1.e-10_wp   !    -       -
-      REAL(wp)  ::   epsi16 = 1.e-16_wp   !    -       -
-      REAL(wp)  ::   epsi20 = 1.e-20_wp   !    -       -
-      REAL(wp)  ::   rzero  = 0._wp       !    -       -
-      REAL(wp)  ::   rone   = 1._wp       !    -       -
-         
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE lim_update
-      !!-------------------------------------------------------------------
-      !!               ***  ROUTINE lim_update  ***
-      !!               
-      !! ** Purpose :  Computes update of sea-ice global variables at 
-      !!               the end of the time step.
-      !!               Address pathological cases
-      !!               This place is very important
-      !!                
-      !! ** Method  :  
-      !!    Ice speed from ice dynamics
-      !!    Ice thickness, Snow thickness, Temperatures, Lead fraction
-      !!      from advection and ice thermodynamics 
-      !!
-      !! ** Action  : - 
-      !!---------------------------------------------------------------------
-      INTEGER  ::   ji, jj, jk, jl, jm    ! dummy loop indices
-      INTEGER  ::   jbnd1, jbnd2
-      INTEGER  ::   i_ice_switch, ind_im, layer             ! local integers
-      REAL(wp) ::   zweight, zesum, zhimax                  ! local scalars
-      REAL(wp) ::   zindb, zindsn, zindic, zacrith          !   -      -
-      REAL(wp) ::   zrtt, zindg, zh, zdvres, zviold         !   -      - 
-      REAL(wp) ::   zbigvalue, zvsold                       !   -      -
-      REAL(wp) ::   z_prescr_hi, zat_i_old, ztmelts, ze_s   !   -      -
-      !
-      INTEGER , POINTER, DIMENSION(:,:,:) ::   internal_melt
-      REAL(wp), POINTER, DIMENSION(:)     ::   zthick0, zqm0   ! thickness of the layers and heat contents for
-      !!-------------------------------------------------------------------
-
-      CALL wrk_alloc( jpi,jpj,jpl, internal_melt )   ! integer
-      CALL wrk_alloc( jkmax, zthick0, zqm0 )
-
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' lim_update '
-         WRITE(numout,*) ' ~~~~~~~~~~ '
-
-         WRITE(numout,*) ' O) Initial values '
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-      ENDIF
-
-      !------------------------------------------------------------------------------
-      ! 1. Update of Global variables                                               |
-      !------------------------------------------------------------------------------
-
-      !---------------------
-      ! Ice dynamics  
-      !---------------------
-      u_ice(:,:) = u_ice(:,:) + d_u_ice_dyn(:,:)
-      v_ice(:,:) = v_ice(:,:) + d_v_ice_dyn(:,:)
-
-      !-----------------------------
-      ! Update ice and snow volumes  
-      !-----------------------------
-      DO jl = 1, jpl
-         v_i(:,:,jl)  = v_i(:,:,jl) + d_v_i_trp(:,:,jl) + d_v_i_thd(:,:,jl) 
-         v_s(:,:,jl)  = v_s(:,:,jl) + d_v_s_trp(:,:,jl) + d_v_s_thd(:,:,jl)
-      END DO
-
-      !---------------------------------
-      ! Classify the pathological cases
-      !---------------------------------
-      ! (1) v_i (new) > 0; d_v_i_thd + v_i(old) > 0 (easy case)
-      ! (2) v_i (new) > 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation)
-      ! (3) v_i (new) < 0; d_v_i_thd + v_i(old) > 0 (combined total ablation)
-      ! (4) v_i (new) < 0; d_v_i_thd + v_i(old) = 0 (total thermodynamic ablation 
-      ! with negative advection, very pathological )
-      ! (5) v_i (old) = 0; d_v_i_trp > 0 (advection of ice in a free-cell)
-      !
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               patho_case(ji,jj,jl) = 1
-               IF( v_i(ji,jj,jl) .GE. 0.0 ) THEN
-                  IF ( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN 
-                     patho_case(ji,jj,jl) = 2
-                  ENDIF
-               ELSE
-                  patho_case(ji,jj,jl) = 3
-                  IF( old_v_i(ji,jj,jl) + d_v_i_thd(ji,jj,jl) .LT. epsi10 ) THEN 
-                     patho_case(ji,jj,jl) = 4
-                  ENDIF
-               ENDIF
-               IF( ( old_v_i(ji,jj,jl) .LE. epsi10 ) .AND. &
-                   ( d_v_i_trp(ji,jj,jl) .GT. epsi06 ) ) THEN
-                  patho_case(ji,jj,jl) = 5 ! advection of ice in an ice-free
-                  ! cell
-                  IF( ln_nicep ) THEN  
-                     WRITE(numout,*) ' ALERTE patho_case still equal to 5 '
-                     WRITE(numout,*) ' ji , jj   : ', ji, jj
-                     WRITE(numout,*) ' old_v_i   : ', old_v_i(ji,jj,jl)
-                     WRITE(numout,*) ' d_v_i_trp : ', d_v_i_trp(ji,jj,jl)
-                  ENDIF
-
-               ENDIF
-            END DO
-         END DO
-      END DO
-
-      !--------------------
-      ! Excessive ablation 
-      !--------------------
-
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF (      ( patho_case(ji,jj,jl) .EQ. 3 ) &
-                  .OR. ( patho_case(ji,jj,jl) .EQ. 4 ) ) THEN
-                  zviold         = old_v_i(ji,jj,jl)
-                  zvsold         = old_v_s(ji,jj,jl)
-                  ! in cases 3 ( combined total ablation )
-                  !      and 4 ( total ablation with negative advection )
-                  ! there is excessive total ablation
-                  ! advection is chosen to be prioritary in order to conserve mass. 
-                  ! dv_i_thd is computed as a residual
-                  ! negative energy has to be kept in memory and to be given to the ocean
-                  ! equivalent salt flux is given to the ocean
-                  !
-                  ! This was the best solution found. Otherwise, mass conservation in advection
-                  ! scheme should have been revised, which could have been a big problem
-                  ! Martin Vancoppenolle (2006, updated 2007)
-
-                  ! is there any ice left ?
-                  zindic = MAX( rzero , SIGN( rone, v_i(ji,jj,jl) - epsi10 ) ) 
-                  !=1 if hi > 1e-3 and 0 if not
-                  zdvres = MAX( 0.0 , -v_i(ji,jj,jl) ) !residual volume if too much ice was molten
-                  !this quantity is positive
-                  v_i(ji,jj,jl) = zindic * v_i(ji,jj,jl)    !ice volume cannot be negative
-                  !correct thermodynamic ablation
-                  d_v_i_thd(ji,jj,jl) = zindic  *  d_v_i_thd(ji,jj,jl) + (1.0-zindic) * (-zviold - d_v_i_trp(ji,jj,jl)) 
-                  ! THIS IS NEW
-                  d_a_i_thd(ji,jj,jl) = zindic  *  d_a_i_thd(ji,jj,jl) + (1.0-zindic) * (-old_a_i(ji,jj,jl)) 
-
-                  !residual salt flux if ice is over-molten
-                  sfx_res(ji,jj)  = sfx_res(ji,jj) - sm_i(ji,jj,jl) * ( rhoic * zdvres * r1_rdtice )
-                  !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhoic * lfus * zdvres * r1_rdtice
-
-                  ! is there any snow left ?
-                  zindsn = MAX(  rzero, SIGN( rone , v_s(ji,jj,jl) - epsi10 )  ) 
-                  zvsold = v_s(ji,jj,jl)
-                  zdvres = MAX( 0.0 , -v_s(ji,jj,jl) )   !residual volume if too much ice was molten
-                  !this quantity is positive
-                  v_s(ji,jj,jl) = zindsn*v_s(ji,jj,jl)    !snow volume cannot be negative
-                  !correct thermodynamic ablation
-                  d_v_s_thd(ji,jj,jl)  = zindsn  *  d_v_s_thd(ji,jj,jl) + & 
-                     (1.0-zindsn) * (-zvsold - d_v_s_trp(ji,jj,jl)) 
-
-                  !no salt flux when snow is over-molten
-                  !             fheat_res(ji,jj)  = fheat_res(ji,jj) + rhosn * lfus * zdvres * r1_rdtice
-               ENDIF
-            END DO !ji
-         END DO !jj
-      END DO !jl
-
-      IF( ln_nicep ) THEN  
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF(  ABS( sfx_res(ji,jj) )  >  1._wp  ) THEN 
-                  WRITE(numout,*) ' ALERTE 1000 : residual salt flux of -> ', sfx_res(ji,jj)
-                  WRITE(numout,*) ' ji, jj : ', ji, jj, ' gphit, glamt : ', gphit(ji,jj), glamt(ji,jj)
-               ENDIF
-            END DO
-         END DO
-
-         WRITE(numout,*) ' 1. Before update of Global variables '
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-      ENDIF
-
-      !---------------------------------------------
-      ! Ice concentration and ice heat content
-      !---------------------------------------------
-
-      a_i (:,:,:) = a_i (:,:,:) + d_a_i_trp(:,:,:) + d_a_i_thd(:,:,:)
-      CALL lim_var_glo2eqv    ! useless, just for debug
-      IF( ln_nicep ) THEN 
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-      ENDIF
-      e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_trp(:,:,:,:)  
-      CALL lim_var_glo2eqv    ! useless, just for debug
-      IF( ln_nicep) THEN
-         WRITE(numout,*) ' After transport update '
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-      ENDIF
-      e_i(:,:,:,:) = e_i(:,:,:,:) + d_e_i_thd(:,:,:,:)  
-      CALL lim_var_glo2eqv ! useless, just for debug
-      IF( ln_nicep ) THEN
-         WRITE(numout,*) ' After thermodyn update '
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-      ENDIF
-
-      at_i(:,:) = a_i(:,:,1)                   ! total ice concentration
-      DO jl = 2, jpl
-         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-      END DO
-
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' 1. After update of Global variables (2) '
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' oa_i : ', oa_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' e_s : ', e_s(jiindx, jjindx, 1, 1:jpl)
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-      ENDIF
-
-      !------------------------------
-      ! Snow temperature and ice age
-      !------------------------------
-      e_s (:,:,:,:) = e_s (:,:,:,:) + d_e_s_trp (:,:,:,:) + d_e_s_thd (:,:,:,:)
-      oa_i(:,:,:)   = oa_i(:,:,:)   + d_oa_i_trp(:,:,:)   + d_oa_i_thd(:,:,:)
-
-      !--------------
-      ! Ice salinity    
-      !--------------
-
-      IF(  num_sal == 2  ) THEN      ! Prognostic salinity [Sice=F(z,t)]
-         !
-         IF( ln_nicep ) THEN  
-            WRITE(numout,*) ' Before everything '
-            WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-            WRITE(numout,*) ' oa_i:  ', oa_i(jiindx, jjindx, 1:jpl)
-            DO jk = 1, nlay_i
-               WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
-            END DO
-            WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-         ENDIF
-         !
-         smv_i(:,:,:) = smv_i(:,:,:) + d_smv_i_thd(:,:,:) + d_smv_i_trp(:,:,:)
-         !
-         IF( ln_nicep ) THEN  
-            WRITE(numout,*) ' After advection   '
-            WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-            WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-         ENDIF
-         !
-      ENDIF
-
-      CALL lim_var_glo2eqv
-
-      !--------------------------------------
-      ! 2. Review of all pathological cases
-      !--------------------------------------
-      zrtt    = 173.15_wp * rone
-      zacrith = 1.e-6_wp
-
-      !-------------------------------------------
-      ! 2.1) Advection of ice in an ice-free cell
-      !-------------------------------------------
-      ! should be removed since it is treated after dynamics now
-
-      zhimax = 5._wp
-      ! first category
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            !--- the thickness of such an ice is often out of bounds
-            !--- thus we recompute a new area while conserving ice volume
-            zat_i_old = SUM(old_a_i(ji,jj,:))
-            zindb          =  MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,1)) - epsi10 ) ) 
-            IF (      ( ABS(d_v_i_trp(ji,jj,1))/MAX(ABS(d_a_i_trp(ji,jj,1)),epsi10)*zindb.GT.zhimax) &
-               .AND.( ( v_i(ji,jj,1)/MAX(a_i(ji,jj,1),epsi10)*zindb).GT.zhimax ) &
-               .AND.( zat_i_old.LT.zacrith ) )  THEN ! new line
-               z_prescr_hi  = hi_max(1) * 0.5_wp
-               a_i(ji,jj,1) = v_i(ji,jj,1) / z_prescr_hi
-            ENDIF
-         END DO
-      END DO
-
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' 2.1 '
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-      ENDIF
-
-      !change this 14h44
-      zhimax = 20.0     ! line added up
-      ! change this also 17 aug
-      zhimax = 30.0     ! line added up
-
-      DO jl = 2, jpl
-         jm = ice_types(jl)
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zindb          =  MAX( rzero, SIGN( rone, ABS(d_a_i_trp(ji,jj,jl)) - epsi10 ) ) 
-               ! this correction is very tricky... sometimes, advection gets wrong i don't know why
-               ! it makes problems when the advected volume and concentration do not seem to be 
-               ! related with each other
-               ! the new thickness is sometimes very big!
-               ! and sometimes d_a_i_trp and d_v_i_trp have different sign
-               ! which of course is plausible
-               ! but fuck! it fucks everything up :)
-               IF ( (ABS(d_v_i_trp(ji,jj,jl))/MAX(ABS(d_a_i_trp(ji,jj,jl)),epsi10)*zindb.GT.zhimax) &
-                  .AND.(v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi10)*zindb).GT.zhimax ) THEN
-                  z_prescr_hi  =  ( hi_max_typ(jl-ice_cat_bounds(jm,1)  ,jm) + &
-                     hi_max_typ(jl-ice_cat_bounds(jm,1)+1,jm) ) / 2.0
-                  a_i (ji,jj,jl) = v_i(ji,jj,jl) / z_prescr_hi
-                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
-               ENDIF
-               zat_i_old = SUM(old_a_i(ji,jj,:))
-
-            END DO ! ji
-         END DO !jj
-      END DO !jl
-
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' 2.1 initial '
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-      ENDIF
-
-      at_i(:,:) = a_i(:,:,1)                   ! total ice concentration
-      DO jl = 2, jpl
-         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-      END DO
-
-      !----------------------------------------------------
-      ! 2.2) Rebin categories with thickness out of bounds
-      !----------------------------------------------------
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' 2.1 before rebinning '
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-      ENDIF
-
-      DO jm = 1, jpm
-         jbnd1 = ice_cat_bounds(jm,1)
-         jbnd2 = ice_cat_bounds(jm,2)
-         IF (ice_ncat_types(jm) .GT. 1 )   CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
-      END DO
-
-
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' 2.1 after rebinning'
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
-            WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-      ENDIF
-
-      at_i(:,:) = a_i(:,:,1)                   ! total ice concentration
-      DO jl = 2, jpl
-         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-      END DO
-
-      !---------------------------------
-      ! 2.3) Melt of an internal layer
-      !---------------------------------
-      internal_melt(:,:,:) = 0
-
-      DO jl = 1, jpl
-         DO jk = 1, nlay_i
-            DO jj = 1, jpj 
-               DO ji = 1, jpi
-                  ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
-                  IF ( ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR. &
-                     ( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) .AND. &
-                     ( v_i(ji,jj,jl) .GT. 0.0 ) .AND. &
-                     ( a_i(ji,jj,jl) .GT. 0.0 ) ) THEN
-                     !                    WRITE(numout,*) ' Internal layer melt : '
-                     !                    WRITE(numout,*) ' ji, jj, jk, jl : ', ji,jj,jk,jl
-                     !                    WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)
-                     !                    WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl)
-                     internal_melt(ji,jj,jl) = 1
-                  ENDIF
-               END DO ! ji
-            END DO ! jj
-         END DO !jk
-      END DO !jl
-
-      DO jl = 1, jpl
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-               IF( internal_melt(ji,jj,jl) == 1 ) THEN
-                  ! initial ice thickness
-                  !-----------------------
-                  ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
-                  !             WRITE(numout,*) ' ji,jj,jl : ', ji,jj,jl
-                  !             WRITE(numout,*) ' old ht_i: ', ht_i(ji,jj,jl)
-                  !             WRITE(numout,*) ' Enthalpy at the beg : ', e_i(ji,jj,1:nlay_i,jl)
-                  !             WRITE(numout,*) ' smv_i   : ', smv_i(ji,jj,jl)
-
-                  ! reduce ice thickness
-                  !-----------------------
-                  ind_im = 0
-                  zesum = 0.0
-                  DO jk = 1, nlay_i
-                     ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
-                     IF ( ( e_i(ji,jj,jk,jl) .LE. 0.0 ) .OR.  & 
-                        ( t_i(ji,jj,jk,jl) .GE. ztmelts ) ) &
-                        ind_im = ind_im + 1
-                     zesum = zesum + e_i(ji,jj,jk,jl)
-                  END DO
-                  IF (ind_im .LT.nlay_i ) smv_i(ji,jj,jl)= smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * & 
-                     ( ht_i(ji,jj,jl) - ind_im*ht_i(ji,jj,jl) / nlay_i )
-                  ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - ind_im*ht_i(ji,jj,jl) / nlay_i
-                  v_i(ji,jj,jl)  = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
-
-                  !             WRITE(numout,*) ' ind_im  : ', ind_im
-                  !             WRITE(numout,*) ' new ht_i: ', ht_i(ji,jj,jl)
-                  !             WRITE(numout,*) ' smv_i   : ', smv_i(ji,jj,jl)
-                  !             WRITE(numout,*) ' zesum   : ', zesum
-
-                  ! redistribute heat
-                  !-----------------------
-                  ! old thicknesses and enthalpies
-                  ind_im = 0
-                  DO jk = 1, nlay_i
-                     ztmelts = - tmut * s_i(ji,jj,jk,jl) + rtt
-                     IF ( ( e_i(ji,jj,jk,jl) .GT. 0.0 ) .AND.  & 
-                        ( t_i(ji,jj,jk,jl) .LT. ztmelts ) ) THEN
-                        ind_im = ind_im + 1
-                        zthick0(ind_im) = ht_i(ji,jj,jl) * ind_im / nlay_i
-                        zqm0   (ind_im) = MAX( e_i(ji,jj,jk,jl) , 0.0 )
-                     ENDIF
-                  END DO
-
-                  !             WRITE(numout,*) ' Old thickness, enthalpy '
-                  !             WRITE(numout,*) ' Number of layer : ind_im ', ind_im
-                  !             WRITE(numout,*) ' zthick0 : ', zthick0(1:ind_im)
-                  !             WRITE(numout,*) ' zqm0    : ', zqm0(1:ind_im)
-
-                  ! Redistributing energy on the new grid
-                  IF ( ind_im .GT. 0 ) THEN
-
-                     DO jk = 1, nlay_i
-                        e_i(ji,jj,jk,jl) = 0.0
-                        DO layer = 1, ind_im
-                           zweight         = MAX (  &
-                              MIN( ht_i(ji,jj,jl) * layer / ind_im , ht_i(ji,jj,jl) * jk / nlay_i ) -       &
-                              MAX( ht_i(ji,jj,jl) * (layer-1) / ind_im , ht_i(ji,jj,jl) * (jk-1) / nlay_i ) , 0.0 ) &
-                              /  ( ht_i(ji,jj,jl) / ind_im )
-
-                           e_i(ji,jj,jk,jl) =  e_i(ji,jj,jk,jl) + zweight*zqm0(layer)
-                        END DO !layer
-                     END DO ! jk
-
-                     zesum = 0.0
-                     DO jk = 1, nlay_i
-                        zesum = zesum + e_i(ji,jj,jk,jl)
-                     END DO
-
-                     !             WRITE(numout,*) ' Enthalpy at the end : ', e_i(ji,jj,1:nlay_i,jl)
-                     !             WRITE(numout,*) ' Volume   at the end : ', v_i(ji,jj,jl)
-                     !             WRITE(numout,*) ' zesum : ', zesum
-
-                  ELSE ! ind_im .EQ. 0, total melt
-                     e_i(ji,jj,jk,jl) = 0.0
-                  ENDIF
-
-               ENDIF ! internal_melt
-
-            END DO ! ji
-         END DO !jj
-      END DO !jl
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' 2.3 after melt of an internal ice layer '
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
-            WRITE(numout,*) ' t_i : ', t_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-      ENDIF
-
-      internal_melt(:,:,:) = 0
-
-      ! Melt of snow
-      !--------------
-      DO jl = 1, jpl
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-               ! snow energy of melting
-               ze_s = e_s(ji,jj,1,jl) * unit_fac / area(ji,jj) /              &
-                  MAX( v_s(ji,jj,jl), 1.0e-6 )  ! snow energy of melting
-
-               ! If snow energy of melting smaller then Lf
-               ! Then all snow melts and meltwater, heat go to the ocean
-               IF ( ze_s .LE. rhosn * lfus ) internal_melt(ji,jj,jl) = 1
-
-               IF( ln_nicep ) THEN  
-                  IF ( (ji.eq.jiindx) .AND. (jj.eq.jjindx) ) THEN
-                     WRITE(numout,*) ' jl    : ', jl
-                     WRITE(numout,*) ' ze_s  : ', ze_s
-                     WRITE(numout,*) ' v_s   : ', v_s(ji,jj,jl)
-                     WRITE(numout,*) ' rhosn : ', rhosn
-                     WRITE(numout,*) ' rhosn : ', lfus 
-                     WRITE(numout,*) ' area  : ', area(ji,jj)
-                     WRITE(numout,*) ' rhosn * lfus : ', rhosn * lfus 
-                     WRITE(numout,*) ' internal_melt : ', internal_melt(ji,jj,jl)
-                  ENDIF
-               ENDIF
-
-            END DO
-         END DO
-      END DO
-
-      DO jl = 1, jpl
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-               IF ( internal_melt(ji,jj,jl) == 1 ) THEN
-                  v_s(ji,jj,jl)   = 0._wp
-                  e_s(ji,jj,1,jl) = 0._wp
-                  !   ! release heat
-                  fheat_res(ji,jj) = fheat_res(ji,jj) + ze_s  * v_s(ji,jj,jl) * r1_rdtice
-                  ! release mass
-                  rdm_snw  (ji,jj) = rdm_snw  (ji,jj) + rhosn * v_s(ji,jj,jl)
-               ENDIF
-            END DO
-         END DO
-      END DO
-
-      zbigvalue      = 1.0d+20
-
-      DO jl = 1, jpl
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-
-               !switches
-               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
-               !switch = 1 if a_i > 1e-06 and 0 if not
-               zindsn        = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi06 ) ) !=1 if hs > 1e-6 and 0 if not
-               zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi04 ) ) !=1 if hi > 1e-3 and 0 if not
-               ! bug fix 25 avril 2007
-               zindb         = zindb*zindic
-
-               !--- 2.3 Correction to ice age 
-               !------------------------------
-               !                IF ((o_i(ji,jj,jl)-1.0)*86400.0.gt.(rdt_ice*float(numit))) THEN
-               !                   o_i(ji,jj,jl) = rdt_ice*FLOAT(numit)/86400.0
-               !                ENDIF
-               IF ((oa_i(ji,jj,jl)-1.0)*86400.0.gt.(rdt_ice*numit*a_i(ji,jj,jl))) THEN
-                  oa_i(ji,jj,jl) = rdt_ice * numit / 86400.0 * a_i(ji,jj,jl)
-               ENDIF
-               oa_i(ji,jj,jl) = zindb*zindic*oa_i(ji,jj,jl)
-
-               !--- 2.4 Correction to snow thickness
-               !-------------------------------------
-               !          ! snow thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hs = 0
-               !             v_s(ji,jj,jl)  = MAX( zindb * v_s(ji,jj,jl), 0.0) 
-               ! snow thickness cannot be smaller than 1e-6
-               v_s(ji,jj,jl)  = zindsn * v_s(ji,jj,jl) * zindb
-               v_s(ji,jj,jl)  = v_s(ji,jj,jl) *  MAX( 0.0 , SIGN( 1.0 , v_s(ji,jj,jl) - epsi06 ) )
-
-               !--- 2.5 Correction to ice thickness
-               !-------------------------------------
-               ! ice thickness has to be greater than 0, and if ice concentration smaller than 1e-6 then hi = 0
-               v_i(ji,jj,jl) = MAX( zindb * v_i(ji,jj,jl), 0.0)
-               ! ice thickness cannot be smaller than 1e-3
-               v_i(ji,jj,jl)  = zindic*v_i(ji,jj,jl)
-
-               !--- 2.6 Snow is transformed into ice if the original ice cover disappears
-               !----------------------------------------------------------------------------
-               zindg          = tms(ji,jj) *  MAX( rzero , SIGN( rone , -v_i(ji,jj,jl) ) )
-               v_i(ji,jj,jl)  = v_i(ji,jj,jl) + zindg * rhosn * v_s(ji,jj,jl) / rau0
-               v_s(ji,jj,jl)  = ( rone - zindg ) * v_s(ji,jj,jl) + & 
-                  zindg * v_i(ji,jj,jl) * ( rau0 - rhoic ) / rhosn
-
-               !--- 2.7 Correction to ice concentrations 
-               !--------------------------------------------
-               ! if greater than 0, ice concentration cannot be smaller than 1e-10
-               a_i(ji,jj,jl) = zindb * MAX(zindsn, zindic) * MAX( a_i(ji,jj,jl), epsi06 )
-               ! then ice volume has to be corrected too...
-               ! instead, zap small areas
-
-               !-------------------------
-               ! 2.8) Snow heat content
-               !-------------------------
-
-               e_s(ji,jj,1,jl) = zindsn *                                &
-                  ( MIN ( MAX ( 0.0, e_s(ji,jj,1,jl) ), zbigvalue ) ) + &
-                  ( 1.0 - zindsn ) * 0.0
-
-            END DO ! ji
-         END DO ! jj
-      END DO ! jl
-
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' 2.8 '
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i: ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-      ENDIF
-
-      !------------------------
-      ! 2.9) Ice heat content 
-      !------------------------
-
-      DO jl = 1, jpl
-         DO jk = 1, nlay_i
-            DO jj = 1, jpj 
-               DO ji = 1, jpi
-                  zindic        = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi06 ) ) 
-                  ! =1 if v_i > 1e-6 and 0 if not
-                  e_i(ji,jj,jk,jl)= zindic * & 
-                     ( MIN ( MAX ( 0.0, e_i(ji,jj,jk,jl) ), zbigvalue ) ) + &
-                     ( 1.0 - zindic ) * 0.0
-               END DO ! ji
-            END DO ! jj
-         END DO !jk
-      END DO !jl
-
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' 2.9 '
-         DO jk = 1, nlay_i
-            WRITE(numout,*) ' e_i : ', e_i(jiindx, jjindx, jk, 1:jpl)
-         END DO
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-      ENDIF
-
-      !---------------------
-      ! 2.11) Ice salinity
-      !---------------------
-
-      IF(  num_sal == 2  ) THEN      ! Prognostic salinity [Sice=F(z,t)]
-         !
-         DO jl = 1, jpl
-            DO jk = 1, nlay_i
-               DO jj = 1, jpj 
-                  DO ji = 1, jpi           ! salinity stays in bounds
-                     smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)), 0.1 * v_i(ji,jj,jl) )
-                     i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl)))
-                     smv_i(ji,jj,jl)  = i_ice_switch*smv_i(ji,jj,jl) + 0.1*(1.0-i_ice_switch)*v_i(ji,jj,jl)
-                  END DO ! ji
-               END DO ! jj
-            END DO !jk
-         END DO !jl
-         !
-      ENDIF
-
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' 2.11 '
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i    ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-      ENDIF
-
-      DO jm = 1, jpm
-         DO jj = 1, jpj 
-            DO ji = 1, jpi
-               jl = ice_cat_bounds(jm,1)
-               !--- 2.12 Constrain the thickness of the smallest category above 5 cm
-               !----------------------------------------------------------------------
-               ! the ice thickness of the smallest category should be higher than 5 cm
-               ! we changed hiclim to 10
-               zindb         = MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
-               ht_i(ji,jj,jl) = zindb*v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl), epsi06)
-               zh            = MAX( rone , zindb * hiclim  / MAX( ht_i(ji,jj,jl) , epsi20 ) )
-               ht_s(ji,jj,jl) = ht_s(ji,jj,jl)* zh
-               !             v_s(ji,jj,jl)  = v_s(ji,jj,jl) * zh
-               ht_i(ji,jj,jl) = ht_i(ji,jj,jl)* zh
-               !             v_i(ji,jj,jl)  = v_i(ji,jj,jl) * zh
-               a_i (ji,jj,jl) = a_i(ji,jj,jl) /zh
-            END DO !ji
-         END DO !jj
-      END DO !jm
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' 2.12 '
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i    ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-      ENDIF
-
-      ! Total ice concentration check : forced to be less than amax
-      !------------------------------
-      at_i(:,:) = a_i(:,:,1)                  ! total ice concentration
-      DO jl = 2, jpl
-         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-      END DO
-      DO jl = 1, jpl                          ! if at_i > amax, set at_i=amax by a rescale of the area of each category
-         a_i(:,:,jl) = a_i(:,:,jl) * amax / MAX( at_i(:,:), amax )
-      END DO
-      at_i(:,:) = MIN( at_i(:,:), amax )      ! at_i becomes less than amax by construction no need to recompute it
-
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' 2.13 '
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl), ' at_i    ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl), ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-      ENDIF
-
-
-      IF( ln_nicep ) THEN  
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               IF (at_i(ji,jj).GT.1.0) THEN
-                  WRITE(numout,*) ' lim_update ! : at_i > 1 -> PAS BIEN -> ALERTE '
-                  WRITE(numout,*) ' ~~~~~~~~~~   at_i     ', at_i(ji,jj)
-                  WRITE(numout,*) ' Point ', ji, jj
-                  WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
-                  DO jl = 1, jpl
-                     WRITE(numout,*) ' a_i ***         ', a_i(ji,jj,jl), ' CAT no ', jl
-                     WRITE(numout,*) ' a_i_old ***     ', old_a_i(ji,jj,jl), ' CAT no ', jl
-                     WRITE(numout,*) ' d_a_i_thd / trp ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
-                  END DO
-               ENDIF
-            END DO
-         END DO
-      ENDIF
-
-      ! Final thickness distribution rebinning
-      ! --------------------------------------
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' rebinning before'
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i    ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-      ENDIF
-      !old version
-      !    CALL lim_itd_th_reb(1,jpl)
-
-      DO jm = 1, jpm
-         jbnd1 = ice_cat_bounds(jm,1)
-         jbnd2 = ice_cat_bounds(jm,2)
-         IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_reb(jbnd1, jbnd2, jm)
-         IF (ice_ncat_types(jm) .EQ. 1 ) THEN
-         ENDIF
-      END DO
-
-      IF( ln_nicep ) THEN  
-         WRITE(numout,*) ' rebinning final'
-         WRITE(numout,*) ' a_i : ', a_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' at_i    ', at_i(jiindx,jjindx)
-         WRITE(numout,*) ' v_i : ', v_i(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' v_s : ', v_s(jiindx, jjindx, 1:jpl)
-         WRITE(numout,*) ' smv_i: ', smv_i(jiindx, jjindx, 1:jpl)
-      ENDIF
-
-      at_i(:,:) = a_i(:,:,1)                   ! total ice concentration
-      DO jl = 2, jpl
-         at_i(:,:) = at_i(:,:) + a_i(:,:,jl)
-      END DO
-
-      !------------------------------------------------------------------------------
-      ! 2) Corrections to avoid wrong values                                        |
-      !------------------------------------------------------------------------------
-      ! Ice drift
-      !------------
-!!gm  BUG ? I don't understand this : it may have a wrong impact on the ice edge advection
-!!gm  and any way there is much faster way to code that...
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            IF ( at_i(ji,jj) .EQ. 0.0 ) THEN ! what to do if there is no ice
-               IF ( at_i(ji+1,jj) .EQ. 0.0 ) u_ice(ji,jj)   = 0.0 ! right side
-               IF ( at_i(ji-1,jj) .EQ. 0.0 ) u_ice(ji-1,jj) = 0.0 ! left side
-               IF ( at_i(ji,jj+1) .EQ. 0.0 ) v_ice(ji,jj)   = 0.0 ! upper side
-               IF ( at_i(ji,jj-1) .EQ. 0.0 ) v_ice(ji,jj-1) = 0.0 ! bottom side
-            ENDIF
-         END DO
-      END DO
-      !mask velocities
-!!gm BUG ?  here the mask are the one of the beginning of the time step, no? 
-!!gm        whereas at this level they should have been updated... To be checked 
-      u_ice(:,:) = u_ice(:,:) * tmu(:,:)
-      v_ice(:,:) = v_ice(:,:) * tmv(:,:)
-      !lateral boundary conditions
-      CALL lbc_lnk( u_ice(:,:), 'U', -1. )
-      CALL lbc_lnk( v_ice(:,:), 'V', -1. )
-
-      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-      ! ALERTES
-      !++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
-
-      IF( ln_nicep ) THEN  
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               DO jl = 1, jpl
-                  IF ( (a_i(ji,jj,jl).GT.1.0).OR.(at_i(ji,jj).GT.1.0) ) THEN
-                     zindb          =  MAX( rzero, SIGN( rone, a_i(ji,jj,jl) - epsi06 ) ) 
-                     WRITE(numout,*) ' lim_update : a_i > 1 '
-                     WRITE(numout,*) ' PAS BIEN ----> ALERTE !!! '
-                     WRITE(numout,*) ' ~~~~~~~~~~   at_i     ', at_i(ji,jj)
-                     WRITE(numout,*) ' Point - category', ji, jj, jl
-                     WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
-                     WRITE(numout,*) ' a_i *** a_i_old ', a_i(ji,jj,jl), old_a_i(ji,jj,jl)
-                     WRITE(numout,*) ' v_i *** v_i_old ', v_i(ji,jj,jl), old_v_i(ji,jj,jl)
-                     WRITE(numout,*) ' ht_i ***        ', v_i(ji,jj,jl)/MAX(a_i(ji,jj,jl),epsi06)*zindb 
-                     WRITE(numout,*) ' hi_max(jl), hi_max(jl-1) ', hi_max(jl), hi_max(jl-1)
-                     WRITE(numout,*) ' d_v_i_thd / trp ', d_v_i_thd(ji,jj,jl), d_v_i_trp(ji,jj,jl)
-                     WRITE(numout,*) ' d_a_i_thd / trp ', d_a_i_thd(ji,jj,jl), d_a_i_trp(ji,jj,jl)
-                  ENDIF
-               END DO
-
-            END DO !jj
-         END DO !ji
-
-         WRITE(numout,*) ' TESTOSC1 ', tio_u(jiindx,jjindx), tio_v(jiindx,jjindx)
-         WRITE(numout,*) ' TESTOSC2 ', u_ice(jiindx,jjindx), v_ice(jiindx,jjindx)
-         WRITE(numout,*) ' TESTOSC3 ', u_oce(jiindx,jjindx), v_oce(jiindx,jjindx)
-         WRITE(numout,*) ' TESTOSC4 ', utau (jiindx,jjindx), vtau (jiindx,jjindx)
-      ENDIF
-
-
-      IF(ln_ctl) THEN   ! Control print
-         CALL prt_ctl_info(' ')
-         CALL prt_ctl_info(' - Cell values : ')
-         CALL prt_ctl_info('   ~~~~~~~~~~~~~ ')
-         CALL prt_ctl(tab2d_1=area       , clinfo1=' lim_update  : cell area   :')
-         CALL prt_ctl(tab2d_1=at_i       , clinfo1=' lim_update  : at_i        :')
-         CALL prt_ctl(tab2d_1=vt_i       , clinfo1=' lim_update  : vt_i        :')
-         CALL prt_ctl(tab2d_1=vt_s       , clinfo1=' lim_update  : vt_s        :')
-         CALL prt_ctl(tab2d_1=strength   , clinfo1=' lim_update  : strength    :')
-         CALL prt_ctl(tab2d_1=u_ice      , clinfo1=' lim_update  : u_ice       :', tab2d_2=v_ice      , clinfo2=' v_ice       :')
-         CALL prt_ctl(tab2d_1=d_u_ice_dyn, clinfo1=' lim_update  : d_u_ice_dyn :', tab2d_2=d_v_ice_dyn, clinfo2=' d_v_ice_dyn :')
-         CALL prt_ctl(tab2d_1=old_u_ice  , clinfo1=' lim_update  : old_u_ice   :', tab2d_2=old_v_ice  , clinfo2=' old_v_ice   :')
-
-         DO jl = 1, jpl
-            CALL prt_ctl_info(' ')
-            CALL prt_ctl_info(' - Category : ', ivar1=jl)
-            CALL prt_ctl_info('   ~~~~~~~~~~')
-            CALL prt_ctl(tab2d_1=ht_i       (:,:,jl)        , clinfo1= ' lim_update  : ht_i        : ')
-            CALL prt_ctl(tab2d_1=ht_s       (:,:,jl)        , clinfo1= ' lim_update  : ht_s        : ')
-            CALL prt_ctl(tab2d_1=t_su       (:,:,jl)        , clinfo1= ' lim_update  : t_su        : ')
-            CALL prt_ctl(tab2d_1=t_s        (:,:,1,jl)      , clinfo1= ' lim_update  : t_snow      : ')
-            CALL prt_ctl(tab2d_1=sm_i       (:,:,jl)        , clinfo1= ' lim_update  : sm_i        : ')
-            CALL prt_ctl(tab2d_1=o_i        (:,:,jl)        , clinfo1= ' lim_update  : o_i         : ')
-            CALL prt_ctl(tab2d_1=a_i        (:,:,jl)        , clinfo1= ' lim_update  : a_i         : ')
-            CALL prt_ctl(tab2d_1=old_a_i    (:,:,jl)        , clinfo1= ' lim_update  : old_a_i     : ')
-            CALL prt_ctl(tab2d_1=d_a_i_trp  (:,:,jl)        , clinfo1= ' lim_update  : d_a_i_trp   : ')
-            CALL prt_ctl(tab2d_1=d_a_i_thd  (:,:,jl)        , clinfo1= ' lim_update  : d_a_i_thd   : ')
-            CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' lim_update  : v_i         : ')
-            CALL prt_ctl(tab2d_1=old_v_i    (:,:,jl)        , clinfo1= ' lim_update  : old_v_i     : ')
-            CALL prt_ctl(tab2d_1=d_v_i_trp  (:,:,jl)        , clinfo1= ' lim_update  : d_v_i_trp   : ')
-            CALL prt_ctl(tab2d_1=d_v_i_thd  (:,:,jl)        , clinfo1= ' lim_update  : d_v_i_thd   : ')
-            CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' lim_update  : v_s         : ')
-            CALL prt_ctl(tab2d_1=old_v_s    (:,:,jl)        , clinfo1= ' lim_update  : old_v_s     : ')
-            CALL prt_ctl(tab2d_1=d_v_s_trp  (:,:,jl)        , clinfo1= ' lim_update  : d_v_s_trp   : ')
-            CALL prt_ctl(tab2d_1=d_v_s_thd  (:,:,jl)        , clinfo1= ' lim_update  : d_v_s_thd   : ')
-            CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)/1.0e9, clinfo1= ' lim_update  : e_i1        : ')
-            CALL prt_ctl(tab2d_1=old_e_i    (:,:,1,jl)/1.0e9, clinfo1= ' lim_update  : old_e_i1    : ')
-            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update  : de_i1_trp   : ')
-            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update  : de_i1_thd   : ')
-            CALL prt_ctl(tab2d_1=e_i        (:,:,2,jl)/1.0e9, clinfo1= ' lim_update  : e_i2        : ')
-            CALL prt_ctl(tab2d_1=old_e_i    (:,:,2,jl)/1.0e9, clinfo1= ' lim_update  : old_e_i2    : ')
-            CALL prt_ctl(tab2d_1=d_e_i_trp  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update  : de_i2_trp   : ')
-            CALL prt_ctl(tab2d_1=d_e_i_thd  (:,:,2,jl)/1.0e9, clinfo1= ' lim_update  : de_i2_thd   : ')
-            CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' lim_update  : e_snow      : ')
-            CALL prt_ctl(tab2d_1=old_e_s    (:,:,1,jl)      , clinfo1= ' lim_update  : old_e_snow  : ')
-            CALL prt_ctl(tab2d_1=d_e_s_trp  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update  : d_e_s_trp   : ')
-            CALL prt_ctl(tab2d_1=d_e_s_thd  (:,:,1,jl)/1.0e9, clinfo1= ' lim_update  : d_e_s_thd   : ')
-            CALL prt_ctl(tab2d_1=smv_i      (:,:,jl)        , clinfo1= ' lim_update  : smv_i       : ')
-            CALL prt_ctl(tab2d_1=old_smv_i  (:,:,jl)        , clinfo1= ' lim_update  : old_smv_i   : ')
-            CALL prt_ctl(tab2d_1=d_smv_i_trp(:,:,jl)        , clinfo1= ' lim_update  : d_smv_i_trp : ')
-            CALL prt_ctl(tab2d_1=d_smv_i_thd(:,:,jl)        , clinfo1= ' lim_update  : d_smv_i_thd : ')
-            CALL prt_ctl(tab2d_1=oa_i       (:,:,jl)        , clinfo1= ' lim_update  : oa_i        : ')
-            CALL prt_ctl(tab2d_1=old_oa_i   (:,:,jl)        , clinfo1= ' lim_update  : old_oa_i    : ')
-            CALL prt_ctl(tab2d_1=d_oa_i_trp (:,:,jl)        , clinfo1= ' lim_update  : d_oa_i_trp  : ')
-            CALL prt_ctl(tab2d_1=d_oa_i_thd (:,:,jl)        , clinfo1= ' lim_update  : d_oa_i_thd  : ')
-            CALL prt_ctl(tab2d_1=REAL(patho_case(:,:,jl))   , clinfo1= ' lim_update  : Path. case  : ')
-
-            DO jk = 1, nlay_i
-               CALL prt_ctl_info(' - Layer : ', ivar1=jk)
-               CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_update  : t_i       : ')
-            END DO
-         END DO
-
-         CALL prt_ctl_info(' ')
-         CALL prt_ctl_info(' - Heat / FW fluxes : ')
-         CALL prt_ctl_info('   ~~~~~~~~~~~~~~~~~~ ')
-         CALL prt_ctl(tab2d_1=fmmec  , clinfo1= ' lim_update : fmmec : ', tab2d_2=fhmec     , clinfo2= ' fhmec     : ')
-         CALL prt_ctl(tab2d_1=sst_m  , clinfo1= ' lim_update : sst   : ', tab2d_2=sss_m     , clinfo2= ' sss       : ')
-         CALL prt_ctl(tab2d_1=fhbri  , clinfo1= ' lim_update : fhbri : ', tab2d_2=fheat_mec , clinfo2= ' fheat_mec : ')
-
-         CALL prt_ctl_info(' ')
-         CALL prt_ctl_info(' - Stresses : ')
-         CALL prt_ctl_info('   ~~~~~~~~~~ ')
-         CALL prt_ctl(tab2d_1=utau       , clinfo1= ' lim_update : utau      : ', tab2d_2=vtau       , clinfo2= ' vtau      : ')
-         CALL prt_ctl(tab2d_1=utau_ice   , clinfo1= ' lim_update : utau_ice  : ', tab2d_2=vtau_ice   , clinfo2= ' vtau_ice  : ')
-         CALL prt_ctl(tab2d_1=u_oce      , clinfo1= ' lim_update : u_oce     : ', tab2d_2=v_oce      , clinfo2= ' v_oce     : ')
-      ENDIF
-
-      CALL wrk_dealloc( jpi,jpj,jpl, internal_melt )   ! integer
-      CALL wrk_dealloc( jkmax, zthick0, zqm0 )
-
-   END SUBROUTINE lim_update
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option         Empty Module               No sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_update     ! Empty routine
-   END SUBROUTINE lim_update
-
-#endif
-
-END MODULE limupdate
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90	(revision 4099)
@@ -26,5 +26,4 @@
    USE limthd
    USE limsbc
-   USE limdia
    USE limdiahsb
    USE limwri
@@ -402,7 +401,8 @@
       ! 2.11) Ice salinity
       !---------------------
-      !clem@bug: smv_i should be updated too: smv_i(:,:,:) = smv_i(:,:,:) + sm_i(:,:,:) * ( v_i(:,:,:) - zviold(:,:,:) )
+      ! clem correct bug on smv_i
+      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
+
       IF (  num_sal == 2  ) THEN ! general case
-
          DO jl = 1, jpl
             !DO jk = 1, nlay_i
@@ -410,12 +410,12 @@
                   DO ji = 1, jpi
                      ! salinity stays in bounds
-                     smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) )
-                     i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl) + epsi20))
-                     smv_i(ji,jj,jl)  = i_ice_switch*smv_i(ji,jj,jl) + 0.1*(1.0-i_ice_switch)*v_i(ji,jj,jl)
+                     !clem smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) )
+                     smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) )
+                     i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20 ) )
+                     smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) + s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl)
                   END DO ! ji
                END DO ! jj
             !END DO !jk
          END DO !jl
-
       ENDIF
 
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90	(revision 4099)
@@ -26,5 +26,4 @@
    USE limthd
    USE limsbc
-   USE limdia
    USE limdiahsb
    USE limwri
@@ -305,7 +304,7 @@
                      zesum = zesum + e_i(ji,jj,jk,jl)
                   END DO
-                  IF (ind_im .LT. nlay_i ) THEN
-                     smv_i(ji,jj,jl) = smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) )
-                  ENDIF
+                  !clem IF (ind_im .LT. nlay_i ) THEN
+                  !clem   smv_i(ji,jj,jl) = smv_i(ji,jj,jl) / ht_i(ji,jj,jl) * ( ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i) )
+                  !clem ENDIF
                   ht_i(ji,jj,jl) = ht_i(ji,jj,jl) - REAL(ind_im)*ht_i(ji,jj,jl) / REAL(nlay_i)
                   v_i(ji,jj,jl)  = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
@@ -532,5 +531,7 @@
       ! 2.11) Ice salinity
       !---------------------
-      !clem@bug: smv_i should be updated too: smv_i(:,:,:) = smv_i(:,:,:) + sm_i(:,:,:) * ( v_i(:,:,:) - zviold(:,:,:) )
+      ! clem correct bug on smv_i
+      smv_i(:,:,:) = sm_i(:,:,:) * v_i(:,:,:)
+
       IF (  num_sal == 2  ) THEN ! general case
          DO jl = 1, jpl
@@ -539,7 +540,8 @@
                   DO ji = 1, jpi
                      ! salinity stays in bounds
-                     smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) )
-                     i_ice_switch    =  1.0-MAX(0.0,SIGN(1.0,-v_i(ji,jj,jl) + epsi20))
-                     smv_i(ji,jj,jl)  = i_ice_switch*smv_i(ji,jj,jl) + 0.1*(1.0-i_ice_switch)*v_i(ji,jj,jl)
+                     !clem smv_i(ji,jj,jl)  =  MAX(MIN((rhoic-rhosn)/rhoic*sss_m(ji,jj),smv_i(ji,jj,jl)),0.1 * v_i(ji,jj,jl) )
+                     smv_i(ji,jj,jl) = MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) )
+                     i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, -v_i(ji,jj,jl) + epsi20 ) )
+                     smv_i(ji,jj,jl) = i_ice_switch * smv_i(ji,jj,jl) + s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl)
                   END DO ! ji
                END DO ! jj
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 4099)
@@ -30,4 +30,8 @@
 
    PUBLIC sbc_ice_alloc ! called in iceini(_2).F90
+
+   CHARACTER (len=8), PUBLIC :: cn_iceflx = 'none'                !: Flux handling over ice categories
+   LOGICAL, PUBLIC :: ln_iceflx_ave    = .FALSE. ! Average heat fluxes over all ice categories
+   LOGICAL, PUBLIC :: ln_iceflx_linear = .FALSE. ! Redistribute mean heat fluxes over all ice categories, using ice temperature and albedo
 
 # if defined  key_lim2
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 4099)
@@ -471,6 +471,6 @@
       ! Coupled case: since cloud cover is not received from atmosphere 
       !               ===> defined as constant value -> definition done in sbc_cpl_init
-      fr1_i0(:,:) = 0.18
-      fr2_i0(:,:) = 0.82
+      IF ( ALLOCATED (fr1_i0)) fr1_i0 (:,:) = 0.18
+      IF ( ALLOCATED (fr2_i0)) fr2_i0 (:,:) = 0.82
       !                                                      ! ------------------------- !
       !                                                      !      10m wind module      !   
@@ -931,5 +931,6 @@
       CALL wrk_alloc( jpi,jpj, ztx, zty )
 
-      IF( srcv(jpr_itx1)%laction ) THEN   ;   itx =  jpr_itx1   
+!AC Pour eviter un stress nul sur la glace dans le cas mixed oce-ice
+      IF( srcv(jpr_itx1)%laction .AND. TRIM( sn_rcv_tau%cldes ) == 'oce and ice') THEN   ;   itx =  jpr_itx1   
       ELSE                                ;   itx =  jpr_otx1
       ENDIF
@@ -938,7 +939,8 @@
       IF(  nrcvinfo(itx) == OASIS_Rcv ) THEN
 
-         !                                                      ! ======================= !
-         IF( srcv(jpr_itx1)%laction ) THEN                      !   ice stress received   !
-            !                                                   ! ======================= !
+         !                                                                                              ! ======================= !
+!AC Pour eviter un stress nul sur la glace dans le cas mixes oce-ice
+         IF( srcv(jpr_itx1)%laction .AND. TRIM( sn_rcv_tau%cldes ) == 'oce and ice') THEN               !   ice stress received   !
+            !                                                                                           ! ======================= !
             !  
             IF( TRIM( sn_rcv_tau%clvref ) == 'cartesian' ) THEN            ! 2 components on the sphere
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_if.F90	(revision 4099)
@@ -99,5 +99,7 @@
          
          fr_i(:,:) = tfreez( sss_m ) * tmask(:,:,1)      ! sea surface freezing temperature [Celcius]
-#if defined key_coupled 
+
+! OM : probleme. a_i pas defini dans les cas lim3 et cice
+#if defined key_coupled && defined key_lim2
          a_i(:,:,1) = fr_i(:,:)         
 #endif
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4099)
@@ -32,4 +32,5 @@
    USE sbcblk_core     ! Surface boundary condition: CORE bulk
    USE sbcblk_clio     ! Surface boundary condition: CLIO bulk
+   USE sbccpl          ! Surface boundary condition: coupled interface
    USE albedo          ! ocean & ice albedo
 
@@ -42,5 +43,4 @@
    USE limitd_me       ! Mechanics on ice thickness distribution
    USE limsbc          ! sea surface boundary condition
-   USE limdia          ! Ice diagnostics
    USE limdiahsb       ! Ice budget diagnostics
    USE limwri          ! Ice outputs
@@ -77,4 +77,35 @@
    !!----------------------------------------------------------------------
 CONTAINS
+
+   FUNCTION fice_cell_ave ( ptab)
+      !!--------------------------------------------------------------------------
+      !! * Compute average over categories, for grid cell (ice covered and free ocean)
+      !!--------------------------------------------------------------------------
+      REAL (wp), DIMENSION (jpi,jpj) :: fice_cell_ave
+      REAL (wp), DIMENSION (jpi,jpj,jpl), INTENT (in) :: ptab
+      INTEGER :: jl ! Dummy loop index
+      
+      fice_cell_ave (:,:) = 0.0_wp
+      
+      DO jl = 1, jpl
+         fice_cell_ave (:,:) = fice_cell_ave (:,:) &
+            &                  + a_i (:,:,jl) * ptab (:,:,jl)
+      END DO
+      
+   END FUNCTION fice_cell_ave
+   
+   FUNCTION fice_ice_ave ( ptab)
+      !!--------------------------------------------------------------------------
+      !! * Compute average over categories, for ice covered part of grid cell
+      !!--------------------------------------------------------------------------
+      REAL (kind=wp), DIMENSION (jpi,jpj) :: fice_ice_ave
+      REAL (kind=wp), DIMENSION (jpi,jpj,jpl), INTENT(in) :: ptab
+
+      fice_ice_ave (:,:) = 0.0_wp
+      WHERE ( at_i (:,:) .GT. 0.0_wp ) fice_ice_ave (:,:) = fice_cell_ave ( ptab (:,:,:)) / at_i (:,:)
+
+   END FUNCTION fice_ice_ave
+
+   !!======================================================================
 
    SUBROUTINE sbc_ice_lim( kt, kblk )
@@ -104,9 +135,29 @@
       REAL(wp) ::   zcoef   ! local scalar
       REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice_os, zalb_ice_cs  ! albedo of the ice under overcast/clear sky
+      REAL(wp), POINTER, DIMENSION(:,:,:)   ::   zalb_ice      ! mean albedo of ice (for coupled)
+
+      REAL(wp), POINTER, DIMENSION(:,:) :: zalb_ice_all    ! Mean albedo over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: ztem_ice_all    ! Mean temperature over all categories
+      
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_qsr_ice_all   ! Mean solar heat flux over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_qns_ice_all   ! Mean non solar heat flux over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_qla_ice_all   ! Mean latent heat flux over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqns_ice_all  ! Mean d(qns)/dT over all categories
+      REAL(wp), POINTER, DIMENSION(:,:) :: z_dqla_ice_all  ! Mean d(qla)/dT over all categories
       !!----------------------------------------------------------------------
 
+      !- O.M. : why do we allocate all these arrays even when MOD( kt-1, nn_fsbc ) /= 0 ?????
+
       IF( nn_timing == 1 )  CALL timing_start('sbc_ice_lim')
 
       CALL wrk_alloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs )
+
+      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+         CALL wrk_alloc( jpi,jpj,jpl, zalb_ice)
+      END IF
+      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+         CALL wrk_alloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
+      ENDIF
+
 
       IF( kt == nit000 ) THEN
@@ -139,4 +190,14 @@
             t_su(:,:,jl) = t_su(:,:,jl) +  rt0 * ( 1. - tmask(:,:,1) )
          END DO
+
+         IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) )
+         
+         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+            !
+            ! Compute mean albedo and temperature
+            zalb_ice_all (:,:) = fice_ice_ave ( zalb_ice (:,:,:) ) 
+            ztem_ice_all (:,:) = fice_ice_ave ( tn_ice   (:,:,:) ) 
+            !
+         ENDIF
                                                      ! Bulk formulea - provides the following fields:
          ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2]
@@ -161,5 +222,49 @@
                &                      tprecip   , sprecip   ,                            &
                &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  )
+            !
+         CASE ( 5 )
+            zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) )
+            
+            CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )
+
+            CALL sbc_cpl_ice_flx( p_frld=ato_i, palbi=zalb_ice, psst=sst_m, pist=tn_ice    )
+
+            ! Latent heat flux is forced to 0 in coupled :
+            !  it is included in qns (non-solar heat flux)
+            qla_ice  (:,:,:) = 0.0e0_wp
+            dqla_ice (:,:,:) = 0.0e0_wp
+            !
          END SELECT
+
+         ! Average over all categories
+         IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+
+            z_qns_ice_all  (:,:) = fice_ice_ave ( qns_ice  (:,:,:) )
+            z_qsr_ice_all  (:,:) = fice_ice_ave ( qsr_ice  (:,:,:) )
+            z_dqns_ice_all (:,:) = fice_ice_ave ( dqns_ice (:,:,:) )
+            z_qla_ice_all  (:,:) = fice_ice_ave ( qla_ice  (:,:,:) )
+            z_dqla_ice_all (:,:) = fice_ice_ave ( dqla_ice (:,:,:) )
+
+            DO jl = 1, jpl
+               dqns_ice (:,:,jl) = z_dqns_ice_all (:,:)
+               dqla_ice (:,:,jl) = z_dqla_ice_all (:,:)
+            END DO
+            !
+            IF ( ln_iceflx_ave ) THEN
+               DO jl = 1, jpl
+                  qns_ice  (:,:,jl) = z_qns_ice_all  (:,:)
+                  qsr_ice  (:,:,jl) = z_qsr_ice_all  (:,:)
+                  qla_ice  (:,:,jl) = z_qla_ice_all  (:,:)
+               END DO
+            END IF
+            !
+            IF ( ln_iceflx_linear ) THEN
+               DO jl = 1, jpl
+                  qns_ice  (:,:,jl) = z_qns_ice_all(:,:) + z_dqns_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:))
+                  qla_ice  (:,:,jl) = z_qla_ice_all(:,:) + z_dqla_ice_all(:,:) * (tn_ice(:,:,jl) - ztem_ice_all(:,:))
+                  qsr_ice  (:,:,jl) = (1.0e0_wp-zalb_ice(:,:,jl)) / (1.0e0_wp-zalb_ice_all(:,:)) * z_qsr_ice_all(:,:)
+               END DO
+            END IF
+         END IF
 
          !                                           !----------------------!
@@ -264,8 +369,8 @@
          !
          !                                           ! Diagnostics and outputs 
-         IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   &
-            &             CALL lim_dia 
-         IF (ln_limdiahsb) CALL lim_diahsb
+         IF (ln_limdiaout) CALL lim_diahsb
+!clem # if ! defined key_iomput
                           CALL lim_wri( 1  )              ! Ice outputs 
+!clem # endif
          IF( kt == nit000 )   CALL iom_close( numrir )  ! clem: close input ice restart file
          IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file 
@@ -287,4 +392,10 @@
       !
       CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice_os, zalb_ice_cs )
+      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+         CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice)
+      END IF
+      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
+         CALL wrk_dealloc( jpi,jpj, ztem_ice_all, zalb_ice_all, z_qsr_ice_all, z_qns_ice_all, z_qla_ice_all, z_dqns_ice_all, z_dqla_ice_all)
+      ENDIF
       !
       IF( nn_timing == 1 )  CALL timing_stop('sbc_ice_lim')
Index: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4098)
+++ branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4099)
@@ -84,5 +84,5 @@
       NAMELIST/namsbc/ nn_fsbc   , ln_ana    , ln_flx,  ln_blk_clio, ln_blk_core, ln_cpl,   &
          &             ln_blk_mfs, ln_apr_dyn, nn_ice,  nn_ice_embd, ln_dm2dc   , ln_rnf,   &
-         &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave , ln_sdw
+         &             ln_ssr    , nn_fwb    , ln_cdgw , ln_wave , ln_sdw, cn_iceflx
       !!----------------------------------------------------------------------
 
@@ -117,4 +117,5 @@
          WRITE(numout,*) '              MFS  bulk  formulation                     ln_blk_mfs  = ', ln_blk_mfs
          WRITE(numout,*) '              coupled    formulation (T if key_sbc_cpl)  ln_cpl      = ', ln_cpl
+         WRITE(numout,*) '              Flux handling over ice categories          cn_iceflx   = ', TRIM (cn_iceflx)
          WRITE(numout,*) '           Misc. options of sbc : '
          WRITE(numout,*) '              Patm gradient added in ocean & ice Eqs.    ln_apr_dyn  = ', ln_apr_dyn
@@ -128,4 +129,19 @@
       ENDIF
 
+      !   Flux handling over ice categories 
+      SELECT CASE ( TRIM (cn_iceflx))
+      CASE ('ave')
+         ln_iceflx_ave    = .TRUE.
+         ln_iceflx_linear = .FALSE.
+      CASE ('linear')
+         ln_iceflx_ave    = .FALSE.
+         ln_iceflx_linear = .TRUE.
+      CASE default
+         ln_iceflx_ave    = .FALSE.
+         ln_iceflx_linear = .FALSE.
+      END SELECT
+      IF(lwp) WRITE(numout,*) '              Fluxes averaged over all ice categories         ln_iceflx_ave    = ', ln_iceflx_ave
+      IF(lwp) WRITE(numout,*) '              Fluxes distributed linearly over ice categories ln_iceflx_linear = ', ln_iceflx_linear
+      !
       !                              ! allocate sbc arrays
       IF( sbc_oce_alloc() /= 0 )   CALL ctl_stop( 'STOP', 'sbc_init : unable to allocate sbc_oce arrays' )
@@ -166,4 +182,10 @@
       IF( ( nn_ice == 3 .OR. nn_ice == 4 ) .AND. nn_ice_embd == 0 )   &
          &   CALL ctl_stop( 'LIM3 and CICE sea-ice models require nn_ice_embd = 1 or 2' )
+
+      IF( ln_iceflx_ave .AND. ln_iceflx_linear ) &
+         &   CALL ctl_stop( ' ln_iceflx_ave and ln_iceflx_linear options are not compatible' )
+
+      IF( ( nn_ice ==3 .AND. lk_cpl) .AND. .NOT. ( ln_iceflx_ave .OR. ln_iceflx_linear ) ) &
+         &   CALL ctl_stop( ' With lim3 coupled, either ln_iceflx_ave or ln_iceflx_linear must be set to .TRUE.' )
       
       IF( ln_dm2dc )   nday_qsr = -1   ! initialisation flag
@@ -295,4 +317,6 @@
       CASE(  2 )   ;         CALL sbc_ice_lim_2( kt, nsbc )          ! LIM-2 ice model
       CASE(  3 )   ;         CALL sbc_ice_lim  ( kt, nsbc )          ! LIM-3 ice model
+      !is it useful?
+      !CASE(  4 )   ;         CALL sbc_ice_cice ( kt, nsbc )          ! CICE ice model
       END SELECT                                              
 
