Index: branches/devmercator2010_1/NEMO/LIM_SRC_2/dom_ice_2.F90
===================================================================
--- branches/devmercator2010_1/NEMO/LIM_SRC_2/dom_ice_2.F90	(revision 2134)
+++ branches/devmercator2010_1/NEMO/LIM_SRC_2/dom_ice_2.F90	(revision 2135)
@@ -4,11 +4,10 @@
    !! LIM 2.0 Sea Ice :   Domain  variables
    !!======================================================================
-   !! History :   2.0  !  03-08  (C. Ethe)  Free form and module
+   !! History :   1.0  !  2003-08  (C. Ethe)  Free form and module
+   !!             3.3  !  2009-05 (G.Garric, C. Bricaud) addition of lim2_evp case
    !!----------------------------------------------------------------------
-#if defined key_lim2
+#if defined   key_lim2
    !!----------------------------------------------------------------------
-   !!   LIM 2.0, UCL-LOCEAN-IPSL (2005)
-   !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!   'key_lim2'                                       LIM2 sea-ice model
    !!----------------------------------------------------------------------
    USE par_ice_2
@@ -20,15 +19,30 @@
 
    INTEGER, PUBLIC ::   njeq , njeqm1          !: j-index of the equator if it is inside the domain
-      !                                        !  (otherwise = jpj+10 (SH) or -10 (SH) )
+   !                                           !  (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
+
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!======================================================================
-#endif
 END MODULE dom_ice_2
Index: branches/devmercator2010_1/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- branches/devmercator2010_1/NEMO/LIM_SRC_2/ice_2.F90	(revision 2134)
+++ branches/devmercator2010_1/NEMO/LIM_SRC_2/ice_2.F90	(revision 2135)
@@ -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,9 +11,4 @@
    !!   '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
 
@@ -20,14 +16,15 @@
    PRIVATE
 
-   !!* Share parameters namelist (namicerun read in iceini) *
+   !                                                                     !!* namicerun namelist *
    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)
+   LOGICAL               , PUBLIC ::   ln_nicep      = .TRUE.             !: flag grid points output (T) or not (F)
+   REAL(wp)              , PUBLIC ::   hsndif        = 0.e0               !: computation of snow temp. (0) or not (9999)
+   REAL(wp)              , PUBLIC ::   hicdif        = 0.e0               !: computation of ice  temp. (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) *
+   !                                       !!* 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
@@ -46,5 +43,7 @@
    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)
-
+   INTEGER , PUBLIC ::   nevp   =   360     !: number of EVP subcycling iterations
+   INTEGER , PUBLIC ::   telast =  3600     !: timescale for EVP elastic waves
+   REAL(wp), PUBLIC ::   alphaevp = 1.e0    !: coefficient for the solution of EVP int. stresses
    REAL(wp), PUBLIC ::   usecc2             !:  = 1.0 / ( ecc * ecc )
    REAL(wp), PUBLIC ::   rhoco              !: = rau0 * cw
@@ -54,9 +53,31 @@
    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) ::   hsnm , hicm   !: mean snow and ice thicknesses
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s                 !: friction velocity
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s         !: friction velocity
 
-   !!* diagnostic quantities
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature (Kelvin)
+# if defined key_lim2_vp
+   !                                                        !!* VP rheology *
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm     !: mean snow and ice thicknesses
+   CHARACTER(len=1), PUBLIC :: cl_grid =   'B'                   !: type of grid used in ice dynamics, 'C' or 'B'
+# else
+   !                                                        !!* EVP rheology *
+   CHARACTER(len=1), PUBLIC             ::   cl_grid = 'C'   !: type of grid used in ice dynamics, 'C' or 'B'
+   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          !:
+   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) ::   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,21 +92,21 @@
    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 bring  ocean surface layer to freezing 
    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       !: Array linked with the max brine pockets heat content
    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         !: Array used to store energy (total lateral ablation case)
    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
 
-   !!* moment used in the advection scheme
+   !                                                                               !!* moment used in the advection
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxice, syice, sxxice, syyice, sxyice   !: for ice  volume
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxsn,  sysn,  sxxsn,  syysn,  sxysn    !: for snow volume
@@ -98,8 +119,12 @@
 #else
    !!----------------------------------------------------------------------
-   !!   Default option         Empty module        NO LIM 2.0 sea-ice model
+   !!   Default option          Empty module         NO LIM-2 sea-ice model
    !!----------------------------------------------------------------------
 #endif
 
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE ice_2
Index: branches/devmercator2010_1/NEMO/LIM_SRC_2/iceini_2.F90
===================================================================
--- branches/devmercator2010_1/NEMO/LIM_SRC_2/iceini_2.F90	(revision 2134)
+++ branches/devmercator2010_1/NEMO/LIM_SRC_2/iceini_2.F90	(revision 2135)
@@ -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
@@ -11,30 +12,30 @@
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
    !!   ice_init_2       : sea-ice model initialization
    !!   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
    PRIVATE
 
-   PUBLIC   ice_init_2               ! called by sbcice_lim_2.F90
+   PUBLIC   ice_init_2   ! called by sbcice_lim_2.F90
+
+   INTEGER, PUBLIC ::   numit   !: iteration number
 
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Id$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -43,24 +44,22 @@
       !!                  ***  ROUTINE ice_init_2  ***
       !!
-      !! ** purpose :   
+      !! ** purpose :   initialisation of LIM-2 domain and variables 
       !!----------------------------------------------------------------------
       !
-      ! Open the namelist 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
+      CALL ice_run_2                  !  read in namelist some run parameters
                  
-      ! Louvain la Neuve Ice model
-      rdt_ice = nn_fsbc * rdttra(1)
+      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
+      !                               ! ice initialisation (start from rest or from a restart)
+      IF( .NOT.ln_rstart ) THEN   ;   CALL lim_istate_2     
+      ELSE                        ;   CALL lim_rst_read_2
       ENDIF
       
-      tn_ice(:,:,1) = sist(:,:)         ! initialisation of ice temperature   
+      tn_ice(:,:,1) = sist(:,:)       ! initialisation of ice temperature   
       fr_i  (:,:) = 1.0 - frld(:,:)   ! initialisation of sea-ice fraction    
       !
@@ -75,5 +74,5 @@
       !!
       !! ** Method  :   Read the namicerun namelist and check the parameter 
-      !!       values called at the first timestep (nit000)
+      !!              values called at the first timestep (nit000)
       !!
       !! ** input   :   Namelist namicerun
@@ -82,8 +81,8 @@
       !!-------------------------------------------------------------------
       !                    
