Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref	(revision 4205)
@@ -49,5 +49,5 @@
    cw          =   5.0e-03 !  drag coefficient for oceanic stress
    angvg       =   0.0     !  turning angle for oceanic stress
-   pstar       =   4.0e+04 !  1st bulk-rheology parameter
+   pstar       =   2.0e+04 !  1st bulk-rheology parameter
    c_rhg       =  20.0     !  2nd bulk-rhelogy parameter
    etamn       =   0.0e+07 !  minimun value for viscosity
@@ -58,5 +58,5 @@
    telast      =9600.0     !  timescale for elastic waves, SB, 720.0
    alphaevp    =   1.0     !  coefficient for the solution of internal ice stresses
-   hminrhg     =   0.05     !  ice thickness (m) below which ice velocity equal ocean velocity
+   hminrhg     =   0.001   !  ice volume (a*h in m) below which ice velocity equal ocean velocity
 /
 !-----------------------------------------------------------------------
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/namelist_ref
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/namelist_ref	(revision 4205)
@@ -212,5 +212,5 @@
                            !  =1 use observed ice-cover      ,
                            !  =2 ice-model used                         ("key_lim3" or "key_lim2)
-   nn_ice_embd = 0         !  =0 levitating ice (no mass exchange, concentration/dilution effect)
+   nn_ice_embd = 1         !  =0 levitating ice (no mass exchange, concentration/dilution effect)
                            !  =1 levitating ice with mass and salt exchange but no presure effect
                            !  =2 embedded sea-ice (full salt and mass exchanges and pressure)
@@ -225,4 +225,5 @@
    ln_cdgw = .false.       !  Neutral drag coefficient read from wave model (T => fill namsbc_wave)
    ln_sdw  = .false.       !  Computation of 3D stokes drift                (T => fill namsbc_wave)
+   cn_iceflx = 'linear'    !  redistribution of solar input into ice categories during coupling ice/atm.
 /
 !-----------------------------------------------------------------------
@@ -283,6 +284,6 @@
    ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data
    rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow)
-   rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity
    rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.)
+   rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity 
                            !  in the calculation of the wind stress (0.=absolute winds or 1.=relative winds)
 /
@@ -532,10 +533,14 @@
                                           !  = 2, use tidal harmonic forcing data from files
                                           !  = 3, use external data AND tidal harmonic forcing
-    nn_dyn3d      =  0                    !  boundary conditions for baroclinic velocities
-    nn_dyn3d_dta  =  0                    !  = 0, bdy data are equal to the initial state
-                           !  = 1, bdy data are read in 'bdydata   .nc' files
-    nn_tra        =  1                    !  boundary conditions for T and S
+    nn_dyn3d      =  1                    !  boundary conditions for baroclinic velocities.....0=None/1=FRS
+    nn_dyn3d_dta  =  1                    !  = 0, bdy data are equal to the initial state
+                                          !  = 1, bdy data are read in 'bdydata   .nc' files
+    nn_tra        =  1                    !  boundary conditions for T and S...................0=None/1=FRS
     nn_tra_dta    =  1                    !  = 0, bdy data are equal to the initial state
-                           !  = 1, bdy data are read in 'bdydata   .nc' files
+                                          !  = 1, bdy data are read in 'bdydata   .nc' files
+!    nn_ice_lim      =  1                  !  boundary conditions for ice.......................0=None/1=FRS
+!    nn_ice_lim_dta  =  1                  !  = 0, bdy data are equal to the initial state
+!                                          !  = 1, bdy data are read in 'bdydata   .nc' files
+
     nn_rimwidth  = 10                      !  width of the relaxation zone
     ln_vol     = .false.                  !  total volume correction (see nn_volctl parameter)
@@ -554,4 +559,7 @@
    bn_tem  =    'amm12_bdyT_tra' ,         24        , 'votemper' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''
    bn_sal  =    'amm12_bdyT_tra' ,         24        , 'vosaline' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''
+!   bn_a_i  =    'amm12_bdyT_ice' ,         24        , 'ileadfra' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''
+!   bn_ht_i =    'amm12_bdyT_ice' ,         24        , 'iicethic' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''
+!   bn_ht_s =    'amm12_bdyT_ice' ,         24        , 'isnowthi' ,     .true.     , .false. ,  'daily'  ,    ''    ,   ''
    cn_dir  =    'bdydta/'
    ln_full_vel = .false.
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 4205)
@@ -396,9 +396,7 @@
    LOGICAL               , PUBLIC ::   ln_limdyn       !: flag for ice dynamics (T) or not (F)
    LOGICAL               , PUBLIC ::   ln_nicep        !: flag for sea-ice points output (T) or not (F)
