Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/dom_ice_2.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/dom_ice_2.F90	(revision 2318)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/dom_ice_2.F90	(revision 2319)
@@ -5,6 +5,9 @@
    !!======================================================================
    !! History :   2.0  !  03-08  (C. Ethe)  Free form and module
+   !!             3.3  !  2009-05 (G.Garric, C. Bricaud) addition of lim2_evp case
    !!----------------------------------------------------------------------
 #if defined key_lim2
+   !!----------------------------------------------------------------------
+   !!   'key_lim2'                                       LIM2 sea-ice model
    !!----------------------------------------------------------------------
    !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
@@ -22,13 +25,24 @@
       !                                        !  (otherwise = jpj+10 (SH) or -10 (SH) )
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor , fcor,   &  !: coriolis factor and coeficient
-      &                                              covrai ,         &  !: sine of geographic latitude
-      &                                              area   ,         &  !: surface of grid cell 
-      &                                              tms    , tmu        !: temperature and velocity points masks
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght   ,         &  !: weight of the 4 neighbours to compute averages
-      &                                              akappa , bkappa     !: first and third group of metric coefficients
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd   !: second group of metric coefficients
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   fs2cor , fcor     !: coriolis factor and coeficient
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   covrai            !: sine of geographic latitude
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   area              !: surface of grid cell 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tms    , tmu      !: temperature and velocity points masks
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   wght              !: weight of the 4 neighbours to compute averages
 
+
+# if defined key_lim2_vp
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2)     ::   akappa , bkappa   !: first and third group of metric coefficients
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,2,2,2,2) ::   alambd            !: second group of metric coefficients
+# else
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tmv    , tmf      !: y-velocity and F-points masks
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj)         ::   tmi               !: ice mask: =1 if ice thick > 0
+# endif
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option          Empty module         NO LIM-2 sea-ice model
+   !!----------------------------------------------------------------------
+#endif
    !!======================================================================
-#endif
 END MODULE dom_ice_2
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 2318)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 2319)
@@ -5,4 +5,5 @@
    !!=====================================================================
    !! History :  2.0  !  03-08  (C. Ethe)  F90: Free form and module
+   !!            3.3  !  2009-05  (G.Garric) addition of the lim2_evp cas
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -25,4 +26,5 @@
    LOGICAL               , PUBLIC ::   ln_limdyn     = .TRUE.             !: flag for ice dynamics (T) or not (F)
    LOGICAL               , PUBLIC ::   ln_limdmp     = .FALSE.            !: Ice damping
+   LOGICAL               , PUBLIC ::   ln_nicep      = .TRUE.             !: flag grid points output (T) or not (F)
    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)
@@ -32,4 +34,6 @@
    INTEGER , PUBLIC ::   nbiter = 1         !: number of sub-time steps for relaxation
    INTEGER , PUBLIC ::   nbitdr = 250       !: maximum number of iterations for relaxation
+   INTEGER , PUBLIC ::   nevp   =   360     !: number of EVP subcycling iterations
+   INTEGER , PUBLIC ::   telast =  3600     !: timescale for EVP elastic waves
    REAL(wp), PUBLIC ::   rdt_ice            !: ice time step
    REAL(wp), PUBLIC ::   epsd   = 1.0e-20   !: tolerance parameter for dynamic
@@ -46,4 +50,5 @@
    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)
+   REAL(wp), PUBLIC ::   alphaevp = 1.e0    !: coefficient for the solution of EVP int. stresses
 
    REAL(wp), PUBLIC ::   usecc2             !:  = 1.0 / ( ecc * ecc )
@@ -54,9 +59,33 @@
    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
+
+   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
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90	(revision 2318)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90	(revision 2319)
@@ -6,4 +6,5 @@
    !! History :   1.0  !  02-08  (G. Madec)  F90: Free form and modules
    !!             2.0  !  03-08  (C. Ethe)  add ice_run
+   !!             3.3  !  09-05  (G.Garric, C. Bricaud) addition of the lim2_evp case
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -15,14 +16,14 @@
    !!   ice_run_2        : Definition some run parameter for ice model
    !!----------------------------------------------------------------------
-   USE dom_oce
-   USE dom_ice_2
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE sbc_ice         ! surface boundary condition: ice
-   USE phycst          ! Define parameters for the routines
-   USE ice_2
-   USE limmsh_2
-   USE limistate_2
-   USE limrst_2   
-   USE in_out_manager
+   USE dom_oce          ! ocean domain
+   USE dom_ice_2        ! LIM2: ice domain
+   USE sbc_oce          ! surface boundary condition: ocean
+   USE sbc_ice          ! surface boundary condition: ice
+   USE phycst           ! Define parameters for the routines
+   USE ice_2            ! LIM2: ice variable
+   USE limmsh_2         ! LIM2: mesh
+   USE limistate_2      ! LIM2: initial state
+   USE limrst_2         ! LIM2: restart
+   USE in_out_manager   ! I/O manager
       
    IMPLICIT NONE
