Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8406)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/ice.F90	(revision 8407)
@@ -74,5 +74,5 @@
    !! in icestp.F90, the routines that compute the changes in the ice     |
    !! state variables are called                                          |
-   !! - lim_dyn : ice dynamics                                            |
+   !! - ice_rhg : ice dynamics                                            |
    !! - lim_trp : ice transport                                           |
    !! - lim_itd_me : mechanical redistribution (ridging and rafting)      |
@@ -208,5 +208,5 @@
    REAL(wp), PUBLIC ::   rn_crhg          !: determines changes in ice strength
    LOGICAL , PUBLIC ::   ln_icestr_bvf    !: use brine volume to diminish ice strength
-                                          ! -- limdyn & limrhg -- !
+                                          ! -- limdyn & icerhg -- !
    REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice
    REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90	(revision 8407)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg.F90	(revision 8407)
@@ -0,0 +1,169 @@
+MODULE icerhg
+   !!======================================================================
+   !!                     ***  MODULE  icerhg  ***
+   !!   Sea-Ice dynamics :  
+   !!======================================================================
+   !! history :  1.0  ! 2002-08  (C. Ethe, G. Madec)  original VP code 
+   !!            3.0  ! 2007-03  (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)  LIM3: EVP-Cgrid
+   !!            3.5  ! 2011-02  (G. Madec) dynamical allocation
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3' :                                 LIM3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!    ice_rhg      : computes ice velocities
+   !!    ice_rhg_init : initialization and namelist read
+   !!----------------------------------------------------------------------
+   USE phycst           ! physical constants
+   USE dom_oce          ! ocean space and time domain
+   USE ice              ! LIM-3 variables
+   USE icerhg_evp       ! EVP rheology
+   USE limcons          ! conservation tests
+   USE limctl           ! control prints
+   USE limvar
+   !
+   USE lbclnk           ! lateral boundary conditions - MPP exchanges
+   USE lib_mpp          ! MPP library
+   USE in_out_manager   ! I/O manager
+   USE lib_fortran      ! glob_sum
+   USE timing           ! Timing
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_rhg        ! routine called by icestp.F90
+   PUBLIC   ice_rhg_init   ! routine called by icestp.F90
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
+   !! $Id: icerhg.F90 8378 2017-07-26 13:55:59Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_rhg( kt )
+      !!-------------------------------------------------------------------
+      !!               ***  ROUTINE ice_rhg  ***
+      !!               
+      !! ** Purpose :   compute ice velocity
+      !!
+      !! ** Action  : comupte - ice velocity (u_ice, v_ice)
+      !!                      - 3 components of the stress tensor (stress1_i, stress2_i, stress12_i)
+      !!                      - shear, divergence and delta (shear_i, divu_i, delta_i)
+      !!--------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   kt     ! number of iteration
+      !!
+      INTEGER  :: jl ! dummy loop indices
+      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
+      !!--------------------------------------------------------------------
+
+      IF( nn_timing == 1 )  CALL timing_start('icerhg')
+
+      CALL lim_var_agg(1)   ! aggregate ice categories
+      !
+      ! conservation test
+      IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+      
+      ! Landfast ice parameterization: define max bottom friction
+      tau_icebfr(:,:) = 0._wp
+      IF( ln_landfast ) THEN
+         DO jl = 1, jpl
+            WHERE( ht_i(:,:,jl) > ht_n(:,:) * rn_gamma )  tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
+         END DO
+      ENDIF
+      
+      ! -----------------------
+      ! Rheology (ice dynamics)
+      ! -----------------------   
+      IF( nn_limdyn /= 0 ) THEN                          ! -- Ice dynamics
+
+         CALL ice_rhg_evp( stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i )
+
+      ELSE
+
+         u_ice(:,:) = rn_uice * umask(:,:,1)             !     or prescribed velocity
+         v_ice(:,:) = rn_vice * vmask(:,:,1)
+         !!CALL RANDOM_NUMBER(u_ice(:,:))
+         !!CALL RANDOM_NUMBER(v_ice(:,:))
+
+      ENDIF
+      !
+      ! conservation test
+      IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'icerhg', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
+
+      ! Control prints
+      IF( ln_ctl )       CALL lim_prt3D( 'icerhg' )
+      !
+      IF( nn_timing == 1 )  CALL timing_stop('icerhg')
+
+   END SUBROUTINE ice_rhg
+
+
+   SUBROUTINE ice_rhg_init
+      !!-------------------------------------------------------------------
+      !!                  ***  ROUTINE ice_rhg_init  ***
+      !!
+      !! ** Purpose : Physical constants and parameters linked to the ice
+      !!      dynamics
+      !!
+      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
+      !!       parameter values called at the first timestep (nit000)
+      !!
+      !! ** input   :   Namelist namicedyn
+      !!-------------------------------------------------------------------
+      INTEGER  ::   ios                 ! Local integer output status for namelist read
+      NAMELIST/namicedyn/ nn_limadv, nn_limadv_ord,                                &
+         &                nn_icestr, rn_pe_rdg, rn_pstar, rn_crhg, ln_icestr_bvf,  &
+         &                rn_ishlat, rn_cio, rn_creepl, rn_ecc,                    &
+         &                nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr, rn_lfrelax
+      !!-------------------------------------------------------------------
+
+      REWIND( numnam_ice_ref )              ! Namelist namicedyn in reference namelist : Ice dynamics
+      READ  ( numnam_ice_ref, namicedyn, IOSTAT = ios, ERR = 901)
+901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in reference namelist', lwp )
+
+      REWIND( numnam_ice_cfg )              ! Namelist namicedyn in configuration namelist : Ice dynamics
+      READ  ( numnam_ice_cfg, namicedyn, IOSTAT = ios, ERR = 902 )
+902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in configuration namelist', lwp )
+      IF(lwm) WRITE ( numoni, namicedyn )
+      
+      IF(lwp) THEN                        ! control print
+         WRITE(numout,*)
+         WRITE(numout,*) 'ice_rhg_init : ice parameters for ice dynamics '
+         WRITE(numout,*) '~~~~~~~~~~~~'
+         ! limtrp
+         WRITE(numout,*)'    choose the advection scheme (-1=Prather, 0=Ulimate-Macho)   nn_limadv     = ', nn_limadv 
+         WRITE(numout,*)'    choose the order of the scheme (if ultimate)                nn_limadv_ord = ', nn_limadv_ord  
+         ! limitd_me
+         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)         nn_icestr     = ', nn_icestr 
+         WRITE(numout,*)'    Ratio of ridging work to PotEner change in ridging          rn_pe_rdg     = ', rn_pe_rdg 
+         WRITE(numout,*) '   first bulk-rheology parameter                               rn_pstar      = ', rn_pstar
+         WRITE(numout,*) '   second bulk-rhelogy parameter                               rn_crhg       = ', rn_crhg
+         WRITE(numout,*)'    Including brine volume in ice strength comp.                ln_icestr_bvf = ', ln_icestr_bvf
+         ! icerhg_evp
+         WRITE(numout,*) '   lateral boundary condition for sea ice dynamics             rn_ishlat     = ', rn_ishlat
+         WRITE(numout,*) '   drag coefficient for oceanic stress                         rn_cio        = ', rn_cio
+         WRITE(numout,*) '   creep limit                                                 rn_creepl     = ', rn_creepl
+         WRITE(numout,*) '   eccentricity of the elliptical yield curve                  rn_ecc        = ', rn_ecc
+         WRITE(numout,*) '   number of iterations for subcycling                         nn_nevp       = ', nn_nevp
+         WRITE(numout,*) '   ratio of elastic timescale over ice time step               rn_relast     = ', rn_relast
+         WRITE(numout,*) '   Landfast: param (T or F)                                    ln_landfast   = ', ln_landfast
+         WRITE(numout,*) '      T: fraction of ocean depth that ice must reach           rn_gamma      = ', rn_gamma
+         WRITE(numout,*) '      T: maximum bottom stress per unit area of contact        rn_icebfr     = ', rn_icebfr
+         WRITE(numout,*) '      T: relax time scale (s-1) to reach static friction       rn_lfrelax    = ', rn_lfrelax
+      ENDIF
+      !
+      IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  free-slip '
+      ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  no-slip '
+      ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  partial-slip '
+      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  strong-slip '
+      ENDIF
+      !
+   END SUBROUTINE ice_rhg_init
+
+#endif 
+
+   !!======================================================================
+END MODULE icerhg
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90	(revision 8407)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icerhg_evp.F90	(revision 8407)
@@ -0,0 +1,765 @@
+MODULE icerhg_evp
+   !!======================================================================
+   !!                     ***  MODULE  icerhg_evp  ***
+   !!   Ice rheology : sea ice rheology
+   !!======================================================================
+   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
+   !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3
+   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy 
+   !!            3.3  !  2009-05  (G.Garric) addition of the evp cas
+   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation 
+   !!            3.5  !  2012-08  (R. Benshila)  AGRIF
+   !!            3.6  !  2016-06  (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013)
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                      LIM-3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   ice_rhg_evp       : computes ice velocities
+   !!----------------------------------------------------------------------
+   USE phycst         ! Physical constant
+   USE par_oce        ! Ocean parameters
+   USE dom_oce        ! Ocean domain
+   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
+   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
+   USE ice            ! ice variables
+   USE limitd_me      ! ice strength
+   !
+   USE lbclnk         ! Lateral Boundary Condition / MPP link
+   USE lib_mpp        ! MPP library
+   USE in_out_manager ! I/O manager
+   USE prtctl         ! Print control
+   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
+#if defined key_agrif
+   USE agrif_lim3_interp
+#endif
+   USE bdy_oce   , ONLY: ln_bdy 
+   USE bdyice 
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_rhg_evp        ! routine called by icerhg.F90
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
+   !! $Id: icerhg_evp.F90 8378 2017-07-26 13:55:59Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_rhg_evp( stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i )
+      !!-------------------------------------------------------------------
+      !!                 ***  SUBROUTINE ice_rhg_evp  ***
+      !!                          EVP-C-grid
+      !!
+      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
+      !!  stress and sea-surface slope. Ice-ice interaction is described by 
+      !!  a non-linear elasto-viscous-plastic (EVP) law including shear 
+      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).	
+      !!
+      !!  The points in the C-grid look like this, dear reader
+      !!
+      !!                              (ji,jj)
+      !!                                 |
+      !!                                 |
+      !!                      (ji-1,jj)  |  (ji,jj)
+      !!                             ---------   
+      !!                            |         |
+      !!                            | (ji,jj) |------(ji,jj)
+      !!                            |         |
+      !!                             ---------   
+      !!                     (ji-1,jj-1)     (ji,jj-1)
+      !!
+      !! ** Inputs  : - wind forcing (stress), oceanic currents
+      !!                ice total volume (vt_i) per unit area
+      !!                snow total volume (vt_s) per unit area
+      !!
+      !! ** Action  : - compute u_ice, v_ice : the components of the 
+      !!                sea-ice velocity vector
+      !!              - compute delta_i, shear_i, divu_i, which are inputs
+      !!                of the ice thickness distribution
+      !!
+      !! ** Steps   : 1) Compute ice snow mass, ice strength 
+      !!              2) Compute wind, oceanic stresses, mass terms and
+      !!                 coriolis terms of the momentum equation
+      !!              3) Solve the momentum equation (iterative procedure)
+      !!              4) Prevent high velocities if the ice is thin
+      !!              5) Recompute invariants of the strain rate tensor
+      !!                 which are inputs of the ITD, store stress
+      !!                 for the next time step
+      !!              6) Control prints of residual (convergence)
+      !!                 and charge ellipse.
+      !!                 The user should make sure that the parameters
+      !!                 nn_nevp, elastic time scale and rn_creepl maintain stress state
+      !!                 on the charge ellipse for plastic flow
+      !!                 e.g. in the Canadian Archipelago
+      !!
+      !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013
+      !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters)
+      !!              but this solution appears very unstable (see Kimmritz et al 2016)
+      !!
+      !! References : Hunke and Dukowicz, JPO97
+      !!              Bouillon et al., Ocean Modelling 2009
+      !!              Bouillon et al., Ocean Modelling 2013
+      !!-------------------------------------------------------------------
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   stress1_i, stress2_i, stress12_i
+      REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   u_ice, v_ice, shear_i, divu_i, delta_i 
+      !!
+      INTEGER ::   ji, jj       ! dummy loop indices
+      INTEGER ::   jter         ! local integers
+      CHARACTER (len=50) ::   charout
+
+      REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio
+      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling
+      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity
+      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013
+      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass
+      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
+      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                                ! temporary scalars
+
+      REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress
+      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
+      REAL(wp) ::   zintb, zintn                                             ! dummy argument
+      REAL(wp) ::   zfac_x, zfac_y
+      
+      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
+      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses
+      
+      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
+      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
+      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
+      REAL(wp), DIMENSION(jpi,jpj) ::   zpice                           ! array used for the calculation of ice surface slope:
+                                                                             !   ocean surface (ssh_m) if ice is not embedded
+                                                                             !   ice top surface if ice is embedded   
+      REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
+
+      REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice
+
+      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
+      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
+      !!-------------------------------------------------------------------
+
+#if defined key_agrif 
+      CALL agrif_interp_lim3( 'U', 0, nn_nevp )   ! First interpolation of coarse values
+      CALL agrif_interp_lim3( 'V', 0, nn_nevp )
+#endif
+      !
+      !------------------------------------------------------------------------------!
+      ! 0) mask at F points for the ice
+      !------------------------------------------------------------------------------!
+      ! ocean/land mask
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1      ! NO vector opt.
+            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
+         END DO
+      END DO
+      CALL lbc_lnk( zfmask, 'F', 1._wp )
+
+      ! Lateral boundary conditions on velocity (modify zfmask)
+      zwf(:,:) = zfmask(:,:)
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1   ! vector opt.
+            IF( zfmask(ji,jj) == 0._wp ) THEN
+               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
+            ENDIF
+         END DO
+      END DO
+      DO jj = 2, jpjm1
+         IF( zfmask(1,jj) == 0._wp ) THEN
+            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
+         ENDIF
+         IF( zfmask(jpi,jj) == 0._wp ) THEN
+            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
+         ENDIF
+      END DO
+      DO ji = 2, jpim1
+         IF( zfmask(ji,1) == 0._wp ) THEN
+            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
+         ENDIF
+         IF( zfmask(ji,jpj) == 0._wp ) THEN
+            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
+         ENDIF
+      END DO
+      CALL lbc_lnk( zfmask, 'F', 1._wp )
+
+      !------------------------------------------------------------------------------!
+      ! 1) define some variables and initialize arrays
+      !------------------------------------------------------------------------------!
+      zrhoco = rau0 * rn_cio 
+
+      ! ecc2: square of yield ellipse eccenticrity
+      ecc2    = rn_ecc * rn_ecc
+      z1_ecc2 = 1._wp / ecc2
+
+      ! Time step for subcycling
+      zdtevp   = rdt_ice / REAL( nn_nevp )
+      z1_dtevp = 1._wp / zdtevp
+
+      ! alpha parameters (Bouillon 2009)
+      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
+      zalph2 = zalph1 * z1_ecc2
+
+      ! alpha and beta parameters (Bouillon 2013)
+      !!zalph1 = 40.
+      !!zalph2 = 40.
+      !!zbeta  = 3000.
+      !!zbeta = REAL( nn_nevp )   ! close to classical EVP of Hunke (2001)
+
+      z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+      z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+
+      ! Initialise stress tensor 
+      zs1 (:,:) = stress1_i (:,:) 
+      zs2 (:,:) = stress2_i (:,:)
+      zs12(:,:) = stress12_i(:,:)
+
+      ! Ice strength
+      CALL lim_itd_me_icestrength( nn_icestr )
+
+      ! scale factors
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
+            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
+         END DO
+      END DO
+            
+      !
+      !------------------------------------------------------------------------------!
+      ! 2) Wind / ocean stress, mass terms, coriolis terms
+      !------------------------------------------------------------------------------!
+
+      IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
+         !                                            
+         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
+         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
+         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
+         !
+         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
+         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
+         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
+         !
+         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
+         !
+      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
+         zpice(:,:) = ssh_m(:,:)
+      ENDIF
+
+      DO jj = 2, jpjm1
+         DO ji = fs_2, fs_jpim1
+
+            ! ice fraction at U-V points
+            zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+            zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+
+            ! Ice/snow mass at U-V points
+            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
+            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
+            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
+            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+
+            ! Ocean currents at U-V points
+            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
+               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
+            
+            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
+               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
+
+            ! Coriolis at T points (m*f)
+            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
+
+            ! m/dt
+            zmU_t(ji,jj)    = zmassU * z1_dtevp
+            zmV_t(ji,jj)    = zmassV * z1_dtevp
+
+            ! Drag ice-atm.
+            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
+            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
+
+            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
+            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
+            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
+
+            ! masks
+            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
+            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
+
+            ! switches
+            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
+            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
+
+         END DO
+      END DO
+      CALL lbc_lnk( zmf, 'T', 1. )
+      !
+      !------------------------------------------------------------------------------!
+      ! 3) Solution of the momentum equation, iterative procedure
+      !------------------------------------------------------------------------------!
+      !
+      !                                               !----------------------!
+      DO jter = 1 , nn_nevp                           !    loop over jter    !
+         !                                            !----------------------!        
+         IF(ln_ctl) THEN   ! Convergence test
+            DO jj = 1, jpjm1
+               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
+               zv_ice(:,jj) = v_ice(:,jj)
+            END DO
+         ENDIF
+
+         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
+         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
+            DO ji = 1, jpim1
+
+               ! shear at F points
+               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+
+            END DO
+         END DO
+         CALL lbc_lnk( zds, 'F', 1. )
+
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1 ! no vector loop
+
+               ! shear**2 at T points (doc eq. A16)
+               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+              
+               ! divergence at T points
+               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+                  &    ) * r1_e1e2t(ji,jj)
+               zdiv2 = zdiv * zdiv
+               
+               ! tension at T points
+               zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+                  &   ) * r1_e1e2t(ji,jj)
+               zdt2 = zdt * zdt
+               
+               ! delta at T points
+               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
+
+               ! P/delta at T points
+               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
+               
+               ! stress at T points
+               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
+             
+            END DO
+         END DO
+         CALL lbc_lnk( zp_delt, 'T', 1. )
+
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+
+               ! P/delta at F points
+               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
+               
+               ! stress at F points
+               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
+
+            END DO
+         END DO
+         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
+ 
+
+         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1               
+
+               ! U points
+               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
+                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
+                  &                    ) * r1_e2u(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
+                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
+                  &                  ) * r1_e1e2u(ji,jj)
+
+               ! V points
+               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
+                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
+                  &                    ) * r1_e1v(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
+                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
+                  &                  ) * r1_e1e2v(ji,jj)
+
+               ! u_ice at V point
+               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
+                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
+               
+               ! v_ice at U point
+               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
+                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
+
+            END DO
+         END DO
+         !
+         ! --- Computation of ice velocity --- !
+         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
+         !  Bouillon et al. 2009 (eq 34-35) => stable
+         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
+            
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+
+                  ! tau_io/(v_oce - v_ice)
+                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
+
+                  ! Ocean-to-Ice stress
+                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+
+                  ! tau_bottom/v_ice
+                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
+                  zTauB = - tau_icebfr(ji,jj) / zvel
+
+                  ! Coriolis at V-points (energy conserving formulation)
+                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
+                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+
+                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
+
+                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
+                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
+                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
+                  ! Bouillon 2013
+                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
+                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
+                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
+
+               END DO
+            END DO
+            CALL lbc_lnk( v_ice, 'V', -1. )
+            
+#if defined key_agrif
+!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
+            CALL agrif_interp_lim3( 'V' )
+#endif
+            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
+
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+                               
+                  ! tau_io/(u_oce - u_ice)
+                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
+
+                  ! Ocean-to-Ice stress
+                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+
+                  ! tau_bottom/u_ice
+                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
+                  zTauB = - tau_icebfr(ji,jj) / zvel
+
+                  ! Coriolis at U-points (energy conserving formulation)
+                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
+                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
+                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+                  
+                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
+
+                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
+                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
+                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
+                     &           ) * zmaskU(ji,jj)
+                  ! Bouillon 2013
+                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
+                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
+                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
+               END DO
+            END DO
+            CALL lbc_lnk( u_ice, 'U', -1. )
+            
+#if defined key_agrif
+!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
+            CALL agrif_interp_lim3( 'U' )
+#endif
+            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
+
+         ELSE ! odd iterations
+
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+
+                  ! tau_io/(u_oce - u_ice)
+                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
+                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
+
+                  ! Ocean-to-Ice stress
+                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
+
+                  ! tau_bottom/u_ice
+                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
+                  zTauB = - tau_icebfr(ji,jj) / zvel
+
+                  ! Coriolis at U-points (energy conserving formulation)
+                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
+                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
+                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
+
+                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
+
+                  ! landfast switch => 0 = static friction ; 1 = sliding friction
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
+                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
+                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskU(ji,jj)
+                  ! Bouillon 2013
+                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
+                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
+                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
+               END DO
+            END DO
+            CALL lbc_lnk( u_ice, 'U', -1. )
+            
+#if defined key_agrif
+!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
+            CALL agrif_interp_lim3( 'U' )
+#endif
+            IF( ln_bdy ) CALL bdy_ice_dyn( 'U' )
+
+            DO jj = 2, jpjm1
+               DO ji = fs_2, fs_jpim1
+         
+                  ! tau_io/(v_oce - v_ice)
+                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
+                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
+
+                  ! Ocean-to-Ice stress
+                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
+
+                  ! tau_bottom/v_ice
+                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
+                  ztauB = - tau_icebfr(ji,jj) / zvel
+                  
+                  ! Coriolis at V-points (energy conserving formulation)
+                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
+                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
+                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
+
+                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
+                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
+
+                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
+                  
+                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
+                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
+                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
+                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
+                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
+                     &           ) * zmaskV(ji,jj)
+                  ! Bouillon 2013
+                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
+                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
+                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
+               END DO
+            END DO
+            CALL lbc_lnk( v_ice, 'V', -1. )
+            
+#if defined key_agrif
+!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
+            CALL agrif_interp_lim3( 'V' )
+#endif
+            IF( ln_bdy ) CALL bdy_ice_dyn( 'V' )
+
+         ENDIF
+         
+         IF(ln_ctl) THEN   ! Convergence test
+            DO jj = 2 , jpjm1
+               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
+            END DO
+            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
+            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
+         ENDIF
+         !
+         !                                                ! ==================== !
+      END DO                                              !  end loop over jter  !
+      !                                                   ! ==================== !
+      !
+      !------------------------------------------------------------------------------!
+      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 
+      !------------------------------------------------------------------------------!
+      DO jj = 1, jpjm1
+         DO ji = 1, jpim1
+
+            ! shear at F points
+            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+
+         END DO
+      END DO           
+      CALL lbc_lnk( zds, 'F', 1. )
+      
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1 ! no vector loop
+            
+            ! tension**2 at T points
+            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+               &   ) * r1_e1e2t(ji,jj)
+            zdt2 = zdt * zdt
+            
+            ! shear**2 at T points (doc eq. A16)
+            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+               &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+            
+            ! shear at T points
+            shear_i(ji,jj) = SQRT( zdt2 + zds2 )
+
+            ! divergence at T points
+            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+               &            ) * r1_e1e2t(ji,jj)
+            
+            ! delta at T points
+            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
+            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
+            delta_i(ji,jj) = zdelta + rn_creepl * rswitch
+
+         END DO
+      END DO
+      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
+      
+      ! --- Store the stress tensor for the next time step --- !
+      stress1_i (:,:) = zs1 (:,:)
+      stress2_i (:,:) = zs2 (:,:)
+      stress12_i(:,:) = zs12(:,:)
+      !
+
+      !------------------------------------------------------------------------------!
+      ! 5) SIMIP diagnostics
+      !------------------------------------------------------------------------------!
+                           
+      DO jj = 2, jpjm1
+         DO ji = 2, jpim1
+             rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
+
+             ! Stress tensor invariants (normal and shear stress N/m)
+             diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
+             diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
+
+             ! Stress terms of the momentum equation (N/m2)
+             diag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term
+             diag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
+
+             diag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term
+             diag_corstry(ji,jj) = zCory(ji,jj) * rswitch
+
+             diag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term
+             diag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
+           
+             diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
+             diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
+
+             ! 2D ice mass, snow mass, area transport arrays (X, Y)
+             zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
+             zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
+
+             diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
+             diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
+
+             diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
+             diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
+
+             diag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component
+             diag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   ''
+
+         END DO
+      END DO
+
+      CALL lbc_lnk_multi(   diag_sig1   , 'T',  1., diag_sig2   , 'T',  1.,   &
+                 &          diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1.,   &
+                 &          diag_corstrx, 'U', -1., diag_corstry, 'V', -1.,   & 
+                 &          diag_intstrx, 'U', -1., diag_intstry, 'V', -1.    )
+
+      CALL lbc_lnk_multi(   diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1.    )
+
+      CALL lbc_lnk_multi(   diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1., &
+                 &          diag_xatrp    , 'U', -1., diag_ymtrp_ice, 'V', -1., &
+                 &          diag_ymtrp_snw, 'V', -1., diag_yatrp    , 'V', -1.  )
+
+      !
+      !------------------------------------------------------------------------------!
+      ! 6) Control prints of residual and charge ellipse
+      !------------------------------------------------------------------------------!
+      !
+      ! print the residual for convergence
+      IF(ln_ctl) THEN
+         WRITE(charout,FMT="('ice_rhg_evp  : res =',D23.16, ' iter =',I4)") zresm, jter
+         CALL prt_ctl_info(charout)
+         CALL prt_ctl(tab2d_1=u_ice, clinfo1=' ice_rhg_evp  : u_ice :', tab2d_2=v_ice, clinfo2=' v_ice :')
+      ENDIF
+
+      ! print charge ellipse
+      ! This can be desactivated once the user is sure that the stress state
+      ! lie on the charge ellipse. See Bouillon et al. 08 for more details
+      IF(ln_ctl) THEN
+         IF( MOD(kt_ice+nn_fsbc-1,nwrite) == 0 ) THEN
+            WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 1000, kt_ice, 0, 0, ' ch. ell. '
+            CALL prt_ctl_info(charout)
+            DO jj = 2, jpjm1
+               DO ji = 2, jpim1
+                  IF (strength(ji,jj) > 1.0) THEN
+                     zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 
+                     zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) )
+                     WRITE(charout,FMT="('ice_rhg_evp  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
+                     CALL prt_ctl_info(charout)
+                  ENDIF
+               END DO
+            END DO
+            WRITE(charout,FMT="('ice_rhg_evp  :', I4, I6, I1, I1, A10)") 2000, kt_ice, 0, 0, ' ch. ell. '
+            CALL prt_ctl_info(charout)
+         ENDIF
+      ENDIF
+      !
+   END SUBROUTINE ice_rhg_evp
+
+#endif
+
+   !!==============================================================================
+END MODULE icerhg_evp
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8406)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90	(revision 8407)
@@ -32,5 +32,5 @@
    USE phycst          ! Define parameters for the routines
    USE eosbn2          ! equation of state
