Index: trunk/NEMO/C1D_SRC/step1d.F90
===================================================================
--- trunk/NEMO/C1D_SRC/step1d.F90	(revision 716)
+++ trunk/NEMO/C1D_SRC/step1d.F90	(revision 717)
@@ -154,5 +154,4 @@
          CALL prt_ctl(tab2d_1=qns    , clinfo1=' qns  -   : ', mask1=tmask, ovlap=1)
          CALL prt_ctl(tab2d_1=qsr    , clinfo1=' qsr  -   : ', mask1=tmask, ovlap=1)
-         CALL prt_ctl(tab2d_1=runoff , clinfo1=' runoff   : ', mask1=tmask, ovlap=1)
          CALL prt_ctl(tab3d_1=tmask  , clinfo1=' tmask    : ', mask1=tmask, ovlap=1, kdim=jpk)
          CALL prt_ctl(tab3d_1=tn     , clinfo1=' sst  -   : ', mask1=tmask, ovlap=1, kdim=1)
Index: trunk/NEMO/LIM_SRC/ice.F90
===================================================================
--- trunk/NEMO/LIM_SRC/ice.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/ice.F90	(revision 717)
@@ -4,16 +4,10 @@
    !! Sea Ice physics:  diagnostics variables of ice defined in memory
    !!=====================================================================
+   !! History :  2.0  !  03-08  (C. Ethe)  F90: Free form and module
+   !!----------------------------------------------------------------------
 #if defined key_ice_lim
    !!----------------------------------------------------------------------
    !!   'key_ice_lim' :                                   LIM sea-ice model
    !!----------------------------------------------------------------------
-   !! History :
-   !!   2.0  !  03-08  (C. Ethe)  F90: Free form and module
-   !!----------------------------------------------------------------------
-   !!  LIM 2.0, UCL-LOCEAN-IPSL (2005)
-   !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
-   !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_ice          ! LIM sea-ice parameters
 
@@ -21,103 +15,85 @@
    PRIVATE
 
-   !! * Share Module variables
-   INTEGER , PUBLIC ::   & !!: ** ice-dynamic namelist (namicedyn) **
-      nbiter = 1      ,  &  !: number of sub-time steps for relaxation
-      nbitdr = 250          !: maximum number of iterations for relaxation
+   !!* ice-dynamic namelist (namicedyn) *
+   INTEGER , PUBLIC ::   nbiter = 1         !: number of sub-time steps for relaxation
+   INTEGER , PUBLIC ::   nbitdr = 250       !: maximum number of iterations for relaxation
+   REAL(wp), PUBLIC ::   epsd   = 1.0e-20   !: tolerance parameter for dynamic
+   REAL(wp), PUBLIC ::   alpha  = 0.5       !: coefficient for semi-implicit coriolis
+   REAL(wp), PUBLIC ::   dm     = 0.6e+03   !: diffusion constant for dynamics
+   REAL(wp), PUBLIC ::   om     = 0.5       !: relaxation constant
+   REAL(wp), PUBLIC ::   resl   = 5.0e-05   !: maximum value for the residual of relaxation
+   REAL(wp), PUBLIC ::   cw     = 5.0e-03   !: drag coefficient for oceanic stress
+   REAL(wp), PUBLIC ::   angvg  = 0.e0      !: turning angle for oceanic stress
+   REAL(wp), PUBLIC ::   pstar  = 1.0e+04   !: first bulk-rheology parameter
+   REAL(wp), PUBLIC ::   c_rhg  = 20.e0     !: second bulk-rhelogy parameter
+   REAL(wp), PUBLIC ::   etamn  = 0.e+07    !: minimun value for viscosity
+   REAL(wp), PUBLIC ::   creepl = 2.e-08    !: creep limit
+   REAL(wp), PUBLIC ::   ecc    = 2.e0      !: eccentricity of the elliptical yield curve
+   REAL(wp), PUBLIC ::   ahi0   = 350.e0    !: sea-ice hor. eddy diffusivity coeff. (m2/s)
 
-   REAL(wp), PUBLIC ::   & !!: ** ice-dynamic namelist (namicedyn) **
-      epsd   = 1.0e-20,  &  !: tolerance parameter for dynamic
-      alpha  = 0.5    ,  &  !: coefficient for semi-implicit coriolis
-      dm     = 0.6e+03,  &  !: diffusion constant for dynamics
-      om     = 0.5    ,  &  !: relaxation constant
-      resl   = 5.0e-05,  &  !: maximum value for the residual of relaxation
-      cw     = 5.0e-03,  &  !: drag coefficient for oceanic stress
-      angvg  = 0.e0   ,  &  !: turning angle for oceanic stress
-      pstar  = 1.0e+04,  &  !: first bulk-rheology parameter
-      c_rhg  = 20.e0  ,  &  !: second bulk-rhelogy parameter
-      etamn  = 0.e+07,   &  !: minimun value for viscosity
-      creepl = 2.e-08,   &  !: creep limit
-      ecc    = 2.e0   ,  &  !: eccentricity of the elliptical yield curve
-      ahi0   = 350.e0       !: sea-ice hor. eddy diffusivity coeff. (m2/s)
+   REAL(wp), PUBLIC ::   usecc2             !:  = 1.0 / ( ecc * ecc )
+   REAL(wp), PUBLIC ::   rhoco              !: = rau0 * cw
+   REAL(wp), PUBLIC ::   sangvg, cangvg     !: sin and cos of the turning angle for ocean stress
+   REAL(wp), PUBLIC ::   pstarh             !: pstar / 2.0
 
-   REAL(wp), PUBLIC ::   &  !:
-      usecc2          ,  &  !:  = 1.0 / ( ecc * ecc )
-      rhoco           ,  &  !: = rau0 * cw
-      sangvg, cangvg  ,  &  !: sin and cos of the turning angle for ocean stress
-      pstarh                !: pstar / 2.0
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ahiu , ahiv   !: hor. diffusivity coeff. at ocean U- and V-points (m2/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pahu , pahv   !: ice hor. eddy diffusivity coef. at ocean U- and V-points
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnm , hicm   !: mean snow and ice thicknesses
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ust2s                 !: friction velocity
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  &  !:
-      u_oce, v_oce,      &  !: surface ocean velocity used in ice dynamics
-      ahiu , ahiv ,      &  !: hor. diffusivity coeff. at ocean U- and V-points (m2/s)
-      pahu , pahv ,      &  !: ice hor. eddy diffusivity coef. at ocean U- and V-points
-      hsnm , hicm ,      &  !: mean snow and ice thicknesses
-      ust2s                 !: friction velocity
+   !!* diagnostic quantities
+!! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   firic         !: IR flux over the ice (only used for outputs)
+!! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fcsic         !: Sensible heat flux over the ice (only used for outputs)
+!! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fleic         !: Latent heat flux over the ice (only used for outputs)
+!! REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qlatic        !: latent flux (only used for outputs)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvosif       !: Variation of volume at surface (only used for outputs)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvobif       !: Variation of ice volume at the bottom ice (only used for outputs)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdvolif       !: Total variation of ice volume (only used for outputs)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdvonif       !: Lateral Variation of ice volume (only used for outputs)
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::  &  !:
-        sst_ini,         &  !: sst read from a file for ice model initialization 
-        sss_ini             !: sss read from a file for ice model initialization 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sist          !: Sea-Ice Surface Temperature (Kelvin ??? degree ??? I don't know)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tfu           !: Freezing/Melting point temperature of sea water at SSS
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicif         !: Ice thickness
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hsnif         !: Snow thickness
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   hicifp        !: Ice production/melting
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   frld          !: Leads fraction = 1-a/totalarea
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   phicif        !: ice thickness  at previous time 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   pfrld         !: Leads fraction at previous time  
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qstoif        !: Energy stored in the brine pockets
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fbif          !: Heat flux at the ice base
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmsnif       !: Variation of snow mass
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   rdmicif       !: Variation of ice mass
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qldif         !: heat balance of the lead (or of the open ocean)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qcmif         !: Energy needed to bring the ocean surface layer until its freezing 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fdtcn         !: net downward heat flux from the ice to the ocean
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qdtcn         !: energy from the ice to the ocean point (at a factor 2)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   thcm          !: part of the solar energy used in the lead heat budget
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fstric        !: Solar flux transmitted trough the ice
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ffltbif       !: Array linked with the max heat contained in brine pockets (?)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fscmbq        !: Linked with the solar flux below the ice (?)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   fsbbq         !: Also linked with the solar flux below the ice (?)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qfvbq         !: Array used to store energy in case of toral lateral ablation (?)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   dmgwi         !: Variation of the mass of snow ice
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
-      firic  ,   &  !: IR flux over the ice (only used for outputs)
-      fcsic  ,   &  !: Sensible heat flux over the ice (only used for outputs)
-      fleic  ,   &  !: Latent heat flux over the ice (only used for outputs)
-      qlatic ,   &  !: latent flux
-      rdvosif,   &  !: Variation of volume at surface (only used for outputs)
-      rdvobif,   &  !: Variation of ice volume at the bottom ice (only used for outputs)
-      fdvolif,   &  !: Total variation of ice volume (only used for outputs)
-      rdvonif,   &  !: Lateral Variation of ice volume (only used for outputs)
-      sist   ,   &  !: Sea-Ice Surface Temperature (Kelvin ??? degree ??? I don't know)
-      tfu    ,   &  !: Melting point temperature of sea water
-      hsnif  ,   &  !: Snow thickness
-      hicif  ,   &  !: Ice thickness
-      hicifp ,   &  !: Ice production/melting
-      frld   ,   &  !: Leads fraction = 1-a/totalarea
-      phicif ,   &  !: ice thickness  at previous time 
-      pfrld  ,   &  !: Leads fraction at previous time  
-      qstoif ,   &  !: Energy stored in the brine pockets
-      fbif   ,   &  !: Heat flux at the ice base
-      rdmsnif,   &  !: Variation of snow mass
-      rdmicif,   &  !: Variation of ice mass
-      qldif  ,   &  !: heat balance of the lead (or of the open ocean)
-      qcmif  ,   &  !: Energy needed to bring the ocean surface layer until its freezing 
-      fdtcn  ,   &  !: net downward heat flux from the ice to the ocean
-      qdtcn  ,   &  !: energy from the ice to the ocean
-      !             !  point (at a factor 2)
-      thcm   ,   &  !: part of the solar energy used in the lead heat budget
-      fstric ,   &  !: Solar flux transmitted trough the ice
-      ffltbif,   &  !: Array linked with the max heat contained in brine pockets (?)
-      fscmbq ,   &  !: Linked with the solar flux below the ice (?)
-      fsbbq  ,   &  !: Also linked with the solar flux below the ice (?)
-      qfvbq  ,   &  !: Array used to store energy in case of toral lateral ablation (?)
-      dmgwi         !: Variation of the mass of snow ice
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   albege        !: Albedo of the snow or ice (only for outputs)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   albecn        !: Albedo of the ocean (only for outputs)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   tauc          !: Cloud optical depth
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
-      albege ,   &  !: Albedo of the snow or ice (only for outputs)
-      albecn ,   &  !: Albedo of the ocean (only for outputs)
-      tauc   ,   &  !: Cloud optical depth
-      sdvt          !: u*^2/(Stress/density)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ui_ice, vi_ice   !: two components of the ice   velocity at I-point (m/s)
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   ui_oce, vi_oce   !: two components of the ocean velocity at I-point (m/s)
 
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpsmax)     ::   scal0   !: ???
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) ::   tbif  !: Temperature inside the ice/snow layer
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
-      u_ice, v_ice,   &  !: two components of the ice velocity (m/s)
-      tio_u, tio_v       !: two components of the ice-ocean stress (N/m2)
+!! REAL(wp), DIMENSION(jpi,jpj,0:jpkmax+1) ::   reslum        !: Relative absorption of solar radiation in each ocean level
 
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jpsmax) ::   &  !:
-      scal0              !: ???
-
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj,jplayersp1) ::   &  !:
-      tbif          !: Temperature inside the ice/snow layer
-
-   REAL(wp), DIMENSION(jpi,jpj,0:jpkmax+1) ::    &  !:
-      reslum        !: Relative absorption of solar radiation in each ocean level
-
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !:
-         sxice, syice, sxxice, syyice, sxyice,      &  !: moments for advection
-         sxsn,  sysn,  sxxsn,  syysn,  sxysn,       &  !:
-         sxa,   sya,   sxxa,   syya,   sxya,        &  !:
-         sxc0,  syc0,  sxxc0,  syyc0,  sxyc0,       &  !:
-         sxc1,  syc1,  sxxc1,  syyc1,  sxyc1,       &  !:
-         sxc2,  syc2,  sxxc2,  syyc2,  sxyc2,       &  !:
-         sxst,  syst,  sxxst,  syyst,  sxyst           !:
+   !!* moment used in the advection scheme
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxice, syice, sxxice, syyice, sxyice   !: for ice  volume
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxsn,  sysn,  sxxsn,  syysn,  sxysn    !: for snow volume
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxa,   sya,   sxxa,   syya,   sxya     !: for ice cover area
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxc0,  syc0,  sxxc0,  syyc0,  sxyc0    !: for heat content of snow
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxc1,  syc1,  sxxc1,  syyc1,  sxyc1    !: for heat content of 1st ice layer
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxc2,  syc2,  sxxc2,  syyc2,  sxyc2    !: for heat content of 2nd ice layer
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   sxst,  syst,  sxxst,  syyst,  sxyst    !: for heat content of brine pockets
 
 #else
@@ -127,4 +103,8 @@
 #endif
 
+   !!----------------------------------------------------------------------
+   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006)
+   !! $Id$
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!======================================================================
 END MODULE ice
Index: trunk/NEMO/LIM_SRC/iceini.F90
===================================================================
--- trunk/NEMO/LIM_SRC/iceini.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/iceini.F90	(revision 717)
@@ -63,11 +63,6 @@
                  
       ! Louvain la Neuve Ice model
-      IF( nacc == 1 ) THEN
-          dtsd2   = nfice * rdtmin * 0.5
-          rdt_ice = nfice * rdtmin
-      ELSE
-          dtsd2   = nfice * rdt * 0.5
-          rdt_ice = nfice * rdt
-      ENDIF
+      dtsd2   = nn_fsbc * rdttra(1) * 0.5
+      rdt_ice = nn_fsbc * rdttra(1)
 
       CALL lim_msh                    ! ice mesh initialization
Index: trunk/NEMO/LIM_SRC/limdia.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limdia.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/limdia.F90	(revision 717)
@@ -19,4 +19,5 @@
    USE par_ice         ! ice parameters
    USE ice_oce         ! ice variables
+   USE sbc_oce         ! surface boundary condition variables
    USE daymod          !
    USE dom_ice         !
@@ -87,5 +88,5 @@
 
       nv = 1 
-      vinfor(nv) = REAL( kt + nfice - 1 )
+      vinfor(nv) = REAL( kt + nn_fsbc - 1 )
       nv = nv + 1
       vinfor(nv) = nyear
@@ -107,6 +108,6 @@
                zicevol = zarea   * hicif(ji,jj)
                zsnwvol = zarea   * hsnif(ji,jj)
-               zicespd = zicevol * ( u_ice(ji,jj) * u_ice(ji,jj)   &
-                  &                + v_ice(ji,jj) * v_ice(ji,jj) )
+               zicespd = zicevol * ( ui_ice(ji,jj) * ui_ice(ji,jj)   &
+                  &                + vi_ice(ji,jj) * vi_ice(ji,jj) )
                vinfor(nv+ 1) = vinfor(nv+ 1) + zarea
                vinfor(nv+ 3) = vinfor(nv+ 3) + zextent15
@@ -133,6 +134,6 @@
                 zicevol = zarea   * hicif(ji,jj)
                 zsnwvol = zarea   * hsnif(ji,jj)
-                zicespd = zicevol * ( u_ice(ji,jj) * u_ice(ji,jj)   &
-                   &                + v_ice(ji,jj) * v_ice(ji,jj) )
+                zicespd = zicevol * ( ui_ice(ji,jj) * ui_ice(ji,jj)   &
+                   &                + vi_ice(ji,jj) * vi_ice(ji,jj) )
                 vinfor(nv+ 1) = vinfor(nv+ 1) + zarea
                 vinfor(nv+ 3) = vinfor(nv+ 3) + zextent15
@@ -154,5 +155,5 @@
     
        ! oututs on file ice_evolu    
-       IF( MOD( kt + nfice - 1, ninfo ) == 0 ) THEN
+       IF( MOD( kt + nn_fsbc - 1, ninfo ) == 0 ) THEN
           WRITE(numevo_ice,fmtw) ( titvar(jv), vinfom(jv)/naveg, jv = 1, nvinfo )
           naveg = 0
@@ -227,8 +228,8 @@
 
        ! Definition et Ecriture de l'entete : nombre d'enregistrements 
-       ndeb   = ( nit000 - 1 + nfice - 1 ) / ninfo
-       IF( nit000 - 1 + nfice == 1 ) ndeb = -1
-
-       nferme = ( nitend + nfice - 1 ) / ninfo ! nit000 - 1 + nfice - 1 + nitend - nit000 + 1
+       ndeb   = ( nit000 - 1 + nn_fsbc - 1 ) / ninfo
+       IF( nit000 - 1 + nn_fsbc == 1 ) ndeb = -1
+
+       nferme = ( nitend + nn_fsbc - 1 ) / ninfo ! nit000 - 1 + nn_fsbc - 1 + nitend - nit000 + 1
        ntot   = nferme - ndeb
        ndeb   = ninfo * ( 1 + ndeb )
Index: trunk/NEMO/LIM_SRC/limdyn.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limdyn.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/limdyn.F90	(revision 717)
@@ -4,4 +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
+   !!---------------------------------------------------------------------
 #if defined key_ice_lim
    !!----------------------------------------------------------------------
@@ -11,34 +16,31 @@
    !!    lim_dyn_init : initialization and namelist read
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE phycst
-   USE in_out_manager  ! I/O manager
-   USE dom_ice
-   USE dom_oce         ! ocean space and time domain
-   USE ice
-   USE ice_oce
-   USE iceini
-   USE limistate
-   USE limrhg          ! ice rheology
-   USE lbclnk
-   USE lib_mpp
-   USE prtctl          ! Print control
+   USE dom_oce        ! ocean space and time domain
+   USE sbc_oce        !
+   USE phycst         !
+   USE ice            !
+   USE ice_oce        !
+   USE dom_ice        !
+   USE iceini         !
+   USE limistate      !
+   USE limrhg         ! ice rheology
+
+   USE lbclnk         !
+   USE lib_mpp        !
+   USE in_out_manager ! I/O manager
+   USE prtctl         ! Print control
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Accessibility
-   PUBLIC lim_dyn  ! routine called by ice_step
-
-   !! * Module variables
+   PUBLIC   lim_dyn   ! routine called by ice_step
+
    REAL(wp)  ::  rone    = 1.e0   ! constant value
 
-
 #  include "vectopt_loop_substitute.h90"
-
-   !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
+   !!----------------------------------------------------------------------
+   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -57,28 +59,31 @@
       !!              - 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
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! number of iteration
-
-      INTEGER ::   ji, jj             ! dummy loop indices
-      INTEGER ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology
-      REAL(wp) ::   &
-         ztairx, ztairy,           &  ! tempory scalars
-         zsang , zrhomod,             &
-         ztglx , ztgly ,           &
-         zt11, zt12, zt21, zt22 ,  &
-         zustm, zsfrld, zsfrldm4,  &
-         zu_ice, zv_ice, ztair2
-      REAL(wp),DIMENSION(jpi,jpj) ::   zmod
-      REAL(wp),DIMENSION(jpj) ::   &
-         zind,                     &  ! i-averaged indicator of sea-ice
-         zmsk                         ! i-averaged of tmask
+      !!
+      INTEGER  ::   ji, jj             ! 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
+!!$      INTEGER ::   ji, jj             ! dummy loop indices
+!!$      INTEGER ::   i_j1, i_jpj        ! Starting/ending j-indices for rheology
+!!$      REAL(wp) ::   &
+!!$         ztairx, ztairy,           &  ! tempory scalars
+!!$         zsang , zrhomod,             &
+!!$         ztglx , ztgly ,           &
+!!$         zt11, zt12, zt21, zt22 ,  &
+!!$         zustm, zsfrld, zsfrldm4,  &
+!!$         zu_ice, zv_ice, ztair2
+!!$      REAL(wp),DIMENSION(jpi,jpj) ::   zmod
+!!$      REAL(wp),DIMENSION(jpj) ::   &
+!!$         zind,                     &  ! i-averaged indicator of sea-ice
+!!$         zmsk                         ! i-averaged of tmask
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000  )   CALL lim_dyn_init   ! Initialization (first time-step only)
+      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only)
       
-      IF ( ln_limdyn ) THEN
+      IF( ln_limdyn ) THEN
 
          ! Mean ice and snow thicknesses.          
@@ -86,7 +91,4 @@
          hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:)
 
-         u_io(:,:)  = u_io(:,:) * tmu(:,:)
-         v_io(:,:)  = v_io(:,:) * tmu(:,:)
-       
          !                                         ! Rheology (ice dynamics)
          !                                         ! ========
@@ -98,7 +100,5 @@
             i_j1 = 1   
             i_jpj = jpj
