Index: /branches/nemo_v3_3_beta/NEMOGCM/ARCH/arch-ifort_osx.fcm
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/ARCH/arch-ifort_osx.fcm	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/ARCH/arch-ifort_osx.fcm	(revision 2370)
@@ -15,6 +15,6 @@
 
 
-%NCDF_INC            -I/usr/local/pub/netcdf/4.0.1-ifort/include 
-%NCDF_LIB            -L /usr/local/pub/netcdf/4.0.1-ifort/lib -lnetcdf   
+%NCDF_INC            -I/usr/local/include 
+%NCDF_LIB            -L /usr/local/lib -lnetcdf   
 %FC                  mpif90
 %FCFLAGS 	     -r8 -O3  -traceback 
Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE/EXP00/namelist
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE/EXP00/namelist	(revision 2370)
@@ -131,7 +131,4 @@
                            !  =1 use observed ice-cover      ,
                            !  =2 ice-model used                             ("key_lim3" or "key_lim2)
-   nn_ico_cpl  = 0         !  ice-ocean coupling : =0 each nn_fsbc 
-                           !                       =1 stresses recomputed each ocean time step ("key_lim3" only)
-                           !                       =2 combination of 0 and 1 cases             ("key_lim3" only)
    ln_dm2dc    = .false.   !  daily mean to diurnal cycle short wave (qsr)
    ln_rnf      = .false.   !  runoffs (T => fill namsbc_rnf)
Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE_LOBSTER/EXP00/namelist
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE_LOBSTER/EXP00/namelist	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/GYRE_LOBSTER/EXP00/namelist	(revision 2370)
@@ -131,7 +131,4 @@
                            !  =1 use observed ice-cover      ,
                            !  =2 ice-model used                             ("key_lim3" or "key_lim2)
-   nn_ico_cpl  = 0         !  ice-ocean coupling : =0 each nn_fsbc 
-                           !                       =1 stresses recomputed each ocean time step ("key_lim3" only)
-                           !                       =2 combination of 0 and 1 cases             ("key_lim3" only)
    ln_dm2dc    = .false.   !  daily mean to diurnal cycle short wave (qsr)
    ln_rnf      = .false.   !  runoffs (T => fill namsbc_rnf)
Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/1_namelist	(revision 2370)
@@ -131,7 +131,4 @@
                            !  =1 use observed ice-cover      ,
                            !  =2 ice-model used                             ("key_lim3" or "key_lim2)
-   nn_ico_cpl  = 0         !  ice-ocean coupling : =0 each nn_fsbc 
-                           !                       =1 stresses recomputed each ocean time step ("key_lim3" only)
-                           !                       =2 combination of 0 and 1 cases             ("key_lim3" only)
    ln_dm2dc    = .false.   !  daily mean to diurnal cycle short wave (qsr)
    ln_rnf      = .false.   !  runoffs (T => fill namsbc_rnf)
Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/namelist	(revision 2370)
@@ -131,7 +131,4 @@
                            !  =1 use observed ice-cover      ,
                            !  =2 ice-model used                             ("key_lim3" or "key_lim2)
-   nn_ico_cpl  = 0         !  ice-ocean coupling : =0 each nn_fsbc 
-                           !                       =1 stresses recomputed each ocean time step ("key_lim3" only)
-                           !                       =2 combination of 0 and 1 cases             ("key_lim3" only)
    ln_dm2dc    = .false.   !  daily mean to diurnal cycle short wave (qsr)
    ln_rnf      = .true.    !  runoffs (T => fill namsbc_rnf)
Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/cpp_ORCA2_LIM.fcm
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/cpp_ORCA2_LIM.fcm	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM/cpp_ORCA2_LIM.fcm	(revision 2370)
@@ -1,1 +1,1 @@
- bld::tool::fppkeys key_trabbl key_vectopt_loop key_orca_r2 key_lim2 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_dtatem key_dtasal key_tradmp key_zdftke key_zdfddm key_iomput key_nproci=1 key_nprocj=1
+ bld::tool::fppkeys  key_trabbl key_vectopt_loop key_orca_r2 key_lim2 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_dtatem key_dtasal key_tradmp key_zdftke key_zdfddm key_iomput key_nproci=1 key_nprocj=1
Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_LIM_PISCES/EXP00/namelist	(revision 2370)
@@ -131,7 +131,4 @@
                            !  =1 use observed ice-cover      ,
                            !  =2 ice-model used                             ("key_lim3" or "key_lim2)
-   nn_ico_cpl  = 0         !  ice-ocean coupling : =0 each nn_fsbc 
-                           !                       =1 stresses recomputed each ocean time step ("key_lim3" only)
-                           !                       =2 combination of 0 and 1 cases             ("key_lim3" only)
    ln_dm2dc    = .false.   !  daily mean to diurnal cycle short wave (qsr)
    ln_rnf      = .true.    !  runoffs (T => fill namsbc_rnf)
Index: /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/CONFIG/ORCA2_OFF_PISCES/EXP00/namelist	(revision 2370)
@@ -131,7 +131,4 @@
                            !  =1 use observed ice-cover      ,
                            !  =2 ice-model used                             ("key_lim3" or "key_lim2)
-   nn_ico_cpl  = 0         !  ice-ocean coupling : =0 each nn_fsbc 
-                           !                       =1 stresses recomputed each ocean time step ("key_lim3" only)
-                           !                       =2 combination of 0 and 1 cases             ("key_lim3" only)
    ln_dm2dc    = .false.   !  daily mean to diurnal cycle short wave (qsr)
    ln_rnf      = .true.    !  runoffs (T => fill namsbc_rnf)
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 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90	(revision 2370)
@@ -4,5 +4,5 @@
    !! 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
    !!----------------------------------------------------------------------
@@ -11,9 +11,4 @@
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_ice_2          ! LIM sea-ice parameters
 
@@ -61,13 +56,15 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s         !: friction velocity
 
-   !!* diagnostic quantities
+   !!* Ice Rheology
 # if defined key_lim2_vp
    !                                                       !!* VP rheology *
+   LOGICAL , PUBLIC ::   lk_lim2_vp = .TRUE.               !: Visco-Plactic reology flag 
+   !
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm   !: mean snow and ice thicknesses
-   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'
+   LOGICAL , PUBLIC::   lk_lim2_vp = .FALSE.               !: Visco-Plactic reology flag 
+   !
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   stress1_i     !: first stress tensor element       
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   stress2_i     !: second stress tensor element
@@ -76,9 +73,8 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   divu_i        !: Divergence of the velocity field [s-1] -> limrhg.F90
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   shear_i       !: Shear of the velocity field [s-1] -> limrhg.F90
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   at_i          !: ice fraction
    !
-   REAL(wp), PUBLIC, DIMENSION(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
 
@@ -131,4 +127,8 @@
 #endif
 
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
+   !! $Id$
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE ice_2
Index: /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 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 2370)
@@ -4,8 +4,9 @@
    !!   Sea-Ice dynamics :  
    !!======================================================================
-   !! History :   1.0  !  01-04  (LIM)  Original code
-   !!             2.0  !  02-08  (C. Ethe, G. Madec)  F90, mpp
-   !!             2.0  !  03-08  (C. Ethe) add lim_dyn_init
-   !!             2.0  !  06-07  (G. Madec)  Surface module
+   !! History :  1.0  ! 2001-04  (LIM)  Original code
+   !!            2.0  ! 2002-08  (C. Ethe, G. Madec)  F90, mpp
+   !!            2.0  ! 2003-08  (C. Ethe) add lim_dyn_init
+   !!            2.0  ! 2006-07  (G. Madec)  Surface module
+   !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
    !!---------------------------------------------------------------------
 #if defined key_lim2
@@ -16,35 +17,30 @@
    !!    lim_dyn_init_2 : initialization and namelist read
    !!----------------------------------------------------------------------
-   USE dom_oce        ! ocean space and time domain
-   USE sbc_oce        !
-   USE phycst         !
-   USE ice_2          !
-   USE dom_ice_2      !
-   USE limistate_2    !
-#if defined key_lim2_vp
-   USE limrhg_2       ! ice rheology
-#else
-   USE limrhg         ! LIM : EVP ice rheology
-#endif
-   USE lbclnk         !
-   USE lib_mpp        !
-   USE in_out_manager ! I/O manager
-   USE prtctl         ! Print control
+   USE dom_oce          ! ocean space and time domain
+   USE sbc_oce          ! ocean surface boundary condition
+   USE phycst           ! physical constant
+   USE ice_2            ! LIM-2: ice variables
+   USE sbc_ice          ! Surface boundary condition: sea-ice fields
+   USE dom_ice_2        ! LIM-2: ice domain
+   USE limistate_2      ! LIM-2: initial state
+   USE limrhg_2         ! LIM-2: VP  ice rheology
+   USE limrhg           ! LIM  : EVP ice rheology
+   USE lbclnk           ! lateral boundary condition - MPP link
+   USE lib_mpp          ! MPP library
+   USE in_out_manager   ! I/O manager
+   USE prtctl           ! Print control
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   lim_dyn_2 ! routine called by sbc_ice_lim
-
-   !! * Module variables
-   REAL(wp)  ::  rone    = 1.e0   ! constant value
-
+   PUBLIC   lim_dyn_2   ! routine called by sbc_ice_lim
+
+   !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
    !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -90,9 +86,7 @@
             i_jpj = jpj
             IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=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
+            IF( lk_lim2_vp )   THEN   ;   CALL lim_rhg_2( i_j1, i_jpj )             !  VP rheology
+            ELSE                      ;   CALL lim_rhg  ( i_j1, i_jpj )             ! EVP rheology
+            ENDIF
             !
          ELSE                                 ! optimization of the computational area
@@ -112,16 +106,13 @@
                   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 )
-               ! 
-#else
-               i_j1 = MAX( 1, i_j1-2 )
+               IF( lk_lim2_vp )   THEN             ! VP  rheology
+                  i_j1 = MAX( 1, i_j1-1 )
+                  CALL lim_rhg_2( i_j1, i_jpj )
+               ELSE                                ! EVP rheology
+                  i_j1 = MAX( 1, i_j1-2 )
+                  CALL lim_rhg( i_j1, i_jpj )
+               ENDIF
                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, 'ij_jpj = ', i_jpj
-               CALL lim_rhg( i_j1, i_jpj )
-               ! 
-#endif
+               !
                ! Southern hemisphere
                i_j1  =  1 
@@ -130,17 +121,13 @@
                   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( lk_lim2_vp )   THEN             ! VP  rheology
+                  i_jpj = MIN( jpj, i_jpj+2 )
+                  CALL lim_rhg_2( i_j1, i_jpj )
+               ELSE                                ! EVP rheology
+                  i_jpj = MIN( jpj, i_jpj+1 )
+                  CALL lim_rhg( i_j1, i_jpj )
+               ENDIF
                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, 'ij_jpj = ', i_jpj
-               CALL lim_rhg( i_j1, i_jpj )
-               ! 
-#endif
-
+               !
             ELSE                                 ! local domain extends over one hemisphere only
                !                                 ! Rheology is computed only over the ice cover
@@ -156,16 +143,15 @@
                   i_jpj = i_jpj - 1
                END DO
-               i_jpj = MIN( jpj, i_jpj+2)
-    
+               i_jpj = MIN( jpj, i_jpj+2 )
+               ! 
+               IF( lk_lim2_vp )   THEN             ! VP  rheology
+                  i_jpj = MIN( jpj, i_jpj+2 )
+                  CALL lim_rhg_2( i_j1, i_jpj )                !  VP rheology
+               ELSE                                ! EVP rheology
+                  i_j1  = MAX( 1  , i_j1-2  )
+                  i_jpj = MIN( jpj, i_jpj+1 )
+                  CALL lim_rhg  ( i_j1, i_jpj )                ! EVP rheology
+               ENDIF
                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
-               ! 
-#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
@@ -177,16 +163,13 @@
          ! computation of friction velocity
          ! --------------------------------
-         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(:,:)  
+         SELECT CASE( cp_ice_msh )           ! ice-ocean relative velocity at u- & v-pts
+         CASE( 'C' )                               ! EVP : C-grid ice dynamics
+            zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)           ! ice-ocean & ice velocity at ocean velocity points
             zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
-            !
-         CASE( '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)
+         CASE( 'I' )                               ! VP  : B-grid ice dynamics (I-point) 
+            DO jj = 1, jpjm1                               ! u_ice v_ice at I-point ; ssu_m, ssv_m at U- & V-points
+               DO ji = 1, jpim1   ! NO vector opt.         !
+                  zu_io(ji,jj) = 0.5_wp * ( u_ice(ji+1,jj+1) + u_ice(ji+1,jj  ) ) - ssu_m(ji,jj)
+                  zv_io(ji,jj) = 0.5_wp * ( v_ice(ji+1,jj+1) + v_ice(ji  ,jj+1) ) - ssv_m(ji,jj)
                END DO
             END DO
@@ -194,9 +177,9 @@
 
          ! frictional velocity at T-point
+         zcoef = 0.5_wp * cw
          DO jj = 2, jpjm1
             DO ji = 2, jpim1   ! NO vector opt. because of zu_io
-               ust2s(ji,jj) = 0.5 * cw                                                          &
-                  &         * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
-                  &            + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
+               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
+                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
             END DO
          END DO
@@ -207,6 +190,6 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ust2s(ji,jj) = zcoef * tms(ji,jj) * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
-                  &                                     + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) )
+               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
+                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1)   ) * tms(ji,jj)
             END DO
          END DO
@@ -217,5 +200,5 @@
       !
       IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
-
+      !
    END SUBROUTINE lim_dyn_2
 
@@ -265,4 +248,10 @@
          WRITE(numout,*) '       coefficient for the solution of int. stresses alphaevp = ', alphaevp
       ENDIF
+      !
+      IF( angvg /= 0._wp .AND. .NOT.lk_lim2_vp ) THEN
+         CALL ctl_warn( 'lim_dyn_init_2: turning angle for oceanic stress not properly coded for EVP ',   &
+            &           '(see limsbc_2 module). We force  angvg = 0._wp'  )
+         angvg = 0._wp
+      ENDIF
 
       !  Initialization
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 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 2370)
@@ -4,4 +4,8 @@
    !! LIM 2.0 ice model :   definition of the ice mesh parameters
    !!======================================================================
+   !! History :   -   ! 2001-04 (LIM) original code
+   !!            1.0  ! 2002-08 (C. Ethe, G. Madec) F90, module
+   !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
+   !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
@@ -10,5 +14,4 @@
    !!   lim_msh_2   : definition of the ice mesh
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE phycst
    USE dom_oce
@@ -20,5 +23,4 @@
    PRIVATE
 
-   !! * Accessibility
    PUBLIC lim_msh_2      ! routine called by ice_ini_2.F90
 
@@ -26,7 +28,6 @@
    !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -43,13 +44,7 @@
       !! 
       !! ** Refer.  : Deleersnijder et al. Ocean Modelling 100, 7-10 
-      !!
-      !! ** History :
-      !!         original    : 01-04 (LIM)
-      !!         addition    : 02-08 (C. Ethe, G. Madec)
       !!--------------------------------------------------------------------- 
-      !! * Local variables
       INTEGER :: ji, jj      ! dummy loop indices
       REAL(wp) ::   zusden   ! local scalars
-
 #if defined key_lim2_vp
       REAL(wp) ::   zusden2           ! local scalars
@@ -267,7 +262,5 @@
          END DO
       END DO
-      
-      !--lateral boundary conditions    
-      CALL lbc_lnk( tmu(:,:), 'I', 1. )
+      CALL lbc_lnk( tmu(:,:), 'I', 1. )      !--lateral boundary conditions    
 #else
       ! EVP rheology : ice velocity point are U- & V-points ; ice vorticity
@@ -278,8 +271,8 @@
       WHERE( fmask(:,:,1) == 1.e0 )   tmf(:,:) = 1.e0
 #endif
-      
+      !
       ! unmasked and masked area of T-grid cell
       area(:,:) = e1t(:,:) * e2t(:,:)
-      
+      !
    END SUBROUTINE lim_msh_2
 
Index: /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 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 2370)
@@ -7,13 +7,12 @@
    !!            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
+   !!             -   !  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'                AND                   LIM-2 sea-ice model
    !!   'key_lim2_vp'                                       VP ice rheology
-   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   lim_rhg_2   : computes ice velocities
@@ -36,6 +35,6 @@
    PUBLIC   lim_rhg_2 ! routine called by lim_dyn
 
-   REAL(wp) ::   rzero   = 0.e0   ! constant value: zero
-   REAL(wp) ::   rone    = 1.e0   !            and  one
+   REAL(wp) ::   rzero   = 0._wp   ! constant value: zero
+   REAL(wp) ::   rone    = 1._wp   !            and  one
 
    !! * Substitutions
@@ -44,7 +43,6 @@
    !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -67,9 +65,9 @@
       INTEGER ::   iter, jter          ! temporary integers
       CHARACTER (len=50) ::   charout
-      REAL(wp) ::   ze11  , ze12  , ze22  , ze21               ! temporary scalars
-      REAL(wp) ::   zt11  , zt12  , zt21  , zt22               !    "         "
-      REAL(wp) ::   zvis11, zvis21, zvis12, zvis22             !    "         "
-      REAL(wp) ::   zgphsx, ztagnx, zgsshx, zunw, zur, zusw    !    "         "
-      REAL(wp) ::   zgphsy, ztagny, zgsshy, zvnw, zvr          !    "         "
+      REAL(wp) ::   ze11  , ze12  , ze22  , ze21               ! local scalars
+      REAL(wp) ::   zt11  , zt12  , zt21  , zt22               !   -      -
+      REAL(wp) ::   zvis11, zvis21, zvis12, zvis22             !   -      -
+      REAL(wp) ::   zgphsx, ztagnx, zgsshx, zunw, zur, zusw    !   -      -
+      REAL(wp) ::   zgphsy, ztagny, zgsshy, zvnw, zvr          !   -      -
       REAL(wp) ::   zresm,  za, zac, zmod
       REAL(wp) ::   zmpzas, zstms, zindu, zusdtp, zmassdt, zcorlal
@@ -91,45 +89,37 @@
       REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zzfrld, zztms
       REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zi1, zi2, zmasst, zpresh
-
       !!-------------------------------------------------------------------
-
-!!bug
-!!    u_oce(:,:) = 0.e0
-!!    v_oce(:,:) = 0.e0
-!!    write(*,*) 'rhg min, max u & v', maxval(u_oce), minval(u_oce), maxval(v_oce), minval(v_oce)
-!!bug
       
       !  Store initial velocities
       !  ----------------