-      REWIND ( numnam_ice )                       ! Read Namelist namicerun 
-      READ   ( numnam_ice , namicerun )
-
-      IF(lwp) THEN
+      REWIND( numnam_ice )                   ! Read namicerun namelist
+      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: branches/devmercator2010_1/NEMO/LIM_SRC_2/limdyn_2.F90
===================================================================
--- branches/devmercator2010_1/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 2134)
+++ branches/devmercator2010_1/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 2135)
@@ -8,4 +8,5 @@
    !!             2.0  !  03-08  (C. Ethe) add lim_dyn_init
    !!             2.0  !  06-07  (G. Madec)  Surface module
+   !!             3.3  !  09-05  (G.Garric, C. Bricaud) addition of the lim2_evp case
    !!---------------------------------------------------------------------
 #if defined key_lim2
@@ -16,32 +17,33 @@
    !!    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 domain
+   USE sbc_oce          ! surface boundary condition variables
+   USE phycst           ! physical constant
+   USE ice_2            ! LIM2: ice variables
+   USE dom_ice_2        ! LIM2: ice domain
+   USE limistate_2      ! LIM2: ice initial state
+#if defined key_lim2_vp
+   USE limrhg_2         ! LIM2:  VP ice rheology
+#else
+   USE limrhg           ! LIM : EVP ice rheology
+#endif
+   USE lbclnk           ! lateral boundary condition - MPP exchanges
+   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 module
+
+   REAL(wp) ::  rone = 1.e0   ! constant value
 
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010)
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -71,57 +73,62 @@
       IF( kt == nit000 )   CALL lim_dyn_init_2   ! Initialization (first time-step only)
       
-      IF( ln_limdyn ) THEN
-         !
-         ! Mean ice and snow thicknesses.          
-         hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:)
+      IF( ln_limdyn ) THEN                     ! Rheology (ice dynamics)
+         !                                     ! ========
+         !         
+         hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:)   ! Mean ice and snow thicknesses
          hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:)
          !
-         !                                     ! Rheology (ice dynamics)
-         !                                     ! ========
+         !                                     ! Define the j-limits where ice rheology is computed
          
-         !  Define the j-limits where ice rheology is computed
-         ! ---------------------------------------------------
-         
-         IF( lk_mpp .OR. nbit_cmp == 1 ) THEN                    ! mpp: compute over the whole domain
+         IF( lk_mpp .OR. nbit_cmp == 1 ) 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 defined key_lim2_vp
+            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  ==!
+            DO jj = 1, jpj
+               zind(jj) = SUM( frld (:,jj  ) )        ! = FLOAT(jpj) if ocean everywhere on a j-line
+               zmsk(jj) = SUM( tmask(:,jj,1) )        ! = 0          if land  everywhere on a j-line
+            END DO
             !
-         ELSE                                 ! optimization of the computational area
-            !
-            DO jj = 1, jpj
-               zind(jj) = SUM( frld (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line
-               zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line
-            END DO
-            !
-            IF( l_jeq ) THEN                     ! local domain include both hemisphere
-               !                                 ! Rheology is computed in each hemisphere
-               !                                 ! only over the ice cover latitude strip
-               ! Northern hemisphere
-               i_j1  = njeq
+            IF( l_jeq ) THEN                          ! local domain include both hemisphere: rheology is computed
+               !                                      ! in each hemisphere only over the ice cover latitude strip
+               i_j1  = njeq                  ! Northern hemisphere
                i_jpj = jpj
                DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
                   i_j1 = i_j1 + 1
                END DO
+#if defined key_lim2_vp
                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 )
-               ! 
-               ! Southern hemisphere
-               i_j1  =  1 
+#else
+               i_j1 = MAX( 1, i_j1-2 )
+               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
+               CALL lim_rhg( i_j1, i_jpj )
+#endif
+               i_j1  =  1                    ! Southern hemisphere
                i_jpj = njeq
                DO WHILE ( i_jpj >= 1 .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
                   i_jpj = i_jpj - 1
                END DO
+#if defined key_lim2_vp
                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 )
+#else
+               i_jpj = MIN( jpj, i_jpj+1 )
+               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
+               CALL lim_rhg( i_j1, i_jpj ) !!!!cbr CALL lim_rhg( i_j1, i_jpj, kt )
+#endif
                ! 
-            ELSE                                 ! local domain extends over one hemisphere only
-               !                                 ! Rheology is computed only over the ice cover
-               !                                 ! latitude strip
+            ELSE                                      ! local domain extends over one hemisphere only: rheology is
+               !                                      ! computed only over the ice cover latitude strip
                i_j1  = 1
                DO WHILE ( i_j1 <= jpj .AND. zind(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 )
@@ -129,33 +136,42 @@
                END DO
                i_j1 = MAX( 1, i_j1-1 )
-    
                i_jpj  = jpj
                DO WHILE ( i_jpj >= 1  .AND. zind(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 )
                   i_jpj = i_jpj - 1
                END DO
+               IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
+#if defined key_lim2_vp
                i_jpj = MIN( jpj, i_jpj+2)
-    
-               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 )
+               CALL lim_rhg_2( i_j1, i_jpj )                !  VP rheology
+#else
+               i_j1 = MAX( 1, i_j1-2 )
+               i_jpj = MIN( jpj, i_jpj+1)
+               CALL lim_rhg  ( i_j1, i_jpj )                ! EVP rheology
+#endif
                !
             ENDIF
             !
          ENDIF
-
+         !
          IF(ln_ctl)   CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
+         ! 
          
-         ! 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
-         ! frictional velocity at T-point
-         DO jj = 2, jpjm1
+         !                                     ! friction velocity
+         !                                     ! =================
+         SELECT CASE( cl_grid )
+         CASE( 'C' )                          ! C-grid ice dynamics (EVP)
+            zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)      ! ice-ocean & ice velocity at ocean velocity points 
+            zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
+            !
+         CASE( 'B' )                          ! B-grid ice dynamics (VP)
+            DO jj = 1, jpjm1                          ! ice velocity at I-point, ice-ocean velocity at ocean points
+               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
+         END SELECT
+         !
+         DO jj = 2, jpjm1                     ! frictional velocity at T-point
             DO ji = 2, jpim1   ! NO vector opt. because of zu_io
                ust2s(ji,jj) = 0.5 * cw                                                          &
@@ -165,6 +181,6 @@
          END DO
          !
-      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
-         !
+      ELSE                                     ! no ice dynamics : transmit directly the atmospheric stress to the ocean
+         !                                     ! ===============
          zcoef = SQRT( 0.5 ) / rau0
          DO jj = 2, jpjm1
@@ -177,8 +193,8 @@
       ENDIF
       !
-      CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
+      CALL lbc_lnk( ust2s, 'T',  1. )   ! lateral boundary condition
       !
       IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
-
+      !
    END SUBROUTINE lim_dyn_2
 
@@ -188,20 +204,17 @@
       !!                  ***  ROUTINE lim_dyn_init_2  ***
       !!
-      !! ** Purpose :   Physical constants and parameters linked to the ice
-      !!              dynamics
-      !!
-      !! ** Method  :   Read the namicedyn namelist and check the ice-dynamic
-      !!              parameter values
+      !! ** Purpose :   initialisation of  the ice dynamics variables
+      !!              
+      !! ** Method  :   Read the namicedyn namelist and check their values
       !!
       !! ** input   :   Namelist namicedyn
       !!-------------------------------------------------------------------
-      NAMELIST/namicedyn/ epsd, alpha,     &
-         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
-         &                c_rhg, etamn, creepl, ecc, ahi0
-      !!-------------------------------------------------------------------
-
+      NAMELIST/namicedyn/ epsd, alpha, dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
+         &                c_rhg, etamn, creepl, ecc, ahi0, nevp, telast, alphaevp
+      !!-------------------------------------------------------------------
+      !
       REWIND ( numnam_ice )                       ! Read Namelist namicedyn
       READ   ( numnam_ice  , namicedyn )
-
+      !
       IF(lwp) THEN                                ! Control print
          WRITE(numout,*)
@@ -223,6 +236,9 @@
          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
-
+      !
       !  Initialization
       usecc2 = 1.0 / ( ecc * ecc )
@@ -240,8 +256,8 @@
 #else
    !!----------------------------------------------------------------------
-   !!   Default option          Empty module       NO LIM 2.0 sea-ice model
+   !!   Default option          Dummy module       NO LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE lim_dyn_2         ! Empty routine
+   SUBROUTINE lim_dyn_2         ! Dummy routine
    END SUBROUTINE lim_dyn_2
 #endif 
Index: branches/devmercator2010_1/NEMO/LIM_SRC_2/limmsh_2.F90
===================================================================
--- branches/devmercator2010_1/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 2134)
+++ branches/devmercator2010_1/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 2135)
@@ -4,5 +4,9 @@
    !! LIM 2.0 ice model :   definition of the ice mesh parameters
    !!======================================================================