@@ -30,4 +31,6 @@
 
    PUBLIC   ice_init_2               ! called by sbcice_lim_2.F90
+
+   INTEGER, PUBLIC ::   numit   !: iteration number
 
    !!----------------------------------------------------------------------
@@ -43,5 +46,5 @@
       !!                  ***  ROUTINE ice_init_2  ***
       !!
-      !! ** purpose :   
+      !! ** purpose :   initialisation of LIM-2 domain and variables  
       !!----------------------------------------------------------------------
       !
@@ -52,4 +55,5 @@
       ! Louvain la Neuve Ice model
       rdt_ice = nn_fsbc * rdttra(1)
+      numit   = nit000 - 1
 
       CALL lim_msh_2                  ! ice mesh initialization
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 2318)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 2319)
@@ -22,6 +22,9 @@
    USE dom_ice_2      !
    USE limistate_2    !
+#if defined key_lim2_vp
    USE limrhg_2       ! ice rheology
-
+#else
+   USE limrhg         ! LIM : EVP ice rheology
+#endif
    USE lbclnk         !
    USE lib_mpp        !
@@ -87,5 +90,9 @@
             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
@@ -105,4 +112,5 @@
                   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
@@ -110,4 +118,10 @@
                CALL lim_rhg_2( i_j1, i_jpj )
                ! 
+#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
                ! Southern hemisphere
                i_j1  =  1 
@@ -116,4 +130,5 @@
                   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
@@ -121,4 +136,11 @@
                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 )
+               ! 
+#endif
+
             ELSE                                 ! local domain extends over one hemisphere only
                !                                 ! Rheology is computed only over the ice cover
@@ -138,5 +160,12 @@
                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 )
+#if defined key_lim2_vp
+               i_jpj = MIN( jpj, i_jpj+2)
+               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
@@ -148,12 +177,20 @@
          ! computation of friction velocity
          ! --------------------------------
-         ! ice-ocean velocity at U & V-points (u_ice v_ice at I-point ; ssu_m, ssv_m at U- & V-points)
-         
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1   ! NO vector opt.
-               zu_io(ji,jj) = 0.5 * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj)
-               zv_io(ji,jj) = 0.5 * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj)
-            END DO
-         END DO
+         SELECT CASE( cl_grid )
+         CASE( 'C' )       ! C-grid ice dynamics (EVP)
+                           ! ice-ocean & ice velocity at  ocean velocity points
+            zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)  
+            zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
+            !
+         CASE( 'B' )       ! B-grid ice dynamics (VP) 
+                           ! 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
+         END SELECT
+
          ! frictional velocity at T-point
          DO jj = 2, jpjm1
@@ -198,5 +235,6 @@
       NAMELIST/namicedyn/ epsd, alpha,     &
          &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   &
-         &                c_rhg, etamn, creepl, ecc, ahi0
+         &                c_rhg, etamn, creepl, ecc, ahi0                   &
+         &                nevp, telast,alphaevp
       !!-------------------------------------------------------------------
 
@@ -223,4 +261,7 @@
          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
 
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 2318)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 2319)
@@ -50,12 +50,12 @@
       !! * Local variables
       INTEGER :: ji, jj      ! dummy loop indices
-
-      REAL(wp), DIMENSION(jpi,jpj) ::  &
-         zd2d1 , zd1d2       ! Derivative of zh2 (resp. zh1) in the x direction
-         !                   ! (resp. y direction) (defined at the center)
-      REAL(wp) ::         &
-         zh1p  , zh2p   , &  ! Idem zh1, zh2 for the bottom left corner of the grid
-         zd2d1p, zd1d2p , &  ! Idem zd2d1, zd1d2 for the bottom left corner of the grid
-         zusden, zusden2     ! temporary scalars
+      REAL(wp) ::   zusden   ! local scalars
+
+#if defined key_lim2_vp
+      REAL(wp) ::   zusden2           ! local scalars
+      REAL(wp) ::   zh1p  , zh2p      !   -      -
+      REAL(wp) ::   zd2d1p, zd1d2p    !   -      -
+      REAL(wp), DIMENSION(jpi,jpj) ::   zd2d1 , zd1d2   ! 2D workspace
+#endif
       !!---------------------------------------------------------------------
 