-      zztms(:,0    ) = 0.e0       ;    zzfrld(:,0    ) = 0.e0
-      zztms(:,jpj+1) = 0.e0       ;    zzfrld(:,jpj+1) = 0.e0
-      zu0(:,0    ) = 0.e0         ;    zv0(:,0    ) = 0.e0
-      zu0(:,jpj+1) = 0.e0         ;    zv0(:,jpj+1) = 0.e0
-      zztms(:,1:jpj) = tms(:,:)   ;    zzfrld(:,1:jpj) = frld(:,:)
-      zu0(:,1:jpj) = u_ice(:,:)   ;    zv0(:,1:jpj) = v_ice(:,:)
-
-      zu_a(:,:)    = zu0(:,:)     ;   zv_a(:,:) = zv0(:,:)
-      zu_n(:,:)    = zu0(:,:)     ;   zv_n(:,:) = zv0(:,:)
+      zztms(:,0    ) = 0._wp        ;   zzfrld(:,0    ) = 0._wp
+      zztms(:,jpj+1) = 0._wp        ;   zzfrld(:,jpj+1) = 0._wp
+      zu0  (:,0    ) = 0._wp        ;   zv0   (:,0    ) = 0._wp
+      zu0  (:,jpj+1) = 0._wp        ;   zv0   (:,jpj+1) = 0._wp
+      zztms(:,1:jpj) = tms  (:,:)   ;   zzfrld(:,1:jpj) = frld (:,:)
+      zu0  (:,1:jpj) = u_ice(:,:)   ;   zv0   (:,1:jpj) = v_ice(:,:)
+      zu_a (:, :   ) = zu0  (:,:)   ;   zv_a  (:, :   ) = zv0  (:,:)
+      zu_n (:, :   ) = zu0  (:,:)   ;   zv_n  (:, :   ) = zv0  (:,:)
 
 !i
-      zi1   (:,:) = 0.e0
-      zi2   (:,:) = 0.e0
-      zpresh(:,:) = 0.e0
-      zmasst(:,:) = 0.e0
+      zi1   (:,:) = 0._wp
+      zi2   (:,:) = 0._wp
+      zpresh(:,:) = 0._wp
+      zmasst(:,:) = 0._wp
 !i
 !!gm violant
-      zfrld(:,:) =0.e0
-      zcorl(:,:) =0.e0
-      zmass(:,:) =0.e0
-      za1ct(:,:) =0.e0
-      za2ct(:,:) =0.e0
+      zfrld(:,:) =0._wp
+      zcorl(:,:) =0._wp
+      zmass(:,:) =0._wp
+      za1ct(:,:) =0._wp
+      za2ct(:,:) =0._wp
 !!gm end
 
-      zviszeta(:,:) = 0.e0
-      zviseta (:,:) = 0.e0
-
-!i    zviszeta(:,0    ) = 0.e0    ;    zviseta(:,0    ) = 0.e0
-!i    zviszeta(:,jpj  ) = 0.e0    ;    zviseta(:,jpj  ) = 0.e0
-!i    zviszeta(:,jpj+1) = 0.e0    ;    zviseta(:,jpj+1) = 0.e0
+      zviszeta(:,:) = 0._wp
+      zviseta (:,:) = 0._wp
+
+!i    zviszeta(:,0    ) = 0._wp    ;    zviseta(:,0    ) = 0._wp
+!i    zviszeta(:,jpj  ) = 0._wp    ;    zviseta(:,jpj  ) = 0._wp
+!i    zviszeta(:,jpj+1) = 0._wp    ;    zviseta(:,jpj+1) = 0._wp
 
 
@@ -143,11 +133,11 @@
          DO ji = 1 , jpi
             ! only the sinus changes its sign with the hemisphere
-            zsang(ji,jj)  = SIGN( 1.e0, fcor(ji,jj) ) * sangvg   ! only the sinus changes its sign with the hemisphere
+            zsang(ji,jj)  = SIGN( 1._wp, fcor(ji,jj) ) * sangvg   ! only the sinus changes its sign with the hemisphere
             !
             zmasst(ji,jj) = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) )
             zpresh(ji,jj) = tms(ji,jj) *  pstarh * hicm(ji,jj) * EXP( -c_rhg * frld(ji,jj) )
 !!gm  :: stress given at I-point (F-point for the ocean) only compute the ponderation with the ice fraction (1-frld)
-            zi1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )
-            zi2(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )
+            zi1(ji,jj)    = tms(ji,jj) * ( 1._wp - frld(ji,jj) )
+            zi2(ji,jj)    = tms(ji,jj) * ( 1._wp - frld(ji,jj) )
          END DO
       END DO
@@ -163,5 +153,5 @@
             zstms = zztms(ji,jj  ) * wght(ji,jj,2,2) + zztms(ji-1,jj  ) * wght(ji,jj,1,2)   &
                &  + zztms(ji,jj-1) * wght(ji,jj,2,1) + zztms(ji-1,jj-1) * wght(ji,jj,1,1)
-            zusw  = 1.0 / MAX( zstms, epsd )
+            zusw  = 1._wp / MAX( zstms, epsd )
 
             zt11 = zztms(ji  ,jj  ) * zzfrld(ji  ,jj  ) 
@@ -201,7 +191,7 @@
             ! Gradient of the sea surface height
             zgsshx =  (   (ssh_m(ji  ,jj  ) - ssh_m(ji-1,jj  ))/e1u(ji-1,jj  )   &
-               &       +  (ssh_m(ji  ,jj-1) - ssh_m(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5
+               &       +  (ssh_m(ji  ,jj-1) - ssh_m(ji-1,jj-1))/e1u(ji-1,jj-1)   ) * 0.5_wp
             zgsshy =  (   (ssh_m(ji  ,jj  ) - ssh_m(ji  ,jj-1))/e2v(ji  ,jj-1)   &
-               &       +  (ssh_m(ji-1,jj  ) - ssh_m(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5
+               &       +  (ssh_m(ji-1,jj  ) - ssh_m(ji-1,jj-1))/e2v(ji-1,jj-1)   ) * 0.5_wp
 
             ! Computation of the velocity field taking into account the ice-ice interaction.                                 
@@ -219,5 +209,5 @@
          !                                                ! ==================== !        
          zindu = MOD( iter , 2 )
-         zusdtp = ( zindu * 2.0 + ( 1.0 - zindu ) * 1.0 )  * REAL( nbiter ) / rdt_ice
+         zusdtp = ( zindu * 2._wp + ( 1._wp - zindu ) * 1._wp )  * REAL( nbiter ) / rdt_ice
 
          ! Computation of free drift field for free slip boundary conditions.
@@ -241,8 +231,8 @@
                zdgi = zt12 + zt21
                ztrace2 = zdgp * zdgp 
-               zdeter  = zt11 * zt22 - 0.25 * zdgi * zdgi
+               zdeter  = zt11 * zt22 - 0.25_wp * zdgi * zdgi
 
                !  Creep limit depends on the size of the grid.
-               zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2 ),  creepl)
+               zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4._wp * zdeter ) * usecc2 ),  creepl)
 
                !-  Computation of viscosities.
@@ -256,8 +246,8 @@
             DO ji = 2, fs_jpim1   ! NO vector opt.
                !* zc1u , zc2v
-               zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm
-               zvis12 =       zviseta (ji-1,jj-1) + dm
-               zvis21 =       zviseta (ji-1,jj-1)
-               zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
+               zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
+               zvis12 =         zviseta (ji-1,jj-1) + dm
+               zvis21 =         zviseta (ji-1,jj-1)
+               zvis22 =         zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
                zdiag  = zvis22 * ( akappa(ji-1,jj-1,1,1) + akappa(ji-1,jj-1,2,1) )
                zs11_11 =  zvis11 * akappa(ji-1,jj-1,1,1) + zdiag
@@ -266,8 +256,8 @@
                zs22_11 =  zvis11 * akappa(ji-1,jj-1,2,1) + zdiag
 
-               zvis11 = 2.0 * zviseta (ji,jj-1) + dm
-               zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
-               zvis12 =       zviseta (ji,jj-1) + dm
-               zvis21 =       zviseta (ji,jj-1)
+               zvis11 = 2._wp * zviseta (ji,jj-1) + dm
+               zvis22 =         zviszeta(ji,jj-1) - zviseta(ji,jj-1)
+               zvis12 =         zviseta (ji,jj-1) + dm
+               zvis21 =         zviseta (ji,jj-1)
                zdiag = zvis22 * ( -akappa(ji,jj-1,1,1) + akappa(ji,jj-1,2,1) )
                zs11_21 = -zvis11 * akappa(ji,jj-1,1,1) + zdiag
@@ -276,9 +266,9 @@
                zs21_21 = -zvis12 * akappa(ji,jj-1,1,2) + zvis21 * akappa(ji,jj-1,2,2)
 
-               zvis11 = 2.0 * zviseta (ji-1,jj) + dm
-               zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
-               zvis12 =       zviseta (ji-1,jj) + dm
-               zvis21 =       zviseta (ji-1,jj)
-               zdiag = zvis22 * ( akappa(ji-1,jj,1,1) + akappa(ji-1,jj,2,1) )
+               zvis11 = 2._wp * zviseta (ji-1,jj) + dm
+               zvis22 =         zviszeta(ji-1,jj) - zviseta(ji-1,jj)
+               zvis12 =         zviseta (ji-1,jj) + dm
+               zvis21 =         zviseta (ji-1,jj)
+               zdiag  = zvis22 * ( akappa(ji-1,jj,1,1) + akappa(ji-1,jj,2,1) )
                zs11_12 =  zvis11 * akappa(ji-1,jj,1,1) + zdiag
                zs12_12 = -zvis12 * akappa(ji-1,jj,2,2) - zvis21 * akappa(ji-1,jj,1,2)
@@ -286,8 +276,8 @@
                zs21_12 = -zvis12 * akappa(ji-1,jj,1,2) - zvis21 * akappa(ji-1,jj,2,2)
 
-               zvis11 = 2.0 * zviseta (ji,jj) + dm
-               zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj)
-               zvis12 =       zviseta (ji,jj) + dm
-               zvis21 =       zviseta (ji,jj)
+               zvis11 = 2._wp * zviseta (ji,jj) + dm
+               zvis22 =         zviszeta(ji,jj) - zviseta(ji,jj)
+               zvis12 =         zviseta (ji,jj) + dm
+               zvis21 =         zviseta (ji,jj)
                zdiag = zvis22 * ( -akappa(ji,jj,1,1) + akappa(ji,jj,2,1) )
                zs11_22 = -zvis11 * akappa(ji,jj,1,1) + zdiag
@@ -315,8 +305,8 @@
 
                !* zc1v , zc2v.
-               zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm
-               zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
-               zvis12 =       zviseta (ji-1,jj-1) + dm
-               zvis21 =       zviseta (ji-1,jj-1)
+               zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
+               zvis22 =         zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
+               zvis12 =         zviseta (ji-1,jj-1) + dm
+               zvis21 =         zviseta (ji-1,jj-1)
                zdiag = zvis22 * ( akappa(ji-1,jj-1,1,2) + akappa(ji-1,jj-1,2,2) )
                zs11_11 =  zvis11 * akappa(ji-1,jj-1,1,2) + zdiag
@@ -325,8 +315,8 @@
                zs21_11 =  zvis12 * akappa(ji-1,jj-1,1,1) - zvis21 * akappa(ji-1,jj-1,2,1)
  
-               zvis11 = 2.0 * zviseta (ji,jj-1) + dm
-               zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
-               zvis12 =       zviseta (ji,jj-1) + dm
-               zvis21 =       zviseta (ji,jj-1)
+               zvis11 = 2._wp * zviseta (ji,jj-1) + dm
+               zvis22 =         zviszeta(ji,jj-1) - zviseta(ji,jj-1)
+               zvis12 =         zviseta (ji,jj-1) + dm
+               zvis21 =         zviseta (ji,jj-1)
                zdiag = zvis22 * ( akappa(ji,jj-1,1,2) + akappa(ji,jj-1,2,2) )
                zs11_21 =  zvis11 * akappa(ji,jj-1,1,2) + zdiag
@@ -335,8 +325,8 @@
                zs21_21 = -zvis12 * akappa(ji,jj-1,1,1) - zvis21 * akappa(ji,jj-1,2,1)
 
-               zvis11 = 2.0 * zviseta (ji-1,jj) + dm
-               zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
-               zvis12 =       zviseta (ji-1,jj) + dm
-               zvis21 =       zviseta (ji-1,jj)
+               zvis11 = 2._wp * zviseta (ji-1,jj) + dm
+               zvis22 =         zviszeta(ji-1,jj) - zviseta(ji-1,jj)
+               zvis12 =         zviseta (ji-1,jj) + dm
+               zvis21 =         zviseta (ji-1,jj)
                zdiag = zvis22 * ( akappa(ji-1,jj,1,2) - akappa(ji-1,jj,2,2) )
                zs11_12 =  zvis11 * akappa(ji-1,jj,1,2) + zdiag
@@ -345,8 +335,8 @@
                zs21_12 =  zvis12 * akappa(ji-1,jj,1,1) - zvis21 * akappa(ji-1,jj,2,1)
 
-               zvis11 = 2.0 * zviseta (ji,jj) + dm
-               zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj)
-               zvis12 =       zviseta (ji,jj) + dm
-               zvis21 =       zviseta (ji,jj)
+               zvis11 = 2._wp * zviseta (ji,jj) + dm
+               zvis22 =         zviszeta(ji,jj) - zviseta(ji,jj)
+               zvis12 =         zviseta (ji,jj) + dm
+               zvis21 =         zviseta (ji,jj)
                zdiag = zvis22 * ( akappa(ji,jj,1,2) - akappa(ji,jj,2,2) )
                zs11_22 =  zvis11 * akappa(ji,jj,1,2) + zdiag
@@ -388,8 +378,8 @@
                   ze22 = + akappa(ji,jj-1,2,2) * zv_a(ji+1,jj) + akappa(ji,jj-1,2,1) * zu_a(ji+1,jj)
                   ze21 =   akappa(ji,jj-1,1,1) * zv_a(ji+1,jj) - akappa(ji,jj-1,1,2) * zu_a(ji+1,jj)
-                  zvis11 = 2.0 * zviseta (ji,jj-1) + dm
-                  zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
-                  zvis12 =       zviseta (ji,jj-1) + dm
-                  zvis21 =       zviseta (ji,jj-1)
+                  zvis11 = 2._wp * zviseta (ji,jj-1) + dm
+                  zvis22 =         zviszeta(ji,jj-1) - zviseta(ji,jj-1)
+                  zvis12 =         zviseta (ji,jj-1) + dm
+                  zvis21 =         zviseta (ji,jj-1)
                   zdiag = zvis22 * ( ze11 + ze22 )
                   zs11_21 =  zvis11 * ze11 + zdiag
@@ -406,8 +396,8 @@
                   ze21 =   akappa(ji-1,jj,1,1) * ( zv_a(ji  ,jj+1) - zv_a(ji-1,jj+1) )   &
                      &   - akappa(ji-1,jj,1,2) * ( zu_a(ji  ,jj+1) + zu_a(ji-1,jj+1) )
-                  zvis11 = 2.0 * zviseta (ji-1,jj) + dm
-                  zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
-                  zvis12 =       zviseta (ji-1,jj) + dm
-                  zvis21 =       zviseta (ji-1,jj)
+                  zvis11 = 2._wp * zviseta (ji-1,jj) + dm
+                  zvis22 =         zviszeta(ji-1,jj) - zviseta(ji-1,jj)
+                  zvis12 =         zviseta (ji-1,jj) + dm
+                  zvis21 =         zviseta (ji-1,jj)
                   zdiag = zvis22 * ( ze11 + ze22 )
                   zs11_12 =  zvis11 * ze11 + zdiag
@@ -424,8 +414,8 @@
                   ze21 =   akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji  ,jj+1) )   &
                      &   - akappa(ji,jj,1,2) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) + zu_a(ji  ,jj+1) )
-                  zvis11 = 2.0 * zviseta (ji,jj) + dm
-                  zvis22 =       zviszeta(ji,jj) - zviseta(ji,jj)
-                  zvis12 =       zviseta (ji,jj) + dm
-                  zvis21 =       zviseta (ji,jj)
+                  zvis11 = 2._wp * zviseta (ji,jj) + dm
+                  zvis22 =         zviszeta(ji,jj) - zviseta(ji,jj)
+                  zvis12 =         zviseta (ji,jj) + dm
+                  zvis21 =         zviseta (ji,jj)
                   zdiag = zvis22 * ( ze11 + ze22 )
                   zs11_22 =  zvis11 * ze11 + zdiag
@@ -443,8 +433,8 @@
                   ze21 =   akappa(ji-1,jj-1,1,1) * ( zv_a(ji  ,jj-1) - zv_a(ji-1,jj-1) - zv_a(ji-1,jj) )   &
                      &  -  akappa(ji-1,jj-1,1,2) * ( zu_a(ji  ,jj-1) + zu_a(ji-1,jj-1) + zu_a(ji-1,jj) )
-                  zvis11 = 2.0 * zviseta (ji-1,jj-1) + dm
-                  zvis22 =       zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
-                  zvis12 =       zviseta (ji-1,jj-1) + dm
-                  zvis21 =       zviseta (ji-1,jj-1)
+                  zvis11 = 2._wp * zviseta (ji-1,jj-1) + dm
+                  zvis22 =         zviszeta(ji-1,jj-1) - zviseta(ji-1,jj-1)
+                  zvis12 =         zviseta (ji-1,jj-1) + dm
+                  zvis21 =         zviseta (ji-1,jj-1)
                   zdiag = zvis22 * ( ze11 + ze22 )
                   zs11_11 =  zvis11 * ze11 + zdiag
@@ -461,8 +451,8 @@
                   ze21 =   akappa(ji,jj-1,1,1) * ( zv_a(ji+1,jj-1) - zv_a(ji  ,jj-1) )   &
                      &   - akappa(ji,jj-1,1,2) * ( zu_a(ji+1,jj-1) + zu_a(ji  ,jj-1) )
-                  zvis11 = 2.0 * zviseta (ji,jj-1) + dm
-                  zvis22 =       zviszeta(ji,jj-1) - zviseta(ji,jj-1)
-                  zvis12 =       zviseta (ji,jj-1) + dm
-                  zvis21 =       zviseta (ji,jj-1)
+                  zvis11 = 2._wp * zviseta (ji,jj-1) + dm
+                  zvis22 =         zviszeta(ji,jj-1) - zviseta(ji,jj-1)
+                  zvis12 =         zviseta (ji,jj-1) + dm
+                  zvis21 =         zviseta (ji,jj-1)
                   zdiag = zvis22 * ( ze11 + ze22 )
                   zs11_21 =  zs11_21 + zvis11 * ze11 + zdiag
@@ -475,8 +465,8 @@
                   ze22 = - akappa(ji-1,jj,2,2) * zv_a(ji-1,jj) + akappa(ji-1,jj,2,1) * zu_a(ji-1,jj)
                   ze21 = - akappa(ji-1,jj,1,1) * zv_a(ji-1,jj) - akappa(ji-1,jj,1,2) * zu_a(ji-1,jj)
-                  zvis11 = 2.0 * zviseta (ji-1,jj) + dm
-                  zvis22 =       zviszeta(ji-1,jj) - zviseta(ji-1,jj)
-                  zvis12 =       zviseta (ji-1,jj) + dm
-                  zvis21 =       zviseta (ji-1,jj)
+                  zvis11 = 2._wp * zviseta (ji-1,jj) + dm
+                  zvis22 =         zviszeta(ji-1,jj) - zviseta(ji-1,jj)
+                  zvis12 =         zviseta (ji-1,jj) + dm
+                  zvis21 =         zviseta (ji-1,jj)
                   zdiag = zvis22 * ( ze11 + ze22 )
                   zs11_12 =  zs11_12 + zvis11 * ze11 + zdiag