-#if defined key_lim2
+   !! History :  LIM  ! 2001-04  (Louvain-la-Neuve) Original code
+   !!            1.0  ! 2002-08  (C. Ethe, G. Madec)
+   !!            3.3  ! 2009-05  (G. Garric, C. Bricaud) addition of the lim2_evp case
+   !!----------------------------------------------------------------------
+#if defined   key_lim2
    !!----------------------------------------------------------------------
    !!   'key_lim2'                                     LIM 2.0sea-ice model
@@ -10,23 +14,20 @@
    !!   lim_msh_2   : definition of the ice mesh
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE phycst
-   USE dom_oce
-   USE dom_ice_2
-   USE lbclnk
-   USE in_out_manager
+   USE phycst           ! physical constants
+   USE dom_oce          ! ocean domain
+   USE dom_ice_2        ! LIM2: ice domain
+   USE lbclnk           ! lateral boundary condition - MPP exchanges
+   USE in_out_manager   ! I/O manager
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC lim_msh_2      ! routine called by ice_ini_2.F90
-
-   !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   PUBLIC   lim_msh_2   ! routine called by ice_ini_2.F90
+
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010)
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -42,22 +43,16 @@
       !!              - Initialization of the ice masks (tmsk, umsk)
       !! 
-      !! ** Refer.  : Deleersnijder et al. Ocean Modelling 100, 7-10 
-      !!
-      !! ** History :
-      !!         original    : 01-04 (LIM)
-      !!         addition    : 02-08 (C. Ethe, G. Madec)
+      !! References : Deleersnijder et al. Ocean Modelling 100, 7-10 
       !!--------------------------------------------------------------------- 
-      !! * 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
       !!---------------------------------------------------------------------
-
+      !
       IF(lwp) THEN
          WRITE(numout,*)
@@ -73,11 +68,6 @@
       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
-!i       zmsk(jj) = SUM( tmask(:,jj,:) )   ! = 0          if land  everywhere on a j-line
-!!ii     write(numout,*) jj, zind(jj)
-!i    END DO
-
       IF( fcor(1,1) * fcor(1,nlcj) < 0.e0 ) THEN   ! local domain include both hemisphere
          l_jeq = .TRUE.
@@ -112,12 +102,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
       
-      
-      ! metric coefficients for sea ice dynamic
+#if defined key_lim2_vp
+      ! metric coefficients for sea ice dynamic (VP rheology)
       !----------------------------------------
       !                                                       ! akappa
@@ -125,11 +120,9 @@
          zd1d2(:,jj) = e1v(:,jj) - e1v(:,jj-1)
       END DO
-      CALL lbc_lnk( zd1d2, 'T', -1. )
-
       DO ji = 2, jpi
          zd2d1(ji,:) = e2u(ji,:) - e2u(ji-1,:)
       END DO
-      CALL lbc_lnk( zd2d1, 'T', -1. )
-
+      CALL lbc_lnk( zd1d2, 'T', -1. )   ;   CALL lbc_lnk( zd2d1, 'T', -1. )      ! lateral boundary condition
+      !
       akappa(:,:,1,1) =        1.0 / ( 2.0 * e1t(:,:) )
       akappa(:,:,1,2) = zd1d2(:,:) / ( 4.0 * e1t(:,:) * e2t(:,:) )
@@ -137,6 +130,6 @@
       akappa(:,:,2,2) =        1.0 / ( 2.0 * e2t(:,:) )
       
