Index: trunk/NEMO/LIM_SRC_3/limsbc.F90
===================================================================
--- trunk/NEMO/LIM_SRC_3/limsbc.F90	(revision 917)
+++ trunk/NEMO/LIM_SRC_3/limsbc.F90	(revision 918)
@@ -4,7 +4,7 @@
    !!           computation of the flux at the sea ice/ocean interface
    !!======================================================================
-   !! History : 00-01 (H. Goosse) Original code
-   !!           02-07 (C. Ethe, G. Madec) re-writing F90
-   !!           06-07 (G. Madec) surface module
+   !! History :   -   !  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
    !!----------------------------------------------------------------------
 #if defined key_lim3
@@ -12,7 +12,7 @@
    !!   'key_lim3'                                    LIM 3.0 sea-ice model
    !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
    !!   lim_sbc  : flux at the ice / ocean interface
    !!----------------------------------------------------------------------
+   USE oce
    USE par_oce          ! ocean parameters
    USE par_ice          ! ice parameters
@@ -35,5 +35,6 @@
    PRIVATE
 
-   PUBLIC lim_sbc       ! called by sbc_ice_lim
+   PUBLIC   lim_sbc_flx   ! called by sbc_ice_lim
+   PUBLIC   lim_sbc_tau   ! called by sbc_ice_lim
 
    REAL(wp)  ::   epsi16 = 1.e-16  ! constant values
@@ -41,9 +42,11 @@
    REAL(wp)  ::   rone   = 1.e0
 
+   REAL(wp), DIMENSION(jpi,jpj) ::   utau_oce   !: atm.-ocean surface i-stress (ocean referential)     [N/m2]
+   REAL(wp), DIMENSION(jpi,jpj) ::   vtau_oce   !: atm.-ocean surface j-stress (ocean referential)     [N/m2]
    !! * Substitutions
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2006) 
-   !! $ Id: $
+   !! NEMO/LIM  3.0 ,  UCL-LOCEAN-IPSL (2008) 
+   !! $Id: $
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -51,14 +54,174 @@
 CONTAINS
 
-   SUBROUTINE lim_sbc( kt )
+   SUBROUTINE lim_sbc_tau( kt, kcpl )
       !!-------------------------------------------------------------------
-      !!                ***  ROUTINE lim_sbc ***
+      !!                ***  ROUTINE lim_sbc_tau ***
       !!  
-      !! ** Purpose : Update surface ocean boundary condition over areas
-      !!      that are at least partially covered by sea-ice
+      !! ** Purpose : Update the ocean surface stresses due to the ice
       !!         
-      !! ** Action  : - comput. of the momentum, heat and freshwater/salt
-      !!      fluxes at the ice-ocean interface.
-      !!              - Update 
+      !! ** 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, zfrldv   ! lead fraction at U- & V-points
+      REAL(wp) ::   ztglx , ztgly
+      REAL(wp) ::   zat_u, zu_ico, zutaui, zu_u, zv_u, zmodu, zmod
+      REAL(wp) ::   zat_v, zv_ico, zvtaui, zu_v, zv_v, zmodv, zsang
+      
+#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
+         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
+
+      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, 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_ico = u_ice(ji,jj) - u_oce(ji,jj)
+               zv_ico = v_ice(ji,jj) - v_oce(ji,jj)
+               zmod  = SQRT( zu_ico * zu_ico + zv_ico * zv_ico ) 
+               ! quadratic drag formulation
+               ztglx = rhoco * zmod * ( cangvg * zu_ico - zsang * zv_ico )
+               ztgly = rhoco * zmod * ( cangvg * zv_ico + zsang * zu_ico )
+               ! 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) = ztglx * exp( - zmod / 0.5 ) 
+               ztio_v(ji,jj) = ztgly * exp( - zmod / 0.5 ) 
+               !
+            END DO
+         END DO
+
+         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
+
+         ! boundary condition on the stress (utau,vtau)
+         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )
+         !
+         !                                        !--------------------------------!
+      CASE( 1 )                                   !  Use the now velocity          !  (at each ocean timestep)
+         !                                        !--------------------------------! 
+         ! ... save the air-ocean stresses at ice time-step
+         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
+            utau_oce(:,:) = utau(:,:)
+            vtau_oce(:,:) = vtau(:,:)
+         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 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
+
+               zu_u   = u_ice(ji,jj) - un(ji,jj,1)                  ! u ice-ocean velocity at U-point
+               zv_v   = v_ice(ji,jj) - vn(ji,jj,1)                  ! v ice-ocean velocity at V-point
+               !                                                    ! u ice-ocean velocity at V-point
+               zu_v   = 0.25 * (  u_ice(ji,jj  ) - un(ji,jj  ,1)  +  u_ice(ji+1,jj  ) - un(ji+1,jj  ,1)    &
+                  &             + u_ice(ji,jj-1) - un(ji,jj-1,1)  +  u_ice(ji+1,jj-1) - un(ji+1,jj-1,1)  ) 
+               !                                                    ! v ice-ocean velocity at U-point
+               zv_u   = 0.25 * (  v_ice(ji-1,jj+1) - vn(ji-1,jj+1,1)  +  v_ice(ji,jj+1) - vn(ji,jj+1,1)    &
+                  &             + v_ice(ji-1,jj  ) - vn(ji-1,jj  ,1)  +  v_ice(ji,jj  ) - vn(ji,jj  ,1)  )
+               zmodu    = SQRT( zu_u * zu_u + zv_u * zv_u )         ! module of the ice-ocean velocity at U-point
+               zmodv    = SQRT( zu_v * zu_v + zv_v * zv_v )         ! module of the ice-ocean velocity at V-point
+               zutaui   = rhoco * zmodu * ( cangvg * zu_u - zsang * zv_u ) 
+               zvtaui   = rhoco * zmodv * ( cangvg * zv_v + zsang * zu_v ) 
+
+               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
+         ! boundary condition on the stress (utau,vtau)
+         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )
+         !
+         !                                        !--------------------------------!
+      CASE( 2 )                                   !  mixed 0 and 2 cases           !  (at each ocean timestep)
+         !                                        !--------------------------------! 
+         ! ... save the air-ocean stresses at ice time-step
+         IF( MOD( kt-1, nn_fsbc ) == 0 ) THEN
+            utau_oce(:,:) = utau(:,:)
+            vtau_oce(:,:) = vtau(:,:)
+         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 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 * ( un(ji,jj,1) - ssu_m(ji,jj) )   ! ice-oce velocity using un and ssu_m
+               zv_ico   = v_ice(ji,jj) - 0.5 * ( vn(ji,jj,1) - ssu_m(ji,jj) )
+               !                                        ! module of the ice-ocean velocity
+               zmod     = SQRT( zu_ico * zu_ico + zv_ico * zv_ico )
+               !                                        ! quadratic drag formulation with rotation
+               zutaui   = rhoco * zmod * ( cangvg * zu_ico - zsang * zv_ico )
+               zvtaui   = rhoco * zmod * ( 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
+         ! boundary condition on the stress (utau,vtau)
+         CALL lbc_lnk( utau, 'U', -1. )   ;   CALL lbc_lnk( vtau, 'V', -1. )
+         !
+      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 )
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE lim_sbc_flx ***
+      !!  
+      !! ** Purpose :   Update the surface ocean boundary condition for heat 
+      !!              salt and mass over areas where sea-ice is non-zero
+      !!         
+      !! ** Action  : - computes the heat and freshwater/salt fluxes
+      !!              at the ice-ocean interface.
+      !!              - Update the ocean sbc
       !!     
       !! ** Outputs : - qsr    : sea heat flux:     solar 
