Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/dom_ice_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/dom_ice_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/dom_ice_2.F90	(revision 2528)
@@ -5,10 +5,13 @@
    !!======================================================================
    !! History :   2.0  !  03-08  (C. Ethe)  Free form and module
+   !!             3.3  !  2009-05 (G.Garric, C. Bricaud) addition of lim2_evp case
    !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
-   !!   LIM 2.0, UCL-LOCEAN-IPSL (2005)
+   !!   'key_lim2'                                       LIM2 sea-ice model
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
    USE par_ice_2
@@ -22,13 +25,24 @@
       !                                        !  (otherwise = jpj+10 (SH) or -10 (SH) )
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor , fcor,   &  !: coriolis factor and coeficient
-      &                                              covrai ,         &  !: sine of geographic latitude
-      &                                              area   ,         &  !: surface of grid cell 
-      &                                              tms    , tmu        !: temperature and velocity points masks
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght   ,         &  !: weight of the 4 neighbours to compute averages
-      &                                              akappa , bkappa     !: first and third group of metric coefficients
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd   !: second group of metric coefficients
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor , fcor     !: coriolis factor and coeficient
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   covrai            !: sine of geographic latitude
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   area              !: surface of grid cell 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tms    , tmu      !: temperature and velocity points masks
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght              !: weight of the 4 neighbours to compute averages
 
+
+# if defined key_lim2_vp
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   akappa , bkappa   !: first and third group of metric coefficients
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd            !: second group of metric coefficients
+# else
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tmv    , tmf      !: y-velocity and F-points masks
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tmi               !: ice mask: =1 if ice thick > 0
+# endif
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option          Empty module         NO LIM-2 sea-ice model
+   !!----------------------------------------------------------------------
+#endif
    !!======================================================================
-#endif
 END MODULE dom_ice_2
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 2528)
@@ -4,5 +4,6 @@
    !! Sea Ice physics:  diagnostics variables of ice defined in memory
    !!=====================================================================
-   !! History :  2.0  !  03-08  (C. Ethe)  F90: Free form and module
+   !! History :  2.0  !  2003-08  (C. Ethe)  F90: Free form and module
+   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -10,53 +11,79 @@
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006)
-   !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_ice_2          ! LIM sea-ice parameters
 
    IMPLICIT NONE
    PRIVATE
+   
+   INTEGER , PUBLIC ::   numit     !: ice iteration index
+   REAL(wp), PUBLIC ::   rdt_ice   !: ice time step
 
-   !!* Share parameters namelist (namicerun read in iceini) *
+   !                                                                     !!* namicerun read in iceini  *
    CHARACTER(len=32)     , PUBLIC ::   cn_icerst_in  = "restart_ice_in"   !: suffix of ice restart name (input)
    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_limdmp     = .FALSE.            !: Ice damping
-   REAL(wp)              , PUBLIC ::   hsndif        = 0.e0               !: computation of temp. in snow (0) or not (9999)
-   REAL(wp)              , PUBLIC ::   hicdif        = 0.e0               !: computation of temp. in ice (0) or not (9999)
-   REAL(wp), DIMENSION(2), PUBLIC ::   acrit = (/ 1.e-06 , 1.e-06 /)      !: minimum fraction for leads in 
-   !                                                                      !: north and south hemisphere
-   !!* ice-dynamic namelist (namicedyn) *
-   INTEGER , PUBLIC ::   nbiter = 1         !: number of sub-time steps for relaxation
-   INTEGER , PUBLIC ::   nbitdr = 250       !: maximum number of iterations for relaxation
-   REAL(wp), PUBLIC ::   rdt_ice            !: ice time step
-   REAL(wp), PUBLIC ::   epsd   = 1.0e-20   !: tolerance parameter for dynamic
-   REAL(wp), PUBLIC ::   alpha  = 0.5       !: coefficient for semi-implicit coriolis
-   REAL(wp), PUBLIC ::   dm     = 0.6e+03   !: diffusion constant for dynamics
-   REAL(wp), PUBLIC ::   om     = 0.5       !: relaxation constant
-   REAL(wp), PUBLIC ::   resl   = 5.0e-05   !: maximum value for the residual of relaxation
-   REAL(wp), PUBLIC ::   cw     = 5.0e-03   !: drag coefficient for oceanic stress
-   REAL(wp), PUBLIC ::   angvg  = 0.e0      !: turning angle for oceanic stress
-   REAL(wp), PUBLIC ::   pstar  = 1.0e+04   !: first bulk-rheology parameter
-   REAL(wp), PUBLIC ::   c_rhg  = 20.e0     !: second bulk-rhelogy parameter
-   REAL(wp), PUBLIC ::   etamn  = 0.e+07    !: minimun value for viscosity
-   REAL(wp), PUBLIC ::   creepl = 2.e-08    !: creep limit
-   REAL(wp), PUBLIC ::   ecc    = 2.e0      !: eccentricity of the elliptical yield curve
-   REAL(wp), PUBLIC ::   ahi0   = 350.e0    !: sea-ice hor. eddy diffusivity coeff. (m2/s)
+   LOGICAL               , PUBLIC ::   ln_nicep      = .TRUE.             !: flag grid points output (T) or not (F)
+   REAL(wp)              , PUBLIC ::   hsndif        = 0._wp              !: snow temp. computation (0) or not (9999)
+   REAL(wp)              , PUBLIC ::   hicdif        = 0._wp              !: ice  temp. computation (0) or not (9999)
+   REAL(wp), DIMENSION(2), PUBLIC ::   acrit = (/ 1.e-6_wp , 1.e-6_wp /)  !: minimum lead fraction in the 2 hemisphere
+   
+   !                                          !!* ice-dynamic namelist (namicedyn) *
+   INTEGER , PUBLIC ::   nbiter =    1         !: number of sub-time steps for relaxation
+   INTEGER , PUBLIC ::   nbitdr =  250         !: maximum number of iterations for relaxation
+   INTEGER , PUBLIC ::   nevp   =  360         !: number of EVP subcycling iterations
+   INTEGER , PUBLIC ::   telast = 3600         !: timescale for EVP elastic waves
+   REAL(wp), PUBLIC ::   epsd   = 1.0e-20_wp   !: tolerance parameter for dynamic
+   REAL(wp), PUBLIC ::   alpha  = 0.5_wp       !: coefficient for semi-implicit coriolis
+   REAL(wp), PUBLIC ::   dm     = 0.6e+03_wp   !: diffusion constant for dynamics
+   REAL(wp), PUBLIC ::   om     = 0.5_wp       !: relaxation constant
+   REAL(wp), PUBLIC ::   resl   = 5.0e-05_wp   !: maximum value for the residual of relaxation
+   REAL(wp), PUBLIC ::   cw     = 5.0e-03_wp   !: drag coefficient for oceanic stress
+   REAL(wp), PUBLIC ::   angvg  = 0._wp        !: turning angle for oceanic stress
+   REAL(wp), PUBLIC ::   pstar  = 1.0e+04_wp   !: first bulk-rheology parameter
+   REAL(wp), PUBLIC ::   c_rhg  = 20._wp       !: second bulk-rhelogy parameter
+   REAL(wp), PUBLIC ::   etamn  = 0._wp        !: minimun value for viscosity
+   REAL(wp), PUBLIC ::   creepl = 2.e-08_wp    !: creep limit
+   REAL(wp), PUBLIC ::   ecc    = 2._wp        !: eccentricity of the elliptical yield curve
+   REAL(wp), PUBLIC ::   ahi0   = 350._wp      !: sea-ice hor. eddy diffusivity coeff. (m2/s)
+   REAL(wp), PUBLIC ::   alphaevp = 1._wp      !: coefficient for the solution of EVP int. stresses
 
-   REAL(wp), PUBLIC ::   usecc2             !:  = 1.0 / ( ecc * ecc )
-   REAL(wp), PUBLIC ::   rhoco              !: = rau0 * cw
-   REAL(wp), PUBLIC ::   sangvg, cangvg     !: sin and cos of the turning angle for ocean stress
-   REAL(wp), PUBLIC ::   pstarh             !: pstar / 2.0
+   REAL(wp), PUBLIC ::   usecc2                !:  = 1.0 / ( ecc * ecc )
+   REAL(wp), PUBLIC ::   rhoco                 !: = rau0 * cw
+   REAL(wp), PUBLIC ::   sangvg, cangvg        !: sin and cos of the turning angle for ocean stress
+   REAL(wp), PUBLIC ::   pstarh                !: pstar / 2.0
 
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ahiu , ahiv   !: hor. diffusivity coeff. at ocean U- and V-points (m2/s)
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pahu , pahv   !: ice hor. eddy diffusivity coef. at ocean U- and V-points
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s         !: friction velocity
+
+   !!* Ice Rheology
+# if defined key_lim2_vp
+   !                                                      !!* VP rheology *
+   LOGICAL , PUBLIC ::   lk_lim2_vp = .TRUE.               !: Visco-Plactic reology flag 
+   !
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm   !: mean snow and ice thicknesses
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s                 !: friction velocity
+   !
+# else
+   !                                                      !!* EVP rheology *
+   LOGICAL , PUBLIC::   lk_lim2_vp = .FALSE.               !: Visco-Plactic reology flag 
+   !
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   stress1_i     !: first stress tensor element       
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   stress2_i     !: second stress tensor element
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   stress12_i    !: diagonal stress tensor element
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   delta_i       !: rheology delta factor (see Flato and Hibler 95) [s-1]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   divu_i        !: Divergence of the velocity field [s-1] -> limrhg.F90
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   shear_i       !: Shear of the velocity field [s-1] -> limrhg.F90
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   at_i          !: ice fraction
+   !
+   REAL(wp), PUBLIC, DIMENSION(:,:)    , POINTER :: vt_s ,vt_i    !: mean snow and ice thicknesses
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj), TARGET  :: hsnm , hicm   !: target vt_s,vt_i pointers 
+#endif
 
-   !!* diagnostic quantities
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature (Kelvin)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvosif       !: ice volume change at ice surface (only used for outputs)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvobif       !: ice volume change at ice bottom  (only used for outputs)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdvolif       !: Total   ice volume change (only used for outputs)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvonif       !: Lateral ice volume change (only used for outputs)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature [Kelvin]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tfu           !: Freezing/Melting point temperature of sea water at SSS
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicif         !: Ice thickness
@@ -71,17 +98,16 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmicif       !: Variation of ice mass
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qldif         !: heat balance of the lead (or of the open ocean)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to bring the ocean surface layer until its freezing 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to freeze the ocean surface layer
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdtcn         !: net downward heat flux from the ice to the ocean
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qdtcn         !: energy from the ice to the ocean point (at a factor 2)
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   thcm          !: part of the solar energy used in the lead heat budget
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fstric        !: Solar flux transmitted trough the ice
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ffltbif       !: Array linked with the max heat contained in brine pockets (?)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ffltbif       !: linked with the max heat contained in brine pockets (?)
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fscmbq        !: Linked with the solar flux below the ice (?)
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fsbbq         !: Also linked with the solar flux below the ice (?)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qfvbq         !: Array used to store energy in case of toral lateral ablation (?)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qfvbq         !: used to store energy in case of toral lateral ablation (?)
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   dmgwi         !: Variation of the mass of snow ice
-
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_ice, v_ice   !: two components of the ice   velocity at I-point (m/s)
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_oce, v_oce   !: two components of the ocean velocity at I-point (m/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_ice, v_ice  !: two components of the ice   velocity at I-point (m/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   u_oce, v_oce  !: two components of the ocean velocity at I-point (m/s)
 
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) ::   tbif  !: Temperature inside the ice/snow layer
@@ -102,4 +128,8 @@
 #endif
 
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE ice_2
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90	(revision 2528)
@@ -6,4 +6,5 @@
    !! History :   1.0  !  02-08  (G. Madec)  F90: Free form and modules
    !!             2.0  !  03-08  (C. Ethe)  add ice_run
+   !!             3.3  !  09-05  (G.Garric, C. Bricaud) addition of the lim2_evp case
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -15,14 +16,14 @@
    !!   ice_run_2        : Definition some run parameter for ice model
    !!----------------------------------------------------------------------
-   USE dom_oce
-   USE dom_ice_2
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE sbc_ice         ! surface boundary condition: ice
-   USE phycst          ! Define parameters for the routines
-   USE ice_2
-   USE limmsh_2
-   USE limistate_2
-   USE limrst_2   
-   USE in_out_manager
+   USE dom_oce          ! ocean domain
+   USE dom_ice_2        ! LIM2: ice domain
+   USE sbc_oce          ! surface boundary condition: ocean
+   USE sbc_ice          ! surface boundary condition: ice
+   USE phycst           ! Define parameters for the routines
+   USE ice_2            ! LIM2: ice variable
+   USE limmsh_2         ! LIM2: mesh
+   USE limistate_2      ! LIM2: initial state
+   USE limrst_2         ! LIM2: restart
+   USE in_out_manager   ! I/O manager
       
    IMPLICIT NONE
@@ -32,9 +33,8 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -43,25 +43,24 @@
       !!                  ***  ROUTINE ice_init_2  ***
       !!
-      !! ** purpose :   
+      !! ** purpose :   initialisation of LIM-2 domain and variables  
       !!----------------------------------------------------------------------
       !
-      ! Open the namelist file 
-      CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp )      
-      CALL ice_run_2                    !  read in namelist some run parameters
-                 
-      ! Louvain la Neuve Ice model
-      rdt_ice = nn_fsbc * rdttra(1)
-
-      CALL lim_msh_2                  ! ice mesh initialization
-     
-      ! Initial sea-ice state
-      IF( .NOT.ln_rstart ) THEN
-         CALL lim_istate_2            ! start from rest: sea-ice deduced from sst
-      ELSE
-         CALL lim_rst_read_2          ! start from a restart file
+      !                                ! Open the namelist file 
+      CALL ctl_opn( numnam_ice, 'namelist_ice', 'OLD', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp ) 
+      !    
+      CALL ice_run_2                   ! read in namelist some run parameters
+      !          
+      rdt_ice = nn_fsbc * rdttra(1)    ! sea-ice time step
+      numit   = nit000 - 1
+      !
+      CALL lim_msh_2                   ! ice mesh initialization
+      !
+      !                                ! Initial sea-ice state
+      IF( .NOT.ln_rstart ) THEN   ;   CALL lim_istate_2     ! start from rest: sea-ice deduced from sst
+      ELSE                        ;   CALL lim_rst_read_2   ! start from a restart file
       ENDIF