@@ -506,12 +496,12 @@
                   zvr     = zv_a(ji,jj) - v_oce(ji,jj)
 !!!!
-                  zmod    = SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) )
+                  zmod    = SQRT( zur*zur + zvr*zvr ) * ( 1._wp - zfrld(ji,jj) )
                   za      = rhoco * zmod
 !!!!
-!!gm chg resul    za      = rhoco * SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) )
+!!gm chg resul    za      = rhoco * SQRT( zur*zur + zvr*zvr ) * ( 1._wp - zfrld(ji,jj) )
                   zac     = za * cangvg
                   zmpzas  = alpha * zcorl(ji,jj) + za * zsang(ji,jj)
                   zmassdt = zusdtp * zmass(ji,jj)
-                  zcorlal = ( 1.0 - alpha ) * zcorl(ji,jj)
+                  zcorlal = ( 1._wp - alpha ) * zcorl(ji,jj)
 
                   za1 =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   &
@@ -527,5 +517,5 @@
                   zunw   = (  ( za1 + zd1 ) * zc2 + ( za2 + zd2 ) * zc1 ) * zden
                   zvnw   = (  ( za2 + zd2 ) * zb1 - ( za1 + zd1 ) * zb2 ) * zden
-                  zmask  = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj)
+                  zmask  = ( 1._wp - MAX( rzero, SIGN( rone , 1._wp - zmass(ji,jj) ) ) ) * tmu(ji,jj)
 
                   zu_n(ji,jj) = ( zu_a(ji,jj) + om * ( zunw - zu_a(ji,jj) ) * tmu(ji,jj) ) * zmask
@@ -582,5 +572,5 @@
 #else
    !!----------------------------------------------------------------------
-   !!   Default option          Dummy module       NO 2.0 LIM sea-ice model
+   !!   Default option        Dummy module      NO VP & LIM-2 sea-ice model
    !!----------------------------------------------------------------------
 CONTAINS
Index: /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 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limrst_2.F90	(revision 2370)
@@ -121,7 +121,7 @@
       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(:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'stress1_i'  , stress1_i (:,:) )    ! EVP rheology
+      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i'  , stress2_i (:,:) )
+      CALL iom_rstput( iter, nitrst, numriw, 'stress12_i' , stress12_i(:,:) )
 #endif
       CALL iom_rstput( iter, nitrst, numriw, 'sxice'      , sxice (:,:)   )
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 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 2370)
@@ -2,10 +2,11 @@
    !!======================================================================
    !!                       ***  MODULE limsbc_2   ***
-   !!           computation of the flux at the sea ice/ocean interface
+   !! LIM-2 :   updates the fluxes at the ocean surface with ice-ocean fluxes
    !!======================================================================
    !! History :  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
+   !!            3.3  ! 2009-05 (G. Garric, C. Bricaud) addition of the lim2_evp case
+   !!             -   ! 2010-11 (G. Madec) ice-ocean stress computed at each ocean time-step
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -13,18 +14,19 @@
    !!   'key_lim2'                                    LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!   lim_sbc_2  : flux at the ice / ocean interface
+   !!   lim_sbc_flx_2  : update mass, heat and salt fluxes at the ocean surface
+   !!   lim_sbc_tau_2  : update i- and j-stresses, and its modulus at the ocean surface
    !!----------------------------------------------------------------------
    USE par_oce          ! ocean parameters
+   USE phycst           ! physical constants
    USE dom_oce          ! ocean domain
+   USE dom_ice_2        ! LIM-2: ice domain
+   USE ice_2            ! LIM-2: ice variables
    USE sbc_ice          ! surface boundary condition: ice
    USE sbc_oce          ! surface boundary condition: ocean
-   USE phycst           ! physical constants
-   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              ! I/O library
    USE albedo           ! albedo parameters
    USE prtctl           ! Print control
@@ -34,12 +36,16 @@
    PRIVATE
 
-   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                !     -      -
+   PUBLIC   lim_sbc_flx_2   ! called by sbc_ice_lim_2
+   PUBLIC   lim_sbc_tau_2   ! called by sbc_ice_lim_2
+
+   REAL(wp)  ::   r1_rdtice            ! = 1. / rdt_ice 
+   REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values
+   REAL(wp)  ::   rzero  = 0._wp       !     -      -
+   REAL(wp)  ::   rone   = 1._wp       !     -      -
    !
-   REAL(wp), DIMENSION(jpi,jpj) ::   soce_r, sice_r   ! constant SSS and ice salinity used in levitating sea-ice case
+   REAL(wp), DIMENSION(jpi,jpj) ::   soce_0, sice_0   ! constant SSS and ice salinity used in levitating sea-ice case
+
+   REAL(wp), DIMENSION(jpi,jpj) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress              [N/m2]
+   REAL(wp), DIMENSION(jpi,jpj) ::   tmod_io              ! modulus of the ice-ocean relative velocity   [m/s]
 
    !! * Substitutions
@@ -48,10 +54,9 @@
    !! NEMO/LIM2 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
-   SUBROUTINE lim_sbc_2( kt )
+   SUBROUTINE lim_sbc_flx_2( kt )
       !!-------------------------------------------------------------------
       !!                ***  ROUTINE lim_sbc_2 ***
@@ -77,48 +82,35 @@
       !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
       !!---------------------------------------------------------------------
-      INTEGER ::   kt    ! number of iteration
+      INTEGER, INTENT(in) ::   kt    ! number of iteration
       !!
-      INTEGER  ::   ji, jj           ! dummy loop indices
+      INTEGER  ::   ji, jj   ! dummy loop indices
       INTEGER  ::   ii0, ii1, ij0, ij1         ! local integers
       INTEGER  ::   ifvt, i1mfr, idfr, iflt    !   -       -
       INTEGER  ::   ial, iadv, ifral, ifrdv    !   -       -
-      REAL(wp) ::   zqsr, zqns, 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) ::   zqsr, zqns, zfm            ! local scalars
+      REAL(wp) ::   zinda, zfons, zemp         !   -      -
+      REAL(wp), DIMENSION(jpi,jpj)   ::   zqnsoce       ! 2D workspace
       REAL(wp), DIMENSION(jpi,jpj,1) ::   zalb, zalbp   ! 2D/3D workspace
-
       !!---------------------------------------------------------------------
      
-      
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-#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,*) '~~~~~~~~~   '
+         IF(lwp) WRITE(numout,*) 'lim_sbc_flx_2 : LIM-2 sea-ice - surface boundary condition - Mass, heat & salt fluxes'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   '
          !
          r1_rdtice = 1. / rdt_ice
          !
-         soce_r(:,:) = soce
-         sice_r(:,:) = sice
-         !
-         IF( cp_cfg == "orca" ) THEN
-           !   ocean/ice salinity in the Baltic sea 
-           DO jj = 1, jpj
-              DO ji = 1, jpi
-                 IF( glamt(ji,jj) >= 14. .AND.  glamt(ji,jj) <= 32. .AND. gphit(ji,jj) >= 54. .AND. gphit(ji,jj) <= 66. ) THEN 
-                   soce_r(ji,jj) = 4.e0 
-                   sice_r(ji,jj) = 2.e0
-                 END IF
-              END DO
-           END DO
-           !
-         END IF
-      END IF
+         soce_0(:,:) = soce                     ! constant SSS and ice salinity used in levitating sea-ice case
+         sice_0(:,:) = sice
+         !
+         IF( cp_cfg == "orca" ) THEN            ! decrease ocean & ice reference salinities in the Baltic sea 
+            WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &
+               &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
+               soce_0(:,:) = 4._wp
+               sice_0(:,:) = 2._wp
+            END WHERE
+         ENDIF
+         !
+      ENDIF
 
       !------------------------------------------!
@@ -126,8 +118,4 @@
       !------------------------------------------!
 
-!!gm
-!!gm CAUTION   
-!!gm re-verifies the non solar expression, especially over open ocen
-!!gm
       zqnsoce(:,:) = qns(:,:)
       DO jj = 1, jpj
@@ -190,5 +178,5 @@
 
             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
@@ -203,12 +191,7 @@
       !      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
             ! freshwater exchanges at the ice-atmosphere / ocean interface (coupled mode)
@@ -223,80 +206,21 @@
                !                                                   !  ice-covered fraction:
 #endif            
-
+            !
             !  computing salt exchanges at the ice/ocean interface
-            zfons =  ( soce_r(ji,jj) - sice_r(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice ) 
-            
+            zfons =  ( soce_0(ji,jj) - sice_0(ji,jj) ) * ( rdmicif(ji,jj) * r1_rdtice ) 
+            !
             !  converting the salt flux from ice to a freshwater flux from ocean
             zfm  = zfons / ( sss_m(ji,jj) + epsi16 )
-            
+            !
             emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution)
             emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution)
-
+            !
          END DO
       END DO
 
-      IF( lk_diaar5 ) THEN
+      IF( lk_diaar5 ) THEN       ! AR5 diagnostics
          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
-
-      !------------------------------------------!
-      !    momentum flux at the ocean surface    !
-      !------------------------------------------!
-
-      IF ( ln_limdyn ) THEN                        ! Update the stress over ice-over area (only in ice-dynamic case)
-         !                                         ! otherwise the atmosphere-ocean stress is used everywhere
-
-         ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point)
-!CDIR NOVERRCHK
-         DO jj = 1, jpj
-!CDIR NOVERRCHK
-            DO ji = 1, jpi
-               ! ... change the cosinus angle sign in the south hemisphere
-               zsang  = SIGN(1.e0, gphif(ji,jj) ) * sangvg
-               ! ... ice velocity relative to the ocean at I-point
-               zu_io  = u_ice(ji,jj) - u_oce(ji,jj)
-               zv_io  = v_ice(ji,jj) - v_oce(ji,jj)
-               zmod   = SQRT( zu_io * zu_io + zv_io * zv_io )
-               zztmp  = rhoco * zmod
-               ! ... components of ice stress over ocean with a ice-ocean rotation angle (at I-point)
-               ztio_u(ji,jj) = zztmp * ( cangvg * zu_io - zsang * zv_io )
-               ztio_v(ji,jj) = zztmp * ( cangvg * zv_io + zsang * zu_io )
-               ! ... module of ice stress over ocean (at I-point)
-               ztiomi(ji,jj) = zztmp * zmod
-               ! 
-            END DO
-         END DO
-
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1   ! NO vector opt.
-               ! ... components of ice-ocean stress at U and V-points  (from I-point values)
-#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  ) )
-               zfrldv = 0.5 * ( frld(ji,jj) + frld(ji  ,jj+1) )
-               ! ... update components of surface ocean stress (ice-cover wheighted)
-               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * zutau
-               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * zvtau
-               ! ... module of ice-ocean stress at T-points (from I-point values)
-               zztmp = 0.25 * ( ztiomi(ji,jj) + ztiomi(ji+1,jj) + ztiomi(ji,jj+1) + ztiomi(ji+1,jj+1) )
-               ! ... update module of surface ocean stress (ice-cover wheighted)
-               taum(ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1. - frld(ji,jj) ) * zztmp
-               !
-            END DO
-         END DO
-
-         ! boundary condition on the stress (utau,vtau,taum)
-         CALL lbc_lnk( utau, 'U', -1. )
-         CALL lbc_lnk( vtau, 'V', -1. )
-         CALL lbc_lnk( taum, 'T',  1. )
-
+         CALL iom_put( 'fsal_virt_cea',   soce_0(:,:) * rdmicif(:,:) * r1_rdtice )
+         CALL iom_put( 'fsal_real_cea', - sice_0(:,:) * rdmicif(:,:) * r1_rdtice )
       ENDIF
 
@@ -305,5 +229,5 @@
       !-----------------------------------------------!
 
-      IF ( lk_cpl ) THEN           
+      IF( lk_cpl ) THEN          ! coupled case
          ! Ice surface temperature 
          tn_ice(:,:,1) = sist(:,:)          ! sea-ice surface temperature       
@@ -314,5 +238,5 @@
       ENDIF
 
-      IF(ln_ctl) THEN
+      IF(ln_ctl) THEN            ! control print
          CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
          CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ')
@@ -321,14 +245,171 @@
          CALL prt_ctl(tab2d_1=fr_i  , clinfo1=' lim_sbc: fr_i   : ', tab2d_2=tn_ice(:,:,1), clinfo2=' tn_ice  : ')
       ENDIF 
-   
-    END SUBROUTINE lim_sbc_2
+      !
+   END SUBROUTINE lim_sbc_flx_2
+
+
+   SUBROUTINE lim_sbc_tau_2( kt , pu_oce, pv_oce )
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE lim_sbc_tau ***
+      !!  
+      !! ** Purpose : Update the ocean surface stresses due to the ice
+      !!         
+      !! ** Action  : * at each ice time step (every nn_fsbc time step):
+      !!                - compute the modulus of ice-ocean relative velocity 
+      !!                  at T-point (C-grid) or I-point (B-grid)
+      !!                      tmod_io = rhoco * | U_ice-U_oce |
+      !!                - update the modulus of stress at ocean surface
+      !!                      taum = frld * taum + (1-frld) * tmod_io * | U_ice-U_oce |
+      !!              * at each ocean time step (each kt): 
+      !!                  compute linearized ice-ocean stresses as
+      !!                      Utau = tmod_io * | U_ice - pU_oce |
+      !!                using instantaneous current ocean velocity (usually before)
+      !!
+      !!    NB: - the averaging operator used depends on the ice dynamics grid (cp_ice_msh='I' or 'C')
+      !!        - ice-ocean rotation angle only allowed in cp_ice_msh='I' case
+      !!        - here we make an approximation: taum is only computed every ice time step
+      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids 
+      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
+      !!
+      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
+      !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
+      !!---------------------------------------------------------------------
+      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
+      !!
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      REAL(wp) ::   zfrldu, zat_u, zu_i, zutau_ice, zu_t, zmodt   ! local scalar
+      REAL(wp) ::   zfrldv, zat_v, zv_i, zvtau_ice, zv_t, zmodi   !   -      -
+      REAL(wp) ::   zsang, zumt                                          !    -         -
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice
+      !!---------------------------------------------------------------------
+      !
+      IF( kt == nit000 .AND. lwp ) THEN         ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'lim_sbc_tau_2 : LIM 2.0 sea-ice - surface ocean momentum fluxes'
+         WRITE(numout,*) '~~~~~~~~~~~~~ '
+         IF( lk_lim2_vp )   THEN   ;   WRITE(numout,*) '                VP  rheology - B-grid case'
+         ELSE                      ;   WRITE(numout,*) '                EVP rheology - C-grid case'
+         ENDIF
+      ENDIF
+      !
+      SELECT CASE( cp_ice_msh )     
+      !                             !-----------------------!
+      CASE( 'I' )                   !  B-grid ice dynamics  !   I-point (i.e. F-point with sea-ice indexation)
+         !                          !--=--------------------!
+         !
+         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==! (i.e. surface module time-step)
+!CDIR NOVERRCHK
+            DO jj = 1, jpj                               !* modulus of ice-ocean relative velocity at I-point
+!CDIR NOVERRCHK
+               DO ji = 1, jpi
+                  zu_i  = u_ice(ji,jj) - u_oce(ji,jj)                   ! ice-ocean relative velocity at I-point
+                  zv_i  = v_ice(ji,jj) - v_oce(ji,jj)
+                  tmod_io(ji,jj) = SQRT( zu_i * zu_i + zv_i * zv_i )    ! modulus of this velocity (at I-point)
+               END DO
+            END DO
+!CDIR NOVERRCHK
+            DO jj = 1, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
+!CDIR NOVERRCHK
+               DO ji = 1, jpim1   ! NO vector opt.
+                  !                                               ! modulus of U_ice-U_oce at T-point
+                  zumt  = 0.25_wp * (  tmod_io(ji+1,jj) + tmod_io(ji+1,jj+1)    &   
+                     &               + tmod_io(ji,jj+1) + tmod_io(ji+1,jj+1)  )
+                  !                                               ! update the modulus of stress at ocean surface
+                  taum(ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1._wp - frld(ji,jj) ) * rhoco * zumt * zumt
+               END DO
+            END DO
+            CALL lbc_lnk( taum, 'T', 1. )
+            !
+            utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
+            vtau_oce(:,:) = vtau(:,:)
+            !
+         ENDIF
+         !
+         !                                        !==  at each ocean time-step  ==!
+         !
+         !                                               !* ice/ocean stress WITH a ice-ocean rotation angle at I-point
+         DO jj = 2, jpj
+            zsang  = SIGN( 1._wp, gphif(1,jj) ) * sangvg          ! change the cosine angle sign in the SH 
+            DO ji = 2, jpi    ! NO vect. opt. possible
+               ! ... ice-ocean relative velocity at I-point using instantaneous surface ocean current at u- & v-pts
+               zu_i = u_ice(ji,jj) - 0.5_wp * ( pu_oce(ji-1,jj  ) + pu_oce(ji-1,jj-1) ) * tmu(ji,jj)
+               zv_i = v_ice(ji,jj) - 0.5_wp * ( pv_oce(ji  ,jj-1) + pv_oce(ji-1,jj-1) ) * tmu(ji,jj)
+               ! ... components of stress with a ice-ocean rotation angle 
+               zmodi = rhoco * tmod_io(ji,jj)                     
+               ztio_u(ji,jj) = zmodi * ( cangvg * zu_i - zsang * zv_i )
+               ztio_v(ji,jj) = zmodi * ( cangvg * zv_i + zsang * zu_i )
+            END DO
+         END DO
+         !                                               !* surface ocean stresses at u- and v-points
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1   ! NO vector opt.
+               !                                   ! ice-ocean stress at U and V-points  (from I-point values)
+               zutau_ice  = 0.5_wp * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
+               zvtau_ice  = 0.5_wp * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
+               !                                   ! open-ocean (lead) fraction at U- & V-points (from T-point values)
+               zfrldu = 0.5_wp * ( frld(ji,jj) + frld(ji+1,jj  ) )
+               zfrldv = 0.5_wp * ( frld(ji,jj) + frld(ji  ,jj+1) )
+               !                                   ! update the surface ocean stress (ice-cover wheighted)
+               utau(ji,jj) = zfrldu * utau_oce(ji,jj) + ( 1._wp - zfrldu ) * zutau_ice
+               vtau(ji,jj) = zfrldv * vtau_oce(ji,jj) + ( 1._wp - zfrldv ) * zvtau_ice
+            END DO
+         END DO
+         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )     ! lateral boundary condition
+         !
+         !
+         !                          !-----------------------!
+      CASE( 'C' )                   !  C-grid ice dynamics  !   U & V-points (same as in the ocean)
+         !                          !--=--------------------!
+         !
+         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==! (i.e. surface module time-step)
+!CDIR NOVERRCHK
+            DO jj = 2, jpjm1                          !* modulus of the ice-ocean velocity at T-point
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1
+                  zu_t  = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   ! 2*(U_ice-U_oce) at T-point
+                  zv_t  = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1)      
+                  zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )                      ! |U_ice-U_oce|^2
+                  !                                               ! update the modulus of stress at ocean surface
+                  taum   (ji,jj) = frld(ji,jj) * taum(ji,jj) + ( 1._wp - frld(ji,jj) ) * rhoco * zmodt
+                  tmod_io(ji,jj) = SQRT( zmodt ) * rhoco          ! rhoco*|Uice-Uoce|
+               END DO
+            END DO
+            CALL lbc_lnk( taum, 'T', 1. )   ;   CALL lbc_lnk( tmod_io, 'T', 1. )
+            !
+            utau_oce(:,:) = utau(:,:)                 !* save the air-ocean stresses at ice time-step
+            vtau_oce(:,:) = vtau(:,:)
+            !
+         ENDIF
+         !
+         !                                        !==  at each ocean time-step  ==!
+         !
+         DO jj = 2, jpjm1                             !* ice stress over ocean WITHOUT a ice-ocean rotation angle
+            DO ji = fs_2, fs_jpim1
+               !                                            ! ocean area at u- & v-points
+               zfrldu  = 0.5_wp * ( frld(ji,jj) + frld(ji+1,jj  ) )
+               zfrldv  = 0.5_wp * ( frld(ji,jj) + frld(ji  ,jj+1) )
+               !                                            ! quadratic drag formulation without rotation
+               !                                            ! using instantaneous surface ocean current
+               zutau_ice = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
+               zvtau_ice = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
+               !                                            ! update the surface ocean stress (ice-cover wheighted)
+               utau(ji,jj) = zfrldu * utau_oce(ji,jj) + ( 1._wp - zfrldu ) * zutau_ice
+               vtau(ji,jj) = zfrldv * vtau_oce(ji,jj) + ( 1._wp - zfrldv ) * zvtau_ice
+            END DO
+         END DO
+         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
+         !
+      END SELECT
+
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
+         &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask )
+      !  
+   END SUBROUTINE lim_sbc_tau_2
 
 #else
    !!----------------------------------------------------------------------