-   REAL(wp)              , PUBLIC ::   hsndif          !: computation of temp. in snow (0) or not (9999)
-   REAL(wp)              , PUBLIC ::   hicdif          !: computation of temp. in ice (0) or not (9999)
    REAL(wp)              , PUBLIC ::   cai             !: atmospheric drag over sea ice
    REAL(wp)              , PUBLIC ::   cao             !: atmospheric drag over ocean
-   REAL(wp), DIMENSION(2), PUBLIC ::   acrit           !: minimum fraction for leads in north and south hemisphere
+   REAL(wp)              , PUBLIC ::   amax            !: maximum ice concentration
    !
    !!--------------------------------------------------------------------------
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90	(revision 4205)
@@ -89,7 +89,4 @@
       CALL lim_itd_ini                 ! ice thickness distribution initialization
       !
-      CALL lim_sbc_init                ! ice surface boundary condition   
-
-
       !                                ! Initial sea-ice state
       IF( .NOT.ln_rstart ) THEN              ! start from rest
@@ -106,5 +103,8 @@
       ENDIF
       !
+      CALL lim_sbc_init                ! ice surface boundary condition   
+      !
       fr_i(:,:) = at_i(:,:)           ! initialisation of sea-ice fraction
+      tn_ice(:,:,:) = t_su(:,:,:)
       !
       nstart = numit  + nn_fsbc      
@@ -248,5 +248,4 @@
       END DO
       !
-      tn_ice(:,:,:) = t_su(:,:,:)
       !
    END SUBROUTINE lim_itd_ini
Index: anches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limdia.F90	(revision 4204)
+++ 	(revision )
@@ -1,599 +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    !: instantaneous values of ice evolution or averaging ntmoy
-   INTEGER, PUBLIC  ::   ninfo    !: 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    ! 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     ! 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
-      INTEGER  ::   ios                         ! Local integer output status for namelist read
-      REAL(wp) ::   zxx0, zxx1    ! local scalars
-      CHARACTER(len=jpchinf) ::   titinf
-      CHARACTER(len=50)      ::   clname
-      !!
-      NAMELIST/namicedia/fmtinf, nfrinf, ninfo, ntmoy
-      !!-------------------------------------------------------------------
-      !
-      REWIND( numnam_ice_ref )              ! Namelist namicedia in reference namelist : Ice diagnostics in ice_evolu
-      READ  ( numnam_ice_ref, namicedia, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedia in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )              ! Namelist namicedia in configuration namelist : Ice diagnostics in ice_evolu 
-      READ  ( numnam_ice_cfg, namicedia, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedia in configuration namelist', lwp )
-      WRITE ( numoni, 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_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90	(revision 4205)
@@ -35,5 +35,4 @@
 
    PUBLIC   lim_istate      ! routine called by lim_init.F90
-
 
    !! * Module variables
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 4205)
@@ -123,5 +123,5 @@
 
       REAL(wp) ::   dtevp              ! time step for subcycling
-      REAL(wp) ::   dtotel, ecc2       ! square of yield ellipse eccenticity
+      REAL(wp) ::   dtotel, ecc2, ecci ! square of yield ellipse eccenticity
       REAL(wp) ::   z0, zr, zcca, zccb ! temporary scalars
       REAL(wp) ::   zu_ice2, zv_ice1   !
@@ -199,5 +199,5 @@
             zc1(ji,jj)    = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )
 #if defined key_lim3
-            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) * 0.5_wp
+            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj)
 #endif
 #if defined key_lim2
@@ -325,4 +325,5 @@
       !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
       ecc2 = ecc * ecc
+      ecci = 1. / ecc2
 
       !-Initialise stress tensor 
