Index: /trunk/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- /trunk/NEMO/LIM_SRC_3/limdyn.F90	(revision 912)
+++ /trunk/NEMO/LIM_SRC_3/limdyn.F90	(revision 913)
@@ -18,4 +18,5 @@
    USE ice
    USE par_ice
+   USE sbc_oce         ! Surface boundary condition: ocean fields
    USE sbc_ice         ! Surface boundary condition: ice fields
    USE ice_oce
@@ -66,19 +67,10 @@
       !!------------------------------------------------------------------------------------
       !! * Local variables
-      INTEGER ::   ji, jj, jl, ja    ! dummy loop indices
-      INTEGER :: i_j1, i_jpj         ! Starting/ending j-indices for rheology
-
-      REAL(wp) ::   &
-         ztairx, ztairy,          &  ! tempory scalars
-         zsang , zmod,            &
-         ztglx , ztgly ,          &
-         zt11, zt12, zt21, zt22 , &
-         zustm,                   &
-         zsfrldmx2, zsfrldmy2,    &
-         zu_ice, zv_ice, ztair2
-
-      REAL(wp),DIMENSION(jpj) ::   &
-           zind,                     &  ! i-averaged indicator of sea-ice
-           zmsk                         ! i-averaged of tmask
+      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), DIMENSION(jpj)     ::   zind           ! i-averaged indicator of sea-ice
+      REAL(wp), DIMENSION(jpj)     ::   zmsk           ! i-averaged of tmask
+      REAL(wp), DIMENSION(jpi,jpj) ::   zu_io, zv_io   ! ice-ocean velocity
       !!---------------------------------------------------------------------
 
@@ -159,43 +151,27 @@
 
          ! computation of friction velocity
+         ! --------------------------------
+         ! ice-ocean velocity at U & V-points (u_ice vi_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points)
+         zu_io(:,:) = u_ice(:,:) - ssu_m(:,:)
+         zv_io(:,:) = v_ice(:,:) - ssv_m(:,:)
+         ! frictional velocity at T-point
+         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)
+            END DO
+         END DO
+         !
+      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
+         !
+         zcoef = SQRT( 0.5 ) / rau0
          DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1
-
-               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
-               zt11  = rhoco * zu_ice * zu_ice
-
-               zu_ice   = u_ice(ji-1,jj) - u_oce(ji-1,jj)
-               zt12  = rhoco * zu_ice * zu_ice
-
-               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
-               zt21  = rhoco * zv_ice * zv_ice
-
-               zv_ice   = v_ice(ji,jj-1) - v_oce(ji,jj-1)
-               zt22  = rhoco * zv_ice * zv_ice
-               ztair2 = ( ( utaui_ice(ji,jj) + utaui_ice(ji-1,jj) ) / 2. )**2 + &
-                        ( ( vtaui_ice(ji,jj) + vtaui_ice(ji,jj-1) ) / 2. )**2
-
-               ! should not be weighted
-               zustm =  ( at_i(ji,jj)       ) * 0.5 * ( zt11 + zt12 + zt21 + zt22 )        &
-                  &  +  ( 1.0 - at_i(ji,jj) ) * SQRT( ztair2 )
-
-               ust2s(ji,jj) = ( zustm / rau0 ) * tms(ji,jj)
-
+            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) )
             END DO
          END DO
-
-       ELSE      ! If no ice dynamics  
-                    
-               ! virer ca (key_lim_cp1)
-          DO jj = 2, jpjm1
-             DO ji = fs_2, fs_jpim1
-                ztair2 = ( ( utaui_ice(ji,jj) + utaui_ice(ji-1,jj) ) / 2. )**2 + &
-                         ( ( vtaui_ice(ji,jj) + vtaui_ice(ji,jj-1) ) / 2. )**2
-                zustm        = SQRT( ztair2  )
-
-                ust2s(ji,jj) = ( zustm / rau0 ) * tms(ji,jj)
-            END DO
-         END DO
-
+         !
       ENDIF
 
Index: /trunk/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- /trunk/NEMO/LIM_SRC_3/limsbc.F90	(revision 912)
+++ /trunk/NEMO/LIM_SRC_3/limsbc.F90	(revision 913)
@@ -77,16 +77,9 @@
       INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
       INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
-      REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux
       REAL(wp) ::   zinda            ! switch for testing the values of ice concentration
       REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface
       REAL(wp) ::   zpme             ! freshwater exchanges at the ice/ocean interface
       REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points