@@ -76,5 +76,5 @@
       njeqm1 = njeq - 1 
 
-      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor
+      fcor(:,:) = 2. * omega * SIN( gphit(:,:) * rad )   !  coriolis factor at T-point
  
 !i    DO jj = 1, jpj
@@ -115,11 +115,17 @@
       !-------------------
 !!ibug ???
-      akappa(:,:,:,:) = 0.e0
       wght(:,:,:,:) = 0.e0
+      tmu(:,:)      = 0.e0
+#if defined key_lim2_vp 
+      akappa(:,:,:,:)     = 0.e0
       alambd(:,:,:,:,:,:) = 0.e0
-      tmu(:,:) = 0.e0
+#else
+      tmv(:,:) = 0.e0
+      tmf(:,:) = 0.e0
+#endif
 !!i
       
-      
+
+#if defined key_lim2_vp      
       ! metric coefficients for sea ice dynamic
       !----------------------------------------
@@ -155,8 +161,25 @@
       CALL lbc_lnk( wght(:,:,2,1), 'I', 1. )      ! but it is never used
       CALL lbc_lnk( wght(:,:,2,2), 'I', 1. )
+#else
+      ! metric coefficients for sea ice dynamic (EVP rheology)
+      !----------------------------------------
+      DO jj = 1, jpjm1                                       ! weights (wght) at F-points
+         DO ji = 1, jpim1
+            zusden = 1. / (  ( e1t(ji+1,jj  ) + e1t(ji,jj) )   &
+               &           * ( e2t(ji  ,jj+1) + e2t(ji,jj) ) ) 
+            wght(ji,jj,1,1) = zusden * e1t(ji+1,jj) * e2t(ji,jj+1)
+            wght(ji,jj,1,2) = zusden * e1t(ji+1,jj) * e2t(ji,jj  )
+            wght(ji,jj,2,1) = zusden * e1t(ji  ,jj) * e2t(ji,jj+1)
+            wght(ji,jj,2,2) = zusden * e1t(ji  ,jj) * e2t(ji,jj  )
+         END DO
+      END DO
+      CALL lbc_lnk( wght(:,:,1,1), 'F', 1. )   ;   CALL lbc_lnk( wght(:,:,1,2),'F', 1. )       ! lateral boundary cond.   
+      CALL lbc_lnk( wght(:,:,2,1), 'F', 1. )   ;   CALL lbc_lnk( wght(:,:,2,2),'F', 1. )
+#endif
     
       ! Coefficients for divergence of the stress tensor
       !-------------------------------------------------
 
+#if defined key_lim2_vp
       DO jj = 2, jpj
          DO ji = 2, jpi   ! NO vector opt.
@@ -226,4 +249,5 @@
       CALL lbc_lnk( alambd(:,:,2,1,1,1), 'I', 1. )      !
       CALL lbc_lnk( alambd(:,:,2,1,1,2), 'I', 1. )      !
+#endif
             
 
@@ -233,4 +257,6 @@
       tms(:,:) = tmask(:,:,1)      ! ice T-point  : use surface tmask
 
+#if defined key_lim2_vp
+      ! VP rheology : ice velocity point is I-point
 !i here we can use umask with a i and j shift of -1,-1
       tmu(:,1) = 0.e0
@@ -244,4 +270,12 @@
       !--lateral boundary conditions    
       CALL lbc_lnk( tmu(:,:), 'I', 1. )
+#else
+      ! EVP rheology : ice velocity point are U- & V-points ; ice vorticity
+      ! point is F-point
+      tmu(:,:) = umask(:,:,1)
+      tmv(:,:) = vmask(:,:,1)
+      tmf(:,:) = 0.e0                        ! used of fmask except its special value along the coast (rn_shlat)
+      WHERE( fmask(:,:,1) == 1.e0 )   tmf(:,:) = 1.e0
+#endif
       
       ! unmasked and masked area of T-grid cell
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limrhg.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limrhg.F90	(revision 2319)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limrhg.F90	(revision 2319)
@@ -0,0 +1,1 @@
+link ../LIM_SRC_3/limrhg.F90
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 2318)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 2319)
@@ -4,13 +4,15 @@
    !!   Ice rheology :  performs sea ice rheology
    !!======================================================================
