Index: trunk/NEMO/LIM_SRC_2/dom_ice_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/dom_ice_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/dom_ice_2.F90	(revision 888)
@@ -1,4 +1,3 @@
 MODULE dom_ice_2
-#if defined key_lim2
    !!======================================================================
    !!                   ***  MODULE  dom_ice  ***
@@ -7,8 +6,8 @@
    !! History :   2.0  !  03-08  (C. Ethe)  Free form and module
    !!----------------------------------------------------------------------
-
+#if defined key_lim2
    !!----------------------------------------------------------------------
    !!   LIM 2.0, UCL-LOCEAN-IPSL (2005)
-   !! $Header$
+   !! $ Id: $
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
Index: trunk/NEMO/LIM_SRC_2/ice_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/ice_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/ice_2.F90	(revision 888)
@@ -4,14 +4,13 @@
    !! 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_lim2
    !!----------------------------------------------------------------------
    !!   'key_lim2' :                                  LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
-   !! History :
-   !!   2.0  !  03-08  (C. Ethe)  F90: Free form and module
-   !!----------------------------------------------------------------------
-   !!  LIM 2.0, UCL-LOCEAN-IPSL (2005)
-   !! $Header$
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
+   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006)
+   !! $ Id: $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
    !! * Modules used
@@ -21,103 +20,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
Index: trunk/NEMO/LIM_SRC_2/iceini_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/iceini_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/iceini_2.F90	(revision 888)
@@ -17,7 +17,7 @@
    USE dom_oce
    USE dom_ice_2
-   USE in_out_manager
    USE ice_oce         ! ice variables
-   USE flx_oce
+   USE sbc_oce         ! surface boundary condition: ocean
+   USE sbc_ice         ! surface boundary condition: ice
    USE phycst          ! Define parameters for the routines
    USE ocfzpt
@@ -27,4 +27,5 @@
    USE limrst_2   
    USE ini1d           ! initialization of the 1D configuration
+   USE in_out_manager
       
    IMPLICIT NONE