@@ -66,6 +229,4 @@
       !!              - emp    : freshwater budget: volume flux 
       !!              - emps   : freshwater budget: concentration/dillution 
-      !!              - 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.
@@ -80,22 +241,15 @@
       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) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity
-      REAL(wp) ::   ztglx , ztgly    ! temporary scalar
-      
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfcm1 , zfcm2    ! solar/non solar heat fluxes
 #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) ::   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
-
       !!---------------------------------------------------------------------
      
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'lim_sbc : LIM 3.0 sea-ice - surface boundary condition'
-         IF(lwp) WRITE(numout,*) '~~~~~~~   '
+         IF(lwp) WRITE(numout,*) 'lim_sbc_flx : LIM 3.0 sea-ice - heat salt and mass ocean surface fluxes'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~ '
       ENDIF
 
@@ -169,5 +323,5 @@
 !                           ! fdtcn : turbulent oceanic heat flux
 
-
+!!gm   this IF prevents the vertorisation of the whole loop
             IF ( ( ji .EQ. jiindx ) .AND. ( jj .EQ. jjindx) ) THEN
                WRITE(numout,*) ' lim_sbc : heat fluxes '
@@ -198,4 +352,5 @@
                WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx)
             ENDIF
+!!gm   end
          END DO
       END DO
@@ -205,4 +360,5 @@
       !------------------------------------------!
 
+!!gm   optimisation: this loop have to be merged with the previous one
       DO jj = 1, jpj
          DO ji = 1, jpi
@@ -254,60 +410,11 @@
       END DO
 
-      emps(:,:) = fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
-      IF (num_sal.eq.2) THEN
-         !In case of variable salinity the salt flux has to be accounted for differently
-         ! Brine drainage has to be added
+      IF( num_sal == 2 ) THEN      ! variable ice salinity: brine drainage included in the salt flux
          emps(:,:) = fsbri(:,:) + fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
+      ELSE                         ! constant ice salinity:
+         emps(:,:) =              fseqv(:,:) + fsalt_res(:,:) + fsalt_rpo(:,:) + emp(:,:)
       ENDIF
       
       IF( lk_dynspg_rl )    emp (:,:) = emps(:,:)      ! rigid-lid formulation : emp = emps
-
-      !------------------------------------------!
-      !    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
-         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 = rhoco * zmod * ( cangvg * zu_io - zsang * zv_io )
-               ztgly = rhoco * zmod * ( cangvg * zv_io + zsang * zu_io )
-               ! 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) = ztglx * exp( - zmod / 0.5 ) 
-               ztio_v(ji,jj) = ztgly * exp( - zmod / 0.5 ) 
-               !
-            END DO
-         END DO
-
-         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
-
-         ! boundary condition on the stress (utau,vtau)
-         CALL lbc_lnk( utau, 'U', -1. )
-         CALL lbc_lnk( vtau, 'V', -1. )
-
-      ENDIF
 
       !-----------------------------------------------!
@@ -331,13 +438,12 @@
 
       IF(ln_ctl) THEN
-         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    : ')
-         CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
-            &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask )
-         CALL prt_ctl(tab2d_1=freeze, clinfo1=' lim_sbc: freeze : ')
-         CALL prt_ctl(tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl)
+         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    : ' )
+         CALL prt_ctl( tab2d_1=freeze, clinfo1=' lim_sbc: freeze : ' )
+         CALL prt_ctl( tab3d_1=tn_ice, clinfo1=' lim_sbc: tn_ice : ', kdim=jpl )
       ENDIF 
-   
-    END SUBROUTINE lim_sbc
+      ! 
+   END SUBROUTINE lim_sbc_flx
+
 
 #else
