Index: /NEMO/branches/2019/dev_r11842_SI3-10_EAP/cfgs/SHARED/namelist_ice_ref
===================================================================
--- /NEMO/branches/2019/dev_r11842_SI3-10_EAP/cfgs/SHARED/namelist_ice_ref	(revision 11866)
+++ /NEMO/branches/2019/dev_r11842_SI3-10_EAP/cfgs/SHARED/namelist_ice_ref	(revision 11867)
@@ -97,4 +97,5 @@
       rn_relast     =   0.333         !     ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 
                                       !        advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300)
+   ln_rhg_EAP       = .false.          !  EVP rheology
 /
 !------------------------------------------------------------------------------
Index: /NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg.F90
===================================================================
--- /NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg.F90	(revision 11866)
+++ /NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg.F90	(revision 11867)
@@ -17,4 +17,5 @@
    USE ice            ! sea-ice: variables
    USE icedyn_rhg_evp ! sea-ice: EVP rheology
+   USE icedyn_rhg_eap ! sea-ice: EAP rheology
    USE icectl         ! sea-ice: control prints
    !
@@ -33,8 +34,9 @@
    !                                        ! associated indices:
    INTEGER, PARAMETER ::   np_rhgEVP = 1   ! EVP rheology
-!! INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology
+   INTEGER, PARAMETER ::   np_rhgEAP = 2   ! EAP rheology
 
    ! ** namelist (namrhg) **
    LOGICAL ::   ln_rhg_EVP       ! EVP rheology
+   LOGICAL ::   ln_rhg_EAP       ! EVP rheology
    !
    !! * Substitutions
@@ -81,8 +83,13 @@
          CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
          !         
+         !                             !----------------------------!
+      CASE( np_rhgEAP )                ! Elasto-Anisotropic-Plastic !
+         !                             !----------------------------!
+         CALL ice_dyn_rhg_eap( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i )
       END SELECT
       !
-      IF( lrst_ice ) THEN                       !* write EVP fields in the restart file
-         IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt )
+      IF( lrst_ice ) THEN                      
+         IF( ln_rhg_EVP )   CALL rhg_evp_rst( 'WRITE', kt ) !* write EVP fields in the restart file
+         IF( ln_rhg_EAP )   CALL rhg_eap_rst( 'WRITE', kt ) !* write EAP fields in the restart file
       ENDIF
       !
@@ -110,5 +117,5 @@
       INTEGER ::   ios, ioptio   ! Local integer output status for namelist read
       !!
-      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast
+      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, ln_rhg_EAP
       !!-------------------------------------------------------------------
       !
@@ -132,4 +139,5 @@
          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,*) '      rheology EAP (icedyn_rhg_eap)                        ln_rhg_EAP = ', ln_rhg_EAP
       ENDIF
       !