@@ -40,5 +41,5 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header$ 
+   !! $ Id: $ 
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    !!----------------------------------------------------------------------
@@ -62,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_2                  ! ice mesh initialization
Index: trunk/NEMO/LIM_SRC_2/icestp_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/icestp_2.F90	(revision 833)
+++ 	(revision )
@@ -1,286 +1,0 @@
-MODULE icestp_2
-   !!======================================================================
-   !!                       ***  MODULE icestp_2   ***
-   !!   Sea-Ice model : LIM Sea ice model time-stepping
-   !!======================================================================
-   !! History :   1.0  !  99-11  (M. Imbard)  Original code
-   !!        !  01-03  (D. Ludicone, E. Durand, G. Madec) free surf.
-   !!   2.0  !  02-09  (G. Madec, C. Ethe)  F90: Free form and module
-   !!----------------------------------------------------------------------
-#if defined key_lim2
-   !!----------------------------------------------------------------------
-   !!   'key_lim2' :                                  LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
-   !!----------------------------------------------------------------------
-   !!   ice_stp_2       : sea-ice model time-stepping
-   !!----------------------------------------------------------------------
-   USE dom_oce
-   USE oce  ! dynamics and tracers variables
-   USE in_out_manager
-   USE ice_oce         ! ice variables
-   USE flx_oce         ! forcings variables
-   USE dom_ice_2
-   USE cpl_oce
-   USE daymod
-   USE phycst          ! Define parameters for the routines
-   USE taumod
-   USE ice_2
-   USE ocesbc
-   USE lbclnk
-   USE limdyn_2
-   USE limtrp_2
-   USE limthd_2
-   USE limflx_2
-   USE limdia_2
-   USE limwri_2
-   USE limrst_2
-   USE limdmp_2        ! Ice damping
-   USE prtctl          ! Print control
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   ice_stp_2  ! called by step.F90
-
-   !! * Substitutions
-#  include "domzgr_substitute.h90"
-#  include "vectopt_loop_substitute.h90"
-   !!-----------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header$ 
-   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
-   !!-----------------------------------------------------
-
-CONTAINS
-
-   SUBROUTINE ice_stp_2 ( kt )
-      !!---------------------------------------------------------------------
-      !!                  ***  ROUTINE ice_stp_2  ***
-      !!                   
-      !! ** Purpose :   Louvain la Neuve Sea Ice Model time stepping 
-      !!
-      !! ** Action  : - call the ice dynamics routine 
-      !!              - call the ice advection/diffusion routine 
-      !!              - call the ice thermodynamics routine 
-      !!              - call the routine that computes mass and 
-      !!                heat fluxes at the ice/ocean interface
-      !!              - save the outputs 
-      !!              - save the outputs for restart when necessary
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT(in)          ::   kt       ! ocean time-step index
-
-      INTEGER                      ::   ji, jj   ! dummy loop indices
-      REAL(wp), DIMENSION(jpi,jpj) ::   zsss_io, zsss2_io, zsss3_io          ! tempory workspaces
-      !!----------------------------------------------------------------------
-
-      IF( kt == nit000 ) THEN
-         IF( lk_cpl ) THEN
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)'
-            IF(lwp) WRITE(numout,*) '~~~~~~~   coupled case'
-         ELSE
-            IF(lwp) WRITE(numout,*)
-            IF(lwp) WRITE(numout,*) 'ice_stp : Louvain la Neuve Ice Model (LIM)' 
-            IF(lwp) WRITE(numout,*) '~~~~~~~   forced case using bulk formulea'
-         ENDIF
-         !  Initialize fluxes fields
-         gtaux(:,:) = 0.e0
-         gtauy(:,:) = 0.e0
-      ENDIF
-
-      ! Temperature , salinity and horizonta wind
-      ! sst_io  and sss_io, u_io and v_io  are initialized at nit000 in limistate.F90 (or limrst.F90) with :
-      !              sst_io = sst_io + (nfice - 1) * (tn(:,:,1)+rt0 )
-      !              sss_io = sss_io + (nfice - 1) * sn(:,:,1)
-      !              u_io  = u_io  + (nfice - 1) * 0.5 * ( un(ji-1,jj  ,1) + un(ji-1,jj-1,1) )
-      !              v_io  = v_io  + (nfice - 1) * 0.5 * ( vn(ji  ,jj-1,1) + vn(ji-1,jj-1,1) )
-      !    cumulate fields
-      !
-      sst_io(:,:) = sst_io(:,:) + tn(:,:,1) + rt0
-      sss_io(:,:) = sss_io(:,:) + sn(:,:,1)
-
- 
-      ! vectors at F-point
-      DO jj = 2, jpj
-         DO ji = fs_2, jpi   ! vector opt.
-            u_io(ji,jj) = u_io(ji,jj) + 0.5 * ( un(ji-1,jj  ,1) + un(ji-1,jj-1,1) )
-            v_io(ji,jj) = v_io(ji,jj) + 0.5 * ( vn(ji  ,jj-1,1) + vn(ji-1,jj-1,1) )
-         END DO
-      END DO
-      
-      
-      IF( MOD( kt-1, nfice ) == 0 ) THEN
-         
-         ! The LIM model is going to be call
-         sst_io(:,:) = sst_io(:,:) / FLOAT( nfice ) * tmask(:,:,1)
-         sss_io(:,:) = sss_io(:,:) / FLOAT( nfice )
-
-         ! stress from ocean U- and V-points to ice U,V point
-         DO jj = 2, jpj
-            DO ji = fs_2, jpi   ! vector opt.
-               gtaux(ji,jj) = 0.5 * ( taux(ji-1,jj  ) + taux(ji-1,jj-1) )
-               gtauy(ji,jj) = 0.5 * ( tauy(ji  ,jj-1) + tauy(ji-1,jj-1) )
-               u_io  (ji,jj) = u_io(ji,jj) / FLOAT( nfice )
-               v_io  (ji,jj) = v_io(ji,jj) / FLOAT( nfice )
-            END DO
-         END DO
-
-         ! lateral boundary condition
-         CALL lbc_lnk( gtaux(:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
-         CALL lbc_lnk( gtauy(:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
-         CALL lbc_lnk( u_io (:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
-         CALL lbc_lnk( v_io (:,:), 'I', -1. )   ! I-point (i.e. ice U-V point)
-
-!!gmbug  in the ocean freezing point computed as :
-!!gm           fzptn (ji,jj) = ( -0.0575 + 1.710523e-3 * SQRT( sn(ji,jj,1) )   &
-!!gm                                     - 2.154996e-4 *       sn(ji,jj,1)   ) * sn(ji,jj,1)   !!   &
-!!gm           !!                        - 7.53e-4 * pressure
-!!gm
-!!gm!bug this is much more accurate and efficient computation
-!!gm       **************************************************
-!!gm freezing point from unesco:
-!!gm     real function tf(s,p)
-!!gm   function to compute the freezing point of seawater
-!!gm
-!!gm   reference: unesco tech. papers in the marine science no. 28. 1978
-!!gm   eighth report jpots
-!!gm   annex 6 freezing point of seawater f.j. millero pp.29-35.
-!!gm
-!!gm  units:
-!!gm         pressure      p          decibars
-!!gm         salinity      s          pss-78
-!!gm         temperature   tf         degrees celsius
-!!gm         freezing pt.
-!!gm************************************************************
-!!gm  checkvalue: tf= -2.588567 deg. c for s=40.0, p=500. decibars
-!!gm     tf=(-.0575+1.710523e-3*sqrt(abs(s))-2.154996e-4*s)*s-7.53e-4*p
-!!gm     return
-!!gm     end
-!!gm!bug
-
-
-!!gm      DO jj = 1, jpj
-!!gm         DO ji = 1, jpi
-!!gm            tfu(ji,jj)  = (  rt0 + ( - 0.0575                              &
-!!gm               &                     + 1.710523e-3 * SQRT( sss_io(ji,jj) )   &
-!!gm               &                     - 2.154996e-4 *       sss_io(ji,jj)   ) * sss_io(ji,jj)  ) * tms(ji,jj)
-!!gm         END DO
-!!gm      END DO
-!!gm
-      zsss_io (:,:) = SQRT( sss_io(:,:) ) 
-      zsss2_io(:,:) =  sss_io(:,:) *  sss_io(:,:)
-      zsss3_io(:,:) = zsss_io(:,:) * zsss_io(:,:) * zsss_io(:,:)
-
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            tfu(ji,jj)  = ABS ( rt0 - 0.0575       *   sss_io(ji,jj)   &
-               &                    + 1.710523e-03 * zsss3_io(ji,jj)   &
-               &                    - 2.154996e-04 * zsss2_io(ji,jj) ) * tms(ji,jj)
-         END DO
-      END DO
-   
-         
-
-         IF(ln_ctl) THEN         ! print mean trends (used for debugging)
-            CALL prt_ctl_info('Ice Forcings ')
-            CALL prt_ctl(tab2d_1=qsr_oce ,clinfo1=' qsr_oce  : ', tab2d_2=qsr_ice , clinfo2=' qsr_ice   : ')
-            CALL prt_ctl(tab2d_1=qnsr_oce,clinfo1=' qnsr_oce : ', tab2d_2=qnsr_ice, clinfo2=' qnsr_ice  : ')
-            CALL prt_ctl(tab2d_1=evap    ,clinfo1=' evap     : ')
-            CALL prt_ctl(tab2d_1=tprecip ,clinfo1=' precip   : ', tab2d_2=sprecip , clinfo2=' Snow      : ')
-            CALL prt_ctl(tab2d_1=gtaux   ,clinfo1=' u-stress : ', tab2d_2=gtauy   , clinfo2=' v-stress  : ')
-            CALL prt_ctl(tab2d_1=sst_io  ,clinfo1=' sst      : ', tab2d_2=sss_io  , clinfo2=' sss       : ')
-            CALL prt_ctl(tab2d_1=u_io    ,clinfo1=' u_io     : ', tab2d_2=v_io    , clinfo2=' v_io      : ')
-            CALL prt_ctl(tab2d_1=hsnif   ,clinfo1=' hsnif  1 : ', tab2d_2=hicif   , clinfo2=' hicif     : ')
-            CALL prt_ctl(tab2d_1=frld    ,clinfo1=' frld   1 : ', tab2d_2=sist    , clinfo2=' sist      : ')
-         ENDIF
-
-         !                                                           !-----------------------!
-         CALL lim_rst_opn_2( kt )                                    ! Open Ice restart file !
-         !                                                           !-----------------------!
-
-         !                                                           !--------------!
-         CALL lim_dyn_2( kt )                                        ! Ice dynamics !   ( rheology/dynamics )
-         !                                                           !--------------!
-         IF(ln_ctl) THEN
-            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  2 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
-            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   2 : ', tab2d_2=sist  , clinfo2=' sist      : ')
-         ENDIF
-
-
-         !                                                           !---------------!
-         CALL lim_trp_2( kt )                                        ! Ice transport !  ( Advection/diffusion )
-         !                                                           !---------------!
-         IF(ln_ctl) THEN
-            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  3 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
-            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   3 : ', tab2d_2=sist  , clinfo2=' sist      : ')
-         ENDIF
-
-         !                                                           !-------------!
-         IF( ln_limdmp ) THEN                                        ! Ice damping !
-         !                                                           !-------------!
-#if defined key_agrif
-            IF( Agrif_Root() )   CALL lim_dmp_2(kt)
-#else
-            CALL lim_dmp_2(kt)                                         
-#endif
-         ENDIF
-         !                                                           !--------------------!
-         CALL lim_thd_2( kt )                                        ! Ice thermodynamics !
-         !                                                           !--------------------!
-         IF(ln_ctl) THEN
-            CALL prt_ctl(tab2d_1=hsnif ,clinfo1=' hsnif  4 : ', tab2d_2=hicif , clinfo2=' hicif     : ')
-            CALL prt_ctl(tab2d_1=frld  ,clinfo1=' frld   4 : ', tab2d_2=sist  , clinfo2=' sist      : ')
-         ENDIF
-
-
-         ! Mass and heat fluxes from ice to ocean
-         !                                                           !------------------------------!
-         CALL lim_flx_2                                              ! Ice/Ocean Mass & Heat fluxes !
-         !                                                           !------------------------------!
-
-         IF( MOD( kt + nfice -1, ninfo ) == 0 .OR. ntmoy == 1 )  THEN        !-----------------!
-            CALL lim_dia_2( kt )                                             ! Ice Diagnostics !
-         ENDIF                                                               !-----------------!
-
-         !                                                           !-------------!
-         CALL lim_wri_2( kt )                                        ! Ice outputs !
-         !                                                           !-------------!
-
-         !                                                           !------------------------!
-         IF( lrst_ice ) CALL lim_rst_write_2( kt )                   ! Write Ice restart file !
-         !                                                           !------------------------!
-
-         ! Re-initialization of forcings
-         qsr_oce (:,:) = 0.e0
-         qsr_ice (:,:) = 0.e0
-         qnsr_oce(:,:) = 0.e0
-         qnsr_ice(:,:) = 0.e0 
-         dqns_ice(:,:) = 0.e0 
-         tprecip (:,:) = 0.e0 
-         sprecip (:,:) = 0.e0
-         fr1_i0  (:,:) = 0.e0
-         fr2_i0  (:,:) = 0.e0
-         evap    (:,:) = 0.e0
-#if defined key_coupled 
-         rrunoff (:,:) = 0.e0
-         calving (:,:) = 0.e0
-#else
-         qla_ice (:,:) = 0.e0
-         dqla_ice(:,:) = 0.e0
-#endif
-      ENDIF
-      !
-   END SUBROUTINE ice_stp_2
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option           Dummy module      NO LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE ice_stp_2 ( kt )     ! Dummy routine
-      WRITE(*,*) 'ice_stp_2: You should not have seen this print! error?', kt
-   END SUBROUTINE ice_stp_2
-#endif
-
-   !!======================================================================
-END MODULE icestp_2
Index: trunk/NEMO/LIM_SRC_2/limadv_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limadv_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limadv_2.F90	(revision 888)
@@ -33,5 +33,5 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header$ 
+   !! $ Id: $ 
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    !!----------------------------------------------------------------------
Index: trunk/NEMO/LIM_SRC_2/limdia_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limdia_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limdia_2.F90	(revision 888)
@@ -19,4 +19,5 @@
    USE par_ice_2       ! ice parameters
    USE ice_oce         ! ice variables
+   USE sbc_oce         ! surface boundary condition variables
    USE daymod          !
    USE dom_ice_2       !
@@ -28,5 +29,5 @@
    PRIVATE
 
-   PUBLIC               lim_dia_2          ! called by ice_step
+   PUBLIC               lim_dia_2          ! called by sbc_ice_lim_2
    INTEGER, PUBLIC ::   ntmoy   = 1 ,   &  !: instantaneous values of ice evolution or averaging ntmoy
       &                 ninfo   = 1        !: frequency of ouputs on file ice_evolu in case of averaging
@@ -58,5 +59,5 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limdia.F90,v 1.9 2007/06/29 17:03:12 opalod Exp $ 
+   !! $ Id: $
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -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_2/limdmp_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limdmp_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limdmp_2.F90	(revision 888)
@@ -40,5 +40,5 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0 , UCL-LOCEAN-IPSL  (2006)
-   !! $Header$
+   !! $ Id: $
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
Index: trunk/NEMO/LIM_SRC_2/limdyn_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limdyn_2.F90	(revision 888)
@@ -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_lim2
    !!----------------------------------------------------------------------
@@ -11,31 +16,32 @@
    !!    lim_dyn_init_2 : initialization and namelist read
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE phycst
-   USE in_out_manager  ! I/O manager
-   USE dom_ice_2
-   USE dom_oce         ! ocean space and time domain
-   USE ice_2
-   USE ice_oce
-   USE iceini_2
-   USE limistate_2
-   USE limrhg_2        ! 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_2          !
+   USE ice_oce        !
+   USE dom_ice_2      !
+   USE iceini_2       !
+   USE limistate_2    !
+   USE limrhg_2       ! 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_2  ! routine called by ice_step
+   PUBLIC   lim_dyn_2 ! routine called by sbc_ice_lim
 
    !! * Module variables
    REAL(wp)  ::  rone    = 1.e0   ! constant value
 
-   !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! $ Id: $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -46,5 +52,5 @@
       !!               ***  ROUTINE lim_dyn_2  ***
       !!               
-      !! ** Purpose :   compute ice velocity and ocean-ice stress
+      !! ** Purpose :   compute ice velocity and ocean-ice friction velocity
       !!                
       !! ** Method  : 
@@ -52,39 +58,27 @@
       !! ** Action  : - Initialisation
       !!              - Call of the dynamic routine for each hemisphere
-      !!              - computation of the stress at the ocean surface         
+      !!              - computation of the friction velocity at the sea-ice base
       !!              - 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 , zmod,             &
-         ztglx , ztgly ,           &
-         zt11, zt12, zt21, zt22 ,  &
-         zustm, zsfrld, zsfrldm4,  &
-         zu_ice, zv_ice, ztair2
-      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
       !!---------------------------------------------------------------------
 
-      IF( kt == nit000  )   CALL lim_dyn_init_2   ! Initialization (first time-step only)
+      IF( kt == nit000 )   CALL lim_dyn_init_2   ! Initialization (first time-step only)
       
-      IF ( ln_limdyn ) THEN
-
+      IF( ln_limdyn ) THEN
+         !
          ! Mean ice and snow thicknesses.          
          hsnm(:,:)  = ( 1.0 - frld(:,:) ) * hsnif(:,:)
          hicm(:,:)  = ( 1.0 - frld(:,:) ) * hicif(:,:)
-
-         u_oce(:,:)  = u_io(:,:) * tmu(:,:)
-         v_oce(:,:)  = v_io(:,:) * tmu(:,:)
-       
-         !                                         ! Rheology (ice dynamics)
-         !                                         ! ========
+         !
+         !                                     ! Rheology (ice dynamics)
+         !                                     ! ========
          
          !  Define the j-limits where ice rheology is computed
@@ -94,17 +88,14 @@
             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_2( i_j1, i_jpj )
-
+            !
          ELSE                                 ! optimization of the computational area
-
+            !
             DO jj = 1, jpj
                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
-
+            END DO
+            !
             IF( l_jeq ) THEN                     ! local domain include both hemisphere
                !                                 ! Rheology is computed in each hemisphere
@@ -118,7 +109,7 @@
                i_j1 = MAX( 1, i_j1-1 )
                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : NH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
-    
+               ! 
                CALL lim_rhg_2( i_j1, i_jpj )
-    
+               ! 
                ! Southern hemisphere
                i_j1  =  1 
@@ -129,7 +120,7 @@
                i_jpj = MIN( jpj, i_jpj+2 )
                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : SH i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
-    
+               ! 
                CALL lim_rhg_2( i_j1, i_jpj )
-    
+               ! 
             ELSE                                 ! local domain extends over one hemisphere only
                !                                 ! Rheology is computed only over the ice cover
@@ -148,105 +139,47 @@
     
                IF(ln_ctl)   WRITE(numout,*) 'lim_dyn : one hemisphere: i_j1 = ', i_j1, ' ij_jpj = ', i_jpj
-    
+               ! 
                CALL lim_rhg_2( i_j1, i_jpj )
-
+               !
             ENDIF
-
+            !
          ENDIF
 
-         IF(ln_ctl)   THEN 
-            CALL prt_ctl(tab2d_1=u_oce , clinfo1=' lim_dyn  : u_oce :', tab2d_2=v_oce , clinfo2=' v_oce :')
-            CALL prt_ctl(tab2d_1=u_ice , clinfo1=' lim_dyn  : u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
-         ENDIF
-         
-         !                                         ! Ice-Ocean stress
-         !                                         ! ================
+         IF(ln_ctl)   CALL prt_ctl(tab2d_1=ui_ice , clinfo1=' lim_dyn  : ui_ice :', tab2d_2=vi_ice , clinfo2=' vi_ice :')
+         
+         ! 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
+         ! frictional velocity at T-point
          DO jj = 2, jpjm1
-            zsang  = SIGN(1.e0, gphif(1,jj-1) ) * sangvg
-            DO ji = 2, jpim1
-               ! computation of wind stress over ocean in X and Y direction
-#if defined key_coupled && defined key_lim_cp1
-               ztairx =  frld(ji-1,jj  ) * gtaux(ji-1,jj  ) + frld(ji,jj  ) * gtaux(ji,jj  )      &
-                  &    + frld(ji-1,jj-1) * gtaux(ji-1,jj-1) + frld(ji,jj-1) * gtaux(ji,jj-1)
-
-               ztairy =  frld(ji-1,jj  ) * gtauy(ji-1,jj  ) + frld(ji,jj  ) * gtauy(ji,jj  )      &
-                  &    + frld(ji-1,jj-1) * gtauy(ji-1,jj-1) + frld(ji,jj-1) * gtauy(ji,jj-1)
-#else
-               zsfrld  = frld(ji,jj) + frld(ji-1,jj) + frld(ji-1,jj-1) + frld(ji,jj-1)
-               ztairx  = zsfrld * gtaux(ji,jj)
-               ztairy  = zsfrld * gtauy(ji,jj)
-#endif
-               zsfrldm4 = 4 - frld(ji,jj) - frld(ji-1,jj) - frld(ji-1,jj-1) - frld(ji,jj-1)
-               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
-               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
-               zmod     = SQRT( zu_ice * zu_ice + zv_ice * zv_ice ) 
-               ztglx   = zsfrldm4 * rhoco * zmod * ( cangvg * zu_ice - zsang * zv_ice ) 
-               ztgly   = zsfrldm4 * rhoco * zmod * ( cangvg * zv_ice + zsang * zu_ice ) 
-
-               tio_u(ji,jj) = - ( ztairx + 1.0 * ztglx ) / ( 4 * rau0 )
-               tio_v(ji,jj) = - ( ztairy + 1.0 * ztgly ) / ( 4 * rau0 )
+            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
-         
-         ! computation of friction velocity
+         !
+      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
+         !
+         zcoef = SQRT( 0.5 ) / rau0
          DO jj = 2, jpjm1
-            DO ji = 2, jpim1
-
-               zu_ice   = u_ice(ji-1,jj-1) - u_oce(ji-1,jj-1)
-               zv_ice   = v_ice(ji-1,jj-1) - v_oce(ji-1,jj-1)
-               zt11  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice )
-
-               zu_ice   = u_ice(ji-1,jj) - u_oce(ji-1,jj)
-               zv_ice   = v_ice(ji-1,jj) - v_oce(ji-1,jj)
-               zt12  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
-
-               zu_ice   = u_ice(ji,jj-1) - u_oce(ji,jj-1)
-               zv_ice   = v_ice(ji,jj-1) - v_oce(ji,jj-1)
-               zt21  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
-
-               zu_ice   = u_ice(ji,jj) - u_oce(ji,jj)
-               zv_ice   = v_ice(ji,jj) - v_oce(ji,jj)
-               zt22  = rhoco * ( zu_ice * zu_ice + zv_ice * zv_ice ) 
-
-               ztair2 = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj)
-
-               zustm =  ( 1 - frld(ji,jj) ) * 0.25 * ( zt11 + zt12 + zt21 + zt22 )        &
-                  &  +        frld(ji,jj)   * SQRT( ztair2 )
-
-               ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
+            DO ji = fs_2, fs_jpim1   ! vector opt.
+               ust2s(ji,jj) = zcoef * tms(ji,jj) * SQRT(  utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj)   &
+                  &                                     + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) )
             END DO
          END DO
-
-       ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean
-                    
-          DO jj = 2, jpjm1
-             DO ji = 2, jpim1
-#if defined key_coupled && defined key_lim_cp1
-                tio_u(ji,jj) = - (  gtaux(ji  ,jj  ) + gtaux(ji-1,jj  )       &
-                   &              + gtaux(ji-1,jj-1) + gtaux(ji  ,jj-1) ) / ( 4 * rau0 )
-
-                tio_v(ji,jj) = - (  gtauy(ji  ,jj )  + gtauy(ji-1,jj  )       &
-                   &              + gtauy(ji-1,jj-1) + gtauy(ji  ,jj-1) ) / ( 4 * rau0 )
-#else
-                tio_u(ji,jj) = - gtaux(ji,jj) / rau0
-                tio_v(ji,jj) = - gtauy(ji,jj) / rau0 
-#endif
-                ztair2       = gtaux(ji,jj) * gtaux(ji,jj) + gtauy(ji,jj) * gtauy(ji,jj)
-                zustm        = SQRT( ztair2  )
-
-                ust2s(ji,jj) = ( zustm / rau0 ) * ( rone + sdvt(ji,jj) ) * tms(ji,jj)
-            END DO
-         END DO
-
+         !
       ENDIF
-
+      !
       CALL lbc_lnk( ust2s, 'T',  1. )   ! T-point
-      CALL lbc_lnk( tio_u, 'I', -1. )   ! I-point (i.e. ice U-V point)
-      CALL lbc_lnk( tio_v, 'I', -1. )   ! I-point (i.e. ice U-V point)
-
-      IF(ln_ctl) THEN 
-            CALL prt_ctl(tab2d_1=tio_u , clinfo1=' lim_dyn  : tio_u :', tab2d_2=tio_v , clinfo2=' tio_v :')
-            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_2
@@ -257,14 +190,11 @@
       !!                  ***  ROUTINE lim_dyn_init_2  ***
       !!
-      !! ** 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,     &
@@ -273,11 +203,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_2: ice parameters for ice dynamics '
@@ -291,5 +218,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
@@ -303,14 +230,12 @@
       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_2
 
Index: trunk/NEMO/LIM_SRC_2/limflx_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limflx_2.F90	(revision 833)
+++ 	(revision )
@@ -1,219 +1,0 @@
-MODULE limflx_2
-   !!======================================================================
-   !!                       ***  MODULE limflx_2   ***
-   !!           computation of the flux at the sea ice/ocean interface
-   !!======================================================================
-#if defined key_lim2
-   !!----------------------------------------------------------------------
-   !!   'key_lim2'                                    LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   lim_flx_2  : flux at the ice / ocean interface
-   !!----------------------------------------------------------------------
-   !! * Modules used
-   USE par_oce          ! ocean parameters
-   USE phycst           ! physical constants
-   USE ocfzpt           ! surface ocean freezing point
-   USE ice_oce          ! sea-ice variable
-   USE flx_oce          ! sea-ice/ocean forcings variables
-   USE ice_2             ! LIM sea-ice variables
-   USE flxblk             ! bulk formulea
-   USE lbclnk           ! ocean lateral boundary condition
-   USE in_out_manager   ! I/O manager
-   USE albedo           ! albedo parameters
-   USE prtctl           ! Print control
-
-   IMPLICIT NONE
-   PRIVATE
-
-   !! * Routine accessibility
-   PUBLIC lim_flx_2     ! called by lim_step_2
-
-   !! * Module variables
-   REAL(wp)  ::           &  ! constant values
-      epsi16 = 1.e-16  ,  &
-      rzero  = 0.e0    ,  &
-      rone   = 1.e0
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE lim_flx_2
-      !!-------------------------------------------------------------------
-      !!                ***  ROUTINE lim_flx_2 ***
-      !!  
-      !! ** Purpose : Computes the mass and heat fluxes to the ocean
-      !!         
-      !! ** Action  : - Initialisation of some variables
-      !!              - comput. of the fluxes at the sea ice/ocean interface
-      !!     
-      !! ** 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
-      !!
-      !! ** References :
-      !!       H. Goosse et al. 1996, Bul. Soc. Roy. Sc. Liege, 65, 87-90
-      !!         original    : 00-01 (LIM)
-      !!         addition    : 02-07 (C. Ethe, G. Madec)
-      !!---------------------------------------------------------------------
-      !! * Local variables
-      INTEGER ::   ji, jj         ! dummy loop indices
-
-      INTEGER ::   &
-         ifvt, i1mfr, idfr ,   &  ! some switches
-         iflt, ial, iadv, ifral, ifrdv
-      
-      REAL(wp) ::   &
-         zinda  ,              &  ! switch for testing the values of ice concentration
-         z1mthcm                  ! 1 - thcm
-!!         zfcm1  ,              &  ! solar  heat fluxes
-!!         zfcm2  ,              &  !  non solar heat fluxes
-#if defined key_lim_fdd   
-      REAL(wp) ::   &
-         zfons,                &  ! salt exchanges at the ice/ocean interface
-         zpme                     ! freshwater exchanges at the ice/ocean interface
-#else
-      REAL(wp) ::   &
-         zprs  , zfons,        &  ! salt exchanges at the ice/ocean interface
-         zpmess                   ! freshwater exchanges at the ice/ocean interface
-#endif
-      REAL(wp), DIMENSION(jpi,jpj) ::  &
-         zfcm1  ,              &  ! solar  heat fluxes
-         zfcm2                    !  non solar heat fluxes      
-#if defined key_coupled    
-      REAL(wp), DIMENSION(jpi,jpj) ::  &
-         zalb  ,               &  ! albedo of ice under overcast sky
-         zalcn ,               &  ! albedo of ocean under overcast sky
-         zalbp ,               &  ! albedo of ice under clear sky
-         zaldum                   ! albedo of ocean under clear sky
-#endif
-      !!---------------------------------------------------------------------
-     
-      !---------------------------------!
-      !      Sea ice/ocean interface    !
-      !---------------------------------!
-       
-       
-      !      heat flux at the ocean surface
-      !-------------------------------------------------------
-       
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-            zinda   = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) )
-            ifvt    = zinda  *  MAX( rzero , SIGN( rone, -phicif(ji,jj) ) )
-            i1mfr   = 1.0 - MAX( rzero , SIGN( rone ,  - ( 1.0 - frld(ji,jj) ) ) )
-            idfr    = 1.0 - MAX( rzero , SIGN( rone , frld(ji,jj) - pfrld(ji,jj) ) )
-            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )
-            ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
-            iadv    = ( 1  - i1mfr ) * zinda
-            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
-            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
-            z1mthcm =   1. - thcm(ji,jj)       
-            !   computation the solar flux at ocean surface
-            zfcm1(ji,jj)   = pfrld(ji,jj) * qsr_oce(ji,jj)  + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj)
-            !  computation the non solar heat flux at ocean surface
-            zfcm2(ji,jj) =  - z1mthcm * zfcm1(ji,jj)   &
-               &           + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            &
-               &           + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) / rdt_ice   &
-               &           + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) / rdt_ice
-
-            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ???
-            
-            fsolar (ji,jj) = zfcm1(ji,jj)                                       ! solar heat flux 
-
-            fnsolar(ji,jj) = zfcm2(ji,jj) - fdtcn(ji,jj)                        ! non solar heat flux
-         END DO
-      END DO
-  
-       
-      !      mass flux at the ocean surface
-      !-------------------------------------------------------
-       
-      DO jj = 1, jpj
-         DO ji = 1, jpi
-#if defined key_lim_fdd
-            !  case of realistic freshwater flux (Tartinville et al., 2001)
-            
-            !  computing freshwater exchanges at the ice/ocean interface
-            zpme = - evap(ji,jj) * frld(ji,jj)            &   !  evaporation over oceanic fraction
-               &   + tprecip(ji,jj)                            &   !  total precipitation
-               &   - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )  &   !  remov. snow precip over ice
-               &   - rdmsnif(ji,jj) / rdt_ice                   !  freshwaterflux due to snow melting 
-            
-            !  computing salt exchanges at the ice/ocean interface
-            zfons =  ( soce - sice ) * ( rdmicif(ji,jj) / rdt_ice ) 
-            
-            !  converting the salt flux from ice to a freshwater flux from ocean
-            fsalt(ji,jj) = zfons / ( sss_io(ji,jj) + epsi16 )
-            
-            !  freshwater masses
-            fmass(ji,jj) = - zpme 
-#else
-            !  case of freshwater flux equivalent as salt flux
-            !  dilution effect due to evaporation and precipitation
-            zprs  = ( tprecip(ji,jj) - sprecip(ji,jj) * ( 1. - pfrld(ji,jj) ) ) * soce  
-            !  freshwater flux
-            zfons = rdmicif(ji,jj) * ( soce - sice )  &  !  fwf : ice formation and melting
-               &   -  dmgwi(ji,jj) * sice             &  !  fwf : salt flx needed to bring the fresh snow to sea/ice salinity
-               &   + rdmsnif(ji,jj) * soce               !  fwf to ocean due to snow melting 
-            !  salt exchanges at the ice/ocean interface
-            zpmess         =  zprs - zfons / rdt_ice - evap(ji,jj) * soce * frld(ji,jj)
-            fsalt(ji,jj) =  - zpmess
-#endif
-         END DO
-      END DO
-
-
-      !-------------------------------------------------------------------!
-      !  computation of others transmitting variables from ice to ocean   !
-      !------------------------------------------ ------------------------!
-
-      !-----------------------------------------------!
-      !   Storing the transmitted variables           !
-      !-----------------------------------------------!
-
-      ftaux (:,:) = - tio_u(:,:) * rau0   ! taux ( ice: N/m2/rau0, ocean: N/m2 )
-      ftauy (:,:) = - tio_v(:,:) * rau0   ! tauy ( ice: N/m2/rau0, ocean: N/m2 )                
-      freeze(:,:) = 1.0 - frld(:,:)       ! Sea ice cover            
-      tn_ice(:,:) = sist(:,:)             ! Ice surface temperature                      
-
-#if defined key_coupled            
-      zalb  (:,:) = 0.e0
-      zalcn (:,:) = 0.e0
-      zalbp (:,:) = 0.e0
-      zaldum(:,:) = 0.e0
-
-      !------------------------------------------------!
-      !  2) Computation of snow/ice and ocean albedo   !
-      !------------------------------------------------!
-      CALL flx_blk_albedo( zalb, zalcn, zalbp, zaldum )
-
-      alb_ice(:,:) =  0.5 * zalbp(:,:) + 0.5 * zalb (:,:)   ! Ice albedo                       
-#endif
-
-      IF(ln_ctl) THEN
-         CALL prt_ctl(tab2d_1=fsolar, clinfo1=' lim_flx: fsolar : ', tab2d_2=fnsolar, clinfo2=' fnsolar : ')
-         CALL prt_ctl(tab2d_1=fmass , clinfo1=' lim_flx: fmass  : ', tab2d_2=fsalt  , clinfo2=' fsalt   : ')
-         CALL prt_ctl(tab2d_1=ftaux , clinfo1=' lim_flx: ftaux  : ', tab2d_2=ftauy  , clinfo2=' ftauy   : ')
-         CALL prt_ctl(tab2d_1=freeze, clinfo1=' lim_flx: freeze : ', tab2d_2=tn_ice , clinfo2=' tn_ice  : ')
-      ENDIF 
-   
-    END SUBROUTINE lim_flx_2
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option :        Dummy module       NO LIM 2.0 sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_flx_2       ! Dummy routine
-   END SUBROUTINE lim_flx_2
-#endif 
-
-   !!======================================================================
-END MODULE limflx_2
Index: trunk/NEMO/LIM_SRC_2/limhdf_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limhdf_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limhdf_2.F90	(revision 888)
@@ -6,5 +6,5 @@
 #if defined key_lim2
    !!----------------------------------------------------------------------