-   !!   Default option :        Dummy module       NO LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_sbc_2         ! Dummy routine
-   END SUBROUTINE lim_sbc_2
+   !!   Default option         Empty module        NO LIM 2.0 sea-ice model
+   !!----------------------------------------------------------------------
 #endif 
 
Index: /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 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 2370)
@@ -7,5 +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
+   !!            3.3  !  2009-05  (G. Garric, C. Bricaud) addition of the lim2_evp case
    !!----------------------------------------------------------------------
 #if defined key_lim2
@@ -41,5 +41,4 @@
    REAL(wp)  ::   rone   = 1.e0
 
-
    !! * Substitution
 #  include "vectopt_loop_substitute.h90"
@@ -88,21 +87,22 @@
          ! ice velocities at ocean U- and V-points (zui_u,zvi_v)
          ! ---------------------------------------
-         zvbord = 1.0 + ( 1.0 - bound )      ! zvbord=2 no-slip, =0 free slip boundary conditions        
-#if defined key_lim2_vp
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1   ! NO vector opt.
-               zui_u(ji,jj) = ( u_ice(ji+1,jj  ) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj  ) + tmu(ji+1,jj+1), zvbord ) )
-               zvi_v(ji,jj) = ( v_ice(ji  ,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji  ,jj+1) + tmu(ji+1,jj+1), zvbord ) )
-            END DO
-         END DO
-         CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )         ! Lateral boundary conditions
-#else
-        zui_u(:,:) = u_ice(:,:)      ! EVP rheology: ice (u,v) at u- and v-points
-        zvi_v(:,:) = v_ice(:,:)
-#endif
+         IF( lk_lim2_vp ) THEN      ! VP rheology : B-grid sea-ice dynamics (I-point ice velocity)
+            zvbord = 1._wp + ( 1._wp - bound )      ! zvbord=2 no-slip, =0 free slip boundary conditions        
+            DO jj = 1, jpjm1
+               DO ji = 1, jpim1   ! NO vector opt.
+                  zui_u(ji,jj) = ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj)+tmu(ji+1,jj+1), zvbord ) )
+                  zvi_v(ji,jj) = ( v_ice(ji,jj+1) + v_ice(ji+1,jj+1) ) / ( MAX( tmu(ji,jj+1)+tmu(ji+1,jj+1), zvbord ) )
+               END DO
+            END DO
+            CALL lbc_lnk( zui_u, 'U', -1. )   ;   CALL lbc_lnk( zvi_v, 'V', -1. )      ! Lateral boundary conditions
+            !
+         ELSE                       ! EVP rheology : C-grid sea-ice dynamics (u- & v-points ice velocity)
+            zui_u(:,:) = u_ice(:,:)      ! EVP rheology: ice (u,v) at u- and v-points
+            zvi_v(:,:) = v_ice(:,:)
+         ENDIF
 
          ! CFL test for stability
          ! ----------------------
-         zcfl  = 0.e0
+         zcfl  = 0._wp
          zcfl  = MAX( zcfl, MAXVAL( ABS( zui_u(1:jpim1, :     ) ) * rdt_ice / e1u(1:jpim1, :     ) ) )
          zcfl  = MAX( zcfl, MAXVAL( ABS( zvi_v( :     ,1:jpjm1) ) * rdt_ice / e2v( :     ,1:jpjm1) ) )
@@ -193,6 +193,6 @@
 !        DO jj = 1, jpjm1          ! NB: has not to be defined on jpj line and jpi row
 !           DO ji = 1 , fs_jpim1   ! vector opt.
-!              IF( MIN( zs0a(ji,jj) , zs0a(ji+1,jj) ) == 0.e0 )   pahu(ji,jj) = 0.e0
-!              IF( MIN( zs0a(ji,jj) , zs0a(ji,jj+1) ) == 0.e0 )   pahv(ji,jj) = 0.e0
+!              IF( MIN( zs0a(ji,jj) , zs0a(ji+1,jj) ) == 0.e0 )   pahu(ji,jj) = 0._wp
+!              IF( MIN( zs0a(ji,jj) , zs0a(ji,jj+1) ) == 0.e0 )   pahv(ji,jj) = 0._wp
 !           END DO
 !        END DO
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 2370)
@@ -4,4 +4,7 @@
    !!   Sea-Ice dynamics :  
    !!======================================================================
+   !! history :  1.0  ! 2002-08 (C. Ethe, G. Madec)  original VP code 
+   !!            3.0  ! 2007-03 (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)  LIM3: EVP-Cgrid
+   !!----------------------------------------------------------------------
 #if defined key_lim3
    !!----------------------------------------------------------------------
@@ -11,5 +14,4 @@
    !!    lim_dyn_init : initialization and namelist read
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE phycst
    USE in_out_manager  ! I/O manager
@@ -30,19 +32,13 @@
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC lim_dyn  ! routine called by ice_step
+   PUBLIC   lim_dyn   ! routine called by ice_step
 
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
-
-   !! * Module variables
-   REAL(wp)  ::  rone    = 1.e0   ! constant value
-
    !!----------------------------------------------------------------------
    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -59,15 +55,10 @@
       !!              - computation of the stress at the ocean surface         
       !!              - treatment of the case if no ice dynamic
-      !! History :
-      !!   1.0  !  01-04   (LIM)  Original code
-      !!   2.0  !  02-08   (C. Ethe, G. Madec)  F90, mpp
-      !!   3.0  !  2007-03 (M.A. Morales Maqueda, S. Bouillon, M. Vancoppenolle) 
-      !!                   LIM3, EVP, C-grid
       !!------------------------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! number of iteration
-      !! * Local variables
+      !!
       INTEGER  ::   ji, jj, jl, ja    ! dummy loop indices
       INTEGER  ::   i_j1, i_jpj       ! Starting/ending j-indices for rheology
-      REAL(wp) ::   zcoef             ! temporary scalar
+      REAL(wp) ::   zcoef             ! local scalar
       REAL(wp), DIMENSION(jpj)     ::   zind           ! i-averaged indicator of sea-ice
       REAL(wp), DIMENSION(jpj)     ::   zmsk           ! i-averaged of tmask
@@ -158,9 +149,9 @@
          zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
          ! frictional velocity at T-point
+         zcoef = 0.5 * cw
          DO jj = 2, jpjm1 
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ust2s(ji,jj) = 0.5 * cw                                                          &
-                  &         * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
-                  &            + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
+               ust2s(ji,jj) = zcoef * (  zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj)   &
+                  &                    + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1)   ) * tms(ji,jj)
             END DO
          END DO
@@ -171,11 +162,10 @@
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ust2s(ji,jj) = zcoef * tms(ji,jj) * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
-                  &                                     + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) )
+               ust2s(ji,jj) = zcoef * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
+                  &                        + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1)   ) * tms(ji,jj)
             END DO
          END DO
          !
       ENDIF
-
       CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
 
@@ -218,6 +208,7 @@
          END DO
       ENDIF
-
+      !
    END SUBROUTINE lim_dyn
+
 
    SUBROUTINE lim_dyn_init
@@ -232,9 +223,4 @@
       !!
       !! ** input   :   Namelist namicedyn
-      !!
-      !! history :
-      !!  8.5  ! 03-08 (C. Ethe) original code
-      !!  9.0  ! 07-03 (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)
-      !!               EVP-Cgrid-LIM3
       !!-------------------------------------------------------------------
       NAMELIST/namicedyn/ epsd, alpha,     &
@@ -244,11 +230,8 @@
       !!-------------------------------------------------------------------
 
-      ! Define the initial parameters
-      ! -------------------------
-
-      ! Read Namelist namicedyn
-      REWIND ( numnam_ice )
-      READ   ( numnam_ice  , namicedyn )
-      IF(lwp) THEN
+      REWIND( numnam_ice )                ! Read Namelist namicedyn
+      READ  ( numnam_ice  , namicedyn )
+      
+      IF(lwp) THEN                        ! control print
          WRITE(numout,*)
          WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
@@ -272,13 +255,18 @@
          WRITE(numout,*) '   timescale for elastic waves                      telast = ', telast
          WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp
-
-      ENDIF
-
-      usecc2 = 1.0 / ( ecc * ecc )
-      rhoco  = rau0 * cw
+      ENDIF
+      !
+      IF( angvg /= 0._wp ) THEN
+         CALL ctl_warn( 'lim_dyn_init: turning angle for oceanic stress not properly coded for EVP ',   &
+            &           '(see limsbc module). We force  angvg = 0._wp'  )
+         angvg = 0._wp
+      ENDIF
+      
+      usecc2 = 1._wp / ( ecc * ecc )
+      rhoco  = rau0  * cw
       angvg  = angvg * rad
       sangvg = SIN( angvg )
       cangvg = COS( angvg )
-      pstarh = pstar / 2.0
+      pstarh = pstar * 0.5_wp
 
       !  Diffusion coefficients.
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 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 2370)
@@ -11,25 +11,24 @@
 #if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
    !!----------------------------------------------------------------------
-   !!   'key_lim3'                                      LIM3 sea-ice model
+   !!   'key_lim3'               OR                     LIM-3 sea-ice model
+   !!   'key_lim2' AND NOT 'key_lim2_vp'             VP LIM-2 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_rhg   : computes ice velocities
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE phycst
-   USE par_oce
-   USE dom_oce
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE sbc_ice         ! Surface boundary condition: ice fields
-   USE lbclnk
-   USE lib_mpp
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
+   USE phycst           ! Physical constant
+   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 link
+   USE lib_mpp          ! MPP library
+   USE in_out_manager   ! I/O manager
+   USE prtctl           ! Print control
 #if defined key_lim3
-   USE ice
-   USE dom_ice
-   USE iceini
-   USE limitd_me
-#endif
-#if defined key_lim2 && ! defined key_lim2_vp
+   USE ice              ! LIM-3: ice variables
+   USE dom_ice          ! LIM-3: ice domain
+   USE iceini           ! LIM-3: ice initialisation
+   USE limitd_me        ! LIM-3: 
+#else
    USE ice_2            ! LIM2: ice variables
    USE dom_ice_2        ! LIM2: ice domain
@@ -37,28 +36,22 @@
 #endif
 
-
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_rhg  ! routine called by lim_dyn
-
+   PUBLIC   lim_rhg   ! routine called by lim_dyn (or lim_dyn_2)
+
+   REAL(wp) ::   rzero   = 0._wp   ! constant values
+   REAL(wp) ::   rone    = 1._wp   ! constant values
+      
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
-
-   !! * Module variables
-   REAL(wp)  ::           &  ! constant values
-      rzero   = 0.e0   ,  &
-      rone    = 1.e0
    !!----------------------------------------------------------------------
    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
    SUBROUTINE lim_rhg( k_j1, k_jpj )
-
       !!-------------------------------------------------------------------
       !!                 ***  SUBROUTINE lim_rhg  ***
@@ -107,34 +100,16 @@
       !!                 e.g. in the Canadian Archipelago
       !!
-      !! ** References : Hunke and Dukowicz, JPO97
-      !!                 Bouillon et al., 08, in prep (update this when
-      !!                 published)
-      !!                 Vancoppenolle et al., OM08
+      !! References : Hunke and Dukowicz, JPO97
+      !!              Bouillon et al., Ocean Modelling 2009
+      !!              Vancoppenolle et al., Ocean Modelling 2008
+      !!-------------------------------------------------------------------
+      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
-
+      INTEGER ::   ji, jj   ! dummy loop indices
+      INTEGER ::   jter     ! local 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) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         !
+      REAL(wp) ::   za, zstms, zsang, zmask   ! local scalars
 
       REAL(wp),DIMENSION(jpi,jpj) :: &
@@ -154,39 +129,34 @@
 
       REAL(wp) :: &
-         dtevp,                      & !: time step for subcycling
-         dtotel,                     & !:
-         ecc2,                       & !: square of yield ellipse eccenticity
-         z0,                         & !: temporary scalar
-         zr,                         & !: temporary scalar
-         zcca, zccb,                 & !: temporary scalars
-         zu_ice2,                    & !: 
-         zv_ice1,                    & !:
-         zddc, zdtc,                 & !: temporary array for delta on corners
-         zdst,                       & !: temporary array for delta on centre
-         zdsshx, zdsshy,             & !: term for the gradient of ocean surface
-         sigma1, sigma2                !: internal ice stress
+         dtevp,                      & ! time step for subcycling
+         dtotel,                     & !
+         ecc2,                       & ! square of yield ellipse eccenticity
+         z0,                         & ! temporary scalar
+         zr,                         & ! temporary scalar
+         zcca, zccb,                 & ! temporary scalars
+         zu_ice2,                    & ! 
+         zv_ice1,                    & !
+         zddc, zdtc,                 & ! temporary array for delta on corners
+         zdst,                       & ! temporary array for delta on centre
+         zdsshx, zdsshy,             & ! term for the gradient of ocean surface
+         sigma1, sigma2                ! internal ice stress
+
+      REAL(wp),DIMENSION(jpi,jpj) ::   zf1, zf2   ! arrays for internal stresses
 
       REAL(wp),DIMENSION(jpi,jpj) :: &
-         zf1, zf2                      !: arrays for internal stresses
-
-      REAL(wp),DIMENSION(jpi,jpj) :: &
-         zdd, zdt,                   & !: Divergence and tension at centre of grid cells
-         zds,                        & !: Shear on northeast corner of grid cells
-         deltat,                     & !: Delta at centre of grid cells
-         deltac,                     & !: Delta on corners
-         zs1, zs2,                   & !: Diagonal stress tensor components zs1 and zs2 
-         zs12                          !: Non-diagonal stress tensor component zs12
+         zdd, zdt,                   & ! Divergence and tension at centre of grid cells
+         zds,                        & ! Shear on northeast corner of grid cells
+         deltat,                     & ! Delta at centre of grid cells
+         deltac,                     & ! Delta on corners
+         zs1, zs2,                   & ! Diagonal stress tensor components zs1 and zs2 
+         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
+         zresm            ,          & ! Maximal error on ice velocity
+         zindb            ,          & ! ice (1) or not (0)      
+         zdummy                        ! dummy argument
+
+      REAL(wp),DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr   ! Local error on velocity
       !!-------------------------------------------------------------------
-
 #if  defined key_lim2 && ! defined key_lim2_vp
 # if defined key_agrif
@@ -205,26 +175,23 @@
       !
       ! Put every vector to 0
-      zpresh(:,:) = 0.0 ; zc1(:,:) = 0.0
-      zpreshc(:,:) = 0.0
-      u_ice2(:,:)  = 0.0 ; v_ice1(:,:)  = 0.0
-      zdd(:,:)     = 0.0 ; zdt(:,:)     = 0.0 ; zds(:,:)     = 0.0
+      zpresh (:,:) = 0._wp   ;   zc1   (:,:) = 0._wp
+      zpreshc(:,:) = 0._wp
+      u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp
+      zdd    (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp
 
 #if defined key_lim3
-      ! Ice strength on T-points
-      CALL lim_itd_me_icestrength(ridge_scheme_swi)
+      CALL lim_itd_me_icestrength( ridge_scheme_swi )      ! LIM-3: Ice strength on T-points
 #endif
 
-      ! Ice mass and temp variables
-!CDIR NOVERRCHK
-      DO jj = k_j1 , k_jpj
+!CDIR NOVERRCHK
+      DO jj = k_j1 , k_jpj       ! Ice mass and temp variables
 !CDIR NOVERRCHK
          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.
+            zpresh(ji,jj) = tms(ji,jj) *  strength(ji,jj) * 0.5_wp
 #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) - &
-               epsd ) ) ) * tms(ji,jj)
+            tmi(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - epsd ) ) ) * tms(ji,jj)
          END DO
       END DO
@@ -271,24 +238,20 @@
          DO ji = fs_2, fs_jpim1
 
-            zt11 = tms(ji,jj)*e1t(ji,jj)
-            zt12 = tms(ji+1,jj)*e1t(ji+1,jj)
-            zt21 = tms(ji,jj)*e2t(ji,jj)
-            zt22 = tms(ji,jj+1)*e2t(ji,jj+1)
+            zt11 = tms(ji  ,jj) * e1t(ji  ,jj)
+            zt12 = tms(ji+1,jj) * e1t(ji+1,jj)
+            zt21 = tms(ji,jj  ) * e2t(ji,jj  )
+            zt22 = tms(ji,jj+1) * e2t(ji,jj+1)
 
             ! Leads area.
-            zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + &
-               &                        zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd )
-            zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + &
-               &                        zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd )
+            zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + epsd )
+            zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + epsd )
 
             ! Mass, coriolis coeff. and currents
             zmass1(ji,jj) = ( zt12*zc1(ji,jj) + zt11*zc1(ji+1,jj) ) / (zt11+zt12+epsd)
             zmass2(ji,jj) = ( zt22*zc1(ji,jj) + zt21*zc1(ji,jj+1) ) / (zt21+zt22+epsd)
-            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + &
-               e1t(ji,jj)*fcor(ji+1,jj) ) &
-               / (e1t(ji,jj) + e1t(ji+1,jj) + epsd )
-            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + &
-               e2t(ji,jj)*fcor(ji,jj+1) ) &
-               / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd )
+            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj)*fcor(ji,jj) + e1t(ji,jj)*fcor(ji+1,jj) )   &
+               &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + epsd )
+            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1)*fcor(ji,jj) + e2t(ji,jj)*fcor(ji,jj+1) )   &
+               &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + epsd )
             !
             u_oce1(ji,jj)  = u_oce(ji,jj)