-      
-      tn_ice(:,:,1) = sist(:,:)         ! initialisation of ice temperature   
-      fr_i  (:,:) = 1.0 - frld(:,:)   ! initialisation of sea-ice fraction    
+      !
+      tn_ice(:,:,1) = sist(:,:)        ! initialisation of ice temperature   
+      fr_i  (:,:) = 1.0 - frld(:,:)    ! initialisation of sea-ice fraction    
       !
    END SUBROUTINE ice_init_2
@@ -82,8 +81,8 @@
       !!-------------------------------------------------------------------
       !                    
-      REWIND ( numnam_ice )                       ! Read Namelist namicerun 
-      READ   ( numnam_ice , namicerun )
-
-      IF(lwp) THEN
+      REWIND( numnam_ice )                      ! Read Namelist namicerun 
+      READ  ( numnam_ice , namicerun )
+      !
+      IF(lwp) THEN                              ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'ice_run : ice share parameters for dynamics/advection/thermo of sea-ice'
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limadv_2.F90	(revision 2528)
@@ -37,7 +37,7 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2009) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limdia_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limdia_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limdia_2.F90	(revision 2528)
@@ -57,7 +57,7 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limdmp_2.F90	(revision 2528)
@@ -2,169 +2,170 @@
    !!======================================================================
    !!                       ***  MODULE limdmp_2   ***
-   !!  Ice model : restoring Ice thickness and Fraction leads
+   !!  LIM-2 ice model : restoring Ice thickness and Fraction leads
    !!======================================================================
-   !! History :   2.0  !  04-04 (S. Theetten) Original code
+   !! History :   2.0  !  2004-04 (S. Theetten) Original code
+   !!             3.3  !  2010-06 (J.-M. Molines) use of fldread
    !!----------------------------------------------------------------------
-#if defined key_lim2   &&   defined key_tradmp
+#if defined key_lim2
    !!----------------------------------------------------------------------
-   !!   'key_lim2'  AND                               LIM 2.0 sea-ice model
-   !!   'key_tradmp'                                                Damping
-   !!----------------------------------------------------------------------
+   !!   'key_lim2'                                    LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_dmp_2      : ice model damping
    !!----------------------------------------------------------------------
    USE in_out_manager  ! I/O manager
-   USE phycst          ! physical constants
-   USE ice_2
-   USE tradmp
-   USE dom_oce
-   USE oce
-   USE iom
+   USE ice_2           ! ice variables 
+   USE sbc_oce, ONLY : nn_fsbc ! for fldread
+   USE dom_oce         ! for mi0; mi1 etc ...
+   USE fldread         ! read input fields
    
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   lim_dmp_2     ! called by ice_step_2
+   PUBLIC   lim_dmp_2     ! called by sbc_ice_lim2
+
+   REAL(wp), DIMENSION(:,:,:), ALLOCATABLE ::   resto_ice   ! restoring coeff. on ICE   [s-1]
+
+   INTEGER, PARAMETER :: jp_hicif = 1 , jp_frld = 2
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_icedmp    ! structure of ice damping input
    
-   INTEGER                        ::   nice1, nice2,  &  ! first and second record used
-      &                                inumice_dmp       ! logical unit for ice variables (damping)
-   REAL(wp), DIMENSION(jpi,jpj)   ::   hicif_dta  ,   &  ! ice thickness at a given time
-      &                                frld_dta          ! fraction lead at a given time
-   REAL(wp), DIMENSION(jpi,jpj,2) ::   hicif_data ,   &  ! ice thickness data at two consecutive times
-      &                                frld_data         ! fraction lead data at two consecutive times
-
    !! * Substitution
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0 , UCL-LOCEAN-IPSL  (2006)
+   !! NEMO/LIM 3.3 , UCL-NEMO-consortium (2010) 
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
-   SUBROUTINE lim_dmp_2(kt)
+   SUBROUTINE lim_dmp_2( kt )
       !!-------------------------------------------------------------------
-      !!                   ***  ROUTINE lim_dmp_2 ***
+      !!                   ***  ROUTINE lim_dmp_2  ***
       !!
-      !! ** purpose : ice model damping : restoring ice thickness and 
-      !!              fraction leads
+      !! ** purpose : ice model damping : restoring ice thickness and fraction leads
       !!
-      !! ** method  : the key_tradmp must be used to compute resto(:,:) coef.
+      !! ** method  : the key_tradmp must be used to compute resto(:,:,1) coef.
       !!---------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt     ! ocean time-step
+      INTEGER, INTENT(in) ::   kt   ! ocean time-step
       !
-      INTEGER             ::   ji, jj         ! dummy loop indices
+      INTEGER  ::   ji, jj         ! dummy loop indices
+      REAL(wp) ::   zfrld, zhice   ! local scalars
       !!---------------------------------------------------------------------
       !
-      CALL dta_lim_2( kt )
-
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            hicif(ji,jj) = hicif(ji,jj) - rdt_ice * resto(ji,jj,1) * ( hicif(ji,jj) - hicif_dta(ji,jj) )
-            frld(ji,jj)  = frld (ji,jj) - rdt_ice * resto(ji,jj,1) * ( frld(ji,jj)  - frld_dta (ji,jj) )  
-         END DO
-      END DO
+      IF (kt == nit000)  THEN 
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'lim_dmp_2 : Ice thickness and ice concentration restoring'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~'
+         !
+         ! ice_resto_init create resto_ice (in 1/s) for restoring ice parameters near open boundaries.
+         ! Double check this routine to verify if it corresponds to your config
+         CALL lim_dmp_init
+      ENDIF
+      !
+      IF( ln_limdmp ) THEN   ! ice restoring in this case
+         !
+         CALL fld_read( kt, nn_fsbc, sf_icedmp )
+         !
+!CDIR COLLAPSE
+         hicif(:,:) = MAX( 0._wp,                     &        ! h >= 0         avoid spurious out of physical range
+            &         hicif(:,:) - rdt_ice * resto_ice(:,:,1) * ( hicif(:,:) - sf_icedmp(jp_hicif)%fnow(:,:,1) )  ) 
+!CDIR COLLAPSE
+         hicif(:,:) = MAX( 0._wp, MIN( 1._wp,         &        ! 0<= frld<=1    values which blow the run up
+            &         frld (:,:) - rdt_ice * resto_ice(:,:,1) * ( frld (:,:) - sf_icedmp(jp_frld )%fnow(:,:,1) )  )  )
+         !
+      ENDIF
       !
    END SUBROUTINE lim_dmp_2
 
 
-   SUBROUTINE dta_lim_2( kt ) 
+   SUBROUTINE lim_dmp_init
       !!----------------------------------------------------------------------
-      !!                   ***  ROUTINE dta_lim_2  ***
+      !!                   ***  ROUTINE lim_dmp_init  ***
       !!
-      !! ** Purpose :   Reads monthly ice thickness and fraction lead  data
+      !! ** Purpose :   Initialization for the ice thickness and concentration 
+      !!                restoring
+      !!              restoring will be used. It is used to mimic ice open
+      !!              boundaries.
       !!
-      !! ** Method  :   Read on unit numicedt the interpolated ice variable
-      !!      onto the model grid.
-      !!      Data begin at january.
-      !!      The value is centered at the middle of month.
-      !!      In the opa model, kt=1 agree with january 1.
-      !!      At each time step, a linear interpolation is applied between
-      !!      two monthly values.
+      !! ** Method  :  ?????
       !!      
-      !! ** Action  :   define hicif_dta and frld_dta arrays at time-step kt
+      !! ** Action  :   define resto_ice(:,:,1)
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt     ! ocean time-step
+      INTEGER  :: ji, jj, jk       ! dummy loop indices
+      INTEGER  :: irelax, ierror   ! error flag for allocation
       !
-      INTEGER  ::   imois, iman, i15          ! temporary integers
-      REAL(wp) ::   zxy
+      REAL(wp) :: zdmpmax, zdmpmin, zfactor, zreltim ! temporary scalar
+      !
+      CHARACTER(len=100)           ::   cn_dir       ! Root directory for location of ssr files
+      TYPE(FLD_N), DIMENSION (2)   ::   sl_icedmp    ! informations about the icedmp  field to be read
+      TYPE(FLD_N)                  ::   sn_hicif     ! 
+      TYPE(FLD_N)                  ::   sn_frld      ! 
+      NAMELIST/namice_dmp/ cn_dir, ln_limdmp, sn_hicif, sn_frld
       !!----------------------------------------------------------------------
+      !
+      ! 1)  initialize fld read structure for input data 
+      !     --------------------------------------------
+      ln_limdmp = .false.                 !* set file information (default values)
+      cn_dir    = './'
+      ! (NB: frequency positive => hours, negative => months)
+      !                !    file     ! frequency ! variable ! time intep !  clim   ! 'yearly' or ! weights  ! rotation !
+      !                !    name     !  (hours)  !  name    !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
+      sn_hicif = FLD_N( 'ice_damping ', -1       , 'hicif'  ,  .true.    , .true.  ,  'yearly'   ,  ''      ,  ''      )
+      sn_frld  = FLD_N( 'ice_damping ', -1       , 'frld'   ,  .true.    , .true.  ,  'yearly'   ,  ''      ,  ''      )
 
-      ! 0. Initialization
-      ! -----------------
-      iman  = INT( raamo )
-!!! better but change the results     i15 = INT( 2*FLOAT( nday ) / ( FLOAT( nobis(nmonth) ) + 0.5 ) )
-      i15   = nday / 16
-      imois = nmonth + i15 - 1
-      IF( imois == 0 ) imois = iman
-      
-      ! 1. First call kt=nit000: Initialization and Open
-      ! -----------------------
-      IF( kt == nit000 ) THEN
-         nice1 = 0
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'dtalim : Ice thickness and lead fraction  monthly fields'
-         IF(lwp) WRITE(numout,*) '~~~~~~'
-         IF(lwp) WRITE(numout,*) '             NetCDF FORMAT'
-         IF(lwp) WRITE(numout,*)
-         ! open file
-         CALL iom_open( 'ice_damping_ATL4.nc', inumice_dmp )
+      REWIND( numnam_ice )                !* read in namelist_ice namicedmp
+      READ  ( numnam_ice, namice_dmp )
+      !
+      IF ( lwp ) THEN                     !* control print
+         WRITE (numout,*)'     lim_dmp_init : lim_dmp initialization ' 
+         WRITE (numout,*)'       Namelist namicedmp read '
+         WRITE (numout,*)'         Ice restoring (T) or not (F) ln_limdmp =', ln_limdmp 
+         WRITE (numout,*)
+         WRITE (numout,*)'     CAUTION : here hard coded ice restoring along northern and southern boundaries'
+         WRITE (numout,*)'               adapt the lim_dmp_init routine to your needs'
       ENDIF
 
+      ! 2)  initialise resto_ice    ==>  config dependant !
+      !     --------------------         ++++++++++++++++
+      !
+      IF( ln_limdmp ) THEN                !* ice restoring is used, follow initialization
+         ! 
+         sl_icedmp ( jp_hicif ) = sn_hicif
+         sl_icedmp ( jp_frld  ) = sn_frld
+         ALLOCATE ( sf_icedmp (2) , resto_ice(jpi,jpj,1), STAT=ierror )
+         IF( ierror > 0 ) THEN
+            CALL ctl_stop( 'lim_dmp_init: unable to allocate sf_icedmp structure or resto_ice array' )   ;   RETURN
+         ENDIF
+         ALLOCATE( sf_icedmp(jp_hicif)%fnow(jpi,jpj,1) , sf_icedmp(jp_hicif)%fdta(jpi,jpj,1,2) )
+         ALLOCATE( sf_icedmp(jp_frld )%fnow(jpi,jpj,1) , sf_icedmp(jp_frld )%fdta(jpi,jpj,1,2) )
+         !                         ! fill sf_icedmp with sn_icedmp and control print
+         CALL fld_fill( sf_icedmp, sl_icedmp, cn_dir, 'lim_dmp_init', 'Ice  restoring input data', 'namicedmp' )
+      
+         resto_ice(:,:,:) = 0._wp
+         !      Re-calculate the North and South boundary restoring term
+         !      because those boundaries may change with the prescribed zoom area.
+         !
+         irelax  = 16                     ! width of buffer zone with respect to close boundary
+         zdmpmax = 10._wp                 ! max restoring time scale  (days) (low restoring)
+         zdmpmin = rdt_ice / 86400._wp    ! min restoring time scale  (days) (high restoring)
+         !                                ! days / grid-point
+         zfactor = ( zdmpmax - zdmpmin ) / REAL( irelax, wp )
 
