Index: /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_2/limrhg.F90
===================================================================
--- /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_2/limrhg.F90	(revision 6859)
+++ /branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/LIM_SRC_2/limrhg.F90	(revision 6859)
@@ -0,0 +1,707 @@
+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 lim2_evp cas
+   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation 
+   !!            3.5  !  2012-08  (R. Benshila)  AGRIF 
+   !!----------------------------------------------------------------------
+#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp )
+   !!----------------------------------------------------------------------
+   !!   'key_lim3'               OR                     LIM-3 sea-ice model
+   !!   'key_lim2' AND NOT 'key_lim2_vp'            EVP LIM-2 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   lim_rhg       : computes ice velocities
+   !!----------------------------------------------------------------------
+   USE phycst         ! Physical constant
+   USE oce     , ONLY :  snwice_mass, snwice_mass_b
+   USE par_oce        ! Ocean parameters
+   USE dom_oce        ! Ocean domain
+   USE sbc_oce        ! Surface boundary condition: ocean fields
+   USE sbc_ice        ! Surface boundary condition: ice fields
+#if defined key_lim3
+   USE ice            ! LIM-3: ice variables
+   USE dom_ice        ! LIM-3: ice domain
+   USE limitd_me      ! LIM-3: 
+#else
+   USE ice_2          ! LIM-2: ice variables
+   USE dom_ice_2      ! LIM-2: ice domain
+#endif
+   USE lbclnk         ! Lateral Boundary Condition / MPP link
+   USE lib_mpp        ! MPP library
+   USE wrk_nemo       ! work arrays
+   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 && defined key_lim2
+   USE agrif_lim2_interp
+#endif
+#if defined key_bdy
+   USE bdyice_lim
+#endif
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   lim_rhg        ! routine called by lim_dyn (or lim_dyn_2)
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011)
+   !! $Id: limrhg.F90 5888 2015-11-13 16:47:32Z clem $
+   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE lim_rhg( k_j1, k_jpj )
+      !!-------------------------------------------------------------------
+      !!                 ***  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
+      !!
+      !! References : Hunke and Dukowicz, JPO97
+      !!              Bouillon et al., Ocean Modelling 2009
+      !!-------------------------------------------------------------------
+      INTEGER, INTENT(in) ::   k_j1    ! southern j-index for ice computation
+      INTEGER, INTENT(in) ::   k_jpj   ! northern j-index for ice computation
+      !!
+      INTEGER ::   ji, jj   ! dummy loop indices
+      INTEGER ::   jter     ! local integers
+      CHARACTER (len=50) ::   charout
+      REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         !
+      REAL(wp) ::   za, zstms          ! local scalars
+      REAL(wp) ::   zc1, zc2, zc3      ! ice mass
+
+      REAL(wp) ::   dtevp , z1_dtevp              ! time step for subcycling
+      REAL(wp) ::   dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity
+      REAL(wp) ::   z0, zr, zcca, zccb            ! temporary scalars
+      REAL(wp) ::   zu_ice2, zv_ice1              !
+      REAL(wp) ::   zddc, zdtc                    ! delta on corners and on centre
+      REAL(wp) ::   zdst                          ! shear at the center of the grid point
+      REAL(wp) ::   zdsshx, zdsshy                ! term for the gradient of ocean surface
+      REAL(wp) ::   sigma1, sigma2                ! internal ice stress
+
+      REAL(wp) ::   zresm         ! Maximal error on ice velocity
+      REAL(wp) ::   zintb, zintn  ! dummy argument
+
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zpreshc          ! Ice strength on grid cell corners (zpreshc)
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld1, zfrld2   ! lead fraction on U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2   ! ice/snow mass on U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays
+      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1           ! ocean u/v component on U points                           
+      REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2           ! ocean u/v component on V points
+      REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zmask            ! mask ocean grid points
+      
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2 
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12
+      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity
+      REAL(wp), POINTER, DIMENSION(:,:) ::   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), PARAMETER               ::   zepsi = 1.0e-20_wp ! tolerance parameter
+      REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity
+      !!-------------------------------------------------------------------
+
+      CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
+      CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               )
+      CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
+      CALL wrk_alloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 )
+
+#if  defined key_lim2 && ! defined key_lim2_vp
+# if defined key_agrif
+      USE ice_2, vt_s => hsnm
+      USE ice_2, vt_i => hicm
+# else
+      vt_s => hsnm
+      vt_i => hicm
+# endif
+      at_i(:,:) = 1. - frld(:,:)
+#endif
+#if defined key_agrif && defined key_lim2 
+      CALL agrif_rhg_lim2_load      ! First interpolation of coarse values
+#endif
+      !
+      !------------------------------------------------------------------------------!
+      ! 1) Ice strength (zpresh)                                !
+      !------------------------------------------------------------------------------!
+      !
+      ! Put every vector to 0
+      delta_i(:,:) = 0._wp   ;
+      zpresh (:,:) = 0._wp   ;  
+      zpreshc(:,:) = 0._wp
+      u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp
+      divu_i (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp
+      shear_i(:,:) = 0._wp
+
+#if defined key_lim3
+      CALL lim_itd_me_icestrength( nn_icestr )      ! LIM-3: Ice strength on T-points
+#endif
+
+      DO jj = k_j1 , k_jpj       ! Ice mass and temp variables
+         DO ji = 1 , jpi
+#if defined key_lim3
+            zpresh(ji,jj) = tmask(ji,jj,1) *  strength(ji,jj)
+#endif
+#if defined key_lim2
+            zpresh(ji,jj) = tmask(ji,jj,1) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )
+#endif
+            ! zmask = 1 where there is ice or on land
+            zmask(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1)
+         END DO
+      END DO
+
+      ! Ice strength on grid cell corners (zpreshc)
+      ! needed for calculation of shear stress 
+      DO jj = k_j1+1, k_jpj-1
+         DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1)
+            zstms          =  tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) +   &
+               &              tmask(ji+1,jj,1)   * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1)   * wght(ji+1,jj+1,1,1)
+            zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) +   &
+               &               zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + zpresh(ji,jj)   * wght(ji+1,jj+1,1,1)     &
+               &             ) / MAX( zstms, zepsi )
+         END DO
+      END DO
+      CALL lbc_lnk( zpreshc(:,:), 'F', 1. )
+      !
+      !------------------------------------------------------------------------------!
+      ! 2) Wind / ocean stress, mass terms, coriolis terms
+      !------------------------------------------------------------------------------!
+      !
+      !  Wind stress, coriolis and mass terms on the sides of the squares        
+      !  zfrld1: lead fraction on U-points                                      
+      !  zfrld2: lead fraction on V-points                                     
+      !  zmass1: ice/snow mass on U-points                                    
+      !  zmass2: ice/snow mass on V-points                                   
+      !  zcorl1: Coriolis parameter on U-points                             
+      !  zcorl2: Coriolis parameter on V-points                            
+      !  (ztagnx,ztagny): wind stress on U/V points                       
+      !  v_oce1: ocean v component on u points                          
+      !  u_oce2: ocean u component on v points                         
+
+      IF( nn_ice_embd == 2 ) 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 = k_j1+1, k_jpj-1
+         DO ji = fs_2, fs_jpim1
+
+            zc1 = tmask(ji  ,jj  ,1) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) )
+            zc2 = tmask(ji+1,jj  ,1) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) )
+            zc3 = tmask(ji  ,jj+1,1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) )
+
+            zt11 = tmask(ji  ,jj,1) * e1t(ji  ,jj)
+            zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj)
+            zt21 = tmask(ji,jj  ,1) * e2t(ji,jj  )
+            zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1)
+
+            ! Leads area.
+            zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi )
+            zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi )
+
+            ! Mass, coriolis coeff. and currents
+            zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi )
+            zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi )
+            zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) )   &
+               &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi )
+            zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) )   &
+               &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi )
+            !
+            ! Ocean has no slip boundary condition
+            v_oce1(ji,jj)  = 0.5 * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji,jj)      &
+               &                   + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) )  &
+               &                   / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)  
+
+            u_oce2(ji,jj)  = 0.5 * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj)      &
+               &                   + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) )  &
+               &                   / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
+
+            ! Wind stress at U,V-point
+            ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj)
+            ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj)
+
+            ! Computation of the velocity field taking into account the ice internal interaction.
+            ! Terms that are independent of the velocity field.
+
+            ! SB On utilise maintenant le gradient de la pente de l'ocean
+            ! include it later
+
+            zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj)
+            zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj)
+
+            za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx
+            za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy
+
+         END DO
+      END DO
+
+      !
+      !------------------------------------------------------------------------------!
+      ! 3) Solution of the momentum equation, iterative procedure
+      !------------------------------------------------------------------------------!
+      !
+      ! Time step for subcycling
+      dtevp  = rdt_ice / nn_nevp
+#if defined key_lim3
+      dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice )
+#else
+      dtotel = dtevp / ( 2._wp * telast )
+#endif
+      z1_dtotel = 1._wp / ( 1._wp + dtotel )
+      z1_dtevp  = 1._wp / dtevp
+      !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)
+      ecc2 = rn_ecc * rn_ecc
+      ecci = 1. / ecc2
+
+      !-Initialise stress tensor 
+      zs1 (:,:) = stress1_i (:,:) 
+      zs2 (:,:) = stress2_i (:,:)
+      zs12(:,:) = stress12_i(:,:)
+
+      !                                               !----------------------!
+      DO jter = 1 , nn_nevp                           !    loop over jter    !
+         !                                            !----------------------!        
+         DO jj = k_j1, k_jpj-1
+            zu_ice(:,jj) = u_ice(:,jj)    ! velocity at previous time step
+            zv_ice(:,jj) = v_ice(:,jj)
+         END DO
+
+         DO jj = k_j1+1, k_jpj-1
+            DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to zmask
+
+               !  
+               !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002)
+               !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells
+               !- zds(:,:): shear on northeast corner of grid cells
+               !
+               !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded, 
+               !                      there are many repeated calculations. 
+               !                      Speed could be improved by regrouping terms. For
+               !                      the moment, however, the stress is on clarity of coding to avoid
+               !                      bugs (Martin, for Miguel).
+               !
+               !- ALSO: arrays zdt, zds and delta could 
+               !  be removed in the future to minimise memory demand.
+               !
+               !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of
+               !              grid cells, exactly as in the B grid case. For simplicity, the indexation on
+               !              the corners is the same as in the B grid.
+               !
+               !
+               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_e12t(ji,jj)
+
+               zdt(ji,jj) = ( ( 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_e12t(ji,jj)
+
+               !
+               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_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) )   &
+                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
+
+
+               v_ice1(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) )   &
+                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
+
+               u_ice2(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  ) )   &
+                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
+            END DO
+         END DO
+
+         CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. )      ! lateral boundary cond.
+         
+         DO jj = k_j1+1, k_jpj-1
+            DO ji = fs_2, fs_jpim1
+
+               !- Calculate Delta at centre of grid cells
+               zdst          = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   &
+                  &            + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   &
+                  &            ) * r1_e12t(ji,jj)
+
+               delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )  
+               delta_i(ji,jj) = delta + rn_creepl
+
+               !- Calculate Delta on corners
+               zddc  = (  ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
+                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
+                  &    ) * r1_e12f(ji,jj)
+
+               zdtc  = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  &
+                  &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  &
+                  &    ) * r1_e12f(ji,jj)
+
+               zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl
+
+               !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide).
+               zs1(ji,jj)  = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj)   * zpresh(ji,jj)   &
+                  &          ) * z1_dtotel
+               zs2(ji,jj)  = ( zs2 (ji,jj) + dtotel *         ecci * zdt(ji,jj) / delta_i(ji,jj)   * zpresh(ji,jj)   &
+                  &          ) * z1_dtotel
+               !-Calculate stress tensor component zs12 at corners
+               zs12(ji,jj) = ( zs12(ji,jj) + dtotel *         ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj)  &
+                  &          ) * z1_dtotel 
+
+            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 = k_j1+1, k_jpj-1
+            DO ji = fs_2, fs_jpim1
+               !- contribution of zs1, zs2 and zs12 to zf1
+               zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  &
+                  &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj)          &
+                  &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj)  &
+                  &                ) * r1_e12u(ji,jj)
+               ! contribution of zs1, zs2 and zs12 to zf2
+               zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  &
+                  &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj)          &
+                  &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj)  &
+                  &               )  * r1_e12v(ji,jj)
+            END DO
+         END DO
+         !
+         ! Computation of ice velocity
+         !
+         ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly.
+         !
+         IF (MOD(jter,2).eq.0) THEN 
+
+            DO jj = k_j1+1, k_jpj-1
+               DO ji = fs_2, fs_jpim1
+                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
+                  z0           = zmass1(ji,jj) * z1_dtevp
+
+                  ! SB modif because ocean has no slip boundary condition
+                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji  ,jj)     &
+                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
+                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
+                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
+                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
+                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
+                  zcca         = z0 + za
+                  zccb         = zcorl1(ji,jj)
+                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
+               END DO
+            END DO
+
+            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
+#if defined key_agrif && defined key_lim2
+            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
+#endif
+#if defined key_bdy
+         CALL bdy_ice_lim_dyn( 'U' )
+#endif         
+
+            DO jj = k_j1+1, k_jpj-1
+               DO ji = fs_2, fs_jpim1
+
+                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
+                  z0           = zmass2(ji,jj) * z1_dtevp
+                  ! SB modif because ocean has no slip boundary condition
+                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
+                     &                 + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
+                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
+                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  & 
+                     &                         ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) )
+                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
+                  zcca         = z0 + za
+                  zccb         = zcorl2(ji,jj)
+                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
+               END DO
+            END DO
+
+            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
+#if defined key_agrif && defined key_lim2
+            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
+#endif
+#if defined key_bdy
+         CALL bdy_ice_lim_dyn( 'V' )
+#endif         
+
+         ELSE 
+            DO jj = k_j1+1, k_jpj-1
+               DO ji = fs_2, fs_jpim1
+                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)
+                  z0           = zmass2(ji,jj) * z1_dtevp
+                  ! SB modif because ocean has no slip boundary condition
+                  zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       &
+                     &                  +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   &
+                     &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)   
+
+                  za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  &
+                     &                         ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) )
+                  zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj)
+                  zcca         = z0 + za
+                  zccb         = zcorl2(ji,jj)
+                  v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch
+               END DO
+            END DO
+
+            CALL lbc_lnk( v_ice(:,:), 'V', -1. )
+#if defined key_agrif && defined key_lim2
+            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )
+#endif
+#if defined key_bdy
+         CALL bdy_ice_lim_dyn( 'V' )
+#endif         
+
+            DO jj = k_j1+1, k_jpj-1
+               DO ji = fs_2, fs_jpim1
+                  rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)
+                  z0           = zmass1(ji,jj) * z1_dtevp
+                  zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji,jj)       &
+                     &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   &
+                     &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
+
+                  za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  &
+                     &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) )
+                  zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj)
+                  zcca         = z0 + za
+                  zccb         = zcorl1(ji,jj)
+                  u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 
+               END DO
+            END DO
+
+            CALL lbc_lnk( u_ice(:,:), 'U', -1. )
+#if defined key_agrif && defined key_lim2
+            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )
+#endif
+#if defined key_bdy
+         CALL bdy_ice_lim_dyn( 'U' )
+#endif         
+
+         ENDIF
+         
+         IF(ln_ctl) THEN
+            !---  Convergence test.
+            DO jj = k_j1+1 , k_jpj-1
+               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) )
+            END DO
+            zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) )
+            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain
+         ENDIF
+
+         !                                                ! ==================== !
+      END DO                                              !  end loop over jter  !
+      !                                                   ! ==================== !
+      !
+      !------------------------------------------------------------------------------!
+      ! 4) Prevent ice velocities when the ice is thin
+      !------------------------------------------------------------------------------!
+      ! If the ice volume is below zvmin then ice velocity should equal the
+      ! ocean velocity. This prevents high velocity when ice is thin
+      DO jj = k_j1+1, k_jpj-1
+         DO ji = fs_2, fs_jpim1
+            IF ( vt_i(ji,jj) <= zvmin ) THEN
+               u_ice(ji,jj) = u_oce(ji,jj)
+               v_ice(ji,jj) = v_oce(ji,jj)
+            ENDIF
+         END DO
+      END DO
+
+      CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. )
+
+#if defined key_agrif && defined key_lim2
+      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' )
+      CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' )
+#endif
+#if defined key_bdy
+      CALL bdy_ice_lim_dyn( 'U' )
+      CALL bdy_ice_lim_dyn( 'V' )
+#endif         
+
+      DO jj = k_j1+1, k_jpj-1 
+         DO ji = fs_2, fs_jpim1
+            IF ( vt_i(ji,jj) <= zvmin ) THEN
+               v_ice1(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) )   &
+                  &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)
+
+               u_ice2(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  ) )   &
+                  &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)
+            ENDIF 
+         END DO
+      END DO
+
+      CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. )
+
+      ! Recompute delta, shear and div, inputs for mechanical redistribution 
+      DO jj = k_j1+1, k_jpj-1
+         DO ji = fs_2, jpim1   !RB bug no vect opt due to zmask
+            !- divu_i(:,:), zdt(:,:): divergence and tension at centre 
+            !- zds(:,:): shear on northeast corner of grid cells
+            IF ( vt_i(ji,jj) <= zvmin ) THEN
+
+               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_e12t(ji,jj)
+
+               zdt(ji,jj) = ( ( 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_e12t(ji,jj)
+               !
+               ! SB modif because ocean has no slip boundary condition 
+               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_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) )                                     &
+                  &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1)
+
+               zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )    &
+                  &   + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1) ) * r1_e12t(ji,jj)
+
+               delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )  
+               delta_i(ji,jj) = delta + rn_creepl
+            
+            ENDIF
+         END DO
+      END DO
+      !
+      !------------------------------------------------------------------------------!
+      ! 5) Store stress tensor and its invariants
+      !------------------------------------------------------------------------------!
+      ! * Invariants of the stress tensor are required for limitd_me
+      !   (accelerates convergence and improves stability)
+      DO jj = k_j1+1, k_jpj-1
+         DO ji = fs_2, fs_jpim1
+            zdst           = (  e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)  &   
+               &              + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj) 
+            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst )
+         END DO
+      END DO
+
+      ! Lateral boundary condition
+      CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1.,  shear_i(:,:), 'T', 1. )
+
+      ! * Store the stress tensor for the next time step
+      stress1_i (:,:) = zs1 (:,:)
+      stress2_i (:,:) = zs2 (:,:)
+      stress12_i(:,:) = zs12(:,:)
+
+      !
+      !------------------------------------------------------------------------------!
+      ! 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
+         CALL prt_ctl_info('lim_rhg  : numit  :',ivar1=numit)
+         CALL prt_ctl_info('lim_rhg  : nwrite :',ivar1=nwrite)
+         CALL prt_ctl_info('lim_rhg  : MOD    :',ivar1=MOD(numit,nwrite))
+         IF( MOD(numit,nwrite) .EQ. 0 ) THEN
+            WRITE(charout,FMT="('lim_rhg  :', I4, I6, I1, I1, A10)") 1000, numit, 0, 0, ' ch. ell. '
+            CALL prt_ctl_info(charout)
+            DO jj = k_j1+1, k_jpj-1
+               DO ji = 2, jpim1
+                  IF (zpresh(ji,jj) > 1.0) THEN
+                     sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
+                     sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(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, numit, 0, 0, ' ch. ell. '
+            CALL prt_ctl_info(charout)
+         ENDIF
+      ENDIF
+      !
+      CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct )
+      CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               )
+      CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  )
+      CALL wrk_dealloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 )
+
+   END SUBROUTINE lim_rhg
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option          Dummy module           NO LIM sea-ice model
+   !!----------------------------------------------------------------------
+CONTAINS
+   SUBROUTINE lim_rhg( k1 , k2 )         ! Dummy routine
+      WRITE(*,*) 'lim_rhg: You should not have seen this print! error?', k1, k2
+   END SUBROUTINE lim_rhg
+#endif
+
+   !!==============================================================================
+END MODULE limrhg