-   !! History :  0.0  !  93-12  (M.A. Morales Maqueda.)  Original code
-   !!            1.0  !  94-12  (H. Goosse) 
-   !!            2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp
-   !!            " "  !  06-08  (G. Madec)  surface module, ice-stress at I-point
-   !!            " "  !  09-09  (G. Madec)  Huge verctor optimisation
-   !!----------------------------------------------------------------------
-#if defined key_lim2
-   !!----------------------------------------------------------------------
-   !!   'key_lim2'                                    LIM 2.0 sea-ice model
+   !! History :  0.0  !  1993-12  (M.A. Morales Maqueda.)  Original code
+   !!            1.0  !  1994-12  (H. Goosse) 
+   !!            2.0  !  2003-12  (C. Ethe, G. Madec)  F90, mpp
+   !!            " "  !  2006-08  (G. Madec)  surface module, ice-stress at I-point
+   !!            " "  !  2009-09  (G. Madec)  Huge verctor optimisation
+   !!            3.3  !  2009-05  (G.Garric, C. Bricaud) addition of the lim2_evp case
+   !!----------------------------------------------------------------------
+#if defined   key_lim2   &&   defined key_lim2_vp
+   !!----------------------------------------------------------------------
+   !!   'key_lim2'                and                 LIM 2.0 sea-ice model
+   !!   'key_lim2_vp'                                       VP ice rheology
    !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
@@ -21,8 +23,8 @@
    USE sbc_oce        ! surface boundary condition: ocean variables
    USE sbc_ice        ! surface boundary condition: ice variables
-   USE dom_ice_2      ! domaine: ice variables
+   USE dom_ice_2      ! LIM2: ice domain
    USE phycst         ! physical constant
-   USE ice_2          ! ice variables
-   USE lbclnk         ! lateral boundary condition
+   USE ice_2          ! LIM2: ice variables
+   USE lbclnk         ! lateral boundary condition - MPP exchanges
    USE lib_mpp        ! MPP library
    USE in_out_manager ! I/O manager
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limrst_2.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limrst_2.F90	(revision 2318)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limrst_2.F90	(revision 2319)
@@ -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/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 2318)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 2319)
@@ -4,7 +4,8 @@
    !!           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
@@ -17,10 +18,10 @@
    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
@@ -35,9 +36,10 @@
    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
+   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
@@ -78,31 +80,29 @@
       !!
       INTEGER  ::   ji, jj           ! dummy loop indices
-      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
-      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
-      REAL(wp) ::   zrdtir           ! 1. / rdt_ice
-      REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux
-      REAL(wp) ::   zinda            ! switch for testing the values of ice concentration
-      REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface
-      REAL(wp) ::   zemp             ! freshwater exchanges at the ice/ocean interface
-      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points
-      REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points
-      REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity
-! interface 2D --> 3D
-      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb     ! albedo of ice under overcast sky
-      REAL(wp), DIMENSION(jpi,jpj,1) ::   zalbp    ! albedo of ice under clear sky
-      REAL(wp) ::   zsang, zmod, zztmp, zfm
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! component of ocean stress below sea-ice at I-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztiomi           ! module    of ocean stress below sea-ice at I-point
-      REAL(wp), DIMENSION(jpi,jpj) ::   zqnsoce          ! save qns before its modification by ice model
+      INTEGER  ::   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
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition'
+#if defined key_lim2_vp
+         IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice (VP rheology)  - surface boundary condition'
+#else
+         IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice (EVP rheology) - surface boundary condition'
+#endif
          IF(lwp) WRITE(numout,*) '~~~~~~~~~   '
-
+         !
+         r1_rdtice = 1. / rdt_ice
+         !
          soce_r(:,:) = soce
          sice_r(:,:) = sice
@@ -186,6 +186,6 @@
             zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads
                &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            &
-               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * zrdtir    &
-               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) * zrdtir
+               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice    &
+               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice 
 
             fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ???
@@ -198,5 +198,5 @@
       CALL iom_put( 'hflx_ice_cea', - fdtcn(:,:) )      
       CALL iom_put( 'qns_io_cea', qns(:,:) - zqnsoce(:,:) * pfrld(:,:) )      
-      CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1. - pfrld(:,:)) )
+      CALL iom_put( 'qsr_io_cea', fstric(:,:) * (1.e0 - pfrld(:,:)) )
 
       !------------------------------------------!
@@ -212,22 +212,18 @@
             
 #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 
-#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 
+            ! freshwater exchanges at the ice-atmosphere / ocean interface (coupled mode)
+            zemp = emp_tot(ji,jj) - emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   ! 
+               &   + rdmsnif(ji,jj) * r1_rdtice                                   !  freshwaterflux due to snow melting 
+#else
             !  computing freshwater exchanges at the ice/ocean interface
             zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction 
                &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean
                &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step