@@ -314,6 +277,6 @@
             ! include it later
 
-            zdsshx =  (ssh_m(ji+1,jj) - ssh_m(ji,jj))/e1u(ji,jj)
-            zdsshy =  (ssh_m(ji,jj+1) - ssh_m(ji,jj))/e2v(ji,jj)
+            zdsshx =  ( ssh_m(ji+1,jj) - ssh_m(ji,jj) ) / e1u(ji,jj)
+            zdsshy =  ( ssh_m(ji,jj+1) - ssh_m(ji,jj) ) / e2v(ji,jj)
 
             za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
@@ -330,19 +293,19 @@
       ! Time step for subcycling
       dtevp  = rdt_ice / nevp
-      dtotel = dtevp / ( 2.0 * telast )
+      dtotel = dtevp / ( 2._wp * telast )
 
       !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
-      ecc2 = ecc*ecc
+      ecc2 = ecc * ecc
 
       !-Initialise stress tensor 
-      zs1(:,:)  = stress1_i(:,:) 
-      zs2(:,:)  = stress2_i(:,:)
+      zs1 (:,:) = stress1_i (:,:) 
+      zs2 (:,:) = stress2_i (:,:)
       zs12(:,:) = stress12_i(:,:)
 
-      !----------------------!
+      !                                               !----------------------!
       DO jter = 1 , nevp                              !    loop over jter    !
-         !----------------------!        
+         !                                            !----------------------!        
          DO jj = k_j1, k_jpj-1
-            zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
+            zu_ice(:,jj) = u_ice(:,jj)    ! velocity at previous time step
             zv_ice(:,jj) = v_ice(:,jj)
          END DO
@@ -409,6 +372,5 @@
             END DO
          END DO
-         CALL lbc_lnk( v_ice1(:,:), 'U', -1. )
-         CALL lbc_lnk( u_ice2(:,:), 'V', -1. )
+         CALL lbc_lnk( v_ice1, 'U', -1. )   ;   CALL lbc_lnk( u_ice2, 'V', -1. )      ! lateral boundary cond.
 
 !CDIR NOVERRCHK
@@ -418,15 +380,13 @@
 
                !- Calculate Delta at centre of grid cells
-               zdst       = (  e2u( ji  , jj   ) * v_ice1(ji,jj)            &
-                  &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)          &
-                  &          + e1v( ji  , jj   ) * u_ice2(ji,jj)            &
-                  &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1)          &
-                  &          )                                              &
+               zdst       = (  e2u(ji  , jj) * v_ice1(ji  ,jj)          &
+                  &          - e2u(ji-1, jj) * v_ice1(ji-1,jj)          &
+                  &          + e1v(ji, jj  ) * u_ice2(ji,jj  )          &
+                  &          - e1v(ji, jj-1) * u_ice2(ji,jj-1)          &
+                  &          )                                          &
                   &         / area(ji,jj)
 
-               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) +                        & 
-                  &                       ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )  
-               deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 +                    &
-                  (zdt(ji,jj)**2 + zdst**2)*usecc2), creepl )
+               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst*zdst ) * usecc2 )  
+               deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zdst**2)*usecc2), creepl )
 
                !-Calculate stress tensor components zs1 and zs2 
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 2370)
@@ -4,7 +4,9 @@
    !!           computation of the flux at the sea ice/ocean interface
    !!======================================================================
-   !! History :   -   !  2006-07 (M. Vancoppelle)  LIM3 original code
-   !!            3.0  !  2008-03 (C. Tallandier)  surface module
-   !!             -   !  2008-04 (C. Tallandier)  split in 2 + new ice-ocean coupling
+   !! History :   -   ! 2006-07 (M. Vancoppelle)  LIM3 original code
+   !!            3.0  ! 2008-03 (C. Tallandier)  surface module
+   !!             -   ! 2008-04 (C. Tallandier)  split in 2 + new ice-ocean coupling
+   !!            3.3  ! 2010-05 (G. Madec) decrease ocean & ice reference salinities in the Baltic sea
+   !!                 !                  + simplification of the ice-ocean stress calculation
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -12,7 +14,7 @@
    !!   'key_lim3'                                    LIM 3.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!   lim_sbc  : flux at the ice / ocean interface
-   !!----------------------------------------------------------------------
-   USE oce
+   !!   lim_sbc_flx  : updates mass, heat and salt fluxes at the ocean surface
+   !!   lim_sbc_tau  : update i- and j-stresses, and its modulus at the ocean surface
+   !!----------------------------------------------------------------------
    USE par_oce          ! ocean parameters
    USE par_ice          ! ice parameters
@@ -35,11 +37,13 @@
    PUBLIC   lim_sbc_tau   ! called by sbc_ice_lim
 
-   REAL(wp)  ::   epsi16 = 1.e-16  ! constant values
-   REAL(wp)  ::   rzero  = 0.e0    
-   REAL(wp)  ::   rone   = 1.e0
-
-   REAL(wp), DIMENSION(jpi,jpj) ::   utau_oce, vtau_oce   !: air-ocean surface i- & j-stress              [N/m2]
-   REAL(wp), DIMENSION(jpi,jpj) ::   tmod_io              !: modulus of the ice-ocean relative velocity   [m/s]
-   REAL(wp), DIMENSION(jpi,jpj) ::   ssu_mb  , ssv_mb     !: before mean ocean surface currents           [m/s]
+   REAL(wp)  ::   r1_rdtice            ! = 1. / rdt_ice 
+   REAL(wp)  ::   epsi16 = 1.e-16_wp   ! constant values
+   REAL(wp)  ::   rzero  = 0._wp    
+   REAL(wp)  ::   rone   = 1._wp
+
+   REAL(wp), DIMENSION(jpi,jpj) ::   utau_oce, vtau_oce   ! air-ocean surface i- & j-stress              [N/m2]
+   REAL(wp), DIMENSION(jpi,jpj) ::   tmod_io              ! modulus of the ice-ocean relative velocity   [m/s]
+
+   REAL(wp), DIMENSION(jpi,jpj) ::   soce_0, sice_0   ! constant SSS and ice salinity used in levitating sea-ice case
 
    !! * Substitutions
@@ -48,168 +52,7 @@
    !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
-
-   SUBROUTINE lim_sbc_tau( kt, kcpl )
-      !!-------------------------------------------------------------------
-      !!                ***  ROUTINE lim_sbc_tau ***
-      !!  
-      !! ** Purpose : Update the ocean surface stresses due to the ice
-      !!         
-      !! ** Action  : - compute the ice-ocean stress depending on kcpl:
-      !!              fluxes at the ice-ocean interface.
-      !!     Case 0  :  old LIM-3 way, computed at ice time-step only
-      !!     Case 1  :  at each ocean time step the stresses are computed
-      !!                using the current ocean velocity (now)
-      !!     Case 2  :  combination of half case 0 + half case 1
-      !!     
-      !! ** Outputs : - utau   : sea surface i-stress (ocean referential)
-      !!              - vtau   : sea surface j-stress (ocean referential)
-      !!
-      !! References : Goosse, H. et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90.
-      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
-      !!---------------------------------------------------------------------
-      INTEGER ::   kt     ! number of ocean iteration
-      INTEGER ::   kcpl   ! ice-ocean coupling indicator: =0  LIM-3 old case
-      !                   !                               =1  stresses computed using now ocean velocity
-      !                   !                               =2  combination of 0 and 1 cases
-      !!
-      INTEGER  ::   ji, jj   ! dummy loop indices
-      REAL(wp) ::   zfrldu, zat_u, zu_ico, zutaui, zu_u, zu_ij, zu_im1j   ! temporary scalar
-      REAL(wp) ::   zfrldv, zat_v, zv_ico, zvtaui, zv_v, zv_ij, zv_ijm1   !    -         -
-      REAL(wp) ::   zsang, zztmp                                          !    -         -
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice
-#if defined key_coupled    
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb     ! albedo of ice under overcast sky
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalbp    ! albedo of ice under clear sky
-#endif
-     !!---------------------------------------------------------------------
-
-      IF( kt == nit000 ) THEN
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'lim_sbc_tau : LIM 3.0 sea-ice - surface ocean momentum fluxes'
-         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
-      ENDIF
-
-      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-!CDIR NOVERRCHK
-         DO jj = 2, jpjm1                         !* modulus of the ice-ocean velocity
-!CDIR NOVERRCHK
-            DO ji = fs_2, fs_jpim1
-               zu_ij   = u_ice(ji  ,jj) - ssu_m(ji  ,jj)               ! (i  ,j)
-               zu_im1j = u_ice(ji-1,jj) - ssu_m(ji-1,jj)               ! (i-1,j)
-               zv_ij   = v_ice(ji,jj  ) - ssv_m(ji,jj  )               ! (i,j  )
-               zv_ijm1 = v_ice(ji,jj-1) - ssv_m(ji,jj-1)               ! (i,j-1)
-               zztmp =  0.5 * ( zu_ij * zu_ij + zu_im1j * zu_im1j   &  ! mean of the square values instead
-                  &           + zv_ij * zv_ij + zv_ijm1 * zv_ijm1 )    ! of the square of the mean values...
-               ! WARNING, here we make an approximation for case 1 and 2: taum is not computed at each time
-               ! step but only every nn_fsbc time step. This avoid mutiple avarage to pass from T -> U,V grids
-               ! and next from U,V grids to T grid. Taum is used in tke, which should not be too sensitive to
-               ! this approximaton...
-               taum(ji,jj) = ( 1. - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * rhoco * zztmp
-               tmod_io(ji,jj) = SQRT( zztmp ) 
-            END DO
-         END DO
-         CALL lbc_lnk( taum, 'T', 1. )   ;   CALL lbc_lnk( tmod_io, 'T', 1. )
-      ENDIF
-
-      SELECT CASE( kcpl )
-         !                                        !--------------------------------!
-      CASE( 0 )                                   !  LIM 3 old stress computation  !  (at ice timestep only)
-         !                                        !--------------------------------! 
-         !                                             !* ice stress over ocean with a ice-ocean rotation angle
-         DO jj = 1, jpjm1
-            zsang  = SIGN( 1.e0, gphif(1,jj) ) * sangvg         ! change the sinus angle sign in the south hemisphere
-            DO ji = 1, fs_jpim1
-               zu_u = u_ice(ji,jj) - u_oce(ji,jj)               ! ice velocity relative to the ocean
-               zv_v = v_ice(ji,jj) - v_oce(ji,jj)
-               !                                                ! quadratic drag formulation with rotation
-!!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1)
-               zutaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_u - zsang * zv_v ) 
-               zvtaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_v + zsang * zu_u ) 
-               !                                                ! bound for too large stress values
-               ! IMPORTANT: the exponential below prevents numerical oscillations in the ice-ocean stress
-               ! This is not physically based. A cleaner solution is offer in CASE kcpl=2
-               ztio_u(ji,jj) = zutaui * EXP( - ( tmod_io(ji,jj) + tmod_io(ji+1,jj) )  )
-               ztio_v(ji,jj) = zvtaui * EXP( - ( tmod_io(ji,jj) + tmod_io(ji,jj+1) )  ) 
-            END DO
-         END DO
-         DO jj = 2, jpjm1                                       ! stress at the surface of the ocean
-            DO ji = fs_2, fs_jpim1   ! vertor opt.
-               zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj  ) )   ! open-ocean fraction at U- & V-points (from T-point values)
-               zfrldv = 0.5 * ( ato_i(ji,jj) + ato_i(ji  ,jj+1) )
-               !                                                    ! update surface ocean stress
-               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * ztio_u(ji,jj)
-               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * ztio_v(ji,jj)
-            END DO
-         END DO
-         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
-
-         !
-         !                                        !--------------------------------!
-      CASE( 1 )                                   !  Use the now velocity          !  (at each ocean timestep)
-         !                                        !--------------------------------! 
-         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-            utau_oce(:,:) = utau(:,:)                 !* save the air-ocean stresses at ice time-step
-            vtau_oce(:,:) = vtau(:,:)
-         ENDIF
-         !
-        DO jj = 2, jpjm1                              !* ice stress over ocean with a ice-ocean rotation angle
-            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg        ! change the sinus angle sign in the south hemisphere
-            DO ji = fs_2, fs_jpim1
-               zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5   ! ice area at u and V-points
-               zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5
-               !                                                ! (u,v) ice-ocean velocity at (U,V)-point, resp.
-               zu_u   = u_ice(ji,jj) - ub(ji,jj,1)
-               zv_v   = v_ice(ji,jj) - vb(ji,jj,1)
-               !                                                ! quadratic drag formulation with rotation
-!!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1)
-               zutaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_u - zsang * zv_v ) 
-               zvtaui   = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_v + zsang * zu_u ) 
-               !                                                   ! stress at the ocean surface
-               utau(ji,jj) = ( 1.- zat_u ) * utau_oce(ji,jj) + zat_u * zutaui
-               vtau(ji,jj) = ( 1.- zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui
-            END DO
-         END DO
-         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
-
-         !
-         !                                        !--------------------------------!
-      CASE( 2 )                                   !  mixed 0 and 2 cases           !  (at each ocean timestep)
-         !                                        !--------------------------------! 
-         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-            utau_oce(:,:) = utau (:,:)               !* save the air-ocean stresses at ice time-step
-            vtau_oce(:,:) = vtau (:,:)
-            ssu_mb  (:,:) = ssu_m(:,:)               !* save the ice-ocean velocity at ice time-step
-            ssv_mb  (:,:) = ssv_m(:,:)
-         ENDIF
-         DO jj = 2, jpjm1                            !* ice stress over ocean with a ice-ocean rotation angle
-            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
-            DO ji = fs_2, fs_jpim1
-               zat_u = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5     ! ice area at u and V-points
-               zat_v = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5 
-               !
-               zu_ico = u_ice(ji,jj) - 0.5 * ( ub(ji,jj,1) + ssu_mb(ji,jj) )   ! ice-oce velocity using ub and ssu_mb
-               zv_ico = v_ice(ji,jj) - 0.5 * ( vb(ji,jj,1) + ssv_mb(ji,jj) )
-               !                                        ! quadratic drag formulation with rotation
-!!gm still an error in the rotation, but usually the angle is zero (zsang=0, cangvg=1)
-               zutaui = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * rhoco * ( cangvg * zu_ico - zsang * zv_ico )
-               zvtaui = 0.5 * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * rhoco * ( cangvg * zv_ico + zsang * zu_ico )
-               !
-               utau(ji,jj) = ( 1.-zat_u ) * utau_oce(ji,jj) + zat_u * zutaui    ! stress at the ocean surface
-               vtau(ji,jj) = ( 1.-zat_v ) * vtau_oce(ji,jj) + zat_v * zvtaui
-            END DO
-         END DO
-         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
-         !
-      END SELECT
-
-      IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
-         &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask )
-      !  
-   END SUBROUTINE lim_sbc_tau
-
 
    SUBROUTINE lim_sbc_flx( kt )
@@ -235,5 +78,5 @@
       !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
       !!---------------------------------------------------------------------
-      INTEGER ::   kt    ! number of iteration
+      INTEGER, INTENT(in) ::   kt    ! number of iteration
       !!
       INTEGER  ::   ji, jj           ! dummy loop indices
@@ -254,4 +97,18 @@
          IF(lwp) WRITE(numout,*) 'lim_sbc_flx : LIM 3.0 sea-ice - heat salt and mass ocean surface fluxes'
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
+         !
+         r1_rdtice = 1. / rdt_ice
+         !
+         soce_0(:,:) = soce
+         sice_0(:,:) = sice
+         !
+         IF( cp_cfg == "orca" ) THEN           ! decrease ocean & ice reference salinities in the Baltic sea 
+            WHERE( 14._wp <= glamt(:,:) .AND. glamt(:,:) <= 32._wp .AND.   &
+               &   54._wp <= gphit(:,:) .AND. gphit(:,:) <= 66._wp         ) 
+               soce_0(:,:) = 4._wp
+               sice_0(:,:) = 2._wp
+            END WHERE
+         ENDIF
+         !
       ENDIF
 
@@ -298,6 +155,6 @@
                ! fscmbq and ffltbif are obsolete
                !              &           + iflt * ffltbif(ji,jj) !!! only if one category is used
-               &           + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) / rdt_ice   &
-               &           + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) / rdt_ice                     &
+               &           + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   &
+               &           + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) )                   * r1_rdtice   &
                &           + fhmec(ji,jj)     & ! new contribution due to snow melt in ridging!!
                &           + fheat_rpo(ji,jj) & ! contribution from ridge formation
@@ -340,11 +197,11 @@
                WRITE(numout,*) ' qcmif     : ', qcmif(jiindx,jjindx)
                WRITE(numout,*) ' qldif     : ', qldif(jiindx,jjindx)
-               WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) / rdt_ice
-               WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) / rdt_ice
+               WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice
+               WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice
                WRITE(numout,*) ' ifrdv     : ', ifrdv
                WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindx,jjindx)
                WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindx,jjindx)
-               WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) / rdt_ice
-               WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) / rdt_ice
+               WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice
+               WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice
                WRITE(numout,*) ' '
                WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx)
@@ -380,5 +237,5 @@
                !              &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )  &   !  remov. snow precip over ice
                &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   !  remov. snow precip over ice
-               &   - rdmsnif(ji,jj) / rdt_ice                &   !  freshwaterflux due to snow melting 
+               &   - rdmsnif(ji,jj) * r1_rdtice                       &   !  freshwaterflux due to snow melting 
                ! new contribution from snow falling when ridging
                &   + fmmec(ji,jj)
@@ -386,7 +243,7 @@
             !  computing salt exchanges at the ice/ocean interface
             !  sice should be the same as computed with the ice model