-            IF(ln_ctl)    THEN
-               CALL prt_ctl_info('lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj)
-            ENDIF
+            IF(ln_ctl)   CALL prt_ctl_info( 'lim_dyn  :    i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj )
             CALL lim_rhg( i_j1, i_jpj )
 
@@ -108,5 +108,4 @@
                zind(jj) = SUM( frld (:,jj  ) )   ! = FLOAT(jpj) if ocean everywhere on a j-line
                zmsk(jj) = SUM( tmask(:,jj,1) )   ! = 0          if land  everywhere on a j-line
-   !!i         write(numout,*) narea, 'limdyn' , jj, zind(jj), zmsk(jj)
             END DO
 
@@ -159,57 +158,62 @@
          ENDIF
 
-         IF(ln_ctl)   THEN 
-            CALL prt_ctl(tab2d_1=u_io , clinfo1=' lim_dyn  : u_io :', tab2d_2=v_io , clinfo2=' v_io :')
-            CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_dyn  : u_ice:', tab2d_2=v_ice, clinfo2=' v_ice:')
-         ENDIF
-         
-         !                                         ! Ice-Ocean stress
-         !                                         ! ================        
-         DO jj = 1, jpj
-            DO ji = 1, jpi
-!!              zsang  = SIGN(1.e0, gphif(ji-1,jj-1) ) * sangvg    ! do the full loop and avoid lbc_lnk
-               zsang  = SIGN(1.e0, gphif(ji,jj) ) * sangvg
-               zu_ice   = u_ice(ji,jj) - u_io(ji,jj)
-               zv_ice   = v_ice(ji,jj) - v_io(ji,jj)
-               zrhomod  = zu_ice * zu_ice + zv_ice * zv_ice 
-               zmod (ji,jj) = zrhomod 
-               zrhomod  = rhoco * SQRT( zrhomod ) 
-               ftaux(ji,jj) = zrhomod * ( cangvg * zu_ice - zsang * zv_ice ) 
-               ftauy(ji,jj) = zrhomod * ( cangvg * zv_ice + zsang * zu_ice ) 
+         IF(ln_ctl)   CALL prt_ctl(tab2d_1=ui_ice , clinfo1=' lim_dyn  : ui_ice :', tab2d_2=vi_ice , clinfo2=' vi_ice :')
+         
+!!$         !                                         ! Ice-Ocean stress
+!!$         !                                         ! ================        
+!!$         DO jj = 1, jpj
+!!$            DO ji = 1, jpi
+!!$!!              zsang  = SIGN(1.e0, gphif(ji-1,jj-1) ) * sangvg    ! do the full loop and avoid lbc_lnk
+!!$               zsang  = SIGN(1.e0, gphif(ji,jj) ) * sangvg
+!!$               zu_ice   = u_ice(ji,jj) - u_io(ji,jj)
+!!$               zv_ice   = v_ice(ji,jj) - v_io(ji,jj)
+!!$               zrhomod  = zu_ice * zu_ice + zv_ice * zv_ice 
+!!$               zmod (ji,jj) = zrhomod 
+!!$               zrhomod  = rhoco * SQRT( zrhomod ) 
+!!$               ftaux(ji,jj) = zrhomod * ( cangvg * zu_ice - zsang * zv_ice ) 
+!!$               ftauy(ji,jj) = zrhomod * ( cangvg * zv_ice + zsang * zu_ice ) 
+!!$            END DO
+!!$         END DO
+
+         ! computation of friction velocity
+         ! --------------------------------
+         ! ice-ocean velocity at U & V-points (ui_ice vi_ice at I-point ; ssu_m, ssv_m at U- & V-points)
+         
+         DO jj = 1, jpjm1
+            DO ji = 1, fs_jpim1   ! vector opt.
+               zu_io(ji,jj) = 0.5 * ( ui_ice(ji+1,jj+1) + ui_ice(ji+1,jj  ) ) - ssu_m(ji,jj)
+               zv_io(ji,jj) = 0.5 * ( vi_ice(ji+1,jj+1) + vi_ice(ji  ,jj+1) ) - ssv_m(ji,jj)
             END DO
          END DO
-
-         ! computation of friction velocity
+         ! frictional velocity at T-point
          DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1     ! vector opt.
-               ust2s(ji,jj) = 0.25 * cw * ( zmod(ji,jj+1) + zmod(ji+1,jj+1) +    &
-                  &                         zmod(ji,jj  ) + zmod(ji+1,jj  ) ) * tms(ji,jj)
+            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
-                    
-          ftaux(:,:) = gtaux(:,:) 
-          ftauy(:,:) = gtauy(:,:) 
-
-          DO jj = 2, jpjm1
-             DO ji = 2, jpim1
-                ztair2 = gtaux(ji  ,jj+1) * gtaux(ji  ,jj+1) + gtauy(ji  ,jj+1) * gtauy(ji  ,jj+1)   &
-                &      + gtaux(ji+1,jj+1) * gtaux(ji+1,jj+1) + gtauy(ji+1,jj+1) * gtauy(ji+1,jj+1)   &
-                &      + gtaux(ji  ,jj  ) * gtaux(ji  ,jj  ) + gtauy(ji  ,jj  ) * gtauy(ji  ,jj  )   &
-                &      + gtaux(ji+1,jj  ) * gtaux(ji+1,jj  ) + gtauy(ji+1,jj  ) * gtauy(ji+1,jj  ) 
-
-                ust2s(ji,jj) = 0.25 / rau0 * SQRT( ztair2 ) * tms(ji,jj)
+!!$         DO jj = 2, jpjm1
+!!$            DO ji = fs_2, fs_jpim1     ! vector opt.
+!!$               ust2s(ji,jj) = 0.25 * cw * ( zmod(ji,jj+1) + zmod(ji+1,jj+1) +    &
+!!$                  &                         zmod(ji,jj  ) + zmod(ji+1,jj  ) ) * 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   ! 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
-
+         !
       ENDIF
 
       CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
 
-      IF(ln_ctl) THEN 
-            CALL prt_ctl(tab2d_1=ftaux , clinfo1=' lim_dyn  : ftaux :', tab2d_2=ftauy , clinfo2=' ftauy :')
-            CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
-      ENDIF
+      IF(ln_ctl)   CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn  : ust2s :')
 
    END SUBROUTINE lim_dyn
@@ -220,14 +224,11 @@
       !!                  ***  ROUTINE lim_dyn_init  ***
       !!
-      !! ** Purpose : Physical constants and parameters linked to the ice
-      !!      dynamics
-      !!
-      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
-      !!       parameter values called at the first timestep (nit000)
+      !! ** Purpose :   Physical constants and parameters linked to the ice
+      !!              dynamics
+      !!
+      !! ** Method  :   Read the namicedyn namelist and check the ice-dynamic
+      !!              parameter values
       !!
       !! ** input   :   Namelist namicedyn
-      !!
-      !! history :
-      !!  8.5  ! 03-08 (C. Ethe) original code
       !!-------------------------------------------------------------------
       NAMELIST/namicedyn/ epsd, alpha,     &
@@ -236,11 +237,8 @@
       !!-------------------------------------------------------------------
 
-      ! Define the initial parameters
-      ! -------------------------
-
-      ! Read Namelist namicedyn
-      REWIND ( numnam_ice )
+      REWIND ( numnam_ice )                       ! Read Namelist namicedyn
       READ   ( numnam_ice  , namicedyn )
-      IF(lwp) THEN
+
+      IF(lwp) THEN                                ! Control print
          WRITE(numout,*)
          WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
@@ -254,5 +252,5 @@
          WRITE(numout,*) '       maximum value for the residual of relaxation     resl   = ', resl
          WRITE(numout,*) '       drag coefficient for oceanic stress              cw     = ', cw
-         WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg
+         WRITE(numout,*) '       turning angle for oceanic stress                 angvg  = ', angvg, ' degrees'
          WRITE(numout,*) '       first bulk-rheology parameter                    pstar  = ', pstar
          WRITE(numout,*) '       second bulk-rhelogy parameter                    c_rhg  = ', c_rhg
@@ -263,17 +261,14 @@
       ENDIF
 
-      !  Initialization
       usecc2 = 1.0 / ( ecc * ecc )
       rhoco  = rau0 * cw
-      angvg  = angvg * rad
+      angvg  = angvg * rad      ! convert angvg from degree to radian
       sangvg = SIN( angvg )
       cangvg = COS( angvg )
       pstarh = pstar / 2.0
-      sdvt(:,:) = 0.e0
-
-      !  Diffusion coefficients.
-      ahiu(:,:) = ahi0 * umask(:,:,1)
+      !
+      ahiu(:,:) = ahi0 * umask(:,:,1)            ! Ice eddy Diffusivity coefficients.
       ahiv(:,:) = ahi0 * vmask(:,:,1)
-
+      !
    END SUBROUTINE lim_dyn_init
 
Index: trunk/NEMO/LIM_SRC/limistate.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limistate.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/limistate.F90	(revision 717)
@@ -4,7 +4,9 @@
    !!              Initialisation of diagnostics ice variables
    !!======================================================================
-   !! History :   2.0  !  01-04  (C. Ethe, G. Madec)  Original code
+   !! History :   1.0  !  01-04  (C. Ethe, G. Madec)  Original code
+   !!             2.0  !  03-08  (G. Madec)  add lim_istate_init
    !!                  !  04-04  (S. Theetten) initialization from a file
    !!                  !  06-07  (S. Masson)  IOM to read the restart
+   !!                  !  07-10  (G. Madec)  surface module
    !!--------------------------------------------------------------------
 #if defined key_ice_lim
@@ -18,11 +20,9 @@
    USE phycst
    USE ocfzpt
-   USE oce             ! dynamics and tracers variables      !!gm used???
-   USE dom_oce                                                     !!gm used???
    USE par_ice         ! ice parameters
    USE ice_oce         ! ice variables
    USE dom_ice
-   USE ice             ! ???
    USE lbclnk
+   USE oce
    USE ice
    USE iom
@@ -32,7 +32,7 @@
    PRIVATE
 
-   PUBLIC lim_istate      ! routine called by lim_init.F90
-
-   REAL(wp) ::           &  !!! ** init namelist (namiceini) **
+   PUBLIC   lim_istate   ! routine called by lim_init.F90
+
+   REAL(wp) ::           & !!! ** init namelist (namiceini) **
       ttest  = 2.0  ,    &  ! threshold water temperature for initial sea ice
       hninn  = 0.5  ,    &  ! initial snow thickness in the north
@@ -67,46 +67,8 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   ztn   ! workspace
       !--------------------------------------------------------------------
-
-       CALL lim_istate_init     !  reading the initials parameters of the ice
-
-      !-- Initialisation of sst,sss,u,v do i=1,jpi
-      u_io(:,:)  = 0.e0       ! ice velocity in x direction
-      v_io(:,:)  = 0.e0       ! ice velocity in y direction
-
-      IF( ln_limini ) THEN    ! 
-        
-         ! Initialisation at tn if no ice or sst_ini if ice
-         ! Idem for salinity
-
-      !--- Criterion for presence (zidto=1.) or absence (zidto=0.) of ice
-         DO jj = 1 , jpj
-            DO ji = 1 , jpi
-               
-               zidto = MAX(zzero, - SIGN(1.,frld(ji,jj) - 1.))
-               
-               sst_io(ji,jj) = ( nfice - 1 ) * (zidto * sst_ini(ji,jj)  + &   ! use the ocean initial values
-                    &          (1.0 - zidto ) * ( tn(ji,jj,1) + rt0 ))        ! tricky trick *(nfice-1) !
-               sss_io(ji,jj) = ( nfice - 1 ) * (zidto * sss_ini(ji,jj) + &
-                    &          (1.0 - zidto ) *  sn(ji,jj,1) )
-
-               ! to avoid the the melting of ice, several layers (mixed layer) should be
-               ! set to sst_ini (sss_ini) if there is ice
-               ! example for one layer 
-               ! tn(ji,jj,1) = zidto * ( sst_ini(ji,jj) - rt0 )  + (1.0 - zidto ) *  tn(ji,jj,1)
-               ! sn(ji,jj,1) = zidto * sss_ini(ji,jj)  + (1.0 - zidto ) *  sn(ji,jj,1)
-               ! tb(ji,jj,1) = tn(ji,jj,1)
-               ! sb(ji,jj,1) = sn(ji,jj,1)
-            END DO
-         END DO
-         
-         
-         !  tfu: Melting point of sea water
-         tfu(:,:)  = ztf   
-         
-         tfu(:,:)  = ABS ( rt0 - 0.0575       * sss_ini(:,:)                               &
-              &                    + 1.710523e-03 * sss_ini(:,:) * SQRT( sss_ini(:,:) )    &
-              &                    - 2.154996e-04 * sss_ini(:,:) * sss_ini(:,:) )
-      ELSE                     !
-
+ 
+      CALL lim_istate_init     !  reading the initials parameters of the ice
+
+      IF( .NOT. ln_limini ) THEN  
          
          ! Initialisation at tn or -2 if ice
@@ -117,18 +79,9 @@
             END DO
          END DO
-         
-         u_io  (:,:) = 0.e0
-         v_io  (:,:) = 0.e0
-         sst_io(:,:) = ( nfice - 1 ) * ( tn(:,:,1) + rt0 )   ! use the ocean initial values
-         sss_io(:,:) = ( nfice - 1 ) *   sn(:,:,1)           ! tricky trick *(nfice-1) !
-         
-         ! reference salinity 34psu
+                  
+         !  tfu: Melting point of sea water [Kelvin]
          zs0 = 34.e0
-         ztf = ABS ( rt0 - 0.0575       * zs0                           &
-              &                    + 1.710523e-03 * zs0 * SQRT( zs0 )   &
-              &                    - 2.154996e-04 * zs0 *zs0          )
-         
-         !  tfu: Melting point of sea water
-         tfu(:,:)  = ztf   
+         ztf = rt0 + ( - 0.0575 + 1.710523e-3 * SQRT( zs0 ) - 2.154996e-4 * zs0 ) * zs0
+         tfu(:,:) = ztf
          
          DO jj = 1, jpj
@@ -153,10 +106,11 @@
          tbif  (:,:,2) = tfu(:,:)
          tbif  (:,:,3) = tfu(:,:)
-      
+
       ENDIF
+     
       fsbbq (:,:)   = 0.e0
       qstoif(:,:)   = 0.e0
-      u_ice (:,:)   = 0.e0
-      v_ice (:,:)   = 0.e0
+      ui_ice(:,:)   = 0.e0
+      vi_ice(:,:)   = 0.e0
 # if defined key_coupled
       albege(:,:)   = 0.8 * tms(:,:)
@@ -192,5 +146,5 @@
 
       CALL lbc_lnk( hsnif, 'T', 1. )
-      CALL lbc_lnk( sist , 'T', 1. )
+      CALL lbc_lnk( sist , 'T', 1. , pval = rt0 )      ! set rt0 on closed boundary (required by bulk formulation)
       DO jk = 1, jplayersp1
          CALL lbc_lnk(tbif(:,:,jk), 'T', 1. )
@@ -198,6 +152,5 @@
       CALL lbc_lnk( fsbbq  , 'T', 1. )
       CALL lbc_lnk( qstoif , 'T', 1. )
-      CALL lbc_lnk( sss_io , 'T', 1. )
-      !
+
    END SUBROUTINE lim_istate
 
@@ -210,8 +163,5 @@
       !!
       !! ** Method  :   Read the namiceini namelist and check the parameter 
-      !!                values called at the first timestep (nit000)
-      !!                or
-      !!                Read 7 variables from a previous restart file
-      !!                sst, sst, hicif, hsnif, frld, ts & tbif
+      !!       values called at the first timestep (nit000)
       !!
       !! ** input   :   Namelist namiceini
@@ -223,43 +173,34 @@
          &                hnins, hgins, alins
       !!-------------------------------------------------------------------
-      
-      ! Read Namelist namiceini 
-      REWIND ( numnam_ice )
+      !
+      REWIND ( numnam_ice )               ! Read Namelist namiceini 
       READ   ( numnam_ice , namiceini )
-      
-      IF(.NOT. ln_limini) THEN 
-         IF(lwp) THEN
-            WRITE(numout,*)
-            WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
-            WRITE(numout,*) '~~~~~~~~~~~~~~~'
-            WRITE(numout,*) '         threshold water temp. for initial sea-ice    ttest      = ', ttest
-            WRITE(numout,*) '         initial snow thickness in the north          hninn      = ', hninn
-            WRITE(numout,*) '         initial ice thickness in the north           hginn      = ', hginn 
-            WRITE(numout,*) '         initial leads area in the north              alinn      = ', alinn            
-            WRITE(numout,*) '         initial snow thickness in the south          hnins      = ', hnins 
-            WRITE(numout,*) '         initial ice thickness in the south           hgins      = ', hgins
-            WRITE(numout,*) '         initial leads area in the south              alins      = ', alins
-         ENDIF
+
+      IF(lwp) THEN
+         WRITE(numout,*)
+         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation '
+         WRITE(numout,*) '~~~~~~~~~~~~~~~'
+         WRITE(numout,*) '         threshold water temp. for initial sea-ice    ttest      = ', ttest
+         WRITE(numout,*) '         initial snow thickness in the north          hninn      = ', hninn
+         WRITE(numout,*) '         initial ice thickness in the north           hginn      = ', hginn 
+         WRITE(numout,*) '         initial leads area in the north              alinn      = ', alinn            
+         WRITE(numout,*) '         initial snow thickness in the south          hnins      = ', hnins 
+         WRITE(numout,*) '         initial ice thickness in the south           hgins      = ', hgins
+         WRITE(numout,*) '         initial leads area in the south              alins      = ', alins
+         WRITE(numout,*) '         Ice state initialization using input file    ln_limini  = ', ln_limini
       ENDIF
 
       IF( ln_limini ) THEN                      ! Ice initialization using input file
-
+         !
          CALL iom_open( 'Ice_initialization.nc', inum_ice )
-
+         !
          IF( inum_ice > 0 ) THEN
-            IF(lwp) THEN
-               WRITE(numout,*) ' '
-               WRITE(numout,*) 'lim_istate_init : ice state initialization with : Ice_initialization.nc'
-               WRITE(numout,*) '~~~~~~~~~~~~~~~'
-               WRITE(numout,*) '         Ice state initialization using input file    ln_limini  = ', ln_limini
-               WRITE(numout,*) ' '
-            ENDIF
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '                  ice state initialization with : Ice_initialization.nc'
             
-            CALL iom_get( inum_ice, jpdom_data, 'sst'  , sst_ini(:,:) )        
-            CALL iom_get( inum_ice, jpdom_data, 'sss'  , sss_ini(:,:) )       
-            CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif  (:,:) )      
-            CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif  (:,:) )      
-            CALL iom_get( inum_ice, jpdom_data, 'frld' , frld   (:,:) )     
-            CALL iom_get( inum_ice, jpdom_data, 'ts'   , sist   (:,:) )
+            CALL iom_get( inum_ice, jpdom_data, 'hicif', hicif )      
+            CALL iom_get( inum_ice, jpdom_data, 'hsnif', hsnif )      
+            CALL iom_get( inum_ice, jpdom_data, 'frld' , frld  )     
+            CALL iom_get( inum_ice, jpdom_data, 'ts'   , sist  )
             CALL iom_get( inum_ice, jpdom_unknown, 'tbif', tbif(1:nlci,1:nlcj,:),   &
                  &        kstart = (/ mig(1),mjg(1),1 /), kcount = (/ nlci,nlcj,jplayersp1 /) )
@@ -269,8 +210,8 @@
 
             CALL iom_close( inum_ice)
-            
+            !
          ENDIF
       ENDIF
-      !
+      !     
    END SUBROUTINE lim_istate_init
 
Index: trunk/NEMO/LIM_SRC/limrhg.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limrhg.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/limrhg.F90	(revision 717)
@@ -4,9 +4,9 @@
    !!   Ice rheology :  performs sea ice rheology
    !!======================================================================
-   !! History :   0.0  !  93-12  (M.A. Morales Maqueda.)  Original code
-   !!             1.0  !  94-12  (H. Goosse) 
-   !!             2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp
-   !!             3.0  !  07-06  (S. Masson, G. Madec)  ice/atmosphere stress
-   !!                             provided at I-point in forced and coupled mode
+   !! History :  0.0  !  93-12  (M.A. Morales Maqueda.)  Original code
+   !!            1.0  !  94-12  (H. Goosse) 
+   !!            2.0  !  03-12  (C. Ethe, G. Madec)  F90, mpp
+   !!            " "  !  06-08  (G. Madec)  surface module, ice-stress at I-point
+   !!            " "  !  09-09  (G. Madec)  Huge verctor optimisation
    !!----------------------------------------------------------------------
 #if defined key_ice_lim
@@ -14,31 +14,32 @@
    !!   'key_ice_lim'                                     LIM sea-ice model
    !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
    !!   lim_rhg   : computes ice velocities
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE phycst
-   USE par_oce
-   USE ice_oce         ! ice variables
-   USE dom_ice
-   USE ice
-   USE lbclnk
-   USE lib_mpp
-   USE in_out_manager  ! I/O manager
-   USE prtctl          ! Print control
+   USE par_oce        ! ocean parameter
+   USE ice_oce        ! ice variables
+   USE sbc_ice        ! surface boundary condition: ice variables
+   USE dom_ice        ! domaine: ice variables
+   USE phycst         ! physical constant
+   USE ice            ! ice variables
+   USE lbclnk         ! lateral boundary condition
+   USE lib_mpp        ! MPP library
+   USE in_out_manager ! I/O manager
+   USE prtctl         ! Print control
 
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_rhg  ! routine called by lim_dyn
-
-   !! * Module variables
-   REAL(wp)  ::           &  ! constant values
-      rzero   = 0.e0   ,  &
-      rone    = 1.e0
-   !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Id$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   PUBLIC   lim_rhg   ! routine called by lim_dyn
+
+   REAL(wp) ::   rzero   = 0.e0   ! constant value: zero
+   REAL(wp) ::   rone    = 1.e0   !            and  one
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! $Header: $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -54,80 +55,150 @@
       !!  viscous-plastic law including shear strength and a bulk rheology.
       !!
-      !! ** Action  : - compute u_ice, v_ice the sea-ice velocity
+      !! ** Action  : - compute ui_ice, vi_ice the sea-ice velocity defined
+      !!              at I-point
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
+      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation
       !!
+      INTEGER ::   ji, jj              ! dummy loop indices
+      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, zunw, zur, zusw    !    "         "
+      REAL(wp) ::   zgphsy, ztagny, zvnw, zvr          !    "         "
+      REAL(wp) ::   zresm,  za, zac, zmod
+      REAL(wp) ::   zmpzas, zstms, zindu, zusdtp, zmassdt, zcorlal
+      REAL(wp) ::   ztrace2, zdeter, zdelta, zmask, zdgp, zdgi, zdiag
+      REAL(wp) ::   za1, zb1, zc1, zd1
+      REAL(wp) ::   za2, zb2, zc2, zd2, zden
+      REAL(wp) ::   zs11_11, zs11_12, zs11_21, zs11_22
+      REAL(wp) ::   zs12_11, zs12_12, zs12_21, zs12_22
+      REAL(wp) ::   zs21_11, zs21_12, zs21_21, zs21_22
+      REAL(wp) ::   zs22_11, zs22_12, zs22_21, zs22_22
+      REAL(wp), DIMENSION(jpi,  jpj  ) ::   zfrld, zmass, zcorl
+      REAL(wp), DIMENSION(jpi,  jpj  ) ::   za1ct, za2ct, zresr
+      REAL(wp), DIMENSION(jpi,  jpj  ) ::   zc1u, zc1v, zc2u, zc2v
+      REAL(wp), DIMENSION(jpi,  jpj  ) ::   zsang
+      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zu0, zv0
+      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zu_n, zv_n
+      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zu_a, zv_a
+      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zviszeta, zviseta
+      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zzfrld, zztms
+      REAL(wp), DIMENSION(jpi,0:jpj+1) ::   zi1, zi2, zmasst, zpresh
+
       !!-------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   k_j1 ,  &  ! southern j-index for ice computation
-         &                     k_jpj      ! northern j-index for ice computation
-
-      INTEGER ::   ji, jj              ! dummy loop indices
-
-      INTEGER  :: &
-         iim1, ijm1, iip1 , ijp1   , & ! temporary integers
-         iter, jter                    !    "          "
-
-      CHARACTER (len=50) :: charout
-
-      REAL(wp) :: &
-         ze11  , ze12  , ze22  , ze21  ,   &  ! temporary scalars
-         zt11  , zt12  , zt21  , zt22  ,   &  !    "         "
-         zvis11, zvis21, zvis12, zvis22,   &  !    "         "
-         zgphsx, ztagnx, zusw  ,           &  !    "         "
-         zgphsy, ztagny                       !    "         "
-      REAL(wp) :: &
-         zresm, zunw, zvnw, zur, zvr, zmod, za, zac, &
-         zmpzas, zstms, zindu, zindu1, zusdtp, zmassdt, zcorlal,  &
-         ztrace2, zdeter, zdelta, zsang, zmask, zdgp, zdgi, zdiag, zic
-      REAL(wp),DIMENSION(jpi,jpj) ::   &
-         zpresh, zfrld, zmass, zcorl,     &
-         zu0, zv0, zviszeta, zviseta,     &
-         zc1u, zc1v, zc2u, zc2v, za1ct, za2ct, za1, za2, zb1, zb2,  &
-         zc1, zc2, zd1, zd2, zden, zu_ice, zv_ice, zresr
-      REAL(wp),DIMENSION(jpi,jpj,2,2) :: &
-         zs11, zs12, zs22, zs21
-      !!-------------------------------------------------------------------
+
+!!bug
+!!    ui_oce(:,:) = 0.e0
+!!    vi_oce(:,:) = 0.e0
+!!    write(*,*) 'rhg min, max u & v', maxval(ui_oce), minval(ui_oce), maxval(vi_oce), minval(vi_oce)
+!!bug
       
       !  Store initial velocities
-      !  ------------------------
-      zu0(:,:) = u_ice(:,:)
-      zv0(:,:) = v_ice(:,:)
-
-      ! masked Ice mass, ice strength, Ice-cover fraction and Lead faction  at T-point 
+      !  ----------------
+      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) = ui_ice(:,:)   ;    zv0(:,1:jpj) = vi_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
+!i
+!!gm violant
+      zfrld(:,:) =0.e0
+      zcorl(:,:) =0.e0
+      zmass(:,:) =0.e0
+      za1ct(:,:) =0.e0
+      za2ct(:,:) =0.e0
+!!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
+
+
+      ! Ice mass, ice strength, and wind stress at the center            |
+      ! of the grid squares.                                             |
       !-------------------------------------------------------------------
 
+!CDIR NOVERRCHK
       DO jj = k_j1 , k_jpj-1
+!CDIR NOVERRCHK
          DO ji = 1 , jpi
-            za1(ji,jj)    = tms(ji,jj) * ( rhosn * hsnm(ji,jj) + rhoic * hicm(ji,jj) )
+            ! 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
+            !
+            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) )
-
-            zb1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )    
-            zb2(ji,jj)    = tms(ji,jj) *         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) )
+
+!!gm   #if defined key_lim_cp1 && defined key_coupled
+!!gm        zi1(ji,jj)    = tms(ji,jj) * gtaux(ji,jj) * ( 1.0 - frld(ji,jj) )
+!!gm        zi2(ji,jj)    = tms(ji,jj) * gtauy(ji,jj) * ( 1.0 - frld(ji,jj) )
+!!gm   #else
+!!gm        zi1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )
+!!gm        zi2(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )
+!!gm   #endif
          END DO
       END DO
 
-      !  Wind stress, coriolis, mass and Gradient of ice strenght at I-point
+
+      !---------------------------------------------------------------------------
+      !  Wind stress, coriolis and mass terms at the corners of the grid squares |
+      !  Gradient of ice strenght.                                               |
       !---------------------------------------------------------------------------
          
       DO jj = k_j1+1, k_jpj-1
-         DO ji = 2, jpi
-            zstms = tms(ji,jj  ) * wght(ji,jj,2,2) + tms(ji-1,jj  ) * wght(ji,jj,1,2)   &
-               &  + tms(ji,jj-1) * wght(ji,jj,2,1) + tms(ji-1,jj-1) * wght(ji,jj,1,1)
+         DO ji = fs_2, jpi
+            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 )
 