@@ -425,14 +426,17 @@
                !-Calculate stress tensor components zs1 and zs2 
                !-at centre of grid cells (see section 3.5 of CICE user's guide).
-               zs1(ji,jj) = ( zs1(ji,jj) &
-                  &          - dtotel*( ( 1.0 - alphaevp) * zs1(ji,jj) +    &
-                  &            ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) &
-                  * zpresh(ji,jj) ) )                          &       
-                  &        / ( 1.0 + alphaevp * dtotel )
-
-               zs2(ji,jj) = ( zs2(ji,jj)   &
-                  &          - dtotel*((1.0-alphaevp)*ecc2*zs2(ji,jj) -  &
-                  zdt(ji,jj)/deltat(ji,jj)*zpresh(ji,jj)) ) &
-                  &        / ( 1.0 + alphaevp*ecc2*dtotel )
+               zs1(ji,jj) = ( zs1(ji,jj) - dtotel*( ( 1._wp - alphaevp) * zs1(ji,jj) +   &
+                  &          ( delta / deltat(ji,jj) - zdd(ji,jj) / deltat(ji,jj) ) * zpresh(ji,jj) ) )  &       
+                  &          / ( 1._wp + alphaevp * dtotel )
+
+               zs2(ji,jj) = ( zs2(ji,jj) - dtotel * ( ( 1._wp - alphaevp ) * ecc2 * zs2(ji,jj) -   &
+                             zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )   &
+                  &          / ( 1._wp + alphaevp * ecc2 * dtotel )
+
+               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
+               !zs1(ji,jj) = ( zs1(ji,jj) + dtotel * ( ( zdd(ji,jj) / deltat(ji,jj) - delta / deltat(ji,jj) )  &
+               !   &         * zpresh(ji,jj) ) ) / ( 1._wp + dtotel )
+               !zs2(ji,jj) = ( zs2(ji,jj) + dtotel * ( ecci * zdt(ji,jj) / deltat(ji,jj) * zpresh(ji,jj) ) )  &
+               !   &         / ( 1._wp + dtotel )
 
             END DO
@@ -468,8 +472,12 @@
 
                !-Calculate stress tensor component zs12 at corners (see section 3.5 of CICE user's guide).
-               zs12(ji,jj) = ( zs12(ji,jj)      &
-                  &        - dtotel*( (1.0-alphaevp)*ecc2*zs12(ji,jj) - zds(ji,jj) / &
-                  &          ( 2.0*deltac(ji,jj) ) * zpreshc(ji,jj))) &
-                  &         / ( 1.0 + alphaevp*ecc2*dtotel ) 
+               zs12(ji,jj) = ( zs12(ji,jj) - dtotel * ( (1.0-alphaevp) * ecc2 * zs12(ji,jj) - zds(ji,jj) /  &
+                  &          ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+                  &          / ( 1._wp + alphaevp * ecc2 * dtotel ) 
+
+               ! new formulation from S. Bouillon to help stabilizing the code (no need of alphaevp)
+               !zs12(ji,jj) = ( zs12(ji,jj) + dtotel *  &
+               !   &          ( ecci * zds(ji,jj) / ( 2._wp * deltac(ji,jj) ) * zpreshc(ji,jj) ) )  &
+               !   &          / ( 1.0 + dtotel ) 
 
             END DO ! ji
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90	(revision 4205)
@@ -17,4 +17,5 @@
    !!----------------------------------------------------------------------
    USE ice            ! sea-ice variables
+   USE oce     , ONLY :  snwice_mass, snwice_mass_b
    USE par_ice        ! sea-ice parameters
    USE dom_oce        ! ocean domain
@@ -159,12 +160,12 @@
       END DO
 
-      CALL iom_rstput( iter, nitrst, numriw, 'u_ice'     , u_ice      )
-      CALL iom_rstput( iter, nitrst, numriw, 'v_ice'     , v_ice      )
-      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq'     , fsbbq      )
-      CALL iom_rstput( iter, nitrst, numriw, 'iatte'     , iatte      ) ! clem modif
-      CALL iom_rstput( iter, nitrst, numriw, 'oatte'     , oatte      ) ! clem modif
-      CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
-      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
-      CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
+      CALL iom_rstput( iter, nitrst, numriw, 'u_ice'        , u_ice      )
+      CALL iom_rstput( iter, nitrst, numriw, 'v_ice'        , v_ice      )
+      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq'        , fsbbq      )
+      CALL iom_rstput( iter, nitrst, numriw, 'stress1_i'    , stress1_i  )
+      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i'    , stress2_i  )
+      CALL iom_rstput( iter, nitrst, numriw, 'stress12_i'   , stress12_i )
+      CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass'  , snwice_mass )   !clem modif
+      CALL iom_rstput( iter, nitrst, numriw, 'snwice_mass_b', snwice_mass_b ) !clem modif
 
       DO jl = 1, jpl 
@@ -371,52 +372,5 @@
          CALL iom_get( numrir, jpdom_autoglo, znam , z2d )
          t_su(:,:,jl) = z2d(:,:)
-         tn_ice (:,:,:) = t_su (:,:,:)
-      END DO
-
-      DO jl = 1, jpl 
-         CALL lbc_lnk( smv_i(:,:,jl) , 'T' ,  1. )
-         CALL lbc_lnk( v_i  (:,:,jl) , 'T' ,  1. )
-         CALL lbc_lnk( a_i  (:,:,jl) , 'T' ,  1. )
-      END DO
-
-      ! we first with bulk ice salinity
-      DO jl = 1, jpl
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-               zindb          = MAX( 0.0 , SIGN( 1.0 , v_i(ji,jj,jl) - 1.0e-4 ) ) 
-               sm_i(ji,jj,jl) = smv_i(ji,jj,jl) / MAX(v_i(ji,jj,jl),1.0e-6) * zindb
-               ht_i(ji,jj,jl) = v_i(ji,jj,jl)   / MAX(a_i(ji,jj,jl),1.0e-6) * zindb
-            END DO
-         END DO
-      END DO
-
-      DO jk = 1, nlay_i
-         s_i(:,:,jk,:) = sm_i(:,:,:)
-      END DO
-
-      IF( num_sal == 2 ) THEN      ! Salinity profile
-         DO jl = 1, jpl
-            DO jk = 1, nlay_i
-               DO jj = 1, jpj
-                  DO ji = 1, jpi
-                     zs_inf        = sm_i(ji,jj,jl)
-                     z_slope_s     = 2._wp * sm_i(ji,jj,jl) / MAX( 0.01_wp , ht_i(ji,jj,jl) )
-                     !- slope of the salinity profile
-                     zs_zero(jk)   = z_slope_s * ( REAL(jk,wp) - 0.5_wp ) * ht_i(ji,jj,jl) / REAL(nlay_i,wp)
-                     zsmax = 4.5_wp
-                     zsmin = 3.5_wp
-                     IF(     sm_i(ji,jj,jl) < zsmin ) THEN
-                        zalpha = 1._wp
-                     ELSEIF( sm_i(ji,jj,jl) < zsmax ) THEN
-                        zalpha = sm_i(ji,jj,jl) / ( zsmin - zsmax ) + zsmax / ( zsmax - zsmin )
-                     ELSE
-                        zalpha = 0._wp
-                     ENDIF
-                     s_i(ji,jj,jk,jl) = zalpha * zs_zero(jk) + ( 1._wp - zalpha ) * zs_inf
-                  END DO
-               END DO
-            END DO
-         END DO
-      ENDIF
+      END DO
 
       DO jl = 1, jpl 
@@ -440,9 +394,9 @@
       CALL iom_get( numrir, jpdom_autoglo, 'v_ice'     , v_ice      )
       CALL iom_get( numrir, jpdom_autoglo, 'fsbbq'     , fsbbq      )
-      CALL iom_get( numrir, jpdom_autoglo, 'iatte'     , iatte      ) ! clem modif
-      CALL iom_get( numrir, jpdom_autoglo, 'oatte'     , oatte      ) ! clem modif
       CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
       CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
       CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
+      CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass'  , snwice_mass )   !clem modif
+      CALL iom_get( numrir, jpdom_autoglo, 'snwice_mass_b', snwice_mass_b ) !clem modif
 
       DO jl = 1, jpl 
@@ -568,5 +522,5 @@
       END DO
       !
-      !clem CALL iom_close( numrir )
+      !CALL iom_close( numrir ) !clem: closed in sbcice_lim.F90
       !
       CALL wrk_dealloc( nlay_i, zs_zero )
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 4205)
@@ -38,5 +38,5 @@
    USE cpl_oasis3, ONLY : lk_cpl
    USE traqsr           ! clem: add penetration of solar flux into the calculation of heat budget
-   USE oce,        ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n
+   USE oce,        ONLY : iatte, oatte, sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n
    USE dom_ice,    ONLY : tms
    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
@@ -249,5 +249,5 @@
             ! mass flux at the ocean/ice interface (sea ice fraction)
             zemp_snw = rdm_snw(ji,jj) * r1_rdtice                         ! snow melting = pure water that enters the ocean
-            zfmm     = rdm_ice(ji,jj) * r1_rdtice                         ! Freezing minus mesting  
+            zfmm     = rdm_ice(ji,jj) * r1_rdtice                         ! Freezing minus melting  
 
             fmmflx(ji,jj) = zfmm                                     ! F/M mass flux save at least for biogeochemical model
@@ -416,51 +416,55 @@
       ENDIF
       ! clem modif
-      iatte(:,:) = 1._wp
-      oatte(:,:) = 1._wp
-      !
-      !                                      ! embedded sea ice
-      IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
-         snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )
-         snwice_mass_b(:,:) = snwice_mass(:,:)
-      ELSE
-         snwice_mass  (:,:) = 0.0_wp         ! no mass exchanges
-         snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges
-      ENDIF
-      IF( nn_ice_embd == 2  .AND.         &  ! full embedment (case 2) & no restart
-         &  .NOT. ln_rstart ) THEN           ! deplete the initial ssh below sea-ice area
-         sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
-         sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
-         !
-         ! Note: Changed the initial values of sshb and sshn=>  need to recompute ssh[u,v,f]_[b,n] 
-         !       which were previously set in domvvl
-         IF ( lk_vvl ) THEN            ! Is this necessary? embd 2 should be restricted to vvl only???
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
-                  zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
-                  zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
-                  zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
-                  sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
-                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
-                  sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
-                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
-                  sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
-                     &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
-                  sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
-                     &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
+      IF( .NOT. ln_rstart ) THEN
+         iatte(:,:) = 1._wp
+         oatte(:,:) = 1._wp
+      ENDIF
+      !
+      ! clem: snwice_mass in the restart file now
+      IF( .NOT. ln_rstart ) THEN
+         !                                      ! embedded sea ice
+         IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass
+            snwice_mass  (:,:) = tms(:,:) * ( rhosn * vt_s(:,:) + rhoic * vt_i(:,:)  )
+            snwice_mass_b(:,:) = snwice_mass(:,:)
+         ELSE
+            snwice_mass  (:,:) = 0.0_wp         ! no mass exchanges
+            snwice_mass_b(:,:) = 0.0_wp         ! no mass exchanges
+         ENDIF
+         IF( nn_ice_embd == 2 ) THEN            ! full embedment (case 2) deplete the initial ssh below sea-ice area
+            sshn(:,:) = sshn(:,:) - snwice_mass(:,:) * r1_rau0
+            sshb(:,:) = sshb(:,:) - snwice_mass(:,:) * r1_rau0
+            !
+            ! Note: Changed the initial values of sshb and sshn=>  need to recompute ssh[u,v,f]_[b,n] 
+            !       which were previously set in domvvl
+            IF ( lk_vvl ) THEN            ! Is this necessary? embd 2 should be restricted to vvl only???
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1                    ! caution: use of Vector Opt. not possible
+                     zcoefu = 0.5  * umask(ji,jj,1) / ( e1u(ji,jj) * e2u(ji,jj) )
+                     zcoefv = 0.5  * vmask(ji,jj,1) / ( e1v(ji,jj) * e2v(ji,jj) )
+                     zcoeff = 0.25 * umask(ji,jj,1) * umask(ji,jj+1,1)
+                     sshu_b(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshb(ji  ,jj)     &
+                        &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshb(ji+1,jj) )
+                     sshv_b(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshb(ji,jj  )     &
+                        &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshb(ji,jj+1) )
+                     sshu_n(ji,jj) = zcoefu * ( e1t(ji  ,jj) * e2t(ji  ,jj) * sshn(ji  ,jj)     &
+                        &                     + e1t(ji+1,jj) * e2t(ji+1,jj) * sshn(ji+1,jj) )
+                     sshv_n(ji,jj) = zcoefv * ( e1t(ji,jj  ) * e2t(ji,jj  ) * sshn(ji,jj  )     &
+                        &                     + e1t(ji,jj+1) * e2t(ji,jj+1) * sshn(ji,jj+1) )
+                  END DO
                END DO
-            END DO
-            CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
-            CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
-            DO jj = 1, jpjm1
-               DO ji = 1, jpim1      ! NO Vector Opt.
-                  sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
-                       &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
-                       &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
-                       &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
+               CALL lbc_lnk( sshu_b, 'U', 1. )   ;   CALL lbc_lnk( sshu_n, 'U', 1. )
+               CALL lbc_lnk( sshv_b, 'V', 1. )   ;   CALL lbc_lnk( sshv_n, 'V', 1. )
+               DO jj = 1, jpjm1
+                  DO ji = 1, jpim1      ! NO Vector Opt.
+                     sshf_n(ji,jj) = 0.5  * umask(ji,jj,1) * umask(ji,jj+1,1)                   &
+                          &               / ( e1f(ji,jj  ) * e2f(ji,jj  ) )                     &
+                          &               * ( e1u(ji,jj  ) * e2u(ji,jj  ) * sshu_n(ji,jj  )     &
+                          &                 + e1u(ji,jj+1) * e2u(ji,jj+1) * sshu_n(ji,jj+1) )
+                  END DO
                END DO
-            END DO
-            CALL lbc_lnk( sshf_n, 'F', 1. )
-          ENDIF
-      ENDIF
+               CALL lbc_lnk( sshf_n, 'F', 1. )
+            ENDIF
+         ENDIF
+      ENDIF ! .NOT. ln_rstart
       !
 !!?      IF( .NOT. ln_rstart ) THEN           ! delete the initial ssh below sea-ice area
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90	(revision 4205)
@@ -22,4 +22,5 @@
    USE phycst         ! physical constants
    USE dom_oce        ! ocean space and time domain variables
+   USE oce     , ONLY :  iatte, oatte
    USE ice            ! LIM: sea-ice variables
    USE par_ice        ! LIM: sea-ice parameters
@@ -220,5 +221,5 @@
             !
             ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx)
-            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) )
+            qcmif  (ji,jj) =  rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) )
             !
             ! oceanic heat flux (limthd_dh)
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90	(revision 4205)
@@ -656,5 +656,5 @@
          DO ji = kideb , kiut
             ! snow temperatures      