-            zfons =  ( soce - sice ) * ( rdmicif(ji,jj) / rdt_ice ) 
+            zfons =  ( soce_0(ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice 
             ! SOCE
-            zfons =  ( sss_m(ji,jj) - sice ) * ( rdmicif(ji,jj) / rdt_ice ) 
+            zfons =  ( sss_m (ji,jj) - sice_0(ji,jj) ) * rdmicif(ji,jj) * r1_rdtice
 
             !CT useless            !  salt flux for constant salinity
@@ -446,4 +303,87 @@
    END SUBROUTINE lim_sbc_flx
 
+
+   SUBROUTINE lim_sbc_tau( kt , pu_oce, pv_oce )
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE lim_sbc_tau ***
+      !!  
+      !! ** Purpose : Update the ocean surface stresses due to the ice
+      !!         
+      !! ** Action  : * at each ice time step (every nn_fsbc time step):
+      !!                - compute the modulus of ice-ocean relative velocity 
+      !!                  (*rho*Cd) at T-point (C-grid) or I-point (B-grid)
+      !!                      tmod_io = rhoco * | U_ice-U_oce |
+      !!                - update the modulus of stress at ocean surface
+      !!                      taum = frld * taum + (1-frld) * tmod_io * | U_ice-U_oce |
+      !!              * at each ocean time step (every kt): 
+      !!                  compute linearized ice-ocean stresses as
+      !!                      Utau = tmod_io * | U_ice - pU_oce |
+      !!                using instantaneous current ocean velocity (usually before)
+      !!
+      !!    NB: - ice-ocean rotation angle no more allowed
+      !!        - here we make an approximation: taum is only computed every ice time step
+      !!          This avoids mutiple average to pass from T -> U,V grids and next from U,V grids 
+      !!          to T grid. taum is used in TKE and GLS, which should not be too sensitive to this approximaton...
+      !!
+      !! ** Outputs : - utau, vtau   : surface ocean i- and j-stress (u- & v-pts) updated with ice-ocean fluxes
+      !!              - taum         : modulus of the surface ocean stress (T-point) updated with ice-ocean fluxes
+      !!---------------------------------------------------------------------
+      INTEGER ,                     INTENT(in) ::   kt               ! ocean time-step index
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(in) ::   pu_oce, pv_oce   ! surface ocean currents
+      !!
+      INTEGER  ::   ji, jj   ! dummy loop indices
+      REAL(wp) ::   zat_u, zutau_ice, zu_t, zmodt   ! local scalar
+      REAL(wp) ::   zat_v, zvtau_ice, zv_t          !   -      -
+     !!---------------------------------------------------------------------
+
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'lim_sbc_tau : LIM-3 sea-ice - surface ocean momentum fluxes'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
+      ENDIF
+
+      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !==  Ice time-step only  ==!   (i.e. surface module time-step)
+!CDIR NOVERRCHK
+         DO jj = 2, jpjm1                             !* update the modulus of stress at ocean surface (T-point)
+!CDIR NOVERRCHK
+            DO ji = fs_2, fs_jpim1
+               !                                               ! 2*(U_ice-U_oce) at T-point
+               zu_t = u_ice(ji,jj) + u_ice(ji-1,jj) - u_oce(ji,jj) - u_oce(ji-1,jj)   
+               zv_t = v_ice(ji,jj) + v_ice(ji,jj-1) - v_oce(ji,jj) - v_oce(ji,jj-1) 
+               !                                              ! |U_ice-U_oce|^2
+               zmodt =  0.25_wp * (  zu_t * zu_t + zv_t * zv_t  )
+               !                                               ! update the ocean stress modulus
+               taum(ji,jj) = ( 1._wp - at_i(ji,jj) ) * taum(ji,jj) + at_i(ji,jj) * rhoco * zmodt
+               tmod_io(ji,jj) = rhoco * SQRT( zmodt )          ! rhoco * |U_ice-U_oce| at T-point
+            END DO
+         END DO
+         CALL lbc_lnk( taum, 'T', 1. )   ;   CALL lbc_lnk( tmod_io, 'T', 1. )
+         !
+         utau_oce(:,:) = utau(:,:)                    !* save the air-ocean stresses at ice time-step
+         vtau_oce(:,:) = vtau(:,:)
+         !
+      ENDIF
+         !
+         !                                     !==  every ocean time-step  ==!
+         !
+      DO jj = 2, jpjm1                                !* update the stress WITHOUT a ice-ocean rotation angle
+         DO ji = fs_2, fs_jpim1   ! Vect. Opt.
+            zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5_wp   ! ice area at u and V-points
+            zat_v  = ( at_i(ji,jj) + at_i(ji,jj+1) ) * 0.5_wp
+            !                                                   ! linearized quadratic drag formulation
+            zutau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji+1,jj) ) * ( u_ice(ji,jj) - pu_oce(ji,jj) )
+            zvtau_ice   = 0.5_wp * ( tmod_io(ji,jj) + tmod_io(ji,jj+1) ) * ( v_ice(ji,jj) - pv_oce(ji,jj) )
+            !                                                   ! stresses at the ocean surface
+            utau(ji,jj) = ( 1._wp - zat_u ) * utau_oce(ji,jj) + zat_u * zutau_ice
+            vtau(ji,jj) = ( 1._wp - zat_v ) * vtau_oce(ji,jj) + zat_v * zvtau_ice
+         END DO
+      END DO
+      CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
+      !
+      IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
+         &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask )
+      !  
+   END SUBROUTINE lim_sbc_tau
+
 #else
    !!----------------------------------------------------------------------
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90	(revision 2370)
@@ -9,5 +9,5 @@
 #if defined key_lim3 || defined key_lim2
    !!----------------------------------------------------------------------
-   !!   'key_lim2' or 'key_lim3' :             LIM 2.0 or 3.0 sea-ice model
+   !!   'key_lim2' or 'key_lim3' :             LIM-2 or LIM-3 sea-ice model
    !!----------------------------------------------------------------------
    USE par_oce          ! ocean parameters
@@ -23,12 +23,16 @@
 
 # if defined  key_lim2
-   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim2        = .TRUE.    !: LIM-2 ice model
-   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3        = .FALSE.   !: no LIM-3
-   CHARACTER(len=1), PUBLIC            ::   cigr_type      = 'I'       !: 'I'-grid ice-velocity (B-grid lower left corner)
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim2    = .TRUE.   !: LIM-2 ice model
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3    = .FALSE.  !: no LIM-3
+#  if defined key_lim2_vp
+   CHARACTER(len=1), PUBLIC, PARAMETER ::   cp_ice_msh = 'I'      !: VP : 'I'-grid ice-velocity (B-grid lower left corner)
+#  else
+   CHARACTER(len=1), PUBLIC, PARAMETER ::   cp_ice_msh = 'C'      !: EVP: 'C'-grid ice-velocity
+#  endif
 # endif
 # if defined  key_lim3
-   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim2        = .FALSE.   !: no LIM-2
-   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3        = .TRUE.    !: LIM-3 ice model
-   CHARACTER(len=1), PUBLIC            ::   cigr_type      = 'C'       !: 'C'-grid ice-velocity
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim2    = .FALSE.  !: no LIM-2
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3    = .TRUE.   !: LIM-3 ice model
+   CHARACTER(len=1), PUBLIC, PARAMETER ::   cp_ice_msh = 'C'      !: 'C'-grid ice-velocity
 # endif
 
@@ -41,8 +45,8 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpl) ::   alb_ice   !: albedo of ice
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau_ice    !: u-stress over ice (I-point for LIM2 or U,V-point for LIM3)   [N/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau_ice    !: v-stress over ice (I-point for LIM2 or U,V-point for LIM3)   [N/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr1_i0      !: 1st fraction of sol. rad. which penetrate inside the ice cover
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr2_i0      !: 2nd fraction of sol. rad. which penetrate inside the ice cover
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau_ice    !: u-stress over ice (I-pt for VP or U,V-pts for EVP)   [N/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   vtau_ice    !: v-stress over ice (I-pt for VP or U,V-pts for EVP)   [N/m2]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr1_i0      !: 1st fraction of Qsr which penetrates inside the ice cover
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fr2_i0      !: 2nd fraction of Qsr which penetrates inside the ice cover
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_ice     !: solid freshwater budget over ice: sublivation - snow
 
@@ -55,7 +59,7 @@
    !!   Default option                      NO LIM 2.0 or 3.0 sea-ice model
    !!----------------------------------------------------------------------
-   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim2        = .FALSE.  !: no LIM-2 ice model
-   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3        = .FALSE.  !: no LIM-3 ice model
-   CHARACTER(len=1), PUBLIC            ::   cigr_type      = '-'      !: no grid ice-velocity
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim2    = .FALSE.  !: no LIM-2 ice model
+   LOGICAL         , PUBLIC, PARAMETER ::   lk_lim3    = .FALSE.  !: no LIM-3 ice model
+   CHARACTER(len=1), PUBLIC, PARAMETER ::   cp_ice_msh = '-'      !: no grid ice-velocity
 #endif
 
@@ -63,7 +67,5 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE sbc_ice
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 2370)
@@ -4,7 +4,8 @@
    !! Surface module :   variables defined in core memory 
    !!======================================================================
-   !! History :  3.0   !  2006-06  (G. Madec)  Original code
-   !!             -    !  2008-08  (G. Madec)  namsbc moved from sbcmod
-   !!            3.3   !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
+   !! History :  3.0  ! 2006-06  (G. Madec)  Original code
+   !!             -   ! 2008-08  (G. Madec)  namsbc moved from sbcmod
+   !!            3.3  ! 2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
+   !!             -   ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
    !!----------------------------------------------------------------------
    USE par_oce          ! ocean parameters
@@ -30,13 +31,9 @@
    !                                             !:  = 1 global mean of e-p-r set to zero at each nn_fsbc time step
    !                                             !:  = 2 annual global mean of e-p-r set to zero
-   INTEGER , PUBLIC ::   nn_ico_cpl  = 0          !: ice-ocean coupling indicator
-   !                                             !:  = 0   LIM-3 old case
-   !                                             !:  = 1   stresses computed using now ocean velocity
-   !                                             !:  = 2   combination of 0 and 1 cases
 
    !!----------------------------------------------------------------------
    !!              Ocean Surface Boundary Condition fields
    !!----------------------------------------------------------------------
-   LOGICAL , PUBLIC ::   lhftau = .FALSE.              !: HF tau used in TKE: mean(stress module) - module(mean stress)
+   LOGICAL , PUBLIC ::   lhftau = .FALSE.        !: HF tau used in TKE: mean(stress module) - module(mean stress)
    !!                                   !!   now    ! before   !!
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   utau   , utau_b   !: sea surface i-stress (ocean referential)     [N/m2]
@@ -53,8 +50,8 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp_tot           !: total E-P over ocean and ice                 [Kg/m2/s]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rnf    , rnf_b    !: river runoff   [Kg/m2/s]  
-   ! - ML - begin
+   !!
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpts) ::  sbc_tsc, sbc_tsc_b  !: sbc content trend                      [K.m/s]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpk) ::   qsr_hc , qsr_hc_b   !: heat content trend due to qsr flux     [K.m/s]
-   ! - ML - end
+   !!
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tprecip           !: total precipitation                          [Kg/m2/s]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sprecip           !: solid precipitation                          [Kg/m2/s]
@@ -63,6 +60,4 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   atm_co2           !: atmospheric pCO2                             [ppm]
 #endif
-!!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rrunoff        !: runoff
-!!$   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   calving        !: calving
 
    !!----------------------------------------------------------------------
@@ -79,5 +74,5 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE sbc_oce
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_clio.F90	(revision 2370)
@@ -83,7 +83,6 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$ 
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -137,13 +136,13 @@
 
          ! (NB: frequency positive => hours, negative => months)
-         !            !    file     ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   !
-         !            !    name     !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      !
-         sn_utau = FLD_N( 'utau'    ,    24     ,  'utau'    ,  .true.    , .false. ,   'yearly'  , ''       , ''         ) 
-         sn_vtau = FLD_N( 'vtau'    ,    24     ,  'vtau'    ,  .true.    , .false. ,   'yearly'  , ''       , ''         ) 
-         sn_wndm = FLD_N( 'mwnd10m' ,    24     ,  'm_10'    ,  .true.    , .false. ,   'yearly'  , ''       , ''         ) 
-         sn_tair = FLD_N( 'tair10m' ,    24     ,  't_10'    ,  .false.   , .false. ,   'yearly'  , ''       , ''         ) 
-         sn_humi = FLD_N( 'humi10m' ,    24     ,  'q_10'    ,  .false.   , .false. ,   'yearly'  , ''       , ''         ) 
-         sn_ccov = FLD_N( 'ccover'  ,    -1     ,  'cloud'   ,  .true.    , .false. ,   'yearly'  , ''       , ''         ) 
-         sn_prec = FLD_N( 'precip'  ,    -1     ,  'precip'  ,  .true.    , .false. ,   'yearly'  , ''       , ''         ) 
+         !            !    file    ! frequency ! variable ! time intep !  clim   ! 'yearly' or ! weights  ! rotation !
+         !            !    name    !  (hours)  !  name    !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
+         sn_utau = FLD_N( 'utau'   ,    24     , 'utau'   ,  .true.    , .false. ,   'yearly'  , ''       , ''       ) 
+         sn_vtau = FLD_N( 'vtau'   ,    24     , 'vtau'   ,  .true.    , .false. ,   'yearly'  , ''       , ''       ) 
+         sn_wndm = FLD_N( 'mwnd10m',    24     , 'm_10'   ,  .true.    , .false. ,   'yearly'  , ''       , ''       ) 
+         sn_tair = FLD_N( 'tair10m',    24     , 't_10'   ,  .false.   , .false. ,   'yearly'  , ''       , ''       ) 
+         sn_humi = FLD_N( 'humi10m',    24     , 'q_10'   ,  .false.   , .false. ,   'yearly'  , ''       , ''       ) 
+         sn_ccov = FLD_N( 'ccover' ,    -1     , 'cloud'  ,  .true.    , .false. ,   'yearly'  , ''       , ''       ) 
+         sn_prec = FLD_N( 'precip' ,    -1     , 'precip' ,  .true.    , .false. ,   'yearly'  , ''       , ''       ) 
 
          REWIND( numnam )                    ! ... read in namlist namsbc_clio
@@ -473,5 +472,5 @@
       INTEGER  ::   ijpl          ! number of ice categories (size of 3rd dim of input arrays)
       !!
-      REAL(wp) ::   zcoef, zmt1, zmt2, zmt3, ztatm3             ! temporary scalars
+      REAL(wp) ::   zcoef, zmt1, zmt2, zmt3, ztatm3     ! temporary scalars
       REAL(wp) ::   ztaevbk, zind1, zind2, zind3, ztamr         !    -         -
       REAL(wp) ::   zesi, zqsati, zdesidt                       !    -         -
@@ -496,6 +495,5 @@
       SELECT CASE( cd_grid )
       CASE( 'C' )                          ! C-grid ice dynamics
-         ! Change from wind speed to wind stress over OCEAN (cao is used)
-         zcoef = cai / cao 
+      zcoef  = cai / cao                         ! Change from air-sea stress to air-ice stress
 !CDIR COLLAPSE
          DO jj = 1 , jpj
@@ -505,15 +503,13 @@
             END DO
          END DO
-      CASE( 'B' )                          ! B-grid ice dynamics
-         ! stress from ocean U- and V-points to ice U,V point
-!CDIR COLLAPSE
-         DO jj = 2, jpj
-            DO ji = 2, jpi   ! B grid : no vector opt.
-               p_taui(ji,jj) = 0.5 * ( utau(ji-1,jj  ) + utau(ji-1,jj-1) )
-               p_tauj(ji,jj) = 0.5 * ( vtau(ji  ,jj-1) + vtau(ji-1,jj-1) )
+      CASE( 'I' )                          ! I-grid ice dynamics:  I-point (i.e. F-point lower-left corner)
+      zcoef  = 0.5_wp * cai / cao                ! Change from air-sea stress to air-ice stress!CDIR COLLAPSE
+         DO jj = 2, jpj         ! stress from ocean U- and V-points to ice U,V point
+            DO ji = 2, jpi   ! I-grid : no vector opt.
+               p_taui(ji,jj) = zcoef * ( utau(ji-1,jj  ) + utau(ji-1,jj-1) )
+               p_tauj(ji,jj) = zcoef * ( vtau(ji  ,jj-1) + vtau(ji-1,jj-1) )
             END DO
          END DO