-            ! Leads area at I-point
-            zfrld(ji,jj) = ( zb2(ji,jj  ) * wght(ji,jj,2,2) + zb2(ji-1,jj  ) * wght(ji,jj,1,2)   &
-               &           + zb2(ji,jj-1) * wght(ji,jj,2,1) + zb2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
-
-            ! Ice cover area at I-point
-            zic          = ( zb1(ji,jj  ) * wght(ji,jj,2,2) + zb1(ji-1,jj  ) * wght(ji,jj,1,2)   &
-               &           + zb1(ji,jj-1) * wght(ji,jj,2,1) + zb1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
+            zt11 = zztms(ji  ,jj  ) * zzfrld(ji  ,jj  ) 
+            zt12 = zztms(ji-1,jj  ) * zzfrld(ji-1,jj  ) 
+            zt21 = zztms(ji  ,jj-1) * zzfrld(ji  ,jj-1) 
+            zt22 = zztms(ji-1,jj-1) * zzfrld(ji-1,jj-1)
+
+            ! Leads area.
+            zfrld(ji,jj) =  (  zt11 * wght(ji,jj,2,2) + zt12 * wght(ji,jj,1,2)   &
+               &             + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw
+
+            ! Mass and coriolis coeff. at I-point
+            zmass(ji,jj) = ( zmasst(ji,jj  ) * wght(ji,jj,2,2) + zmasst(ji-1,jj  ) * wght(ji,jj,1,2)   &
+               &           + zmasst(ji,jj-1) * wght(ji,jj,2,1) + zmasst(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
+            zcorl(ji,jj) = zmass(ji,jj) * fcor(ji,jj)
 
             ! Wind stress.
-            ztagnx = zic * gtaux(ji,jj)
-            ztagny = zic * gtauy(ji,jj)
-
-            ! Mass and coriolis coeff.
-            zmass(ji,jj) = ( za1(ji,jj  ) * wght(ji,jj,2,2) + za1(ji-1,jj  ) * wght(ji,jj,1,2)   &
-               &           + za1(ji,jj-1) * wght(ji,jj,2,1) + za1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
-            zcorl(ji,jj) = zmass(ji,jj) * fcor(ji,jj)
+!!gm   #if defined key_lim_cp1 && defined key_coupled
+!!gm   ztagnx = ( zi1(ji,jj  ) * wght(ji,jj,2,2) + zi1(ji-1,jj  ) * wght(ji,jj,1,2)   &
+!!gm      &     + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
+!!gm   ztagny = ( zi2(ji,jj  ) * wght(ji,jj,2,2) + zi2(ji-1,jj  ) * wght(ji,jj,1,2)   &
+!!gm      &     + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw
+!!gm   #else
+!!gm   ztagnx = ( zi1(ji,jj  ) * wght(ji,jj,2,2) + zi1(ji-1,jj  ) * wght(ji,jj,1,2)   &
+!!gm      &     + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtaux(ji,jj)
+!!gm   ztagny = ( zi2(ji,jj  ) * wght(ji,jj,2,2) + zi2(ji-1,jj  ) * wght(ji,jj,1,2)   &
+!!gm      &     + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * gtauy(ji,jj)
+!!gm   #endif
+            !!gm  always stress provided at I-point (ocean F-point)
+            ztagnx = ( zi1(ji,jj  ) * wght(ji,jj,2,2) + zi1(ji-1,jj  ) * wght(ji,jj,1,2)   &
+               &     + zi1(ji,jj-1) * wght(ji,jj,2,1) + zi1(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * utaui_ice(ji,jj)
+            ztagny = ( zi2(ji,jj  ) * wght(ji,jj,2,2) + zi2(ji-1,jj  ) * wght(ji,jj,1,2)   &
+               &     + zi2(ji,jj-1) * wght(ji,jj,2,1) + zi2(ji-1,jj-1) * wght(ji,jj,1,1) ) * zusw * vtaui_ice(ji,jj)
 
             ! Gradient of ice strength
@@ -143,13 +214,9 @@
 
             ! Computation of the velocity field taking into account the ice-ice interaction.                                 
-            ! Terms that are independent of the velocity field.
-            za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * v_io(ji,jj) - zgphsx
-            za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u_io(ji,jj) - zgphsy
+            ! Terms that are independent of the ice velocity field.
+            za1ct(ji,jj) = ztagnx - zcorl(ji,jj) * vi_oce(ji,jj) - zgphsx
+            za2ct(ji,jj) = ztagny + zcorl(ji,jj) * ui_oce(ji,jj) - zgphsy
          END DO
       END DO
-
-!! inutile!!
-!!??    CALL lbc_lnk( za1ct, 'I', -1. )
-!!??    CALL lbc_lnk( za2ct, 'I', -1. )
 
 
@@ -164,470 +231,357 @@
          ! Computation of free drift field for free slip boundary conditions.
 
-           DO jj = k_j1, k_jpj-1
-              DO ji = 1, jpim1
-                 !- Rate of strain tensor.
-                 zt11 =   akappa(ji,jj,1,1) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) - u_ice(ji,jj  ) - u_ice(ji  ,jj+1) )  &
-                    &   + akappa(ji,jj,1,2) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) + v_ice(ji,jj  ) + v_ice(ji  ,jj+1) )
-                 zt12 = - akappa(ji,jj,2,2) * ( u_ice(ji  ,jj) + u_ice(ji+1,jj  ) - u_ice(ji,jj+1) - u_ice(ji+1,jj+1) )  &
-                    &   - akappa(ji,jj,2,1) * ( v_ice(ji  ,jj) + v_ice(ji+1,jj  ) + v_ice(ji,jj+1) + v_ice(ji+1,jj+1) )
-                 zt22 = - akappa(ji,jj,2,2) * ( v_ice(ji  ,jj) + v_ice(ji+1,jj  ) - v_ice(ji,jj+1) - v_ice(ji+1,jj+1) )  &
-                    &   + akappa(ji,jj,2,1) * ( u_ice(ji  ,jj) + u_ice(ji+1,jj  ) + u_ice(ji,jj+1) + u_ice(ji+1,jj+1) )
-                 zt21 =   akappa(ji,jj,1,1) * ( v_ice(ji+1,jj) + v_ice(ji+1,jj+1) - v_ice(ji,jj  ) - v_ice(ji  ,jj+1) )  &
-                    &   - akappa(ji,jj,1,2) * ( u_ice(ji+1,jj) + u_ice(ji+1,jj+1) + u_ice(ji,jj  ) + u_ice(ji  ,jj+1) )
-
-                 !- Rate of strain tensor. 
-                 zdgp = zt11 + zt22
-                 zdgi = zt12 + zt21
-                 ztrace2 = zdgp * zdgp 
-                 zdeter  = zt11 * zt22 - 0.25 * zdgi * zdgi
-
-                 !  Creep limit depends on the size of the grid.
-                 zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2),  creepl)
-
-                 !-  Computation of viscosities.
-                 zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn )
-                 zviseta (ji,jj) = zviszeta(ji,jj) * usecc2
-              END DO
-           END DO
-!!??       CALL lbc_lnk( zviszeta, 'I', -1. )  ! or T point???   semble reellement inutile
-!!??       CALL lbc_lnk( zviseta , 'I', -1. )
-
-
-           !-  Determination of zc1u, zc2u, zc1v and zc2v.
-           DO jj = k_j1+1, k_jpj-1
-              DO ji = 2, jpim1
-                 ze11   =  akappa(ji-1,jj-1,1,1)
-                 ze12   = +akappa(ji-1,jj-1,2,2)
-                 ze22   =  akappa(ji-1,jj-1,2,1)
-                 ze21   = -akappa(ji-1,jj-1,1,2)
-                 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)
-
-                 zdiag = zvis22 * ( ze11 + ze22 )
-                 zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag
-                 zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21
-                 zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag
-                 zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12
-
-                 ze11   = -akappa(ji,jj-1,1,1)
-                 ze12   = +akappa(ji,jj-1,2,2)
-                 ze22   =  akappa(ji,jj-1,2,1)
-                 ze21   = -akappa(ji,jj-1,1,2)
-                 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)
-
-                 zdiag = zvis22 * ( ze11 + ze22 )
-                 zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag
-                 zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21
-                 zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag
-                 zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12
-
-                 ze11   =  akappa(ji-1,jj,1,1)
-                 ze12   = -akappa(ji-1,jj,2,2)
-                 ze22   =  akappa(ji-1,jj,2,1)
-                 ze21   = -akappa(ji-1,jj,1,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 * ( ze11 + ze22 ) 
-                 zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag
-                 zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21
-                 zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag
-                 zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12
-
-                 ze11   = -akappa(ji,jj,1,1)
-                 ze12   = -akappa(ji,jj,2,2)
-                 ze22   =  akappa(ji,jj,2,1)
-                 ze21   = -akappa(ji,jj,1,2)
-                 zvis11 = 2.0 * 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(ji,jj,2,2) =  zvis11 * ze11 + zdiag
-                 zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21
-                 zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag 
-                 zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12
-              END DO
-           END DO
-
-           DO jj = k_j1+1, k_jpj-1
-              DO ji = 2, jpim1
-                 zc1u(ji,jj) =   &
-                    + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)   &
-                    - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)   &
-                    - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)   &
-                    + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)   &
-                    + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)   &
-                    + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)   &
-                    - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)   &
-                    - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2)
-                 
-                 zc2u(ji,jj) =   &
-                    + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)   &
-                    - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)   &
-                    - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)   &
-                    + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)   &
-                    - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)   &
-                    - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)   &
-                    + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)   &
-                    + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2)
-             END DO
-           END DO
-
-           DO jj = k_j1+1, k_jpj-1
-              DO ji = 2, jpim1
-                 !  zc1v , zc2v.
-                 ze11   =  akappa(ji-1,jj-1,1,2)
-                 ze12   = -akappa(ji-1,jj-1,2,1)
-                 ze22   = +akappa(ji-1,jj-1,2,2)
-                 ze21   =  akappa(ji-1,jj-1,1,1)
-                 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)
-
-                 zdiag = zvis22 * ( ze11 + ze22 )
-                 zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag
-                 zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21
-                 zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag
-                 zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12
+!CDIR NOVERRCHK
+         DO jj = k_j1, k_jpj-1
+!CDIR NOVERRCHK
+            DO ji = 1, fs_jpim1
+               !- Rate of strain tensor.
+               zt11 =   akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji,jj  ) - zu_a(ji  ,jj+1) )  &
+                  &   + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji,jj  ) + zv_a(ji  ,jj+1) )
+               zt12 = - akappa(ji,jj,2,2) * ( zu_a(ji  ,jj) + zu_a(ji+1,jj  ) - zu_a(ji,jj+1) - zu_a(ji+1,jj+1) )  &
+                  &   - akappa(ji,jj,2,1) * ( zv_a(ji  ,jj) + zv_a(ji+1,jj  ) + zv_a(ji,jj+1) + zv_a(ji+1,jj+1) )
+               zt22 = - akappa(ji,jj,2,2) * ( zv_a(ji  ,jj) + zv_a(ji+1,jj  ) - zv_a(ji,jj+1) - zv_a(ji+1,jj+1) )  &
+                  &   + akappa(ji,jj,2,1) * ( zu_a(ji  ,jj) + zu_a(ji+1,jj  ) + zu_a(ji,jj+1) + zu_a(ji+1,jj+1) )
+               zt21 =   akappa(ji,jj,1,1) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) - zv_a(ji,jj  ) - 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  ) + zu_a(ji  ,jj+1) )
+
+               !- Rate of strain tensor. 
+               zdgp = zt11 + zt22
+               zdgi = zt12 + zt21
+               ztrace2 = zdgp * zdgp 
+               zdeter  = zt11 * zt22 - 0.25 * zdgi * zdgi
+
+               !  Creep limit depends on the size of the grid.
+               zdelta = MAX( SQRT( ztrace2 + ( ztrace2 - 4.0 * zdeter ) * usecc2 ),  creepl)
+
+               !-  Computation of viscosities.
+               zviszeta(ji,jj) = MAX( zpresh(ji,jj) / zdelta, etamn )
+               zviseta (ji,jj) = zviszeta(ji,jj) * usecc2
+            END DO
+         END DO
+
+         !-  Determination of zc1u, zc2u, zc1v and zc2v.
+         DO jj = k_j1+1, k_jpj-1
+            DO ji = fs_2, fs_jpim1
+               !* 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)
+               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
+               zs12_11 =  zvis12 * akappa(ji-1,jj-1,2,2) - zvis21 * akappa(ji-1,jj-1,1,2)
+               zs21_11 = -zvis12 * akappa(ji-1,jj-1,1,2) + zvis21 * akappa(ji-1,jj-1,2,2)
+               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)
+               zdiag = zvis22 * ( -akappa(ji,jj-1,1,1) + akappa(ji,jj-1,2,1) )
+               zs11_21 = -zvis11 * akappa(ji,jj-1,1,1) + zdiag
+               zs12_21 =  zvis12 * akappa(ji,jj-1,2,2) - zvis21 * akappa(ji,jj-1,1,2)
+               zs22_21 =  zvis11 * akappa(ji,jj-1,2,1) + zdiag
+               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) )
+               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)
+               zs22_12 =  zvis11 * akappa(ji-1,jj,2,1) + zdiag
+               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)
+               zdiag = zvis22 * ( -akappa(ji,jj,1,1) + akappa(ji,jj,2,1) )
+               zs11_22 = -zvis11 * akappa(ji,jj,1,1) + zdiag
+               zs12_22 = -zvis12 * akappa(ji,jj,2,2) - zvis21 * akappa(ji,jj,1,2)
+               zs22_22 =  zvis11 * akappa(ji,jj,2,1) + zdiag
+               zs21_22 = -zvis12 * akappa(ji,jj,1,2) - zvis21 * akappa(ji,jj,2,2)
+
+               zc1u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22   &
+                  &          - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12   &
+                  &          - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11   &
+                  &          + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12   &
+                  &          + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21   &
+                  &          + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22   &
+                  &          - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21   &
+                  &          - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22
+
+               zc2u(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22   &
+                  &          - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12   &
+                  &          - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11   &
+                  &          + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12   &
+                  &          - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21   &
+                  &          - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22   &
+                  &          + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21   &
+                  &          + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22
+
+               !* 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)
+               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
+               zs12_11 = -zvis12 * akappa(ji-1,jj-1,2,1) + zvis21 * akappa(ji-1,jj-1,1,1)
+               zs22_11 =  zvis11 * akappa(ji-1,jj-1,2,2) + zdiag
+               zs21_11 =  zvis12 * akappa(ji-1,jj-1,1,1) - zvis21 * akappa(ji-1,jj-1,2,1)
  
-                 ze11   =  akappa(ji,jj-1,1,2)
-                 ze12   = -akappa(ji,jj-1,2,1)
-                 ze22   = +akappa(ji,jj-1,2,2)
-                 ze21   = -akappa(ji,jj-1,1,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)
-
-                 zdiag = zvis22 * ( ze11 + ze22 )
-                 zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag
-                 zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21
-                 zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag
-                 zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12
-
-                 ze11   =  akappa(ji-1,jj,1,2)
-                 ze12   = -akappa(ji-1,jj,2,1)
-                 ze22   = -akappa(ji-1,jj,2,2)
-                 ze21   =  akappa(ji-1,jj,1,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)
-
-                 zdiag = zvis22 * ( ze11 + ze22 )
-                 zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag
-                 zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21
-                 zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag
-                 zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12
-
-                 ze11   =  akappa(ji,jj,1,2)
-                 ze12   = -akappa(ji,jj,2,1)
-                 ze22   = -akappa(ji,jj,2,2)
-                 ze21   = -akappa(ji,jj,1,1)
-                 zvis11 = 2.0 * 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(ji,jj,2,2) =  zvis11 * ze11 + zdiag
-                 zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21
-                 zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag
-                 zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12
-
-              END DO
-           END DO
-
-           DO jj = k_j1+1, k_jpj-1
-              DO ji = 2, jpim1
-                 zc1v(ji,jj) =   &
-                    + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)   &
-                    - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)   &
-                    - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)   &
-                    + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)   &
-                    + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)   &
-                    + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)   &
-                    - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)   &
-                    - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2)
-                 zc2v(ji,jj) =   &
-                    + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)   &
-                    - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)   &
-                    - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)   &
-                    + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)   &
-                    - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)   &
-                    - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)   &
-                    + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)   &
-                    + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2)
-              END DO
-           END DO
-
-         ! Relaxation.
-           
-iflag:   DO jter = 1 , nbitdr
-
-            !  Store previous drift field.   
-            DO jj = k_j1, k_jpj-1
-               zu_ice(:,jj) = u_ice(:,jj)
-               zv_ice(:,jj) = v_ice(:,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)
+               zdiag = zvis22 * ( akappa(ji,jj-1,1,2) + akappa(ji,jj-1,2,2) )
+               zs11_21 =  zvis11 * akappa(ji,jj-1,1,2) + zdiag
+               zs12_21 = -zvis12 * akappa(ji,jj-1,2,1) - zvis21 * akappa(ji,jj-1,1,1)
+               zs22_21 =  zvis11 * akappa(ji,jj-1,2,2) + zdiag
+               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)
+               zdiag = zvis22 * ( akappa(ji-1,jj,1,2) - akappa(ji-1,jj,2,2) )
+               zs11_12 =  zvis11 * akappa(ji-1,jj,1,2) + zdiag
+               zs12_12 = -zvis12 * akappa(ji-1,jj,2,1) + zvis21 * akappa(ji-1,jj,1,1)
+               zs22_12 = -zvis11 * akappa(ji-1,jj,2,2) + zdiag
+               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)
+               zdiag = zvis22 * ( akappa(ji,jj,1,2) - akappa(ji,jj,2,2) )
+               zs11_22 =  zvis11 * akappa(ji,jj,1,2) + zdiag
+               zs12_22 = -zvis12 * akappa(ji,jj,2,1) - zvis21 * akappa(ji,jj,1,1)
+               zs22_22 = -zvis11 * akappa(ji,jj,2,2) + zdiag
+               zs21_22 = -zvis12 * akappa(ji,jj,1,1) - zvis21 * akappa(ji,jj,2,1)
+
+               zc1v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22   &
+                  &          - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12   &
+                  &          - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11   &
+                  &          + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12   &
+                  &          + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21   &
+                  &          + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22   &
+                  &          - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21   &
+                  &          - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22
+
+               zc2v(ji,jj) = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22   &
+                  &          - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12   &
+                  &          - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11   &
+                  &          + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12   &
+                  &          - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21   &
+                  &          - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22   &
+                  &          + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21   &
+                  &          + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22
             END DO
-
+         END DO
+
+         ! GAUSS-SEIDEL method
+         !                                                      ! ================ !
+iflag:   DO jter = 1 , nbitdr                                   !    Relaxation    !
+            !                                                   ! ================ !
+!CDIR NOVERRCHK
             DO jj = k_j1+1, k_jpj-1
-               zsang  = SIGN( 1.e0, fcor(1,jj) ) * sangvg   ! only the sinus changes its sign with the hemisphere
-               DO ji = 2, jpim1
-                 zur     = u_ice(ji,jj) - u_io(ji,jj)
-                 zvr     = v_ice(ji,jj) - v_io(ji,jj)
-                 zmod    = SQRT( zur * zur + zvr * zvr) * ( 1.0 - zfrld(ji,jj) )
-                 za      = rhoco * zmod
-                 zac     = za * cangvg
-                  zmpzas  = alpha * zcorl(ji,jj) + za * zsang
+!CDIR NOVERRCHK
+               DO ji = fs_2, fs_jpim1
+                  !
+                  ze11 =   akappa(ji,jj-1,1,1) * zu_a(ji+1,jj) + akappa(ji,jj-1,1,2) * zv_a(ji+1,jj)
+                  ze12 = + akappa(ji,jj-1,2,2) * zu_a(ji+1,jj) - akappa(ji,jj-1,2,1) * zv_a(ji+1,jj)
+                  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)
+                  zdiag = zvis22 * ( ze11 + ze22 )
+                  zs11_21 =  zvis11 * ze11 + zdiag
+                  zs12_21 =  zvis12 * ze12 + zvis21 * ze21
+                  zs22_21 =  zvis11 * ze22 + zdiag
+                  zs21_21 =  zvis12 * ze21 + zvis21 * ze12
+
+                  ze11 =   akappa(ji-1,jj,1,1) * ( zu_a(ji  ,jj+1) - zu_a(ji-1,jj+1) )   &
+                     &   + akappa(ji-1,jj,1,2) * ( zv_a(ji  ,jj+1) + zv_a(ji-1,jj+1) )
+                  ze12 = + akappa(ji-1,jj,2,2) * ( zu_a(ji-1,jj+1) + zu_a(ji  ,jj+1) )   &
+                     &   - akappa(ji-1,jj,2,1) * ( zv_a(ji-1,jj+1) + zv_a(ji  ,jj+1) )
+                  ze22 = + akappa(ji-1,jj,2,2) * ( zv_a(ji-1,jj+1) + zv_a(ji  ,jj+1) )   &
+                     &   + akappa(ji-1,jj,2,1) * ( zu_a(ji-1,jj+1) + zu_a(ji  ,jj+1) )
+                  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)
+                  zdiag = zvis22 * ( ze11 + ze22 )
+                  zs11_12 =  zvis11 * ze11 + zdiag
+                  zs12_12 =  zvis12 * ze12 + zvis21 * ze21
+                  zs22_12 =  zvis11 * ze22 + zdiag
+                  zs21_12 =  zvis12 * ze21 + zvis21 * ze12
+
+                  ze11 =   akappa(ji,jj,1,1) * ( zu_a(ji+1,jj) + zu_a(ji+1,jj+1) - zu_a(ji  ,jj+1) )   &
+                     &   + akappa(ji,jj,1,2) * ( zv_a(ji+1,jj) + zv_a(ji+1,jj+1) + zv_a(ji  ,jj+1) )
+                  ze12 = - akappa(ji,jj,2,2) * ( zu_a(ji+1,jj) - zu_a(ji  ,jj+1) - zu_a(ji+1,jj+1) )   &
+                     &   - akappa(ji,jj,2,1) * ( zv_a(ji+1,jj) + zv_a(ji  ,jj+1) + zv_a(ji+1,jj+1) )
+                  ze22 = - akappa(ji,jj,2,2) * ( zv_a(ji+1,jj) - zv_a(ji  ,jj+1) - zv_a(ji+1,jj+1) )   &
+                     &   + akappa(ji,jj,2,1) * ( zu_a(ji+1,jj) + zu_a(ji  ,jj+1) + zu_a(ji+1,jj+1) )
+                  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)
+                  zdiag = zvis22 * ( ze11 + ze22 )
+                  zs11_22 =  zvis11 * ze11 + zdiag
+                  zs12_22 =  zvis12 * ze12 + zvis21 * ze21
+                  zs22_22 =  zvis11 * ze22 + zdiag
+                  zs21_22 =  zvis12 * ze21 + zvis21 * ze12
+
+            ! 2nd part
+                  ze11 =   akappa(ji-1,jj-1,1,1) * ( zu_a(ji  ,jj-1) - zu_a(ji-1,jj-1) - zu_a(ji-1,jj) )   &
+                     &   + akappa(ji-1,jj-1,1,2) * ( zv_a(ji  ,jj-1) + zv_a(ji-1,jj-1) + zv_a(ji-1,jj) )
+                  ze12 = - akappa(ji-1,jj-1,2,2) * ( zu_a(ji-1,jj-1) + zu_a(ji  ,jj-1) - zu_a(ji-1,jj) )   &
+                     &   - akappa(ji-1,jj-1,2,1) * ( zv_a(ji-1,jj-1) + zv_a(ji  ,jj-1) + zv_a(ji-1,jj) )
+                  ze22 = - akappa(ji-1,jj-1,2,2) * ( zv_a(ji-1,jj-1) + zv_a(ji  ,jj-1) - zv_a(ji-1,jj) )   &
+                     &   + akappa(ji-1,jj-1,2,1) * ( zu_a(ji-1,jj-1) + zu_a(ji  ,jj-1) + zu_a(ji-1,jj) )
+                  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)
+                  zdiag = zvis22 * ( ze11 + ze22 )
+                  zs11_11 =  zvis11 * ze11 + zdiag
+                  zs12_11 =  zvis12 * ze12 + zvis21 * ze21
+                  zs22_11 =  zvis11 * ze22 + zdiag
+                  zs21_11 =  zvis12 * ze21 + zvis21 * ze12
+
+                  ze11 =   akappa(ji,jj-1,1,1) * ( zu_a(ji+1,jj-1) - zu_a(ji  ,jj-1) )   &
+                     &   + akappa(ji,jj-1,1,2) * ( zv_a(ji+1,jj-1) + zv_a(ji  ,jj-1) )
+                  ze12 = - akappa(ji,jj-1,2,2) * ( zu_a(ji  ,jj-1) + zu_a(ji+1,jj-1) )   &
+                     &   - akappa(ji,jj-1,2,1) * ( zv_a(ji  ,jj-1) + zv_a(ji+1,jj-1) )
+                  ze22 = - akappa(ji,jj-1,2,2) * ( zv_a(ji  ,jj-1) + zv_a(ji+1,jj-1) )   &
+                     &   + akappa(ji,jj-1,2,1) * ( zu_a(ji  ,jj-1) + zu_a(ji+1,jj-1) )
+                  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)
+                  zdiag = zvis22 * ( ze11 + ze22 )
+                  zs11_21 =  zs11_21 + zvis11 * ze11 + zdiag
+                  zs12_21 =  zs12_21 + zvis12 * ze12 + zvis21 * ze21
+                  zs22_21 =  zs22_21 + zvis11 * ze22 + zdiag
+                  zs21_21 =  zs21_21 + zvis12 * ze21 + zvis21 * ze12
+
+                  ze11 = - akappa(ji-1,jj,1,1) * zu_a(ji-1,jj) + akappa(ji-1,jj,1,2) * zv_a(ji-1,jj)
+                  ze12 = - akappa(ji-1,jj,2,2) * zu_a(ji-1,jj) - akappa(ji-1,jj,2,1) * zv_a(ji-1,jj)
+                  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)
+                  zdiag = zvis22 * ( ze11 + ze22 )
+                  zs11_12 =  zs11_12 + zvis11 * ze11 + zdiag
+                  zs12_12 =  zs12_12 + zvis12 * ze12 + zvis21 * ze21
+                  zs22_12 =  zs22_12 + zvis11 * ze22 + zdiag
+                  zs21_12 =  zs21_12 + zvis12 * ze21 + zvis21 * ze12
+
+                  zd1 = + alambd(ji,jj,2,2,2,1) * zs11_21 + alambd(ji,jj,2,2,2,2) * zs11_22  &
+                     &  - alambd(ji,jj,2,2,1,1) * zs11_11 - alambd(ji,jj,2,2,1,2) * zs11_12  &
+                     &  - alambd(ji,jj,1,1,2,1) * zs12_21 - alambd(ji,jj,1,1,1,1) * zs12_11  &
+                     &  + alambd(ji,jj,1,1,2,2) * zs12_22 + alambd(ji,jj,1,1,1,2) * zs12_12  &
+                     &  + alambd(ji,jj,1,2,1,1) * zs21_11 + alambd(ji,jj,1,2,2,1) * zs21_21  &
+                     &  + alambd(ji,jj,1,2,1,2) * zs21_12 + alambd(ji,jj,1,2,2,2) * zs21_22  &
+                     &  - alambd(ji,jj,2,1,1,1) * zs22_11 - alambd(ji,jj,2,1,2,1) * zs22_21  &
+                     &  - alambd(ji,jj,2,1,1,2) * zs22_12 - alambd(ji,jj,2,1,2,2) * zs22_22
+
+                  zd2 = + alambd(ji,jj,2,2,2,1) * zs21_21 + alambd(ji,jj,2,2,2,2) * zs21_22  &
+                     &  - alambd(ji,jj,2,2,1,1) * zs21_11 - alambd(ji,jj,2,2,1,2) * zs21_12  &
+                     &  - alambd(ji,jj,1,1,2,1) * zs22_21 - alambd(ji,jj,1,1,1,1) * zs22_11  &
+                     &  + alambd(ji,jj,1,1,2,2) * zs22_22 + alambd(ji,jj,1,1,1,2) * zs22_12  &
+                     &  - alambd(ji,jj,1,2,1,1) * zs11_11 - alambd(ji,jj,1,2,2,1) * zs11_21  &
+                     &  - alambd(ji,jj,1,2,1,2) * zs11_12 - alambd(ji,jj,1,2,2,2) * zs11_22  &
+                     &  + alambd(ji,jj,2,1,1,1) * zs12_11 + alambd(ji,jj,2,1,2,1) * zs12_21  &
+                     &  + alambd(ji,jj,2,1,1,2) * zs12_12 + alambd(ji,jj,2,1,2,2) * zs12_22
+
+                  zur     = zu_a(ji,jj) - ui_oce(ji,jj)
+                  zvr     = zv_a(ji,jj) - vi_oce(ji,jj)
+!!!!
+                  zmod    = SQRT( zur*zur + zvr*zvr ) * ( 1.0 - zfrld(ji,jj) )
+                  za      = rhoco * zmod
+!!!!
+!!gm chg resul    za      = rhoco * SQRT( zur*zur + zvr*zvr ) * ( 1.0 - 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)
 