-      !                                                      ! weights (wght)
-      DO jj = 2, jpj
+      !             
+      DO jj = 2, jpj                                         ! weights (wght) at I-points
          DO ji = 2, jpi
             zusden = 1. / (  ( e1t(ji,jj) + e1t(ji-1,jj  ) )   &
@@ -152,49 +145,51 @@
       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
       !-------------------------------------------------
-
-      DO jj = 2, jpj
+#if defined key_lim2_vp
+      DO jj = 2, jpj                                     ! VP rheology
          DO ji = 2, jpi   ! NO vector opt.
-            zh1p  =  e1t(ji  ,jj  ) * wght(ji,jj,2,2)   &
-               &   + e1t(ji-1,jj  ) * wght(ji,jj,1,2)   &
-               &   + e1t(ji  ,jj-1) * wght(ji,jj,2,1)   &
-               &   + e1t(ji-1,jj-1) * wght(ji,jj,1,1)
-
-            zh2p  =  e2t(ji  ,jj  ) * wght(ji,jj,2,2)   &
-               &   + e2t(ji-1,jj  ) * wght(ji,jj,1,2)   &
-               &   + e2t(ji  ,jj-1) * wght(ji,jj,2,1)   &
-               &   + e2t(ji-1,jj-1) * wght(ji,jj,1,1)
-
-! better written but change the last digit and thus solver in less than 100 timestep
-!           zh1p  = e1t(ji-1,jj  ) * wght(ji,jj,1,2) + e1t(ji,jj  ) * wght(ji,jj,2,2)   &
-!              &  + e1t(ji-1,jj-1) * wght(ji,jj,1,1) + e1t(ji,jj-1) * wght(ji,jj,2,1) 
-
-!           zh2p  = e2t(ji-1,jj  ) * wght(ji,jj,1,2) + e2t(ji,jj  ) * wght(ji,jj,2,2)   &
-!              &  + e2t(ji-1,jj-1) * wght(ji,jj,1,1) + e2t(ji,jj-1) * wght(ji,jj,2,1)
-
-!!ibug =0   zusden = 1.0 / ( zh1p * zh2p * 4.e0 )
-            zusden = 1.0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 )
+            zh1p  =  e1t(ji  ,jj  ) * wght(ji,jj,2,2) + e1t(ji-1,jj  ) * wght(ji,jj,1,2)   &
+               &   + e1t(ji  ,jj-1) * wght(ji,jj,2,1) + e1t(ji-1,jj-1) * wght(ji,jj,1,1)
+            zh2p  =  e2t(ji  ,jj  ) * wght(ji,jj,2,2) + e2t(ji-1,jj  ) * wght(ji,jj,1,2)   &
+               &   + e2t(ji  ,jj-1) * wght(ji,jj,2,1) + e2t(ji-1,jj-1) * wght(ji,jj,1,1)
+            !
+            zusden  = 1.e0 / MAX( zh1p * zh2p * 4.e0 , 1.e-20 )
             zusden2 = zusden * 2.0 
-
-            zd1d2p = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj  ) - e1t(ji,jj-1) + e1t(ji  ,jj)   )
-            zd2d1p = zusden * 0.5 * (  e2t(ji  ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj  ) - e2t(ji-1,jj)   )
-
+            zd1d2p  = zusden * 0.5 * ( -e1t(ji-1,jj-1) + e1t(ji-1,jj  ) - e1t(ji,jj-1) + e1t(ji  ,jj)   )
+            zd2d1p  = zusden * 0.5 * (  e2t(ji  ,jj-1) - e2t(ji-1,jj-1) + e2t(ji,jj  ) - e2t(ji-1,jj)   )
+            !
             alambd(ji,jj,2,2,2,1) = zusden2 * e2t(ji  ,jj-1)
             alambd(ji,jj,2,2,2,2) = zusden2 * e2t(ji  ,jj  )
             alambd(ji,jj,2,2,1,1) = zusden2 * e2t(ji-1,jj-1)
             alambd(ji,jj,2,2,1,2) = zusden2 * e2t(ji-1,jj  )
-
+            !
             alambd(ji,jj,1,1,2,1) = zusden2 * e1t(ji  ,jj-1)
             alambd(ji,jj,1,1,2,2) = zusden2 * e1t(ji  ,jj  )
             alambd(ji,jj,1,1,1,1) = zusden2 * e1t(ji-1,jj-1)
             alambd(ji,jj,1,1,1,2) = zusden2 * e1t(ji-1,jj  )
-
+            !
             alambd(ji,jj,1,2,2,1) = zd1d2p
             alambd(ji,jj,1,2,2,2) = zd1d2p
             alambd(ji,jj,1,2,1,1) = zd1d2p
             alambd(ji,jj,1,2,1,2) = zd1d2p
-
+            !
             alambd(ji,jj,2,1,2,1) = zd2d1p
             alambd(ji,jj,2,1,2,2) = zd2d1p
@@ -203,46 +198,46 @@
          END DO
       END DO
-
-      CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. )      ! CAUTION: even with the lbc_lnk at ice U-V point
-      CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. )      ! the value of wght at jpj is wrong
-      CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. )      ! but it is never used
-      CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. )      ! 
-
-      CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. )      ! CAUTION: idem
-      CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. )      ! 
-      CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. )      !
-      CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. )      !
-
-      CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. )      ! CAUTION: idem
-      CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. )      !
-      CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. )      !
-      CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. )      !
-
-      CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. )      ! CAUTION: idem
-      CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. )      !
-      CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )      !
-      CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )      !
+      !
+      ! lateral boundary conditions
+      ! CAUTION: even with the lbc_lnk at ice U-V point, the value of wght at jpj is wrong but it is never used
+      CALL lbc_lnk( alambd(:,:,2,2,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,2,2,1), 'I', 1. ) 
+      CALL lbc_lnk( alambd(:,:,2,2,1,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,2,2,2), 'I', 1. ) 
+      !
+      CALL lbc_lnk( alambd(:,:,1,1,2,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,1,2,1), 'I', 1. )
+      CALL lbc_lnk( alambd(:,:,1,1,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,1,1,2), 'I', 1. )
+      !
+      CALL lbc_lnk( alambd(:,:,1,2,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,2,2,1), 'I', 1. )
+      CALL lbc_lnk( alambd(:,:,1,2,1,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,1,2,2,2), 'I', 1. )
+      !
+      CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,1,2,2), 'I', 1. )
+      CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )   ;   CALL lbc_lnk( alambd(:,:,2,1,2,1), 'I', 1. )
+#endif
             
 
       ! Initialization of ice masks
       !----------------------------
-      
+      !      
       tms(:,:) = tmask(:,:,1)      ! ice T-point  : use surface tmask
-
-!i here we can use umask with a i and j shift of -1,-1
-      tmu(:,1) = 0.e0
+      !
+#if defined key_lim2_vp
+      ! VP rheology : ice velocity point is I-point
+      tmu(:,1) = 0.e0              ! 
       tmu(1,:) = 0.e0
-      DO jj = 2, jpj               ! ice U.V-point: computed from ice T-point mask
+      DO jj = 2, jpj
          DO ji = 2, jpim1   ! NO vector opt.
             tmu(ji,jj) =  tms(ji,jj) * tms(ji-1,jj) * tms(ji,jj-1) * tms(ji-1,jj-1)            
          END DO
       END DO
-      
-      !--lateral boundary conditions    
-      CALL lbc_lnk( tmu(:,:), 'I', 1. )
-      
-      ! unmasked and masked area of T-grid cell
-      area(:,:) = e1t(:,:) * e2t(:,:)
-      
+      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
+      !
+      area(:,:) = e1t(:,:) * e2t(:,:)        ! unmasked and masked area of T-grid cell
+      !
    END SUBROUTINE lim_msh_2
 