-   !!   'key_lim2'  i                                 LIM 2.0 sea-ice model
+   !!   'key_lim2'                                    LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
    !!   lim_hdf_2  : diffusion trend on sea-ice variable
@@ -34,5 +34,5 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header$ 
+   !! $ Id: $
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    !!----------------------------------------------------------------------
Index: trunk/NEMO/LIM_SRC_2/limistate_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limistate_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limistate_2.F90	(revision 888)
@@ -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_lim2
@@ -18,10 +20,9 @@
    USE phycst
    USE ocfzpt
-   USE oce             ! dynamics and tracers variables      !!gm used???
-   USE dom_oce                                                     !!gm used???
    USE par_ice_2       ! ice parameters
    USE ice_oce         ! ice variables
    USE dom_ice_2
    USE lbclnk
+   USE oce
    USE ice_2
    USE iom
@@ -47,5 +48,5 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
-   !! $Header$ 
+   !! $ Id: $
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -66,46 +67,8 @@
       REAL(wp), DIMENSION(jpi,jpj) ::   ztn   ! workspace
       !--------------------------------------------------------------------
-
-       CALL lim_istate_init_2   !  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_2     !  reading the initials parameters of the ice
+
+      IF( .NOT. ln_limini ) THEN  
          
          ! Initialisation at tn or -2 if ice