-                  za1(ji,jj) =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   &
-                     &        + za * ( cangvg * u_io(ji,jj) - zsang * v_io(ji,jj) )
-
-                  za2(ji,jj) =  zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj)   &
-                     &        + za * ( cangvg * v_io(ji,jj) + zsang * u_io(ji,jj) )
-
-                  zb1(ji,jj)  = zmassdt + zac - zc1u(ji,jj)
-                  zb2(ji,jj)  = zmpzas  - zc2u(ji,jj)
-                  zc1(ji,jj)  = zmpzas  + zc1v(ji,jj)
-                  zc2(ji,jj)  = zmassdt + zac  - zc2v(ji,jj) 
-                  zdeter      = zc1(ji,jj) * zb2(ji,jj) + zc2(ji,jj) * zb1(ji,jj)
-                  zden(ji,jj) = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) )
+                  za1 =  zmassdt * zu0(ji,jj) + zcorlal * zv0(ji,jj) + za1ct(ji,jj)   &
+                     &        + za * ( cangvg * ui_oce(ji,jj) - zsang(ji,jj) * vi_oce(ji,jj) )
+                  za2 =  zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj)   &
+                     &        + za * ( cangvg * vi_oce(ji,jj) + zsang(ji,jj) * ui_oce(ji,jj) )
+                  zb1    = zmassdt + zac - zc1u(ji,jj)
+                  zb2    = zmpzas        - zc2u(ji,jj)
+                  zc1    = zmpzas        + zc1v(ji,jj)
+                  zc2    = zmassdt + zac - zc2v(ji,jj)
+                  zdeter = zc1 * zb2 + zc2 * zb1
+                  zden   = SIGN( rone, zdeter) / MAX( epsd , ABS( zdeter ) )
+                  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)
+
+                  zu_n(ji,jj) = ( zu_a(ji,jj) + om * ( zunw - zu_a(ji,jj) ) * tmu(ji,jj) ) * zmask
+                  zv_n(ji,jj) = ( zv_a(ji,jj) + om * ( zvnw - zv_a(ji,jj) ) * tmu(ji,jj) ) * zmask
                END DO
             END DO
 
-            ! The computation of ice interaction term is splitted into two parts
-            !-------------------------------------------------------------------------
-
-            ! Terms that do not involve already up-dated velocities.
-         
-            DO jj = k_j1+1, k_jpj-1
-               DO ji = 2, jpim1
-                  iim1 = ji
-                  ijm1 = jj - 1
-                  iip1 = ji + 1
-                  ijp1 = jj
-                  ze11 =   akappa(iim1,ijm1,1,1) * u_ice(iip1,ijp1) + akappa(iim1,ijm1,1,2) * v_ice(iip1,ijp1)
-                  ze12 = + akappa(iim1,ijm1,2,2) * u_ice(iip1,ijp1) - akappa(iim1,ijm1,2,1) * v_ice(iip1,ijp1)
-                  ze22 = + akappa(iim1,ijm1,2,2) * v_ice(iip1,ijp1) + akappa(iim1,ijm1,2,1) * u_ice(iip1,ijp1)
-                  ze21 =   akappa(iim1,ijm1,1,1) * v_ice(iip1,ijp1) - akappa(iim1,ijm1,1,2) * u_ice(iip1,ijp1)
-                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
-                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
-                  zvis12 =       zviseta (iim1,ijm1) + dm
-                  zvis21 =       zviseta (iim1,ijm1)
-                  zdiag = zvis22 * ( ze11 + ze22 )
-                  zs11(ji,jj,2,1) =  zvis11 * ze11 + zdiag
-                  zs12(ji,jj,2,1) =  zvis12 * ze12 + zvis21 * ze21
-                  zs22(ji,jj,2,1) =  zvis11 * ze22 + zdiag
-                  zs21(ji,jj,2,1) =  zvis12 * ze21 + zvis21 * ze12
-
-
-                  iim1 = ji - 1
-                  ijm1 = jj
-                  iip1 = ji
-                  ijp1 = jj + 1                   
-                  ze11 =   akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijp1) - u_ice(iim1,ijp1) )   &
-                     &   + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijp1) + v_ice(iim1,ijp1) )
-                  ze12 = + akappa(iim1,ijm1,2,2) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )   &
-                     &   - akappa(iim1,ijm1,2,1) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )
-                  ze22 = + akappa(iim1,ijm1,2,2) * ( v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )   &
-                     &   + akappa(iim1,ijm1,2,1) * ( u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )
-                  ze21 =   akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijp1) - v_ice(iim1,ijp1) )   &
-                     &   - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijp1) + u_ice(iim1,ijp1) )
-                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
-                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
-                  zvis12 =       zviseta (iim1,ijm1) + dm
-                  zvis21 =       zviseta (iim1,ijm1)
-                  zdiag = zvis22 * ( ze11 + ze22 )
-                  zs11(ji,jj,1,2) =  zvis11 * ze11 + zdiag
-                  zs12(ji,jj,1,2) =  zvis12 * ze12 + zvis21 * ze21
-                  zs22(ji,jj,1,2) =  zvis11 * ze22 + zdiag
-                  zs21(ji,jj,1,2) =  zvis12 * ze21 + zvis21 * ze12
-
-                  iim1 = ji
-                  ijm1 = jj
-                  iip1 = ji + 1
-                  ijp1 = jj + 1
-                  ze11 =   akappa(iim1,ijm1,1,1) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) - u_ice(iim1,ijp1) )   &
-                     &   + akappa(iim1,ijm1,1,2) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) + v_ice(iim1,ijp1) )
-                  ze12 = - akappa(iim1,ijm1,2,2) * ( u_ice(iip1,ijm1) - u_ice(iim1,ijp1) - u_ice(iip1,ijp1) )   &
-                     &   - akappa(iim1,ijm1,2,1) * ( v_ice(iip1,ijm1) + v_ice(iim1,ijp1) + v_ice(iip1,ijp1) )
-                  ze22 = - akappa(iim1,ijm1,2,2) * ( v_ice(iip1,ijm1) - v_ice(iim1,ijp1) - v_ice(iip1,ijp1) )   &
-                     &   + akappa(iim1,ijm1,2,1) * ( u_ice(iip1,ijm1) + u_ice(iim1,ijp1) + u_ice(iip1,ijp1) )
-                  ze21 =   akappa(iim1,ijm1,1,1) * ( v_ice(iip1,ijm1) + v_ice(iip1,ijp1) - v_ice(iim1,ijp1) )   &
-                     &   - akappa(iim1,ijm1,1,2) * ( u_ice(iip1,ijm1) + u_ice(iip1,ijp1) + u_ice(iim1,ijp1) ) 
-                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
-                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
-                  zvis12 =       zviseta (iim1,ijm1) + dm
-                  zvis21 =       zviseta (iim1,ijm1)
-
-                  zdiag = zvis22 * ( ze11 + ze22 )
-                  zs11(ji,jj,2,2) =  zvis11 * ze11 + zdiag
-                  zs12(ji,jj,2,2) =  zvis12 * ze12 + zvis21 * ze21
-                  zs22(ji,jj,2,2) =  zvis11 * ze22 + zdiag
-                  zs21(ji,jj,2,2) =  zvis12 * ze21 + zvis21 * ze12
-
-               END DO
+            CALL lbc_lnk( zu_n(:,1:jpj), 'I', -1. )
+            CALL lbc_lnk( zv_n(:,1:jpj), 'I', -1. )
+
+            ! Test of Convergence
+            DO jj = k_j1+1 , k_jpj-1
+               zresr(:,jj) = MAX( ABS( zu_a(:,jj) - zu_n(:,jj) ) , ABS( zv_a(:,jj) - zv_n(:,jj) ) )
             END DO
-
-            ! Terms involving already up-dated velocities.
-            !-Using the arrays zu_ice and zv_ice in the computation of the terms ze leads to JACOBI's method; 
-            ! Using arrays u and v in the computation of the terms ze leads to GAUSS-SEIDEL method.
-             
-            DO jj = k_j1+1, k_jpj-1
-               DO ji = 2, jpim1
-                  iim1 = ji - 1
-                  ijm1 = jj - 1
-                  iip1 = ji
-                  ijp1 = jj
-                  ze11 =   akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) - zu_ice(iim1,ijp1) )   &
-                     &   + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) + zv_ice(iim1,ijp1) )
-                  ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) - zu_ice(iim1,ijp1) )   &
-                     &   - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) + zv_ice(iim1,ijp1) )
-                  ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) - zv_ice(iim1,ijp1) )   &
-                     &   + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) + zu_ice(iim1,ijp1) )
-                  ze21 =   akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) - zv_ice(iim1,ijp1) )   &
-                     &  -  akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) + zu_ice(iim1,ijp1) )
-                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
-                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
-                  zvis12 =       zviseta (iim1,ijm1) + dm
-                  zvis21 =       zviseta (iim1,ijm1)
-
-                  zdiag = zvis22 * ( ze11 + ze22 )
-                  zs11(ji,jj,1,1) =  zvis11 * ze11 + zdiag
-                  zs12(ji,jj,1,1) =  zvis12 * ze12 + zvis21 * ze21
-                  zs22(ji,jj,1,1) =  zvis11 * ze22 + zdiag
-                  zs21(ji,jj,1,1) =  zvis12 * ze21 + zvis21 * ze12
-
-#if defined key_agrif
-             END DO
-          END DO
-
-          DO jj = k_j1+1, k_jpj-1
-             DO ji = 2, jpim1
-#endif
-
-                  iim1 = ji
-                  ijm1 = jj - 1
-                  iip1 = ji + 1
-                  ze11 =   akappa(iim1,ijm1,1,1) * ( zu_ice(iip1,ijm1) - zu_ice(iim1,ijm1) )   &
-                     &   + akappa(iim1,ijm1,1,2) * ( zv_ice(iip1,ijm1) + zv_ice(iim1,ijm1) )
-                  ze12 = - akappa(iim1,ijm1,2,2) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) )   &
-                     &   - akappa(iim1,ijm1,2,1) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) )
-                  ze22 = - akappa(iim1,ijm1,2,2) * ( zv_ice(iim1,ijm1) + zv_ice(iip1,ijm1) )   &
-                     &   + akappa(iim1,ijm1,2,1) * ( zu_ice(iim1,ijm1) + zu_ice(iip1,ijm1) )
-                  ze21 =   akappa(iim1,ijm1,1,1) * ( zv_ice(iip1,ijm1) - zv_ice(iim1,ijm1) )   &
-                     &   - akappa(iim1,ijm1,1,2) * ( zu_ice(iip1,ijm1) + zu_ice(iim1,ijm1) )
-                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
-                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
-                  zvis12 =       zviseta (iim1,ijm1) + dm
-                  zvis21 =       zviseta (iim1,ijm1)
-
-                  zdiag = zvis22 * ( ze11 + ze22 )
-                  zs11(ji,jj,2,1) =  zs11(ji,jj,2,1) + zvis11 * ze11 + zdiag
-                  zs12(ji,jj,2,1) =  zs12(ji,jj,2,1) + zvis12 * ze12 + zvis21 * ze21
-                  zs22(ji,jj,2,1) =  zs22(ji,jj,2,1) + zvis11 * ze22 + zdiag
-                  zs21(ji,jj,2,1) =  zs21(ji,jj,2,1) + zvis12 * ze21 + zvis21 * ze12
-
-
-                  iim1 = ji - 1
-                  ijm1 = jj 
-                  ze11 = - akappa(iim1,ijm1,1,1) * zu_ice(iim1,ijm1) + akappa(iim1,ijm1,1,2) * zv_ice(iim1,ijm1)
-                  ze12 = - akappa(iim1,ijm1,2,2) * zu_ice(iim1,ijm1) - akappa(iim1,ijm1,2,1) * zv_ice(iim1,ijm1)
-                  ze22 = - akappa(iim1,ijm1,2,2) * zv_ice(iim1,ijm1) + akappa(iim1,ijm1,2,1) * zu_ice(iim1,ijm1)
-                  ze21 = - akappa(iim1,ijm1,1,1) * zv_ice(iim1,ijm1) - akappa(iim1,ijm1,1,2) * zu_ice(iim1,ijm1)
-                  zvis11 = 2.0 * zviseta (iim1,ijm1) + dm
-                  zvis22 =       zviszeta(iim1,ijm1) - zviseta(iim1,ijm1)
-                  zvis12 =       zviseta (iim1,ijm1) + dm
-                  zvis21 =       zviseta (iim1,ijm1)
-
-                  zdiag = zvis22 * ( ze11 + ze22 )
-                  zs11(ji,jj,1,2) =  zs11(ji,jj,1,2) + zvis11 * ze11 + zdiag 
-                  zs12(ji,jj,1,2) =  zs12(ji,jj,1,2) + zvis12 * ze12 + zvis21 * ze21
-                  zs22(ji,jj,1,2) =  zs22(ji,jj,1,2) + zvis11 * ze22 + zdiag
-                  zs21(ji,jj,1,2) =  zs21(ji,jj,1,2) + zvis12 * ze21 + zvis21 * ze12
-
-#if defined key_agrif
-             END DO
-          END DO
-
-          DO jj = k_j1+1, k_jpj-1
-             DO ji = 2, jpim1
-#endif
-                  zd1(ji,jj) =   &
-                     + alambd(ji,jj,2,2,2,1) * zs11(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs11(ji,jj,2,2)  &
-                     - alambd(ji,jj,2,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs11(ji,jj,1,2)  &
-                     - alambd(ji,jj,1,1,2,1) * zs12(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs12(ji,jj,1,1)  &
-                     + alambd(ji,jj,1,1,2,2) * zs12(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs12(ji,jj,1,2)  &
-                     + alambd(ji,jj,1,2,1,1) * zs21(ji,jj,1,1) + alambd(ji,jj,1,2,2,1) * zs21(ji,jj,2,1)  &
-                     + alambd(ji,jj,1,2,1,2) * zs21(ji,jj,1,2) + alambd(ji,jj,1,2,2,2) * zs21(ji,jj,2,2)  &
-                     - alambd(ji,jj,2,1,1,1) * zs22(ji,jj,1,1) - alambd(ji,jj,2,1,2,1) * zs22(ji,jj,2,1)  &
-                     - alambd(ji,jj,2,1,1,2) * zs22(ji,jj,1,2) - alambd(ji,jj,2,1,2,2) * zs22(ji,jj,2,2)
-                  zd2(ji,jj) =   &
-                     + alambd(ji,jj,2,2,2,1) * zs21(ji,jj,2,1) + alambd(ji,jj,2,2,2,2) * zs21(ji,jj,2,2)  &
-                     - alambd(ji,jj,2,2,1,1) * zs21(ji,jj,1,1) - alambd(ji,jj,2,2,1,2) * zs21(ji,jj,1,2)  &
-                     - alambd(ji,jj,1,1,2,1) * zs22(ji,jj,2,1) - alambd(ji,jj,1,1,1,1) * zs22(ji,jj,1,1)  &
-                     + alambd(ji,jj,1,1,2,2) * zs22(ji,jj,2,2) + alambd(ji,jj,1,1,1,2) * zs22(ji,jj,1,2)  &
-                     - alambd(ji,jj,1,2,1,1) * zs11(ji,jj,1,1) - alambd(ji,jj,1,2,2,1) * zs11(ji,jj,2,1)  &
-                     - alambd(ji,jj,1,2,1,2) * zs11(ji,jj,1,2) - alambd(ji,jj,1,2,2,2) * zs11(ji,jj,2,2)  &
-                     + alambd(ji,jj,2,1,1,1) * zs12(ji,jj,1,1) + alambd(ji,jj,2,1,2,1) * zs12(ji,jj,2,1)  &
-                     + alambd(ji,jj,2,1,1,2) * zs12(ji,jj,1,2) + alambd(ji,jj,2,1,2,2) * zs12(ji,jj,2,2)
-               END DO
+            zresm = MAXVAL( zresr(1:jpi,k_j1+1:k_jpj-1) )
+!!!! this should be faster on scalar processor
+!           zresm = MAXVAL(  MAX( ABS( zu_a(1:jpi,k_j1+1:k_jpj-1) - zu_n(1:jpi,k_j1+1:k_jpj-1) ),   &
+!              &                  ABS( zv_a(1:jpi,k_j1+1:k_jpj-1) - zv_n(1:jpi,k_j1+1:k_jpj-1) ) )  )
+!!!!
+            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
+
+            DO jj = k_j1, k_jpj
+               zu_a(:,jj) = zu_n(:,jj)
+               zv_a(:,jj) = zv_n(:,jj)
             END DO
 
-            DO jj = k_j1+1, k_jpj-1
-               DO ji = 2, jpim1
-                  zunw = (  ( za1(ji,jj) + zd1(ji,jj) ) * zc2(ji,jj)        &
-                     &    + ( za2(ji,jj) + zd2(ji,jj) ) * zc1(ji,jj) ) * zden(ji,jj)
-
-                  zvnw = (  ( za2(ji,jj) + zd2(ji,jj) ) * zb1(ji,jj)        &
-                     &    - ( za1(ji,jj) + zd1(ji,jj) ) * zb2(ji,jj) ) * zden(ji,jj)
-
-                  zmask = ( 1.0 - MAX( rzero, SIGN( rone , 1.0 - zmass(ji,jj) ) ) ) * tmu(ji,jj)
-
-                  u_ice(ji,jj) = ( u_ice(ji,jj) + om * ( zunw - u_ice(ji,jj) ) * tmu(ji,jj) ) * zmask
-                  v_ice(ji,jj) = ( v_ice(ji,jj) + om * ( zvnw - v_ice(ji,jj) ) * tmu(ji,jj) ) * zmask
-               END DO
+            IF( zresm <= resl )   EXIT   iflag
+
+            !                                                   ! ================ !
+         END DO    iflag                                        !  end Relaxation  !
+         !                                                      ! ================ !
+
+         IF( zindu == 0 ) THEN      ! even iteration
+            DO jj = k_j1 , k_jpj-1
+               zu0(:,jj) = zu_a(:,jj)
+               zv0(:,jj) = zv_a(:,jj)
             END DO
-
-            CALL lbc_lnk( u_ice, 'I', -1. )
-            CALL lbc_lnk( v_ice, 'I', -1. )
-
-            !---  5.2.5.4. Convergence test.
-            DO jj = k_j1+1 , k_jpj-1
-               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ) , ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
-            END DO
-            zresm = MAXVAL( zresr( 1:jpi , k_j1+1:k_jpj-1 ) )
-            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
-
-            IF ( zresm <= resl) EXIT iflag
-
-         END DO iflag
-
-         zindu1 = 1.0 - zindu
-         DO jj = k_j1 , k_jpj-1
-            zu0(:,jj) = zindu * zu0(:,jj) + zindu1 * u_ice(:,jj)
-            zv0(:,jj) = zindu * zv0(:,jj) + zindu1 * v_ice(:,jj)
-         END DO
-      !                                                   ! ==================== !
+         ENDIF
+         !                                                ! ==================== !
       END DO                                              !  end loop over iter  !
       !                                                   ! ==================== !
+
+      ui_ice(:,:) = zu_a(:,1:jpj)
+      vi_ice(:,:) = zv_a(:,1:jpj)
 
       IF(ln_ctl) THEN
          WRITE(charout,FMT="('lim_rhg  : res =',D23.16, ' iter =',I4)") zresm, jter
          CALL prt_ctl_info(charout)
-         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' lim_rhg  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
+         CALL prt_ctl(tab2d_1=ui_ice, clinfo1=' lim_rhg  : ui_ice :', tab2d_2=vi_ice, clinfo2=' vi_ice :')
       ENDIF
 
Index: trunk/NEMO/LIM_SRC/limrst.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limrst.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/limrst.F90	(revision 717)
@@ -17,6 +17,6 @@
    !!----------------------------------------------------------------------
    USE ice
-   USE dom_oce
-   USE ice_oce         ! ice variables
+   USE sbc_oce
+   USE sbc_ice
    USE daymod
 
@@ -28,6 +28,6 @@
 
    PUBLIC   lim_rst_opn     ! routine called by ??? module
-   PUBLIC   lim_rst_write   ! routine called by ??? module
-   PUBLIC   lim_rst_read    ! routine called by ??? module
+   PUBLIC   lim_rst_write   ! routine called by lim_step.F90
+   PUBLIC   lim_rst_read    ! routine called by lim_init.F90
 
    LOGICAL, PUBLIC ::   lrst_ice         !: logical to control the ice restart write 
@@ -56,6 +56,6 @@
       IF( kt == nit000 )   lrst_ice = .FALSE.
       
-      IF( kt == nitrst - 2*nfice + 1 .OR.  nitend - nit000 + 1 <= nfice ) THEN
-         ! beware if model runs less than nfice + 1 time step
+      IF( kt == nitrst - 2*nn_fsbc + 1 .OR.  nitend - nit000 + 1 <= nn_fsbc ) THEN
+         ! beware if model runs less than nn_fsbc + 1 time step
          ! beware of the format used to write kt (default is i8.8, that should be large enough)
          IF( nitrst > 1.0e9 ) THEN   
@@ -74,4 +74,5 @@
    END SUBROUTINE lim_rst_opn
 
+
    SUBROUTINE lim_rst_write( kt )
       !!----------------------------------------------------------------------
@@ -80,22 +81,21 @@
       !! ** purpose  :   output of sea-ice variable in a netcdf file
       !!----------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt     ! number of iteration
-      !!
-      INTEGER ::   iter     ! kt + nfice -1
-      !!----------------------------------------------------------------------
-
-      iter = kt + nfice -1
+      INTEGER, INTENT(in) ::   kt   ! number of iteration
+      !
+      INTEGER ::   iter   ! kt + nn_fsbc -1
+      !!----------------------------------------------------------------------
+
+      iter = kt + nn_fsbc - 1
 
       IF( iter == nitrst ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'lim_rst_write : write ice restart.output NetCDF file  kt =', kt
-         IF(lwp) WRITE(numout,*) '~~~~~~~'         
-      ENDIF
-
-      ! Write in numriw (if iter == nitrst)
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~'
+      ENDIF
+
+      ! Write in numriw (if iter  == nitrst)
       ! ------------------ 
       !                                                                     ! calendar control
-      CALL iom_rstput( iter, nitrst, numriw, 'nfice' , REAL( nfice, wp) )      ! time-step 
-      CALL iom_rstput( iter, nitrst, numriw, 'kt_ice', REAL( iter, wp) )      ! date
+      CALL iom_rstput( iter, nitrst, numriw, 'kt_ice', REAL( iter, wp) ) 
       
       CALL iom_rstput( iter, nitrst, numriw, 'hicif' , hicif (:,:)   )      ! prognostic variables 
@@ -109,8 +109,6 @@
       CALL iom_rstput( iter, nitrst, numriw, 'tbif2' , tbif  (:,:,2) )
       CALL iom_rstput( iter, nitrst, numriw, 'tbif3' , tbif  (:,:,3) )
-      CALL iom_rstput( iter, nitrst, numriw, 'u_ice' , u_ice (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'v_ice' , v_ice (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'gtaux' , gtaux (:,:)   )
-      CALL iom_rstput( iter, nitrst, numriw, 'gtauy' , gtauy (:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'ui_ice', ui_ice(:,:)   )
+      CALL iom_rstput( iter, nitrst, numriw, 'vi_ice', vi_ice(:,:)   )
       CALL iom_rstput( iter, nitrst, numriw, 'qstoif', qstoif(:,:)   )
       CALL iom_rstput( iter, nitrst, numriw, 'fsbbq' , fsbbq (:,:)   )
@@ -165,6 +163,5 @@
       !! ** purpose  :   read of sea-ice variable restart in a netcdf file
       !!----------------------------------------------------------------------
-      !
-      REAL(wp) ::   zfice, ziter
+      REAL(wp) ::   ziter
       !!----------------------------------------------------------------------
 
@@ -172,12 +169,11 @@
          WRITE(numout,*)
          WRITE(numout,*) 'lim_rst_read : read ice NetCDF restart file'
-         WRITE(numout,*) '~~~~~~~~'
+         WRITE(numout,*) '~~~~~~~~~~~~'
       ENDIF
 
       CALL iom_open ( 'restart_ice_in', numrir, kiolib = jprstlib )
 
-      CALL iom_get( numrir, 'nfice' , zfice )
-      CALL iom_get( numrir, 'kt_ice', ziter )    
-      IF(lwp) WRITE(numout,*) '   read ice restart file at time step    : ', ziter
+      CALL iom_get( numrir, 'kt_ice' , ziter )
+      IF(lwp) WRITE(numout,*) '   read ice restart file at time step    : ', INT( ziter )
       IF(lwp) WRITE(numout,*) '   in any case we force it to nit000 - 1 : ', nit000 - 1
 
@@ -186,8 +182,4 @@
       IF( ( nit000 - INT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 )   &
          &     CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nit000 in ice restart',  &
-         &                   '   verify the file or rerun with the value 0 for the',        &
-         &                   '   control of time parameter  nrstdt' )
-      IF( INT(zfice) /= nfice          .AND. ABS( nrstdt ) == 1 )   &
-         &     CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nfice in ice restart',  &
          &                   '   verify the file or rerun with the value 0 for the',        &
          &                   '   control of time parameter  nrstdt' )
@@ -203,8 +195,6 @@
       CALL iom_get( numrir, jpdom_autoglo, 'tbif2' , tbif(:,:,2) )    
       CALL iom_get( numrir, jpdom_autoglo, 'tbif3' , tbif(:,:,3) )    
-      CALL iom_get( numrir, jpdom_autoglo, 'u_ice' , u_ice  )    
-      CALL iom_get( numrir, jpdom_autoglo, 'v_ice' , v_ice  )    
-      CALL iom_get( numrir, jpdom_autoglo, 'gtaux' , gtaux  )    
-      CALL iom_get( numrir, jpdom_autoglo, 'gtauy' , gtauy  )    
+      CALL iom_get( numrir, jpdom_autoglo, 'ui_ice', ui_ice )    
+      CALL iom_get( numrir, jpdom_autoglo, 'vi_ice', vi_ice )    
       CALL iom_get( numrir, jpdom_autoglo, 'qstoif', qstoif )    
       CALL iom_get( numrir, jpdom_autoglo, 'fsbbq' , fsbbq  )    
Index: trunk/NEMO/LIM_SRC/limsbc.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limsbc.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/limsbc.F90	(revision 717)
@@ -61,8 +61,10 @@
       !!              - Update 
       !!     
-      !! ** Outputs : - fsolar  : solar heat flux at sea ice/ocean interface
-      !!              - fnsolar : non solar heat flux 
-      !!              - fsalt   : salt flux at sea ice/ocean interface
-      !!              - fmass   : freshwater flux at sea ice/ocean interface
+      !! ** Outputs : - qsr    : sea heat flux:     solar 
+      !!              - qns    : sea heat flux: non solar
+      !!              - 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.
@@ -183,12 +185,13 @@
             DO ji = 1, jpi
                ! ... change the cosinus angle sign in the south hemisphere
-               zsang  = SIGN(1.e0, gphif(ji-1,jj-1) ) * sangvg 
+!!              zsang  = SIGN(1.e0, gphif(ji-1,jj-1) ) * sangvg    ! do the full loop and avoid lbc_lnk
+               zsang  = SIGN(1.e0, gphif(ji,jj) ) * sangvg
                ! ... 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 )
+               zu_io  = ui_ice(ji,jj) - ui_oce(ji,jj)
+               zv_io  = vi_ice(ji,jj) - vi_oce(ji,jj)
+               zmod   = rhoco * SQRT( zu_io * zu_io + zv_io * zv_io )
                ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point)
-               ztio_u(ji,jj) = cw * zmod * ( cangvg * zu_io - zsang * zv_io )
-               ztio_v(ji,jj) = cw * zmod * ( cangvg * zv_io + zsang * zu_io )
+               ztio_u(ji,jj) = zmod * ( cangvg * zu_io - zsang * zv_io )
+               ztio_v(ji,jj) = zmod * ( cangvg * zv_io + zsang * zu_io )
                ! 
             END DO
@@ -196,5 +199,5 @@
 
          DO jj = 2, jpjm1
-            DO ji = 2, jpim1
+            DO ji = fs_2, fs_jpim1   ! vertor opt.
                ! ... ice-cover wheighted ice-ocean stress at U and V-points  (from I-point values)
                zutau  = 0.5 * ( ztio_u(ji+1,jj) + ztio_u(ji+1,jj+1) )
Index: trunk/NEMO/LIM_SRC/limthd.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limthd.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/limthd.F90	(revision 717)
@@ -4,4 +4,8 @@
    !!              LIM thermo ice model : ice thermodynamic
    !!======================================================================
+   !! History :  1.0  !  00-01 (LIM)
+   !!            2.0  !  02-07 (C. Ethe, G. Madec) F90
+   !!            2.0  !  03-08 (C. Ethe)  add lim_thd_init
+   !!---------------------------------------------------------------------
 #if defined key_ice_lim
    !!----------------------------------------------------------------------
@@ -11,33 +15,30 @@
    !!   lim_thd_init : initialisation of sea-ice thermodynamic
    !!----------------------------------------------------------------------
-   !! * Modules used
+   USE phycst          ! physical constants
    USE dom_oce         ! ocean space and time domain variables
-   USE sbc_oce         ! surface boundary condition: ocean
+   USE lbclnk
+   USE in_out_manager  ! I/O manager
+   USE ice             ! LIM sea-ice variables
    USE ice_oce         ! sea-ice/ocean variables
+   USE sbc_oce         ! 
+   USE sbc_ice         ! 
    USE thd_ice         ! LIM thermodynamic sea-ice variables
    USE dom_ice         ! LIM sea-ice domain
-   USE ice             ! LIM sea-ice variables
    USE iceini
    USE limthd_zdf
    USE limthd_lac
    USE limtab
-   USE phycst          ! physical constants
-   USE in_out_manager  ! I/O manager
    USE prtctl          ! Print control
-   USE lbclnk
       
    IMPLICIT NONE
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_thd       ! called by lim_step
-
-   !! * Module variables
-   REAL(wp)  ::            &  ! constant values
-      epsi20 = 1.e-20   ,  &
-      epsi16 = 1.e-16   ,  &
-      epsi04 = 1.e-04   ,  &
-      zzero  = 0.e0     ,  &
-      zone   = 1.e0
+   PUBLIC   lim_thd    ! called by lim_step
+
+   REAL(wp)  ::   epsi20 = 1.e-20   ,  &  ! constant values
+      &           epsi16 = 1.e-16   ,  &
+      &           epsi04 = 1.e-04   ,  &
+      &           zzero  = 0.e0     ,  &
+      &           zone   = 1.e0
 
    !! * Substitutions
@@ -47,5 +48,5 @@
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -68,13 +69,8 @@
       !!             - back to the geographic grid
       !!
-      !! ** References :
-      !!       H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
-      !!
-      !! History :
-      !!   1.0  !  00-01 (LIM)
-      !!   2.0  !  02-07 (C. Ethe, G. Madec) F90
+      !! References :   Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
       !!---------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! number of iteration
-
+      !!
       INTEGER  ::   ji, jj,    &   ! dummy loop indices
          nbpb  ,               &   ! nb of icy pts for thermo. cal.
@@ -92,12 +88,10 @@
          zfontn             ,  &   ! heat flux from snow thickness
          zfntlat, zpareff          ! test. the val. of lead heat budget
-      REAL(wp), DIMENSION(jpi,jpj) :: &
-         zhicifp            ,  &   ! ice thickness for outputs
-         zqlbsbq                   ! link with lead energy budget qldif
-      REAL(wp), DIMENSION(jpi,jpj,jpk) :: &
-         zmsk                      ! working array
+      REAL(wp), DIMENSION(jpi,jpj) ::   zhicifp,   &  ! ice thickness for outputs
+         &                              zqlbsbq       ! link with lead energy budget qldif
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zmsk      ! working array
       !!-------------------------------------------------------------------
 
-      IF( kt == nit000  )   CALL lim_thd_init  ! Initialization (first time-step only)
+      IF( kt == nit000 )   CALL lim_thd_init  ! Initialization (first time-step only)
    
       !-------------------------------------------!
@@ -188,11 +182,12 @@
             !  temperature and turbulent mixing (McPhee, 1992)
             zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )  ! friction velocity
-            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_io(ji,jj) - tfu(ji,jj) ) 
+!!gm old    fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_io(ji,jj) - tfu(ji,jj) ) 
+            fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( sst_m(ji,jj) - tfu(ji,jj) ) 
             qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * frld(ji,jj) * rdt_ice
                         
             !  partial computation of the lead energy budget (qldif)
             zfontn         = ( sprecip(ji,jj) / rhosn ) * xlsn  !   energy for melting
-            zfnsol         = qnsr_oce(ji,jj)  !  total non solar flux
-            qldif(ji,jj)   = tms(ji,jj) * ( qsr_oce(ji,jj) * ( 1.0 - thcm(ji,jj) )   &
+            zfnsol         = qns(ji,jj)                         !  total non solar flux over the ocean
+            qldif(ji,jj)   = tms(ji,jj) * ( qsr(ji,jj) * ( 1.0 - thcm(ji,jj) )   &
                &                               + zfnsol + fdtcn(ji,jj) - zfontn     &
                &                               + ( 1.0 - zindb ) * fsbbq(ji,jj) )   &
@@ -206,5 +201,6 @@
             
             !  energy needed to bring ocean surface layer until its freezing
-            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_io(ji,jj) ) * ( 1 - zinda )
+!!gm old    qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_io(ji,jj) ) * ( 1 - zinda )
+            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda )
             
             !  calculate oceanic heat flux.