-   USE limdyn          ! Ice dynamics
+   USE icerhg          ! Ice dynamics
    USE limtrp          ! Ice transport
    USE limthd          ! Ice thermodynamics
@@ -60,5 +60,5 @@
 
    USE bdy_oce   , ONLY: ln_bdy
-   USE bdyice_lim      ! unstructured open boundary data  (bdy_ice_lim routine)
+   USE bdyice          ! unstructured open boundary data
 # if defined key_agrif
    USE agrif_ice
@@ -106,8 +106,5 @@
                                        !                    3 = bulk formulation,
                                        !                    4 = Pure Coupled formulation)
-      !!
       INTEGER  ::   jl                 ! dummy loop index
-      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_os, zalb_cs  ! ice albedo under overcast/clear sky
-      REAL(wp), DIMENSION(jpi,jpj)     ::   zutau_ice, zvtau_ice 
       !!----------------------------------------------------------------------
 
@@ -158,5 +155,5 @@
             !
             IF( nn_limdyn /= 0 ) THEN                          ! -- Ice dynamics
-                                      CALL lim_dyn( kt )       !     rheology  
+                                      CALL ice_rhg( kt )       !     rheology  
             ELSE
                u_ice(:,:) = rn_uice * umask(:,:,1)             !     or prescribed velocity