-         CALL lbc_lnk( p_taui(:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
-         CALL lbc_lnk( p_tauj(:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
+         CALL lbc_lnk( p_taui(:,:), 'I', -1. )   ;   CALL lbc_lnk( p_tauj(:,:), 'I', -1. )   ! I-point
       END SELECT
 
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk_core.F90	(revision 2370)
@@ -75,5 +75,5 @@
    !! NEMO/OPA 3.3 , NEMO-consortium (2010) 
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
@@ -133,15 +133,15 @@
          !
          ! (NB: frequency positive => hours, negative => months)
-         !            !    file     ! frequency !  variable  ! time intep !  clim   ! 'yearly' or ! weights  ! rotation   !
-         !            !    name     !  (hours)  !   name     !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs      !
-         sn_wndi = FLD_N( 'uwnd10m' ,    24     ,  'u_10'    ,  .false.   , .false. ,   'yearly'  , ''       , ''         )
-         sn_wndj = FLD_N( 'vwnd10m' ,    24     ,  'v_10'    ,  .false.   , .false. ,   'yearly'  , ''       , ''         )
-         sn_qsr  = FLD_N( 'qsw'     ,    24     ,  'qsw'     ,  .false.   , .false. ,   'yearly'  , ''       , ''         )
-         sn_qlw  = FLD_N( 'qlw'     ,    24     ,  'qlw'     ,  .false.   , .false. ,   'yearly'  , ''       , ''         )
-         sn_tair = FLD_N( 'tair10m' ,    24     ,  't_10'    ,  .false.   , .false. ,   'yearly'  , ''       , ''         )
-         sn_humi = FLD_N( 'humi10m' ,    24     ,  'q_10'    ,  .false.   , .false. ,   'yearly'  , ''       , ''         )
-         sn_prec = FLD_N( 'precip'  ,    -1     ,  'precip'  ,  .true.    , .false. ,   'yearly'  , ''       , ''         )
-         sn_snow = FLD_N( 'snow'    ,    -1     ,  'snow'    ,  .true.    , .false. ,   'yearly'  , ''       , ''         )
-         sn_tdif = FLD_N( 'taudif'  ,    24     ,  'taudif'  ,  .true.    , .false. ,   'yearly'  , ''       , ''         )
+         !            !    file    ! frequency ! variable ! time intep !  clim   ! 'yearly' or ! weights  ! rotation !
+         !            !    name    !  (hours)  !  name    !   (T/F)    !  (T/F)  !  'monthly'  ! filename ! pairs    !
+         sn_wndi = FLD_N( 'uwnd10m',    24     , 'u_10'   ,  .false.   , .false. ,   'yearly'  , ''       , ''       )
+         sn_wndj = FLD_N( 'vwnd10m',    24     , 'v_10'   ,  .false.   , .false. ,   'yearly'  , ''       , ''       )
+         sn_qsr  = FLD_N( 'qsw'    ,    24     , 'qsw'    ,  .false.   , .false. ,   'yearly'  , ''       , ''       )
+         sn_qlw  = FLD_N( 'qlw'    ,    24     , 'qlw'    ,  .false.   , .false. ,   'yearly'  , ''       , ''       )
+         sn_tair = FLD_N( 'tair10m',    24     , 't_10'   ,  .false.   , .false. ,   'yearly'  , ''       , ''       )
+         sn_humi = FLD_N( 'humi10m',    24     , 'q_10'   ,  .false.   , .false. ,   'yearly'  , ''       , ''       )
+         sn_prec = FLD_N( 'precip' ,    -1     , 'precip' ,  .true.    , .false. ,   'yearly'  , ''       , ''       )
+         sn_snow = FLD_N( 'snow'   ,    -1     , 'snow'   ,  .true.    , .false. ,   'yearly'  , ''       , ''       )
+         sn_tdif = FLD_N( 'taudif' ,    24     , 'taudif' ,  .true.    , .false. ,   'yearly'  , ''       , ''       )
          !
          REWIND( numnam )                          ! read in namlist namsbc_core
@@ -396,21 +396,21 @@
       !! caution : the net upward water flux has with mm/day unit
       !!---------------------------------------------------------------------
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)      ::   pst      ! ice surface temperature (>0, =rt0 over land) [Kelvin]
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)    ::   pui      ! ice surface velocity (i- and i- components      [m/s]
-      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj)    ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid)
-      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)      ::   palb     ! ice albedo (clear sky) (alb_ice_cs)               [%]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_taui   ! i- & j-components of surface ice stress        [N/m2]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid)
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_qns    ! non solar heat flux over ice (T-point)         [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_qsr    !     solar heat flux over ice (T-point)         [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_qla    ! latent    heat flux over ice (T-point)         [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_dqns   ! non solar heat sensistivity  (T-point)         [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(:,:,:)      ::   p_dqla   ! latent    heat sensistivity  (T-point)         [W/m2]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_tpr    ! total precipitation          (T-point)      [Kg/m2/s]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_spr    ! solid precipitation          (T-point)      [Kg/m2/s]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_fr1    ! 1sr fraction of qsr penetration in ice (T-point)  [%]
-      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj)    ::   p_fr2    ! 2nd fraction of qsr penetration in ice (T-point)  [%]
-      CHARACTER(len=1), INTENT(in   )                ::   cd_grid  ! ice grid ( C or B-grid)
-      INTEGER, INTENT(in   )                         ::   pdim     ! number of ice categories
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   pst      ! ice surface temperature (>0, =rt0 over land) [Kelvin]
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pui      ! ice surface velocity (i- and i- components      [m/s]
+      REAL(wp), INTENT(in   ), DIMENSION(jpi,jpj) ::   pvi      !    at I-point (B-grid) or U & V-point (C-grid)
+      REAL(wp), INTENT(in   ), DIMENSION(:,:,:)   ::   palb     ! ice albedo (clear sky) (alb_ice_cs)               [%]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_taui   ! i- & j-components of surface ice stress        [N/m2]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_tauj   !   at I-point (B-grid) or U & V-point (C-grid)
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qns    ! non solar heat flux over ice (T-point)         [W/m2]
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qsr    !     solar heat flux over ice (T-point)         [W/m2]
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_qla    ! latent    heat flux over ice (T-point)         [W/m2]
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_dqns   ! non solar heat sensistivity  (T-point)         [W/m2]
+      REAL(wp), INTENT(  out), DIMENSION(:,:,:)   ::   p_dqla   ! latent    heat sensistivity  (T-point)         [W/m2]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_tpr    ! total precipitation          (T-point)      [Kg/m2/s]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_spr    ! solid precipitation          (T-point)      [Kg/m2/s]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_fr1    ! 1sr fraction of qsr penetration in ice (T-point)  [%]
+      REAL(wp), INTENT(  out), DIMENSION(jpi,jpj) ::   p_fr2    ! 2nd fraction of qsr penetration in ice (T-point)  [%]
+      CHARACTER(len=1), INTENT(in   )             ::   cd_grid  ! ice grid ( C or B-grid)
+      INTEGER, INTENT(in   )                      ::   pdim     ! number of ice categories
       !!
       INTEGER  ::   ji, jj, jl    ! dummy loop indices
@@ -449,5 +449,5 @@
       ! ----------------------------------------------------------------------------- !
       SELECT CASE( cd_grid )
-      CASE( 'B' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation)
+      CASE( 'I' )                  ! B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation)
          !                           and scalar wind at T-point ( = | U10m - U_ice | ) (masked)
 #if defined key_vectopt_loop
@@ -605,6 +605,4 @@
       !!   9.0  !  05-08  (L. Brodeau) Rewriting and optimization
       !!----------------------------------------------------------------------
-      !! * Arguments
-
       REAL(wp), INTENT(in) :: zu                 ! altitude of wind measurement       [m]
       REAL(wp), INTENT(in), DIMENSION(jpi,jpj) ::  &
@@ -646,5 +644,5 @@
          grav   = 9.8,          &  ! gravity                       
          kappa  = 0.4              ! von Karman s constant
-
+      !!----------------------------------------------------------------------
       !! * Start
       !! Air/sea differences
@@ -770,5 +768,5 @@
          grav   = 9.8,      &  ! gravity                       
          kappa  = 0.4          ! von Karman's constant
-
+      !!----------------------------------------------------------------------
       !!  * Start
 
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90	(revision 2370)
@@ -4,7 +4,7 @@
    !! Surface Boundary Condition :  momentum, heat and freshwater fluxes in coupled mode
    !!======================================================================
-   !! History :  2.0  !  06-2007  (R. Redler, N. Keenlyside, W. Park) Original code split into flxmod & taumod
-   !!            3.0  !  02-2008  (G. Madec, C Talandier)  surface module
-   !!            3.1  !  02-2009  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
+   !! History :  2.0  ! 2007-06  (R. Redler, N. Keenlyside, W. Park) Original code split into flxmod & taumod
+   !!            3.0  ! 2008-02  (G. Madec, C Talandier)  surface module
+   !!            3.1  ! 2009_02  (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface
    !!----------------------------------------------------------------------
 #if defined key_oasis3 || defined key_oasis4
@@ -788,7 +788,7 @@
       !! ** Method  :   transform the received stress from the atmosphere into
       !!             an atmosphere-ice stress in the (i,j) ocean referencial
-      !!             and at the velocity point of the sea-ice model (cigr_type):
+      !!             and at the velocity point of the sea-ice model (cp_ice_msh):
       !!                'C'-grid : i- (j-) components given at U- (V-) point 
-      !!                'B'-grid : both components given at I-point 
+      !!                'I'-grid : B-grid lower-left corner: both components given at I-point 
       !!
       !!                The received stress are :
@@ -803,5 +803,5 @@
       !!                 first  as  2 components on the sphere 
       !!                 second as  2 components oriented along the local grid
-      !!                 third  as  2 components on the cigr_type point 
+      !!                 third  as  2 components on the cp_ice_msh point 
       !!
       !!                In 'oce and ice' case, only one vector stress field 
@@ -809,7 +809,7 @@
       !!             so that it is now defined as (i,j) components given at U-
       !!             and V-points, respectively. Therefore, here only the third
-      !!             transformation is done and only if the ice-grid is a 'B'-grid. 
-      !!
-      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice at cigr_type point
+      !!             transformation is done and only if the ice-grid is a 'I'-grid. 
+      !!
+      !! ** Action  :   return ptau_i, ptau_j, the stress over the ice at cp_ice_msh point
       !!----------------------------------------------------------------------
       REAL(wp), INTENT(out), DIMENSION(jpi,jpj) ::   p_taui   ! i- & j-components of atmos-ice stress [N/m2]
@@ -872,5 +872,5 @@
          !    
          !                                                  j+1   j     -----V---F
-         ! ice stress on ice velocity point (cigr_type)                  !       |
+         ! ice stress on ice velocity point (cp_ice_msh)                 !       |
          ! (C-grid ==>(U,V) or B-grid ==> I or F)                 j      |   T   U
          !                                                               |       |
@@ -879,5 +879,5 @@
          !                                                              i-1  i   i
          !                                                               i      i+1 (for I)
-         SELECT CASE ( cigr_type )
+         SELECT CASE ( cp_ice_msh )
             !
          CASE( 'I' )                                         ! B-grid ==> I
@@ -1258,5 +1258,5 @@
             END DO
          CASE( 'weighted oce and ice' )   
-            SELECT CASE ( cigr_type )
+            SELECT CASE ( cp_ice_msh )
             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
                DO jj = 2, jpjm1
@@ -1293,5 +1293,5 @@
             CALL lbc_lnk( zitx1, 'T', -1. )   ;   CALL lbc_lnk( zity1, 'T', -1. )
          CASE( 'mixed oce-ice'        )
-            SELECT CASE ( cigr_type )
+            SELECT CASE ( cp_ice_msh )
             CASE( 'C' )                      ! Ocean and Ice on C-grid ==> T
                DO jj = 2, jpjm1
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 2370)
@@ -5,8 +5,9 @@
    !!       &           covered area using LIM sea-ice model
    !! Sea-Ice model  :  LIM 3.0 Sea ice model time-stepping
-   !!======================================================================
-   !! History :  2.0   !  2006-12  (M. Vancoppenolle) Original code
-   !!            3.0   !  2008-02  (C. Talandier)  Surface module from icestp.F90
-   !!            9.0   !  2008-04  (G. Madec)  sltyle and lim_ctl routine
+   !!=====================================================================
+   !! History :  2.0  ! 2006-12  (M. Vancoppenolle) Original code
+   !!            3.0  ! 2008-02  (C. Talandier)  Surface module from icestp.F90
+   !!             -   ! 2008-04  (G. Madec)  sltyle and lim_ctl routine
+   !!            3.3  ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -20,9 +21,9 @@
    USE oce             ! ocean dynamics and tracers
    USE dom_oce         ! ocean space and time domain
-   USE lib_mpp
+   USE lib_mpp         ! MPP library
    USE par_ice         ! sea-ice parameters
-   USE ice
-   USE iceini
-   USE dom_ice
+   USE ice             ! LIM-3: ice variables
+   USE iceini          ! LIM-3: ice initialisation
+   USE dom_ice         ! LIM-3: ice domain
 
    USE sbc_oce         ! Surface boundary condition: ocean fields
@@ -30,5 +31,5 @@
    USE sbcblk_core     ! Surface boundary condition: CORE bulk
    USE sbcblk_clio     ! Surface boundary condition: CLIO bulk
-   USE albedo
+   USE albedo          ! ocean & ice albedo
 
    USE phycst          ! Define parameters for the routines
@@ -46,5 +47,5 @@
    USE limvar          ! Ice variables switch
 
-   USE lbclnk
+   USE lbclnk          ! lateral boundary condition - MPP link
    USE iom             ! I/O manager library
    USE in_out_manager  ! I/O manager
@@ -56,6 +57,4 @@
    PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90
    
-   CHARACTER(len=1) ::   cl_grid = 'C'     ! type of grid used in ice dynamics
-
    !! * Substitutions
 #  include "domzgr_substitute.h90"
@@ -64,10 +63,9 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
-   SUBROUTINE sbc_ice_lim( kt, kblk, kico )
+   SUBROUTINE sbc_ice_lim( kt, kblk )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE sbc_ice_lim  ***
@@ -91,5 +89,4 @@
       INTEGER, INTENT(in) ::   kt      ! ocean time step
       INTEGER, INTENT(in) ::   kblk    ! type of bulk (=3 CLIO, =4 CORE)
-      INTEGER, INTENT(in) ::   kico    ! ice-ocean stress treatment
       !!
       INTEGER  ::   jl                 ! loop index
@@ -142,5 +139,5 @@
                &                      qla_ice   , dqns_ice  , dqla_ice  ,               &
                &                      tprecip   , sprecip   ,                           &
-               &                      fr1_i0    , fr2_i0    , cl_grid, jpl  )
+               &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  )
             !         
          CASE( 4 )                                       ! CORE bulk formulation
@@ -149,5 +146,5 @@
                &                      qla_ice   , dqns_ice  , dqla_ice  ,               &
                &                      tprecip   , sprecip   ,                           &
-               &                      fr1_i0    , fr2_i0    , cl_grid, jpl  )
+               &                      fr1_i0    , fr2_i0    , cp_ice_msh, jpl  )
          END SELECT
 
@@ -160,34 +157,34 @@
          !                                           ! Store previous ice values
 !!gm : remark   old_...   should becomes ...b  as tn versus tb  
-         old_a_i(:,:,:)   = a_i(:,:,:)     ! ice area
-         old_e_i(:,:,:,:) = e_i(:,:,:,:)   ! ice thermal energy
-         old_v_i(:,:,:)   = v_i(:,:,:)     ! ice volume
-         old_v_s(:,:,:)   = v_s(:,:,:)     ! snow volume 
-         old_e_s(:,:,:,:) = e_s(:,:,:,:)   ! snow thermal energy
-         old_smv_i(:,:,:) = smv_i(:,:,:)   ! salt content
-         old_oa_i(:,:,:)  = oa_i(:,:,:)    ! areal age content
+         old_a_i  (:,:,:)   = a_i  (:,:,:)     ! ice area
+         old_e_i  (:,:,:,:) = e_i  (:,:,:,:)   ! ice thermal energy
+         old_v_i  (:,:,:)   = v_i  (:,:,:)     ! ice volume
+         old_v_s  (:,:,:)   = v_s  (:,:,:)     ! snow volume 
+         old_e_s  (:,:,:,:) = e_s  (:,:,:,:)   ! snow thermal energy
+         old_smv_i(:,:,:)   = smv_i(:,:,:)     ! salt content
+         old_oa_i (:,:,:)   = oa_i (:,:,:)     ! areal age content
 
          !                                           ! intialisation to zero    !!gm is it truly necessary ???
-         d_a_i_thd(:,:,:)   = 0.e0 ; d_a_i_trp(:,:,:)   = 0.e0
-         d_v_i_thd(:,:,:)   = 0.e0 ; d_v_i_trp(:,:,:)   = 0.e0
-         d_e_i_thd(:,:,:,:) = 0.e0 ; d_e_i_trp(:,:,:,:) = 0.e0
-         d_v_s_thd(:,:,:)   = 0.e0 ; d_v_s_trp(:,:,:)   = 0.e0
-         d_e_s_thd(:,:,:,:) = 0.e0 ; d_e_s_trp(:,:,:,:) = 0.e0
-         d_smv_i_thd(:,:,:) = 0.e0 ; d_smv_i_trp(:,:,:) = 0.e0
-         d_oa_i_thd(:,:,:)  = 0.e0 ; d_oa_i_trp(:,:,:)  = 0.e0
-         !
-         fseqv(:,:)     = 0.e0
-         fsbri(:,:)     = 0.e0     ; fsalt_res(:,:) = 0.e0
+         d_a_i_thd  (:,:,:)   = 0.e0   ;   d_a_i_trp  (:,:,:)   = 0.e0
+         d_v_i_thd  (:,:,:)   = 0.e0   ;   d_v_i_trp  (:,:,:)   = 0.e0
+         d_e_i_thd  (:,:,:,:) = 0.e0   ;   d_e_i_trp  (:,:,:,:) = 0.e0
+         d_v_s_thd  (:,:,:)   = 0.e0   ;   d_v_s_trp  (:,:,:)   = 0.e0
+         d_e_s_thd  (:,:,:,:) = 0.e0   ;   d_e_s_trp  (:,:,:,:) = 0.e0
+         d_smv_i_thd(:,:,:)   = 0.e0   ;   d_smv_i_trp(:,:,:)   = 0.e0
+         d_oa_i_thd (:,:,:)   = 0.e0   ;   d_oa_i_trp (:,:,:)   = 0.e0
+         !
+         fseqv    (:,:) = 0.e0
+         fsbri    (:,:) = 0.e0     ;   fsalt_res(:,:) = 0.e0
          fsalt_rpo(:,:) = 0.e0
-         fhmec(:,:)     = 0.e0     ; fhbri(:,:)     = 0.e0
-         fmmec(:,:)     = 0.e0     ; fheat_res(:,:) = 0.e0
-         fheat_rpo(:,:) = 0.e0     ; focea2D(:,:)   = 0.e0
-         fsup2D(:,:)    = 0.e0
+         fhmec    (:,:) = 0.e0     ;   fhbri    (:,:) = 0.e0
+         fmmec    (:,:) = 0.e0     ;   fheat_res(:,:) = 0.e0
+         fheat_rpo(:,:) = 0.e0     ;   focea2D  (:,:) = 0.e0
+         fsup2D   (:,:) = 0.e0
          ! 
-         diag_sni_gr(:,:) = 0.e0   ; diag_lat_gr(:,:) = 0.e0
-         diag_bot_gr(:,:) = 0.e0   ; diag_dyn_gr(:,:) = 0.e0
-         diag_bot_me(:,:) = 0.e0   ; diag_sur_me(:,:) = 0.e0
+         diag_sni_gr(:,:) = 0.e0   ;   diag_lat_gr(:,:) = 0.e0
+         diag_bot_gr(:,:) = 0.e0   ;   diag_dyn_gr(:,:) = 0.e0
+         diag_bot_me(:,:) = 0.e0   ;   diag_sur_me(:,:) = 0.e0
          ! dynamical invariants
-         delta_i(:,:) = 0.e0       ; divu_i (:,:) = 0.e0       ;    shear_i(:,:) = 0.e0
+         delta_i(:,:) = 0.e0       ;   divu_i(:,:) = 0.e0       ;   shear_i(:,:) = 0.e0
 
                           CALL lim_rst_opn( kt )     ! Open Ice restart file
@@ -196,5 +193,5 @@
          !
 #if ! defined key_c1d
-         ! Ice dynamics & transport (not in 1D case)
+                                                     ! Ice dynamics & transport (not in 1D case)
                           CALL lim_dyn( kt )              ! Ice dynamics    ( rheology/dynamics )
                           CALL lim_trp( kt )              ! Ice transport   ( Advection/diffusion )
@@ -204,5 +201,4 @@
                           CALL lim_itd_me                 ! Mechanical redistribution ! (ridging/rafting)
 #endif
-         !
          !                                           ! Ice thermodynamics 
                           CALL lim_var_glo2eqv            ! equivalent variables
@@ -216,5 +212,5 @@
                           CALL lim_itd_th( kt )           !  Remap ice categories, lateral accretion  !
          !
-         !                                           ! Global variables update |
+         !                                           ! Global variables update
                           CALL lim_var_agg( 1 )           ! requested by limupdate
                           CALL lim_update                 ! Global variables update
@@ -223,8 +219,5 @@
          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 2, ' - Final state - ' )   ! control print
          !
-         !                                           ! Fluxes of mass and heat to the ocean |
-                         CALL lim_sbc_flx( kt )           ! Ice/Ocean heat freshwater/salt fluxes
-         IF( ln_limdyn .AND. kico == 0 )   &              ! Ice/Ocean stresses (only in ice-dynamic case)
-            &            CALL lim_sbc_tau( kt, kico )     ! otherwise the atm.-ocean stresses are used everywhere
+                          CALL lim_sbc_flx( kt )     ! Update surface ocean mass, heat and salt fluxes
          !
          IF( ln_nicep )   CALL lim_prt_state( jiindx, jjindx, 3, ' - Final state lim_sbc - ' )   ! control print
@@ -239,12 +232,15 @@
          IF( ln_nicep )   CALL lim_ctl               ! alerts in case of model crash
          !