@@ -258,5 +254,5 @@
          CALL tab_2d_1d( nbpb, fr1_i0_1d  (1:nbpb)     , fr1_i0     , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d( nbpb, fr2_i0_1d  (1:nbpb)     , fr2_i0     , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d( nbpb, qnsr_ice_1d(1:nbpb)     , qnsr_ice   , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d( nbpb, qns_ice_1d (1:nbpb)     , qns_ice    , jpi, jpj, npb(1:nbpb) )
 #if ! defined key_coupled
          CALL tab_2d_1d( nbpb, qla_ice_1d (1:nbpb)     , qla_ice    , jpi, jpj, npb(1:nbpb) )
@@ -394,19 +390,17 @@
       IF(ln_ctl) THEN
          CALL prt_ctl_info(' lim_thd  end  ')
-         CALL prt_ctl(tab2d_1=hicif , clinfo1=' lim_thd: hicif   : ', tab2d_2=hsnif , clinfo2=' hsnif  : ')
-         CALL prt_ctl(tab2d_1=frld  , clinfo1=' lim_thd: frld    : ', tab2d_2=hicifp, clinfo2=' hicifp : ')
-         CALL prt_ctl(tab2d_1=phicif, clinfo1=' lim_thd: phicif  : ', tab2d_2=pfrld , clinfo2=' pfrld  : ')
-         CALL prt_ctl(tab2d_1=sist  , clinfo1=' lim_thd: sist    : ')
-         CALL prt_ctl(tab2d_1=tbif(:,:,1), clinfo1=' lim_thd: tbif 1  : ')
-         CALL prt_ctl(tab2d_1=tbif(:,:,2), clinfo1=' lim_thd: tbif 2  : ')
-         CALL prt_ctl(tab2d_1=tbif(:,:,3), clinfo1=' lim_thd: tbif 3  : ')
-         CALL prt_ctl(tab2d_1=fdtcn , clinfo1=' lim_thd: fdtcn   : ', tab2d_2=qdtcn , clinfo2=' qdtcn  : ')
-         CALL prt_ctl(tab2d_1=qstoif, clinfo1=' lim_thd: qstoif  : ', tab2d_2=fsbbq , clinfo2=' fsbbq  : ')
-      ENDIF
-
-    END SUBROUTINE lim_thd
-
-
-    SUBROUTINE lim_thd_init
+         CALL prt_ctl(tab2d_1=hicif      , clinfo1=' hicif   : ', tab2d_2=hsnif      , clinfo2=' hsnif  : ')
+         CALL prt_ctl(tab2d_1=frld       , clinfo1=' frld    : ', tab2d_2=hicifp     , clinfo2=' hicifp : ')
+         CALL prt_ctl(tab2d_1=phicif     , clinfo1=' phicif  : ', tab2d_2=pfrld      , clinfo2=' pfrld  : ')
+         CALL prt_ctl(tab2d_1=sist       , clinfo1=' sist    : ', tab2d_2=tbif(:,:,1), clinfo2=' tbif 1 : ')
+         CALL prt_ctl(tab2d_1=tbif(:,:,2), clinfo1=' tbif 2  : ', tab2d_2=tbif(:,:,3), clinfo2=' tbif 3 : ')
+         CALL prt_ctl(tab2d_1=fdtcn      , clinfo1=' fdtcn   : ', tab2d_2=qdtcn      , clinfo2=' qdtcn  : ')
+         CALL prt_ctl(tab2d_1=qstoif     , clinfo1=' qstoif  : ', tab2d_2=fsbbq      , clinfo2=' fsbbq  : ')
+      ENDIF
+       !
+   END SUBROUTINE lim_thd
+
+
+   SUBROUTINE lim_thd_init
       !!-------------------------------------------------------------------
       !!                   ***  ROUTINE lim_thd_init *** 
@@ -419,7 +413,4 @@
       !!
       !! ** input   :   Namelist namicether
-      !!
-      !! history :
-      !!  8.5  ! 03-08 (C. Ethe) original code
       !!-------------------------------------------------------------------
       NAMELIST/namicethd/ hmelt , hiccrit, hicmin, hiclim, amax  ,        &
Index: trunk/NEMO/LIM_SRC/limthd_zdf.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limthd_zdf.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/limthd_zdf.F90	(revision 717)
@@ -4,11 +4,14 @@
    !!                thermodynamic growth and decay of the ice 
    !!======================================================================
+   !! History :  1.0  !  01-04 (LIM) Original code
+   !!            2.0  !  02-08 (C. Ethe, G. Madec) F90
+   !!----------------------------------------------------------------------
 #if defined key_ice_lim
    !!----------------------------------------------------------------------
    !!   'key_ice_lim'                                     LIM sea-ice model
    !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
    !!   lim_thd_zdf  : vertical accr./abl. and lateral ablation of sea ice
    !!----------------------------------------------------------------------
-   !! * Modules used
    USE par_oce          ! ocean parameters
    USE phycst           ! ???
@@ -22,18 +25,16 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_thd_zdf        ! called by lim_thd
-
-   !! * Module variables
-   REAL(wp)  ::           &  ! constant values
-      epsi20 = 1.e-20  ,  &
-      epsi13 = 1.e-13  ,  &
-      zzero  = 0.e0    ,  &
-      zone   = 1.e0
+   PUBLIC   lim_thd_zdf        ! called by lim_thd
+
+   REAL(wp) ::   epsi20 = 1.e-20  ,  &  ! constant values
+      &          epsi13 = 1.e-13  ,  &
+      &          zzero  = 0.e0    ,  &
+      &          zone   = 1.e0
    !!----------------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
    !! $Id$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
 CONTAINS
 
@@ -64,23 +65,12 @@
       !!              - Performs lateral ablation
       !!
-      !! References :
-      !!   Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
-      !!   Fichefet T. and M. Maqueda 1999, Clim. Dyn, 15(4), 251-268  
+      !! References : Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646   
+      !!              Fichefet T. and M. Maqueda 1999, Clim. Dyn, 15(4), 251-268  
+      !!------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kideb    ! Start point on which the  the computation is applied
+      INTEGER, INTENT(in) ::   kiut     ! End point on which the  the computation is applied
       !!
-      !! History :
-      !!   original    : 01-04 (LIM)
-      !!   addition    : 02-08 (C. Ethe, G. Madec)
-      !!------------------------------------------------------------------
-      !! * Arguments
-      INTEGER , INTENT (in) ::  &
-         kideb ,  &  ! Start point on which the  the computation is applied
-         kiut        ! End point on which the  the computation is applied
-
-      !! * Local variables
       INTEGER ::   ji       ! dummy loop indices
-
-      REAL(wp) , DIMENSION(jpij,2) ::  &
-         zqcmlt        ! energy due to surface( /1 ) and bottom melting( /2 )
-
+      REAL(wp), DIMENSION(jpij,2) ::   zqcmlt        ! energy due to surface( /1 ) and bottom melting( /2 )
       REAL(wp), DIMENSION(jpij) ::  &
          ztsmlt      &    ! snow/ice surface melting temperature
@@ -97,6 +87,5 @@
          , zts_old   &    ! previous surface temperature
          , zidsn , z1midsn , zidsnic ! tempory variables
-
-      REAL(wp), DIMENSION(jpij) :: &
+      REAL(wp), DIMENSION(jpij) ::   &
           zfnet       &  ! net heat flux at the top surface( incl. conductive heat flux)
           , zsprecip  &    ! snow accumulation
@@ -109,5 +98,4 @@
           , zfrld_1d    &    ! new sea/ice fraction
           , zep            ! internal temperature of the 2nd layer of the snow/ice system
-
        REAL(wp), DIMENSION(3) :: & 
           zplediag  &    ! principle diagonal, subdiag. and supdiag. of the 
@@ -115,5 +103,4 @@
           , zsupdiag  &    ! of the temperatures inside the snow-ice system
           , zsmbr          ! second member
-
        REAL(wp) :: & 
           zhsu     &     ! thickness of surface layer
@@ -130,5 +117,4 @@
           , zb2 , zd2 , zb3 , zd3 &
           , ztint          ! equivalent temperature at the snow-ice interface
-
        REAL(wp) :: & 
           zexp      &     ! exponential function of the ice thickness
@@ -148,5 +134,4 @@
           , zdtic        &  ! ice internal temp. increment
           , zqnes          ! conductive energy due to ice melting in the first ice layer
-
        REAL(wp) :: & 
           ztbot     &      ! temperature at the bottom surface
@@ -162,5 +147,4 @@
           , zc1, zpc1, zc2, zpc2, zp1, zp2 & ! tempory variables
           , ztb2, ztb3
-
        REAL(wp) :: & 
           zdrmh         &   ! change in snow/ice thick. after snow-ice formation
@@ -181,5 +165,4 @@
        !     Computation of energies due to surface and bottom melting 
        !-----------------------------------------------------------------------
-       
        
        DO ji = kideb , kiut
@@ -201,5 +184,4 @@
        END DO
 
-
        !-------------------------------------------
        !  2. Calculate some intermediate variables.  
@@ -265,5 +247,4 @@
        !     - qstbif_1d, total energy stored in brine pockets (updating)
        !-------------------------------------------------------------------
-
 
        DO ji = kideb , kiut
@@ -288,5 +269,4 @@
        END DO
 
-
        !--------------------------------------------------------------------------------
        !  4. Computation of the surface temperature : determined by considering the 
@@ -333,6 +313,6 @@
           !---computation of the energy balance function 
           zfts    = - z1mi0 (ji) * qsr_ice_1d(ji)   & ! net absorbed solar radiation
-             &      - qnsr_ice_1d(ji)                & ! total non solar radiation
-             &      - zfcsu (ji)                  ! conductive heat flux from the surface
+             &      - qns_ice_1d(ji)                & ! total non solar radiation
+             &      - zfcsu (ji)                      ! conductive heat flux from the surface
           !---computation of surface temperature increment  
           zdts    = -zfts / zdfts
@@ -360,5 +340,5 @@
           sist_1d(ji) = MIN( ztsmlt(ji) , sist_1d(ji) )
 #if ! defined key_coupled
-          qnsr_ice_1d(ji) = qnsr_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) )
+          qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) )
           qla_ice_1d (ji) = qla_ice_1d (ji) + dqla_ice_1d(ji) * ( sist_1d(ji) - zts_old(ji) )
 #endif
@@ -366,11 +346,9 @@
        END DO
 
-     
-
        !     5.2. Calculate available heat for surface ablation. 
        !---------------------------------------------------------------------
 
        DO ji = kideb, kiut
-          zfnet(ji) = qnsr_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)          
+          zfnet(ji) = qns_ice_1d(ji) + z1mi0(ji) * qsr_ice_1d(ji) + zfcsu(ji)          
           zfnet(ji) = MAX( zzero , zfnet(ji) )
           zfnet(ji) = zfnet(ji) * MAX( zzero , SIGN( zone , sist_1d(ji) - ztsmlt(ji) ) )
@@ -730,18 +708,11 @@
           dvnbq_1d(ji) = ( 1.0 - frld_1d(ji) ) * ( zhicnew - h_ice_1d(ji) )
           dmgwi_1d(ji) = dmgwi_1d(ji) + ( 1.0 -frld_1d(ji) ) * ( h_snow_1d(ji) - zhsnnew ) * rhosn
-          !  case of natural freshwater flux
-#if defined key_lim_fdd  
-          rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhicnew - h_ice_1d(ji) )  * rhoic
+          !---  volume change of ice and snow (used for ocean-ice freshwater flux computation)
+          rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhicnew - h_ice_1d (ji) ) * rhoic
           rdmsnif_1d(ji) = rdmsnif_1d(ji) + ( 1.0 - frld_1d(ji) )   * ( zhsnnew - h_snow_1d(ji) ) * rhosn
 
-#else
-          rdmicif_1d(ji) = rdmicif_1d(ji) + ( 1.0 - frld_1d(ji) ) * (  ( zhicnew - h_ice_1d(ji) ) * rhoic   &
-             &                                                    + ( zhsnnew - h_snow_1d(ji) ) * rhosn )
-#endif
-
           !---  Actualize new snow and ice thickness.
-
           h_snow_1d(ji)  = zhsnnew
-          h_ice_1d(ji)  = zhicnew
+          h_ice_1d (ji)  = zhicnew
 
        END DO
@@ -799,16 +770,18 @@
           qstbif_1d(ji) = zdrfrl2 * qstbif_1d(ji)
           frld_1d(ji)    = zfrld_1d(ji)
-
-       END DO
-       
+          !
+       END DO
+       ! 
     END SUBROUTINE lim_thd_zdf
+
 #else
-   !!======================================================================
-   !!                       ***  MODULE limthd_zdf   ***
-   !!                              no sea ice model  
-   !!======================================================================
+   !!----------------------------------------------------------------------
+   !!   Default Option                                     NO sea-ice model
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE lim_thd_zdf          ! Empty routine
    END SUBROUTINE lim_thd_zdf
 #endif
- END MODULE limthd_zdf
+
+   !!======================================================================
+END MODULE limthd_zdf
Index: trunk/NEMO/LIM_SRC/limtrp.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limtrp.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/limtrp.F90	(revision 717)
@@ -112,6 +112,6 @@
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
-               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 ) )
+               zui_u(ji,jj) = ( ui_ice(ji+1,jj  ) + ui_ice(ji+1,jj+1) ) / ( MAX( tmu(ji+1,jj  ) + tmu(ji+1,jj+1), zvbord ) )
+               zvi_v(ji,jj) = ( vi_ice(ji  ,jj+1) + vi_ice(ji+1,jj+1) ) / ( MAX( tmu(ji  ,jj+1) + tmu(ji+1,jj+1), zvbord ) )
             END DO
          END DO
Index: trunk/NEMO/LIM_SRC/limwri.F90
===================================================================
--- trunk/NEMO/LIM_SRC/limwri.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/limwri.F90	(revision 717)
@@ -15,14 +15,16 @@
    !!   lim_wri_init : initialization and namelist read
    !!----------------------------------------------------------------------
+   USE phycst
    USE dom_oce
+   USE daymod
    USE ice_oce         ! ice variables
+   USE sbc_oce
+   USE sbc_ice
    USE dom_ice
-   USE sbc_oce         ! surface boundary condition: ocean
-   USE phycst
-   USE daymod
    USE ice
+
+   USE lbclnk
+   USE dianam    ! build name of file (routine)
    USE in_out_manager
-   USE dianam    ! build name of file (routine)
-   USE lbclnk
    USE ioipsl
 
@@ -30,5 +32,5 @@
    PRIVATE
 
-   PUBLIC   lim_wri        ! routine called by lim_step.F90
+   PUBLIC   lim_wri        ! routine called by sbcice_lim module
 
    INTEGER, PARAMETER                       ::   jpnoumax = 40   ! maximum number of variable for ice output
@@ -49,7 +51,9 @@
       zone   = 1.e0
 
-   !!----------------------------------------------------------------------
-   !!  LIM 2.0, UCL-LOCEAN-IPSL (2005)
-   !! $Id$
+   !! * Substitutions
+#   include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006)
+   !! $Header: $
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -79,5 +83,5 @@
       !!-------------------------------------------------------------------
       INTEGER, INTENT(in) ::   kt     ! number of iteration
-
+      !!
       INTEGER  ::   ji, jj, jf                      ! dummy loop indices
       CHARACTER(len = 40)  ::   clhstnam, clop
@@ -90,5 +94,5 @@
 
       !                                          !--------------------!
-      IF ( kt == nit000 ) THEN                !   Initialisation   !
+      IF( kt == nit000 ) THEN                    !   Initialisation   !
          !                                       !--------------------!
          CALL lim_wri_init 
@@ -97,5 +101,5 @@
 !!Chris         clop     = "ave(only(x))"      !ibug  namelist parameter a ajouter
          clop     = "ave(x)"
-         zout     = nwrite * rdt_ice / nfice
+         zout     = nwrite * rdt_ice / nn_fsbc
          zsec     = 0.
          niter    = 0
@@ -110,9 +114,9 @@
          
          DO jf = 1, noumef
-            IF ( nc(jf) == 1 )   CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj   &
-                  &                                , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout )
+            IF( nc(jf) == 1 )   CALL histdef( nice, nam(jf), titn(jf), uni(jf), jpi, jpj   &
+               &                                  , nhorid, 1, 1, 1, -99, 32, clop, zsto, zout )
          END DO
          CALL histend( nice )