-               &   + rdmsnif(ji,jj) * zrdtir                       !  freshwaterflux due to snow melting 
+               &   + rdmsnif(ji,jj) * r1_rdtice                    !  freshwaterflux due to snow melting 
                !                                                   !  ice-covered fraction:
 #endif            
 
             !  computing salt exchanges at the ice/ocean interface
-            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * zrdtir ) 
+            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice ) 
             
             !  converting the salt flux from ice to a freshwater flux from ocean
@@ -241,7 +237,7 @@
 
       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
 
@@ -277,6 +273,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  ) )
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 2318)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 2319)
@@ -7,4 +7,5 @@
    !!            2.0  !  2001-05 (G. Madec, R. Hordoir) opa norm
    !!             -   !  2004-01 (G. Madec, C. Ethe)  F90, mpp
+   !!            3.3  !  2009-05  (G.Garric, C. Bricaud) addition of the lim2_evp case
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -32,12 +33,12 @@
    PUBLIC   lim_trp_2   ! called by sbc_ice_lim_2
 
-   REAL(wp), PUBLIC  ::   bound  = 0.e0   !: boundary condit. (0.0 no-slip, 1.0 free-slip)
-
-   REAL(wp)  ::           &  ! constant values
-      epsi06 = 1.e-06  ,  &
-      epsi03 = 1.e-03  ,  &
-      epsi16 = 1.e-16  ,  &
-      rzero  = 0.e0    ,  &
-      rone   = 1.e0
+   REAL(wp), PUBLIC ::   bound  = 0.e0          !: boundary condit. (0.0 no-slip, 1.0 free-slip)
+
+   REAL(wp)  ::   epsi06 = 1.e-06   ! constant values
+   REAL(wp)  ::   epsi03 = 1.e-03  
+   REAL(wp)  ::   epsi16 = 1.e-16  
+   REAL(wp)  ::   rzero  = 0.e0   
+   REAL(wp)  ::   rone   = 1.e0
+
 
    !! * Substitution
@@ -88,4 +89,5 @@
          ! ---------------------------------------
          zvbord = 1.0 + ( 1.0 - bound )      ! zvbord=2 no-slip, =0 free slip boundary conditions        
+#if defined key_lim2_vp
          DO jj = 1, jpjm1
             DO ji = 1, jpim1   ! NO vector opt.
@@ -95,5 +97,8 @@
          END DO
          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
@@ -109,6 +114,6 @@
          ! content of properties
          ! ---------------------
-         zs0sn (:,:) =  hsnm(:,:) * area(:,:)                 ! Snow volume.
-         zs0ice(:,:) =  hicm(:,:) * area(:,:)                 ! Ice volume.
+         zs0sn (:,:) =  hsnm(:,:)              * area  (:,:)  ! Snow volume.
+         zs0ice(:,:) =  hicm(:,:)              * area  (:,:)  ! Ice volume.
          zs0a  (:,:) =  ( 1.0 - frld(:,:) )    * area  (:,:)  ! Surface covered by ice.
          zs0c0 (:,:) =  tbif(:,:,1) / rt0_snow * zs0sn (:,:)  ! Heat content of the snow layer.
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 2318)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 2319)
@@ -7,6 +7,7 @@
    !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3
    !!             -   !  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
@@ -18,9 +19,6 @@
    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
@@ -28,4 +26,14 @@
    USE limitd_me
    USE prtctl          ! Print control
+#if defined key_lim3
+   USE ice
+   USE dom_ice
+   USE iceini
+#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
 
 
@@ -179,5 +187,11 @@
          zv_ice           ,          &
          zresr                         !: Local error on velocity
-
+      !!-------------------------------------------------------------------
+
+#if  defined key_lim2 && ! defined key_lim2_vp
+     vt_s => hsnm
+     vt_i => hicm
+     at_i(:,:) = 1. - frld(:,:)
+#endif
       !
       !------------------------------------------------------------------------------!
@@ -191,6 +205,8 @@
       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 +216,7 @@
          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.
+#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) - &
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 2318)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 2319)
@@ -53,6 +53,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,5 +170,7 @@
          !  Ice model step  !
          ! ---------------- !
-                                        CALL lim_rst_opn_2  ( kt )      ! Open Ice restart file
+         numit = numit + nn_fsbc                                             !Ice model time step
+
+                                        CALL lim_rst_opn_2  ( kt )           ! Open Ice restart file
 #if ! defined key_c1d
             ! Ice dynamics & transport (not in 1D case)