@@ -176,5 +173,5 @@
          IF( .NOT. Agrif_Root() )     CALL agrif_interp_lim3('T')
 #endif
-         IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice_lim( kt )   ! -- bdy ice thermo 
+         IF( ln_limthd .AND. ln_bdy ) CALL bdy_ice( kt )       ! -- bdy ice thermo 
          ! previous lead fraction and ice volume for flux calculations
                                       CALL lim_var_glo2eqv     ! ht_i and ht_s for ice albedo calculation
@@ -278,5 +275,5 @@
       IF( ierr /= 0 )   CALL ctl_stop('STOP', 'ice_init : unable to allocate ice arrays')
       !
-      CALL lim_dyn_init                ! set ice dynamics parameters
+      CALL ice_rhg_init                ! set ice dynamics parameters
       !
       CALL ice_itd_init                ! ice thickness distribution initialization
Index: anches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90	(revision 8406)
+++ 	(revision )
@@ -1,169 +1,0 @@
-MODULE limdyn
-   !!======================================================================
-   !!                     ***  MODULE  limdyn  ***
-   !!   Sea-Ice dynamics :  
-   !!======================================================================
-   !! history :  1.0  ! 2002-08  (C. Ethe, G. Madec)  original VP code 
-   !!            3.0  ! 2007-03  (MA Morales Maqueda, S. Bouillon, M. Vancoppenolle)  LIM3: EVP-Cgrid
-   !!            3.5  ! 2011-02  (G. Madec) dynamical allocation
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3' :                                 LIM3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!    lim_dyn      : computes ice velocities
-   !!    lim_dyn_init : initialization and namelist read
-   !!----------------------------------------------------------------------
-   USE phycst           ! physical constants
-   USE dom_oce          ! ocean space and time domain
-   USE ice              ! LIM-3 variables
-   USE limrhg           ! LIM-3 rheology
-   USE limcons          ! conservation tests
-   USE limctl           ! control prints
-   USE limvar
-   !
-   USE lbclnk           ! lateral boundary conditions - MPP exchanges
-   USE lib_mpp          ! MPP library
-   USE in_out_manager   ! I/O manager
-   USE lib_fortran      ! glob_sum
-   USE timing           ! Timing
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   lim_dyn        ! routine called by icestp.F90
-   PUBLIC   lim_dyn_init   ! routine called by icestp.F90
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE lim_dyn( kt )
-      !!-------------------------------------------------------------------
-      !!               ***  ROUTINE lim_dyn  ***
-      !!               
-      !! ** Purpose :   compute ice velocity
-      !!                
-      !! ** Method  : 
-      !!
-      !! ** Action  : - Initialisation
-      !!              - Call of the dynamic routine for each hemisphere
-      !!------------------------------------------------------------------------------------
-      INTEGER, INTENT(in) ::   kt     ! number of iteration
-      !!
-      INTEGER  :: jl, jk ! dummy loop indices
-      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 
-     !!---------------------------------------------------------------------
-
-      IF( nn_timing == 1 )  CALL timing_start('limdyn')
-
-      CALL lim_var_agg(1)                      ! aggregate ice categories
-      !
-      ! conservation test
-      IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-      
-      ! ice velocities before rheology
-      u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1)
-      v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1)
-      
-      ! Landfast ice parameterization: define max bottom friction
-      tau_icebfr(:,:) = 0._wp
-      IF( ln_landfast ) THEN
-         DO jl = 1, jpl
-            WHERE( ht_i(:,:,jl) > ht_n(:,:) * rn_gamma )  tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr
-         END DO
-      ENDIF
-      
-      ! Rheology (ice dynamics)
-      ! ========     
-      CALL lim_rhg
-      !
-      ! conservation test
-      IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)
-
-      ! Control prints
-      IF( ln_ctl )       CALL lim_prt3D( 'limdyn' )
-      !
-      IF( nn_timing == 1 )  CALL timing_stop('limdyn')
-
-   END SUBROUTINE lim_dyn
-
-
-   SUBROUTINE lim_dyn_init
-      !!-------------------------------------------------------------------
-      !!                  ***  ROUTINE lim_dyn_init  ***
-      !!
-      !! ** Purpose : Physical constants and parameters linked to the ice
-      !!      dynamics
-      !!
-      !! ** Method  :  Read the namicedyn namelist and check the ice-dynamic
-      !!       parameter values called at the first timestep (nit000)
-      !!
-      !! ** input   :   Namelist namicedyn
-      !!-------------------------------------------------------------------
-      INTEGER  ::   ios                 ! Local integer output status for namelist read
-      NAMELIST/namicedyn/ nn_limadv, nn_limadv_ord,                                &
-         &                nn_icestr, rn_pe_rdg, rn_pstar, rn_crhg, ln_icestr_bvf,  &
-         &                rn_ishlat, rn_cio, rn_creepl, rn_ecc,                    &
-         &                nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr, rn_lfrelax
-      !!-------------------------------------------------------------------
-
-      REWIND( numnam_ice_ref )              ! Namelist namicedyn in reference namelist : Ice dynamics
-      READ  ( numnam_ice_ref, namicedyn, IOSTAT = ios, ERR = 901)
-901   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in reference namelist', lwp )
-
-      REWIND( numnam_ice_cfg )              ! Namelist namicedyn in configuration namelist : Ice dynamics
-      READ  ( numnam_ice_cfg, namicedyn, IOSTAT = ios, ERR = 902 )
-902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namicedyn in configuration namelist', lwp )
-      IF(lwm) WRITE ( numoni, namicedyn )
-      
-      IF(lwp) THEN                        ! control print
-         WRITE(numout,*)
-         WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics '
-         WRITE(numout,*) '~~~~~~~~~~~~'
-         ! limtrp
-         WRITE(numout,*)'    choose the advection scheme (-1=Prather, 0=Ulimate-Macho)   nn_limadv     = ', nn_limadv 
-         WRITE(numout,*)'    choose the order of the scheme (if ultimate)                nn_limadv_ord = ', nn_limadv_ord  
-         ! limitd_me
-         WRITE(numout,*)'    ice strength parameterization (0=Hibler 1=Rothrock)         nn_icestr     = ', nn_icestr 
-         WRITE(numout,*)'    Ratio of ridging work to PotEner change in ridging          rn_pe_rdg     = ', rn_pe_rdg 
-         WRITE(numout,*) '   first bulk-rheology parameter                               rn_pstar      = ', rn_pstar
-         WRITE(numout,*) '   second bulk-rhelogy parameter                               rn_crhg       = ', rn_crhg
-         WRITE(numout,*)'    Including brine volume in ice strength comp.                ln_icestr_bvf = ', ln_icestr_bvf
-         ! limrhg
-         WRITE(numout,*) '   lateral boundary condition for sea ice dynamics             rn_ishlat     = ', rn_ishlat
-         WRITE(numout,*) '   drag coefficient for oceanic stress                         rn_cio        = ', rn_cio
-         WRITE(numout,*) '   creep limit                                                 rn_creepl     = ', rn_creepl
-         WRITE(numout,*) '   eccentricity of the elliptical yield curve                  rn_ecc        = ', rn_ecc
-         WRITE(numout,*) '   number of iterations for subcycling                         nn_nevp       = ', nn_nevp
-         WRITE(numout,*) '   ratio of elastic timescale over ice time step               rn_relast     = ', rn_relast
-         WRITE(numout,*) '   Landfast: param (T or F)                                    ln_landfast   = ', ln_landfast
-         WRITE(numout,*) '      T: fraction of ocean depth that ice must reach           rn_gamma      = ', rn_gamma
-         WRITE(numout,*) '      T: maximum bottom stress per unit area of contact        rn_icebfr     = ', rn_icebfr
-         WRITE(numout,*) '      T: relax time scale (s-1) to reach static friction       rn_lfrelax    = ', rn_lfrelax
-      ENDIF
-      !
-      IF     (      rn_ishlat == 0.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  free-slip '
-      ELSEIF (      rn_ishlat == 2.                ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  no-slip '
-      ELSEIF ( 0. < rn_ishlat .AND. rn_ishlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  partial-slip '
-      ELSEIF ( 2. < rn_ishlat                      ) THEN   ;   IF(lwp) WRITE(numout,*) '   ice lateral  strong-slip '
-      ENDIF
-      !
-   END SUBROUTINE lim_dyn_init
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option          Empty module           NO LIM sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_dyn         ! Empty routine
-   END SUBROUTINE lim_dyn
-#endif 
-
-   !!======================================================================
-END MODULE limdyn
Index: anches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90	(revision 8406)
+++ 	(revision )
@@ -1,770 +1,0 @@
-MODULE limrhg
-   !!======================================================================
-   !!                     ***  MODULE  limrhg  ***
-   !!   Ice rheology : sea ice rheology
-   !!======================================================================
-   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
-   !!            3.0  !  2008-03  (M. Vancoppenolle) LIM3
-   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy 
-   !!            3.3  !  2009-05  (G.Garric) addition of the evp cas
-   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation 
-   !!            3.5  !  2012-08  (R. Benshila)  AGRIF
-   !!            3.6  !  2016-06  (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013)
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                      LIM-3 sea-ice model
-   !!----------------------------------------------------------------------
-   !!   lim_rhg       : computes ice velocities
-   !!----------------------------------------------------------------------
-   USE phycst         ! Physical constant
-   USE par_oce        ! Ocean parameters
-   USE dom_oce        ! Ocean domain
-   USE sbc_oce , ONLY : ln_ice_embd, nn_fsbc, ssh_m
-   USE sbc_ice , ONLY : utau_ice, vtau_ice, snwice_mass, snwice_mass_b
-   USE ice            ! ice variables
-   USE limitd_me      ! ice strength
-   !
-   USE lbclnk         ! Lateral Boundary Condition / MPP link
-   USE lib_mpp        ! MPP library
-   USE in_out_manager ! I/O manager
-   USE prtctl         ! Print control
-   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)  
-#if defined key_agrif
-   USE agrif_lim3_interp
-#endif
-   USE bdy_oce   , ONLY: ln_bdy 
-   USE bdyice_lim 
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   lim_rhg        ! routine called by lim_dyn
-
-   !! * Substitutions
-#  include "vectopt_loop_substitute.h90"
-   !!----------------------------------------------------------------------
-   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
-   !! $Id$
-   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE lim_rhg
-      !!-------------------------------------------------------------------
-      !!                 ***  SUBROUTINE lim_rhg  ***
-      !!                          EVP-C-grid
-      !!
-      !! ** purpose : determines sea ice drift from wind stress, ice-ocean
-      !!  stress and sea-surface slope. Ice-ice interaction is described by 
-      !!  a non-linear elasto-viscous-plastic (EVP) law including shear 
-      !!  strength and a bulk rheology (Hunke and Dukowicz, 2002).	
-      !!
-      !!  The points in the C-grid look like this, dear reader
-      !!
-      !!                              (ji,jj)
-      !!                                 |
-      !!                                 |
-      !!                      (ji-1,jj)  |  (ji,jj)
-      !!                             ---------   
-      !!                            |         |
-      !!                            | (ji,jj) |------(ji,jj)
-      !!                            |         |
-      !!                             ---------   
-      !!                     (ji-1,jj-1)     (ji,jj-1)
-      !!
-      !! ** Inputs  : - wind forcing (stress), oceanic currents
-      !!                ice total volume (vt_i) per unit area
-      !!                snow total volume (vt_s) per unit area
-      !!
-      !! ** Action  : - compute u_ice, v_ice : the components of the 
-      !!                sea-ice velocity vector
-      !!              - compute delta_i, shear_i, divu_i, which are inputs
-      !!                of the ice thickness distribution
-      !!
-      !! ** Steps   : 1) Compute ice snow mass, ice strength 
-      !!              2) Compute wind, oceanic stresses, mass terms and
-      !!                 coriolis terms of the momentum equation
-      !!              3) Solve the momentum equation (iterative procedure)
-      !!              4) Prevent high velocities if the ice is thin
-      !!              5) Recompute invariants of the strain rate tensor
-      !!                 which are inputs of the ITD, store stress
-      !!                 for the next time step
-      !!              6) Control prints of residual (convergence)
-      !!                 and charge ellipse.
-      !!                 The user should make sure that the parameters
-      !!                 nn_nevp, elastic time scale and rn_creepl maintain stress state
-      !!                 on the charge ellipse for plastic flow
-      !!                 e.g. in the Canadian Archipelago
-      !!
-      !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013
-      !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters)
-      !!              but this solution appears very unstable (see Kimmritz et al 2016)
-      !!
-      !! References : Hunke and Dukowicz, JPO97
-      !!              Bouillon et al., Ocean Modelling 2009
-      !!              Bouillon et al., Ocean Modelling 2013
-      !!-------------------------------------------------------------------
-      INTEGER ::   ji, jj       ! dummy loop indices
-      INTEGER ::   jter         ! local integers
-      CHARACTER (len=50) ::   charout
-
-      REAL(wp) ::   zrhoco                                                   ! rau0 * rn_cio
-      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling
-      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity
-      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013
-      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass
-      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars
-      REAL(wp) ::   zTauO, zTauB, zTauE, zvel                                ! temporary scalars
-
-      REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress
-      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity
-      REAL(wp) ::   zintb, zintn                                             ! dummy argument
-      REAL(wp) ::   zfac_x, zfac_y
-      
-      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors
-      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
-      !
-      REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points
-      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points
-      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfU   , zfV                     ! internal stresses
-      
-      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
-      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components
-      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence
-      REAL(wp), DIMENSION(jpi,jpj) ::   zpice                           ! array used for the calculation of ice surface slope:
-                                                                             !   ocean surface (ssh_m) if ice is not embedded
-                                                                             !   ice top surface if ice is embedded   
-      REAL(wp), DIMENSION(jpi,jpj) ::   zCorx, zCory                    ! Coriolis stress array
-      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! Ocean-to-ice stress array
-
-      REAL(wp), DIMENSION(jpi,jpj) ::   zswitchU, zswitchV              ! dummy arrays
-      REAL(wp), DIMENSION(jpi,jpj) ::   zmaskU, zmaskV                  ! mask for ice presence
-      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice
-
-      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter
-      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity
-      !!-------------------------------------------------------------------
-
-#if defined key_agrif 
-      CALL agrif_interp_lim3( 'U', 0, nn_nevp )   ! First interpolation of coarse values
-      CALL agrif_interp_lim3( 'V', 0, nn_nevp )
-#endif
-      !
-      !------------------------------------------------------------------------------!
-      ! 0) mask at F points for the ice
-      !------------------------------------------------------------------------------!
-      ! ocean/land mask
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1      ! NO vector opt.
-            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1)
-         END DO
-      END DO
-      CALL lbc_lnk( zfmask, 'F', 1._wp )
-
-      ! Lateral boundary conditions on velocity (modify zfmask)
-      zwf(:,:) = zfmask(:,:)
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1   ! vector opt.
-            IF( zfmask(ji,jj) == 0._wp ) THEN
-               zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) )
-            ENDIF
-         END DO
-      END DO
-      DO jj = 2, jpjm1
-         IF( zfmask(1,jj) == 0._wp ) THEN
-            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) )
-         ENDIF
-         IF( zfmask(jpi,jj) == 0._wp ) THEN
-            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) )
-         ENDIF
-      END DO
-      DO ji = 2, jpim1
-         IF( zfmask(ji,1) == 0._wp ) THEN
-            zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) )
-         ENDIF
-         IF( zfmask(ji,jpj) == 0._wp ) THEN
-            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) )
-         ENDIF
-      END DO
-      CALL lbc_lnk( zfmask, 'F', 1._wp )
-
-      !------------------------------------------------------------------------------!
-      ! 1) define some variables and initialize arrays
-      !------------------------------------------------------------------------------!
-      zrhoco = rau0 * rn_cio 
-
-      ! ecc2: square of yield ellipse eccenticrity
-      ecc2    = rn_ecc * rn_ecc
-      z1_ecc2 = 1._wp / ecc2
-
-      ! Time step for subcycling
-      zdtevp   = rdt_ice / REAL( nn_nevp )
-      z1_dtevp = 1._wp / zdtevp
-
-      ! alpha parameters (Bouillon 2009)
-      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
-      zalph2 = zalph1 * z1_ecc2
-
-      ! alpha and beta parameters (Bouillon 2013)
-      !!zalph1 = 40.
-      !!zalph2 = 40.
-      !!zbeta  = 3000.
-      !!zbeta = REAL( nn_nevp )   ! close to classical EVP of Hunke (2001)
-
-      z1_alph1 = 1._wp / ( zalph1 + 1._wp )
-      z1_alph2 = 1._wp / ( zalph2 + 1._wp )
-
-      ! Initialise stress tensor 
-      zs1 (:,:) = stress1_i (:,:) 
-      zs2 (:,:) = stress2_i (:,:)
-      zs12(:,:) = stress12_i(:,:)
-
-      ! Ice strength
-      CALL lim_itd_me_icestrength( nn_icestr )
-
-      ! scale factors
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) )
-            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) )
-         END DO
-      END DO
-            
-      !
-      !------------------------------------------------------------------------------!
-      ! 2) Wind / ocean stress, mass terms, coriolis terms
-      !------------------------------------------------------------------------------!
-
-      IF( ln_ice_embd ) THEN             !== embedded sea ice: compute representative ice top surface ==!
-         !                                            
-         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[n/nn_fsbc], n=0,nn_fsbc-1}
-         !                                               = (1/nn_fsbc)^2 * {SUM[n], n=0,nn_fsbc-1}
-         zintn = REAL( nn_fsbc - 1 ) / REAL( nn_fsbc ) * 0.5_wp     
-         !
-         ! average interpolation coeff as used in dynspg = (1/nn_fsbc) * {SUM[1-n/nn_fsbc], n=0,nn_fsbc-1}
-         !                                               = (1/nn_fsbc)^2 * (nn_fsbc^2 - {SUM[n], n=0,nn_fsbc-1})
-         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp
-         !
-         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0
-         !
-      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==!
-         zpice(:,:) = ssh_m(:,:)
-      ENDIF
-
-      DO jj = 2, jpjm1
-         DO ji = fs_2, fs_jpim1
-
-            ! ice fraction at U-V points
-            zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
-            zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
-
-            ! Ice/snow mass at U-V points
-            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
-            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
-            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
-            zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
-            zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
-
-            ! Ocean currents at U-V points
-            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    &
-               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
-            
-            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    &
-               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
-
-            ! Coriolis at T points (m*f)
-            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
-
-            ! m/dt
-            zmU_t(ji,jj)    = zmassU * z1_dtevp
-            zmV_t(ji,jj)    = zmassV * z1_dtevp
-
-            ! Drag ice-atm.
-            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
-            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj)
-
-            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points
-            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
-            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
-
-            ! masks
-            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
-            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
-
-            ! switches
-            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin
-            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin
-
-         END DO
-      END DO
-      CALL lbc_lnk( zmf, 'T', 1. )
-      !
-      !------------------------------------------------------------------------------!
-      ! 3) Solution of the momentum equation, iterative procedure
-      !------------------------------------------------------------------------------!
-      !
-      !                                               !----------------------!
-      DO jter = 1 , nn_nevp                           !    loop over jter    !
-         !                                            !----------------------!        
-         IF(ln_ctl) THEN   ! Convergence test
-            DO jj = 1, jpjm1
-               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step
-               zv_ice(:,jj) = v_ice(:,jj)
-            END DO
-         ENDIF
-
-         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
-         DO jj = 1, jpjm1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points
-            DO ji = 1, jpim1
-
-               ! shear at F points
-               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
-                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
-
-            END DO
-         END DO
-         CALL lbc_lnk( zds, 'F', 1. )
-
-         DO jj = 2, jpjm1
-            DO ji = 2, jpim1 ! no vector loop
-
-               ! shear**2 at T points (doc eq. A16)
-               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
-                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
-                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
-              
-               ! divergence at T points
-               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
-                  &    ) * r1_e1e2t(ji,jj)
-               zdiv2 = zdiv * zdiv
-               
-               ! tension at T points
-               zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
-                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
-                  &   ) * r1_e1e2t(ji,jj)
-               zdt2 = zdt * zdt
-               
-               ! delta at T points
-               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
-
-               ! P/delta at T points
-               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
-               
-               ! stress at T points
-               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1
-               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2
-             
-            END DO
-         END DO
-         CALL lbc_lnk( zp_delt, 'T', 1. )
-
-         DO jj = 1, jpjm1
-            DO ji = 1, jpim1
-
-               ! P/delta at F points
-               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
-               
-               ! stress at F points
-               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2
-
-            END DO
-         END DO
-         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
- 
-
-         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
-         DO jj = 2, jpjm1
-            DO ji = fs_2, fs_jpim1               
-
-               ! U points
-               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
-                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
-                  &                    ) * r1_e2u(ji,jj)                                                                      &
-                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
-                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
-                  &                  ) * r1_e1e2u(ji,jj)
-
-               ! V points
-               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
-                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
-                  &                    ) * r1_e1v(ji,jj)                                                                      &
-                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
-                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
-                  &                  ) * r1_e1e2v(ji,jj)
-
-               ! u_ice at V point
-               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     &
-                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1)
-               
-               ! v_ice at U point
-               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     &
-                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1)
-
-            END DO
-         END DO
-         !
-         ! --- Computation of ice velocity --- !
-         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp
-         !  Bouillon et al. 2009 (eq 34-35) => stable
-         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations
-            
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-
-                  ! tau_io/(v_oce - v_ice)
-                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
-                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
-
-                  ! Ocean-to-Ice stress
-                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
-
-                  ! tau_bottom/v_ice
-                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
-                  zTauB = - tau_icebfr(ji,jj) / zvel
-
-                  ! Coriolis at V-points (energy conserving formulation)
-                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
-                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
-                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
-
-                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-                  
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
-                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
-                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
-                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
-                     &           ) * zmaskV(ji,jj)
-                  ! Bouillon 2013
-                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
-                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
-                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
-
-               END DO
-            END DO
-            CALL lbc_lnk( v_ice, 'V', -1. )
-            
-#if defined key_agrif
-!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
-            CALL agrif_interp_lim3( 'V' )
-#endif
-            IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'V' )
-
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-                               
-                  ! tau_io/(u_oce - u_ice)
-                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
-                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
-
-                  ! Ocean-to-Ice stress
-                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
-
-                  ! tau_bottom/u_ice
-                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
-                  zTauB = - tau_icebfr(ji,jj) / zvel
-
-                  ! Coriolis at U-points (energy conserving formulation)
-                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
-                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
-                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
-                  
-                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
-                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
-                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
-                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin 
-                     &           ) * zmaskU(ji,jj)
-                  ! Bouillon 2013
-                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
-                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
-                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
-               END DO
-            END DO
-            CALL lbc_lnk( u_ice, 'U', -1. )
-            
-#if defined key_agrif
-!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
-            CALL agrif_interp_lim3( 'U' )
-#endif
-            IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'U' )
-
-         ELSE ! odd iterations
-
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-
-                  ! tau_io/(u_oce - u_ice)
-                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
-                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
-
-                  ! Ocean-to-Ice stress
-                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
-
-                  ! tau_bottom/u_ice
-                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) )
-                  zTauB = - tau_icebfr(ji,jj) / zvel
-
-                  ! Coriolis at U-points (energy conserving formulation)
-                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  &
-                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  &
-                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) )
-
-                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction ; 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
-                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                   &  ! previous velocity
-                     &                                     + zTauE + zTauO * u_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
-                     &             + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &             ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
-                     &           ) * zmaskU(ji,jj)
-                  ! Bouillon 2013
-                  !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  &
-                  !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  &
-                  !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj)
-               END DO
-            END DO
-            CALL lbc_lnk( u_ice, 'U', -1. )
-            
-#if defined key_agrif
-!!            CALL agrif_interp_lim3( 'U', jter, nn_nevp )
-            CALL agrif_interp_lim3( 'U' )
-#endif
-            IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'U' )
-
-            DO jj = 2, jpjm1
-               DO ji = fs_2, fs_jpim1
-         
-                  ! tau_io/(v_oce - v_ice)
-                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  &
-                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) )
-
-                  ! Ocean-to-Ice stress
-                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) )
-
-                  ! tau_bottom/v_ice
-                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) )
-                  ztauB = - tau_icebfr(ji,jj) / zvel
-                  
-                  ! Coriolis at V-points (energy conserving formulation)
-                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  &
-                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  &
-                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) )
-
-                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io
-                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
-
-                  ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction
-                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) )
-                  
-                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009)
-                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                   &  ! previous velocity
-                     &                                     + zTauE + zTauO * v_ice(ji,jj)                  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
-                     &                                     ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )  &  ! m/dt + tau_io(only ice part) + landfast
-                     &             + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0
-                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )        &  ! v_ice = v_oce if mass < zmmin
-                     &           ) * zmaskV(ji,jj)
-                  ! Bouillon 2013
-                  !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  &
-                  !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  &
-                  !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj)
-               END DO
-            END DO
-            CALL lbc_lnk( v_ice, 'V', -1. )
-            
-#if defined key_agrif
-!!            CALL agrif_interp_lim3( 'V', jter, nn_nevp )
-            CALL agrif_interp_lim3( 'V' )
-#endif
-            IF( ln_bdy ) CALL bdy_ice_lim_dyn( 'V' )
-
-         ENDIF
-         
-         IF(ln_ctl) THEN   ! Convergence test
-            DO jj = 2 , jpjm1
-               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
-            END DO
-            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) )
-            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
-         ENDIF
-         !
-         !                                                ! ==================== !
-      END DO                                              !  end loop over jter  !
-      !                                                   ! ==================== !
-      !
-      !------------------------------------------------------------------------------!
-      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 
-      !------------------------------------------------------------------------------!
-      DO jj = 1, jpjm1
-         DO ji = 1, jpim1
-
-            ! shear at F points
-            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
-               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
-               &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
-
-         END DO
-      END DO           
-      CALL lbc_lnk( zds, 'F', 1. )
-      
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1 ! no vector loop
-            
-            ! tension**2 at T points
-            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
-               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
-               &   ) * r1_e1e2t(ji,jj)
-            zdt2 = zdt * zdt
-            
-            ! shear**2 at T points (doc eq. A16)
-            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
-               &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
-               &   ) * 0.25_wp * r1_e1e2t(ji,jj)
-            
-            ! shear at T points
-            shear_i(ji,jj) = SQRT( zdt2 + zds2 )
-
-            ! divergence at T points
-            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
-               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
-               &            ) * r1_e1e2t(ji,jj)
-            
-            ! delta at T points
-            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
-            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
-            delta_i(ji,jj) = zdelta + rn_creepl * rswitch
-
-         END DO
-      END DO
-      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. )
-      
-      ! --- Store the stress tensor for the next time step --- !
-      stress1_i (:,:) = zs1 (:,:)
-      stress2_i (:,:) = zs2 (:,:)
-      stress12_i(:,:) = zs12(:,:)
-      !
-
-      !------------------------------------------------------------------------------!
-      ! 5) SIMIP diagnostics
-      !------------------------------------------------------------------------------!
-                           
-      DO jj = 2, jpjm1
-         DO ji = 2, jpim1
-             rswitch  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
-
-             ! Stress tensor invariants (normal and shear stress N/m)
-             diag_sig1(ji,jj) = ( zs1(ji,jj) + zs2(ji,jj) ) * rswitch                                 ! normal stress
-             diag_sig2(ji,jj) = SQRT( ( zs1(ji,jj) - zs2(ji,jj) )**2 + 4*zs12(ji,jj)**2 ) * rswitch   ! shear stress
-
-             ! Stress terms of the momentum equation (N/m2)
-             diag_dssh_dx(ji,jj) = zspgU(ji,jj) * rswitch    ! sea surface slope stress term
-             diag_dssh_dy(ji,jj) = zspgV(ji,jj) * rswitch
-
-             diag_corstrx(ji,jj) = zCorx(ji,jj) * rswitch    ! Coriolis stress term
-             diag_corstry(ji,jj) = zCory(ji,jj) * rswitch
-
-             diag_intstrx(ji,jj) = zfU(ji,jj)   * rswitch    ! internal stress term
-             diag_intstry(ji,jj) = zfV(ji,jj)   * rswitch
-           
-             diag_utau_oi(ji,jj) = ztaux_oi(ji,jj) * rswitch  ! oceanic stress
-             diag_vtau_oi(ji,jj) = ztauy_oi(ji,jj) * rswitch
-
-             ! 2D ice mass, snow mass, area transport arrays (X, Y)
-             zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * rswitch
-             zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * rswitch
-
-             diag_xmtrp_ice(ji,jj) = rhoic * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
-             diag_ymtrp_ice(ji,jj) = rhoic * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
-
-             diag_xmtrp_snw(ji,jj) = rhosn * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
-             diag_ymtrp_snw(ji,jj) = rhosn * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
-
-             diag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )         ! area transport,      X-component
-             diag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )         !        ''            Y-   ''
-
-         END DO
-      END DO
-
-      CALL lbc_lnk_multi(   diag_sig1   , 'T',  1., diag_sig2   , 'T',  1.,   &
-                 &          diag_dssh_dx, 'U', -1., diag_dssh_dy, 'V', -1.,   &
-                 &          diag_corstrx, 'U', -1., diag_corstry, 'V', -1.,   & 
-                 &          diag_intstrx, 'U', -1., diag_intstry, 'V', -1.    )
-
-      CALL lbc_lnk_multi(   diag_utau_oi, 'U', -1., diag_vtau_oi, 'V', -1.    )
-
-      CALL lbc_lnk_multi(   diag_xmtrp_ice, 'U', -1., diag_xmtrp_snw, 'U', -1., &
-                 &          diag_xatrp    , 'U', -1., diag_ymtrp_ice, 'V', -1., &
-                 &          diag_ymtrp_snw, 'V', -1., diag_yatrp    , 'V', -1.  )
-
-      !
-      !------------------------------------------------------------------------------!
-      ! 6) Control prints of residual and charge ellipse
-      !------------------------------------------------------------------------------!
-      !
-      ! print the residual for convergence
-      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 :')
-      ENDIF
-
-      ! print charge ellipse
-      ! This can be desactivated once the user is sure that the stress state
-      ! lie on the charge ellipse. See Bouillon et al. 08 for more details
-      IF(ln_ctl) THEN
-         IF( MOD(kt_ice+nn_fsbc-1,nwrite) == 0 ) THEN
-            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, kt_ice, 0, 0, ' ch. ell. '
-            CALL prt_ctl_info(charout)
-            DO jj = 2, jpjm1
-               DO ji = 2, jpim1
-                  IF (strength(ji,jj) > 1.0) THEN
-                     zsig1 = ( zs1(ji,jj) + SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) ) 
-                     zsig2 = ( zs1(ji,jj) - SQRT(zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 ) ) / ( 2*strength(ji,jj) )
-                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)")
-                     CALL prt_ctl_info(charout)
-                  ENDIF
-               END DO
-            END DO
-            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 2000, kt_ice, 0, 0, ' ch. ell. '
-            CALL prt_ctl_info(charout)
-         ENDIF
-      ENDIF
-      !
-   END SUBROUTINE lim_rhg
-
-#else
-   !!----------------------------------------------------------------------
-   !!   Default option          Dummy module           NO LIM sea-ice model
-   !!----------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE lim_rhg         ! Dummy routine
-      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?'
-   END SUBROUTINE lim_rhg
-#endif
-
-   !!==============================================================================
-END MODULE limrhg
Index: /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90	(revision 8407)
+++ /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90	(revision 8407)
@@ -0,0 +1,366 @@
+MODULE bdyice
+   !!======================================================================
+   !!                       ***  MODULE  bdyice  ***
+   !! Unstructured Open Boundary Cond. :  Open boundary conditions for sea-ice (LIM3)
+   !!======================================================================
+   !!  History :  3.3  !  2010-09 (D. Storkey)  Original code
+   !!             3.4  !  2011    (D. Storkey)  rewrite in preparation for OBC-BDY merge
+   !!              -   !  2012-01 (C. Rousset)  add lim3 and remove useless jk loop 
+   !!----------------------------------------------------------------------
+#if defined key_lim3
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'                                                 LIM-3 sea ice model
+   !!----------------------------------------------------------------------
+   !!   bdy_ice        : Application of open boundaries to ice
+   !!   bdy_ice_frs        : Application of Flow Relaxation Scheme
+   !!----------------------------------------------------------------------
+   USE timing          ! Timing
+   USE phycst          ! physical constant
+   USE eosbn2          ! equation of state
+   USE oce             ! ocean dynamics and tracers variables
+   USE ice             ! LIM_3 ice variables
+   USE limvar
+   USE limctl
+   USE par_oce         ! ocean parameters
+   USE dom_oce         ! ocean space and time domain variables 
+   USE sbc_oce         ! Surface boundary condition: ocean fields
+   USE bdy_oce         ! ocean open boundary conditions
+   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
+   USE in_out_manager  ! write to numout file
+   USE lib_mpp         ! distributed memory computing
+   USE lib_fortran     ! to use key_nosignedzero
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   bdy_ice     ! routine called in sbcmod
+   PUBLIC   bdy_ice_dyn ! routine called in limrhg
+
+   !!----------------------------------------------------------------------
+   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
+   !! $Id: bdyice.F90 8306 2017-07-10 10:18:03Z clem $
+   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE bdy_ice( kt )
+      !!----------------------------------------------------------------------
+      !!                  ***  SUBROUTINE bdy_ice  ***
+      !!
+      !! ** Purpose : - Apply open boundary conditions for ice (LIM3)
+      !!
+      !!----------------------------------------------------------------------
+      INTEGER, INTENT( in ) ::   kt   ! Main time step counter
+      !
+      INTEGER ::   ib_bdy   ! Loop index
+      !!----------------------------------------------------------------------
+      !
+      CALL lim_var_glo2eqv
+      !
+      DO ib_bdy=1, nb_bdy
+         !
+         SELECT CASE( cn_ice_lim(ib_bdy) )
+         CASE('none')
+            CYCLE
+         CASE('frs')
+            CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
+         CASE DEFAULT
+            CALL ctl_stop( 'bdy_ice : unrecognised option for open boundaries for ice fields' )
+         END SELECT
+         !
+      END DO
+      !
+                        CALL lim_var_zapsmall
+                        CALL lim_var_agg(1)
+      IF( ln_limctl )   CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )
+      !
+   END SUBROUTINE bdy_ice
+
+
+   SUBROUTINE bdy_ice_frs( idx, dta, kt, ib_bdy )
+      !!------------------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_ice_frs  ***
+      !!                    
+      !! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields in the case 
+      !!              of unstructured open boundaries.
+      !! 
+      !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in a three-
+      !!             dimensional baroclinic ocean model with realistic topography. Tellus, 365-382.
+      !!------------------------------------------------------------------------------
+      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices
+      TYPE(OBC_DATA),  INTENT(in) ::   dta     ! OBC external data
+      INTEGER,         INTENT(in) ::   kt      ! main time-step counter
+      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
+      !
+      INTEGER  ::   jpbound            ! 0 = incoming ice
+      !                                ! 1 = outgoing ice
+      INTEGER  ::   jb, jk, jgrd, jl   ! dummy loop indices
+      INTEGER  ::   ji, jj, ii, ij     ! local scalar
+      REAL(wp) ::   zwgt, zwgt1        ! local scalar
+      REAL(wp) ::   ztmelts, zdh
+      !!------------------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 )   CALL timing_start('bdy_ice_frs')
+      !
+      jgrd = 1      ! Everything is at T-points here
+      !
+      DO jl = 1, jpl
+         DO jb = 1, idx%nblenrim(jgrd)
+            ji    = idx%nbi(jb,jgrd)
+            jj    = idx%nbj(jb,jgrd)
+            zwgt  = idx%nbw(jb,jgrd)
+            zwgt1 = 1.e0 - idx%nbw(jb,jgrd)
+            a_i (ji,jj,jl) = ( a_i (ji,jj,jl) * zwgt1 + dta%a_i (jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Leads fraction 
+            ht_i(ji,jj,jl) = ( ht_i(ji,jj,jl) * zwgt1 + dta%ht_i(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice depth 
+            ht_s(ji,jj,jl) = ( ht_s(ji,jj,jl) * zwgt1 + dta%ht_s(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow depth
+
+            ! -----------------
+            ! Pathological case
+            ! -----------------
+            ! In case a) snow load would be in excess or b) ice is coming into a warmer environment that would lead to 
+            ! very large transformation from snow to ice (see limthd_dh.F90)
+
+            ! Then, a) transfer the snow excess into the ice (different from limthd_dh)
+            zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 )
+            ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment)
+            !zdh = MAX( 0._wp, ht_s(ji,jj,jl) * rhosn / rhoic )
+
+            ! recompute ht_i, ht_s
+            ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
+            ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn ) 
+
+         ENDDO
+         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy )
+      ENDDO
+      ! retrieve at_i
+      at_i(:,:) = 0._wp
+      DO jl = 1, jpl
+         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
+      END DO
+
+      DO jl = 1, jpl
+         DO jb = 1, idx%nblenrim(jgrd)
+            ji    = idx%nbi(jb,jgrd)
+            jj    = idx%nbj(jb,jgrd)
+
+            ! condition on ice thickness depends on the ice velocity
+            ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values
+            jpbound = 0   ;   ii = ji   ;   ij = jj
+            !
+            IF( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. ) jpbound = 1; ii = ji+1; ij = jj
+            IF( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. ) jpbound = 1; ii = ji-1; ij = jj
+            IF( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj+1
+            IF( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj-1
+            !
+            IF( nn_ice_lim_dta(ib_bdy) == 0 ) jpbound = 0; ii = ji; ij = jj   ! case ice boundaries = initial conditions
+            !                                                                 !      do not make state variables dependent on velocity
+            !
+            rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ii,ij) - 0.01 ) ) ! 0 if no ice
+            !
+            ! concentration and thickness
+            a_i (ji,jj,jl) = a_i (ii,ij,jl) * rswitch
+            ht_i(ji,jj,jl) = ht_i(ii,ij,jl) * rswitch
+            ht_s(ji,jj,jl) = ht_s(ii,ij,jl) * rswitch
+            !
+            ! Ice and snow volumes
+            v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
+            v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl)
+            !
+            SELECT CASE( jpbound )
+            !
+            CASE( 0 )   ! velocity is inward
+               !
+               ! Ice salinity, age, temperature
+               sm_i(ji,jj,jl)   = rswitch * rn_ice_sal(ib_bdy)  + ( 1.0 - rswitch ) * rn_simin
+               oa_i(ji,jj,jl)   = rswitch * rn_ice_age(ib_bdy) * a_i(ji,jj,jl)
+               t_su(ji,jj,jl)   = rswitch * rn_ice_tem(ib_bdy)  + ( 1.0 - rswitch ) * rn_ice_tem(ib_bdy)
+               DO jk = 1, nlay_s
+                  t_s(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rt0
+               END DO
+               DO jk = 1, nlay_i
+                  t_i(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rt0 
+                  s_i(ji,jj,jk,jl) = rswitch * rn_ice_sal(ib_bdy) + ( 1.0 - rswitch ) * rn_simin
+               END DO
+               !
+            CASE( 1 )   ! velocity is outward
+               !
+               ! Ice salinity, age, temperature
+               sm_i(ji,jj,jl)   = rswitch * sm_i(ii,ij,jl)  + ( 1.0 - rswitch ) * rn_simin
+               oa_i(ji,jj,jl)   = rswitch * oa_i(ii,ij,jl)
+               t_su(ji,jj,jl)   = rswitch * t_su(ii,ij,jl)  + ( 1.0 - rswitch ) * rt0
+               DO jk = 1, nlay_s
+                  t_s(ji,jj,jk,jl) = rswitch * t_s(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rt0
+               END DO
+               DO jk = 1, nlay_i
+                  t_i(ji,jj,jk,jl) = rswitch * t_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rt0
+                  s_i(ji,jj,jk,jl) = rswitch * s_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rn_simin
+               END DO
+               !
+            END SELECT
+            !
+            IF( nn_icesal == 1 ) THEN     ! constant salinity : overwrite rn_icesal
+               sm_i(ji,jj  ,jl) = rn_icesal
+               s_i (ji,jj,:,jl) = rn_icesal
+            ENDIF
+            !
+            ! contents
+            smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
+            DO jk = 1, nlay_s
+               ! Snow energy of melting
+               e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
+               ! Multiply by volume, so that heat content in J/m2
+               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
+            END DO
+            DO jk = 1, nlay_i
+               ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K                  
+               ! heat content per unit volume
+               e_i(ji,jj,jk,jl) = rswitch * rhoic * &
+                  (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
+                  +   lfus    * ( 1.0 - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
+                  - rcp      * ( ztmelts - rt0 ) )
+               ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
+               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) * r1_nlay_i
+            END DO
+            !
+         END DO
+         !
+         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk(  v_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk(  v_s(:,:,jl), 'T', 1., ib_bdy )
+         !
+         CALL lbc_bdy_lnk( smv_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk(  sm_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk(  oa_i(:,:,jl), 'T', 1., ib_bdy )
+         CALL lbc_bdy_lnk(  t_su(:,:,jl), 'T', 1., ib_bdy )
+         DO jk = 1, nlay_s
+            CALL lbc_bdy_lnk(t_s(:,:,jk,jl), 'T', 1., ib_bdy )
+            CALL lbc_bdy_lnk(e_s(:,:,jk,jl), 'T', 1., ib_bdy )
+         END DO
+         DO jk = 1, nlay_i
+            CALL lbc_bdy_lnk(t_i(:,:,jk,jl), 'T', 1., ib_bdy )
+            CALL lbc_bdy_lnk(e_i(:,:,jk,jl), 'T', 1., ib_bdy )
+         END DO
+         !
+      END DO !jl
+      !
+      !      
+      IF( nn_timing == 1 )   CALL timing_stop('bdy_ice_frs')
+      !
+   END SUBROUTINE bdy_ice_frs
+
+
+   SUBROUTINE bdy_ice_dyn( cd_type )
+      !!------------------------------------------------------------------------------
+      !!                 ***  SUBROUTINE bdy_ice_dyn  ***
+      !!                    
+      !! ** Purpose : Apply dynamics boundary conditions for sea-ice in the cas of unstructured open boundaries.
+      !!              u_ice and v_ice are equal to the value of the adjacent grid point if this latter is not ice free
+      !!              if adjacent grid point is ice free, then u_ice and v_ice are equal to ocean velocities
+      !!
+      !! 2013-06 : C. Rousset
+      !!------------------------------------------------------------------------------
+      CHARACTER(len=1), INTENT(in)  ::   cd_type   ! nature of velocity grid-points
+      !
+      INTEGER  ::   jb, jgrd           ! dummy loop indices
+      INTEGER  ::   ji, jj             ! local scalar
+      INTEGER  ::   ib_bdy             ! Loop index
+      REAL(wp) ::   zmsk1, zmsk2, zflag
+      !!------------------------------------------------------------------------------
+      !
+      IF( nn_timing == 1 ) CALL timing_start('bdy_ice_dyn')
+      !
+      DO ib_bdy=1, nb_bdy
+         !
+         SELECT CASE( cn_ice_lim(ib_bdy) )
+         !
+         CASE('none')
+            CYCLE
+            !
+         CASE('frs')
+            !
+            IF( nn_ice_lim_dta(ib_bdy) == 0 ) CYCLE            ! case ice boundaries = initial conditions 
+            !                                                  !      do not change ice velocity (it is only computed by rheology)
+            SELECT CASE ( cd_type )
+            !     
+            CASE ( 'U' )  
+               jgrd = 2      ! u velocity
+               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(jgrd)
+                  ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
+                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
+                  zflag = idx_bdy(ib_bdy)%flagu(jb,jgrd)
+                  !
+                  IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries
+                     ! one of the two zmsk is always 0 (because of zflag)
+                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice
+                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice
+                     !  
+                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)
+                     u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + &
+                        &            u_ice(ji-1,jj) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + &
+                        &            u_oce(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
+                  ELSE                             ! everywhere else
+                     !u_ice(ji,jj) = u_oce(ji,jj)
+                     u_ice(ji,jj) = 0._wp
+                  ENDIF
+                  ! mask ice velocities
+                  rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ji,jj) - 0.01_wp ) ) ! 0 if no ice
+                  u_ice(ji,jj) = rswitch * u_ice(ji,jj)
+                  !
+               END DO
+               CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy )
+               !
+            CASE ( 'V' )
+               jgrd = 3      ! v velocity
+               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(jgrd)
+                  ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
+                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
+                  zflag = idx_bdy(ib_bdy)%flagv(jb,jgrd)
+                  !
+                  IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern boundaries
+                     ! one of the two zmsk is always 0 (because of zflag)
+                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice
+                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice
+                     !  
+                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)
+                     v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + &
+                        &            v_ice(ji,jj-1) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + &
+                        &            v_oce(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
+                  ELSE                             ! everywhere else
+                     !v_ice(ji,jj) = v_oce(ji,jj)
+                     v_ice(ji,jj) = 0._wp
+                  ENDIF
+                  ! mask ice velocities
+                  rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ji,jj) - 0.01 ) ) ! 0 if no ice
+                  v_ice(ji,jj) = rswitch * v_ice(ji,jj)
+                  !
+               END DO
+               CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy )
+               !
+            END SELECT
+            !
+         CASE DEFAULT
+            CALL ctl_stop( 'bdy_ice_dyn : unrecognised option for open boundaries for ice fields' )
+         END SELECT
+         !
+      END DO
+      !
+      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_dyn')
+      !
+    END SUBROUTINE bdy_ice_dyn
+
+#else
+   !!---------------------------------------------------------------------------------
+   !!   Default option                                                    Empty module
+   !!---------------------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE bdy_ice( kt )      ! Empty routine
+      WRITE(*,*) 'bdy_ice: You should not have seen this print! error?', kt
+   END SUBROUTINE bdy_ice
+#endif
+
+   !!=================================================================================
+END MODULE bdyice
Index: anches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90
===================================================================
--- /branches/2017/dev_r8183_ICEMODEL/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice_lim.F90	(revision 8406)
+++ 	(revision )
@@ -1,375 +1,0 @@
-MODULE bdyice_lim
-   !!======================================================================
-   !!                       ***  MODULE  bdyice_lim  ***
-   !! Unstructured Open Boundary Cond. :  Open boundary conditions for sea-ice (LIM3)
-   !!======================================================================
-   !!  History :  3.3  !  2010-09 (D. Storkey)  Original code
-   !!             3.4  !  2011    (D. Storkey)  rewrite in preparation for OBC-BDY merge
-   !!              -   !  2012-01 (C. Rousset)  add lim3 and remove useless jk loop 
-   !!----------------------------------------------------------------------
-#if defined key_lim3
-   !!----------------------------------------------------------------------
-   !!   'key_lim3'                                                 LIM-3 sea ice model
-   !!----------------------------------------------------------------------
-   !!   bdy_ice_lim        : Application of open boundaries to ice
-   !!   bdy_ice_frs        : Application of Flow Relaxation Scheme
-   !!----------------------------------------------------------------------
-   USE timing          ! Timing
-   USE phycst          ! physical constant
-   USE eosbn2          ! equation of state
-   USE oce             ! ocean dynamics and tracers variables
-#if defined key_lim3
-   USE ice             ! LIM_3 ice variables
-   USE limvar
-   USE limctl
-#endif 
-   USE par_oce         ! ocean parameters
-   USE dom_oce         ! ocean space and time domain variables 
-   USE sbc_oce         ! Surface boundary condition: ocean fields
-   USE bdy_oce         ! ocean open boundary conditions
-   USE lbclnk          ! ocean lateral boundary conditions (or mpp link)
-   USE in_out_manager  ! write to numout file
-   USE lib_mpp         ! distributed memory computing
-   USE lib_fortran     ! to use key_nosignedzero
-
-   IMPLICIT NONE
-   PRIVATE
-
-   PUBLIC   bdy_ice_lim     ! routine called in sbcmod
-   PUBLIC   bdy_ice_lim_dyn ! routine called in limrhg
-
-   !!----------------------------------------------------------------------
-   !! NEMO/OPA 3.3 , NEMO Consortium (2010)
-   !! $Id$
-   !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
-   !!----------------------------------------------------------------------
-CONTAINS
-
-   SUBROUTINE bdy_ice_lim( kt )
-      !!----------------------------------------------------------------------
-      !!                  ***  SUBROUTINE bdy_ice_lim  ***
-      !!
-      !! ** Purpose : - Apply open boundary conditions for ice (LIM3)
-      !!
-      !!----------------------------------------------------------------------
-      INTEGER, INTENT( in ) ::   kt   ! Main time step counter
-      !
-      INTEGER ::   ib_bdy   ! Loop index
-      !!----------------------------------------------------------------------
-      !
-#if defined key_lim3
-      CALL lim_var_glo2eqv
-#endif
-      !
-      DO ib_bdy=1, nb_bdy
-         !
-         SELECT CASE( cn_ice_lim(ib_bdy) )
-         CASE('none')
-            CYCLE
-         CASE('frs')
-            CALL bdy_ice_frs( idx_bdy(ib_bdy), dta_bdy(ib_bdy), kt, ib_bdy )
-         CASE DEFAULT
-            CALL ctl_stop( 'bdy_ice_lim : unrecognised option for open boundaries for ice fields' )
-         END SELECT
-         !
-      END DO
-      !
-#if defined key_lim3
-                        CALL lim_var_zapsmall
-                        CALL lim_var_agg(1)
-      IF( ln_limctl )   CALL lim_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' )
-#endif
-      !
-   END SUBROUTINE bdy_ice_lim
-
-
-   SUBROUTINE bdy_ice_frs( idx, dta, kt, ib_bdy )
-      !!------------------------------------------------------------------------------
-      !!                 ***  SUBROUTINE bdy_ice_frs  ***
-      !!                    
-      !! ** Purpose : Apply the Flow Relaxation Scheme for sea-ice fields in the case 
-      !!              of unstructured open boundaries.
-      !! 
-      !! Reference : Engedahl H., 1995: Use of the flow relaxation scheme in a three-
-      !!             dimensional baroclinic ocean model with realistic topography. Tellus, 365-382.
-      !!------------------------------------------------------------------------------
-      TYPE(OBC_INDEX), INTENT(in) ::   idx     ! OBC indices
-      TYPE(OBC_DATA),  INTENT(in) ::   dta     ! OBC external data
-      INTEGER,         INTENT(in) ::   kt      ! main time-step counter
-      INTEGER,         INTENT(in) ::   ib_bdy  ! BDY set index
-      !
-      INTEGER  ::   jpbound            ! 0 = incoming ice
-      !                                ! 1 = outgoing ice
-      INTEGER  ::   jb, jk, jgrd, jl   ! dummy loop indices
-      INTEGER  ::   ji, jj, ii, ij     ! local scalar
-      REAL(wp) ::   zwgt, zwgt1        ! local scalar
-      REAL(wp) ::   ztmelts, zdh
-      !!------------------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 )   CALL timing_start('bdy_ice_frs')
-      !
-      jgrd = 1      ! Everything is at T-points here
-      !
-#if defined key_lim3
-
-      DO jl = 1, jpl
-         DO jb = 1, idx%nblenrim(jgrd)
-            ji    = idx%nbi(jb,jgrd)
-            jj    = idx%nbj(jb,jgrd)
-            zwgt  = idx%nbw(jb,jgrd)
-            zwgt1 = 1.e0 - idx%nbw(jb,jgrd)
-            a_i (ji,jj,jl) = ( a_i (ji,jj,jl) * zwgt1 + dta%a_i (jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Leads fraction 
-            ht_i(ji,jj,jl) = ( ht_i(ji,jj,jl) * zwgt1 + dta%ht_i(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Ice depth 
-            ht_s(ji,jj,jl) = ( ht_s(ji,jj,jl) * zwgt1 + dta%ht_s(jb,jl) * zwgt ) * tmask(ji,jj,1)  ! Snow depth
-
-            ! -----------------
-            ! Pathological case
-            ! -----------------
-            ! In case a) snow load would be in excess or b) ice is coming into a warmer environment that would lead to 
-            ! very large transformation from snow to ice (see limthd_dh.F90)
-
-            ! Then, a) transfer the snow excess into the ice (different from limthd_dh)
-            zdh = MAX( 0._wp, ( rhosn * ht_s(ji,jj,jl) + ( rhoic - rau0 ) * ht_i(ji,jj,jl) ) * r1_rau0 )
-            ! Or, b) transfer all the snow into ice (if incoming ice is likely to melt as it comes into a warmer environment)
-            !zdh = MAX( 0._wp, ht_s(ji,jj,jl) * rhosn / rhoic )
-
-            ! recompute ht_i, ht_s
-            ht_i(ji,jj,jl) = MIN( hi_max(jl), ht_i(ji,jj,jl) + zdh )
-            ht_s(ji,jj,jl) = MAX( 0._wp, ht_s(ji,jj,jl) - zdh * rhoic / rhosn ) 
-
-         ENDDO
-         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )
-         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy )
-         CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy )
-      ENDDO
-      ! retrieve at_i
-      at_i(:,:) = 0._wp
-      DO jl = 1, jpl
-         at_i(:,:) = a_i(:,:,jl) + at_i(:,:)
-      END DO
-
-      DO jl = 1, jpl
-         DO jb = 1, idx%nblenrim(jgrd)
-            ji    = idx%nbi(jb,jgrd)
-            jj    = idx%nbj(jb,jgrd)
-
-            ! condition on ice thickness depends on the ice velocity
-            ! if velocity is outward (strictly), then ice thickness, volume... must be equal to adjacent values
-            jpbound = 0   ;   ii = ji   ;   ij = jj
-            !
-            IF( u_ice(ji+1,jj  ) < 0. .AND. umask(ji-1,jj  ,1) == 0. ) jpbound = 1; ii = ji+1; ij = jj
-            IF( u_ice(ji-1,jj  ) > 0. .AND. umask(ji+1,jj  ,1) == 0. ) jpbound = 1; ii = ji-1; ij = jj
-            IF( v_ice(ji  ,jj+1) < 0. .AND. vmask(ji  ,jj-1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj+1
-            IF( v_ice(ji  ,jj-1) > 0. .AND. vmask(ji  ,jj+1,1) == 0. ) jpbound = 1; ii = ji  ; ij = jj-1
-            !
-            IF( nn_ice_lim_dta(ib_bdy) == 0 ) jpbound = 0; ii = ji; ij = jj   ! case ice boundaries = initial conditions
-            !                                                                 !      do not make state variables dependent on velocity
-            !
-            rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ii,ij) - 0.01 ) ) ! 0 if no ice
-            !
-            ! concentration and thickness
-            a_i (ji,jj,jl) = a_i (ii,ij,jl) * rswitch
-            ht_i(ji,jj,jl) = ht_i(ii,ij,jl) * rswitch
-            ht_s(ji,jj,jl) = ht_s(ii,ij,jl) * rswitch
-            !
-            ! Ice and snow volumes
-            v_i(ji,jj,jl) = ht_i(ji,jj,jl) * a_i(ji,jj,jl)
-            v_s(ji,jj,jl) = ht_s(ji,jj,jl) * a_i(ji,jj,jl)
-            !
-            SELECT CASE( jpbound )
-            !
-            CASE( 0 )   ! velocity is inward
-               !
-               ! Ice salinity, age, temperature
-               sm_i(ji,jj,jl)   = rswitch * rn_ice_sal(ib_bdy)  + ( 1.0 - rswitch ) * rn_simin
-               oa_i(ji,jj,jl)   = rswitch * rn_ice_age(ib_bdy) * a_i(ji,jj,jl)
-               t_su(ji,jj,jl)   = rswitch * rn_ice_tem(ib_bdy)  + ( 1.0 - rswitch ) * rn_ice_tem(ib_bdy)
-               DO jk = 1, nlay_s
-                  t_s(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rt0
-               END DO
-               DO jk = 1, nlay_i
-                  t_i(ji,jj,jk,jl) = rswitch * rn_ice_tem(ib_bdy) + ( 1.0 - rswitch ) * rt0 
-                  s_i(ji,jj,jk,jl) = rswitch * rn_ice_sal(ib_bdy) + ( 1.0 - rswitch ) * rn_simin
-               END DO
-               !
-            CASE( 1 )   ! velocity is outward
-               !
-               ! Ice salinity, age, temperature
-               sm_i(ji,jj,jl)   = rswitch * sm_i(ii,ij,jl)  + ( 1.0 - rswitch ) * rn_simin
-               oa_i(ji,jj,jl)   = rswitch * oa_i(ii,ij,jl)
-               t_su(ji,jj,jl)   = rswitch * t_su(ii,ij,jl)  + ( 1.0 - rswitch ) * rt0
-               DO jk = 1, nlay_s
-                  t_s(ji,jj,jk,jl) = rswitch * t_s(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rt0
-               END DO
-               DO jk = 1, nlay_i
-                  t_i(ji,jj,jk,jl) = rswitch * t_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rt0
-                  s_i(ji,jj,jk,jl) = rswitch * s_i(ii,ij,jk,jl) + ( 1.0 - rswitch ) * rn_simin
-               END DO
-               !
-            END SELECT
-            !
-            IF( nn_icesal == 1 ) THEN     ! constant salinity : overwrite rn_icesal
-               sm_i(ji,jj  ,jl) = rn_icesal
-               s_i (ji,jj,:,jl) = rn_icesal
-            ENDIF
-            !
-            ! contents
-            smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl)
-            DO jk = 1, nlay_s
-               ! Snow energy of melting
-               e_s(ji,jj,jk,jl) = rswitch * rhosn * ( cpic * ( rt0 - t_s(ji,jj,jk,jl) ) + lfus )
-               ! Multiply by volume, so that heat content in J/m2
-               e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * v_s(ji,jj,jl) * r1_nlay_s
-            END DO
-            DO jk = 1, nlay_i
-               ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rt0 !Melting temperature in K                  
-               ! heat content per unit volume
-               e_i(ji,jj,jk,jl) = rswitch * rhoic * &
-                  (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) &
-                  +   lfus    * ( 1.0 - (ztmelts-rt0) / MIN((t_i(ji,jj,jk,jl)-rt0),-epsi20) ) &
-                  - rcp      * ( ztmelts - rt0 ) )
-               ! Mutliply by ice volume, and divide by number of layers to get heat content in J/m2
-               e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) * r1_nlay_i
-            END DO
-            !
-         END DO
-         !
-         CALL lbc_bdy_lnk(  a_i(:,:,jl), 'T', 1., ib_bdy )
-         CALL lbc_bdy_lnk( ht_i(:,:,jl), 'T', 1., ib_bdy )
-         CALL lbc_bdy_lnk( ht_s(:,:,jl), 'T', 1., ib_bdy )
-         CALL lbc_bdy_lnk(  v_i(:,:,jl), 'T', 1., ib_bdy )
-         CALL lbc_bdy_lnk(  v_s(:,:,jl), 'T', 1., ib_bdy )
-         !
-         CALL lbc_bdy_lnk( smv_i(:,:,jl), 'T', 1., ib_bdy )
-         CALL lbc_bdy_lnk(  sm_i(:,:,jl), 'T', 1., ib_bdy )
-         CALL lbc_bdy_lnk(  oa_i(:,:,jl), 'T', 1., ib_bdy )
-         CALL lbc_bdy_lnk(  t_su(:,:,jl), 'T', 1., ib_bdy )
-         DO jk = 1, nlay_s
-            CALL lbc_bdy_lnk(t_s(:,:,jk,jl), 'T', 1., ib_bdy )
-            CALL lbc_bdy_lnk(e_s(:,:,jk,jl), 'T', 1., ib_bdy )
-         END DO
-         DO jk = 1, nlay_i
-            CALL lbc_bdy_lnk(t_i(:,:,jk,jl), 'T', 1., ib_bdy )
-            CALL lbc_bdy_lnk(e_i(:,:,jk,jl), 'T', 1., ib_bdy )
-         END DO
-         !
-      END DO !jl
-      !
-#endif
-      !      
-      IF( nn_timing == 1 )   CALL timing_stop('bdy_ice_frs')
-      !
-   END SUBROUTINE bdy_ice_frs
-
-
-   SUBROUTINE bdy_ice_lim_dyn( cd_type )
-      !!------------------------------------------------------------------------------
-      !!                 ***  SUBROUTINE bdy_ice_lim_dyn  ***
-      !!                    
-      !! ** Purpose : Apply dynamics boundary conditions for sea-ice in the cas of unstructured open boundaries.
-      !!              u_ice and v_ice are equal to the value of the adjacent grid point if this latter is not ice free
-      !!              if adjacent grid point is ice free, then u_ice and v_ice are equal to ocean velocities
-      !!
-      !! 2013-06 : C. Rousset
-      !!------------------------------------------------------------------------------
-      CHARACTER(len=1), INTENT(in)  ::   cd_type   ! nature of velocity grid-points
-      !
-      INTEGER  ::   jb, jgrd           ! dummy loop indices
-      INTEGER  ::   ji, jj             ! local scalar
-      INTEGER  ::   ib_bdy             ! Loop index
-      REAL(wp) ::   zmsk1, zmsk2, zflag
-      !!------------------------------------------------------------------------------
-      !
-      IF( nn_timing == 1 ) CALL timing_start('bdy_ice_lim_dyn')
-      !
-      DO ib_bdy=1, nb_bdy
-         !
-         SELECT CASE( cn_ice_lim(ib_bdy) )
-         !
-         CASE('none')
-            CYCLE
-            !
-         CASE('frs')
-            !
-            IF( nn_ice_lim_dta(ib_bdy) == 0 ) CYCLE            ! case ice boundaries = initial conditions 
-            !                                                  !      do not change ice velocity (it is only computed by rheology)
-            SELECT CASE ( cd_type )
-            !     
-            CASE ( 'U' )  
-               jgrd = 2      ! u velocity
-               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(jgrd)
-                  ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
-                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
-                  zflag = idx_bdy(ib_bdy)%flagu(jb,jgrd)
-                  !
-                  IF ( ABS( zflag ) == 1. ) THEN  ! eastern and western boundaries
-                     ! one of the two zmsk is always 0 (because of zflag)
-                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji+1,jj) ) ) ! 0 if no ice
-                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji-1,jj) ) ) ! 0 if no ice
-                     !  
-                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)
-                     u_ice (ji,jj) = u_ice(ji+1,jj) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + &
-                        &            u_ice(ji-1,jj) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + &
-                        &            u_oce(ji  ,jj) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
-                  ELSE                             ! everywhere else
-                     !u_ice(ji,jj) = u_oce(ji,jj)
-                     u_ice(ji,jj) = 0._wp
-                  ENDIF
-                  ! mask ice velocities
-                  rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ji,jj) - 0.01_wp ) ) ! 0 if no ice
-                  u_ice(ji,jj) = rswitch * u_ice(ji,jj)
-                  !
-               END DO
-               CALL lbc_bdy_lnk( u_ice(:,:), 'U', -1., ib_bdy )
-               !
-            CASE ( 'V' )
-               jgrd = 3      ! v velocity
-               DO jb = 1, idx_bdy(ib_bdy)%nblenrim(jgrd)
-                  ji    = idx_bdy(ib_bdy)%nbi(jb,jgrd)
-                  jj    = idx_bdy(ib_bdy)%nbj(jb,jgrd)
-                  zflag = idx_bdy(ib_bdy)%flagv(jb,jgrd)
-                  !
-                  IF ( ABS( zflag ) == 1. ) THEN  ! northern and southern boundaries
-                     ! one of the two zmsk is always 0 (because of zflag)
-                     zmsk1 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj+1) ) ) ! 0 if no ice
-                     zmsk2 = 1._wp - MAX( 0.0_wp, SIGN ( 1.0_wp , - vt_i(ji,jj-1) ) ) ! 0 if no ice
-                     !  
-                     ! u_ice = u_ice of the adjacent grid point except if this grid point is ice-free (then u_ice = u_oce)
-                     v_ice (ji,jj) = v_ice(ji,jj+1) * 0.5_wp * ABS( zflag + 1._wp ) * zmsk1 + &
-                        &            v_ice(ji,jj-1) * 0.5_wp * ABS( zflag - 1._wp ) * zmsk2 + &
-                        &            v_oce(ji,jj  ) * ( 1._wp - MIN( 1._wp, zmsk1 + zmsk2 ) )
-                  ELSE                             ! everywhere else
-                     !v_ice(ji,jj) = v_oce(ji,jj)
-                     v_ice(ji,jj) = 0._wp
-                  ENDIF
-                  ! mask ice velocities
-                  rswitch = MAX( 0.0_wp , SIGN ( 1.0_wp , at_i(ji,jj) - 0.01 ) ) ! 0 if no ice
-                  v_ice(ji,jj) = rswitch * v_ice(ji,jj)
-                  !
-               END DO
-               CALL lbc_bdy_lnk( v_ice(:,:), 'V', -1., ib_bdy )
-               !
-            END SELECT
-            !
-         CASE DEFAULT
-            CALL ctl_stop( 'bdy_ice_lim_dyn : unrecognised option for open boundaries for ice fields' )
-         END SELECT
-         !
-      END DO
-      !
-      IF( nn_timing == 1 ) CALL timing_stop('bdy_ice_lim_dyn')
-      !
-    END SUBROUTINE bdy_ice_lim_dyn
-
-#else
-   !!---------------------------------------------------------------------------------
-   !!   Default option                                                    Empty module
-   !!---------------------------------------------------------------------------------
-CONTAINS
-   SUBROUTINE bdy_ice_lim( kt )      ! Empty routine
-      WRITE(*,*) 'bdy_ice_lim: You should not have seen this print! error?', kt
-   END SUBROUTINE bdy_ice_lim
-#endif
-
-   !!=================================================================================
-END MODULE bdyice_lim