@@ -116,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
@@ -152,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(:,:)
@@ -191,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. )
@@ -197,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_2
 
@@ -209,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
@@ -222,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_2 : 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 /) )
@@ -268,8 +210,8 @@
 
             CALL iom_close( inum_ice)
-            
+            !
          ENDIF
       ENDIF
-      !
+      !     
    END SUBROUTINE lim_istate_init_2
 
Index: trunk/NEMO/LIM_SRC_2/limmsh_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limmsh_2.F90	(revision 888)
@@ -25,5 +25,5 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header$ 
+   !! $ Id: $
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    !!----------------------------------------------------------------------
Index: trunk/NEMO/LIM_SRC_2/limrhg_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limrhg_2.F90	(revision 888)
@@ -4,35 +4,42 @@
    !!   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
+   !!            " "  !  06-08  (G. Madec)  surface module, ice-stress at I-point
+   !!            " "  !  09-09  (G. Madec)  Huge verctor optimisation
+   !!----------------------------------------------------------------------
 #if defined key_lim2
    !!----------------------------------------------------------------------
    !!   'key_lim2'                                    LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
    !!   lim_rhg_2   : computes ice velocities
    !!----------------------------------------------------------------------
-   !! * Modules used
-   USE phycst
-   USE par_oce
-   USE ice_oce         ! ice variables
-   USE dom_ice_2
-   USE ice_2
-   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_2      ! domaine: ice variables
+   USE phycst         ! physical constant
+   USE ice_2          ! 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_2  ! routine called by lim_dyn_2
-
-   !! * Module variables
-   REAL(wp)  ::           &  ! constant values
-      rzero   = 0.e0   ,  &
-      rone    = 1.e0
-   !!----------------------------------------------------------------------
-   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   PUBLIC   lim_rhg_2 ! 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) 
+   !! $ Id: $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -48,47 +55,79 @@
       !!  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
       !!