Index: branches/devmercator2010_1/NEMO/LIM_SRC_2/limrhg_2.F90
===================================================================
--- branches/devmercator2010_1/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 2134)
+++ branches/devmercator2010_1/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 2135)
@@ -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 :  LIM  !  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.0 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 phycst         ! physical constant
-   USE ice_2          ! ice variables
-   USE lbclnk         ! lateral boundary condition
+   USE dom_ice_2      ! LIM2: ice domain
+   USE phycst         ! physical constants
+   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
@@ -32,5 +33,5 @@
    PRIVATE
 
-   PUBLIC   lim_rhg_2 ! routine called by lim_dyn
+   PUBLIC   lim_rhg_2   ! routine called by lim_dyn
 
    REAL(wp) ::   rzero   = 0.e0   ! constant value: zero
@@ -40,5 +41,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010)
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -89,12 +90,5 @@
       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
@@ -106,5 +100,5 @@
       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(:,:)
@@ -126,8 +120,4 @@
       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
 
 
@@ -372,8 +362,7 @@
             END DO
          END DO
-
-         ! GAUSS-SEIDEL method
+         ! 
          !                                                      ! ================ !
-iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    !
+iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    ! GAUSS-SEIDEL method
             !                                                   ! ================ !
 !CDIR NOVERRCHK
Index: branches/devmercator2010_1/NEMO/LIM_SRC_2/limrst_2.F90
===================================================================
--- branches/devmercator2010_1/NEMO/LIM_SRC_2/limrst_2.F90	(revision 2134)
+++ branches/devmercator2010_1/NEMO/LIM_SRC_2/limrst_2.F90	(revision 2135)
@@ -6,4 +6,5 @@
    !! 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
@@ -108,50 +109,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 (:,:)   )
+      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 +224,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: branches/devmercator2010_1/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- branches/devmercator2010_1/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 2134)
+++ branches/devmercator2010_1/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 2135)
@@ -4,12 +4,12 @@
    !!           computation of the flux at the sea ice/ocean interface
    !!======================================================================
-   !! 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
    !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
    !!   'key_lim2'                                    LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   lim_sbc_2  : flux at the ice / ocean interface
@@ -17,13 +17,13 @@
    USE par_oce          ! ocean parameters
    USE dom_oce          ! ocean domain
-   USE sbc_ice          ! surface boundary condition
-   USE sbc_oce          ! surface boundary condition
+   USE sbc_ice          ! surface boundary condition: ice
+   USE sbc_oce          ! surface boundary condition: ocean
    USE phycst           ! physical constants
-   USE ice_2            ! LIM sea-ice variables
-
-   USE lbclnk           ! ocean lateral boundary condition
+   USE ice_2            ! LIM2: ice variables
+
+   USE lbclnk           ! ocean lateral boundary condition - MPP exchanges
    USE in_out_manager   ! I/O manager
    USE diaar5, ONLY :   lk_diaar5
-   USE iom              ! 
+   USE iom              ! IOM library
    USE albedo           ! albedo parameters
    USE prtctl           ! Print control
@@ -33,20 +33,20 @@
    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_2   ! called by sbc_ice_lim_2
+
+   REAL(wp)  ::   r1_rdtice                    ! constant values
+   REAL(wp)  ::   epsi16 = 1.e-16              !     -      -
+   REAL(wp)  ::   rzero  = 0.e0                !     -      -
+   REAL(wp)  ::   rone   = 1.e0                !     -      -
+   !
+   REAL(wp), DIMENSION(jpi,jpj) ::   soce_r, sice_r   ! constant SSS and ice salinity used in levitating sea-ice case
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! NEMO/LIM2 3.3, UCL-LOCEAN-IPSL (2010)
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -78,26 +78,16 @@
       !!
       INTEGER  ::   ji, jj           ! dummy loop indices
-      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
-      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
-      INTEGER  ::   ii0, ii1, ij0, ij1  ! temporary integers
-      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  ::   ii0, ii1, ij0, ij1         ! local integers
+      INTEGER  ::   ifvt, i1mfr, idfr, iflt    !   -       -
+      INTEGER  ::   ial, iadv, ifral, ifrdv    !   -       -
+      REAL(wp) ::   zqsr, zqns, zsang, zmod, zfm   ! local scalars
+      REAL(wp) ::   zinda, zfons, zemp, zztmp      !   -      -
+      REAL(wp) ::   zfrldu, zutau, zu_io           !   -      -
+      REAL(wp) ::   zfrldv, zvtau, zv_io           !   -      -
+      REAL(wp), DIMENSION(jpi,jpj)   ::   ztio_u, ztio_v    ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj)   ::   ztiomi, zqnsoce   !  -     -
+      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb, zalbp   ! 2D/3D workspace
       !!---------------------------------------------------------------------
      
-      zrdtir = 1. / rdt_ice
       
       IF( kt == nit000 ) THEN
@@ -105,12 +95,11 @@
          IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~   '
-
+         !
+         r1_rdtice = 1. / rdt_ice
+         !
          soce_r(:,:) = soce
          sice_r(:,:) = sice
          !
-         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN
-            !                                        ! =======================
-            !                                        !  ORCA_R2 configuration
-            !                                        ! =======================
+         IF( cp_cfg == "orca"  .AND. jp_cfg == 2 ) THEN     !  ORCA_R2 configuration
             ii0 = 145   ;   ii1 = 180        ! Baltic Sea
             ij0 = 113   ;   ij1 = 130   ;   soce_r(mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 4.e0
@@ -175,5 +164,5 @@
 !!$
 
-            !   computation the solar flux at ocean surface
+            ! solar flux at ocean surface
 #if defined key_coupled 
             zqsr = qsr_tot(ji,jj) + ( fstric(ji,jj) - qsr_ice(ji,jj,1) ) * ( 1.0 - pfrld(ji,jj) )
@@ -181,12 +170,12 @@
             zqsr = pfrld(ji,jj) * qsr(ji,jj)  + ( 1.  - pfrld(ji,jj) ) * fstric(ji,jj)
 #endif            
-            !  computation the non solar heat flux at ocean surface
+            ! non solar heat flux at ocean surface
             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
-
+               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                              &
+               &       + 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
@@ -194,52 +183,39 @@
       END DO
 
-      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( 'hflx_ice_cea', - fdtcn(:,:)                         )      
+      CALL iom_put( 'qns_io_cea'  , qns(:,:) - zqnsoce(:,:) * 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) )   &   ! atmosphere-ocean freshwater flux
+               &                  + rdmsnif(ji,jj) * r1_rdtice                   ! freshwater flux 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 
-               !                                                   !  ice-covered fraction:
+            ! freshwater exchanges at the ice-atmosphere / ocean interface (forced mode)
+            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) )   &   ! (account for change in ice cover within the timestep
+               &   + rdmsnif(ji,jj) * r1_rdtice                   ! freshwaterflux due to snow melting 
 #endif            
-
-            !  computing salt exchanges at the ice/ocean interface
-            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * zrdtir ) 
-            
-            !  converting the salt flux from ice to a freshwater flux from ocean
+            ! salt exchanges at the ice/ocean interface
+            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice ) 
+            !
+            ! convert the salt flux from ice into 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 )
+         CALL iom_put( 'isnwmlt_cea'  ,                 rdmsnif(:,:) * r1_rdtice )
+         CALL iom_put( 'fsal_virt_cea',   soce_r(:,:) * rdmicif(:,:) * r1_rdtice )
+         CALL iom_put( 'fsal_real_cea', - sice_r(:,:) * rdmicif(:,:) * r1_rdtice )
       ENDIF
 