-      ! 2. Read monthly file
-      ! --------------------
-      IF( ( kt == nit000 ) .OR. imois /= nice1 ) THEN
-         !
-         ! Calendar computation
-         nice1 = imois        ! first file record used 
-         nice2 = nice1 + 1    ! last  file record used
-         nice1 = MOD( nice1, iman )
-         nice2 = MOD( nice2, iman )
-         IF( nice1 == 0 )   nice1 = iman
-         IF( nice2 == 0 )   nice2 = iman
-         IF(lwp) WRITE(numout,*) 'first record file used nice1 ', nice1
-         IF(lwp) WRITE(numout,*) 'last  record file used nice2 ', nice2
-         
-         ! Read monthly ice thickness Levitus 
-         CALL iom_get( inumice_dmp, jpdom_data, 'iicethic', hicif_data(:,:,1), nice1 ) 
-         CALL iom_get( inumice_dmp, jpdom_data, 'iicethic', hicif_data(:,:,2), nice2 ) 
-         
-         ! Read monthly ice thickness Levitus 
-         CALL iom_get( inumice_dmp, jpdom_data, 'ileadfra', frld_data(:,:,1), nice1 ) 
-         CALL iom_get( inumice_dmp, jpdom_data, 'ileadfra', frld_data(:,:,2), nice2 ) 
-         
-         ! The fraction lead read in the file is in fact the 
-         ! ice concentration which is 1 - the fraction lead
-         frld_data = 1 - frld_data          
-         
-         IF(lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) ' Ice thickness month ', nice1,' and ', nice2
-            WRITE(numout,*)
-            WRITE(numout,*) ' Ice thickness month = ', nice1
-            CALL prihre( hicif_data(1,1,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
-            WRITE(numout,*)
-            WRITE(numout,*) ' Fraction lead months ', nice1,' and ', nice2
-            WRITE(numout,*)
-            WRITE(numout,*) ' Fraction lead month = ', nice1
-            CALL prihre( frld_data(1,1,1), jpi, jpj, 1, jpi, 20, 1, jpj, 20, 1., numout )
-         ENDIF
+         !    South boundary restoring term
+         ! REM: if there is no ice in the model and in the data, 
+         !      no restoring even with non zero resto_ice
+         DO jj = mj0(jpjzoom - 1 + 1), mj1(jpjzoom -1 + irelax)
+            zreltim = zdmpmin + zfactor * ( mjg(jj) - jpjzoom + 1 )
+            resto_ice(:,jj,:) = 1._wp / ( zreltim * 86400._wp )
+         END DO
 
-         CALL FLUSH(numout)
-
+         ! North boundary restoring term
+         DO jj =  mj0(jpjzoom -1 + jpjglo - irelax), mj1(jpjzoom - 1 + jpjglo)
+            zreltim = zdmpmin + zfactor * (jpjglo - ( mjg(jj) - jpjzoom + 1 ))
+            resto_ice(:,jj,:) = 1.e0 / ( zreltim * 86400 )
+         END DO
       ENDIF
-        
-      ! 3. At every time step compute ice thickness and fraction lead data
-      ! ------------------------------------------------------------------
-      zxy = FLOAT( nday + 15 - 30 * i15 ) / 30.
-      hicif_dta(:,:) = (1.-zxy) * hicif_data(:,:,1) + zxy * hicif_data(:,:,2)
-      frld_dta (:,:) = (1.-zxy) * frld_data (:,:,1) + zxy * frld_data (:,:,2)
-
-      IF( kt == nitend )   CALL iom_close( inumice_dmp )
       !
-   END SUBROUTINE dta_lim_2
-
+   END SUBROUTINE lim_dmp_init
+   
 #else
    !!----------------------------------------------------------------------
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 2528)
@@ -4,8 +4,9 @@
    !!   Sea-Ice dynamics :  
    !!======================================================================
-   !! History :   1.0  !  01-04  (LIM)  Original code
-   !!             2.0  !  02-08  (C. Ethe, G. Madec)  F90, mpp
-   !!             2.0  !  03-08  (C. Ethe) add lim_dyn_init
-   !!             2.0  !  06-07  (G. Madec)  Surface module
+   !! History :  1.0  ! 2001-04  (LIM)  Original code
+   !!            2.0  ! 2002-08  (C. Ethe, G. Madec)  F90, mpp
+   !!            2.0  ! 2003-08  (C. Ethe) add lim_dyn_init
+   !!            2.0  ! 2006-07  (G. Madec)  Surface module
+   !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
    !!---------------------------------------------------------------------
 #if defined key_lim2
@@ -16,32 +17,30 @@
    !!    lim_dyn_init_2 : initialization and namelist read
    !!----------------------------------------------------------------------
-   USE dom_oce        ! ocean space and time domain
-   USE sbc_oce        !
-   USE phycst         !
-   USE ice_2          !
-   USE dom_ice_2      !
-   USE limistate_2    !
-   USE limrhg_2       ! ice rheology
-
-   USE lbclnk         !
-   USE lib_mpp        !
-   USE in_out_manager ! I/O manager
-   USE prtctl         ! Print control
+   USE dom_oce          ! ocean space and time domain
+   USE sbc_oce          ! ocean surface boundary condition
+   USE phycst           ! physical constant
+   USE ice_2            ! LIM-2: ice variables
+   USE sbc_ice          ! Surface boundary condition: sea-ice fields
+   USE dom_ice_2        ! LIM-2: ice domain
+   USE limistate_2      ! LIM-2: initial state
+   USE limrhg_2         ! LIM-2: VP  ice rheology
+   USE limrhg           ! LIM  : EVP ice rheology
+   USE lbclnk           ! lateral boundary condition - MPP link
+   USE lib_mpp          ! MPP library
+   USE in_out_manager   ! I/O manager
+   USE prtctl           ! Print control
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   lim_dyn_2 ! routine called by sbc_ice_lim
-
-   !! * Module variables
-   REAL(wp)  ::  rone    = 1.e0   ! constant value
-
+   PUBLIC   lim_dyn_2   ! routine called by sbc_ice_lim
+
+   !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -83,9 +82,11 @@
          ! ---------------------------------------------------
          
-         IF( lk_mpp .OR. nbit_cmp == 1 ) THEN                    ! mpp: compute over the whole domain
+         IF( lk_mpp .OR. lk_mpp_rep ) THEN                    ! mpp: compute over the whole domain
             i_j1 = 1   
             i_jpj = jpj
             IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
-            CALL lim_rhg_2( i_j1, i_jpj )
+            IF( lk_lim2_vp )   THEN   ;   CALL lim_rhg_2( i_j1, i_jpj )             !  VP rheology
+            ELSE                      ;   CALL lim_rhg  ( i_j1, i_jpj )             ! EVP rheology
+            ENDIF
             !
          ELSE                                 ! optimization of the computational area
@@ -105,9 +106,13 @@
                   i_j1 = i_j1 + 1
                END DO
-               i_j1 = MAX( 1, i_j1-1 )
-               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
-               ! 
-               CALL lim_rhg_2( i_j1, i_jpj )
-               ! 
+               IF( lk_lim2_vp )   THEN             ! VP  rheology
+                  i_j1 = MAX( 1, i_j1-1 )
+                  CALL lim_rhg_2( i_j1, i_jpj )
+               ELSE                                ! EVP rheology
+                  i_j1 = MAX( 1, i_j1-2 )
+                  CALL lim_rhg( i_j1, i_jpj )
+               ENDIF
+               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, 'ij_jpj = ', i_jpj
+               !
                ! Southern hemisphere
                i_j1  =  1 
@@ -116,9 +121,13 @@
                   i_jpj = i_jpj - 1
                END DO
-               i_jpj = MIN( jpj, i_jpj+2 )
-               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
-               ! 
-               CALL lim_rhg_2( i_j1, i_jpj )
-               ! 
+               IF( lk_lim2_vp )   THEN             ! VP  rheology
+                  i_jpj = MIN( jpj, i_jpj+2 )
+                  CALL lim_rhg_2( i_j1, i_jpj )
+               ELSE                                ! EVP rheology
+                  i_jpj = MIN( jpj, i_jpj+1 )
+                  CALL lim_rhg( i_j1, i_jpj )
+               ENDIF
+               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, 'ij_jpj = ', i_jpj
+               !
             ELSE                                 ! local domain extends over one hemisphere only
                !                                 ! Rheology is computed only over the ice cover
@@ -134,9 +143,15 @@
                   i_jpj = i_jpj - 1
                END DO
-               i_jpj = MIN( jpj, i_jpj+2)
-    
+               i_jpj = MIN( jpj, i_jpj+2 )
+               ! 
+               IF( lk_lim2_vp )   THEN             ! VP  rheology
+                  i_jpj = MIN( jpj, i_jpj+2 )
+                  CALL lim_rhg_2( i_j1, i_jpj )                !  VP rheology
+               ELSE                                ! EVP rheology
+                  i_j1  = MAX( 1  , i_j1-2  )
+                  i_jpj = MIN( jpj, i_jpj+1 )
+                  CALL lim_rhg  ( i_j1, i_jpj )                ! EVP rheology
+               ENDIF
                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
-               ! 
-               CALL lim_rhg_2( i_j1, i_jpj )
                !
             ENDIF
@@ -148,18 +163,23 @@
          ! computation of friction velocity
          ! --------------------------------
-         ! ice-ocean velocity at U & V-points (u_ice v_ice at I-point ; ssu_m, ssv_m at U- & V-points)
-         
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1   ! NO vector opt.
-               zu_io(ji,jj) = 0.5 * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj)
-               zv_io(ji,jj) = 0.5 * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj)
-            END DO
-         END DO
+         SELECT CASE( cp_ice_msh )           ! ice-ocean relative velocity at u- & v-pts
+         CASE( 'C' )                               ! EVP : C-grid ice dynamics
+            zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)           ! ice-ocean & ice velocity at ocean velocity points
+            zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
+         CASE( 'I' )                               ! VP  : B-grid ice dynamics (I-point) 
+            DO jj = 1, jpjm1                               ! u_ice v_ice at I-point ; ssu_m, ssv_m at U- & V-points
+               DO ji = 1, jpim1   ! NO vector opt.         !
+                  zu_io(ji,jj) = 0.5_wp * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj)
+                  zv_io(ji,jj) = 0.5_wp * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj)
+               END DO
+            END DO
+         END SELECT
+
          ! frictional velocity at T-point
+         zcoef = 0.5_wp * cw
          DO jj = 2, jpjm1
             DO ji = 2, jpim1   ! NO vector opt. because of zu_io
-               ust2s(ji,jj) = 0.5 * cw                                                          &
-                  &         * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
-                  &            + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
+               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
+                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
             END DO
          END DO
@@ -170,6 +190,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ust2s(ji,jj) = zcoef * tms(ji,jj) * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
-                  &                                     + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) )
+               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
+                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1)   ) * tms(ji,jj)
             END DO
          END DO
@@ -180,5 +200,5 @@
       !
       IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
-
+      !
    END SUBROUTINE lim_dyn_2
 
@@ -198,5 +218,6 @@
       NAMELIST/namicedyn/ epsd, alpha,     &
          &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
-         &                c_rhg, etamn, creepl, ecc, ahi0
+         &                c_rhg, etamn, creepl, ecc, ahi0,                  &
+         &                nevp, telast,alphaevp
       !!-------------------------------------------------------------------
 
@@ -223,4 +244,13 @@
          WRITE(numout,*) '       eccentricity of the elliptical yield curve       ecc    = ', ecc
          WRITE(numout,*) '       horizontal diffusivity coeff. for sea-ice        ahi0   = ', ahi0
+         WRITE(numout,*) '       number of iterations for subcycling nevp   = ', nevp
+         WRITE(numout,*) '       timescale for elastic waves telast = ', telast
+         WRITE(numout,*) '       coefficient for the solution of int. stresses alphaevp = ', alphaevp
+      ENDIF
+      !
+      IF( angvg /= 0._wp .AND. .NOT.lk_lim2_vp ) THEN
+         CALL ctl_warn( 'lim_dyn_init_2: turning angle for oceanic stress not properly coded for EVP ',   &
+            &           '(see limsbc_2 module). We force  angvg = 0._wp'  )
+         angvg = 0._wp
       ENDIF
 
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limhdf_2.F90	(revision 2528)
@@ -32,7 +32,7 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limistate_2.F90	(revision 2528)
@@ -46,7 +46,7 @@
    REAL(wp) ::   zone      = 1.e0     ! constant value = 1
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 2528)
@@ -4,4 +4,8 @@
    !! LIM 2.0 ice model :   definition of the ice mesh parameters
    !!======================================================================
+   !! History :   -   ! 2001-04 (LIM) original code
+   !!            1.0  ! 2002-08 (C. Ethe, G. Madec) F90, module
+   !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
+   !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
@@ -10,5 +14,4 @@
    !!   lim_msh_2   : definition of the ice mesh
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE phycst
    USE dom_oce
@@ -20,13 +23,11 @@
    PRIVATE
 
-   !! * Accessibility
    PUBLIC lim_msh_2      ! routine called by ice_ini_2.F90
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -43,19 +44,13 @@
       !! 
       !! ** Refer.  : Deleersnijder et al. Ocean Modelling 100, 7-10 
-      !!
-      !! ** History :
-      !!         original    : 01-04 (LIM)
-      !!         addition    : 02-08 (C. Ethe, G. Madec)
       !!--------------------------------------------------------------------- 
-      !! * Local variables
       INTEGER :: ji, jj      ! dummy loop indices
-
-      REAL(wp), DIMENSION(jpi,jpj) ::  &
-         zd2d1 , zd1d2       ! Derivative of zh2 (resp. zh1) in the x direction
-         !                   ! (resp. y direction) (defined at the center)
-      REAL(wp) ::         &
-         zh1p  , zh2p   , &  ! Idem zh1, zh2 for the bottom left corner of the grid
-         zd2d1p, zd1d2p , &  ! Idem zd2d1, zd1d2 for the bottom left corner of the grid
-         zusden, zusden2     ! temporary scalars
+      REAL(wp) ::   zusden   ! local scalars
+#if defined key_lim2_vp
+      REAL(wp) ::   zusden2           ! local scalars
+      REAL(wp) ::   zh1p  , zh2p      !   -      -
+      REAL(wp) ::   zd2d1p, zd1d2p    !   -      -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zd2d1 , zd1d2   ! 2D workspace
+#endif
       !!---------------------------------------------------------------------
 