-      !! 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
+      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
+
       !!-------------------------------------------------------------------
-      ! * Arguments
-      INTEGER, INTENT(in) ::   k_j1 ,  &  ! southern j-index for ice computation
-         &                     k_jpj      ! northern j-index for ice computation
-
-      ! * Local variables
-      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
-      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(:,:)
+      !  ----------------
+      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            |
@@ -96,15 +135,16 @@
       !-------------------------------------------------------------------
 
+!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) )
-#if defined key_lim_cp1 && defined key_coupled
-            zb1(ji,jj)    = tms(ji,jj) * gtaux(ji,jj) * ( 1.0 - frld(ji,jj) )
-            zb2(ji,jj)    = tms(ji,jj) * gtauy(ji,jj) * ( 1.0 - frld(ji,jj) )
-#else
-            zb1(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )
-            zb2(ji,jj)    = tms(ji,jj) * ( 1.0 - frld(ji,jj) )
-#endif
+!!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) )
          END DO
       END DO
@@ -117,13 +157,13 @@
          
       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 )
 
-            zt11 = tms(ji  ,jj  ) * frld(ji  ,jj  ) 
-            zt12 = tms(ji-1,jj  ) * frld(ji-1,jj  ) 
-            zt21 = tms(ji  ,jj-1) * frld(ji  ,jj-1) 
-            zt22 = tms(ji-1,jj-1) * frld(ji-1,jj-1)
+            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.
@@ -131,21 +171,15 @@
                &             + zt21 * wght(ji,jj,2,1) + zt22 * wght(ji,jj,1,1) ) * zusw
 