@@ -275,6 +251,11 @@
             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) )
+#if defined key_lim2_vp
+               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )      ! VP rheology
                zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
+#else
+               zutau  = ztio_u(ji,jj)                                      ! EVP rheology
+               zvtau  = ztio_v(ji,jj)
+#endif
                ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
                zfrldu = 0.5 * ( frld(ji,jj) + frld(ji+1,jj  ) )
@@ -290,20 +271,15 @@
             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( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )     ! lateral boundary condition 
          CALL lbc_lnk( taum, 'T',  1. )
 
       ENDIF
 
+      IF( lk_cpl ) THEN           
       !-----------------------------------------------!
       !   Coupling variables                          !
       !-----------------------------------------------!
-
-      IF ( lk_cpl ) THEN           
-         ! Ice surface temperature 
-         tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature       
-         ! Computation of snow/ice and ocean albedo
+         tn_ice(:,:,1) = sist(:,:)           ! sea-ice surface temperature       
+         !                                   ! snow/ice and ocean albedo
          CALL albedo_ice( tn_ice, reshape( hicif, (/jpi,jpj,1/) ), reshape( hsnif, (/jpi,jpj,1/) ), zalbp, zalb )
          alb_ice(:,:,1) =  0.5 * ( zalbp(:,:,1) + zalb (:,:,1) )   ! Ice albedo (mean clear and overcast skys)
@@ -318,6 +294,6 @@
          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_2
 
 #else
Index: branches/devmercator2010_1/NEMO/LIM_SRC_2/limtrp_2.F90
===================================================================
--- branches/devmercator2010_1/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 2134)
+++ branches/devmercator2010_1/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 2135)
@@ -4,4 +4,9 @@
    !! LIM 2.0 transport ice model : sea-ice advection/diffusion
    !!======================================================================
+   !! History :  LIM  !  2000-01  (LIM)  Original code
+   !!            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
    !!----------------------------------------------------------------------
@@ -11,42 +16,35 @@
    !!   lim_trp_init_2 : initialization and namelist read
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE phycst
-   USE dom_oce
+   USE phycst          ! physical constants
+   USE dom_oce         ! ocean domain
    USE in_out_manager  ! I/O manager
-   USE dom_ice_2
-   USE ice_2
-   USE limistate_2
-   USE limadv_2
-   USE limhdf_2
-   USE lbclnk
-   USE lib_mpp
+   USE dom_ice_2       ! LIM2 sea-ice domain
+   USE ice_2           ! LIM2 variables
+   USE limistate_2     ! LIM2 initial state
+   USE limadv_2        ! LIM2 advection
+   USE limhdf_2        ! LIM2 horizontal diffusion
+   USE lbclnk          ! Lateral Boundary condition - MPP exchanges
+   USE lib_mpp         ! MPP library
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_trp_2     ! called by sbc_ice_lim_2
-
-   !! * Shared module variables
-   REAL(wp), PUBLIC  ::   &  !:
-      bound  = 0.e0          !: boundary condit. (0.0 no-slip, 1.0 free-slip)
-
-   !! * Module variables
-   REAL(wp)  ::           &  ! constant values
-      epsi06 = 1.e-06  ,  &
-      epsi03 = 1.e-03  ,  &
-      epsi16 = 1.e-16  ,  &
-      rzero  = 0.e0    ,  &
-      rone   = 1.e0
+   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)  ::   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"
    !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !! NEMO/LIM2 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -62,34 +60,17 @@
       !!
       !! ** action :
-      !!
-      !! History :
-      !!   1.0  !  00-01 (LIM)  Original code
-      !!        !  01-05 (G. Madec, R. Hordoir) opa norm
-      !!   2.0  !  04-01 (G. Madec, C. Ethe)  F90, mpp
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! number of iteration
-
-      INTEGER  ::   ji, jj, jk,   &  ! dummy loop indices
-         &          initad           ! number of sub-timestep for the advection
-
-      REAL(wp) ::  &                              
-         zindb  ,  &
-         zacrith, &
-         zindsn , &
-         zindic , &
-         zusvosn, &
-         zusvoic, &
-         zignm  , &
-         zindhe , &
-         zvbord , &
-         zcfl   , &
-         zusnit , &
-         zrtt, ztsn, ztic1, ztic2
-
-      REAL(wp), DIMENSION(jpi,jpj)  ::   &  ! temporary workspace
-         zui_u , zvi_v , zsm   ,         &
-         zs0ice, zs0sn , zs0a  ,         &
-         zs0c0 , zs0c1 , zs0c2 ,         &
-         zs0st
+      !!
+      INTEGER  ::   ji, jj, jk   ! dummy loop indices
+      INTEGER  ::   initad       ! number of sub-timestep for the advection
+
+      REAL(wp) ::   zindb  , zacrith, zindsn , zindic , zusvosn   ! local scalars
+      REAL(wp) ::   zusvoic, zignm  , zindhe , zvbord , zcfl      !   -      -
+      REAL(wp) ::   zusnit , zrtt   , ztsn   , ztic1  , ztic2     !   -      -
+
+      REAL(wp), DIMENSION(jpi,jpj) ::   zui_u , zvi_v , zsm             ! 2D workspace
+      REAL(wp), DIMENSION(jpi,jpj) ::   zs0ice, zs0sn , zs0a            !  -      -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zs0c0 , zs0c1 , zs0c2 , zs0st   !  -      -
       !---------------------------------------------------------------------
 
@@ -107,13 +88,16 @@
          ! zvbord factor between 1 and 2 to take into account slip or no-slip boundary conditions.        
          zvbord = 1.0 + ( 1.0 - bound )
-         DO jj = 1, jpjm1
+#if defined key_lim2_vp
+         DO jj = 1, jpjm1     ! VP rheology: ice (u,v) at I-point
             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 ) )