@@ -76,5 +71,5 @@
       njeqm1 = njeq - 1 
 
-      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor
+      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor at T-point
  
 !i    DO jj = 1, jpj
@@ -115,11 +110,17 @@
       !-------------------
 !!ibug ???
-      akappa(:,:,:,:) = 0.e0
       wght(:,:,:,:) = 0.e0
+      tmu(:,:)      = 0.e0
+#if defined key_lim2_vp 
+      akappa(:,:,:,:)     = 0.e0
       alambd(:,:,:,:,:,:) = 0.e0
-      tmu(:,:) = 0.e0
+#else
+      tmv(:,:) = 0.e0
+      tmf(:,:) = 0.e0
+#endif
 !!i
       
-      
+
+#if defined key_lim2_vp      
       ! metric coefficients for sea ice dynamic
       !----------------------------------------
@@ -155,8 +156,25 @@
       CALL lbc_lnk( wght(:,:,2,1), 'I', 1. )      ! but it is never used
       CALL lbc_lnk( wght(:,:,2,2), 'I', 1. )
+#else
+      ! metric coefficients for sea ice dynamic (EVP rheology)
+      !----------------------------------------
+      DO jj = 1, jpjm1                                       ! weights (wght) at F-points
+         DO ji = 1, jpim1
+            zusden = 1. / (  ( e1t(ji+1,jj  ) + e1t(ji,jj) )   &
+               &           * ( e2t(ji  ,jj+1) + e2t(ji,jj) ) ) 
+            wght(ji,jj,1,1) = zusden * e1t(ji+1,jj) * e2t(ji,jj+1)
+            wght(ji,jj,1,2) = zusden * e1t(ji+1,jj) * e2t(ji,jj  )
+            wght(ji,jj,2,1) = zusden * e1t(ji  ,jj) * e2t(ji,jj+1)
+            wght(ji,jj,2,2) = zusden * e1t(ji  ,jj) * e2t(ji,jj  )
+         END DO
+      END DO
+      CALL lbc_lnk( wght(:,:,1,1), 'F', 1. )   ;   CALL lbc_lnk( wght(:,:,1,2),'F', 1. )       ! lateral boundary cond.   
+      CALL lbc_lnk( wght(:,:,2,1), 'F', 1. )   ;   CALL lbc_lnk( wght(:,:,2,2),'F', 1. )
+#endif
     
       ! Coefficients for divergence of the stress tensor
       !-------------------------------------------------
 
+#if defined key_lim2_vp
       DO jj = 2, jpj
          DO ji = 2, jpi   ! NO vector opt.
@@ -226,4 +244,5 @@
       CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )      !
       CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )      !
+#endif
             
 
@@ -233,4 +252,6 @@
       tms(:,:) = tmask(:,:,1)      ! ice T-point  : use surface tmask
 
+#if defined key_lim2_vp
+      ! VP rheology : ice velocity point is I-point
 !i here we can use umask with a i and j shift of -1,-1
       tmu(:,1) = 0.e0
@@ -241,11 +262,17 @@
          END DO
       END DO
-      
-      !--lateral boundary conditions    
-      CALL lbc_lnk( tmu(:,:), 'I', 1. )
-      
+      CALL lbc_lnk( tmu(:,:), 'I', 1. )      !--lateral boundary conditions    
+#else
+      ! EVP rheology : ice velocity point are U- & V-points ; ice vorticity
+      ! point is F-point
+      tmu(:,:) = umask(:,:,1)
+      tmv(:,:) = vmask(:,:,1)
+      tmf(:,:) = 0.e0                        ! used of fmask except its special value along the coast (rn_shlat)
+      WHERE( fmask(:,:,1) == 1.e0 )   tmf(:,:) = 1.e0
+#endif
+      !
       ! unmasked and masked area of T-grid cell
       area(:,:) = e1t(:,:) * e2t(:,:)
-      
+      !
    END SUBROUTINE lim_msh_2
 
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limrhg.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limrhg.F90	(revision 2528)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limrhg.F90	(revision 2528)
@@ -0,0 +1,1 @@
+link ../LIM_SRC_3/limrhg.F90
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 2528)
@@ -4,14 +4,15 @@
    !!   Ice rheology :  performs sea ice rheology
    !!======================================================================
-   !! History :  0.0  !  93-12  (M.A. Morales Maqueda.)  Original code
-   !!            1.0  !  94-12  (H. Goosse) 
-   !!            2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp
-   !!            " "  !  06-08  (G. Madec)  surface module, ice-stress at I-point
-   !!            " "  !  09-09  (G. Madec)  Huge verctor optimisation
-   !!----------------------------------------------------------------------
-#if defined key_lim2
-   !!----------------------------------------------------------------------
-   !!   'key_lim2'                                    LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
+   !! History :  0.0  !  1993-12  (M.A. Morales Maqueda.)  Original code
+   !!            1.0  !  1994-12  (H. Goosse) 
+   !!            2.0  !  2003-12  (C. Ethe, G. Madec)  F90, mpp
+   !!             -   !  2006-08  (G. Madec)  surface module, ice-stress at I-point
+   !!             -   !  2009-09  (G. Madec)  Huge verctor optimisation
+   !!            3.3  !  2009-05  (G.Garric, C. Bricaud) addition of the lim2_evp case
+   !!----------------------------------------------------------------------
+#if defined   key_lim2   &&   defined key_lim2_vp
+   !!----------------------------------------------------------------------
+   !!   'key_lim2'                AND                   LIM-2 sea-ice model
+   !!   'key_lim2_vp'                                       VP ice rheology
    !!----------------------------------------------------------------------
    !!   lim_rhg_2   : computes ice velocities
@@ -21,8 +22,8 @@
    USE sbc_oce        ! surface boundary condition: ocean variables
    USE sbc_ice        ! surface boundary condition: ice variables
-   USE dom_ice_2      ! domaine: ice variables
+   USE dom_ice_2      ! LIM2: ice domain
    USE phycst         ! physical constant
-   USE ice_2          ! ice variables
-   USE lbclnk         ! lateral boundary condition
+   USE ice_2          ! LIM2: ice variables
+   USE lbclnk         ! lateral boundary condition - MPP exchanges
    USE lib_mpp        ! MPP library
    USE in_out_manager ! I/O manager
@@ -34,15 +35,14 @@
    PUBLIC   lim_rhg_2 ! routine called by lim_dyn
 
-   REAL(wp) ::   rzero   = 0.e0   ! constant value: zero
-   REAL(wp) ::   rone    = 1.e0   !            and  one
+   REAL(wp) ::   rzero   = 0._wp   ! constant value: zero
+   REAL(wp) ::   rone    = 1._wp   !            and  one
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -65,9 +65,9 @@
       INTEGER ::   iter, jter          ! temporary integers
       CHARACTER (len=50) ::   charout
-      REAL(wp) ::   ze11  , ze12  , ze22  , ze21               ! temporary scalars
-      REAL(wp) ::   zt11  , zt12  , zt21  , zt22               !    "         "
-      REAL(wp) ::   zvis11, zvis21, zvis12, zvis22             !    "         "
-      REAL(wp) ::   zgphsx, ztagnx, zgsshx, zunw, zur, zusw    !    "         "
-      REAL(wp) ::   zgphsy, ztagny, zgsshy, zvnw, zvr          !    "         "
+      REAL(wp) ::   ze11  , ze12  , ze22  , ze21               ! local scalars
+      REAL(wp) ::   zt11  , zt12  , zt21  , zt22               !   -      -
+      REAL(wp) ::   zvis11, zvis21, zvis12, zvis22             !   -      -
+      REAL(wp) ::   zgphsx, ztagnx, zgsshx, zunw, zur, zusw    !   -      -
+      REAL(wp) ::   zgphsy, ztagny, zgsshy, zvnw, zvr          !   -      -
       REAL(wp) ::   zresm,  za, zac, zmod
       REAL(wp) ::   zmpzas, zstms, zindu, zusdtp, zmassdt, zcorlal
@@ -89,45 +89,37 @@
       REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zzfrld, zztms
       REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zi1, zi2, zmasst, zpresh
-
       !!-------------------------------------------------------------------
-
-!!bug
-!!    u_oce(:,:) = 0.e0
-!!    v_oce(:,:) = 0.e0
-!!    write(*,*) 'rhg min, max u & v', maxval(u_oce), minval(u_oce), maxval(v_oce), minval(v_oce)
-!!bug
       
       !  Store initial velocities
       !  ----------------
-      zztms(:,0    ) = 0.e0       ;    zzfrld(:,0    ) = 0.e0
-      zztms(:,jpj+1) = 0.e0       ;    zzfrld(:,jpj+1) = 0.e0
-      zu0(:,0    ) = 0.e0         ;    zv0(:,0    ) = 0.e0
-      zu0(:,jpj+1) = 0.e0         ;    zv0(:,jpj+1) = 0.e0
-      zztms(:,1:jpj) = tms(:,:)   ;    zzfrld(:,1:jpj) = frld(:,:)
-      zu0(:,1:jpj) = u_ice(:,:)   ;    zv0(:,1:jpj) = v_ice(:,:)
-
-      zu_a(:,:)    = zu0(:,:)     ;   zv_a(:,:) = zv0(:,:)
-      zu_n(:,:)    = zu0(:,:)     ;   zv_n(:,:) = zv0(:,:)
+      zztms(:,0    ) = 0._wp        ;   zzfrld(:,0    ) = 0._wp
+      zztms(:,jpj+1) = 0._wp        ;   zzfrld(:,jpj+1) = 0._wp
+      zu0  (:,0    ) = 0._wp        ;   zv0   (:,0    ) = 0._wp
+      zu0  (:,jpj+1) = 0._wp        ;   zv0   (:,jpj+1) = 0._wp
+      zztms(:,1:jpj) = tms  (:,:)   ;   zzfrld(:,1:jpj) = frld (:,:)
+      zu0  (:,1:jpj) = u_ice(:,:)   ;   zv0   (:,1:jpj) = v_ice(:,:)
+      zu_a (:, :   ) = zu0  (:,:)   ;   zv_a  (:, :   ) = zv0  (:,:)
+      zu_n (:, :   ) = zu0  (:,:)   ;   zv_n  (:, :   ) = zv0  (:,:)
 
 !i
-      zi1   (:,:) = 0.e0
-      zi2   (:,:) = 0.e0
-      zpresh(:,:) = 0.e0
-      zmasst(:,:) = 0.e0
+      zi1   (:,:) = 0._wp
+      zi2   (:,:) = 0._wp
+      zpresh(:,:) = 0._wp
+      zmasst(:,:) = 0._wp
 !i
 !!gm violant
-      zfrld(:,:) =0.e0
-      zcorl(:,:) =0.e0
-      zmass(:,:) =0.e0
-      za1ct(:,:) =0.e0
-      za2ct(:,:) =0.e0
+      zfrld(:,:) =0._wp
+      zcorl(:,:) =0._wp
+      zmass(:,:) =0._wp
+      za1ct(:,:) =0._wp
+      za2ct(:,:) =0._wp
 !!gm end
 
-      zviszeta(:,:) = 0.e0
-      zviseta (:,:) = 0.e0
-
-!i    zviszeta(:,0    ) = 0.e0    ;    zviseta(:,0    ) = 0.e0
-!i    zviszeta(:,jpj  ) = 0.e0    ;    zviseta(:,jpj  ) = 0.e0
-!i    zviszeta(:,jpj+1) = 0.e0    ;    zviseta(:,jpj+1) = 0.e0
+      zviszeta(:,:) = 0._wp
+      zviseta (:,:) = 0._wp
+
+!i    zviszeta(:,0    ) = 0._wp    ;    zviseta(:,0    ) = 0._wp
+!i    zviszeta(:,jpj  ) = 0._wp    ;    zviseta(:,jpj  ) = 0._wp
+!i    zviszeta(:,jpj+1) = 0._wp    ;    zviseta(:,jpj+1) = 0._wp
 
 
@@ -141,11 +133,11 @@
          DO ji = 1 , jpi
             ! only the sinus changes its sign with the hemisphere
-            zsang(ji,jj)  = SIGN( 1.e0, fcor(ji,jj) ) * sangvg   ! only the sinus changes its sign with the hemisphere
+            zsang(ji,jj)  = SIGN( 1._wp, fcor(ji,jj) ) * sangvg   ! only the sinus changes its sign with the hemisphere
             !
             zmasst(ji,jj) = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) )
             zpresh(ji,jj) = tms(ji,jj) *  pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) )
 !!gm  :: stress given at I-point (F-point for the ocean) only compute the ponderation with the ice fraction (1-frld)
-            zi1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )
-            zi2(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )
+            zi1(ji,jj)    = tms(ji,jj) * ( 1._wp - frld(ji,jj) )
+            zi2(ji,jj)    = tms(ji,jj) * ( 1._wp - frld(ji,jj) )
          END DO
       END DO
@@ -161,5 +153,5 @@
             zstms = zztms(ji,jj  ) * wght(ji,jj,2,2) + zztms(ji-1,jj  ) * wght(ji,jj,1,2)   &
                &  + zztms(ji,jj-1) * wght(ji,jj,2,1) + zztms(ji-1,jj-1) * wght(ji,jj,1,1)
-            zusw  = 1.0 / MAX( zstms, epsd )
+            zusw  = 1._wp / MAX( zstms, epsd )
 
             zt11 = zztms(ji  ,jj  ) * zzfrld(ji  ,jj  ) 