-            ! 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
+            ! 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.
-#if defined key_lim_cp1 && defined key_coupled
-            ztagnx = ( 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
-            ztagny = ( 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
-#else
-            ztagnx = ( 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 * gtaux(ji,jj)
-            ztagny = ( 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 * gtauy(ji,jj)
-#endif
+            ! always provide stress 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
@@ -161,13 +195,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_oce(ji,jj) - zgphsx
-            za2ct(ji,jj) = ztagny + zcorl(ji,jj) * u_oce(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. )
 
 
@@ -182,470 +212,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_oce(ji,jj)
-                 zvr     = v_ice(ji,jj) - v_oce(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_oce(ji,jj) - zsang * v_oce(ji,jj) )
-
-                  za2(ji,jj) =  zmassdt * zv0(ji,jj) - zcorlal * zu0(ji,jj) + za2ct(ji,jj)   &
-                     &        + za * ( cangvg * v_oce(ji,jj) + zsang * u_oce(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_2/limrst_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limrst_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limrst_2.F90	(revision 888)
@@ -17,6 +17,6 @@
    !!----------------------------------------------------------------------
    USE ice_2
-   USE dom_oce
-   USE ice_oce         ! ice variables
+   USE sbc_oce
+   USE sbc_ice
    USE daymod
 
@@ -27,7 +27,7 @@
    PRIVATE
 
-   PUBLIC   lim_rst_opn_2     ! routine called by ??? module
-   PUBLIC   lim_rst_write_2   ! routine called by ??? module
-   PUBLIC   lim_rst_read_2    ! routine called by ??? module
+   PUBLIC   lim_rst_opn_2     ! routine called by sbcice_lim_2.F90
+   PUBLIC   lim_rst_write_2   ! routine called by sbcice_lim_2.F90
+   PUBLIC   lim_rst_read_2    ! routine called by iceini_2.F90
 
    LOGICAL, PUBLIC ::   lrst_ice         !: logical to control the ice restart write 
@@ -36,5 +36,5 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
-   !! $Header: /home/opalod/NEMOCVSROOT/NEMO/LIM_SRC/limrst.F90,v 1.15 2007/06/29 14:54:06 opalod Exp $ 
+   !! $ Id: $
    !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
@@ -57,7 +57,7 @@
       
       ! to get better performances with NetCDF format:
-      ! we open and define the ice restart file one ice time step before writing the data (-> at nitrst - 2*nfice + 1)
-      ! except if we write ice restart files every ice time step or if an ice restart file was writen at nitend - 2*nfice + 1
-      IF( kt == nitrst - 2*nfice + 1 .OR. nstock == nfice .OR. ( kt == nitend - nfice + 1 .AND. .NOT. lrst_ice ) ) THEN
+      ! we open and define the ice restart file one ice time step before writing the data (-> at nitrst - 2*nn_fsbc + 1)
+      ! except if we write ice restart files every ice time step or if an ice restart file was writen at nitend - 2*nn_fsbc + 1
+      IF( kt == nitrst - 2*nn_fsbc + 1 .OR. nstock == nn_fsbc .OR. ( kt == nitend - nn_fsbc + 1 .AND. .NOT. lrst_ice ) ) THEN
          ! beware of the format used to write kt (default is i8.8, that should be large enough...)
          IF( nitrst > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst
@@ -72,7 +72,7 @@
             CASE DEFAULT         ;   WRITE(numout,*) '             open ice restart NetCDF file: '//clname
             END SELECT
-            IF( kt == nitrst - 2*nfice + 1 ) THEN   
-               WRITE(numout,*)         '             kt = nitrst - 2*nfice + 1 = ', kt,' date= ', ndastp
-            ELSE   ;   WRITE(numout,*) '             kt = '                       , kt,' date= ', ndastp
+            IF( kt == nitrst - 2*nn_fsbc + 1 ) THEN   
+               WRITE(numout,*)         '             kt = nitrst - 2*nn_fsbc + 1 = ', kt,' date= ', ndastp
+            ELSE   ;   WRITE(numout,*) '             kt = '                         , kt,' date= ', ndastp
             ENDIF
          ENDIF
@@ -90,15 +90,15 @@
       !! ** 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   ! ice restarts are written at kt == nitrst - nfice + 1
+      INTEGER, INTENT(in) ::   kt   ! number of iteration
+      !
+      INTEGER ::   iter   ! kt + nn_fsbc -1
+      !!----------------------------------------------------------------------
+
+      iter = kt + nn_fsbc - 1   ! ice restarts are written at kt == nitrst - nn_fsbc + 1
 
       IF( iter == nitrst ) THEN
          IF(lwp) WRITE(numout,*)
          IF(lwp) WRITE(numout,*) 'lim_rst_write_2 : write ice restart file  kt =', kt
-         IF(lwp) WRITE(numout,*) '~~~~~~~'         
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~'
       ENDIF
 
@@ -106,6 +106,5 @@
       ! ------------------ 
       !                                                                     ! 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 
@@ -119,8 +118,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 (:,:)   )
@@ -175,6 +172,5 @@
       !! ** purpose  :   read of sea-ice variable restart in a netcdf file
       !!----------------------------------------------------------------------
-      !
-      REAL(wp) ::   zfice, ziter
+      REAL(wp) ::   ziter
       !!----------------------------------------------------------------------
 
@@ -182,12 +178,11 @@
          WRITE(numout,*)
          WRITE(numout,*) 'lim_rst_read_2 : 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
 
@@ -196,8 +191,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' )
@@ -213,8 +204,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_2/limsbc_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 888)
+++ trunk/NEMO/LIM_SRC_2/limsbc_2.F90	(revision 888)
@@ -0,0 +1,251 @@
+MODULE limsbc_2
+   !!======================================================================
+   !!                       ***  MODULE limsbc_2   ***
+   !!           computation of the flux at the sea ice/ocean interface
+   !!======================================================================
+   !! History : 00-01 (H. Goosse) Original code
+   !!           02-07 (C. Ethe, G. Madec) re-writing F90
+   !!           06-07 (G. Madec) surface module
+   !!----------------------------------------------------------------------
+#if defined key_lim2
+   !!----------------------------------------------------------------------
+   !!   'key_lim2'                                    LIM 2.0 sea-ice model
+   !!----------------------------------------------------------------------
+   !!----------------------------------------------------------------------
+   !!   lim_sbc_2  : flux at the ice / ocean interface
+   !!----------------------------------------------------------------------
+   USE par_oce          ! ocean parameters
+   USE dom_oce          ! ocean domain
+   USE sbc_ice          ! surface boundary condition
+   USE sbc_oce          ! surface boundary condition
+   USE phycst           ! physical constants
+   USE ocfzpt           ! surface ocean freezing point
+   USE ice_oce          ! sea-ice variable
+   USE ice_2            ! LIM sea-ice variables
+   USE iceini_2         ! ???
+   USE dynspg_oce       ! choice of the surface pressure gradient scheme 
+
+   USE lbclnk           ! ocean lateral boundary condition
+   USE in_out_manager   ! I/O manager
+   USE albedo           ! albedo parameters
+   USE prtctl           ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC lim_sbc_2     ! called by sbc_ice_lim_2
+
+   REAL(wp)  ::   epsi16 = 1.e-16  ! constant values
+   REAL(wp)  ::   rzero  = 0.e0    
+   REAL(wp)  ::   rone   = 1.e0
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!   LIM 2.0,  UCL-LOCEAN-IPSL (2006) 
+   !! $ Id: $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
+CONTAINS
+
+   SUBROUTINE lim_sbc_2( kt )
+      !!-------------------------------------------------------------------
+      !!                ***  ROUTINE lim_sbc_2 ***
+      !!  
+      !! ** Purpose : Update surface ocean boundary condition over areas
+      !!      that are at least partially covered by sea-ice
+      !!         
+      !! ** Action  : - comput. of the momentum, heat and freshwater/salt
+      !!      fluxes at the ice-ocean interface.
+      !!              - Update 
+      !!     
+      !! ** 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.
+      !!              Tartinville et al. 2001 Ocean Modelling, 3, 95-108.
+      !!---------------------------------------------------------------------
+      INTEGER ::   kt    ! number of iteration
+      !!
+      INTEGER  ::   ji, jj           ! dummy loop indices
+      INTEGER  ::   ifvt, i1mfr, idfr               ! some switches
+      INTEGER  ::   iflt, ial, iadv, ifral, ifrdv
+      REAL(wp) ::   zqsr  , zqns     ! solar & non solar heat flux
+      REAL(wp) ::   zinda            ! switch for testing the values of ice concentration
+      REAL(wp) ::   zfons            ! salt exchanges at the ice/ocean interface
+      REAL(wp) ::   zemp             ! freshwater exchanges at the ice/ocean interface
+      REAL(wp) ::   zfrldu, zfrldv   ! lead fraction at U- & V-points
+      REAL(wp) ::   zutau , zvtau    ! lead fraction at U- & V-points
+      REAL(wp) ::   zu_io , zv_io    ! 2 components of the ice-ocean velocity
+#if defined key_coupled    
+      REAL(wp), DIMENSION(jpi,jpj) ::   zalb     ! albedo of ice under overcast sky
+      REAL(wp), DIMENSION(jpi,jpj) ::   zalbp    ! albedo of ice under clear sky
+#endif
+      REAL(wp) ::   zsang, zmod, zfm
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztio_u, ztio_v   ! ocean stress below sea-ice
+
+      !!---------------------------------------------------------------------
+     
+      IF( kt == nit000 ) THEN
+         IF(lwp) WRITE(numout,*)
+         IF(lwp) WRITE(numout,*) 'lim_sbc_2 : LIM 2.0 sea-ice - surface boundary condition'
+         IF(lwp) WRITE(numout,*) '~~~~~~~~~   '
+      ENDIF
+
+      !------------------------------------------!
+      !      heat flux at the ocean surface      !
+      !------------------------------------------!
+
+!!gm
+!!gm CAUTION   
+!!gm re-verifies the non solar expression, especially over open ocen
+!!gm
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zinda   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - pfrld(ji,jj) )   ) )
+            ifvt    = zinda * MAX( rzero , SIGN( rone,  - phicif(ji,jj)           ) )
+            i1mfr   = 1.0   - MAX( rzero , SIGN( rone, - ( 1.0 - frld(ji,jj) )    ) )
+            idfr    = 1.0   - MAX( rzero , SIGN( rone, frld(ji,jj) - pfrld(ji,jj) ) )
+            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt )
+            ial     = ifvt   * i1mfr + ( 1 - ifvt ) * idfr
+            iadv    = ( 1  - i1mfr ) * zinda
+            ifral   = ( 1  - i1mfr * ( 1 - ial ) )   
+            ifrdv   = ( 1  - ifral * ( 1 - ial ) ) * iadv 
+            !   computation the solar flux at ocean surface
+            zqsr    = pfrld(ji,jj) * qsr(ji,jj)  + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj)
+            !  computation the non solar heat flux at ocean surface
+            zqns    =  - ( 1. - thcm(ji,jj) ) * zqsr   &   ! part of the solar energy used in leads
+               &       + iflt    * ( fscmbq(ji,jj) + ffltbif(ji,jj) )                            &
+               &       + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) / rdt_ice   &
+               &       + ifrdv   * ( qfvbq(ji,jj) + qdtcn(ji,jj) ) / rdt_ice
+
+            fsbbq(ji,jj) = ( 1.0 - ( ifvt + iflt ) ) * fscmbq(ji,jj)     ! ???
+            
+            qsr  (ji,jj) = zqsr                                          ! solar heat flux 
+            qns  (ji,jj) = zqns - fdtcn(ji,jj)                           ! non solar heat flux
+         END DO
+      END DO
+  
+      !------------------------------------------!
+      !      mass flux at the ocean surface      !
+      !------------------------------------------!
+
+!!gm
+!!gm CAUTION   
+!!gm re-verifies the emp & emps expression, especially the absence of 1-frld on zfm
+!!gm
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            
+            !  computing freshwater exchanges at the ice/ocean interface
+            zemp = + emp(ji,jj)     *         frld(ji,jj)      &   !  e-p budget over open ocean fraction 
+               &   - tprecip(ji,jj) * ( 1. -  frld(ji,jj) )    &   !  liquid precipitation reaches directly the ocean
+               &   + sprecip(ji,jj) * ( 1. - pfrld(ji,jj) )    &   !  taking into account change in ice cover within the time step
+               &   + rdmsnif(ji,jj) / rdt_ice                      !  freshwaterflux due to snow melting 
+               !                                                   !  ice-covered fraction:
+
+            !  computing salt exchanges at the ice/ocean interface
+            zfons =  ( soce - sice ) * ( rdmicif(ji,jj) / rdt_ice ) 
+            
+            !  converting the salt flux from ice to a freshwater flux from ocean
+            zfm  = zfons / ( sss_m(ji,jj) + epsi16 )
+            
+            emps(ji,jj) = zemp + zfm      ! surface ocean concentration/dilution effect (use on SSS evolution)
+            emp (ji,jj) = zemp            ! surface ocean volume flux (use on sea-surface height evolution)
+
+         END DO
+      END DO
+
+      IF( lk_dynspg_rl )    emp (:,:) = emps(:,:)      ! rigid-lid formulation : emp = emps
+
+      !------------------------------------------!
+      !    momentum flux at the ocean surface    !
+      !------------------------------------------!
+
+      IF ( ln_limdyn ) THEN                        ! Update the stress over ice-over area (only in ice-dynamic case)
+         !                                         ! otherwise the atmosphere-ocean stress is used everywhere
+
+         ! ... ice stress over ocean with a ice-ocean rotation angle (at I-point)
+!CDIR NOVERRCHK
+         DO jj = 1, jpj
+!CDIR NOVERRCHK
+            DO ji = 1, jpi
+               ! ... change the cosinus angle sign in the south hemisphere
+               zsang  = SIGN(1.e0, gphif(ji,jj) ) * sangvg
+               ! ... ice velocity relative to the ocean
+               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) = zmod * ( cangvg * zu_io - zsang * zv_io )
+               ztio_v(ji,jj) = zmod * ( cangvg * zv_io + zsang * zu_io )
+               ! 
+            END DO
+         END DO
+
+         DO jj = 2, jpjm1
+            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) )
+               zvtau  = 0.5 * ( ztio_v(ji,jj+1) + ztio_v(ji+1,jj+1) )
+               ! ... open-ocean (lead) fraction at U- & V-points (from T-point values)
+               zfrldu = 0.5 * ( frld (ji,jj) + frld (ji+1,jj  ) )
+               zfrldv = 0.5 * ( frld (ji,jj) + frld (ji  ,jj+1) )
+               ! update surface ocean stress
+               utau(ji,jj) = zfrldu * utau(ji,jj) + ( 1. - zfrldu ) * zutau
+               vtau(ji,jj) = zfrldv * vtau(ji,jj) + ( 1. - zfrldv ) * zvtau
+               !
+            END DO
+         END DO
+
+         ! boundary condition on the stress (utau,vtau)
+         CALL lbc_lnk( utau, 'U', -1. )
+         CALL lbc_lnk( vtau, 'V', -1. )
+
+      ENDIF
+
+      !-----------------------------------------------!
+      !   Storing the transmitted variables           !
+      !-----------------------------------------------!
+
+      freeze(:,:) = 1.0 - frld(:,:)       ! Sea ice cover            
+      tn_ice(:,:) = sist(:,:)             ! Ice surface temperature                      
+
+#if defined key_coupled            
+      !------------------------------------------------!
+      !    Computation of snow/ice and ocean albedo    !
+      !------------------------------------------------!
+      zalb  (:,:) = 0.e0
+      zalbp (:,:) = 0.e0
+
+      CALL albedo_ice( sist, hicif, hsnif, zalbp, zalb )
+
+      alb_ice(:,:) =  0.5 * zalbp(:,:) + 0.5 * zalb (:,:)   ! Ice albedo (mean clear and overcast skys)
+#endif
+
+      IF(ln_ctl) THEN
+         CALL prt_ctl(tab2d_1=qsr   , clinfo1=' lim_sbc: qsr    : ', tab2d_2=qns   , clinfo2=' qns     : ')
+         CALL prt_ctl(tab2d_1=emp   , clinfo1=' lim_sbc: emp    : ', tab2d_2=emps  , clinfo2=' emps    : ')
+         CALL prt_ctl(tab2d_1=utau  , clinfo1=' lim_sbc: utau   : ', mask1=umask,   &
+            &         tab2d_2=vtau  , clinfo2=' vtau    : '        , mask2=vmask )
+         CALL prt_ctl(tab2d_1=freeze, clinfo1=' lim_sbc: freeze : ', tab2d_2=tn_ice, clinfo2=' tn_ice  : ')
+      ENDIF 
+   
+    END SUBROUTINE lim_sbc_2
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option :        Dummy module       NO LIM 2.0 sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE lim_sbc_2         ! Dummy routine
+   END SUBROUTINE lim_sbc_2
+#endif 
+
+   !!======================================================================
+END MODULE limsbc_2
Index: trunk/NEMO/LIM_SRC_2/limtab_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limtab_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limtab_2.F90	(revision 888)
@@ -21,6 +21,6 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! $ Id: $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 CONTAINS
Index: trunk/NEMO/LIM_SRC_2/limthd_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limthd_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limthd_2.F90	(revision 888)
@@ -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_lim2
    !!----------------------------------------------------------------------