+               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
-         ! Lateral boundary conditions on zui_u, zvi_v
-         CALL lbc_lnk( zui_u, 'U', -1. )
-         CALL lbc_lnk( zvi_v, 'V', -1. )
+         CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )      ! Lateral boundary conditions 
+#else
+        zui_u(:,:) = u_ice(:,:)      ! EVP rheology: ice (u,v) at u- and v-points
+        zvi_v(:,:) = v_ice(:,:)
+#endif
 
          ! CFL test for stability
@@ -123,18 +107,17 @@
          zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) )
 
-         IF (lk_mpp ) CALL mpp_max(zcfl)
-
-         IF ( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl
+         IF(lk_mpp ) CALL mpp_max(zcfl)
+
+         IF( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : cfl criterion violation. day ',nday,' cfl = ',zcfl
 
          ! content of properties
          ! ---------------------
-         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.
-         zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)  ! Heat content of the first ice layer.
-         zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)  ! Heat content of the second ice layer.
-         zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a(:,:)    ! Heat reservoir for brine pockets.
-         
+         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.
+         zs0c1 (:,:) =  tbif(:,:,2) / rt0_ice  * zs0ice(:,:)   ! Heat content of the first ice layer.
+         zs0c2 (:,:) =  tbif(:,:,3) / rt0_ice  * zs0ice(:,:)   ! Heat content of the second ice layer.
+         zs0st (:,:) =  qstoif(:,:) / xlic     * zs0a  (:,:)   ! Heat reservoir for brine pockets.
  
          ! Advection 
@@ -144,5 +127,7 @@
          zusnit = 1.0 / REAL( initad ) 
          
-         IF ( MOD( nday , 2 ) == 0) THEN
+!!gm this has been changed in the reference to become odd and even ice time step
+
+         IF( MOD( nday , 2 ) == 0) THEN
             DO jk = 1,initad
                CALL lim_adv_x_2( zusnit, zui_u, rone , zsm, zs0ice, sxice, sxxice, syice, syyice, sxyice )
@@ -231,5 +216,4 @@
          ! Up-dating and limitation of sea ice properties after transport.
          DO jj = 1, jpj
-!!!iii      zindhe = REAL( MAX( 0, isign(1, jj - njeq ) ) )              !ibug mpp  !!bugmpp  njeq!
             zindhe = MAX( 0.e0, SIGN( 1.e0, fcor(1,jj) ) )              ! = 0 for SH, =1 for NH
             DO ji = 1, jpi
@@ -272,7 +256,7 @@
             END DO
          END DO
-         
+         !
       ENDIF
-      
+      !
    END SUBROUTINE lim_trp_2
 
@@ -288,11 +272,8 @@
       !!
       !! ** input   :   Namelist namicetrp
-      !!
-      !! history :
-      !!   2.0  !  03-08 (C. Ethe)  Original code
       !!-------------------------------------------------------------------
       NAMELIST/namicetrp/ bound
       !!-------------------------------------------------------------------
-
+      !
       ! Read Namelist namicetrp
       REWIND ( numnam_ice )
@@ -304,5 +285,5 @@
          WRITE(numout,*) '   boundary conditions (0. no-slip, 1. free-slip) bound  = ', bound
       ENDIF
-            
+      !     
    END SUBROUTINE lim_trp_init_2
 
Index: branches/devmercator2010_1/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- branches/devmercator2010_1/NEMO/LIM_SRC_3/limrhg.F90	(revision 2134)
+++ branches/devmercator2010_1/NEMO/LIM_SRC_3/limrhg.F90	(revision 2135)
@@ -4,9 +4,10 @@
    !!   Ice rheology : sea ice rheology
    !!======================================================================
-   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
-   !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3
+   !! History :  LIM  !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
+   !!            3.0  !  2008-03  (M. Vancoppenolle) LIM 3
    !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy 
+   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas
    !!----------------------------------------------------------------------
-#if defined key_lim3
+#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
    !!----------------------------------------------------------------------
    !!   'key_lim3'                                      LIM3 sea-ice model
@@ -14,28 +15,29 @@
    !!   lim_rhg   : computes ice velocities
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE phycst
-   USE par_oce
-   USE dom_oce
-   USE dom_ice
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE sbc_ice         ! Surface boundary condition: ice fields
-   USE ice
-   USE iceini
-   USE lbclnk
-   USE lib_mpp
-   USE in_out_manager  ! I/O manager
-   USE limitd_me
-   USE prtctl          ! Print control
-
+   USE phycst           ! physical constants
+   USE par_oce          ! ocean parameters
+   USE dom_oce          ! ocean domain
+   USE sbc_oce          ! Surface boundary condition: ocean fields
+   USE sbc_ice          ! Surface boundary condition: ice fields
+   USE lbclnk           ! lateral boundary condition - MPP exchanges
+   USE lib_mpp          ! MPP library
+   USE in_out_manager   ! I/O manager
+   USE limitd_me        ! LIM3: 
+   USE prtctl           ! control print
+#if defined key_lim3
+   USE ice              ! LIM3: ice variables
+   USE dom_ice          ! LIM3: ice domain
+   USE iceini           ! LIM3: ice initialisation
+#endif
+#if defined key_lim2 && ! defined key_lim2_vp
+   USE ice_2            ! LIM2: ice variables
+   USE dom_ice_2        ! LIM2: ice domain
+   USE iceini_2         ! LIM2: ice initialisation
+#endif
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_rhg  ! routine called by lim_dyn
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
+   PUBLIC   lim_rhg   ! routine called by lim_dyn module
 
    !! * Module variables
@@ -43,14 +45,15 @@
       rzero   = 0.e0   ,  &
       rone    = 1.e0
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008) 
+   !! NEMO/LIM-3 3.3,  UCL-LOCEAN-IPSL (2010) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
    SUBROUTINE lim_rhg( k_j1, k_jpj )
-
       !!-------------------------------------------------------------------
       !!                 ***  SUBROUTINE lim_rhg  ***
@@ -100,48 +103,23 @@
       !!
       !! ** References : Hunke and Dukowicz, JPO97
-      !!                 Bouillon et al., 08, in prep (update this when
-      !!                 published)
-      !!                 Vancoppenolle et al., OM08
+      !!                 Bouillon et al., 2009, Ocean. Modelling, 27, 174-184.
+      !!                 Vancoppenolle et al. 2009, Ocean Modelling, 27, 33-53.
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
+      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation
       !!