-            IF (ht_s_b(ji).GT.0) &
+            IF (ht_s_b(ji).GT.0._wp) &
                t_s_b(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) &
                *  t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) &
@@ -692,5 +692,5 @@
             DO ji = kideb , kiut
                ztmelt_i        = -tmut * s_i_b(ji,layer) + rtt 
-               t_i_b(ji,layer) =  MAX(MIN(t_i_b(ji,layer),ztmelt_i),190.0)
+               t_i_b(ji,layer) =  MAX(MIN(t_i_b(ji,layer),ztmelt_i), 190._wp)
                zerrit(ji)      =  MAX(zerrit(ji),ABS(t_i_b(ji,layer) - ztitemp(ji,layer)))
             END DO
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90	(revision 4205)
@@ -74,5 +74,5 @@
       REAL(wp) ::   zindb  , zindsn , zindic, zindh, zinda      ! local scalar
       REAL(wp) ::   zusvosn, zusvoic, zbigval     !   -      -
-      REAL(wp) ::   zcfl , zusnit , zrtt          !   -      -
+      REAL(wp) ::   zcfl , zusnit                 !   -      -
       REAL(wp) ::   ze   , zsal   , zage          !   -      -
       !
@@ -450,5 +450,4 @@
                   zusvosn         = 1.0/MAX( v_s(ji,jj,jl) , epsi16 )
                   zusvoic         = 1.0/MAX( v_i(ji,jj,jl) , epsi16 )