@@ -18,5 +22,6 @@
    USE ice_2           ! LIM sea-ice variables
    USE ice_oce         ! sea-ice/ocean variables
-   USE flx_oce         ! sea-ice/ocean forcings variables 
+   USE sbc_oce         ! 
+   USE sbc_ice         ! 
    USE thd_ice_2       ! LIM thermodynamic sea-ice variables
    USE dom_ice_2       ! LIM sea-ice domain
@@ -30,14 +35,11 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_thd_2       ! called by lim_step_2
-
-   !! * 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_2  ! 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
@@ -46,6 +48,6 @@
    !!-------- -------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
+   !! $ Id: $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
    !!----------------------------------------------------------------------
 
@@ -68,13 +70,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 +89,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_2  ! Initialization (first time-step only)
+      IF( kt == nit000 )   CALL lim_thd_init_2  ! Initialization (first time-step only)
    
       !-------------------------------------------!
@@ -173,4 +168,6 @@
       !--------------------------------------------------------------------------
 
+      sst_m(:,:) = sst_m(:,:) + rt0
+
       !CDIR NOVERRCHK
       DO jj = 1, jpj
@@ -188,11 +185,11 @@
             !  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) ) 
+            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 +203,5 @@
             
             !  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 )
+            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( tfu(ji,jj) - sst_m(ji,jj) ) * ( 1 - zinda )
             
             !  calculate oceanic heat flux.
@@ -216,4 +213,5 @@
       END DO
       
+      sst_m(:,:) = sst_m(:,:) - rt0
       
       !         Select icy points and fulfill arrays for the vectorial grid.
@@ -258,5 +256,5 @@
          CALL tab_2d_1d_2( nbpb, fr1_i0_1d  (1:nbpb)     , fr1_i0     , jpi, jpj, npb(1:nbpb) )
          CALL tab_2d_1d_2( nbpb, fr2_i0_1d  (1:nbpb)     , fr2_i0     , jpi, jpj, npb(1:nbpb) )
-         CALL tab_2d_1d_2( nbpb, qnsr_ice_1d(1:nbpb)     , qnsr_ice   , jpi, jpj, npb(1:nbpb) )
+         CALL tab_2d_1d_2( nbpb, qns_ice_1d (1:nbpb)     , qns_ice    , jpi, jpj, npb(1:nbpb) )
 #if ! defined key_coupled
          CALL tab_2d_1d_2( nbpb, qla_ice_1d (1:nbpb)     , qla_ice    , jpi, jpj, npb(1:nbpb) )
@@ -404,5 +402,5 @@
          CALL prt_ctl(tab2d_1=qstoif, clinfo1=' lim_thd: qstoif  : ', tab2d_2=fsbbq , clinfo2=' fsbbq  : ')
       ENDIF
-
+       !
     END SUBROUTINE lim_thd_2
 
@@ -419,7 +417,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_2/limthd_lac_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limthd_lac_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limthd_lac_2.F90	(revision 888)
@@ -30,5 +30,5 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header$ 
+   !! $ Id: $
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    !!----------------------------------------------------------------------
Index: trunk/NEMO/LIM_SRC_2/limthd_zdf_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limthd_zdf_2.F90	(revision 888)
@@ -4,7 +4,11 @@
    !!                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_lim2
    !!----------------------------------------------------------------------
    !!   'key_lim2'                                    LIM 2.0 sea-ice model
+   !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
    !!   lim_thd_zdf_2 : vertical accr./abl. and lateral ablation of sea ice
@@ -22,18 +26,16 @@
    PRIVATE
 
-   !! * Routine accessibility
-   PUBLIC lim_thd_zdf_2      ! called by lim_thd_2
-
-   !! * Module variables
-   REAL(wp)  ::           &  ! constant values
-      epsi20 = 1.e-20  ,  &
-      epsi13 = 1.e-13  ,  &
-      zzero  = 0.e0    ,  &
-      zone   = 1.e0
+   PUBLIC   lim_thd_zdf_2        ! called by lim_thd_2
+
+   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) 
-   !! $Header$ 
-   !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
-   !!----------------------------------------------------------------------
+   !! $ Id: $
+   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+
 CONTAINS
 
@@ -64,23 +66,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 +88,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 +99,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 +104,4 @@
           , zsupdiag  &    ! of the temperatures inside the snow-ice system
           , zsmbr          ! second member