@@ -137,8 +145,9 @@
       ioptio = 0 
       IF( ln_rhg_EVP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEVP    ;   ENDIF
-!!    IF( ln_rhg_EAP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEAP    ;   ENDIF
+      IF( ln_rhg_EAP ) THEN   ;   ioptio = ioptio + 1   ;   nice_rhg = np_rhgEAP    ;   ENDIF
       IF( ioptio /= 1 )   CALL ctl_stop( 'ice_dyn_rhg_init: choose one and only one ice rheology' )
       !
       IF( ln_rhg_EVP  )   CALL rhg_evp_rst( 'READ' )  !* read or initialize all required files
+      IF( ln_rhg_EAP  )   CALL rhg_eap_rst( 'READ' )  !* read or initialize all required files
       !
    END SUBROUTINE ice_dyn_rhg_init
Index: /NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90
===================================================================
--- /NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90	(revision 11867)
+++ /NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90	(revision 11867)
@@ -0,0 +1,919 @@
+MODULE icedyn_rhg_eap
+   !!======================================================================
+   !!                     ***  MODULE  icedyn_rhg_eap  ***
+   !!   Sea-Ice dynamics : rheology Elasto-Viscous-Plastic
+   !!======================================================================
+   !! History :   -   !  2007-03  (M.A. Morales Maqueda, S. Bouillon) Original code
+   !!            3.0  !  2008-03  (M. Vancoppenolle) adaptation to new model
+   !!             -   !  2008-11  (M. Vancoppenolle, S. Bouillon, Y. Aksenov) add surface tilt in ice rheolohy 
+   !!            3.3  !  2009-05  (G.Garric)    addition of the evp case
+   !!            3.4  !  2011-01  (A. Porter)   dynamical allocation 
+   !!            3.5  !  2012-08  (R. Benshila) AGRIF
+   !!            3.6  !  2016-06  (C. Rousset)  Rewriting + landfast ice + mEVP (Bouillon 2013)
+   !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017)
+   !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube]
+   !!----------------------------------------------------------------------
+#if defined key_si3
+   !!----------------------------------------------------------------------
+   !!   'key_si3'                                       SI3 sea-ice model
+   !!----------------------------------------------------------------------
+   !!   ice_dyn_rhg_eap : computes ice velocities from EVP rheology
+   !!   rhg_eap_rst     : read/write EVP fields in ice restart
+   !!----------------------------------------------------------------------
+   USE phycst         ! Physical constant
+   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            ! sea-ice: ice variables
+   USE icevar         ! ice_var_sshdyn
+   USE icedyn_rdgrft  ! sea-ice: ice strength
+   USE bdy_oce , ONLY : ln_bdy 
+   USE bdyice 
+#if defined key_agrif
+   USE agrif_ice_interp
+#endif
+   !
+   USE in_out_manager ! I/O manager
+   USE iom            ! I/O manager library
+   USE lib_mpp        ! MPP library
+   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero)
+   USE lbclnk         ! lateral boundary conditions (or mpp links)
+   USE prtctl         ! Print control
+
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_dyn_rhg_eap   ! called by icedyn_rhg.F90
+   PUBLIC   rhg_eap_rst       ! called by icedyn_rhg.F90
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+   !!----------------------------------------------------------------------
+   !! NEMO/ICE 4.0 , NEMO Consortium (2018)
+   !! $Id: icedyn_rhg_eap.F90 11536 2019-09-11 13:54:18Z smasson $
+   !! Software governed by the CeCILL license (see ./LICENSE)
+   !!----------------------------------------------------------------------
+CONTAINS
+
+   SUBROUTINE ice_dyn_rhg_eap( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i )
+      !!-------------------------------------------------------------------
+      !!                 ***  SUBROUTINE ice_dyn_rhg_eap  ***
+      !!                             EAP-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-anisotropic-plastic (EAP) law including shear 
+      !!  strength and a bulk rheology .	
+      !!
+      !!  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   : 0) compute mask at F point
+      !!              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) Recompute delta, shear and divergence
+      !!                 (which are inputs of the ITD) & store stress
+      !!                 for the next time step
+      !!              5) Diagnostics including charge ellipse
+      !!
+      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017)
+      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters).
+      !!              This is an upgraded version of mEVP from Bouillon et al. 2013
+      !!              (i.e. more stable and better convergence)
+      !!
+      !! References : Hunke and Dukowicz, JPO97
+      !!              Bouillon et al., Ocean Modelling 2009
+      !!              Bouillon et al., Ocean Modelling 2013
+      !!              Kimmritz et al., Ocean Modelling 2016 & 2017
+      !!-------------------------------------------------------------------
+      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i   !
+      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i  , pdivu_i   , pdelta_i      !
+      !!
+      INTEGER ::   ji, jj       ! dummy loop indices
+      INTEGER ::   jter         ! local integers
+      !
+      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) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
+      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
+      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
+      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars
+      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
+      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast
+      !
+      REAL(wp) ::   zresm                                               ! Maximal error on ice velocity
+      REAL(wp) ::   zintb, zintn                                        ! dummy argument
+      REAL(wp) ::   zfac_x, zfac_y
+      REAL(wp) ::   zshear, zdum1, zdum2
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on 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) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                           
+      !
+      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) ::   zsshdyn                         ! array used for the calculation of ice surface slope:
+      !                                                                 !    ocean surface (ssh_m) if ice is not embedded
+      !                                                                 !    ice bottom surface if ice is embedded   
+      REAL(wp), DIMENSION(jpi,jpj) ::   zfU  , zfV                      ! internal stresses
+      REAL(wp), DIMENSION(jpi,jpj) ::   zspgU, zspgV                    ! surface pressure gradient at U/V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zCorU, zCorV                    ! Coriolis stress array
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_ai, ztauy_ai              ! ice-atm. stress at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_oi, ztauy_oi              ! ice-ocean stress at U-V points
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_bi, ztauy_bi              ! ice-OceanBottom stress at U-V points (landfast)
+      REAL(wp), DIMENSION(jpi,jpj) ::   ztaux_base, ztauy_base          ! ice-bottom stress at U-V points (landfast)
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! 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 becomes very small
+      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small
+      !! --- diags
+      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3
+      !! --- SIMIP diags
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_ice ! X-component of ice mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_ice ! Y-component of ice mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xmtrp_snw ! X-component of snow mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_ymtrp_snw ! Y-component of snow mass transport (kg/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_xatrp     ! X-component of area transport (m2/s)
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zdiag_yatrp     ! Y-component of area transport (m2/s)      
+      !!-------------------------------------------------------------------
+
+      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_eap: EAP sea-ice rheology'
+      !
+!!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
+      !------------------------------------------------------------------------------!
+      ! 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( 'icedyn_rhg_eap', 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( 'icedyn_rhg_eap', 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)
+      IF( .NOT. ln_aEVP ) THEN
+         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
+         zalph2 = zalph1 * z1_ecc2
+
+         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+      ENDIF
+         
+      ! Initialise stress tensor 
+      zs1 (:,:) = pstress1_i (:,:) 
+      zs2 (:,:) = pstress2_i (:,:)
+      zs12(:,:) = pstress12_i(:,:)
+
+      ! Ice strength
+      CALL ice_strength
+
+      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010)
+      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile
+      ELSE                         ;   zkt = 0._wp
+      ENDIF
+      !
+      !------------------------------------------------------------------------------!
+      ! 2) Wind / ocean stress, mass terms, coriolis terms
+      !------------------------------------------------------------------------------!
+      ! sea surface height
+      !    embedded sea ice: compute representative ice top surface
+      !    non-embedded sea ice: use ocean surface for slope calculation
+      zsshdyn(:,:) = ice_var_sshdyn( ssh_m, snwice_mass, snwice_mass_b)
+
+      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 = ( rhos * vt_s(ji  ,jj  ) + rhoi * vt_i(ji  ,jj  ) )
+            zm2 = ( rhos * vt_s(ji+1,jj  ) + rhoi * vt_i(ji+1,jj  ) )
+            zm3 = ( rhos * vt_s(ji  ,jj+1) + rhoi * 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.25_wp * ( v_oce(ji,jj) + v_oce(ji,jj-1) + v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * umask(ji,jj,1)
+            u_oceV(ji,jj)   = 0.25_wp * ( u_oce(ji,jj) + u_oce(ji-1,jj) + u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * vmask(ji,jj,1)
+
+            ! Coriolis at T points (m*f)
+            zmf(ji,jj)      = zm1 * ff_t(ji,jj)
+
+            ! dt/m at T points (for alpha and beta coefficients)
+            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin )
+            
+            ! m/dt
+            zmU_t(ji,jj)    = zmassU * z1_dtevp
+            zmV_t(ji,jj)    = zmassV * z1_dtevp
+            
+            ! Drag ice-atm.
+            ztaux_ai(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj)
+            ztauy_ai(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 * ( zsshdyn(ji+1,jj) - zsshdyn(ji,jj) ) * r1_e1u(ji,jj)
+            zspgV(ji,jj)    = - zmassV * grav * ( zsshdyn(ji,jj+1) - zsshdyn(ji,jj) ) * r1_e2v(ji,jj)
+
+            ! masks
+            zmsk00x(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice
+            zmsk00y(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice
+
+            ! switches
+            IF( zmassU <= zmmin .AND. zaU(ji,jj) <= zamin ) THEN   ;   zmsk01x(ji,jj) = 0._wp
+            ELSE                                                   ;   zmsk01x(ji,jj) = 1._wp   ;   ENDIF
+            IF( zmassV <= zmmin .AND. zaV(ji,jj) <= zamin ) THEN   ;   zmsk01y(ji,jj) = 0._wp
+            ELSE                                                   ;   zmsk01y(ji,jj) = 1._wp   ;   ENDIF
+
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'icedyn_rhg_eap', zmf, 'T', 1., zdt_m, 'T', 1. )
+      !
+      !                                  !== Landfast ice parameterization ==!
+      !
+      IF( ln_landfast_L16 ) THEN         !-- Lemieux 2016
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               ! ice thickness at U-V points
+               zvU = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1)
+               zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1)
+               ! ice-bottom stress at U points
+               zvCr = zaU(ji,jj) * rn_depfra * hu_n(ji,jj)
+               ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
+               ! ice-bottom stress at V points
+               zvCr = zaV(ji,jj) * rn_depfra * hv_n(ji,jj)
+               ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) )
+               ! ice_bottom stress at T points
+               zvCr = at_i(ji,jj) * rn_depfra * ht_n(ji,jj)
+               tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) )
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_rhg_eap', tau_icebfr(:,:), 'T', 1. )
+         !
+      ELSE                               !-- no landfast
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1
+               ztaux_base(ji,jj) = 0._wp
+               ztauy_base(ji,jj) = 0._wp
+            END DO
+         END DO
+      ENDIF
+
+      !------------------------------------------------------------------------------!
+      ! 3) Solution of the momentum equation, iterative procedure
+      !------------------------------------------------------------------------------!
+      !
+      !                                               ! ==================== !
+      DO jter = 1 , nn_nevp                           !    loop over jter    !
+         !                                            ! ==================== !        
+         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1
+         !
+!!$         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( 'icedyn_rhg_eap', zds, 'F', 1. )
+
+         DO jj = 2, jpj    ! loop to jpi,jpj to avoid making a communication for zs1,zs2,zs12
+            DO ji = 2, jpi ! 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 )
+
+               ! alpha & beta for aEVP
+               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
+               !   alpha = beta = sqrt(4*gamma)
+               IF( ln_aEVP ) THEN
+                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+                  zalph2   = zalph1
+                  z1_alph2 = z1_alph1
+               ENDIF
+               
+               ! stress at T points (zkt/=0 if landfast)
+               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
+             
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. )
+
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+
+               ! alpha & beta for aEVP
+               IF( ln_aEVP ) THEN
+                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+                  zbeta(ji,jj) = zalph2
+               ENDIF
+               
+               ! 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 (zkt/=0 if landfast)
+               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
+
+            END DO
+         END DO
+
+         ! --- 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)
+               !
+               !                !--- ice currents at U-V point
+               v_iceU(ji,jj) = 0.25_wp * ( v_ice(ji,jj) + v_ice(ji,jj-1) + v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * umask(ji,jj,1)
+               u_iceV(ji,jj) = 0.25_wp * ( u_ice(ji,jj) + u_ice(ji-1,jj) + u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * vmask(ji,jj,1)
+               !
+            END DO
+         END DO
+         !
+         ! --- Computation of ice velocity --- !
+         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017
+         !  Bouillon et al. 2009 (eq 34-35) => stable
+         IF( MOD(jter,2) == 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  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
+                  zTauB = ztauy_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
+                  !
+                  !                 !--- Coriolis at V-points (energy conserving formulation)
+                  zCorV(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
+                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + 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
+                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &           )   * zmsk00y(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zRHS + 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
+                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &            )   * zmsk00y(ji,jj)
+                  ENDIF
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1. )
+            !
+#if defined key_agrif
+!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
+            CALL agrif_interp_ice( '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  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
+                  zTauB = ztaux_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
+                  !
+                  !                 !--- Coriolis at U-points (energy conserving formulation)
+                  zCorU(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
+                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + 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
+                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                        &           )   * zmsk00x(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zRHS + 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
+                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                        &            )   * zmsk00x(ji,jj)
+                  ENDIF
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1. )
+            !
+#if defined key_agrif
+!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
+            CALL agrif_interp_ice( '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  = 5.e-05_wp + SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) )
+                  zTauB = ztaux_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztaux_bi(ji,jj) = zTauB * u_ice(ji,jj)
+                  !
+                  !                 !--- Coriolis at U-points (energy conserving formulation)
+                  zCorU(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
+                  zRHS = zfU(ji,jj) + ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztaux_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                     u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + 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
+                        &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
+                        &           )   * zmsk00x(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                     u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zRHS + 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
+                        &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &            )   * zmsk00x(ji,jj)
+                  ENDIF
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_rhg_eap', u_ice, 'U', -1. )
+            !
+#if defined key_agrif
+!!            CALL agrif_interp_ice( 'U', jter, nn_nevp )
+            CALL agrif_interp_ice( '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  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
+                  zTauB = ztauy_base(ji,jj) / zvel
+                  !                 !--- OceanBottom-to-Ice stress
+                  ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
+                  !
+                  !                 !--- Coriolis at v-points (energy conserving formulation)
+                  zCorV(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
+                  zRHS = zfV(ji,jj) + ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)
+                  !
+                  !                 !--- landfast switch => 0 = static  friction : TauB > RHS & sign(TauB) /= sign(RHS)
+                  !                                         1 = sliding friction : TauB < RHS
+                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zRHS + ztauy_base(ji,jj) ) - SIGN( 1._wp, zRHS ) ) )
+                  !
+                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017)
+                     v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity
+                        &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part)
+                        &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + 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
+                        &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &           )   * zmsk00y(ji,jj)
+                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009)
+                     v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                                       & ! previous velocity
+                        &                                     + zRHS + 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
+                        &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin
+                        &            )   * zmsk00y(ji,jj)
+                  ENDIF
+               END DO
+            END DO
+            CALL lbc_lnk( 'icedyn_rhg_eap', v_ice, 'V', -1. )
+            !
+#if defined key_agrif
+!!            CALL agrif_interp_ice( 'V', jter, nn_nevp )
+            CALL agrif_interp_ice( '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 ) )
+!!$            CALL mpp_max( 'icedyn_rhg_eap', 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           
+      
+      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
+            pshear_i(ji,jj) = SQRT( zdt2 + zds2 )
+
+            ! divergence at T points
+            pdivu_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( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )  
+            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0
+            pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch
+
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1. )
+      
+      ! --- Store the stress tensor for the next time step --- !
+      CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
+      pstress1_i (:,:) = zs1 (:,:)
+      pstress2_i (:,:) = zs2 (:,:)
+      pstress12_i(:,:) = zs12(:,:)
+      !
+
+      !------------------------------------------------------------------------------!
+      ! 5) diagnostics
+      !------------------------------------------------------------------------------!
+      DO jj = 1, jpj
+         DO ji = 1, jpi
+            zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice
+         END DO
+      END DO
+
+      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- !
+      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. &
+         & iom_use('utau_bi') .OR. iom_use('vtau_bi') ) THEN
+         !
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', ztaux_oi, 'U', -1., ztauy_oi, 'V', -1., ztaux_ai, 'U', -1., ztauy_ai, 'V', -1., &
+            &                                  ztaux_bi, 'U', -1., ztauy_bi, 'V', -1. )
+         !
+         CALL iom_put( 'utau_oi' , ztaux_oi * zmsk00 )
+         CALL iom_put( 'vtau_oi' , ztauy_oi * zmsk00 )
+         CALL iom_put( 'utau_ai' , ztaux_ai * zmsk00 )
+         CALL iom_put( 'vtau_ai' , ztauy_ai * zmsk00 )
+         CALL iom_put( 'utau_bi' , ztaux_bi * zmsk00 )
+         CALL iom_put( 'vtau_bi' , ztauy_bi * zmsk00 )
+      ENDIF
+       
+      ! --- divergence, shear and strength --- !
+      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * zmsk00 )   ! divergence
+      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * zmsk00 )   ! shear
+      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
+
+      ! --- stress tensor --- !
+      IF( iom_use('isig1') .OR. iom_use('isig2') .OR. iom_use('isig3') .OR. iom_use('normstr') .OR. iom_use('sheastr') ) THEN
+         !
+         ALLOCATE( zsig1(jpi,jpj) , zsig2(jpi,jpj) , zsig3(jpi,jpj) )
+         !         
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               zdum1 = ( zmsk00(ji-1,jj) * pstress12_i(ji-1,jj) + zmsk00(ji  ,jj-1) * pstress12_i(ji  ,jj-1) +  &  ! stress12_i at T-point
+                  &      zmsk00(ji  ,jj) * pstress12_i(ji  ,jj) + zmsk00(ji-1,jj-1) * pstress12_i(ji-1,jj-1) )  &
+                  &    / MAX( 1._wp, zmsk00(ji-1,jj) + zmsk00(ji,jj-1) + zmsk00(ji,jj) + zmsk00(ji-1,jj-1) )
+
+               zshear = SQRT( pstress2_i(ji,jj) * pstress2_i(ji,jj) + 4._wp * zdum1 * zdum1 ) ! shear stress  
+
+               zdum2 = zmsk00(ji,jj) / MAX( 1._wp, strength(ji,jj) )
+
+!!               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) + zshear ) ! principal stress (y-direction, see Hunke & Dukowicz 2002)
+!!               zsig2(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) - zshear ) ! principal stress (x-direction, see Hunke & Dukowicz 2002)
+!!               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) ! quadratic relation linking compressive stress to shear stress
+!!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
+               zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
+               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
+               zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zsig1, 'T', 1., zsig2, 'T', 1., zsig3, 'T', 1. )
+         !
+         CALL iom_put( 'isig1' , zsig1 )
+         CALL iom_put( 'isig2' , zsig2 )
+         CALL iom_put( 'isig3' , zsig3 )
+         !
+         ! Stress tensor invariants (normal and shear stress N/m)
+         IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress
+         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress
+
+         DEALLOCATE( zsig1 , zsig2 , zsig3 )
+      ENDIF
+      
+      ! --- SIMIP --- !
+      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
+         & iom_use('corstrx') .OR. iom_use('corstry') .OR. iom_use('intstrx') .OR. iom_use('intstry') ) THEN
+         !
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zspgU, 'U', -1., zspgV, 'V', -1., &
+            &                                  zCorU, 'U', -1., zCorV, 'V', -1., zfU, 'U', -1., zfV, 'V', -1. )
+
+         CALL iom_put( 'dssh_dx' , zspgU * zmsk00 )   ! Sea-surface tilt term in force balance (x)
+         CALL iom_put( 'dssh_dy' , zspgV * zmsk00 )   ! Sea-surface tilt term in force balance (y)
+         CALL iom_put( 'corstrx' , zCorU * zmsk00 )   ! Coriolis force term in force balance (x)
+         CALL iom_put( 'corstry' , zCorV * zmsk00 )   ! Coriolis force term in force balance (y)
+         CALL iom_put( 'intstrx' , zfU   * zmsk00 )   ! Internal force term in force balance (x)
+         CALL iom_put( 'intstry' , zfV   * zmsk00 )   ! Internal force term in force balance (y)
+      ENDIF
+
+      IF(  iom_use('xmtrpice') .OR. iom_use('ymtrpice') .OR. &
+         & iom_use('xmtrpsnw') .OR. iom_use('ymtrpsnw') .OR. iom_use('xatrp') .OR. iom_use('yatrp') ) THEN
+         !
+         ALLOCATE( zdiag_xmtrp_ice(jpi,jpj) , zdiag_ymtrp_ice(jpi,jpj) , &
+            &      zdiag_xmtrp_snw(jpi,jpj) , zdiag_ymtrp_snw(jpi,jpj) , zdiag_xatrp(jpi,jpj) , zdiag_yatrp(jpi,jpj) )
+         !
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1
+               ! 2D ice mass, snow mass, area transport arrays (X, Y)
+               zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * zmsk00(ji,jj)
+               zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * zmsk00(ji,jj)
+
+               zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
+               zdiag_ymtrp_ice(ji,jj) = rhoi * zfac_y * ( vt_i(ji,jj+1) + vt_i(ji,jj) ) !        ''           Y-   ''
+
+               zdiag_xmtrp_snw(ji,jj) = rhos * zfac_x * ( vt_s(ji+1,jj) + vt_s(ji,jj) ) ! snow mass transport, X-component
+               zdiag_ymtrp_snw(ji,jj) = rhos * zfac_y * ( vt_s(ji,jj+1) + vt_s(ji,jj) ) !          ''          Y-   ''
+
+               zdiag_xatrp(ji,jj)     = zfac_x * ( at_i(ji+1,jj) + at_i(ji,jj) )        ! area transport,      X-component
+               zdiag_yatrp(ji,jj)     = zfac_y * ( at_i(ji,jj+1) + at_i(ji,jj) )        !        ''            Y-   ''
+
+            END DO
+         END DO
+
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zdiag_xmtrp_ice, 'U', -1., zdiag_ymtrp_ice, 'V', -1., &
+            &                                  zdiag_xmtrp_snw, 'U', -1., zdiag_ymtrp_snw, 'V', -1., &
+            &                                  zdiag_xatrp    , 'U', -1., zdiag_yatrp    , 'V', -1. )
+
+         CALL iom_put( 'xmtrpice' , zdiag_xmtrp_ice )   ! X-component of sea-ice mass transport (kg/s)
+         CALL iom_put( 'ymtrpice' , zdiag_ymtrp_ice )   ! Y-component of sea-ice mass transport 
+         CALL iom_put( 'xmtrpsnw' , zdiag_xmtrp_snw )   ! X-component of snow mass transport (kg/s)
+         CALL iom_put( 'ymtrpsnw' , zdiag_ymtrp_snw )   ! Y-component of snow mass transport
+         CALL iom_put( 'xatrp'    , zdiag_xatrp     )   ! X-component of ice area transport
+         CALL iom_put( 'yatrp'    , zdiag_yatrp     )   ! Y-component of ice area transport
+
+         DEALLOCATE( zdiag_xmtrp_ice , zdiag_ymtrp_ice , &
+            &        zdiag_xmtrp_snw , zdiag_ymtrp_snw , zdiag_xatrp , zdiag_yatrp )
+
+      ENDIF
+      !
+   END SUBROUTINE ice_dyn_rhg_eap
+
+
+   SUBROUTINE rhg_eap_rst( cdrw, kt )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE rhg_eap_rst  ***
+      !!                     
+      !! ** Purpose :   Read or write RHG file in restart file
+      !!
+      !! ** Method  :   use of IOM library
+      !!----------------------------------------------------------------------
+      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
+      !
+      INTEGER  ::   iter            ! local integer
+      INTEGER  ::   id1, id2, id3   ! local integers
+      !!----------------------------------------------------------------------
+      !
+      IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialize
+         !                                   ! ---------------
+         IF( ln_rstart ) THEN                   !* Read the restart file
+            !
+            id1 = iom_varid( numrir, 'stress1_i' , ldstop = .FALSE. )
+            id2 = iom_varid( numrir, 'stress2_i' , ldstop = .FALSE. )
+            id3 = iom_varid( numrir, 'stress12_i', ldstop = .FALSE. )
+            !
+            IF( MIN( id1, id2, id3 ) > 0 ) THEN      ! fields exist
+               CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  )
+               CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  )
+               CALL iom_get( numrir, jpdom_autoglo, 'stress12_i', stress12_i )
+            ELSE                                     ! start rheology from rest
+               IF(lwp) WRITE(numout,*)
+               IF(lwp) WRITE(numout,*) '   ==>>>   previous run without rheology, set stresses to 0'
+               stress1_i (:,:) = 0._wp
+               stress2_i (:,:) = 0._wp
+               stress12_i(:,:) = 0._wp
+            ENDIF
+         ELSE                                   !* Start from rest
+            IF(lwp) WRITE(numout,*)
+            IF(lwp) WRITE(numout,*) '   ==>>>   start from rest: set stresses to 0'
+            stress1_i (:,:) = 0._wp
+            stress2_i (:,:) = 0._wp
+            stress12_i(:,:) = 0._wp
+         ENDIF
+         !
+      ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
+         !                                   ! -------------------
+         IF(lwp) WRITE(numout,*) '---- rhg-rst ----'
+         iter = kt + nn_fsbc - 1             ! ice restarts are written at kt == nitrst - nn_fsbc + 1
+         !
+         CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  )
+         CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  )
+         CALL iom_rstput( iter, nitrst, numriw, 'stress12_i', stress12_i )
+         !
+      ENDIF
+      !
+   END SUBROUTINE rhg_eap_rst
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option         Empty module           NO SI3 sea-ice model
+   !!----------------------------------------------------------------------
+#endif
+
+   !!==============================================================================
+END MODULE icedyn_rhg_eap