-         
+         !
       ENDIF
       !                                          !--------------------!
@@ -120,16 +124,9 @@
       !                                          !--------------------!
 
-!!gm  change the print below to have it only at output time step or when nitend =< 100
-      IF(lwp) THEN
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_wri : write ice outputs in NetCDF files at time : ', nyear, nmonth, nday, kt + nfice - 1
-         WRITE(numout,*) '~~~~~~~ '
-      ENDIF
-
-      !-- calculs des valeurs instantanees
+      !-- Store instantaneous values in zcmo
       
       zcmo(:,:, 1:jpnoumax ) = 0.e0 
       DO jj = 2 , jpjm1
-         DO ji = 2 , jpim1
+         DO ji = fs_2 , fs_jpim1
             zindh  = MAX( zzero , SIGN( zone , hicif(ji,jj) * (1.0 - frld(ji,jj) ) - 0.10 ) )
             zinda  = MAX( zzero , SIGN( zone , ( 1.0 - frld(ji,jj) ) - 0.10 ) )
@@ -142,23 +139,27 @@
             zcmo(ji,jj,5)  = sist  (ji,jj)
             zcmo(ji,jj,6)  = fbif  (ji,jj)
-            zcmo(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * tmu(ji,jj  ) + u_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
-                                      + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
+            zcmo(ji,jj,7)  = zindb * (  ui_ice(ji,jj  ) * tmu(ji,jj  ) + ui_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
+                                      + ui_ice(ji,jj+1) * tmu(ji,jj+1) + ui_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
                                   / ztmu 
 
-            zcmo(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * tmu(ji,jj  ) + v_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
-                                      + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
+            zcmo(ji,jj,8)  = zindb * (  vi_ice(ji,jj  ) * tmu(ji,jj  ) + vi_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
+                                      + vi_ice(ji,jj+1) * tmu(ji,jj+1) + vi_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
                                   / ztmu
-            zcmo(ji,jj,9)  = sst_io(ji,jj)
-            zcmo(ji,jj,10) = sss_io(ji,jj)
-
-            zcmo(ji,jj,11) = fnsolar(ji,jj) + fsolar(ji,jj)
-            zcmo(ji,jj,12) = fsolar (ji,jj)
-            zcmo(ji,jj,13) = fnsolar(ji,jj)
+            zcmo(ji,jj,9)  = sst_m(ji,jj)
+            zcmo(ji,jj,10) = sss_m(ji,jj)
+
+!!gm        zcmo(ji,jj,11) = fnsolar(ji,jj) + fsolar(ji,jj)
+!!gm        zcmo(ji,jj,12) = fsolar (ji,jj)
+!!gm        zcmo(ji,jj,13) = fnsolar(ji,jj)
+            zcmo(ji,jj,11) = qns(ji,jj) + qsr(ji,jj)
+            zcmo(ji,jj,12) = qsr(ji,jj)
+            zcmo(ji,jj,13) = qns(ji,jj)
             ! See thersf for the coefficient
-            zcmo(ji,jj,14) = - fsalt(ji,jj) * rday * ( sss_io(ji,jj) + epsi16 ) / soce
-            zcmo(ji,jj,15) = gtaux(ji,jj)
-            zcmo(ji,jj,16) = gtauy(ji,jj)
-            zcmo(ji,jj,17) = ( 1.0 - frld(ji,jj) ) * qsr_ice (ji,jj) + frld(ji,jj) * qsr_oce (ji,jj)
-            zcmo(ji,jj,18) = ( 1.0 - frld(ji,jj) ) * qnsr_ice(ji,jj) + frld(ji,jj) * qnsr_oce(ji,jj)
+!!gm        zcmo(ji,jj,14) = - fsalt(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce
+            zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce    !!gm ???
+            zcmo(ji,jj,15) = utaui_ice(ji,jj)
+            zcmo(ji,jj,16) = vtaui_ice(ji,jj)
+            zcmo(ji,jj,17) = ( 1.0 - frld(ji,jj) ) * qsr_ice(ji,jj) + frld(ji,jj) * qsr(ji,jj)   !!gm a revoir
+            zcmo(ji,jj,18) = ( 1.0 - frld(ji,jj) ) * qns_ice(ji,jj) + frld(ji,jj) * qns(ji,jj)   !!gm a revoir
             zcmo(ji,jj,19) = sprecip(ji,jj)
          END DO
@@ -175,5 +176,6 @@
          END DO
          
-         IF ( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN 
+         IF( jf == 7  .OR. jf == 8  .OR. jf == 11 .OR. jf == 12 .OR. jf == 15 .OR.   &
+            &                            jf == 23 .OR. jf == 24 .OR. jf == 16 ) THEN 
             CALL lbc_lnk( zfield, 'T', -1. )
          ELSE 
@@ -181,11 +183,9 @@
          ENDIF
          
-         IF ( nc(jf) == 1 ) CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 )
+         IF( nc(jf) == 1 )   CALL histwrite( nice, nam(jf), niter, zfield, ndim, ndex51 )
          
       END DO
       
-      IF ( ( nfice * niter + nit000 - 1 ) >= nitend ) THEN
-         CALL histclo( nice ) 
-      ENDIF
+      IF( ( nn_fsbc * niter + nit000 - 1 ) >= nitend )   CALL histclo( nice ) 
       !
    END SUBROUTINE lim_wri
@@ -225,24 +225,18 @@
          field_13, field_14, field_15, field_16, field_17, field_18,   &
          field_19
-!!gm      NAMELIST/namiceout/ noumef, &
-!!           zfield( 1), zfield( 2), zfield( 3), zfield( 4), zfield( 5),   &
-!!           zfield( 6), zfield( 7), zfield( 8), zfield( 9), zfield(10),   &
-!!           zfield(11), zfield(12), zfield(13), zfield(14), zfield(15),   &
-!!gm         zfield(16), zfield(17), zfield(18), zfield(19)
-      !!-------------------------------------------------------------------
-
-      ! Read Namelist namicewri
-      REWIND ( numnam_ice )
+      !!-------------------------------------------------------------------
+
+      REWIND ( numnam_ice )                ! Read Namelist namicewri
       READ   ( numnam_ice  , namiceout )
       
-      zfield(1)  = field_1
-      zfield(2)  = field_2
-      zfield(3)  = field_3
-      zfield(4)  = field_4
-      zfield(5)  = field_5
-      zfield(6)  = field_6
-      zfield(7)  = field_7
-      zfield(8)  = field_8
-      zfield(9)  = field_9
+      zfield( 1) = field_1
+      zfield( 2) = field_2
+      zfield( 3) = field_3
+      zfield( 4) = field_4
+      zfield( 5) = field_5
+      zfield( 6) = field_6
+      zfield( 7) = field_7
+      zfield( 8) = field_8
+      zfield( 9) = field_9
       zfield(10) = field_10
       zfield(11) = field_11
@@ -274,5 +268,5 @@
          DO nf = 1 , noumef         
             WRITE(numout,*) '   ', titn(nf), '   ', nam(nf),'      ', uni(nf),'  ', nc(nf),'        ', cmulti(nf),   &
-               '        ', cadd(nf)
+               &       '        ', cadd(nf)
          END DO
       ENDIF
Index: trunk/NEMO/LIM_SRC/limwri_dimg.h90
===================================================================
--- trunk/NEMO/LIM_SRC/limwri_dimg.h90	(revision 716)
+++ trunk/NEMO/LIM_SRC/limwri_dimg.h90	(revision 717)
@@ -82,5 +82,5 @@
 
        zsto     = rdt_ice
-       zout     = nwrite * rdt_ice / nfice
+       zout     = nwrite * rdt_ice / nn_fsbc
        zsec     = 0.
        niter    = 0
@@ -106,10 +106,10 @@
           zcmo(ji,jj,5)  = sist  (ji,jj)
           zcmo(ji,jj,6)  = fbif  (ji,jj)
-          zcmo(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * tmu(ji,jj  ) + u_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
-               + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
+          zcmo(ji,jj,7)  = zindb * (  ui_ice(ji,jj  ) * tmu(ji,jj  ) + ui_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
+             &                      + ui_ice(ji,jj+1) * tmu(ji,jj+1) + ui_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
                / ztmu 
 
-          zcmo(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * tmu(ji,jj  ) + v_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
-               + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
+          zcmo(ji,jj,8)  = zindb * (  vi_ice(ji,jj  ) * tmu(ji,jj  ) + vi_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
+             &                      + vi_ice(ji,jj+1) * tmu(ji,jj+1) + vi_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
                / ztmu
           zcmo(ji,jj,9)  = sst_io(ji,jj)
@@ -132,8 +132,8 @@
     nmoyice = nmoyice + 1 
     ! compute mean value if it is time to write on file
-    IF ( MOD(kt+nfice-1-nit000+1,nwrite) == 0 ) THEN
+    IF ( MOD(kt+nn_fsbc-1-nit000+1,nwrite) == 0 ) THEN
        rcmoy(:,:,:) = rcmoy(:,:,:) / FLOAT(nmoyice)
 #else  
-       IF ( MOD(kt-nfice-1-nit000+1,nwrite) == 0 ) THEN 
+       IF ( MOD(kt-nn_fsbc-1-nit000+1,nwrite) == 0 ) THEN 
           !  case of instantaneaous output rcmoy(:,:, 1:jpnoumax ) = 0.e0
           DO jj = 2 , jpjm1
@@ -149,10 +149,10 @@
                 rcmoy(ji,jj,5)  = sist  (ji,jj)
                 rcmoy(ji,jj,6)  = fbif  (ji,jj)
-                rcmoy(ji,jj,7)  = zindb * (  u_ice(ji,jj  ) * tmu(ji,jj  ) + u_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
-                     + u_ice(ji,jj+1) * tmu(ji,jj+1) + u_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
+                rcmoy(ji,jj,7)  = zindb * (  ui_ice(ji,jj  ) * tmu(ji,jj  ) + ui_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
+                   &                       + ui_ice(ji,jj+1) * tmu(ji,jj+1) + ui_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
                      / ztmu
 
-                rcmoy(ji,jj,8)  = zindb * (  v_ice(ji,jj  ) * tmu(ji,jj  ) + v_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
-                     + v_ice(ji,jj+1) * tmu(ji,jj+1) + v_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
+                rcmoy(ji,jj,8)  = zindb * (  vi_ice(ji,jj  ) * tmu(ji,jj  ) + vi_ice(ji+1,jj  ) * tmu(ji+1,jj  )   &
+                   &                       + vi_ice(ji,jj+1) * tmu(ji,jj+1) + vi_ice(ji+1,jj+1) * tmu(ji+1,jj+1) ) &
                      / ztmu
                 rcmoy(ji,jj,9)  = sst_io(ji,jj)
@@ -201,5 +201,5 @@
           rcmoy(:,:,:) = 0.0
           nmoyice = 0 
-       END IF     !  MOD(kt+nfice-1-nit000+1, nwrite == 0 ) !
+       END IF     !  MOD(kt+nn_fsbc-1-nit000+1, nwrite == 0 ) !
 
      END SUBROUTINE lim_wri
Index: trunk/NEMO/LIM_SRC/thd_ice.F90
===================================================================
--- trunk/NEMO/LIM_SRC/thd_ice.F90	(revision 716)
+++ trunk/NEMO/LIM_SRC/thd_ice.F90	(revision 717)
@@ -56,5 +56,5 @@
       fr1_i0_1d   ,     &  !:    "                  "      fr1_i0
       fr2_i0_1d   ,     &  !:    "                  "      fr2_i0
-      qnsr_ice_1d ,     &  !:    "                  "      qns_ice
+      qns_ice_1d  ,     &  !:    "                  "      qns_ice
       qfvbq_1d    ,     &  !:    "                  "      qfvbq
       sist_1d     ,     &  !:    "                  "      sist
Index: trunk/NEMO/OPA_SRC/DIA/diawri.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DIA/diawri.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/DIA/diawri.F90	(revision 717)
@@ -16,4 +16,5 @@
    USE sbc_oce         ! surface boundary condition: ocean
    USE sbc_ice         ! surface boundary condition: ice
+   USE sbcssr          ! restoring term toward SST/SSS climatology
    USE phycst          ! physical constants
    USE ocfzpt          ! ocean freezing point
@@ -243,17 +244,17 @@
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
 #endif
-#if ! defined key_dynspg_rl && defined key_ice_lim
-         ! sowaflup = sowaflep + sorunoff + sowafldp + a term associated to
-         !    internal damping to Levitus that can be diagnosed from others
-         ! sowaflcd = sowaflep + sorunoff + sowafldp + iowaflup
-         CALL histdef( nid_T, "iowaflup", "Ice=>ocean net freshwater"          , "kg/m2/s",   &  ! fsalt
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sowaflep", "atmos=>ocean net freshwater"        , "kg/m2/s",   &  ! fmass
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-#endif
+!!$#if ! defined key_dynspg_rl && defined key_ice_lim
+!!$         ! sowaflup = sowaflep + sorunoff + sowafldp + a term associated to
+!!$         !    internal damping to Levitus that can be diagnosed from others
+!!$         ! sowaflcd = sowaflep + sorunoff + sowafldp + iowaflup
+!!$         CALL histdef( nid_T, "iowaflup", "Ice=>ocean net freshwater"          , "kg/m2/s",   &  ! fsalt
+!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+!!$         CALL histdef( nid_T, "sowaflep", "atmos=>ocean net freshwater"        , "kg/m2/s",   &  ! fmass
+!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+!!$#endif
          CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux"              , "Kg/m2/s",   &  ! emp
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
-         CALL histdef( nid_T, "sorunoff", "Runoffs"                            , "Kg/m2/s",   &  ! runoffs
-            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
+!!$         CALL histdef( nid_T, "sorunoff", "Runoffs"                            , "Kg/m2/s",   &  ! runoffs
+!!$            &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
          CALL histdef( nid_T, "sowaflcd", "concentration/dilution water flux"  , "kg/m2/s",   &  ! emps
             &          jpi, jpj, nh_T, 1  , 1, 1  , -99 , 32, clop, zsto, zout )
@@ -421,10 +422,10 @@
       CALL histwrite( nid_T, "sossheig", it, sshn          , ndim_hT, ndex_hT )   ! sea surface height
 #endif
-#if ! defined key_dynspg_rl && defined key_ice_lim
-      CALL histwrite( nid_T, "iowaflup", it, fsalt(:,:)    , ndim_hT, ndex_hT )   ! ice=>ocean water flux
-      CALL histwrite( nid_T, "sowaflep", it, fmass(:,:)    , ndim_hT, ndex_hT )   ! atmos=>ocean water flux
-#endif
+!!$#if ! defined key_dynspg_rl && defined key_ice_lim
+!!$      CALL histwrite( nid_T, "iowaflup", it, fsalt(:,:)    , ndim_hT, ndex_hT )   ! ice=>ocean water flux
+!!$      CALL histwrite( nid_T, "sowaflep", it, fmass(:,:)    , ndim_hT, ndex_hT )   ! atmos=>ocean water flux
+!!$#endif
       CALL histwrite( nid_T, "sowaflup", it, emp           , ndim_hT, ndex_hT )   ! upward water flux
-      CALL histwrite( nid_T, "sorunoff", it, runoff        , ndim_hT, ndex_hT )   ! runoff
+!!$      CALL histwrite( nid_T, "sorunoff", it, runoff        , ndim_hT, ndex_hT )   ! runoff
       CALL histwrite( nid_T, "sowaflcd", it, emps          , ndim_hT, ndex_hT )   ! c/d water flux
       zw2d(:,:) = emps(:,:) * sn(:,:,1) * tmask(:,:,1)
Index: trunk/NEMO/OPA_SRC/DIA/diawri_dimg.h90
===================================================================
--- trunk/NEMO/OPA_SRC/DIA/diawri_dimg.h90	(revision 716)
+++ trunk/NEMO/OPA_SRC/DIA/diawri_dimg.h90	(revision 717)
@@ -76,5 +76,4 @@
     !! * modules used
     USE lib_mpp
-    USE dtasst, ONLY : sst
 
     !! * Arguments
Index: trunk/NEMO/OPA_SRC/DOM/domain.F90
===================================================================
--- trunk/NEMO/OPA_SRC/DOM/domain.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/DOM/domain.F90	(revision 717)
@@ -143,9 +143,9 @@
       NAMELIST/namrun/ no    , cexper   , ln_rstart , nrstdt , nit000,         &
          &             nitend, ndate0   , nleapy   , ninist , nstock,          &
-         &             nwrite, nrunoff  , ln_dimgnnn
+         &             nwrite, ln_dimgnnn
 
       NAMELIST/namdom/ ntopo , e3zps_min, e3zps_rat, ngrid  , nmsh  ,   &
          &             nacc  , atfp     , rdt      , rdtmin , rdtmax,   &
-         &             rdth  , rdtbt    , nfice    , nfbulk , nclosea
+         &             rdth  , rdtbt    , nclosea
       NAMELIST/namcla/ n_cla
       !!----------------------------------------------------------------------
@@ -174,5 +174,4 @@
          WRITE(numout,*) '           frequency of restart file       nstock    = ', nstock
          WRITE(numout,*) '           frequency of output file        nwrite    = ', nwrite
-         WRITE(numout,*) '           runoff option                   nrunoff   = ', nrunoff
          WRITE(numout,*) '           multi file dimgout           ln_dimgnnn   = ', ln_dimgnnn
       ENDIF
@@ -258,16 +257,4 @@
       ENDIF
 
-      IF( lk_ice_lim ) THEN
-         IF(lwp) WRITE(numout,*) '           ice model coupling frequency      nfice  = ', nfice
-         nfbulk = nfice
-         IF( MOD( rday, nfice*rdt ) /= 0 ) THEN 
-            IF(lwp) WRITE(numout,*) ' '
-            IF(lwp) WRITE(numout,*) 'W A R N I N G :  nfice is NOT a multiple of the number of time steps in a day'
-            IF(lwp) WRITE(numout,*) ' '
-         ENDIF
-         IF(lwp) WRITE(numout,*) '           bulk computation frequency       nfbulk  = ', nfbulk, ' = nfice if ice model used'
-         IF(lwp) WRITE(numout,*) '           flag closed sea or not           nclosea = ', nclosea
-      ENDIF
-
       ! Default values
       n_cla = 0
Index: trunk/NEMO/OPA_SRC/LDF/ldfeiv.F90
===================================================================
--- trunk/NEMO/OPA_SRC/LDF/ldfeiv.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/LDF/ldfeiv.F90	(revision 717)
@@ -17,4 +17,5 @@
    USE dom_oce         ! ocean space and time domain
    USE sbc_oce         ! surface boundary condition: ocean
+   USE sbcrnf          ! river runoffs
    USE ldftra_oce      ! ocean tracer   lateral physics
    USE phycst          ! physical constants
@@ -177,5 +178,5 @@
             DO ji = 1, jpi
                zaht      = ( 1. -  MIN( 1., ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min )  &
-                  &      + aht0 * upsrnfh(ji,jj)                          ! enhanced near river mouths
+                  &      + aht0 * rnfmsk(ji,jj)                          ! enhanced near river mouths
                ahtu(ji,jj) = MAX( MAX( zaht_min, aeiu(ji,jj) ) + zaht, aht0 )
                ahtv(ji,jj) = MAX( MAX( zaht_min, aeiv(ji,jj) ) + zaht, aht0 )
@@ -352,5 +353,5 @@
             DO ji = 1, jpi
                zaht      = ( 1. -  MIN( 1., ABS( ff(ji,jj) / zf20 ) ) ) * ( aht0 - zaht_min )  &
-                  &      + aht0 * upsrnfh(ji,jj)                          ! enhanced near river mouths
+                  &      + aht0 * rnfmsk(ji,jj)                          ! enhanced near river mouths
                ahtu(ji,jj) = MAX( MAX( zaht_min, aeiu(ji,jj) ) + zaht, aht0 )
                ahtv(ji,jj) = MAX( MAX( zaht_min, aeiv(ji,jj) ) + zaht, aht0 )
Index: trunk/NEMO/OPA_SRC/SBC/sbc_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/SBC/sbc_oce.F90	(revision 717)
@@ -18,8 +18,6 @@
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qns       !: sea heat flux: non solar                     [W/m2]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qsr       !: sea heat flux:     solar                     [W/m2]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qrp       !: heat flux damping                            [w/m2]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emp       !: freshwater budget: volume flux               [Kg/m2/s]
    REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   emps      !: freshwater budget: concentration/dillution   [Kg/m2/s]
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   erp       !: evaporation damping                          [kg/m2/s]
 
    !!----------------------------------------------------------------------
Index: trunk/NEMO/OPA_SRC/SBC/sbcice_lim.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/SBC/sbcice_lim.F90	(revision 717)
@@ -104,10 +104,10 @@
          DO jj = 2, jpj
             DO ji = fs_2, jpi   ! vector opt.
-               u_oce(ji,jj) = 0.5 * ( ssu_m(ji-1,jj  ) + ssu_m(ji-1,jj-1) )
-               v_oce(ji,jj) = 0.5 * ( ssv_m(ji  ,jj-1) + ssv_m(ji-1,jj-1) )
+               ui_oce(ji,jj) = 0.5 * ( ssu_m(ji-1,jj  ) + ssu_m(ji-1,jj-1) )
+               vi_oce(ji,jj) = 0.5 * ( ssv_m(ji  ,jj-1) + ssv_m(ji-1,jj-1) )
             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( ui_oce, 'I', -1. )   ! I-point (i.e. F-point with ice indices)
+         CALL lbc_lnk( vi_oce, 'I', -1. )   ! I-point (i.e. F-point with ice indices)
 
          ! ... masked sea surface freezing temperature [Kelvin] (set to rt0 over land)
@@ -124,5 +124,5 @@
             CALL blk_ice_clio()
          CASE( 4 )           ! CORE bulk formulation
-            CALL blk_ice_core( sist     , u_ice    , v_ice    , alb_ice_cs,                      &
+            CALL blk_ice_core( sist     , ui_ice   , vi_ice   , alb_ice_cs,                      &
                &                                     utaui_ice, vtaui_ice , qns_ice , qsr_ice,   &
                &                                     qla_ice  , dqns_ice  , dqla_ice,            &
@@ -148,5 +148,5 @@
          ;                              CALL lim_thd      ( kt )      ! Ice thermodynamics 
          ;                              CALL lim_sbc      ( kt )      ! Ice/Ocean Mass & Heat fluxes 
-         IF( MOD( kt+nfice-1, ninfo ) == 0 .OR.   &
+         IF( MOD( kt+nn_fsbc-1, ninfo ) == 0 .OR.   &
             &  ntmoy == 1 )             CALL lim_dia      ( kt )      ! Ice Diagnostics 
          ;                              CALL lim_wri      ( kt )      ! Ice outputs 
Index: trunk/NEMO/OPA_SRC/SBC/sbcssr.F90
===================================================================
--- trunk/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/SBC/sbcssr.F90	(revision 717)
@@ -33,4 +33,7 @@
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sst   ! structure of input SST (file informations, fields read)
    TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   sf_sss   ! structure of input SSS (file informations, fields read)
+
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   erp       !: evaporation damping                          [kg/m2/s]
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   qrp       !: heat flux damping                            [w/m2]
 
    !! * Namelist namsbc_ssr
@@ -99,5 +102,5 @@
             WRITE(numout,*) '                           (Yes=2, volume flux) '
             WRITE(numout,*) '          dQ/dT (restoring magnitude on SST)     dqdt     = ', dqdt, ' W/m2/K'
-            WRITE(numout,*) '          dE/dS (restoring magnitude on SST)     deds     = ', deds, ' ???'
+            WRITE(numout,*) '          dE/dS (restoring magnitude on SST)     deds     = ', deds, ' ...'
          ENDIF
 
Index: trunk/NEMO/OPA_SRC/TRA/traadv_cen2.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv_cen2.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/TRA/traadv_cen2.F90	(revision 717)
@@ -1,11 +1,13 @@
 MODULE traadv_cen2
-   !!==============================================================================
-   !!                       ***  MODULE  traadv_cen2  ***
+   !!======================================================================
+   !!                     ***  MODULE  traadv_cen2  ***
    !! Ocean active tracers:  horizontal & vertical advective trend
-   !!==============================================================================
-   !! History :  8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv 
-   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
-   !!            9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-   !!            " "  !  06-04  (R. Benshila, G. Madec) Step reorganization
+   !!======================================================================
+   !! History :   8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad=traadv 
+   !!             8.5  !  02-06  (G. Madec)  F90: Free form and module
+   !!             9.0  !  04-08  (C. Talandier) New trends organization
+   !!             " "  !  05-11  (V. Garnier) Surface pressure gradient organization
+   !!             " "  !  06-04  (R. Benshila, G. Madec) Step reorganization
+   !!             " "  !  06-07  (G. madec)  add ups_orca_set routine
    !!----------------------------------------------------------------------
 
@@ -14,4 +16,6 @@
    !!                  vertical advection trends using a seconder order
    !!                  centered scheme. (k-j-i loops)
+   !!   ups_orca_set : allow mixed upstream/centered scheme in specific
+   !!                  area (set for orca 2 and 4 only)
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and active tracers
@@ -21,6 +25,8 @@
    USE trdmod_oce      ! ocean variables trends
    USE trdmod          ! ocean active tracers trends 
+   USE closea          ! closed sea
    USE trabbl          ! advective term in the BBL
    USE ocfzpt          !
+   USE sbcrnf          ! river runoffs
    USE in_out_manager  ! I/O manager
    USE lib_mpp
@@ -32,5 +38,9 @@
    PRIVATE
 
-   PUBLIC   tra_adv_cen2   ! routine called by step.F90
+   PUBLIC   tra_adv_cen2    ! routine called by step.F90
+   PUBLIC   ups_orca_set    ! routine used by traadv_cen2_jki.F90
+
+   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   upsmsk    !: mixed upstream/centered scheme near some straits 
+   !                                                   !  and in closed seas (orca 2 and 4 configurations)
 
    REAL(wp), DIMENSION(jpi,jpj) ::   btr2   ! inverse of T-point surface [1/(e1t*e2t)]
@@ -40,5 +50,5 @@
 #  include "vectopt_loop_substitute.h90"
    !!----------------------------------------------------------------------
-   !!   OPA 9.0 , LOCEAN-IPSL (2005) 
+   !!   OPA 9.0 , LOCEAN-IPSL (2006) 
    !! $Id$
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
@@ -119,14 +129,13 @@
       !!
       INTEGER  ::   ji, jj, jk                           ! dummy loop indices
-      REAL(wp) ::                           &
-         zbtr, zta, zsa, zfui, zfvj,        &  ! temporary scalars
-         zhw, ze3tr, zcofi, zcofj,          &  !    "         "
-         zupsut, zupsvt, zupsus, zupsvs,    &  !    "         "
-         zfp_ui, zfp_vj, zfm_ui, zfm_vj,    &  !    "         "
-         zcofk, zupst, zupss, zcent,        &  !    "         "
-         zcens, zfp_w, zfm_w,               &  !    "         "
-         zcenut, zcenvt, zcenus, zcenvs,    &  !    "         "
-         z_hdivn_x, z_hdivn_y, z_hdivn
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, ztrdt, zind   ! 3D workspace 
+      REAL(wp) ::   zta, zsa, zbtr, zhw, ze3tr,       &  ! temporary scalars
+         &          zfp_ui, zfp_vj, zfp_w , zfui  ,   &  !    "         "
+         &          zfm_ui, zfm_vj, zfm_w , zfvj  ,   &  !    "         "
+         &          zcofi , zcofj , zcofk ,           &  !    "         "
+         &          zupsut, zupsus, zcenut, zcenus,   &  !    "         "
+         &          zupsvt, zupsvs, zcenvt, zcenvs,   &  !    "         "
+         &          zupst , zupss , zcent , zcens ,   &  !    "         "
+         &          z_hdivn_x, z_hdivn_y, z_hdivn 
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, ztrdt, zind   ! 3D workspace
       REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zww, ztrds         !  "      "
       !!----------------------------------------------------------------------
@@ -137,6 +146,10 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~   Vector optimization case'
          IF(lwp) WRITE(numout,*)
-         ! 
-         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
+         !
+         upsmsk(:,:) = 0.e0                              ! not upstream by default
+         IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits
+         !                                               ! and in closed seas (orca2 and orca4 only)
+         !   
+         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )        ! inverse of T-point surface
       ENDIF
 
@@ -146,10 +159,10 @@
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zind(ji,jj,jk) =  MAX ( upsrnfh(ji,jj) * upsrnfz(jk),     &  ! changing advection scheme near runoff
-                  &                    upsadv(ji,jj)                     &  ! in the vicinity of some straits
+               zind(ji,jj,jk) = MAX (   &
+                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
+                  upsmsk(ji,jj)                      &  ! some of some straits
 #if defined key_ice_lim
-                  &                  , tmask(ji,jj,jk)                   &  ! half upstream tracer fluxes
-                  &                  * MAX( 0., SIGN( 1., fzptn(ji,jj)   &  ! if tn < ("freezing"+0.1 )
-                  &                                +0.1-tn(ji,jj,jk) ) ) &
+                  !                                     ! below ice covered area (if tn < "freezing"+0.1 )
+                  , MAX(  0., SIGN( 1., fzptn(ji,jj) + 0.1 - tn(ji,jj,jk) )  ) * tmask(ji,jj,jk)   &
 #endif
                   &                  )
@@ -158,7 +171,8 @@
       END DO
 
-
-      !  Horizontal advective fluxes
-      ! -----------------------------
+      ! I. Horizontal advective fluxes
+      ! ------------------------------
+      !  Second order centered tracer flux at u and v-points
+      ! -----------------------------------------------------
       !                                                ! ===============
       DO jk = 1, jpkm1                                 ! Horizontal slab
@@ -208,11 +222,9 @@
                zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
 #endif
-               ! horizontal advective trends
-               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
+               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &    ! horizontal advective trends
                   &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
                zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   &
                   &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  )
-               ! add it to the general tracer trends
-               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
+               ta(ji,jj,jk) = ta(ji,jj,jk) + zta                          ! add it to the general tracer trends
                sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
             END DO
@@ -279,7 +291,5 @@
          &                       tab3d_2=sa, clinfo2=            ' Sa: ', mask2=tmask, clinfo3='tra' )
 
-      ! 4. "zonal" mean advective heat and salt transport 
-      ! -------------------------------------------------
-
+      ! "zonal" mean advective heat and salt transport
       IF( ln_diaptr .AND. ( MOD( kt, nf_ptr ) == 0 ) ) THEN
          IF( lk_zco ) THEN
@@ -313,24 +323,18 @@
       ENDIF
 
-      ! 1. Vertical advective fluxes
+      ! 1. Vertical advective fluxes (Second order centered tracer flux at w-point)
       ! ----------------------------
-      ! Second order centered tracer flux at w-point
       DO jk = 2, jpk
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
-               ! upstream indicator
-               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
-               ! velocity * 1/2
-               zhw = 0.5 * pwn(ji,jj,jk)
-               ! upstream scheme
-               zfp_w = zhw + ABS( zhw )
+               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )         ! upstream indicator
+               zhw = 0.5 * pwn(ji,jj,jk)                               ! velocity * 1/2
+               zfp_w = zhw + ABS( zhw )                                ! upstream scheme
                zfm_w = zhw - ABS( zhw )
                zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1)
                zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1)
-               ! centered scheme
-               zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )
+               zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )         ! centered scheme
                zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) )
-               ! mixed centered / upstream scheme
-               zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
+               zwx(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent      ! mixed centered / upstream scheme
                zwy(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens
             END DO
@@ -344,9 +348,7 @@
             DO ji = fs_2, fs_jpim1   ! vector opt.
                ze3tr = 1. / fse3t(ji,jj,jk)
-               ! vertical advective trends
-               zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )
+               zta = - ze3tr * ( zwx(ji,jj,jk) - zwx(ji,jj,jk+1) )     ! vertical advective trends
                zsa = - ze3tr * ( zwy(ji,jj,jk) - zwy(ji,jj,jk+1) )
-               ! add it to the general tracer trends
-               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
+               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta                      ! add it to the general tracer trends
                sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
             END DO
@@ -387,4 +389,65 @@
       !
    END SUBROUTINE tra_adv_cen2
+   
+   
+   SUBROUTINE ups_orca_set
+      !!----------------------------------------------------------------------
+      !!                  ***  ROUTINE ups_orca_set  ***
+      !!       
+      !! ** Purpose :   add a portion of upstream scheme in area where the
+      !!                centered scheme generates too strong overshoot
+      !!
+      !! ** Method  :   orca (R4 and R2) confiiguration setting. Set upsmsk
+      !!                array to nozero value in some straith. 
+      !!
+      !! ** Action : - upsmsk set to 1 at some strait, 0 elsewhere for orca
+      !!----------------------------------------------------------------------
+      INTEGER  ::   ii0, ii1, ij0, ij1      ! temporary integers
+      !!----------------------------------------------------------------------
+      
+      ! mixed upstream/centered scheme near river mouths
+      ! ------------------------------------------------
+      SELECT CASE ( jp_cfg )
+      !                                        ! =======================
+      CASE ( 4 )                               !  ORCA_R4 configuration 
+         !                                     ! =======================
+         !                                          ! Gibraltar Strait
+         ii0 =  70   ;   ii1 =  71
+         ij0 =  52   ;   ij1 =  53   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
+         !
+         !                                     ! =======================
+      CASE ( 2 )                               !  ORCA_R2 configuration 
+         !                                     ! =======================
+         !                                          ! Gibraltar Strait
+         ij0 = 102   ;   ij1 = 102
+         ii0 = 138   ;   ii1 = 138   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.20
+         ii0 = 139   ;   ii1 = 139   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
+         ii0 = 140   ;   ii1 = 140   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
+         ij0 = 101   ;   ij1 = 102
+         ii0 = 141   ;   ii1 = 141   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
+         !                                          ! Bab el Mandeb Strait
+         ij0 =  87   ;   ij1 =  88
+         ii0 = 164   ;   ii1 = 164   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.10
+         ij0 =  88   ;   ij1 =  88
+         ii0 = 163   ;   ii1 = 163   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
+         ii0 = 162   ;   ii1 = 162   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.40
+         ii0 = 160   ;   ii1 = 161   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
+         ij0 =  89   ;   ij1 =  89
+         ii0 = 158   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
+         ij0 =  90   ;   ij1 =  90
+         ii0 = 160   ;   ii1 = 160   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
+         !                                          ! Sound Strait
+         ij0 = 116   ;   ij1 = 116
+         ii0 = 144   ;   ii1 = 144   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
+         ii0 = 145   ;   ii1 = 147   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.50
+         ii0 = 148   ;   ii1 = 148   ;   upsmsk( mi0(ii0):mi1(ii1) , mj0(ij0):mj1(ij1) ) = 0.25
+         !
+      END SELECT 
+      
+      ! mixed upstream/centered scheme over closed seas
+      ! -----------------------------------------------
+      CALL clo_ups( upsmsk(:,:) )
+      !
+   END SUBROUTINE ups_orca_set
 
    !!======================================================================
Index: trunk/NEMO/OPA_SRC/TRA/traadv_cen2_jki.F90
===================================================================
--- trunk/NEMO/OPA_SRC/TRA/traadv_cen2_jki.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/TRA/traadv_cen2_jki.F90	(revision 717)
@@ -4,15 +4,17 @@
    !! Ocean active tracers:  horizontal & vertical advective trend
    !!==============================================================================
-   !! History :  8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv
-   !!            8.5  !  02-06  (G. Madec)  F90: Free form and module
-   !!            9.0  !  05-11  (V. Garnier) Surface pressure gradient organization
-   !!            " "  !  06-04  (R. Benshila, G. Madec) Step reorganization
+   !! History :   8.2  !  01-08  (G. Madec, E. Durand)  trahad+trazad = traadv
+   !!             8.5  !  02-06  (G. Madec)  F90: Free form and module
+   !!             9.0  !  04-08  (C. Talandier) New trends organization
+   !!             " "  !  05-11  (V. Garnier) Surface pressure gradient organization
+   !!             " "  !  06-04  (R. Benshila, G. Madec) Step reorganization
+   !!             " "  !  06-07  (G. madec)  add ups_orca_set routine
    !!----------------------------------------------------------------------
 
    !!----------------------------------------------------------------------
    !!   tra_adv_cen2_jki : update the tracer trend with the horizontal and
-   !!                  vertical advection trends using a seconder order
-   !!                  centered scheme. Auto-tasking case, k-slab for
-   !!                  hor. adv., j-slab for vert. adv. (j-k-i loops)
+   !!                      vertical advection trends using a seconder order
+   !!                      centered scheme. Auto-tasking case, k-slab for
+   !!                      hor. adv., j-slab for vert. adv. (j-k-i loops)
    !!----------------------------------------------------------------------
    USE oce             ! ocean dynamics and active tracers
@@ -26,14 +28,15 @@
    USE lib_mpp
    USE lbclnk          ! ocean lateral boundary condition (or mpp link)
+   USE sbcrnf          ! river runoffs
    USE in_out_manager  ! I/O manager
    USE diaptr          ! poleward transport diagnostics
    USE prtctl          ! Print control
-
+   
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   tra_adv_cen2_jki   ! routine called by step.F90
-
-   REAL(wp), DIMENSION(jpi,jpj) ::   btr2   ! inverse of T-point surface [1/(e1t*e2t)]
+   PUBLIC   tra_adv_cen2_jki    ! routine called by step.F90
+
+   REAL(wp), DIMENSION(jpi,jpj) ::   btr2    ! inverse of T-point surface [1/(e1t*e2t)]
 
    !! * Substitutions
@@ -87,4 +90,6 @@
       !!         zta+tn*divn=1/(e1t*e2t*e3t) { mi-1( e2u*e3u  un  di[tn] )
       !!                                      +mj-1( e1v*e3v  vn  mj[tn] )  }
+      !!         C A U T I O N : the trend saved is the centered trend only.
+      !!      It does not take into account the upstream part of the scheme.
       !!         NB:in z-coordinate - full step (ln_zco=T) e3u=e3v=e3t, so 
       !!      they vanish from the expression of the flux and divergence.
@@ -109,28 +114,26 @@
       !!
       !! ** Action :  - update (ta,sa) with the now advective tracer trends
-      !!              - save trends in (ztrdt,ztrds) ('key_trdtra')
       !!----------------------------------------------------------------------
-      USE oce, ONLY :   zwx => ua   ! use ua as workspace
-      USE oce, ONLY :   zwy => va   ! use va as workspace
-      !!
-      INTEGER , INTENT(in)                         ::   kt    ! ocean time-step index
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun   ! ocean velocity u-component
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pvn   ! ocean velocity v-component
-      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pwn   ! ocean velocity w-component
-      !!
+      USE oce               ,   zwx => ua          ! use ua as workspace
+      USE oce               ,   zwy => va          ! use va as workspace
+      USE traadv_cen2, ONLY :   ups_orca_set,   &  ! upstream indicator near some straits 
+         &                      upsmsk             ! and over closed sea (orca 2 and 4)
+
+      INTEGER , INTENT(in)                         ::   kt              ! ocean time-step index
+      REAL(wp), INTENT(in), DIMENSION(jpi,jpj,jpk) ::   pun, pvn, pwn   ! now ocean velocity fields
+
       INTEGER  ::   ji, jj, jk                           ! dummy loop indices
-      REAL(wp) ::                           &
-         zbtr, zta, zsa, zfui, zfvj,        &  ! temporary scalars
-         zhw, ze3tr, zcofi, zcofj,          &  !    "         "
-         zupsut, zupsvt, zupsus, zupsvs,    &  !    "         "
-         zfp_ui, zfp_vj, zfm_ui, zfm_vj,    &  !    "         "
-         zcofk, zupst, zupss, zcent,        &  !    "         "
-         zcens, zfp_w, zfm_w,               &  !    "         "
-         zcenut, zcenvt, zcenus, zcenvs,    &  !    "         "
-         z_hdivn_x, z_hdivn_y, z_hdivn
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, ztrdt, zind   ! 3D workspace 
-      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zww, ztrds         !  "      "
+      REAL(wp) ::   zta, zsa, zbtr, zhw, ze3tr,       &  ! temporary scalars
+         &          zfp_ui, zfp_vj, zfp_w , zfui  ,   &  !    "         "
+         &          zfm_ui, zfm_vj, zfm_w , zfvj  ,   &  !    "         "
+         &          zcofi , zcofj , zcofk ,           &  !    "         "
+         &          zupsut, zupsus, zcenut, zcenus,   &  !    "         "
+         &          zupsvt, zupsvs, zcenvt, zcenvs,   &  !    "         "
+         &          zupst , zupss , zcent , zcens ,   &  !    "         "
+         &          z_hdivn_x, z_hdivn_y, z_hdivn 
+      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   zwz, ztrdt, zind,   &  ! temporary workspace arrays
+         &                                  zww, ztrds             !    "         "
       !!----------------------------------------------------------------------
-
+      !
       IF( kt == nit000 ) THEN
          IF(lwp) WRITE(numout,*)
@@ -138,6 +141,10 @@
          IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~   Auto-tasking case'
          IF(lwp) WRITE(numout,*)
-         !
-         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )
+         !         
+         upsmsk(:,:) = 0.e0                              ! not upstream by default
+         IF( cp_cfg == "orca" )   CALL ups_orca_set      ! set mixed Upstream/centered scheme near some straits
+         !                                               ! and in closed seas (orca2 and orca4 only)
+
+         btr2(:,:) = 1. / ( e1t(:,:) * e2t(:,:) )        ! inverse of T-point surface
       ENDIF
 
@@ -147,37 +154,36 @@
       DO jk = 1, jpkm1                                 ! Horizontal slab
          !                                             ! ===============
-         !  Upstream / centered scheme indicator
-         ! --------------------------------------
+
+         ! 0. Upstream / centered scheme indicator
+         ! ---------------------------------------
          DO jj = 1, jpj
             DO ji = 1, jpi
-               zind(ji,jj,jk) =  MAX (              &
-                  upsrnfh(ji,jj) * upsrnfz(jk),     &  ! changing advection scheme near runoff
-                  upsadv(ji,jj)                     &  ! in the vicinity of some straits
+               zind(ji,jj,jk) = MAX (   &
+                  rnfmsk(ji,jj) * rnfmsk_z(jk),      &  ! near runoff mouths (& closed sea outflows)
+                  upsmsk(ji,jj)                      &  ! some of some straits
 #if defined key_ice_lim
-                  , tmask(ji,jj,jk)                 &  ! half upstream tracer fluxes if tn < ("freezing"+0.1 )
-                  * MAX( 0., SIGN( 1., fzptn(ji,jj)+0.1-tn(ji,jj,jk) ) )   &
-#endif
-               )
-            END DO
-         END DO
-
-         !  Horizontal advective fluxes
-         ! -----------------------------
+                  !                                     ! below ice covered area (if tn < "freezing"+0.1 )
+                  , MAX(  0., SIGN( 1., fzptn(ji,jj) + 0.1 - tn(ji,jj,jk) )  ) * tmask(ji,jj,jk)   &
+#endif
+                  &                  )
+            END DO
+         END DO
+
+
+         ! I. Horizontal advective fluxes
+         ! ------------------------------
          ! Second order centered tracer flux at u and v-points
          DO jj = 1, jpjm1
             DO ji = 1, fs_jpim1   ! vector opt.
-               ! upstream indicator
-               zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )
+               zcofi = MAX( zind(ji+1,jj,jk), zind(ji,jj,jk) )               ! upstream indicator
                zcofj = MAX( zind(ji,jj+1,jk), zind(ji,jj,jk) )
-               ! volume fluxes * 1/2
-#if defined key_zco
-               zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)
+#if defined key_zco
+               zfui = 0.5 * e2u(ji,jj) * pun(ji,jj,jk)                       ! volume fluxes * 1/2 (zco)
                zfvj = 0.5 * e1v(ji,jj) * pvn(ji,jj,jk)
 #else
-               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)
+               zfui = 0.5 * e2u(ji,jj) * fse3u(ji,jj,jk) * pun(ji,jj,jk)     ! volume fluxes * 1/2
                zfvj = 0.5 * e1v(ji,jj) * fse3v(ji,jj,jk) * pvn(ji,jj,jk)
 #endif
-               ! upstream scheme
-               zfp_ui = zfui + ABS( zfui )
+               zfp_ui = zfui + ABS( zfui )                                   ! upstream scheme
                zfp_vj = zfvj + ABS( zfvj )
                zfm_ui = zfui - ABS( zfui )
@@ -187,11 +193,9 @@
                zupsus = zfp_ui * sb(ji,jj,jk) + zfm_ui * sb(ji+1,jj  ,jk)
                zupsvs = zfp_vj * sb(ji,jj,jk) + zfm_vj * sb(ji  ,jj+1,jk)
-               ! centered scheme
-               zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) )
+               zcenut = zfui * ( tn(ji,jj,jk) + tn(ji+1,jj  ,jk) )           ! centered scheme
                zcenvt = zfvj * ( tn(ji,jj,jk) + tn(ji  ,jj+1,jk) )
                zcenus = zfui * ( sn(ji,jj,jk) + sn(ji+1,jj  ,jk) )
                zcenvs = zfvj * ( sn(ji,jj,jk) + sn(ji  ,jj+1,jk) )
-               ! mixed centered / upstream scheme
-               zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut
+               zwx(ji,jj,jk) = zcofi * zupsut + (1.-zcofi) * zcenut          ! mixed centered / upstream scheme
                zwy(ji,jj,jk) = zcofj * zupsvt + (1.-zcofj) * zcenvt
                zww(ji,jj,jk) = zcofi * zupsus + (1.-zcofi) * zcenus
@@ -200,19 +204,19 @@
          END DO
 
-         ! Tracer flux divergence at t-point added to the general trend
+         ! 2. Tracer flux divergence at t-point added to the general trend
+         ! ---------------------------------------------------------------
+
          DO jj = 2, jpjm1
             DO ji = fs_2, fs_jpim1   ! vector opt.
 #if defined key_zco
-               zbtr = btr2(ji,jj) 
-#else
-               zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)
-#endif
-               ! horizontal advective trends
-               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &
+               zbtr = btr2(ji,jj)                                           ! inverse of the volume (zco)
+#else
+               zbtr = btr2(ji,jj) / fse3t(ji,jj,jk)                         ! inverse of the volume
+#endif
+               zta = - zbtr * (  zwx(ji,jj,jk) - zwx(ji-1,jj  ,jk)   &      ! horizontal advective trends
                   &            + zwy(ji,jj,jk) - zwy(ji  ,jj-1,jk)  )
                zsa = - zbtr * (  zww(ji,jj,jk) - zww(ji-1,jj  ,jk)   &
                   &            + zwz(ji,jj,jk) - zwz(ji  ,jj-1,jk)  )
-               ! add it to the general tracer trends
-               ta(ji,jj,jk) = ta(ji,jj,jk) + zta
+               ta(ji,jj,jk) = ta(ji,jj,jk) + zta                            ! add it to the general tracer trends
                sa(ji,jj,jk) = sa(ji,jj,jk) + zsa
             END DO
@@ -225,5 +229,5 @@
          &                       tab3d_2=sa, clinfo2=            ' Sa: ', mask2=tmask, clinfo3='tra' )
 
-      ! Save the horizontal advective trends for diagnostics
+      ! 3. Save the horizontal advective trends for diagnostic
       IF( l_trdtra )   THEN
          ztrdt(:,:,:) = 0.e0   ;   ztrds(:,:,:) = 0.e0
@@ -298,6 +302,6 @@
       ENDIF
 
-      !  Vertical advection
-      ! --------------------
+      ! II. Vertical advection
+      ! ----------------------
 !CDIR PARALLEL DO
 !$OMP PARALLEL DO
@@ -312,39 +316,34 @@
          IF( lk_dynspg_rl ) THEN                          ! rigid lid : flux set to zero
             zwz(:,jj, 1 ) = 0.e0    ;    zww(:,jj, 1 ) = 0.e0
-         ELSE                                             ! free surface-constant volume
+         ELSE                                             ! free surface-constant volume : advection across the surface
             zwz(:,jj, 1 ) = pwn(:,jj,1) * tn(:,jj,1)
             zww(:,jj, 1 ) = pwn(:,jj,1) * sn(:,jj,1)
          ENDIF
          
-         ! Vertical advective fluxes at w-point
+         ! 1. Vertical advective fluxes (Second order centered tracer flux at w-point)
+         ! ----------------------------
          DO jk = 2, jpk
             DO ji = 2, jpim1
-               ! upstream indicator
-               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )
-               ! velocity * 1/2
-               zhw = 0.5 * pwn(ji,jj,jk)
-               ! upstream scheme
-               zfp_w = zhw + ABS( zhw )
+               zcofk = MAX( zind(ji,jj,jk-1), zind(ji,jj,jk) )          ! upstream indicator
+               zhw = 0.5 * pwn(ji,jj,jk)                                ! velocity * 1/2
+               zfp_w = zhw + ABS( zhw )                                 ! upstream scheme
                zfm_w = zhw - ABS( zhw )
                zupst = zfp_w * tb(ji,jj,jk) + zfm_w * tb(ji,jj,jk-1)
                zupss = zfp_w * sb(ji,jj,jk) + zfm_w * sb(ji,jj,jk-1)
-               ! centered scheme
-               zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )
+               zcent = zhw * ( tn(ji,jj,jk) + tn(ji,jj,jk-1) )          ! centered scheme
                zcens = zhw * ( sn(ji,jj,jk) + sn(ji,jj,jk-1) )
-               ! mixed centered / upstream scheme
-               zwz(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent
+               zwz(ji,jj,jk) = zcofk * zupst + (1.-zcofk) * zcent       ! mixed centered / upstream scheme
                zww(ji,jj,jk) = zcofk * zupss + (1.-zcofk) * zcens
             END DO
          END DO
 
-         ! Tracer flux divergence at t-point added to the general trend
+         ! 2. Tracer flux divergence at t-point added to the general trend
+         ! -------------------------
          DO jk = 1, jpkm1
             DO ji = 2, jpim1
                ze3tr = 1. / fse3t(ji,jj,jk)
-               ! vertical advective trends
-               zta = - ze3tr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )
+               zta = - ze3tr * ( zwz(ji,jj,jk) - zwz(ji,jj,jk+1) )      ! vertical advective trends
                zsa = - ze3tr * ( zww(ji,jj,jk) - zww(ji,jj,jk+1) )
-               ! add it to the general tracer trends
-               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta
+               ta(ji,jj,jk) =  ta(ji,jj,jk) + zta                       ! add it to the general tracer trends
                sa(ji,jj,jk) =  sa(ji,jj,jk) + zsa
             END DO
@@ -354,5 +353,6 @@
       !                                                ! ===============
 
-      ! Save the vertical advective trends for diagnostic
+      ! 3. Save the vertical advective trends for diagnostic
+      ! ----------------------------------------------------
       IF( l_trdtra )   THEN
          ! Recompute the vertical advection zta & zsa trends computed 
Index: trunk/NEMO/OPA_SRC/ice_oce.F90
===================================================================
--- trunk/NEMO/OPA_SRC/ice_oce.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/ice_oce.F90	(revision 717)
@@ -39,12 +39,4 @@
       fcalving            !: Iceberg calving 
 # endif
-
-   REAL(wp), PUBLIC, DIMENSION(jpi,jpj) ::   &  !: field exchanges with ice model to ocean
-      sst_io, sss_io , &  !: sea surface temperature (C) and salinity (PSU)
-      u_io  , v_io   , &  !: velocity at ice surface (m/s)
-      fsolar, fnsolar, &  !: solar and non-solar heat fluxes (W/m2)
-      fsalt , fmass  , &  !: salt and freshwater fluxes
-      ftaux , ftauy  , &  !: wind stresses
-      gtaux , gtauy       !: wind stresses
    
    REAL(wp), PUBLIC ::   &  !:
@@ -59,7 +51,4 @@
 #endif
 
-   INTEGER, PUBLIC ::   &  !: namdom : space/time domain (namlist)
-      nfice =  5           !: coupling frequency OPA ICELLN  nfice 
-
    !!----------------------------------------------------------------------
 END MODULE ice_oce
Index: trunk/NEMO/OPA_SRC/lbclnk.F90
===================================================================
--- trunk/NEMO/OPA_SRC/lbclnk.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/lbclnk.F90	(revision 717)
@@ -329,5 +329,5 @@
 
 
-   SUBROUTINE lbc_lnk_3d( pt3d, cd_type, psgn, cd_mpp )
+   SUBROUTINE lbc_lnk_3d( pt3d, cd_type, psgn, cd_mpp, pval )
       !!---------------------------------------------------------------------
       !!                  ***  ROUTINE lbc_lnk_3d  ***
@@ -355,8 +355,10 @@
       CHARACTER(len=3), INTENT( in ), OPTIONAL ::    &
          cd_mpp        ! fill the overlap area only (here do nothing)
+      REAL(wp)        , INTENT(in   ), OPTIONAL           ::   pval      ! background value (used at closed boundaries)
 
       !! * Local declarations
       INTEGER  ::   ji, jk
       INTEGER  ::   ijt, iju
+      REAL(wp) ::   zland
       !!----------------------------------------------------------------------
       !!  OPA 9.0 , LOCEAN-IPSL (2005) 
@@ -365,5 +367,12 @@
       !!----------------------------------------------------------------------
 
-      IF (PRESENT(cd_mpp)) THEN
+      IF( PRESENT( pval ) ) THEN      ! set land value (zero by default)
+         zland = pval
+      ELSE
+         zland = 0.e0
+      ENDIF
+
+
+      IF( PRESENT( cd_mpp ) ) THEN
          ! only fill the overlap area and extra allows 
          ! this is in mpp case. In this module, just do nothing
@@ -385,8 +394,8 @@
             SELECT CASE ( cd_type )
             CASE ( 'T' , 'U' , 'V' , 'W' )             ! T-, U-, V-, W-points
-               pt3d( 1 ,:,jk) = 0.e0
-               pt3d(jpi,:,jk) = 0.e0
-            CASE ( 'F' )                               ! F-point
-               pt3d(jpi,:,jk) = 0.e0
+               pt3d( 1 ,:,jk) = zland
+               pt3d(jpi,:,jk) = zland
+            CASE ( 'F' )                               ! F-point
+               pt3d(jpi,:,jk) = zland
             END SELECT
 
@@ -402,14 +411,14 @@
             CASE ( 'T' , 'U' , 'W' )                   ! T-, U-, W-points
                pt3d(:, 1 ,jk) = pt3d(:,3,jk)
-               pt3d(:,jpj,jk) = 0.e0
+               pt3d(:,jpj,jk) = zland
             CASE ( 'V' , 'F' )                         ! V-, F-points
                pt3d(:, 1 ,jk) = psgn * pt3d(:,2,jk)
-               pt3d(:,jpj,jk) = 0.e0
+               pt3d(:,jpj,jk) = zland
             END SELECT
 
          CASE ( 3 , 4 )                        ! *  North fold  T-point pivot
 
-            pt3d( 1 ,jpj,jk) = 0.e0
-            pt3d(jpi,jpj,jk) = 0.e0
+            pt3d( 1 ,jpj,jk) = zland
+            pt3d(jpi,jpj,jk) = zland
 
             SELECT CASE ( cd_type )
@@ -417,5 +426,5 @@
                DO ji = 2, jpi
                   ijt = jpi-ji+2
-                  pt3d(ji, 1 ,jk) = 0.e0
+                  pt3d(ji, 1 ,jk) = zland
                   pt3d(ji,jpj,jk) = psgn * pt3d(ijt,jpj-2,jk)
                END DO
@@ -427,5 +436,5 @@
                DO ji = 1, jpi-1
                   iju = jpi-ji+1
-                  pt3d(ji, 1 ,jk) = 0.e0
+                  pt3d(ji, 1 ,jk) = zland
                   pt3d(ji,jpj,jk) = psgn * pt3d(iju,jpj-2,jk)
                END DO
@@ -437,5 +446,5 @@
                   DO ji = 2, jpi
                      ijt = jpi-ji+2
-                     pt3d(ji,  1  ,jk) = 0.e0
+                     pt3d(ji,  1  ,jk) = zland
                      pt3d(ji,jpj-1,jk) = psgn * pt3d(ijt,jpj-2,jk)
                      pt3d(ji,jpj  ,jk) = psgn * pt3d(ijt,jpj-3,jk)
@@ -451,6 +460,6 @@
          CASE ( 5 , 6 )                        ! *  North fold  F-point pivot
 
-            pt3d( 1 ,jpj,jk) = 0.e0
-            pt3d(jpi,jpj,jk) = 0.e0
+            pt3d( 1 ,jpj,jk) = zland
+            pt3d(jpi,jpj,jk) = zland
 
             SELECT CASE ( cd_type )
@@ -458,5 +467,5 @@
                DO ji = 1, jpi
                   ijt = jpi-ji+1
-                  pt3d(ji, 1 ,jk) = 0.e0
+                  pt3d(ji, 1 ,jk) = zland
                   pt3d(ji,jpj,jk) = psgn * pt3d(ijt,jpj-1,jk)
                END DO
@@ -464,5 +473,5 @@
                   DO ji = 1, jpi-1
                      iju = jpi-ji
-                     pt3d(ji, 1 ,jk) = 0.e0
+                     pt3d(ji, 1 ,jk) = zland
                      pt3d(ji,jpj,jk) = psgn * pt3d(iju,jpj-1,jk)
                   END DO
@@ -470,5 +479,5 @@
                   DO ji = 1, jpi
                      ijt = jpi-ji+1
-                     pt3d(ji, 1 ,jk) = 0.e0
+                     pt3d(ji, 1 ,jk) = zland
                      pt3d(ji,jpj,jk) = psgn * pt3d(ijt,jpj-2,jk)
                   END DO
@@ -492,8 +501,8 @@
             SELECT CASE ( cd_type )
             CASE ( 'T' , 'U' , 'V' , 'W' )             ! T-, U-, V-, W-points
-               pt3d(:, 1 ,jk) = 0.e0
-               pt3d(:,jpj,jk) = 0.e0
-            CASE ( 'F' )                               ! F-point
-               pt3d(:,jpj,jk) = 0.e0
+               pt3d(:, 1 ,jk) = zland
+               pt3d(:,jpj,jk) = zland
+            CASE ( 'F' )                               ! F-point
+               pt3d(:,jpj,jk) = zland
             END SELECT
 
@@ -506,5 +515,5 @@
 
 
-   SUBROUTINE lbc_lnk_2d( pt2d, cd_type, psgn, cd_mpp )
+   SUBROUTINE lbc_lnk_2d( pt2d, cd_type, psgn, cd_mpp, pval )
       !!---------------------------------------------------------------------
       !!                 ***  ROUTINE lbc_lnk_2d  ***
@@ -532,11 +541,17 @@
       CHARACTER(len=3), INTENT( in ), OPTIONAL ::    &
          cd_mpp        ! fill the overlap area only (here do nothing)
+      REAL(wp)        , INTENT(in   ), OPTIONAL           ::   pval      ! background value (used at closed boundaries)
 
       !! * Local declarations
       INTEGER  ::   ji
       INTEGER  ::   ijt, iju
+      REAL(wp) ::   zland
       !!----------------------------------------------------------------------
-      !!  OPA 8.5, LODYC-IPSL (2002)
-      !!----------------------------------------------------------------------
+
+      IF( PRESENT( pval ) ) THEN      ! set land value (zero by default)
+         zland = pval
+      ELSE
+         zland = 0.e0
+      ENDIF
 
       IF (PRESENT(cd_mpp)) THEN
@@ -556,11 +571,11 @@
          SELECT CASE ( cd_type )
          CASE ( 'T' , 'U' , 'V' , 'W' )                ! T-, U-, V-, W-points
-            pt2d( 1 ,:) = 0.e0
-            pt2d(jpi,:) = 0.e0
+            pt2d( 1 ,:) = zland
+            pt2d(jpi,:) = zland
          CASE ( 'F' )                                  ! F-point, ice U-V point
-            pt2d(jpi,:) = 0.e0 
+            pt2d(jpi,:) = zland
          CASE ( 'I' )                                  ! F-point, ice U-V point
-            pt2d( 1 ,:) = 0.e0 
-            pt2d(jpi,:) = 0.e0 
+            pt2d( 1 ,:) = zland
+            pt2d(jpi,:) = zland
          END SELECT
 
@@ -576,15 +591,15 @@
          CASE ( 'T' , 'U' , 'W' )                      ! T-, U-, W-points
             pt2d(:, 1 ) = pt2d(:,3)
-            pt2d(:,jpj) = 0.e0
+            pt2d(:,jpj) = zland
          CASE ( 'V' , 'F' , 'I' )                      ! V-, F-points, ice U-V point
             pt2d(:, 1 ) = psgn * pt2d(:,2)
-            pt2d(:,jpj) = 0.e0
+            pt2d(:,jpj) = zland
          END SELECT
 
       CASE ( 3 , 4 )                           ! * North fold  T-point pivot
 
-         pt2d( 1 , 1 ) = 0.e0        !!!!!  bug gm ??? !Edmee
-         pt2d( 1 ,jpj) = 0.e0
-         pt2d(jpi,jpj) = 0.e0
+         pt2d( 1 , 1 ) = zland       !!!!!  bug gm ??? !Edmee
+         pt2d( 1 ,jpj) = zland
+         pt2d(jpi,jpj) = zland
 
          SELECT CASE ( cd_type )
@@ -593,5 +608,5 @@
             DO ji = 2, jpi
                ijt = jpi-ji+2
-               pt2d(ji, 1 ) = 0.e0
+               pt2d(ji, 1 ) = zland
                pt2d(ji,jpj) = psgn * pt2d(ijt,jpj-2)
             END DO
@@ -604,5 +619,5 @@
             DO ji = 1, jpi-1
                iju = jpi-ji+1
-               pt2d(ji, 1 ) = 0.e0
+               pt2d(ji, 1 ) = zland
                pt2d(ji,jpj) = psgn * pt2d(iju,jpj-2)
             END DO
@@ -615,5 +630,5 @@
             DO ji = 2, jpi
                ijt = jpi-ji+2
-               pt2d(ji, 1   ) = 0.e0
+               pt2d(ji, 1   ) = zland
                pt2d(ji,jpj-1) = psgn * pt2d(ijt,jpj-2)
                pt2d(ji,jpj  ) = psgn * pt2d(ijt,jpj-3)
@@ -628,5 +643,5 @@
 
          CASE ( 'I' )                                  ! ice U-V point
-            pt2d(:, 1 ) = 0.e0
+            pt2d(:, 1 ) = zland
             pt2d(2,jpj) = psgn * pt2d(3,jpj-1)
             DO ji = 3, jpi
@@ -639,7 +654,7 @@
       CASE ( 5 , 6 )                           ! * North fold  F-point pivot
 
-         pt2d( 1 , 1 ) = 0.e0           !!bug  ???
-         pt2d( 1 ,jpj) = 0.e0
-         pt2d(jpi,jpj) = 0.e0
+         pt2d( 1 , 1 ) = zland          !!bug  ???
+         pt2d( 1 ,jpj) = zland
+         pt2d(jpi,jpj) = zland
 
          SELECT CASE ( cd_type )
@@ -648,5 +663,5 @@
             DO ji = 1, jpi
                ijt = jpi-ji+1
-               pt2d(ji, 1 ) = 0.e0
+               pt2d(ji, 1 ) = zland
                pt2d(ji,jpj) = psgn * pt2d(ijt,jpj-1)
             END DO
@@ -655,5 +670,5 @@
             DO ji = 1, jpi-1
                iju = jpi-ji
-               pt2d(ji, 1 ) = 0.e0
+               pt2d(ji, 1 ) = zland
                pt2d(ji,jpj) = psgn * pt2d(iju,jpj-1)
             END DO
@@ -662,5 +677,5 @@
             DO ji = 1, jpi
                ijt = jpi-ji+1
-               pt2d(ji, 1 ) = 0.e0
+               pt2d(ji, 1 ) = zland
                pt2d(ji,jpj) = psgn * pt2d(ijt,jpj-2)
             END DO
@@ -681,6 +696,6 @@
 
          CASE ( 'I' )                                  ! ice U-V point
-            pt2d( : , 1 ) = 0.e0
-            pt2d( 2 ,jpj) = 0.e0
+            pt2d( : , 1 ) = zland
+            pt2d( 2 ,jpj) = zland
             DO ji = 2 , jpim1
                ijt = jpi - ji + 2
@@ -694,11 +709,11 @@
          SELECT CASE ( cd_type )
          CASE ( 'T' , 'U' , 'V' , 'W' )                ! T-, U-, V-, W-points
-            pt2d(:, 1 ) = 0.e0
-            pt2d(:,jpj) = 0.e0
+            pt2d(:, 1 ) = zland
+            pt2d(:,jpj) = zland
          CASE ( 'F' )                                  ! F-point
-            pt2d(:,jpj) = 0.e0
+            pt2d(:,jpj) = zland
          CASE ( 'I' )                                  ! ice U-V point
-            pt2d(:, 1 ) = 0.e0
-            pt2d(:,jpj) = 0.e0
+            pt2d(:, 1 ) = zland
+            pt2d(:,jpj) = zland
          END SELECT
 
Index: trunk/NEMO/OPA_SRC/lib_mpp.F90
===================================================================
--- trunk/NEMO/OPA_SRC/lib_mpp.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/lib_mpp.F90	(revision 717)
@@ -597,5 +597,5 @@
 #endif
 
-   SUBROUTINE mpp_lnk_3d( ptab, cd_type, psgn, cd_mpp )
+   SUBROUTINE mpp_lnk_3d( ptab, cd_type, psgn, cd_mpp, pval )
       !!----------------------------------------------------------------------
       !!                  ***  routine mpp_lnk_3d  ***
@@ -632,4 +632,5 @@
       CHARACTER(len=3), INTENT( in ), OPTIONAL ::    &
          cd_mpp        ! fill the overlap area only 
+      REAL(wp)        , INTENT(in   ), OPTIONAL           ::   pval      ! background value (used at closed boundaries)
 
       !! * Local variables
@@ -638,8 +639,15 @@
       INTEGER ::   ml_req1, ml_req2, ml_err     ! for key_mpi_isend
       INTEGER ::   ml_stat(MPI_STATUS_SIZE)     ! for key_mpi_isend
+      REAL(wp) ::   zland
       !!----------------------------------------------------------------------
 
       ! 1. standard boundary treatment
       ! ------------------------------
+
+      IF( PRESENT( pval ) ) THEN      ! set land value (zero by default)
+         zland = pval
+      ELSE
+         zland = 0.e0
+      ENDIF
 
       IF( PRESENT( cd_mpp ) ) THEN
@@ -662,8 +670,8 @@
             SELECT CASE ( cd_type )
             CASE ( 'T', 'U', 'V', 'W' )
-               ptab(     1       :jpreci,:,:) = 0.e0
-               ptab(nlci-jpreci+1:jpi   ,:,:) = 0.e0
+               ptab(     1       :jpreci,:,:) = zland
+               ptab(nlci-jpreci+1:jpi   ,:,:) = zland
             CASE ( 'F' )
-               ptab(nlci-jpreci+1:jpi   ,:,:) = 0.e0
+               ptab(nlci-jpreci+1:jpi   ,:,:) = zland
             END SELECT 
          ENDIF
@@ -673,8 +681,8 @@
          SELECT CASE ( cd_type )
          CASE ( 'T', 'U', 'V', 'W' )
-            ptab(:,     1       :jprecj,:) = 0.e0
-            ptab(:,nlcj-jprecj+1:jpj   ,:) = 0.e0
+            ptab(:,     1       :jprecj,:) = zland
+            ptab(:,nlcj-jprecj+1:jpj   ,:) = zland
          CASE ( 'F' )
-            ptab(:,nlcj-jprecj+1:jpj   ,:) = 0.e0
+            ptab(:,nlcj-jprecj+1:jpj   ,:) = zland
          END SELECT
      
@@ -1050,5 +1058,5 @@
 
 
-   SUBROUTINE mpp_lnk_2d( pt2d, cd_type, psgn, cd_mpp )
+   SUBROUTINE mpp_lnk_2d( pt2d, cd_type, psgn, cd_mpp, pval )
       !!----------------------------------------------------------------------
       !!                  ***  routine mpp_lnk_2d  ***
@@ -1084,4 +1092,5 @@
       CHARACTER(len=3), INTENT( in ), OPTIONAL ::    &
          cd_mpp        ! fill the overlap area only 
+      REAL(wp)        , INTENT(in   ), OPTIONAL           ::   pval      ! background value (used at closed boundaries)
 
       !! * Local variables
@@ -1092,5 +1101,12 @@
       INTEGER  ::   ml_req1, ml_req2, ml_err     ! for key_mpi_isend
       INTEGER  ::   ml_stat(MPI_STATUS_SIZE)     ! for key_mpi_isend
+      REAL(wp) ::   zland
       !!----------------------------------------------------------------------
+
+      IF( PRESENT( pval ) ) THEN      ! set land value (zero by default)
+         zland = pval
+      ELSE
+         zland = 0.e0
+      ENDIF
 
       ! 1. standard boundary treatment
@@ -1115,8 +1131,8 @@
             SELECT CASE ( cd_type )
             CASE ( 'T', 'U', 'V', 'W' , 'I' )
-               pt2d(     1       :jpreci,:) = 0.e0
-               pt2d(nlci-jpreci+1:jpi   ,:) = 0.e0
+               pt2d(     1       :jpreci,:) = zland
+               pt2d(nlci-jpreci+1:jpi   ,:) = zland
             CASE ( 'F' )
-               pt2d(nlci-jpreci+1:jpi   ,:) = 0.e0
+               pt2d(nlci-jpreci+1:jpi   ,:) = zland
             END SELECT
          ENDIF
@@ -1126,8 +1142,8 @@
          SELECT CASE ( cd_type )
          CASE ( 'T', 'U', 'V', 'W' , 'I' )
-            pt2d(:,     1       :jprecj) = 0.e0
-            pt2d(:,nlcj-jprecj+1:jpj   ) = 0.e0
+            pt2d(:,     1       :jprecj) = zland
+            pt2d(:,nlcj-jprecj+1:jpj   ) = zland
          CASE ( 'F' )
-            pt2d(:,nlcj-jprecj+1:jpj   ) = 0.e0
+            pt2d(:,nlcj-jprecj+1:jpj   ) = zland
          END SELECT
 
@@ -1394,5 +1410,5 @@
   
             CASE ( 'I' )                                  ! ice U-V point
-               pt2d( 2 ,nlcj) = 0.e0
+               pt2d( 2 ,nlcj) = zland
                DO ji = 2 , nlci-1
                   ijt = iloc - ji + 2
Index: trunk/NEMO/OPA_SRC/opa.F90
===================================================================
--- trunk/NEMO/OPA_SRC/opa.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/opa.F90	(revision 717)
@@ -309,8 +309,4 @@
 #endif
 
-      CALL flx_init                         ! Thermohaline forcing initialization
-
-      CALL flx_fwb_init                     ! FreshWater Budget correction
-
       CALL dia_ptr_init                     ! Poleward TRansports initialization
 
Index: trunk/NEMO/OPA_SRC/restart.F90
===================================================================
--- trunk/NEMO/OPA_SRC/restart.F90	(revision 716)
+++ trunk/NEMO/OPA_SRC/restart.F90	(revision 717)
@@ -19,5 +19,4 @@
    USE phycst          ! physical constants
    USE daymod          ! calendar
-   USE ice_oce         ! ice variables
    USE cpl_oce, ONLY : lk_cpl              !
    USE in_out_manager  ! I/O manager
@@ -138,19 +137,4 @@
       CALL iom_rstput( kt, nitrst, numrow, 'rotn'   , rotn    )
       CALL iom_rstput( kt, nitrst, numrow, 'hdivn'  , hdivn   )
-
-#if defined key_ice_lim        
-      CALL iom_rstput( kt, nitrst, numrow, 'nfice'  , REAL( nfice, wp) )   !  ice computation frequency
-      CALL iom_rstput( kt, nitrst, numrow, 'sst_io' , sst_io  )
-      CALL iom_rstput( kt, nitrst, numrow, 'sss_io' , sss_io  )
-      CALL iom_rstput( kt, nitrst, numrow, 'u_io'   , u_io    )
-      CALL iom_rstput( kt, nitrst, numrow, 'v_io'   , v_io    )
-# if defined key_coupled
-      CALL iom_rstput( kt, nitrst, numrow, 'alb_ice', alb_ice )
-# endif
-#endif
-#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_core 
-      CALL iom_rstput( kt, nitrst, numrow, 'nfbulk' , REAL( nfbulk, wp) )   !  bulk computation frequency
-      CALL iom_rstput( kt, nitrst, numrow, 'gsst'   , gsst    )
-#endif
 
       IF( nn_dynhpg_rst == 1 .OR. lk_vvl ) THEN
@@ -204,5 +188,5 @@
       !!                    has been stored in the restart file.
       !!----------------------------------------------------------------------
-      REAL(wp) ::   zcoef, zkt, zrdt, zrdttra1, zndastp, znfice, znfbulk
+      REAL(wp) ::   zcoef, zkt, zrdt, zrdttra1, zndastp
 #if defined key_ice_lim
       INTEGER  ::   ji, jj
@@ -299,62 +283,4 @@
       ENDIF
 
-      !!sm: TO BE MOVED IN NEW SURFACE MODULE...
-
-#if defined key_ice_lim
-      ! Louvain La Neuve Sea Ice Model
-      IF( iom_varid( numror, 'nfice' ) > 0 ) then 
-         CALL iom_get( numror             , 'nfice'  , znfice  )   ! ice computation frequency
-         CALL iom_get( numror, jpdom_autoglo, 'sst_io' , sst_io  )
-         CALL iom_get( numror, jpdom_autoglo, 'sss_io' , sss_io  )
-         CALL iom_get( numror, jpdom_autoglo, 'u_io'   , u_io    )
-         CALL iom_get( numror, jpdom_autoglo, 'v_io'   , v_io    )
-# if defined key_coupled
-         CALL iom_get( numror, jpdom_autoglo, 'alb_ice', alb_ice )
-# endif
-         IF( znfice /= REAL( nfice, wp ) ) THEN      ! if nfice changed between 2 runs
-            zcoef = REAL( nfice-1, wp ) / znfice
-            sst_io(:,:) = zcoef * sst_io(:,:)
-            sss_io(:,:) = zcoef * sss_io(:,:)
-            u_io  (:,:) = zcoef * u_io  (:,:)
-            v_io  (:,:) = zcoef * v_io  (:,:)
-         ENDIF
-      ELSE
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'rst_read :  LLN sea Ice Model => Ice initialization'
-         IF(lwp) WRITE(numout,*)
-         zcoef = REAL( nfice-1, wp )
-         sst_io(:,:) = zcoef *( tn(:,:,1) + rt0 )          !!bug a explanation is needed here!
-         sss_io(:,:) = zcoef *  sn(:,:,1)
-         zcoef = 0.5 * REAL( nfice-1, wp )
-         DO jj = 2, jpj
-            DO ji = fs_2, jpi   ! vector opt.
-               u_io(ji,jj) = zcoef * ( un(ji-1,jj  ,1) + un(ji-1,jj-1,1) )
-               v_io(ji,jj) = zcoef * ( vn(ji  ,jj-1,1) + vn(ji-1,jj-1,1) )
-            END DO
-         END DO
-# if defined key_coupled
-         alb_ice(:,:) = 0.8 * tmask(:,:,1)
-# endif
-      ENDIF
-#endif
-#if defined key_flx_bulk_monthly || defined key_flx_bulk_daily || defined key_flx_core 
-      ! Louvain La Neuve Sea Ice Model
-      IF( iom_varid( numror, 'nfbulk' ) > 0 ) THEN 
-         CALL iom_get( numror             , 'nfbulk', znfbulk )   ! bulk computation frequency
-         CALL iom_get( numror, jpdom_autoglo, 'gsst'  , gsst    )
-         IF( znfbulk /= REAL(nfbulk, wp) ) THEN      ! if you change nfbulk between 2 runs
-            zcoef = REAL( nfbulk-1, wp ) / znfbulk
-            gsst(:,:) = zcoef * gsst(:,:)
-         ENDIF
-      ELSE
-         IF(lwp) WRITE(numout,*)
-         IF(lwp) WRITE(numout,*) 'rst_read :  LLN sea Ice Model => Ice initialization'
-         IF(lwp) WRITE(numout,*)
-         gsst(:,:) = REAL( nfbulk - 1, wp )*( tn(:,:,1) + rt0 )
-      ENDIF
-#endif
-
-      !!sm: end of TO BE MOVED IN NEW SURFACE MODULE...
-
       IF( iom_varid( numror, 'rhd' ) > 0 ) THEN
          CALL iom_get( numror, jpdom_autoglo, 'rhd' , rhd  )