-
        REAL(wp) :: & 
           zhsu     &     ! thickness of surface layer
@@ -130,5 +118,4 @@
           , zb2 , zd2 , zb3 , zd3 &
           , ztint          ! equivalent temperature at the snow-ice interface
-
        REAL(wp) :: & 
           zexp      &     ! exponential function of the ice thickness
@@ -148,5 +135,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 +148,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 +166,4 @@
        !     Computation of energies due to surface and bottom melting 
        !-----------------------------------------------------------------------
-       
        
        DO ji = kideb , kiut
@@ -201,5 +185,4 @@
        END DO
 
-
        !-------------------------------------------
        !  2. Calculate some intermediate variables.  
@@ -265,5 +248,4 @@
        !     - qstbif_1d, total energy stored in brine pockets (updating)
        !-------------------------------------------------------------------
-
 
        DO ji = kideb , kiut
@@ -288,5 +270,4 @@
        END DO
 
-
        !--------------------------------------------------------------------------------
        !  4. Computation of the surface temperature : determined by considering the 
@@ -333,6 +314,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 +341,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 +347,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 +709,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 +771,18 @@
           qstbif_1d(ji) = zdrfrl2 * qstbif_1d(ji)
           frld_1d(ji)    = zfrld_1d(ji)
-
-       END DO
-       
+          !
+       END DO
+       ! 
     END SUBROUTINE lim_thd_zdf_2
+
 #else
-   !!======================================================================
-   !!                       ***  MODULE limthd_zdf_2   ***
-   !!                              no sea ice model  
-   !!======================================================================
+   !!----------------------------------------------------------------------
+   !!   Default Option                                     NO sea-ice model
+   !!----------------------------------------------------------------------
 CONTAINS
    SUBROUTINE lim_thd_zdf_2          ! Empty routine
    END SUBROUTINE lim_thd_zdf_2
 #endif
- END MODULE limthd_zdf_2
+
+   !!======================================================================
+END MODULE limthd_zdf_2
Index: trunk/NEMO/LIM_SRC_2/limtrp_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limtrp_2.F90	(revision 888)
@@ -30,5 +30,5 @@
 
    !! * Routine accessibility
-   PUBLIC lim_trp_2     ! called by ice_step_2
+   PUBLIC lim_trp_2     ! called by sbc_ice_lim_2
 
    !! * Shared module variables
@@ -48,5 +48,5 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0,  UCL-LOCEAN-IPSL (2005) 
-   !! $Header$ 
+   !! $ Id: $
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt 
    !!----------------------------------------------------------------------
@@ -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
@@ -128,5 +128,5 @@
          IF (lk_mpp ) CALL mpp_max(zcfl)
 
-         IF ( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl
+         IF ( zcfl > 0.5 .AND. lwp )   WRITE(numout,*) 'lim_trp_2 : violation of cfl criterion the ',nday,'th day, cfl = ',zcfl
 
          ! content of properties
Index: trunk/NEMO/LIM_SRC_2/limwri_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limwri_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limwri_2.F90	(revision 888)
@@ -9,5 +9,5 @@
 #if defined key_lim2
    !!----------------------------------------------------------------------
-   !!   'key_lim2'  i                                 LIM 2.0 sea-ice model
+   !!   'key_lim2'                                    LIM 2.0 sea-ice model
    !!----------------------------------------------------------------------
    !!----------------------------------------------------------------------
@@ -15,20 +15,22 @@
    !!   lim_wri_init_2 : initialization and namelist read
    !!----------------------------------------------------------------------
-   USE ioipsl
-   USE dianam    ! build name of file (routine)
    USE phycst
    USE dom_oce
    USE daymod
-   USE in_out_manager
    USE ice_oce         ! ice variables
-   USE flx_oce
+   USE sbc_oce
+   USE sbc_ice
    USE dom_ice_2
    USE ice_2
+
    USE lbclnk
+   USE dianam          ! build name of file (routine)
+   USE in_out_manager
+   USE ioipsl
 
    IMPLICIT NONE
    PRIVATE
 
-   PUBLIC   lim_wri_2        ! routine called by lim_step_2.F90
+   PUBLIC   lim_wri_2      ! routine called by sbc_ice_lim_2
 
    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)
-   !! $Header$
+   !! * Substitutions
+#   include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !!  LIM 2.0, UCL-LOCEAN-IPSL (2006)
+   !! $ Id: $
    !! 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_2 
@@ -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_2 : 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,22 @@
             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)
+            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)
+            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) = qsr_ice(ji,jj)
+            zcmo(ji,jj,18) = qns_ice(ji,jj)
             zcmo(ji,jj,19) = sprecip(ji,jj)
          END DO
@@ -175,5 +171,5 @@
          END DO
          
-         IF ( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN 
+         IF( jf == 7  .OR. jf == 8  .OR. jf == 15 .OR. jf == 16 ) THEN
             CALL lbc_lnk( zfield, 'T', -1. )
          ELSE 
@@ -181,11 +177,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_2
@@ -225,24 +219,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 +262,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_2/limwri_dimg_2.h90
===================================================================
--- trunk/NEMO/LIM_SRC_2/limwri_dimg_2.h90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/limwri_dimg_2.h90	(revision 888)
@@ -2,5 +2,5 @@
    !!----------------------------------------------------------------------
    !!  LIM 2.0, UCL-LOCEAN-IPSL (2005)
-   !! $Header$
+   !! $ Id: $
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
    !!----------------------------------------------------------------------
@@ -82,5 +82,5 @@
 
        zsto     = rdt_ice
-       zout     = nwrite * rdt_ice / nfice
+       zout     = nwrite * rdt_ice / nn_fsbc
        zsec     = 0.
        niter    = 0
@@ -106,23 +106,23 @@
           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)
+
+          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)
+          zcmo(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce
+          zcmo(ji,jj,15) = utaui_ice(ji,jj)
+          zcmo(ji,jj,16) = vtaui_ice(ji,jj)
+          zcmo(ji,jj,17) = qsr_ice(ji,jj)
+          zcmo(ji,jj,18) = qns_ice(ji,jj)
           zcmo(ji,jj,19) = sprecip(ji,jj)
        END DO
@@ -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,23 +149,23 @@
                 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)
-                rcmoy(ji,jj,10) = sss_io(ji,jj)
-
-                rcmoy(ji,jj,11) = fnsolar(ji,jj) + fsolar(ji,jj)
-                rcmoy(ji,jj,12) = fsolar (ji,jj)
-                rcmoy(ji,jj,13) = fnsolar(ji,jj)
+                rcmoy(ji,jj,9)  = sst_m(ji,jj)
+                rcmoy(ji,jj,10) = sss_m(ji,jj)
+
+                rcmoy(ji,jj,11) = qns(ji,jj) + qsr(ji,jj)
+                rcmoy(ji,jj,12) = qsr(ji,jj)
+                rcmoy(ji,jj,13) = qns(ji,jj)
                 ! See thersf for the coefficient
-                rcmoy(ji,jj,14) = - fsalt(ji,jj) * rday * ( sss_io(ji,jj) + epsi16 ) / soce
-                rcmoy(ji,jj,15) = gtaux(ji,jj)
-                rcmoy(ji,jj,16) = gtauy(ji,jj)
-                rcmoy(ji,jj,17) = ( 1.0 - frld(ji,jj) ) * qsr_ice (ji,jj) + frld(ji,jj) * qsr_oce (ji,jj)
-                rcmoy(ji,jj,18) = ( 1.0 - frld(ji,jj) ) * qnsr_ice(ji,jj) + frld(ji,jj) * qnsr_oce(ji,jj)
+                rcmoy(ji,jj,14) = - emps(ji,jj) * rday * ( sss_m(ji,jj) + epsi16 ) / soce
+                rcmoy(ji,jj,15) = utaui_ice(ji,jj)
+                rcmoy(ji,jj,16) = vtaui_ice(ji,jj)
+                rcmoy(ji,jj,17) = qsr_ice(ji,jj)
+                rcmoy(ji,jj,18) = qns_ice(ji,jj)
                 rcmoy(ji,jj,19) = sprecip(ji,jj)
              END DO
@@ -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_2
Index: trunk/NEMO/LIM_SRC_2/par_ice_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/par_ice_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/par_ice_2.F90	(revision 888)
@@ -7,5 +7,5 @@
    !!----------------------------------------------------------------------
    !!  LIM 2.0, UCL-LOCEAN-IPSL (2005)
-   !! $Header$
+   !! $ Id: $
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
    !!----------------------------------------------------------------------
Index: trunk/NEMO/LIM_SRC_2/thd_ice_2.F90
===================================================================
--- trunk/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 833)
+++ trunk/NEMO/LIM_SRC_2/thd_ice_2.F90	(revision 888)
@@ -9,5 +9,5 @@
    !!----------------------------------------------------------------------
    !!   LIM 2.0, UCL-LOCEAN-IPSL (2005)
-   !! $Header$
+   !! $ Id: $
    !! This software is governed by the CeCILL licence see modipsl/doc/NEMO_CeCILL.txt
    !!----------------------------------------------------------------------
@@ -57,5 +57,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