@@ -199,7 +191,7 @@
             ! Gradient of the sea surface height
             zgsshx =  (   (ssh_m(ji  ,jj  ) - ssh_m(ji-1,jj  ))/e1u(ji-1,jj  )   &
-               &       +  (ssh_m(ji  ,jj-1) - ssh_m(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5
+               &       +  (ssh_m(ji  ,jj-1) - ssh_m(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5_wp
             zgsshy =  (   (ssh_m(ji  ,jj  ) - ssh_m(ji  ,jj-1))/e2v(ji  ,jj-1)   &
-               &       +  (ssh_m(ji-1,jj  ) - ssh_m(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5
+               &       +  (ssh_m(ji-1,jj  ) - ssh_m(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5_wp
 
             ! Computation of the velocity field taking into account the ice-ice interaction.                                 
@@ -217,5 +209,5 @@
          !                                                ! ==================== !        
          zindu = MOD( iter , 2 )
-         zusdtp = ( zindu * 2.0 + ( 1.0 - zindu ) * 1.0 )  * REAL( nbiter ) / rdt_ice
+         zusdtp = ( zindu * 2._wp + ( 1._wp - zindu ) * 1._wp )  * REAL( nbiter ) / rdt_ice
 
          ! Computation of free drift field for free slip boundary conditions.
@@ -239,8 +231,8 @@
                zdgi = zt12 + zt21
                ztrace2 = zdgp * zdgp 
-               zdeter  = zt11 * zt22 - 0.25 * zdgi * zdgi
+               zdeter  = zt11 * zt22 - 0.25_wp * zdgi * zdgi
 
                !  Creep limit depends on the size of the grid.
-               zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2 ),  creepl)
+               zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4._wp * zdeter ) * usecc2 ),  creepl)
 
                !-  Computation of viscosities.
@@ -254,8 +246,8 @@
             DO ji = 2, fs_jpim1   ! NO vector opt.
                !* zc1u , zc2v
-               zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm
-               zvis12 =       zviseta (ji-1,jj-1) + dm
-               zvis21 =       zviseta (ji-1,jj-1)
-               zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
+               zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
+               zvis12 =         zviseta (ji-1,jj-1) + dm
+               zvis21 =         zviseta (ji-1,jj-1)
+               zvis22 =         zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
                zdiag  = zvis22 * ( akappa(ji-1,jj-1,1,1) + akappa(ji-1,jj-1,2,1) )
                zs11_11 =  zvis11 * akappa(ji-1,jj-1,1,1) + zdiag
@@ -264,8 +256,8 @@
                zs22_11 =  zvis11 * akappa(ji-1,jj-1,2,1) + zdiag
 
-               zvis11 = 2.0 * zviseta (ji,jj-1) + dm
-               zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
-               zvis12 =       zviseta (ji,jj-1) + dm
-               zvis21 =       zviseta (ji,jj-1)
+               zvis11 = 2._wp * zviseta (ji,jj-1) + dm
+               zvis22 =         zviszeta(ji,jj-1) - zviseta(ji,jj-1)
+               zvis12 =         zviseta (ji,jj-1) + dm
+               zvis21 =         zviseta (ji,jj-1)
                zdiag = zvis22 * ( -akappa(ji,jj-1,1,1) + akappa(ji,jj-1,2,1) )
                zs11_21 = -zvis11 * akappa(ji,jj-1,1,1) + zdiag
@@ -274,9 +266,9 @@
                zs21_21 = -zvis12 * akappa(ji,jj-1,1,2) + zvis21 * akappa(ji,jj-1,2,2)
 
-               zvis11 = 2.0 * zviseta (ji-1,jj) + dm
-               zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
-               zvis12 =       zviseta (ji-1,jj) + dm
-               zvis21 =       zviseta (ji-1,jj)
-               zdiag = zvis22 * ( akappa(ji-1,jj,1,1) + akappa(ji-1,jj,2,1) )
+               zvis11 = 2._wp * zviseta (ji-1,jj) + dm
+               zvis22 =         zviszeta(ji-1,jj) - zviseta(ji-1,jj)
+               zvis12 =         zviseta (ji-1,jj) + dm
+               zvis21 =         zviseta (ji-1,jj)
+               zdiag  = zvis22 * ( akappa(ji-1,jj,1,1) + akappa(ji-1,jj,2,1) )
                zs11_12 =  zvis11 * akappa(ji-1,jj,1,1) + zdiag
                zs12_12 = -zvis12 * akappa(ji-1,jj,2,2) - zvis21 * akappa(ji-1,jj,1,2)
@@ -284,8 +276,8 @@
                zs21_12 = -zvis12 * akappa(ji-1,jj,1,2) - zvis21 * akappa(ji-1,jj,2,2)
 
-               zvis11 = 2.0 * zviseta (ji,jj) + dm
-               zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj)
-               zvis12 =       zviseta (ji,jj) + dm
-               zvis21 =       zviseta (ji,jj)
+               zvis11 = 2._wp * zviseta (ji,jj) + dm
+               zvis22 =         zviszeta(ji,jj) - zviseta(ji,jj)
+               zvis12 =         zviseta (ji,jj) + dm
+               zvis21 =         zviseta (ji,jj)
                zdiag = zvis22 * ( -akappa(ji,jj,1,1) + akappa(ji,jj,2,1) )
                zs11_22 = -zvis11 * akappa(ji,jj,1,1) + zdiag
@@ -313,8 +305,8 @@
 
                !* zc1v , zc2v.
-               zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm
-               zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
-               zvis12 =       zviseta (ji-1,jj-1) + dm
-               zvis21 =       zviseta (ji-1,jj-1)
+               zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
+               zvis22 =         zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
+               zvis12 =         zviseta (ji-1,jj-1) + dm
+               zvis21 =         zviseta (ji-1,jj-1)
                zdiag = zvis22 * ( akappa(ji-1,jj-1,1,2) + akappa(ji-1,jj-1,2,2) )
                zs11_11 =  zvis11 * akappa(ji-1,jj-1,1,2) + zdiag
@@ -323,8 +315,8 @@
                zs21_11 =  zvis12 * akappa(ji-1,jj-1,1,1) - zvis21 * akappa(ji-1,jj-1,2,1)
  
-               zvis11 = 2.0 * zviseta (ji,jj-1) + dm
-               zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
-               zvis12 =       zviseta (ji,jj-1) + dm
-               zvis21 =       zviseta (ji,jj-1)
+               zvis11 = 2._wp * zviseta (ji,jj-1) + dm
+               zvis22 =         zviszeta(ji,jj-1) - zviseta(ji,jj-1)
+               zvis12 =         zviseta (ji,jj-1) + dm
+               zvis21 =         zviseta (ji,jj-1)
                zdiag = zvis22 * ( akappa(ji,jj-1,1,2) + akappa(ji,jj-1,2,2) )
                zs11_21 =  zvis11 * akappa(ji,jj-1,1,2) + zdiag
@@ -333,8 +325,8 @@
                zs21_21 = -zvis12 * akappa(ji,jj-1,1,1) - zvis21 * akappa(ji,jj-1,2,1)
 
-               zvis11 = 2.0 * zviseta (ji-1,jj) + dm
-               zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
-               zvis12 =       zviseta (ji-1,jj) + dm
-               zvis21 =       zviseta (ji-1,jj)
+               zvis11 = 2._wp * zviseta (ji-1,jj) + dm
+               zvis22 =         zviszeta(ji-1,jj) - zviseta(ji-1,jj)
+               zvis12 =         zviseta (ji-1,jj) + dm
+               zvis21 =         zviseta (ji-1,jj)
                zdiag = zvis22 * ( akappa(ji-1,jj,1,2) - akappa(ji-1,jj,2,2) )
                zs11_12 =  zvis11 * akappa(ji-1,jj,1,2) + zdiag
@@ -343,8 +335,8 @@
                zs21_12 =  zvis12 * akappa(ji-1,jj,1,1) - zvis21 * akappa(ji-1,jj,2,1)
 
-               zvis11 = 2.0 * zviseta (ji,jj) + dm
-               zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj)
-               zvis12 =       zviseta (ji,jj) + dm
-               zvis21 =       zviseta (ji,jj)
+               zvis11 = 2._wp * zviseta (ji,jj) + dm
+               zvis22 =         zviszeta(ji,jj) - zviseta(ji,jj)
+               zvis12 =         zviseta (ji,jj) + dm
+               zvis21 =         zviseta (ji,jj)
                zdiag = zvis22 * ( akappa(ji,jj,1,2) - akappa(ji,jj,2,2) )
                zs11_22 =  zvis11 * akappa(ji,jj,1,2) + zdiag
@@ -386,8 +378,8 @@
                   ze22 = + akappa(ji,jj-1,2,2) * zv_a(ji+1,jj) + akappa(ji,jj-1,2,1) * zu_a(ji+1,jj)
                   ze21 =   akappa(ji,jj-1,1,1) * zv_a(ji+1,jj) - akappa(ji,jj-1,1,2) * zu_a(ji+1,jj)
-                  zvis11 = 2.0 * zviseta (ji,jj-1) + dm
-                  zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
-                  zvis12 =       zviseta (ji,jj-1) + dm
-                  zvis21 =       zviseta (ji,jj-1)
+                  zvis11 = 2._wp * zviseta (ji,jj-1) + dm
+                  zvis22 =         zviszeta(ji,jj-1) - zviseta(ji,jj-1)
+                  zvis12 =         zviseta (ji,jj-1) + dm
+                  zvis21 =         zviseta (ji,jj-1)
                   zdiag = zvis22 * ( ze11 + ze22 )
                   zs11_21 =  zvis11 * ze11 + zdiag
@@ -404,8 +396,8 @@
                   ze21 =   akappa(ji-1,jj,1,1) * ( zv_a(ji  ,jj+1) - zv_a(ji-1,jj+1) )   &
                      &   - akappa(ji-1,jj,1,2) * ( zu_a(ji  ,jj+1) + zu_a(ji-1,jj+1) )
-                  zvis11 = 2.0 * zviseta (ji-1,jj) + dm
-                  zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
-                  zvis12 =       zviseta (ji-1,jj) + dm
-                  zvis21 =       zviseta (ji-1,jj)
+                  zvis11 = 2._wp * zviseta (ji-1,jj) + dm
+                  zvis22 =         zviszeta(ji-1,jj) - zviseta(ji-1,jj)
+                  zvis12 =         zviseta (ji-1,jj) + dm
+                  zvis21 =         zviseta (ji-1,jj)
                   zdiag = zvis22 * ( ze11 + ze22 )
                   zs11_12 =  zvis11 * ze11 + zdiag
@@ -422,8 +414,8 @@
                   ze21 =   akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji  ,jj+1) )   &
                      &   - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji  ,jj+1) )
-                  zvis11 = 2.0 * zviseta (ji,jj) + dm
-                  zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj)
-                  zvis12 =       zviseta (ji,jj) + dm
-                  zvis21 =       zviseta (ji,jj)
+                  zvis11 = 2._wp * zviseta (ji,jj) + dm
+                  zvis22 =         zviszeta(ji,jj) - zviseta(ji,jj)
+                  zvis12 =         zviseta (ji,jj) + dm
+                  zvis21 =         zviseta (ji,jj)
                   zdiag = zvis22 * ( ze11 + ze22 )
                   zs11_22 =  zvis11 * ze11 + zdiag
@@ -441,8 +433,8 @@
                   ze21 =   akappa(ji-1,jj-1,1,1) * ( zv_a(ji  ,jj-1) - zv_a(ji-1,jj-1) - zv_a(ji-1,jj) )   &
                      &  -  akappa(ji-1,jj-1,1,2) * ( zu_a(ji  ,jj-1) + zu_a(ji-1,jj-1) + zu_a(ji-1,jj) )
-                  zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm
-                  zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
-                  zvis12 =       zviseta (ji-1,jj-1) + dm
-                  zvis21 =       zviseta (ji-1,jj-1)
+                  zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
+                  zvis22 =         zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
+                  zvis12 =         zviseta (ji-1,jj-1) + dm
+                  zvis21 =         zviseta (ji-1,jj-1)
                   zdiag = zvis22 * ( ze11 + ze22 )
                   zs11_11 =  zvis11 * ze11 + zdiag
@@ -459,8 +451,8 @@
                   ze21 =   akappa(ji,jj-1,1,1) * ( zv_a(ji+1,jj-1) - zv_a(ji  ,jj-1) )   &
                      &   - akappa(ji,jj-1,1,2) * ( zu_a(ji+1,jj-1) + zu_a(ji  ,jj-1) )
-                  zvis11 = 2.0 * zviseta (ji,jj-1) + dm
-                  zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
-                  zvis12 =       zviseta (ji,jj-1) + dm
-                  zvis21 =       zviseta (ji,jj-1)
+                  zvis11 = 2._wp * zviseta (ji,jj-1) + dm
+                  zvis22 =         zviszeta(ji,jj-1) - zviseta(ji,jj-1)
+                  zvis12 =         zviseta (ji,jj-1) + dm
+                  zvis21 =         zviseta (ji,jj-1)
                   zdiag = zvis22 * ( ze11 + ze22 )
                   zs11_21 =  zs11_21 + zvis11 * ze11 + zdiag
@@ -473,8 +465,8 @@
                   ze22 = - akappa(ji-1,jj,2,2) * zv_a(ji-1,jj) + akappa(ji-1,jj,2,1) * zu_a(ji-1,jj)
                   ze21 = - akappa(ji-1,jj,1,1) * zv_a(ji-1,jj) - akappa(ji-1,jj,1,2) * zu_a(ji-1,jj)
-                  zvis11 = 2.0 * zviseta (ji-1,jj) + dm
-                  zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
-                  zvis12 =       zviseta (ji-1,jj) + dm
-                  zvis21 =       zviseta (ji-1,jj)
+                  zvis11 = 2._wp * zviseta (ji-1,jj) + dm
+                  zvis22 =         zviszeta(ji-1,jj) - zviseta(ji-1,jj)
+                  zvis12 =         zviseta (ji-1,jj) + dm
+                  zvis21 =         zviseta (ji-1,jj)
                   zdiag = zvis22 * ( ze11 + ze22 )
                   zs11_12 =  zs11_12 + zvis11 * ze11 + zdiag