-      ENDIF                                          ! End sea-ice time step only
-
-      !                                              !--------------------------!
-      ! Ice/Ocean stresses (nn_ico_cpl=1 or 2 cases) !  at all ocean time step  !
-      !                                              !--------------------------!
-      IF( ln_limdyn .AND. kico /= 0 )   &
-         &                CALL lim_sbc_tau( kt, kico ) 
-!!gm   remark, in this case the ocean-ice stress is not saved in diag call above .....  find a solution!!!
+      ENDIF                                    ! End sea-ice time step only
+
+      !                                        !--------------------------!
+      !                                        !  at all ocean time step  !
+      !                                        !--------------------------!
+      !                                               
+      !                                              ! Update surface ocean stresses (only in ice-dynamic case)
+      !                                                   ! otherwise the atm.-ocean stresses are used everywhere
+      IF( ln_limdyn )     CALL lim_sbc_tau( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents
+      
+!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!!
       !
    END SUBROUTINE sbc_ice_lim
@@ -664,6 +660,6 @@
    !!----------------------------------------------------------------------
 CONTAINS
-   SUBROUTINE sbc_ice_lim ( kt, kblk, kico )     ! Dummy routine
-      WRITE(*,*) 'sbc_ice_lim: You should not have seen this print! error?', kt, kblk, kico
+   SUBROUTINE sbc_ice_lim ( kt, kblk )     ! Dummy routine
+      WRITE(*,*) 'sbc_ice_lim: You should not have seen this print! error?', kt, kblk
    END SUBROUTINE sbc_ice_lim
 #endif
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 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim_2.F90	(revision 2370)
@@ -2,7 +2,6 @@
    !!======================================================================
    !!                       ***  MODULE  sbcice_lim_2  ***
-   !! Surface module :  update surface ocean boundary condition over ice
-   !!                   covered area using LIM sea-ice model
-   !! Sea-Ice model  :  LIM 2.0 Sea ice model time-stepping
+   !! Surface module :  update surface ocean boundary condition over ice covered area using LIM sea-ice model
+   !! Sea-Ice model  :  LIM-2 Sea ice model time-stepping
    !!======================================================================
    !! History :  1.0   !  06-2006  (G. Madec)  from icestp_2.F90
@@ -12,12 +11,10 @@
 #if defined key_lim2
    !!----------------------------------------------------------------------
-   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   sbc_ice_lim_2  : sea-ice model time-stepping and
-   !!                    update ocean sbc over ice-covered area
-   !!----------------------------------------------------------------------
-   USE oce             ! ocean dynamics and tracers
-   USE dom_oce         ! ocean space and time domain
-   USE lib_mpp
+   !!   'key_lim2' :                                    LIM-2 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   sbc_ice_lim_2   : sea-ice model time-stepping and update ocean sbc over ice-covered area
+   !!----------------------------------------------------------------------
+   USE oce              ! ocean dynamics and tracers
+   USE dom_oce          ! ocean space and time domain
    USE ice_2
    USE par_ice_2
@@ -25,26 +22,27 @@
    USE dom_ice_2
 
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE sbc_ice         ! Surface boundary condition: ice   fields
-   USE sbcblk_core     ! Surface boundary condition: CORE bulk
-   USE sbcblk_clio     ! Surface boundary condition: CLIO bulk
-   USE sbccpl          ! Surface boundary condition: coupled interface
+   USE sbc_oce          ! Surface boundary condition: ocean fields
+   USE sbc_ice          ! Surface boundary condition: ice   fields
+   USE sbcblk_core      ! Surface boundary condition: CORE bulk
+   USE sbcblk_clio      ! Surface boundary condition: CLIO bulk
+   USE sbccpl           ! Surface boundary condition: coupled interface
    USE albedo
 
-   USE phycst          ! Define parameters for the routines
-   USE eosbn2          ! equation of state
+   USE phycst           ! Define parameters for the routines
+   USE eosbn2           ! equation of state
    USE limdyn_2
    USE limtrp_2
    USE limdmp_2
    USE limthd_2
-   USE limsbc_2        ! sea surface boundary condition
+   USE limsbc_2         ! sea surface boundary condition
    USE limdia_2
    USE limwri_2
    USE limrst_2
 
-   USE lbclnk
-   USE iom             ! I/O manager library
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
+   USE lbclnk           ! lateral boundary condition - MPP link
+   USE lib_mpp          ! MPP library
+   USE iom              ! I/O manager library
+   USE in_out_manager   ! I/O manager
+   USE prtctl           ! Print control
 
    IMPLICIT NONE
@@ -59,7 +57,6 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -97,21 +94,30 @@
          IF(lwp) WRITE(numout,*) 'sbc_ice_lim_2 : update ocean surface boudary condition' 
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   via Louvain la Neuve Ice Model (LIM) time stepping'
-
+         !
          CALL ice_init_2
-
       ENDIF
 
-      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-         !
+      !                                        !----------------------!
+      IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN     !  Ice time-step only  !
+         !                                     !----------------------!
+         !  Bulk Formulea !
+         !----------------!
          ! ... mean surface ocean current at ice dynamics point
-         !     B-grid dynamics :  I-point 
-         DO jj = 2, jpj
-            DO ji = 2, jpi   ! B grid : no vector opt.
-               u_oce(ji,jj) = 0.5 * ( ssu_m(ji-1,jj  ) + ssu_m(ji-1,jj-1) ) * tmu(ji,jj)
-               v_oce(ji,jj) = 0.5 * ( ssv_m(ji  ,jj-1) + ssv_m(ji-1,jj-1) ) * tmu(ji,jj)
+         SELECT CASE( cp_ice_msh )
+         CASE( 'I' )                  !== B-grid ice dynamics :   I-point (i.e. F-point with sea-ice indexation)
+            DO jj = 2, jpj
+               DO ji = 2, jpi   ! NO vector opt. possible
+                  u_oce(ji,jj) = 0.5_wp * ( ssu_m(ji-1,jj  ) + ssu_m(ji-1,jj-1) ) * tmu(ji,jj)
+                  v_oce(ji,jj) = 0.5_wp * ( ssv_m(ji  ,jj-1) + ssv_m(ji-1,jj-1) ) * tmu(ji,jj)
+               END DO
             END DO
-         END DO
-         CALL lbc_lnk( u_oce, 'I', -1. )   ! I-point (i.e. F-point with ice indices)
-         CALL lbc_lnk( v_oce, 'I', -1. )   ! I-point (i.e. F-point with ice indices)
+            CALL lbc_lnk( u_oce, 'I', -1. )   ! I-point (i.e. F-point with ice indices)
+            CALL lbc_lnk( v_oce, 'I', -1. )   ! I-point (i.e. F-point with ice indices)
+            !
+         CASE( 'C' )                  !== C-grid ice dynamics :   U & V-points (same as ocean)
+            u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point
+            v_oce(:,:) = ssv_m(:,:)
+            !
+         END SELECT
 
          ! ... masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
@@ -142,5 +148,5 @@
                &                      qla_ice    , dqns_ice   , dqla_ice   ,            &
                &                      tprecip    , sprecip    ,                         &
-               &                      fr1_i0     , fr2_i0     , cl_grid    , jpl  )
+               &                      fr1_i0     , fr2_i0     , cp_ice_msh , jpl  )
 
          CASE( 4 )           ! CORE bulk formulation
@@ -149,5 +155,5 @@
                &                      qla_ice    , dqns_ice   , dqla_ice   ,            &
                &                      tprecip    , sprecip    ,                         &
-               &                      fr1_i0     , fr2_i0     , cl_grid    , jpl  )
+               &                      fr1_i0     , fr2_i0     , cp_ice_msh , jpl  )
          CASE( 5 )           ! Coupled formulation : atmosphere-ice stress only (fluxes provided after ice dynamics)
             CALL sbc_cpl_ice_tau( utau_ice , vtau_ice )
@@ -170,31 +176,40 @@
          !  Ice model step  !
          ! ---------------- !
-         numit = numit + nn_fsbc                                             !Ice model time 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)
-                                        CALL lim_dyn_2      ( kt )           ! Ice dynamics    ( rheology/dynamics )
-                                        CALL lim_trp_2      ( kt )           ! Ice transport   ( Advection/diffusion )
-            IF( ln_limdmp )             CALL lim_dmp_2      ( kt )           ! Ice damping 
+                                                       ! Ice dynamics & transport (except in 1D case)
+                           CALL lim_dyn_2      ( kt )      ! Ice dynamics    ( rheology/dynamics )
+                           CALL lim_trp_2      ( kt )      ! Ice transport   ( Advection/diffusion )
+         IF( ln_limdmp )   CALL lim_dmp_2      ( kt )      ! Ice damping 
 #endif
 #if defined key_coupled
-         IF( ksbc == 5    )             CALL sbc_cpl_ice_flx( reshape( frld, (/jpi,jpj,1/) ),        &
-      &                                                       qns_tot, qns_ice, qsr_tot , qsr_ice,   &
-      &                                                       emp_tot, emp_ice, dqns_ice, sprecip,   &
+         !                                             ! Ice surface fluxes in coupled mode 
+         IF( ksbc == 5 )   CALL sbc_cpl_ice_flx( reshape( frld, (/jpi,jpj,1/) ),                 &
+      &                                                   qns_tot, qns_ice, qsr_tot , qsr_ice,   &
+      &                                                   emp_tot, emp_ice, dqns_ice, sprecip,   &
       !                                      optional arguments, used only in 'mixed oce-ice' case
-      &                                                       palbi = zalb_ice_cs, psst = sst_m, pist = sist )
+      &                                                   palbi = zalb_ice_cs, psst = sst_m, pist = sist )
 #endif
-                                        CALL lim_thd_2      ( kt )      ! Ice thermodynamics 
-                                        CALL lim_sbc_2      ( kt )      ! Ice/Ocean Mass & Heat fluxes 
+                           CALL lim_thd_2      ( kt )      ! Ice thermodynamics 
+                           CALL lim_sbc_flx_2  ( kt )      ! update surface ocean mass, heat & salt fluxes 
 
          IF( ( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR. ntmoy == 1 ) .AND. .NOT. lk_mpp )   &
-            &                           CALL lim_dia_2      ( kt )      ! Ice Diagnostics
+            &              CALL lim_dia_2      ( kt )      ! Ice Diagnostics
 # if ! defined key_iomput
-                                        CALL lim_wri_2      ( kt )      ! Ice outputs
+                           CALL lim_wri_2      ( kt )      ! Ice outputs
 # endif
-         IF( lrst_ice )                 CALL lim_rst_write_2( kt )      ! Ice restart file
+         IF( lrst_ice  )   CALL lim_rst_write_2( kt )      ! Ice restart file
          !
-      ENDIF
+      ENDIF                                    ! End sea-ice time step only
+      !
+      !                                        !--------------------------!
+      !                                        !  at all ocean time step  !
+      !                                        !--------------------------!
+      !                                               
+      !                                              ! Update surface ocean stresses (only in ice-dynamic case)
+      !                                                   ! otherwise the atm.-ocean stresses are used everywhere
+      IF( ln_limdyn    )   CALL lim_sbc_tau_2( kt, ub(:,:,1), vb(:,:,1) )  ! using before instantaneous surf. currents
       !
    END SUBROUTINE sbc_ice_lim_2
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90	(revision 2370)
@@ -4,9 +4,10 @@
    !! Surface module :  provide to the ocean its surface boundary condition
    !!======================================================================
-   !! History :  3.0  !  2006-07  (G. Madec)  Original code
-   !!            3.1  !  2008-08  (S. Masson, A. Caubel, E. Maisonnave, G. Madec) coupled interface
-   !!            3.3  !  2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
-   !!            3.3  !  2010-10  (S. Masson)  add diurnal cycle
-   !!            3.3   !  09-2010  (D. Storkey) add ice boundary conditions (BDY)
+   !! History :  3.0  ! 2006-07  (G. Madec)  Original code
+   !!            3.1  ! 2008-08  (S. Masson, A. Caubel, E. Maisonnave, G. Madec) coupled interface
+   !!            3.3  ! 2010-04  (M. Leclair, G. Madec)  Forcing averaged over 2 time steps
+   !!            3.3  ! 2010-10  (S. Masson)  add diurnal cycle
+   !!            3.3  ! 2010-09  (D. Storkey) add ice boundary conditions (BDY)
+   !!             -   ! 2010-11  (G. Madec) ice-ocean stress always computed at each ocean time-step
    !!----------------------------------------------------------------------
 
@@ -74,5 +75,5 @@
       !!
       NAMELIST/namsbc/ nn_fsbc, ln_ana  , ln_flx, ln_blk_clio, ln_blk_core, ln_cpl    ,   &
-         &             nn_ice , ln_dm2dc, ln_rnf, ln_ssr     , nn_fwb     , nn_ico_cpl
+         &             nn_ice , ln_dm2dc, ln_rnf, ln_ssr     , nn_fwb
       !!----------------------------------------------------------------------
 
@@ -107,5 +108,4 @@
          WRITE(numout,*) '           Misc. options of sbc : '
          WRITE(numout,*) '              ice management in the sbc (=0/1/2/3)       nn_ice      = ', nn_ice 
-         WRITE(numout,*) '              ice-ocean stress computation (=0/1/2)      nn_ico_cpl  = ', nn_ico_cpl
          WRITE(numout,*) '              daily mean to diurnal cycle qsr            ln_dm2dc    = ', ln_dm2dc 
          WRITE(numout,*) '              runoff / runoff mouths                     ln_rnf      = ', ln_rnf
@@ -238,10 +238,10 @@
       
       SELECT CASE( nn_ice )                                     ! Update heat and freshwater fluxes over sea-ice areas
-      CASE(  1 )   ;       CALL sbc_ice_if   ( kt )                   ! Ice-cover climatology ("Ice-if" model)
+      CASE(  1 )   ;       CALL sbc_ice_if   ( kt )                  ! Ice-cover climatology ("Ice-if" model)
          !                                                      
-      CASE(  2 )   ;       CALL sbc_ice_lim_2( kt, nsbc )             ! LIM 2.0 ice model
-                           IF( lk_bdy ) CALL bdy_ice_frs( kt )
+      CASE(  2 )   ;       CALL sbc_ice_lim_2( kt, nsbc )            ! LIM-2 ice model
+         IF( lk_bdy )      CALL bdy_ice_frs  ( kt )                  ! BDY boundary condition
          !                                                     
-      CASE(  3 )   ;       CALL sbc_ice_lim  ( kt, nsbc, nn_ico_cpl)  ! LIM 3.0 ice model
+      CASE(  3 )   ;       CALL sbc_ice_lim  ( kt, nsbc )            ! LIM-3 ice model
       END SELECT                                              
 
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/SBC/sbcrnf.F90	(revision 2370)
@@ -4,8 +4,9 @@
    !! Ocean forcing:  river runoff
    !!=====================================================================
-   !! History :  OPA  !  2000-11  (R. Hordoir, E. Durand)  NetCDF FORMAT
-   !!   NEMO     1.0  !  2002-09  (G. Madec)  F90: Free form and module
-   !!            3.0  !  2006-07  (G. Madec)  Surface module 
-   !!            3.2  !  2009-04  (B. Lemaire)  Introduce iom_put
+   !! History :  OPA  ! 2000-11  (R. Hordoir, E. Durand)  NetCDF FORMAT
+   !!   NEMO     1.0  ! 2002-09  (G. Madec)  F90: Free form and module
+   !!            3.0  ! 2006-07  (G. Madec)  Surface module 
+   !!            3.2  ! 2009-04  (B. Lemaire)  Introduce iom_put
+   !!            3.3  ! 2010-10  (R. Furner, G. Madec) runoff distributed over ocean levels
    !!----------------------------------------------------------------------
 
@@ -49,14 +50,13 @@
    REAL(wp), PUBLIC, DIMENSION(jpk)     ::   rnfmsk_z     !: river mouth mask (vert.)
 
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnf       !: structure of input river runoff (file information, fields read)
-
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_s_rnf     !: structure of input river runoff salinity (file information, fields read)  
-   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     !: structure of input river runoff temperature (file information, fields read)  
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_rnf       !: structure: river runoff (file information, fields read)
+
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_s_rnf     !: structure: river runoff salinity (file information, fields read)  
+   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_t_rnf     !: structure: river runoff temperature (file information, fields read)  
  
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   h_rnf        !: depth of runoff in m
    INTEGER,  PUBLIC, DIMENSION(jpi,jpj) ::   nk_rnf       !: depth of runoff in model levels
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpts) :: rnf_tsc_b, rnf_tsc  !: before and now 
-   !                                                                 !: temp. & sal. content of river runoffs   [K.m/s & PSU.m/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpts) :: rnf_tsc_b, rnf_tsc  !: before and now T & S contents of runoffs  [K.m/s & PSU.m/s]
 
    !! * Substitutions  
@@ -65,7 +65,6 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
-
 CONTAINS
 
@@ -171,7 +170,7 @@
          CALL iom_rstput( kt, nitrst, numrow, 'rnf_sc_b', rnf_tsc(:,:,jp_sal) )
       ENDIF
-
       !
    END SUBROUTINE sbc_rnf
+
 
    SUBROUTINE sbc_rnf_div( phdivn )
Index: /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/step.F90
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/NEMO/OPA_SRC/step.F90	(revision 2370)
@@ -26,14 +26,14 @@
 
    !!----------------------------------------------------------------------
-   !!   stp            : OPA system time-stepping
+   !!   stp             : OPA system time-stepping
    !!----------------------------------------------------------------------
    USE step_oce         ! time stepping definition modules 
 #if defined key_top
-   USE trcstp          ! passive tracer time-stepping      (trc_stp routine)
-#endif
-   USE asminc          ! assimilation increments    (tra_asm_inc, dyn_asm_inc routines)
-   USE stpctl          ! time stepping control            (stp_ctl routine)
-   USE restart         ! ocean restart                    (rst_wri routine)
-   USE prtctl          ! Print control                    (prt_ctl routine)
+   USE trcstp           ! passive tracer time-stepping      (trc_stp routine)
+#endif
+   USE asminc           ! assimilation increments    (tra_asm_inc, dyn_asm_inc routines)
+   USE stpctl           ! time stepping control            (stp_ctl routine)
+   USE restart          ! ocean restart                    (rst_wri routine)
+   USE prtctl           ! Print control                    (prt_ctl routine)
 
 #if defined key_agrif
@@ -44,5 +44,5 @@
    PRIVATE
 
-   PUBLIC   stp        ! called by opa.F90
+   PUBLIC   stp   ! called by opa.F90
 
    !! * Substitutions
@@ -52,7 +52,6 @@
    !! NEMO/OPA 3.3 , NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
 
@@ -149,10 +148,9 @@
          IF( ln_zps )    CALL zps_hde( kstp, jpts, tsb, gtsu, gtsv,  &    ! Partial steps: before horizontal gradient
             &                                      rhd, gru , grv  )      ! of t, s, rd at the last ocean level
-         IF ( .NOT. ln_traldf_grif ) THEN
-                         CALL ldf_slp( kstp, rhd, rn2b )        ! before slope of the lateral mixing
+         IF( ln_traldf_grif ) THEN                           ! before slope for Griffies operator
+                         CALL ldf_slp_grif( kstp )
          ELSE
-                         CALL ldf_slp_grif( kstp )
-                         IF ( ln_dynldf_bilap .OR. ln_dynldf_iso ) CALL ldf_slp( kstp, rhd, rn2b )
-      ENDIF
+                         CALL ldf_slp( kstp, rhd, rn2b )     ! before slope for Madec operator
+         ENDIF
       ENDIF
 #if defined key_traldf_c2d
Index: /branches/nemo_v3_3_beta/NEMOGCM/TOOLS/COMPILE/cfg.txt
===================================================================
--- /branches/nemo_v3_3_beta/NEMOGCM/TOOLS/COMPILE/cfg.txt	(revision 2369)
+++ /branches/nemo_v3_3_beta/NEMOGCM/TOOLS/COMPILE/cfg.txt	(revision 2370)
@@ -4,3 +4,4 @@
 ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC
 POMME OPA_SRC NST_SRC
+ORCA2_LIM3 OPA_SRC LIM_SRC_3
 ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC C1D_SRC