-      REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points
       REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity
-      REAL(wp) ::   ztairx, ztairy   
-      REAL(wp) ::   zsfrldmx2, zsfrldmy2
-      REAL(wp) ::   zu_ice, zv_ice
-      REAL(wp) ::   ztglx , ztgly
-      REAL(wp) ::   ztaux , ztauy
       
 #if defined key_coupled    
@@ -94,5 +87,5 @@
       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalbp    ! albedo of ice under clear sky
 #endif
-      REAL(wp) ::   zsang, zmod, zfm
+      REAL(wp) ::   zsang, zmod
       REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice
       REAL(wp), DIMENSION(jpi,jpj) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes
@@ -277,62 +270,30 @@
 
          ! ... ice stress over ocean with a ice-ocean rotation angle
-         DO jj = 2, jpjm1
-            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
-
-            DO ji = fs_2, fs_jpim1
-               ! computation of wind stress over ocean in X and Y direction
-               ztairx =  ( 2.0 - at_i(ji,jj) - at_i(ji+1,jj) ) * utau(ji,jj)
-               ztairy =  ( 2.0 - at_i(ji,jj) - at_i(ji,jj+1) ) * vtau(ji,jj)
-
-               zsfrldmx2 = at_i(ji,jj) + at_i(ji+1,jj)
-               zsfrldmy2 = at_i(ji,jj) + at_i(ji,jj+1)
-
-               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
-               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
-               zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 
-
+         DO jj = 1, jpj
+            ! ... change the sinus angle sign in the south hemisphere
+            zsang  = SIGN(1.e0, gphif(1,jj) ) * sangvg
+            DO ji = 1, jpi
+               ! ... ice velocity relative to the ocean
+               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 ) 
                ! quadratic drag formulation
-               ztglx   = zsfrldmx2 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 
-               ztgly   = zsfrldmy2 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 
-!
-!              ! IMPORTANT
-!              ! these lines are bound to prevent numerical oscillations
-!              ! in the ice-ocean stress
-!              ! They are physically ill-based. There is a cleaner solution
-!              ! to try (remember discussion in Paris Gurvan)
-!
-               ztglx   = ztglx * exp( - zmod / 0.5 )
-               ztgly   = ztglx * exp( - zmod / 0.5 )
-
-               tio_u(ji,jj) = ( ztairx + 1.0 * ztglx ) / 2.
-               tio_v(ji,jj) = ( ztairy + 1.0 * ztgly ) / 2.
+               ztio_u(ji,jj) = rhoco * zmod * ( cangvg * zu_io - zsang * zv_io ) 
+               ztio_v(ji,jj) = rhoco * zmod * ( cangvg * zv_io + zsang * zu_io ) 
+               !
             END DO
          END DO
 
-         DO jj = 1, jpjm1
-            DO ji = 1, fs_jpim1   ! vertor opt.
-               zfrldu = MAX( freezn(ji,jj), freezn(ji,jj+1) )   ! ice/ocean indicator at U- and V-points
-               zfrldv = MAX( freezn(ji,jj), freezn(ji+1,jj) )
-               ztaux = tio_u(ji,jj)
-               ztauy = tio_v(ji,jj)
-               utau(ji,jj) = (1.-zfrldu) * utau(ji,jj) + zfrldu * ztaux    ! stress at the ocean surface
-               vtau(ji,jj) = (1.-zfrldv) * vtau(ji,jj) + zfrldv * ztauy
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1   ! vertor opt.
+               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
+               zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj  ) )
+               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
-
-!new fashion         DO jj = 2, jpjm1
-!new fashion            DO ji = fs_2, fs_jpim1   ! vertor opt.
-!new fashion               ! ... ice-cover wheighted ice-ocean stress at U and V-points  (from I-point values)
-!new fashion               zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
-!new fashion               zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
-!new fashion               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
-!new fashion               zfrldu = 0.5 * ( frld (ji,jj) + frld (ji+1,jj  ) )
-!new fashion               zfrldv = 0.5 * ( frld (ji,jj) + frld (ji  ,jj+1) )
-!new fashion               ! update surface ocean stress
-!new fashion               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * zutau
-!new fashion               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * zvtau
-!new fashion               !
-!new fashion            END DO
-!new fashion         END DO
 
          ! boundary condition on the stress (utau,vtau)