-                  zrtt            = 173.15 * rone 
                   zindsn          = MAX( rzero, SIGN( rone, v_s(ji,jj,jl) - epsi10 ) )
                   zindic          = MAX( rzero, SIGN( rone, v_i(ji,jj,jl) - epsi10 ) )
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90	(revision 4205)
@@ -180,5 +180,4 @@
                ht_s(ji,jj,jl) = v_s (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb
                o_i(ji,jj,jl)  = oa_i(ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi10 ) * zindb
-               a_i(ji,jj,jl) = a_i (ji,jj,jl) * zindb ! clem correction
             END DO
          END DO
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90	(revision 4205)
@@ -22,5 +22,4 @@
    REAL(wp), PUBLIC ::   hicmin      !: (REMOVE)
    REAL(wp), PUBLIC ::   hiclim      !: minimum ice thickness
-   REAL(wp), PUBLIC ::   amax        !: maximum lead fraction
    REAL(wp), PUBLIC ::   sbeta       !: numerical scheme for diffusion in ice  (REMOVE)
    REAL(wp), PUBLIC ::   parlat      !: (REMOVE)
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90	(revision 4205)
@@ -260,5 +260,4 @@
             CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )         
          ENDIF
-                                                    
          !
          CALL dia_hsb_rst( nit000, 'READ' )  !* read or initialize all required files
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/IOM/iom.F90	(revision 4205)
@@ -406,5 +406,5 @@
 
 
-   FUNCTION iom_varid ( kiomid, cdvar, kdimsz, ldstop )  
+   FUNCTION iom_varid ( kiomid, cdvar, kdimsz, kndims, ldstop )  
       !!-----------------------------------------------------------------------
       !!                  ***  FUNCTION  iom_varid  ***