@@ -504,12 +496,12 @@
                   zvr     = zv_a(ji,jj) - v_oce(ji,jj)
 !!!!
-                  zmod    = SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) )
+                  zmod    = SQRT( zur*zur + zvr*zvr ) * ( 1._wp - zfrld(ji,jj) )
                   za      = rhoco * zmod
 !!!!
-!!gm chg resul    za      = rhoco * SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) )
+!!gm chg resul    za      = rhoco * SQRT( zur*zur + zvr*zvr ) * ( 1._wp - zfrld(ji,jj) )
                   zac     = za * cangvg
                   zmpzas  = alpha * zcorl(ji,jj) + za * zsang(ji,jj)
                   zmassdt = zusdtp * zmass(ji,jj)
-                  zcorlal = ( 1.0 - alpha ) * zcorl(ji,jj)
+                  zcorlal = ( 1._wp - alpha ) * zcorl(ji,jj)
 
                   za1 =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   &
@@ -525,5 +517,5 @@
                   zunw   = (  ( za1 + zd1 ) * zc2 + ( za2 + zd2 ) * zc1 ) * zden
                   zvnw   = (  ( za2 + zd2 ) * zb1 - ( za1 + zd1 ) * zb2 ) * zden
-                  zmask  = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj)
+                  zmask  = ( 1._wp - MAX( rzero, SIGN( rone , 1._wp - zmass(ji,jj) ) ) ) * tmu(ji,jj)
 
                   zu_n(ji,jj) = ( zu_a(ji,jj) + om * ( zunw - zu_a(ji,jj) ) * tmu(ji,jj) ) * zmask
@@ -580,5 +572,5 @@
 #else
    !!----------------------------------------------------------------------
-   !!   Default option          Dummy module       NO 2.0 LIM sea-ice model
+   !!   Default option        Dummy module      NO VP & LIM-2 sea-ice model
    !!----------------------------------------------------------------------
 CONTAINS
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limrst_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limrst_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limrst_2.F90	(revision 2528)
@@ -6,9 +6,9 @@
    !! History :  2.0  !  01-04  (C. Ethe, G. Madec)  Original code
    !!                 !  06-07  (S. Masson)  use IOM for restart read/write
+   !!            3.3  !  09-05  (G.Garric) addition of the lim2_evp case
    !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   lim_rst_opn_2   : open ice restart file
@@ -16,11 +16,10 @@
    !!   lim_rst_read_2  : read  the ice restart file 
    !!----------------------------------------------------------------------
-   USE dom_oce         ! ocean space and time domain
-   USE ice_2
-   USE sbc_oce
-   USE sbc_ice
-
-   USE in_out_manager
-   USE iom
+   USE dom_oce          ! ocean domain
+   USE ice_2            ! LIM-2: sea-ice variables
+   USE sbc_oce          ! Surface Boundary Condition: ocean
+   USE sbc_ice          ! Surface Boundary Condition: sea-ice
+   USE in_out_manager   ! I/O manager
+   USE iom              ! I/O library
 
    IMPLICIT NONE
@@ -35,9 +34,8 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -84,4 +82,5 @@
    END SUBROUTINE lim_rst_opn_2
 
+
    SUBROUTINE lim_rst_write_2( kt )
       !!----------------------------------------------------------------------
@@ -108,50 +107,55 @@
       CALL iom_rstput( iter, nitrst, numriw, 'kt_ice', REAL( iter, wp) ) 
       