-      !!-------------------------------------------------------------------
-      ! * Arguments
-      !
-      INTEGER, INTENT(in) :: &
-         k_j1 ,                      & !: southern j-index for ice computation
-         k_jpj                         !: northern j-index for ice computation
-
-      ! * Local variables
-      INTEGER ::   ji, jj              !: dummy loop indices
-
-      INTEGER  :: &
-         jter                          !: temporary integers
-
-      CHARACTER (len=50) ::   charout
-
-      REAL(wp) :: &
-         zt11, zt12, zt21, zt22,     & !: temporary scalars
-         ztagnx, ztagny,             & !: wind stress on U/V points                       
-         delta                         !
-
-      REAL(wp) :: &
-         za,                         & !:
-         zstms,                      & !: temporary scalar for ice strength
-         zsang,                      & !: temporary scalar for coriolis term
-         zmask                         !: mask for the computation of ice mass
-
-      REAL(wp),DIMENSION(jpi,jpj) :: &
-         zpresh        ,             & !: temporary array for ice strength
-         zpreshc       ,             & !: Ice strength on grid cell corners (zpreshc)
-         zfrld1, zfrld2,             & !: lead fraction on U/V points                                    
-         zmass1, zmass2,             & !: ice/snow mass on U/V points                                    
-         zcorl1, zcorl2,             & !: coriolis parameter on U/V points
-         za1ct, za2ct  ,             & !: temporary arrays
-         zc1           ,             & !: ice mass
-         zusw          ,             & !: temporary weight for the computation
-                                !: of ice strength
-         u_oce1, v_oce1,             & !: ocean u/v component on U points                           
-         u_oce2, v_oce2,             & !: ocean u/v component on V points
-         u_ice2,                     & !: ice u component on V point
-         v_ice1                        !: ice v component on U point
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      INTEGER  ::   jter     ! local integers
+      CHARACTER (len=50) ::   charout   ! local character
+      REAL(wp) ::   zt11, zt12, zt21, zt22   ! local scalars
+      REAL(wp) ::   ztagnx, ztagny, delta    !   -      -
+      REAL(wp) ::   za, zstms, zsang, zmask  !   -      -
+      REAL(wp) ::   zresm, zindb, zdummy     !   -    -
+
+      REAL(wp),DIMENSION(jpi,jpj) ::   zpresh , zfrld1, zmass1, zcorl1, za1ct    ! 2D workspace
+      REAL(wp),DIMENSION(jpi,jpj) ::   zpreshc, zfrld2, zmass2, zcorl2, za2ct    !  -      -
+      REAL(wp),DIMENSION(jpi,jpj) ::   u_oce1, v_oce1, u_ice2, zc1               !  -      -                          
+      REAL(wp),DIMENSION(jpi,jpj) ::   u_oce2, v_oce2, v_ice1, zusw              !  -      -
+      REAL(wp),DIMENSION(jpi,jpj) ::   zf1, zf2                                  !  -      -
 
       REAL(wp) :: &
@@ -159,6 +137,4 @@
          sigma1, sigma2                !: internal ice stress
 
-      REAL(wp),DIMENSION(jpi,jpj) :: &
-         zf1, zf2                      !: arrays for internal stresses
 
       REAL(wp),DIMENSION(jpi,jpj) :: &
@@ -170,14 +146,13 @@
          zs12                          !: Non-diagonal stress tensor component zs12
 
-      REAL(wp) :: &
-         zresm            ,          & !: Maximal error on ice velocity
-         zindb            ,          & !: ice (1) or not (0)      
-         zdummy                        !: dummy argument
-
-      REAL(wp),DIMENSION(jpi,jpj) :: &
-         zu_ice           ,          & !: Ice velocity on previous time step
-         zv_ice           ,          &
-         zresr                         !: Local error on velocity
-
+
+      REAL(wp),DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr   !
+      !!-------------------------------------------------------------------
+
+#if  defined key_lim2 && ! defined key_lim2_vp
+     vt_s => hsnm
+     vt_i => hicm
+     at_i(:,:) = 1. - frld(:,:)
+#endif
       !
       !------------------------------------------------------------------------------!
@@ -190,7 +165,8 @@
       u_ice2(:,:)  = 0.0 ; v_ice1(:,:)  = 0.0
       zdd(:,:)     = 0.0 ; zdt(:,:)     = 0.0 ; zds(:,:)     = 0.0
-
+#if defined key_lim3
       ! Ice strength on T-points
       CALL lim_itd_me_icestrength(ridge_scheme_swi)
+#endif
 
       ! Ice mass and temp variables
@@ -200,5 +176,9 @@
          DO ji = 1 , jpi
             zc1(ji,jj)    = tms(ji,jj) * ( rhosn * vt_s(ji,jj) + rhoic * vt_i(ji,jj) )
+#if defined key_lim3
             zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) / 2.
+#else
+            zpresh(ji,jj) = tms(ji,jj) *  2. * pstar * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) )
+#endif
             ! tmi = 1 where there is ice or on land
             tmi(ji,jj)    = 1.0 - ( 1.0 - MAX( 0.0 , SIGN ( 1.0 , vt_i(ji,jj) - &
@@ -269,4 +249,5 @@
                / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd )
             !
+            ! Mass, coriolis coeff. and currents
             u_oce1(ji,jj)  = u_oce(ji,jj)
             v_oce2(ji,jj)  = v_oce(ji,jj)
Index: branches/devmercator2010_1/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90
===================================================================
--- branches/devmercator2010_1/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 2134)
+++ branches/devmercator2010_1/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 2135)
@@ -8,4 +8,5 @@
    !! History :  1.0   !  06-2006  (G. Madec)  from icestp_2.F90
    !!            3.0   !  08-2008  (S. Masson, E. .... ) coupled interface
+   !!            3.3   !  05-2009  (G.Garric) addition of the lim2_evp case
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -53,6 +54,4 @@
    PUBLIC sbc_ice_lim_2 ! routine called by sbcmod.F90
    
-   CHARACTER(len=1) ::   cl_grid = 'B'     ! type of grid used in ice dynamics
-
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -172,7 +171,8 @@
          !  Ice model step  !
          ! ---------------- !
-
-                                        CALL lim_rst_opn_2  ( kt )      ! Open Ice restart file
-         IF( .NOT. lk_c1d ) THEN                                        ! Ice dynamics & transport (not in 1D case)
+         numit = numit + nn_fsbc                                             ! Ice model time step
+
+                                        CALL lim_rst_opn_2  ( kt )           ! Open Ice restart file
+         IF( .NOT. lk_c1d ) THEN                                             ! Ice dynamics & transport (not in 1D case)
                                         CALL lim_dyn_2      ( kt )           ! Ice dynamics    ( rheology/dynamics )
                                         CALL lim_trp_2      ( kt )           ! Ice transport   ( Advection/diffusion )