@@ -415,4 +415,5 @@
       CHARACTER(len=*)     , INTENT(in   )           ::   cdvar    ! name of the variable
       INTEGER, DIMENSION(:), INTENT(  out), OPTIONAL ::   kdimsz   ! size of the dimensions
+      INTEGER,               INTENT(  out), OPTIONAL ::   kndims   ! size of the dimensions
       LOGICAL              , INTENT(in   ), OPTIONAL ::   ldstop   ! stop if looking for non-existing variable (default = .TRUE.)
       !
@@ -446,5 +447,5 @@
                   SELECT CASE (iom_file(kiomid)%iolib)
                   CASE (jpioipsl )   ;   iom_varid = iom_ioipsl_varid( kiomid, cdvar, iiv, kdimsz )
-                  CASE (jpnf90   )   ;   iom_varid = iom_nf90_varid  ( kiomid, cdvar, iiv, kdimsz )
+                  CASE (jpnf90   )   ;   iom_varid = iom_nf90_varid  ( kiomid, cdvar, iiv, kdimsz, kndims )
                   CASE (jprstdimg)   ;   iom_varid = -1   ! all variables are listed in iom_file
                   CASE DEFAULT   
@@ -467,4 +468,5 @@
                   ENDIF
                ENDIF
+               IF( PRESENT(kndims) )  kndims = iom_file(kiomid)%ndims(iiv)
             ENDIF
          ENDIF
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/IOM/iom_def.F90	(revision 4205)
@@ -44,5 +44,5 @@
 
    INTEGER, PARAMETER, PUBLIC ::   jpmax_files  = 100   !: maximum number of simultaneously opened file
-   INTEGER, PARAMETER, PUBLIC ::   jpmax_vars   = 360  !: maximum number of variables in one file
+   INTEGER, PARAMETER, PUBLIC ::   jpmax_vars   = 600  !: maximum number of variables in one file
    INTEGER, PARAMETER, PUBLIC ::   jpmax_dims   =  4   !: maximum number of dimensions for one variable
    INTEGER, PARAMETER, PUBLIC ::   jpmax_digits =  5   !: maximum number of digits for the cpu number in the file name
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/IOM/restart.F90	(revision 4205)
@@ -132,4 +132,6 @@
                      CALL iom_rstput( kt, nitrst, numrow, 'rhd'    , rhd       )
 #endif