-      CALL iom_rstput( iter, nitrst, numriw, 'hicif' , hicif (:,:)   )      ! prognostic variables 
-      CALL iom_rstput( iter, nitrst, numriw, 'hsnif' , hsnif (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'frld'  , frld  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sist'  , sist  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'tbif1' , tbif  (:,:,1) )
-      CALL iom_rstput( iter, nitrst, numriw, 'tbif2' , tbif  (:,:,2) )
-      CALL iom_rstput( iter, nitrst, numriw, 'tbif3' , tbif  (:,:,3) )
-      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, 'qstoif', qstoif(:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq' , fsbbq (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxice' , sxice (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'syice' , syice (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxxice', sxxice(:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'syyice', syyice(:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxyice', sxyice(:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxsn'  , sxsn  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sysn'  , sysn  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxxsn' , sxxsn (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'syysn' , syysn (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxysn' , sxysn (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxa'   , sxa   (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sya'   , sya   (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxxa'  , sxxa  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'syya'  , syya  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxya'  , sxya  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxc0'  , sxc0  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'syc0'  , syc0  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxxc0' , sxxc0 (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'syyc0' , syyc0 (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxyc0' , sxyc0 (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxc1'  , sxc1  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'syc1'  , syc1  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxxc1' , sxxc1 (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'syyc1' , syyc1 (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxyc1' , sxyc1 (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxc2'  , sxc2  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'syc2'  , syc2  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxxc2' , sxxc2 (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'syyc2' , syyc2 (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxyc2' , sxyc2 (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxst'  , sxst  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'syst'  , syst  (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxxst' , sxxst (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'syyst' , syyst (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'sxyst' , sxyst (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'hicif'      , hicif (:,:)   )      ! prognostic variables 
+      CALL iom_rstput( iter, nitrst, numriw, 'hsnif'      , hsnif (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'frld'       , frld  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sist'       , sist  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'tbif1'      , tbif  (:,:,1) )
+      CALL iom_rstput( iter, nitrst, numriw, 'tbif2'      , tbif  (:,:,2) )
+      CALL iom_rstput( iter, nitrst, numriw, 'tbif3'      , tbif  (:,:,3) )
+      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, 'qstoif'     , qstoif(:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq'      , fsbbq (:,:)   )
+#if ! defined key_lim2_vp
+      CALL iom_rstput( iter, nitrst, numriw, 'stress1_i'  , stress1_i (:,:) )    ! EVP rheology
+      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i'  , stress2_i (:,:) )
+      CALL iom_rstput( iter, nitrst, numriw, 'stress12_i' , stress12_i(:,:) )
+#endif
+      CALL iom_rstput( iter, nitrst, numriw, 'sxice'      , sxice (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'syice'      , syice (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxxice'     , sxxice(:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'syyice'     , syyice(:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxyice'     , sxyice(:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxsn'       , sxsn  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sysn'       , sysn  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxxsn'      , sxxsn (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'syysn'      , syysn (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxysn'      , sxysn (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxa'        , sxa   (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sya'        , sya   (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxxa'       , sxxa  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'syya'       , syya  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxya'       , sxya  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxc0'       , sxc0  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'syc0'       , syc0  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxxc0'      , sxxc0 (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'syyc0'      , syyc0 (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxyc0'      , sxyc0 (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxc1'       , sxc1  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'syc1'       , syc1  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxxc1'      , sxxc1 (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'syyc1'      , syyc1 (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxyc1'      , sxyc1 (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxc2'       , sxc2  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'syc2'       , syc2  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxxc2'      , sxxc2 (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'syyc2'      , syyc2 (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxyc2'      , sxyc2 (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxst'       , sxst  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'syst'       , syst  (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxxst'      , sxxst (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'syyst'      , syyst (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'sxyst'      , sxyst (:,:)   )
       
       IF( iter == nitrst ) THEN
@@ -218,41 +222,46 @@
       ENDIF
 
-      CALL iom_get( numrir, jpdom_autoglo, 'qstoif', qstoif )    
-      CALL iom_get( numrir, jpdom_autoglo, 'fsbbq' , fsbbq  )    
-      CALL iom_get( numrir, jpdom_autoglo, 'sxice' , sxice  )
-      CALL iom_get( numrir, jpdom_autoglo, 'syice' , syice  )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxxice', sxxice )
-      CALL iom_get( numrir, jpdom_autoglo, 'syyice', syyice )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxyice', sxyice )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxsn'  , sxsn   )
-      CALL iom_get( numrir, jpdom_autoglo, 'sysn'  , sysn   )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxxsn' , sxxsn  )
-      CALL iom_get( numrir, jpdom_autoglo, 'syysn' , syysn  )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxysn' , sxysn  )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxa'   , sxa    )
-      CALL iom_get( numrir, jpdom_autoglo, 'sya'   , sya    )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxxa'  , sxxa   )
-      CALL iom_get( numrir, jpdom_autoglo, 'syya'  , syya   )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxya'  , sxya   )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxc0'  , sxc0   )
-      CALL iom_get( numrir, jpdom_autoglo, 'syc0'  , syc0   )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxxc0' , sxxc0  )
-      CALL iom_get( numrir, jpdom_autoglo, 'syyc0' , syyc0  )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxyc0' , sxyc0  )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxc1'  , sxc1   )
-      CALL iom_get( numrir, jpdom_autoglo, 'syc1'  , syc1   )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxxc1' , sxxc1  )
-      CALL iom_get( numrir, jpdom_autoglo, 'syyc1' , syyc1  )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxyc1' , sxyc1  )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxc2'  , sxc2   )
-      CALL iom_get( numrir, jpdom_autoglo, 'syc2'  , syc2   )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxxc2' , sxxc2  )
-      CALL iom_get( numrir, jpdom_autoglo, 'syyc2' , syyc2  )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxyc2' , sxyc2  )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxst'  , sxst   )
-      CALL iom_get( numrir, jpdom_autoglo, 'syst'  , syst   )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxxst' , sxxst  )
-      CALL iom_get( numrir, jpdom_autoglo, 'syyst' , syyst  )
-      CALL iom_get( numrir, jpdom_autoglo, 'sxyst' , sxyst  )
+      CALL iom_get( numrir, jpdom_autoglo, 'qstoif'     , qstoif )    
+      CALL iom_get( numrir, jpdom_autoglo, 'fsbbq'      , fsbbq  )    
+#if ! defined key_lim2_vp
+      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 )
+#endif
+      CALL iom_get( numrir, jpdom_autoglo, 'sxice'      , sxice  )
+      CALL iom_get( numrir, jpdom_autoglo, 'syice'      , syice  )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxxice'     , sxxice )
+      CALL iom_get( numrir, jpdom_autoglo, 'syyice'     , syyice )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxyice'     , sxyice )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxsn'       , sxsn   )
+      CALL iom_get( numrir, jpdom_autoglo, 'sysn'       , sysn   )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxxsn'      , sxxsn  )
+      CALL iom_get( numrir, jpdom_autoglo, 'syysn'      , syysn  )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxysn'      , sxysn  )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxa'        , sxa    )
+      CALL iom_get( numrir, jpdom_autoglo, 'sya'        , sya    )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxxa'       , sxxa   )
+      CALL iom_get( numrir, jpdom_autoglo, 'syya'       , syya   )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxya'       , sxya   )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxc0'       , sxc0   )
+      CALL iom_get( numrir, jpdom_autoglo, 'syc0'       , syc0   )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxxc0'      , sxxc0  )
+      CALL iom_get( numrir, jpdom_autoglo, 'syyc0'      , syyc0  )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxyc0'      , sxyc0  )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxc1'       , sxc1   )
+      CALL iom_get( numrir, jpdom_autoglo, 'syc1'       , syc1   )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxxc1'      , sxxc1  )
+      CALL iom_get( numrir, jpdom_autoglo, 'syyc1'      , syyc1  )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxyc1'      , sxyc1  )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxc2'       , sxc2   )
+      CALL iom_get( numrir, jpdom_autoglo, 'syc2'       , syc2   )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxxc2'      , sxxc2  )
+      CALL iom_get( numrir, jpdom_autoglo, 'syyc2'      , syyc2  )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxyc2'      , sxyc2  )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxst'       , sxst   )
+      CALL iom_get( numrir, jpdom_autoglo, 'syst'       , syst   )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxxst'      , sxxst  )
+      CALL iom_get( numrir, jpdom_autoglo, 'syyst'      , syyst  )
+      CALL iom_get( numrir, jpdom_autoglo, 'sxyst'      , sxyst  )
       
       CALL iom_close( numrir )
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 2528)
@@ -2,9 +2,11 @@
    !!======================================================================
    !!                       ***  MODULE limsbc_2   ***
-   !!           computation of the flux at the sea ice/ocean interface
+   !! LIM-2 :   updates the fluxes at the ocean surface with ice-ocean fluxes
    !!======================================================================
-   !! History : 00-01 (H. Goosse) Original code
-   !!           02-07 (C. Ethe, G. Madec) re-writing F90
-   !!           06-07 (G. Madec) surface module
+   !! History :  LIM  ! 2000-01 (H. Goosse) Original code
+   !!            1.0  ! 2002-07 (C. Ethe, G. Madec) re-writing F90
+   !!            3.0  ! 2006-07 (G. Madec) surface module
+   !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
+   !!             -   ! 2010-11 (G. Madec) ice-ocean stress computed at each ocean time-step
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -12,18 +14,19 @@
    !!   'key_lim2'                                    LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!   lim_sbc_2  : flux at the ice / ocean interface
+   !!   lim_sbc_flx_2  : update mass, heat and salt fluxes at the ocean surface
+   !!   lim_sbc_tau_2  : update i- and j-stresses, and its modulus at the ocean surface
    !!----------------------------------------------------------------------
    USE par_oce          ! ocean parameters
+   USE phycst           ! physical constants
    USE dom_oce          ! ocean domain
-   USE sbc_ice          ! surface boundary condition
-   USE sbc_oce          ! surface boundary condition
-   USE phycst           ! physical constants
-   USE ice_2            ! LIM sea-ice variables
-
-   USE lbclnk           ! ocean lateral boundary condition
+   USE dom_ice_2        ! LIM-2: ice domain
+   USE ice_2            ! LIM-2: ice variables
+   USE sbc_ice          ! surface boundary condition: ice
+   USE sbc_oce          ! surface boundary condition: ocean
+
+   USE lbclnk           ! ocean lateral boundary condition - MPP exchanges
    USE in_out_manager   ! I/O manager
    USE diaar5, ONLY :   lk_diaar5
-   USE iom              ! 
+   USE iom              ! I/O library
    USE albedo           ! albedo parameters
    USE prtctl           ! Print control
@@ -33,23 +36,27 @@
    PRIVATE
 
-   PUBLIC lim_sbc_2     ! called by sbc_ice_lim_2
-
-   REAL(wp)  ::   epsi16 = 1.e-16  ! constant values
-   REAL(wp)  ::   rzero  = 0.e0    
-   REAL(wp)  ::   rone   = 1.e0
-   REAL(wp), DIMENSION(jpi,jpj)  ::   soce_r
-   REAL(wp), DIMENSION(jpi,jpj)  ::   sice_r
+   PUBLIC   lim_sbc_flx_2   ! called by sbc_ice_lim_2
+   PUBLIC   lim_sbc_tau_2   ! called by sbc_ice_lim_2
+
+   REAL(wp)  ::   r1_rdtice            ! = 1. / rdt_ice 
+   REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values
+   REAL(wp)  ::   rzero  = 0._wp       !     -      -
+   REAL(wp)  ::   rone   = 1._wp       !     -      -
+   !
+   REAL(wp), DIMENSION(jpi,jpj) ::   soce_0, sice_0   ! constant SSS and ice salinity used in levitating sea-ice case
+
+   REAL(wp), DIMENSION(jpi,jpj) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress              [N/m2]
+   REAL(wp), DIMENSION(jpi,jpj) ::   tmod_io              ! modulus of the ice-ocean relative velocity   [m/s]
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
-   SUBROUTINE lim_sbc_2( kt )
+   SUBROUTINE lim_sbc_flx_2( kt )
       !!-------------------------------------------------------------------
       !!                ***  ROUTINE lim_sbc_2 ***
@@ -75,50 +82,35 @@
       !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
       !!---------------------------------------------------------------------
-      INTEGER ::   kt    ! number of iteration
+      INTEGER, INTENT(in) ::   kt    ! number of iteration
       !!
-      INTEGER  ::   ji, jj           ! dummy loop indices
-      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
-      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
-      REAL(wp) ::   zrdtir           ! 1. / rdt_ice
-      REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux
-      REAL(wp) ::   zinda            ! switch for testing the values of ice concentration
-      REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface
-      REAL(wp) ::   zemp             ! freshwater exchanges at the ice/ocean interface
-      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points
-      REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points
-      REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity
-! interface 2D --> 3D
-      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb     ! albedo of ice under overcast sky
-      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalbp    ! albedo of ice under clear sky
-      REAL(wp) ::   zsang, zmod, zztmp, zfm
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! component of ocean stress below sea-ice at I-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztiomi           ! module    of ocean stress below sea-ice at I-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   zqnsoce          ! save qns before its modification by ice model
-
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      INTEGER  ::   ii0, ii1, ij0, ij1         ! local integers
+      INTEGER  ::   ifvt, i1mfr, idfr, iflt    !   -       -
+      INTEGER  ::   ial, iadv, ifral, ifrdv    !   -       -
+      REAL(wp) ::   zqsr, zqns, zfm            ! local scalars
+      REAL(wp) ::   zinda, zfons, zemp         !   -      -
+      REAL(wp), DIMENSION(jpi,jpj)   ::   zqnsoce       ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb, zalbp   ! 2D/3D workspace
       !!---------------------------------------------------------------------
      
-      zrdtir = 1. / rdt_ice
-      
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~   '
-
-         soce_r(:,:) = soce
-         sice_r(:,:) = sice
-         !
-         IF( cp_cfg == "orca" ) THEN
-           !   ocean/ice salinity in the Baltic sea 
-           DO jj = 1, jpj
-              DO ji = 1, jpi
-                 IF( glamt(ji,jj) >= 14. .AND.  glamt(ji,jj) <= 32. .AND. gphit(ji,jj) >= 54. .AND. gphit(ji,jj) <= 66. ) THEN 
-                   soce_r(ji,jj) = 4.e0 
-                   sice_r(ji,jj) = 2.e0
-                 END IF
-              END DO
-           END DO
-           !
-         END IF
-      END IF
+         IF(lwp) WRITE(numout,*) 'lim_sbc_flx_2 : LIM-2 sea-ice - surface boundary condition - Mass, heat & salt fluxes'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   '
+         !
+         r1_rdtice = 1. / rdt_ice
+         !
+         soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating sea-ice case
+         sice_0(:,:) = sice
+         !
+         IF( cp_cfg == "orca" ) THEN            ! decrease ocean & ice reference salinities in the Baltic sea 
+            WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &
+               &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
+               soce_0(:,:) = 4._wp
+               sice_0(:,:) = 2._wp
+            END WHERE
+         ENDIF
+         !
+      ENDIF
 
       !------------------------------------------!
@@ -126,8 +118,4 @@
       !------------------------------------------!
 
-!!gm
-!!gm CAUTION   
-!!gm re-verifies the non solar expression, especially over open ocen
-!!gm
       zqnsoce(:,:) = qns(:,:)
       DO jj = 1, jpj
@@ -186,9 +174,9 @@
             zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads
                &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            &
-               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * zrdtir    &
-               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * zrdtir
+               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice    &
+               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice 
 
             fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ???
-            
+            !
             qsr  (ji,jj) = zqsr                                          ! solar heat flux 
             qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux
@@ -198,104 +186,41 @@
       CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )      
       CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )      
-      CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1. - pfrld(:,:)) )
+      CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1.e0 - pfrld(:,:)) )
 
       !------------------------------------------!
       !      mass flux at the ocean surface      !
       !------------------------------------------!
-
-!!gm
-!!gm CAUTION   
-!!gm re-verifies the emp & emps expression, especially the absence of 1-frld on zfm
-!!gm
       DO jj = 1, jpj
          DO ji = 1, jpi
-            
+            !
 #if defined key_coupled
-          zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   ! 
-             &   + rdmsnif(ji,jj) * zrdtir                                      !  freshwaterflux due to snow melting 
+            ! freshwater exchanges at the ice-atmosphere / ocean interface (coupled mode)
+            zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   ! 
+               &   + rdmsnif(ji,jj) * r1_rdtice                                   !  freshwaterflux due to snow melting 
 #else
-!!$            !  computing freshwater exchanges at the ice/ocean interface
-!!$            zpme = - evap(ji,jj)    *   frld(ji,jj)           &   !  evaporation over oceanic fraction
-!!$               &   + tprecip(ji,jj)                           &   !  total precipitation
-!!$               &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )   &   !  remov. snow precip over ice
-!!$               &   - rdmsnif(ji,jj) / rdt_ice                     !  freshwaterflux due to snow melting 
             !  computing freshwater exchanges at the ice/ocean interface
             zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction 
                &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean
                &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step
-               &   + rdmsnif(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting 
+               &   + rdmsnif(ji,jj) * r1_rdtice                    !  freshwaterflux due to snow melting 
                !                                                   !  ice-covered fraction:
 #endif            
-
+            !
             !  computing salt exchanges at the ice/ocean interface
-            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * zrdtir ) 
-            
+            zfons =  ( soce_0(ji,jj) - sice_0(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice ) 
+            !
             !  converting the salt flux from ice to a freshwater flux from ocean
             zfm  = zfons / ( sss_m(ji,jj) + epsi16 )
-            
+            !
             emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution)
             emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution)
-
+            !
          END DO
       END DO
 
-      IF( lk_diaar5 ) THEN
-         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * zrdtir )
-         CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdmicif(:,:) * zrdtir )
-         CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdmicif(:,:) * zrdtir )
-      ENDIF
-
-      !------------------------------------------!
-      !    momentum flux at the ocean surface    !
-      !------------------------------------------!
-
-      IF ( ln_limdyn ) THEN                        ! Update the stress over ice-over area (only in ice-dynamic case)
-         !                                         ! otherwise the atmosphere-ocean stress is used everywhere
-
-         ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point)
-!CDIR NOVERRCHK
-         DO jj = 1, jpj
-!CDIR NOVERRCHK
-            DO ji = 1, jpi
-               ! ... change the cosinus angle sign in the south hemisphere
-               zsang  = SIGN(1.e0, gphif(ji,jj) ) * sangvg
-               ! ... ice velocity relative to the ocean at I-point
-               zu_io  = u_ice(ji,jj) - u_oce(ji,jj)
-               zv_io  = v_ice(ji,jj) - v_oce(ji,jj)
-               zmod   = SQRT( zu_io * zu_io + zv_io * zv_io )
-               zztmp  = rhoco * zmod
-               ! ... components of ice stress over ocean with a ice-ocean rotation angle (at I-point)
-               ztio_u(ji,jj) = zztmp * ( cangvg * zu_io - zsang * zv_io )
-               ztio_v(ji,jj) = zztmp * ( cangvg * zv_io + zsang * zu_io )
-               ! ... module of ice stress over ocean (at I-point)
-               ztiomi(ji,jj) = zztmp * zmod
-               ! 
-            END DO
-         END DO
-
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1   ! NO vector opt.
-               ! ... components of ice-ocean stress at U and V-points  (from I-point values)
-               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
-               zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
-               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
-               zfrldu = 0.5 * ( frld(ji,jj) + frld(ji+1,jj  ) )
-               zfrldv = 0.5 * ( frld(ji,jj) + frld(ji  ,jj+1) )
-               ! ... update components of surface ocean stress (ice-cover wheighted)
-               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * zutau
-               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * zvtau
-               ! ... module of ice-ocean stress at T-points (from I-point values)
-               zztmp = 0.25 * ( ztiomi(ji,jj) + ztiomi(ji+1,jj) + ztiomi(ji,jj+1) + ztiomi(ji+1,jj+1) )
-               ! ... update module of surface ocean stress (ice-cover wheighted)
-               taum(ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1. - frld(ji,jj) ) * zztmp
-               !
-            END DO
-         END DO
-
-         ! boundary condition on the stress (utau,vtau,taum)
-         CALL lbc_lnk( utau, 'U', -1. )
-         CALL lbc_lnk( vtau, 'V', -1. )
-         CALL lbc_lnk( taum, 'T',  1. )
-
+      IF( lk_diaar5 ) THEN       ! AR5 diagnostics
+         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * r1_rdtice )
+         CALL iom_put( 'fsal_virt_cea',   soce_0(:,:) * rdmicif(:,:) * r1_rdtice )
+         CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdmicif(:,:) * r1_rdtice )
       ENDIF
 
@@ -304,5 +229,5 @@
       !-----------------------------------------------!
 
-      IF ( lk_cpl ) THEN           
+      IF( lk_cpl ) THEN          ! coupled case
          ! Ice surface temperature 
          tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature       
@@ -313,5 +238,5 @@
       ENDIF
 
-      IF(ln_ctl) THEN
+      IF(ln_ctl) THEN            ! control print
          CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
          CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ')
@@ -320,14 +245,171 @@
          CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ')
       ENDIF 
-   
-    END SUBROUTINE lim_sbc_2
+      !
+   END SUBROUTINE lim_sbc_flx_2
+
+
+   SUBROUTINE lim_sbc_tau_2( kt , pu_oce, pv_oce )
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE lim_sbc_tau ***
+      !!  
+      !! ** Purpose : Update the ocean surface stresses due to the ice
+      !!         
+      !! ** Action  : * at each ice time step (every nn_fsbc time step):
+      !!                - compute the modulus of ice-ocean relative velocity 
+      !!                  at T-point (C-grid) or I-point (B-grid)
+      !!                      tmod_io = rhoco * | U_ice-U_oce |
+      !!                - update the modulus of stress at ocean surface
+      !!                      taum = frld * taum + (1-frld) * tmod_io * | U_ice-U_oce |
+      !!              * at each ocean time step (each kt): 
+      !!                  compute linearized ice-ocean stresses as
+      !!                      Utau = tmod_io * | U_ice - pU_oce |
+      !!                using instantaneous current ocean velocity (usually before)
+      !!
+      !!    NB: - the averaging operator used depends on the ice dynamics grid (cp_ice_msh='I' or 'C')
+      !!        - ice-ocean rotation angle only allowed in cp_ice_msh='I' case
+      !!        - here we make an approximation: taum is only computed every ice time step
+      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids 
+      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
+      !!
+      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
+      !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
+      !!---------------------------------------------------------------------
+      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
+      !!
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      REAL(wp) ::   zfrldu, zat_u, zu_i, zutau_ice, zu_t, zmodt   ! local scalar
+      REAL(wp) ::   zfrldv, zat_v, zv_i, zvtau_ice, zv_t, zmodi   !   -      -
+      REAL(wp) ::   zsang, zumt                                          !    -         -
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice
+      !!---------------------------------------------------------------------
+      !
+      IF( kt == nit000 .AND. lwp ) THEN         ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'lim_sbc_tau_2 : LIM 2.0 sea-ice - surface ocean momentum fluxes'
+         WRITE(numout,*) '~~~~~~~~~~~~~ '
+         IF( lk_lim2_vp )   THEN   ;   WRITE(numout,*) '                VP  rheology - B-grid case'
+         ELSE                      ;   WRITE(numout,*) '                EVP rheology - C-grid case'
+         ENDIF
+      ENDIF
+      !
+      SELECT CASE( cp_ice_msh )     
+      !                             !-----------------------!
+      CASE( 'I' )                   !  B-grid ice dynamics  !   I-point (i.e. F-point with sea-ice indexation)
+         !                          !--=--------------------!
+         !
+         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==! (i.e. surface module time-step)
+!CDIR NOVERRCHK
+            DO jj = 1, jpj                               !* modulus of ice-ocean relative velocity at I-point
+!CDIR NOVERRCHK
+               DO ji = 1, jpi
+                  zu_i  = u_ice(ji,jj) - u_oce(ji,jj)                   ! ice-ocean relative velocity at I-point
+                  zv_i  = v_ice(ji,jj) - v_oce(ji,jj)
+                  tmod_io(ji,jj) = SQRT( zu_i * zu_i + zv_i * zv_i )    ! modulus of this velocity (at I-point)
+               END DO
+            END DO
+!CDIR NOVERRCHK
+            DO jj = 1, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
+!CDIR NOVERRCHK
+               DO ji = 1, jpim1   ! NO vector opt.
+                  !                                               ! modulus of U_ice-U_oce at T-point
+                  zumt  = 0.25_wp * (  tmod_io(ji+1,jj) + tmod_io(ji+1,jj+1)    &   
+                     &               + tmod_io(ji,jj+1) + tmod_io(ji+1,jj+1)  )
+                  !                                               ! update the modulus of stress at ocean surface
+                  taum(ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1._wp - frld(ji,jj) ) * rhoco * zumt * zumt
+               END DO
+            END DO
+            CALL lbc_lnk( taum, 'T', 1. )
+            !
+            utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
+            vtau_oce(:,:) = vtau(:,:)
+            !
+         ENDIF
+         !
+         !                                        !==  at each ocean time-step  ==!
+         !
+         !                                               !* ice/ocean stress WITH a ice-ocean rotation angle at I-point
+         DO jj = 2, jpj
+            zsang  = SIGN( 1._wp, gphif(1,jj) ) * sangvg          ! change the cosine angle sign in the SH 
+            DO ji = 2, jpi    ! NO vect. opt. possible
+               ! ... ice-ocean relative velocity at I-point using instantaneous surface ocean current at u- & v-pts
+               zu_i = u_ice(ji,jj) - 0.5_wp * ( pu_oce(ji-1,jj  ) + pu_oce(ji-1,jj-1) ) * tmu(ji,jj)
+               zv_i = v_ice(ji,jj) - 0.5_wp * ( pv_oce(ji  ,jj-1) + pv_oce(ji-1,jj-1) ) * tmu(ji,jj)
+               ! ... components of stress with a ice-ocean rotation angle 
+               zmodi = rhoco * tmod_io(ji,jj)                     
+               ztio_u(ji,jj) = zmodi * ( cangvg * zu_i - zsang * zv_i )
+               ztio_v(ji,jj) = zmodi * ( cangvg * zv_i + zsang * zu_i )
+            END DO
+         END DO
+         !                                               !* surface ocean stresses at u- and v-points
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1   ! NO vector opt.
+               !                                   ! ice-ocean stress at U and V-points  (from I-point values)
+               zutau_ice  = 0.5_wp * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
+               zvtau_ice  = 0.5_wp * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
+               !                                   ! open-ocean (lead) fraction at U- & V-points (from T-point values)
+               zfrldu = 0.5_wp * ( frld(ji,jj) + frld(ji+1,jj  ) )
+               zfrldv = 0.5_wp * ( frld(ji,jj) + frld(ji  ,jj+1) )
+               !                                   ! update the surface ocean stress (ice-cover wheighted)
+               utau(ji,jj) = zfrldu * utau_oce(ji,jj) + ( 1._wp - zfrldu ) * zutau_ice
+               vtau(ji,jj) = zfrldv * vtau_oce(ji,jj) + ( 1._wp - zfrldv ) * zvtau_ice
+            END DO
+         END DO
+         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )     ! lateral boundary condition
+         !
+         !
+         !                          !-----------------------!
+      CASE( 'C' )                   !  C-grid ice dynamics  !   U & V-points (same as in the ocean)
+         !                          !--=--------------------!
+         !
+         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==! (i.e. surface module time-step)
+!CDIR NOVERRCHK
+            DO jj = 2, jpjm1                          !* modulus of the ice-ocean velocity at T-point
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1
+                  zu_t  = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   ! 2*(U_ice-U_oce) at T-point
+                  zv_t  = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)      
+                  zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )                      ! |U_ice-U_oce|^2
+                  !                                               ! update the modulus of stress at ocean surface
+                  taum   (ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1._wp - frld(ji,jj) ) * rhoco * zmodt
+                  tmod_io(ji,jj) = SQRT( zmodt ) * rhoco          ! rhoco*|Uice-Uoce|
+               END DO
+            END DO
+            CALL lbc_lnk( taum, 'T', 1. )   ;   CALL lbc_lnk( tmod_io, 'T', 1. )
+            !
+            utau_oce(:,:) = utau(:,:)                 !* save the air-ocean stresses at ice time-step
+            vtau_oce(:,:) = vtau(:,:)
+            !
+         ENDIF
+         !
+         !                                        !==  at each ocean time-step  ==!
+         !
+         DO jj = 2, jpjm1                             !* ice stress over ocean WITHOUT a ice-ocean rotation angle
+            DO ji = fs_2, fs_jpim1
+               !                                            ! ocean area at u- & v-points
+               zfrldu  = 0.5_wp * ( frld(ji,jj) + frld(ji+1,jj  ) )
+               zfrldv  = 0.5_wp * ( frld(ji,jj) + frld(ji  ,jj+1) )
+               !                                            ! quadratic drag formulation without rotation
+               !                                            ! using instantaneous surface ocean current
+               zutau_ice = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
+               zvtau_ice = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
+               !                                            ! update the surface ocean stress (ice-cover wheighted)
+               utau(ji,jj) = zfrldu * utau_oce(ji,jj) + ( 1._wp - zfrldu ) * zutau_ice
+               vtau(ji,jj) = zfrldv * vtau_oce(ji,jj) + ( 1._wp - zfrldv ) * zvtau_ice
+            END DO
+         END DO
+         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
+         !
+      END SELECT
+
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
+         &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask )
+      !  
+   END SUBROUTINE lim_sbc_tau_2
 
 #else
    !!----------------------------------------------------------------------
