Index: trunk/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- trunk/NEMO/LIM_SRC_3/limsbc.F90	(revision 1525)
+++ trunk/NEMO/LIM_SRC_3/limsbc.F90	(revision 1526)
@@ -40,11 +40,12 @@
    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) ::   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), DIMENSION(jpi,jpj) ::   ssu_mb  , ssv_mb     !: before mean ocean surface currents           [m/s]
+
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !! NEMO/LIM  3.0 ,  UCL-LOCEAN-IPSL (2008) 
+   !! NEMO/LIM  3.2 ,  UCL-LOCEAN-IPSL (2009) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -77,16 +78,14 @@
       !                   !                               =2  combination of 0 and 1 cases
       !!
-      INTEGER  ::   ji, jj           ! dummy loop indices
-      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points
-      REAL(wp) ::   zat_u, zu_ico, zutaui, zu_u
-      REAL(wp) ::   zat_v, zv_ico, zvtaui, zv_v, zsang
-      REAL(wp) ::   zu_ij, zu_im1j, zv_ij, zv_ijm1
-
+      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                                                 !    -         -
+      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
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice
-      !!---------------------------------------------------------------------
+     !!---------------------------------------------------------------------
 
       IF( kt == nit000 ) THEN
@@ -100,5 +99,5 @@
       CASE( 0 )                                   !  LIM 3 old stress computation  !  (at ice timestep only)
          !                                        !--------------------------------! 
-         DO jj = 2, jpjm1                      ! ... modulus of the ice-ocean velocity
+         DO jj = 2, jpjm1                             !* modulus of the ice-ocean velocity
             DO ji = fs_2, fs_jpim1
                zu_ij   = u_ice(ji  ,jj) - ssu_m(ji  ,jj)               ! (i  ,j)
@@ -110,42 +109,34 @@
             END DO
          END DO
-         CALL lbc_lnk( tmod_io, 'T', 1. )
-         ! ... ice stress over ocean with a ice-ocean rotation angle
+         CALL lbc_lnk( tmod_io, 'T', 1. )   ! lateral boundary condition
+         !
+         !                                             !* ice stress over ocean with a ice-ocean rotation angle
          DO jj = 1, jpjm1
-            ! ... change the sinus angle sign in the south hemisphere
-            zsang  = SIGN(1.e0, gphif(1,jj) ) * sangvg
+            zsang  = SIGN( 1.e0, gphif(1,jj) ) * sangvg         ! change the sinus angle sign in the south hemisphere
             DO ji = 1, fs_jpim1
-               ! ... ice velocity relative to the ocean
-               zu_u = u_ice(ji,jj) - u_oce(ji,jj)
+               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
+               !                                                ! 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 ) 
-               ! 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)
-               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
+               !                                                ! 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.
-               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
-               zfrldu = 0.5 * ( ato_i(ji,jj) + ato_i(ji+1,jj  ) )
+               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
+               !                                                    ! 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
-
-         ! boundary condition on the stress (utau,vtau)
-         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )
+            END DO
+         END DO
+         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
+
          !
          !                                        !--------------------------------!
@@ -153,7 +144,7 @@
          !                                        !--------------------------------! 
          IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-            utau_oce(:,:) = utau(:,:)             ! ... save the air-ocean stresses at ice time-step
+            utau_oce(:,:) = utau(:,:)                 !* save the air-ocean stresses at ice time-step
             vtau_oce(:,:) = vtau(:,:)
-            DO jj = 2, jpjm1                      ! ... modulus of the ice-ocean velocity
+            DO jj = 2, jpjm1                          !* modulus of the ice-ocean velocity
                DO ji = fs_2, fs_jpim1
                   zu_ij   = u_ice(ji  ,jj) - ssu_m(ji  ,jj)               ! (i  ,j)
@@ -165,28 +156,26 @@
                END DO
             END DO
-         CALL lbc_lnk( tmod_io, 'T', 1. )
+            CALL lbc_lnk( tmod_io, 'T', 1. )          ! lateral boundary condition
          ENDIF
-         ! ... ice stress over ocean with a ice-ocean rotation angle
-         DO jj = 2, jpjm1
-            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
+         !
+        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
-               ! computation of wind stress over ocean in X and Y direction
-               zat_u  = ( at_i(ji,jj) + at_i(ji+1,jj) ) * 0.5      ! ice area at u and V-points
+               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_u   = u_ice(ji,jj) - ub(ji,jj,1)                  ! u ice-ocean velocity at U-point
-               zv_v   = v_ice(ji,jj) - vb(ji,jj,1)                  ! v ice-ocean velocity at V-point
-               !
+               !                                                ! (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 ) 
-
-               utau(ji,jj) = ( 1.- zat_u ) * utau_oce(ji,jj) + zat_u * zutaui    ! stress at the ocean surface
+               !                                                   ! 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
-               !                                                    ! u ice-ocean velocity at V-point
-            END DO
-         END DO
-         ! boundary condition on the stress (utau,vtau)
-         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )
+            END DO
+         END DO
+         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
+
          !
          !                                        !--------------------------------!
@@ -194,9 +183,9 @@
          !                                        !--------------------------------! 
          IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
-            utau_oce(:,:) = utau (:,:)            ! ... save the air-ocean stresses at ice time-step
+            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
+            ssu_mb  (:,:) = ssu_m(:,:)               !* save the ice-ocean velocity at ice time-step
             ssv_mb  (:,:) = ssv_m(:,:)
-            DO jj = 2, jpjm1                      ! ... modulus of the ice-ocean velocity
+            DO jj = 2, jpjm1                         !* modulus of the ice-ocean velocity
                DO ji = fs_2, fs_jpim1
                   zu_ij   = u_ice(ji  ,jj) - ssu_m(ji  ,jj)               ! (i  ,j)
@@ -208,16 +197,13 @@
                END DO
             END DO
-         CALL lbc_lnk( tmod_io, 'T', 1. )
+            CALL lbc_lnk( tmod_io, 'T', 1. )
          ENDIF
-         ! ... ice stress over ocean with a ice-ocean rotation angle
-         DO jj = 2, jpjm1
+         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
-               ! computation of wind stress over ocean in X and Y direction
                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 
-
-               !!gm bug mixing U and V points value below     ====>>> to be corrected
-               zu_ico = u_ice(ji,jj) - 0.5 * ( ub(ji,jj,1) - ssu_mb(ji,jj) )   ! ice-oce velocity using un and ssu_mb
+               !
+               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
@@ -230,8 +216,8 @@
             END DO
          END DO
-         ! boundary condition on the stress (utau,vtau)
-         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )
+         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )   ! lateral boundary condition
          !
       END SELECT
+
       IF(ln_ctl)   CALL prt_ctl( tab2d_1=utau, clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
          &                       tab2d_2=vtau, clinfo2=' vtau    : '        , mask2=vmask )
@@ -389,5 +375,5 @@
       !------------------------------------------!
 
-      !!gm   optimisation: this loop have to be merged with the previous one
+!!gm   optimisation: this loop have to be merged with the previous one
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -475,5 +461,4 @@
    END SUBROUTINE lim_sbc_flx
 
-
 #else
    !!----------------------------------------------------------------------