+                     CALL iom_rstput( kt, nitrst, numrow, 'iatte'  , iatte     ) !clem modif
+                     CALL iom_rstput( kt, nitrst, numrow, 'oatte'  , oatte     ) !clem modif
       IF( kt == nitrst ) THEN
          CALL iom_close( numrow )     ! close the restart file (only at last time step)
@@ -232,4 +234,7 @@
       ENDIF
       !
+      CALL iom_get( numror, jpdom_autoglo, 'iatte' , iatte ) ! clem modif
+      CALL iom_get( numror, jpdom_autoglo, 'oatte' , oatte ) ! clem modif
+      !
    END SUBROUTINE rst_read
 
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 4205)
@@ -82,5 +82,4 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::  qsr_hc , qsr_hc_b   !: heat content trend due to qsr flux     [K.m/s] jpi,jpj,jpk
    !!
-   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   oatte, iatte      !: clem attenuation coef of the input solar flux [unitless]
    !!
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tprecip           !: total precipitation                          [Kg/m2/s]
@@ -127,6 +126,5 @@
          !
       ALLOCATE( rnf  (jpi,jpj) , sbc_tsc  (jpi,jpj,jpts) , qsr_hc  (jpi,jpj,jpk) ,     &
-         &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) ,     &
-         &      iatte(jpi,jpj) , oatte    (jpi,jpj)                              , STAT=ierr(3) )
+         &      rnf_b(jpi,jpj) , sbc_tsc_b(jpi,jpj,jpts) , qsr_hc_b(jpi,jpj,jpk) , STAT=ierr(3) )
          !
       ALLOCATE( tprecip(jpi,jpj) , sprecip(jpi,jpj) , fr_i(jpi,jpj) ,     &
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 4205)
@@ -153,11 +153,9 @@
       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 defined key_coupled
+      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_alloc( jpi,jpj,jpl, zalb_ice)
+      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
+         &   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
@@ -169,5 +167,5 @@
          !
          IF( ln_nicep ) THEN      ! control print at a given point
-            jiindx = 15   ;   jjindx = 46
+            jiindx = 6   ;   jjindx = 47
             WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx
          ENDIF
@@ -193,4 +191,5 @@
          IF ( ln_cpl ) zalb_ice (:,:,:) = 0.5 * ( zalb_ice_cs (:,:,:) +  zalb_ice_os (:,:,:) )
          
+#if defined key_coupled
          IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
             !
@@ -200,5 +199,6 @@
             !
          ENDIF
-                                                     ! Bulk formulea - provides the following fields:
+#endif
+                                               ! Bulk formulea - provides the following fields:
          ! utau_ice, vtau_ice : surface ice stress                     (U- & V-points)   [N/m2]
          ! qsr_ice , qns_ice  : solar & non solar heat flux over ice   (T-point)         [W/m2]
@@ -238,4 +238,5 @@
 
          ! Average over all categories
+#if defined key_coupled
          IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) THEN
 
@@ -267,5 +268,5 @@
             END IF
          END IF
-
+#endif
          !                                           !----------------------!
          !                                           ! LIM-3  time-stepping !
@@ -283,5 +284,7 @@
          old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content
          old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content
-
+         !
+         old_u_ice(:,:) = u_ice(:,:)
+         old_v_ice(:,:) = v_ice(:,:)
          !                                           ! intialisation to zero    !!gm is it truly necessary ???
          d_a_i_thd  (:,:,:)   = 0._wp   ;   d_a_i_trp  (:,:,:)   = 0._wp
@@ -292,4 +295,7 @@
          d_smv_i_thd(:,:,:)   = 0._wp   ;   d_smv_i_trp(:,:,:)   = 0._wp
          d_oa_i_thd (:,:,:)   = 0._wp   ;   d_oa_i_trp (:,:,:)   = 0._wp
+         !
+         d_u_ice_dyn(:,:) = 0._wp
+         d_v_ice_dyn(:,:) = 0._wp
          !
          sfx    (:,:) = 0._wp   ;   sfx_thd  (:,:) = 0._wp
@@ -373,5 +379,7 @@
                           CALL lim_wri( 1  )              ! Ice outputs 
 !clem # endif
-         IF( kt == nit000 )   CALL iom_close( numrir )  ! clem: close input ice restart file
+         IF( kt == nit000 .AND. ln_rstart )   &
+            &             CALL iom_close( numrir )        ! clem: close input ice restart file
+         !
          IF( lrst_ice )   CALL lim_rst_write( kt )        ! Ice restart file 
                           CALL lim_var_glo2eqv            ! ???