-   !!   Default option :        Dummy module       NO LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_sbc_2         ! Dummy routine
-   END SUBROUTINE lim_sbc_2
+   !!   Default option         Empty module        NO LIM 2.0 sea-ice model
+   !!----------------------------------------------------------------------
 #endif 
 
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limtab_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limtab_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limtab_2.F90	(revision 2528)
@@ -20,7 +20,7 @@
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limthd_2.F90	(revision 2528)
@@ -50,7 +50,7 @@
 #  include "vectopt_loop_substitute.h90"
    !!-------- -------------------------------------------------------------
-   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2009) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limthd_lac_2.F90	(revision 2528)
@@ -28,7 +28,7 @@
       zone   = 1.e0
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 2528)
@@ -33,7 +33,7 @@
       &          zone   = 1.e0
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 2528)
@@ -7,4 +7,5 @@
    !!            2.0  !  2001-05 (G. Madec, R. Hordoir) opa norm
    !!             -   !  2004-01 (G. Madec, C. Ethe)  F90, mpp
+   !!            3.3  !  2009-05  (G. Garric, C. Bricaud) addition of the lim2_evp case
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -32,19 +33,18 @@
    PUBLIC   lim_trp_2   ! called by sbc_ice_lim_2
 
-   REAL(wp), PUBLIC  ::   bound  = 0.e0   !: boundary condit. (0.0 no-slip, 1.0 free-slip)
-
-   REAL(wp)  ::           &  ! constant values
-      epsi06 = 1.e-06  ,  &
-      epsi03 = 1.e-03  ,  &
-      epsi16 = 1.e-16  ,  &
-      rzero  = 0.e0    ,  &
-      rone   = 1.e0
+   REAL(wp), PUBLIC ::   bound  = 0.e0          !: boundary condit. (0.0 no-slip, 1.0 free-slip)
+
+   REAL(wp)  ::   epsi06 = 1.e-06   ! constant values
+   REAL(wp)  ::   epsi03 = 1.e-03  
+   REAL(wp)  ::   epsi16 = 1.e-16  
+   REAL(wp)  ::   rzero  = 0.e0   
+   REAL(wp)  ::   rone   = 1.e0
 
    !! * Substitution
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/LIM 3.2,  UCL-LOCEAN-IPSL (2010) 
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -87,17 +87,22 @@
          ! ice velocities at ocean U- and V-points (zui_u,zvi_v)
          ! ---------------------------------------
-         zvbord = 1.0 + ( 1.0 - bound )      ! zvbord=2 no-slip, =0 free slip boundary conditions        
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1   ! NO vector opt.
-               zui_u(ji,jj) = ( u_ice(ji+1,jj  ) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj  ) + tmu(ji+1,jj+1), zvbord ) )
-               zvi_v(ji,jj) = ( v_ice(ji  ,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji  ,jj+1) + tmu(ji+1,jj+1), zvbord ) )
-            END DO
-         END DO
-         CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )         ! Lateral boundary conditions
-
+         IF( lk_lim2_vp ) THEN      ! VP rheology : B-grid sea-ice dynamics (I-point ice velocity)
+            zvbord = 1._wp + ( 1._wp - bound )      ! zvbord=2 no-slip, =0 free slip boundary conditions        
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1   ! NO vector opt.
+                  zui_u(ji,jj) = ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj)+tmu(ji+1,jj+1), zvbord ) )
+                  zvi_v(ji,jj) = ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji,jj+1)+tmu(ji+1,jj+1), zvbord ) )
+               END DO
+            END DO
+            CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )      ! Lateral boundary conditions
+            !
+         ELSE                       ! EVP rheology : C-grid sea-ice dynamics (u- & v-points ice velocity)
+            zui_u(:,:) = u_ice(:,:)      ! EVP rheology: ice (u,v) at u- and v-points
+            zvi_v(:,:) = v_ice(:,:)
+         ENDIF
 
          ! CFL test for stability
          ! ----------------------
-         zcfl  = 0.e0
+         zcfl  = 0._wp
          zcfl  = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, :     ) ) * rdt_ice / e1u(1:jpim1, :     ) ) )
          zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) )
@@ -109,6 +114,6 @@
          ! content of properties
          ! ---------------------
-         zs0sn (:,:) =  hsnm(:,:) * area(:,:)                 ! Snow volume.
-         zs0ice(:,:) =  hicm(:,:) * area(:,:)                 ! Ice volume.
+         zs0sn (:,:) =  hsnm(:,:)              * area  (:,:)  ! Snow volume.
+         zs0ice(:,:) =  hicm(:,:)              * area  (:,:)  ! Ice volume.
          zs0a  (:,:) =  ( 1.0 - frld(:,:) )    * area  (:,:)  ! Surface covered by ice.
          zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn (:,:)  ! Heat content of the snow layer.
@@ -188,6 +193,6 @@
 !        DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
 !           DO ji = 1 , fs_jpim1   ! vector opt.
-!              IF( MIN( zs0a(ji,jj) , zs0a(ji+1,jj) ) == 0.e0 )   pahu(ji,jj) = 0.e0
-!              IF( MIN( zs0a(ji,jj) , zs0a(ji,jj+1) ) == 0.e0 )   pahv(ji,jj) = 0.e0
+!              IF( MIN( zs0a(ji,jj) , zs0a(ji+1,jj) ) == 0.e0 )   pahu(ji,jj) = 0._wp
+!              IF( MIN( zs0a(ji,jj) , zs0a(ji,jj+1) ) == 0.e0 )   pahv(ji,jj) = 0._wp
 !           END DO
 !        END DO
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limwri_2.F90	(revision 2528)
@@ -60,7 +60,7 @@
 #   include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006)
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -112,5 +112,5 @@
          CALL dia_nam ( clhstnam, nwrite, 'icemod' )
          CALL histbeg ( clhstnam, jpi, glamt, jpj, gphit,    &
-            &           1, jpi, 1, jpj, niter, zjulian, rdt_ice, nhorid, nice , domain_id=nidom)
+            &           1, jpi, 1, jpj, niter, zjulian, rdt_ice, nhorid, nice , domain_id=nidom, snc4chunks=snc4set)
          CALL histvert( nice, "deptht", "Vertical T levels", "m", 1, zdept, ndepid, "down")
          CALL wheneq  ( jpij , tmask(:,:,1), 1, 1., ndex51, ndim)
@@ -120,5 +120,5 @@
                &                                  , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout )
          END DO
-         CALL histend( nice )
+         CALL histend( nice, snc4set )
          !
       ENDIF
@@ -305,5 +305,5 @@
       CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s", jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
 
-      CALL histend( kid )   ! end of the file definition
+      CALL histend( kid, snc4set )   ! end of the file definition
 
       CALL histwrite( kid, "isnowthi", kt, hsnif          , jpi*jpj, (/1/) )   
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/limwri_dimg_2.h90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/limwri_dimg_2.h90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/limwri_dimg_2.h90	(revision 2528)
@@ -1,7 +1,7 @@
     SUBROUTINE lim_wri_2(kt)
    !!----------------------------------------------------------------------
-   !!  LIM 2.0, UCL-LOCEAN-IPSL (2005)
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
     !!-------------------------------------------------------------------
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/par_ice_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/par_ice_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/par_ice_2.F90	(revision 2528)
@@ -24,7 +24,7 @@
 #endif
    !!----------------------------------------------------------------------
-   !! NEMO/LIM 2.0, UCL-LOCEAN-IPSL (2008)
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE par_ice_2
Index: trunk/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 2528)
@@ -8,7 +8,7 @@
    !!   2.0  !  02-11  (C. Ethe)  F90: Free form and module
    !!----------------------------------------------------------------------
-   !!   LIM 2.0, UCL-LOCEAN-IPSL (2005)
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
    !! * Modules used
