Index: trunk/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- trunk/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 2524)
+++ trunk/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90	(revision 2528)
@@ -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
@@ -21,9 +23,8 @@
    USE sbc_oce          ! Surface boundary condition: ocean fields
    USE phycst           ! physical constants
+   USE albedo           ! albedo parameters
    USE ice              ! LIM sea-ice variables
-
    USE lbclnk           ! ocean lateral boundary condition
    USE in_out_manager   ! I/O manager
-   USE albedo           ! albedo parameters
    USE prtctl           ! Print control
 
@@ -34,181 +35,22 @@
    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
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/LIM  3.2 ,  UCL-LOCEAN-IPSL (2009) 
+   !! NEMO/LIM3 3.3 , UCL - NEMO Consortium (2010)
    !! $Id$
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
 CONTAINS
-
-   SUBROUTINE lim_sbc_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 )
@@ -234,5 +76,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
@@ -253,4 +95,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
 
@@ -297,6 +153,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
@@ -339,11 +195,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)
@@ -374,18 +230,16 @@
 
             !  computing freshwater exchanges at the ice/ocean interface
-            zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj) )  &   !  evaporation over oceanic fraction
-               &   + tprecip(ji,jj) *         at_i(ji,jj)    &   !  total precipitation
-               ! old fashioned way               
-               !              &   - 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 
-               ! new contribution from snow falling when ridging
-               &   + fmmec(ji,jj)
+            zpme = - emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction
+               &   + tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean
+               &   - sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice
+               &   - rdmsnif(ji,jj) * r1_rdtice                       &   ! freshwaterflux due to snow melting 
+               &   + fmmec(ji,jj)                                         ! snow falling when ridging
+
 
             !  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
@@ -445,4 +299,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
    !!----------------------------------------------------------------------