@@ -392,10 +400,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 defined key_coupled
+      IF ( ln_cpl .OR. ln_iceflx_ave .OR. ln_iceflx_linear ) CALL wrk_dealloc( jpi,jpj,jpl, zalb_ice)
+      IF ( ln_iceflx_ave .OR. ln_iceflx_linear ) &
+         &    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_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 4205)
@@ -136,5 +136,6 @@
       ENDIF
 
-      !   Flux handling over ice categories 
+      !   Flux handling over ice categories
+#if defined key_coupled 
       SELECT CASE ( TRIM (cn_iceflx))
       CASE ('ave')
@@ -150,4 +151,5 @@
       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
+#endif
       !
       !                              ! allocate sbc arrays
@@ -189,11 +191,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 defined key_coupled
       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.' )
-      
+#endif      
       IF( ln_dm2dc )   nday_qsr = -1   ! initialisation flag
 
@@ -247,5 +248,4 @@
          IF( nsbc ==  6 )   WRITE(numout,*) '              MFS Bulk formulation'
       ENDIF
-
       !
                           CALL sbc_ssm_init               ! Sea-surface mean fields initialisation
@@ -255,5 +255,4 @@
       IF( nn_ice == 4 )   CALL cice_sbc_init( nsbc )      ! CICE initialisation
       !
-
    END SUBROUTINE sbc_init
 
@@ -329,5 +328,5 @@
       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
+      CASE(  4 )   ;         CALL sbc_ice_cice ( kt, nsbc )          ! CICE ice model
       END SELECT                                              
 
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/TRA/traqsr.F90	(revision 4205)
@@ -45,10 +45,9 @@
    LOGICAL , PUBLIC ::   ln_qsr_2bd   !: 2 band         light absorption flag
    LOGICAL , PUBLIC ::   ln_qsr_bio   !: bio-model      light absorption flag
+   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem)
    INTEGER , PUBLIC ::   nn_chldta    !: use Chlorophyll data (=1) or not (=0)
    REAL(wp), PUBLIC ::   rn_abs       !: fraction absorbed in the very near surface (RGB & 2 bands)
    REAL(wp), PUBLIC ::   rn_si0       !: very near surface depth of extinction      (RGB & 2 bands)
    REAL(wp), PUBLIC ::   rn_si1       !: deepest depth of extinction (water type I)       (2 bands)
-   LOGICAL , PUBLIC ::   ln_qsr_ice   !: light penetration for ice-model LIM3 (clem)
-
    
    ! Module variables
@@ -395,6 +394,8 @@
       !
       ! clem init for oatte and iatte
-      oatte(:,:) = 1._wp
-      iatte(:,:) = 1._wp
+      IF( .NOT. ln_rstart ) THEN
+         oatte(:,:) = 1._wp
+         iatte(:,:) = 1._wp
+      ENDIF
       !
       CALL wrk_alloc( jpi, jpj,      zekb, zekg, zekr        ) 
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/oce.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/oce.F90	(revision 4205)
@@ -57,4 +57,7 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   snwice_fmass       !: time evolution of mass of snow+ice               [Kg/m2/s]
 
+   !! arrays related to penetration of solar fluxes to calculate the heat budget for sea ice
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   oatte, iatte       !: attenuation coef of the input solar flux [unitless]
+
    !!----------------------------------------------------------------------
    !! NEMO/OPA 4.0 , NEMO Consortium (2011)
@@ -68,5 +71,5 @@
       !!                   ***  FUNCTION oce_alloc  ***
       !!----------------------------------------------------------------------
-      INTEGER :: ierr(3)
+      INTEGER :: ierr(4)
       !!----------------------------------------------------------------------
       !
@@ -90,6 +93,7 @@
          &      gru(jpi,jpj)      , grv(jpi,jpj)                      , STAT=ierr(2) )
          !
-      ALLOCATE( snwice_mass(jpi,jpj)  , snwice_mass_b(jpi,jpj),             &
-         &      snwice_fmass(jpi,jpj), STAT= ierr(3) )
+      ALLOCATE( snwice_mass(jpi,jpj) , snwice_mass_b(jpi,jpj), snwice_fmass(jpi,jpj) , STAT=ierr(3) )
+         !
+      ALLOCATE( iatte(jpi,jpj) , oatte(jpi,jpj) , STAT=ierr(4) )
          !
       oce_alloc = MAXVAL( ierr )
Index: /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/par_oce.F90
===================================================================
--- /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/par_oce.F90	(revision 4204)
+++ /branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/OPA_SRC/par_oce.F90	(revision 4205)
@@ -56,5 +56,5 @@
    !!  Values set to pp_to_be_computed  indicates that variables will be computed in domzgr
    REAL(wp)      ::   pp_not_used       = 999999._wp   !: vertical grid parameter
-   REAL(wp)      ::   pp_to_be_computed = 999998._wp   !:    -      -       -
+   REAL(wp)      ::   pp_to_be_computed = 999999._wp   !:    -      -       -
 
 
