Index: /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/cfgs/SHARED/field_def_nemo-ice.xml
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/cfgs/SHARED/field_def_nemo-ice.xml	(revision 15569)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/cfgs/SHARED/field_def_nemo-ice.xml	(revision 15570)
@@ -81,4 +81,5 @@
           <field id="sig1_pnorm"   long_name="P-normalized 1st principal stress component"                                                                       unit=""     />
           <field id="sig2_pnorm"   long_name="P-normalized 2nd principal stress component"                                                                       unit=""     />
+          <field id="icedlt"       long_name="delta"                                                   standard_name="delta"                                     unit=""     />
           <field id="normstr"      long_name="Average normal stress in sea ice"                        standard_name="average_normal_stress"                     unit="N/m"  />
           <field id="sheastr"      long_name="Maximum shear stress in sea ice"                         standard_name="maximum_shear_stress"                      unit="N/m"  />
@@ -86,4 +87,8 @@
           <field id="icediv"       long_name="Divergence of the sea-ice velocity field"                standard_name="divergence_of_sea_ice_velocity"            unit="s-1"  />
           <field id="iceshe"       long_name="Maximum shear of sea-ice velocity field"                 standard_name="maximum_shear_of_sea_ice_velocity"         unit="s-1"  />
+          <field id="aniso"        long_name="anisotropy of sea ice floe orientation (0.5 - 1)"        standard_name="anisotropy"                                unit=""     />
+          <field id="yield11"      long_name="yield surface tensor component 11"                       standard_name="yield11"                                   unit="N/m"  />
+          <field id="yield22"      long_name="yield surface tensor component 22"                       standard_name="yield22"                                   unit="N/m"  />
+          <field id="yield12"      long_name="yield surface tensor component 12"                       standard_name="yield12"                                   unit="N/m"  />
           <field id="beta_evp"     long_name="Relaxation parameter of ice rheology (beta)"             standard_name="relaxation_parameter_of_ice_rheology"      unit=""  />	 
  
@@ -401,4 +406,7 @@
 	  <field field_ref="normstr"          name="normstr" />
 	  <field field_ref="sheastr"          name="sheastr" />
+          <field field_ref="sig1_pnorm"       name="sig1_pnorm"/>
+          <field field_ref="sig2_pnorm"       name="sig2_pnorm"/>
+	  <field field_ref="icedlt"           name="sidelta" />
 	  
 	  <!-- heat fluxes -->
Index: /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/ice.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/ice.F90	(revision 15569)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/ice.F90	(revision 15570)
@@ -150,4 +150,6 @@
    !
    !                                     !!** ice-rheology namelist (namdyn_rhg) **
+   LOGICAL , PUBLIC ::   ln_rhg_EVP       ! EVP rheology switch, used for rdgrft and rheology
+   LOGICAL , PUBLIC ::   ln_rhg_EAP       ! EAP rheology switch, used for rdgrft and rheology
    LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F) 
    REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9
@@ -246,4 +248,6 @@
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   divu_i          !: Divergence of the velocity field             [s-1]
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   shear_i         !: Shear of the velocity field                  [s-1]
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   aniso_11, aniso_12   !: structure tensor elements
+   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rdg_conv
    !
    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   t_bo            !: Sea-Ice bottom temperature [Kelvin]     
@@ -436,5 +440,6 @@
       ALLOCATE( u_oce    (jpi,jpj) , v_oce    (jpi,jpj) , ht_i_new  (jpi,jpj) , strength(jpi,jpj) ,  &
          &      stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) ,                      &
-         &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) , STAT=ierr(ii) )
+         &      delta_i  (jpi,jpj) , divu_i   (jpi,jpj) , shear_i   (jpi,jpj) ,                      &
+         &      aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv  (jpi,jpj) , STAT=ierr(ii) )
 
       ii = ii + 1
Index: /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rdgrft.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rdgrft.F90	(revision 15569)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rdgrft.F90	(revision 15570)
@@ -138,4 +138,5 @@
       INTEGER , DIMENSION(jpij) ::   iptidx        ! compute ridge/raft or not
       REAL(wp), DIMENSION(jpij) ::   zdivu, zdelt  ! 1D divu_i & delta_i
+      REAL(wp), DIMENSION(jpij) ::   zconv         ! 1D rdg_conv (if EAP rheology)
       !
       INTEGER, PARAMETER ::   jp_itermax = 20    
@@ -174,5 +175,6 @@
         
          ! just needed here
-         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zdelt   (1:npti)      , delta_i  )
+         CALL tab_2d_1d( npti, nptidx(1:npti), zconv   (1:npti)      , rdg_conv )
          ! needed here and in the iteration loop
          CALL tab_2d_1d( npti, nptidx(1:npti), zdivu   (1:npti)      , divu_i) ! zdivu is used as a work array here (no change in divu_i)
@@ -184,5 +186,6 @@
             ! closing_net = rate at which open water area is removed + ice area removed by ridging 
             !                                                        - ice area added in new ridges
-            closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
+            IF( ln_rhg_EVP )  closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp )
+            IF( ln_rhg_EAP )  closing_net(ji) = zconv(ji)
             !
             IF( zdivu(ji) < 0._wp )   closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) )   ! make sure the closing rate is large enough
Index: /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg.F90	(revision 15569)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg.F90	(revision 15570)
@@ -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,7 @@
    !                                        ! 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
    !
    !! * Substitutions
@@ -81,8 +81,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, aniso_11, aniso_12, rdg_conv )
       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 +115,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, nn_rhg_chkcvg
+      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, ln_rhg_EAP, rn_creepl, rn_ecc , nn_nevp, rn_relast, nn_rhg_chkcvg
       !!-------------------------------------------------------------------
       !
@@ -137,4 +142,5 @@
          ELSEIF( nn_rhg_chkcvg == 2 ) THEN   ;   WRITE(numout,*) '         check cvg at both main and rheology time steps'
          ENDIF
+         WRITE(numout,*) '      rheology EAP (icedyn_rhg_eap)                        ln_rhg_EAP = ', ln_rhg_EAP
       ENDIF
       !
@@ -142,8 +148,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/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg_eap.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg_eap.F90	(revision 15570)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg_eap.F90	(revision 15570)
@@ -0,0 +1,2037 @@
+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]
+   !!                 !  2019     (S. Rynders, Y. Aksenov, C. Rousset)  change into eap rheology from
+   !!                                           CICE code (Tsamados, Heorton) 
+   !!----------------------------------------------------------------------
+#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
+
+   USE netcdf         ! NetCDF library for convergence test
+   IMPLICIT NONE
+   PRIVATE
+
+   PUBLIC   ice_dyn_rhg_eap   ! called by icedyn_rhg.F90
+   PUBLIC   rhg_eap_rst       ! called by icedyn_rhg.F90
+
+   REAL(wp), PARAMETER ::   pphi = 3.141592653589793_wp/12._wp    ! diamond shaped floe smaller angle (default phi = 30 deg)
+
+   ! look-up table for calculating structure tensor
+   INTEGER, PARAMETER ::   nx_yield = 41
+   INTEGER, PARAMETER ::   ny_yield = 41
+   INTEGER, PARAMETER ::   na_yield = 21
+
+   REAL(wp), DIMENSION(nx_yield, ny_yield, na_yield) ::   s11r, s12r, s22r, s11s, s12s, s22s
+
+   !! * Substitutions
+#  include "vectopt_loop_substitute.h90"
+
+   !! for convergence tests
+   INTEGER ::   ncvgid   ! netcdf file id
+   INTEGER ::   nvarid   ! netcdf variable id
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
+   !!----------------------------------------------------------------------
+   !! 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, paniso_11, paniso_12, prdg_conv )
+      !!-------------------------------------------------------------------
+      !!                 ***  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      !
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   paniso_11 , paniso_12                 ! structure tensor components
+      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   prdg_conv                             ! for ridging
+      !!
+      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) ::   zbetau, zbetav
+      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
+      REAL(wp) ::   zds2, zdt, zdt2, zdiv, zdiv2, zdsT                  ! 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) ::   zintb, zintn                                        ! dummy argument
+      REAL(wp) ::   zfac_x, zfac_y
+      REAL(wp) ::   zshear, zdum1, zdum2
+      REAL(wp) ::   zstressptmp, zstressmtmp, zstress12tmpF             ! anisotropic stress tensor components
+      REAL(wp) ::   zalphar, zalphas                                    ! for mechanical redistribution
+      REAL(wp) ::   zmresult11, zmresult12, z1dtevpkth, zp5kth, z1_dtevp_A  ! for structure tensor evolution
+      !
+      REAL(wp), DIMENSION(jpi,jpj) ::   zstress12tmp                    ! anisotropic stress tensor component for regridding
+      REAL(wp), DIMENSION(jpi,jpj) ::   zyield11, zyield22, zyield12    ! yield surface tensor for history
+      REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points
+      REAL(wp), DIMENSION(jpi,jpj) ::   zten_i                          ! tension
+      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) ::   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                          ! 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
+      !! --- check convergence
+      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice
+      !! --- diags
+      REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength
+      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p         
+      !! --- 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'
+      !
+      ! for diagnostics and convergence tests
+      ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) )
+      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
+            zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
+         END DO
+      END DO
+      !
+!!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)
+      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( umask(ji,jj,1), umask(ji,jj+1,1), &
+                  &                                          vmask(ji,jj,1), vmask(ji+1,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( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) )
+         ENDIF
+         IF( zfmask(jpi,jj) == 0._wp ) THEN
+            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) )
+         ENDIF
+      END DO
+      DO ji = 2, jpim1
+         IF( zfmask(ji,1) == 0._wp ) THEN
+            zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) )
+         ENDIF
+         IF( zfmask(ji,jpj) == 0._wp ) THEN
+            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) )
+         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
+
+      ! alpha parameters (Bouillon 2009)
+      IF( .NOT. ln_aEVP ) THEN
+         zdtevp   = rdt_ice / REAL( nn_nevp )
+         zalph1 =   2._wp * rn_relast * REAL( nn_nevp )
+         zalph2 = zalph1 * z1_ecc2
+
+         z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+      ELSE
+         zdtevp   = rdt_ice
+         ! zalpha parameters set later on adaptatively
+      ENDIF
+      z1_dtevp = 1._wp / zdtevp
+         
+      ! Initialise stress tensor 
+      zs1 (:,:) = pstress1_i (:,:) 
+      zs2 (:,:) = pstress2_i (:,:)
+      zs12(:,:) = pstress12_i(:,:)
+
+      ! constants for structure tensor
+      z1_dtevp_A = z1_dtevp/10.0_wp
+      z1dtevpkth = 1._wp / (z1_dtevp_A + 0.00002_wp)
+      zp5kth = 0.5_wp * 0.00002_wp
+
+      ! 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_lf_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_lf_depfra * hu_n(ji,jj)
+               ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) )
+               ! ice-bottom stress at V points
+               zvCr = zaV(ji,jj) * rn_lf_depfra * hv_n(ji,jj)
+               ztauy_base(ji,jj) = - rn_lf_bfr * 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_lf_depfra * ht_n(ji,jj)
+               tau_icebfr(ji,jj) = - rn_lf_bfr * 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
+         !
+         ! convergence test
+         IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2  ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step
+                  zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1)
+               END DO
+            END DO
+         ENDIF
+
+         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- !
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+
+               ! shear at F points
+               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   &
+                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   &
+                  &         ) * r1_e1e2f(ji,jj) * zfmask(ji,jj)
+
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_rhg_eap', zds, 'F', 1._wp )
+
+         DO jj = 2, jpjm1
+            DO ji = 2, jpim1 ! no vector loop
+
+               ! shear**2 at T points (doc eq. A16)
+               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+              
+               ! divergence at T points
+               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   &
+                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   &
+                  &    ) * r1_e1e2t(ji,jj)
+               zdiv2 = zdiv * zdiv
+               
+               ! tension at T points
+               zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   &
+                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   &
+                  &   ) * r1_e1e2t(ji,jj)
+               zdt2 = zdt * zdt
+
+               ! delta at T points
+               zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
+
+            END DO
+         END DO
+         CALL lbc_lnk( 'icedyn_rhg_eap', zdelta, 'T', 1._wp )
+               
+         ! P/delta at T points
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl )
+            END DO
+         END DO
+
+         DO jj = 2, jpj    ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
+            DO ji = 2, jpi ! no vector loop
+
+                ! shear at T points 
+               zdsT = ( zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+                  &   + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+               !WRITE(numout,*) 'shear output', ji, jj, zdsT
+               
+              ! divergence at T points (duplication to avoid communications)
+               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)
+               
+               ! tension at T points (duplication to avoid communications)
+               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)
+
+               ! --- anisotropic stress calculation --- !
+               CALL update_stress_rdg (jter, nn_nevp, zdiv, zdt, zdsT, paniso_11(ji,jj), paniso_12(ji,jj), &
+                                       zstressptmp, zstressmtmp, zstress12tmp(ji,jj), strength(ji,jj), zalphar, zalphas)
+
+               ! structure tensor update
+!!$               IF (mod(jter,10) == 0) THEN
+                  CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), paniso_11(ji,jj), paniso_12(ji,jj), zmresult11,  zmresult12)
+
+!!$                  paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A  + zp5kth + zmresult11) * z1dtevpkth ! implicit
+!!$                  paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A  + zmresult12) * z1dtevpkth ! implicit
+                  paniso_11(ji,jj) = (paniso_11(ji,jj)  + 0.5*2.e-5*zdtevp + zmresult11*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit
+                  paniso_12(ji,jj) = (paniso_12(ji,jj)                     + zmresult12*zdtevp) / (1. + 2.e-5*zdtevp) ! implicit
+!!$               ENDIF
+
+               IF (jter == nn_nevp) THEN
+                  zyield11(ji,jj) = 0.5_wp * (zstressptmp + zstressmtmp)
+                  zyield22(ji,jj) = 0.5_wp * (zstressptmp - zstressmtmp)
+                  zyield12(ji,jj) = zstress12tmp(ji,jj)
+                  prdg_conv(ji,jj) = -min( zalphar, 0._wp)    
+               ENDIF
+
+               ! alpha 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
+                  ! explicit:
+                  ! z1_alph1 = 1._wp / zalph1
+                  ! z1_alph2 = 1._wp / zalph1
+                  ! zalph1 = zalph1 - 1._wp
+                  ! zalph2 = zalph1
+               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(ji,jj)*(1._wp - zkt) ) ) * z1_alph1
+               !zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
+!!$               zs1(ji,jj) = ( zs1(ji,jj) + zstressptmp * zalph1  ) * z1_alph1
+!!$               zs2(ji,jj) = ( zs2(ji,jj) + zstressmtmp * zalph1  ) * z1_alph1
+               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zstressptmp ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj) * zalph1 + zstressmtmp ) * z1_alph1
+               !zs1(ji,jj) = ( stressptmp * zs1(ji,jj) + zalph1  ) * z1_alph1
+               !zs2(ji,jj) = ( stressmtmp * zs2(ji,jj) + zalph1  ) * z1_alph1
+               !WRITE(numout,*) 'stress output', ji, jj, zs1(ji,jj)
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zstress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.)
+
+        ! Save beta at T-points for further computations
+         IF( ln_aEVP ) THEN
+            DO jj = 1, jpj
+               DO ji = 1, jpi
+                  zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+               END DO
+            END DO
+         ENDIF
+
+         DO jj = 1, jpjm1
+            DO ji = 1, jpim1
+               ! stress12tmp at F points 
+               zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  &
+                  &            + zstress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + zstress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  &
+                  &            ) * 0.25_wp * r1_e1e2f(ji,jj)
+
+               ! alpha for aEVP
+               IF( ln_aEVP ) THEN
+                  zalph2   = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) )
+                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
+                  ! explicit:
+                  ! z1_alph2 = 1._wp / zalph2
+                  ! zalph2 = zalph2 - 1._wp
+               ENDIF
+               
+               ! 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
+!!$               zs12(ji,jj) = ( zs12(ji,jj) + zstress12tmpF * zalph1  ) * z1_alph1
+               zs12(ji,jj) = ( zs12(ji,jj) * zalph1 + zstress12tmpF ) * z1_alph1
+               !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1  ) * z1_alph1
+
+            END DO
+         END DO
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. )
+
+         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- !
+         DO jj = 2, jpjm1
+            DO ji = fs_2, fs_jpim1               
+               !                   !--- U points
+               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             &
+                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    &
+                  &                    ) * r1_e2u(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  &
+                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              &
+                  &                  ) * r1_e1e2u(ji,jj)
+               !
+               !                !--- V points
+               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             &
+                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    &
+                  &                    ) * r1_e1v(ji,jj)                                                                      &
+                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  &
+                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              &
+                  &                  ) * r1_e1e2v(ji,jj)
+               !
+               !                !--- 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)
+                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
+                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * 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) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   & 
+                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
+                        &                                    ) / ( zbetav + 1._wp )                                              &
+                        &             ) * 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_lf_relax )         & ! 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)
+                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
+                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * 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) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
+                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
+                        &                                    ) / ( zbetau + 1._wp )                                              &
+                        &             ) * 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_lf_relax )         & ! 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)
+                     zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) )
+                     u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * 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) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   &
+                        &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
+                        &                                    ) / ( zbetau + 1._wp )                                              &
+                        &             ) * 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_lf_relax )         & ! 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)
+                     zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) )
+                     v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * 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) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast
+                        &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &
+                        &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0
+                        &                                    ) / ( zbetav + 1._wp )                                              & 
+                        &             ) * 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_lf_relax )         & ! 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
+
+         ! convergence test
+         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg_eap( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
+         !
+         !                                                ! ==================== !
+      END DO                                              !  end loop over jter  !
+      !                                                   ! ==================== !
+      IF( ln_aEVP )   CALL iom_put( 'beta_evp' , zbeta )
+      !
+      CALL lbc_lnk( 'icedyn_rhg_eap', prdg_conv, 'T', 1. )  ! only need this in ridging module after rheology completed
+      !
+      !------------------------------------------------------------------------------!
+      ! 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
+            
+            zten_i(ji,jj) = 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
+            zfac            = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) ! delta  
+            rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zfac ) ) ! 0 if delta=0
+            pdelta_i(ji,jj) = zfac + rn_creepl * rswitch ! delta+creepl
+
+         END DO
+      END DO
+      CALL lbc_lnk_multi( 'icedyn_rhg_eap', pshear_i, 'T', 1., pdivu_i, 'T', 1., pdelta_i, 'T', 1., zten_i, 'T', 1., &
+         &                                  zs1     , 'T', 1., zs2    , 'T', 1., zs12    , 'F', 1. )
+      
+      ! --- Store the stress tensor for the next time step --- !
+      pstress1_i (:,:) = zs1 (:,:)
+      pstress2_i (:,:) = zs2 (:,:)
+      pstress12_i(:,:) = zs12(:,:)
+      !
+
+      !------------------------------------------------------------------------------!
+      ! 5) diagnostics
+      !------------------------------------------------------------------------------!
+      ! --- 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('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * zmsk00 )   ! delta
+      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * zmsk00 )   ! strength
+
+      ! --- Stress tensor invariants (SIMIP diags) --- !
+      IF( iom_use('normstr') .OR. iom_use('sheastr') ) THEN
+         !
+         ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
+         !         
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+            
+               ! Ice stresses
+               ! sigma1, sigma2, sigma12 are some useful recombination of the stresses (Hunke and Dukowicz MWR 2002, Bouillon et al., OM2013)
+               ! These are NOT stress tensor components, neither stress invariants, neither stress principal components
+               ! I know, this can be confusing...
+               zfac             =   strength(ji,jj) / ( pdelta_i(ji,jj) + rn_creepl ) 
+               zsig1            =   zfac * ( pdivu_i(ji,jj) - pdelta_i(ji,jj) )
+               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
+               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
+               
+               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008)
+               zsig_I (ji,jj)   =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure
+               zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress
+               
+            END DO
+         END DO         
+         !
+         ! Stress tensor invariants (normal and shear stress N/m) - SIMIP diags - definitions following Coon (1974) and Feltham (2008)
+         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * zmsk00(:,:) ) ! Normal stress
+         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * zmsk00(:,:) ) ! Maximum shear stress
+         
+         DEALLOCATE ( zsig_I, zsig_II )
+         
+      ENDIF
+
+      ! --- Normalized stress tensor principal components --- !
+      ! This are used to plot the normalized yield curve, see Lemieux & Dupont, 2020
+      ! Recommendation 1 : we use ice strength, not replacement pressure
+      ! Recommendation 2 : need to use deformations at PREVIOUS iterate for viscosities
+      IF( iom_use('sig1_pnorm') .OR. iom_use('sig2_pnorm') ) THEN
+         !
+         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )         
+         !         
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+            
+               ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 
+               !                        and **deformations** at current iterates
+               !                        following Lemieux & Dupont (2020)
+               zfac             =   zp_delt(ji,jj)
+               zsig1            =   zfac * ( pdivu_i(ji,jj) - ( zdelta(ji,jj) + rn_creepl ) )
+               zsig2            =   zfac * z1_ecc2 * zten_i(ji,jj)
+               zsig12           =   zfac * z1_ecc2 * pshear_i(ji,jj)
+               
+               ! Stress invariants (sigma_I, sigma_II, Coon 1974, Feltham 2008), T-point
+               zsig_I(ji,jj)    =   zsig1 * 0.5_wp                                           ! 1st stress invariant, aka average normal stress, aka negative pressure
+               zsig_II(ji,jj)   =   SQRT ( MAX( 0._wp, zsig2 * zsig2 * 0.25_wp + zsig12 ) )  ! 2nd  ''       '', aka maximum shear stress
+      
+               ! Normalized  principal stresses (used to display the ellipse)
+               z1_strength      =   1._wp / MAX( 1._wp, strength(ji,jj) )
+               zsig1_p(ji,jj)   =   ( zsig_I(ji,jj) + zsig_II(ji,jj) ) * z1_strength
+               zsig2_p(ji,jj)   =   ( zsig_I(ji,jj) - zsig_II(ji,jj) ) * z1_strength
+            END DO
+         END DO               
+         !
+         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
+         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
+      
+         DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
+         
+      ENDIF
+
+      ! --- yieldcurve --- !
+      IF( iom_use('yield11') .OR. iom_use('yield12') .OR. iom_use('yield22')) THEN
+
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zyield11, 'T', 1., zyield22, 'T', 1., zyield12, 'T', 1. )
+
+         CALL iom_put( 'yield11', zyield11 * zmsk00 )
+         CALL iom_put( 'yield22', zyield22 * zmsk00 )
+         CALL iom_put( 'yield12', zyield12 * zmsk00 )
+      ENDIF
+
+      ! --- anisotropy tensor --- !
+      IF( iom_use('aniso') ) THEN       
+         CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1. ) 
+         CALL iom_put( 'aniso' , paniso_11 * zmsk00 )
+      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
+      !
+      ! --- convergence tests --- !
+      IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN
+         IF( iom_use('uice_cvg') ) THEN
+            IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
+               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , &
+                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) )
+            ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) )
+               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , &
+                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) )
+            ENDIF
+         ENDIF
+      ENDIF      
+      !
+      DEALLOCATE( zmsk00, zmsk15 )
+      !
+   END SUBROUTINE ice_dyn_rhg_eap
+
+   
+   SUBROUTINE rhg_cvg_eap( kt, kiter, kitermax, pu, pv, pub, pvb )
+      !!----------------------------------------------------------------------
+      !!                    ***  ROUTINE rhg_cvg_eap  ***
+      !!                     
+      !! ** Purpose :   check convergence of oce rheology
+      !!
+      !! ** Method  :   create a file ice_cvg.nc containing the convergence of ice velocity
+      !!                during the sub timestepping of rheology so as:
+      !!                  uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) )
+      !!                This routine is called every sub-iteration, so it is cpu expensive
+      !!
+      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)   
+      !!----------------------------------------------------------------------
+      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index
+      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities
+      !!
+      INTEGER           ::   it, idtime, istatus
+      INTEGER           ::   ji, jj          ! dummy loop indices
+      REAL(wp)          ::   zresm           ! local real 
+      CHARACTER(len=20) ::   clname
+      REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence
+      !!----------------------------------------------------------------------
+
+      ! create file
+      IF( kt == nit000 .AND. kiter == 1 ) THEN
+         !
+         IF( lwp ) THEN
+            WRITE(numout,*)
+            WRITE(numout,*) 'rhg_cvg_eap : ice rheology convergence control'
+            WRITE(numout,*) '~~~~~~~'
+         ENDIF
+         !
+         IF( lwm ) THEN
+            clname = 'ice_cvg.nc'
+            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname)
+            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid )
+            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime )
+            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid )
+            istatus = NF90_ENDDEF(ncvgid)
+         ENDIF
+         !
+      ENDIF
+
+      ! time
+      it = ( kt - 1 ) * kitermax + kiter
+      
+      ! convergence
+      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large)
+         zresm = 0._wp
+      ELSE
+         DO jj = 1, jpj
+            DO ji = 1, jpi
+               zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), &
+                  &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj)
+            END DO
+         END DO
+         zresm = MAXVAL( zres )
+         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
+      ENDIF
+
+      IF( lwm ) THEN
+         ! write variables
+         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) )
+         ! close file
+         IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid)
+      ENDIF
+      
+   END SUBROUTINE rhg_cvg_eap
+
+
+   SUBROUTINE update_stress_rdg( ksub, kndte, pdivu, ptension, pshear, pa11, pa12, &
+      &                          pstressp,  pstressm, pstress12, pstrength, palphar, palphas )
+      !!---------------------------------------------------------------------
+      !!                   ***  SUBROUTINE update_stress_rdg  ***
+      !!                     
+      !! ** Purpose :   Updates the stress depending on values of strain rate and structure
+      !!                tensor and for the last subcycle step it computes closing and sliding rate
+      !!---------------------------------------------------------------------
+      INTEGER,  INTENT(in   ) ::   ksub, kndte
+      REAL(wp), INTENT(in   ) ::   pstrength
+      REAL(wp), INTENT(in   ) ::   pdivu, ptension, pshear
+      REAL(wp), INTENT(in   ) ::   pa11, pa12
+      REAL(wp), INTENT(  out) ::   pstressp, pstressm, pstress12
+      REAL(wp), INTENT(  out) ::   palphar, palphas
+
+      INTEGER  ::   kx ,ky, ka
+
+      REAL(wp) ::   zstemp11r, zstemp12r, zstemp22r   
+      REAL(wp) ::   zstemp11s, zstemp12s, zstemp22s 
+      REAL(wp) ::   za22, zQd11Qd11, zQd11Qd12, zQd12Qd12 
+      REAL(wp) ::   zQ11Q11, zQ11Q12, zQ12Q12 
+      REAL(wp) ::   zdtemp11, zdtemp12, zdtemp22 
+      REAL(wp) ::   zrotstemp11r, zrotstemp12r, zrotstemp22r   
+      REAL(wp) ::   zrotstemp11s, zrotstemp12s, zrotstemp22s
+      REAL(wp) ::   zsig11, zsig12, zsig22 
+      REAL(wp) ::   zsgprm11, zsgprm12, zsgprm22 
+      REAL(wp) ::   zinvstressconviso 
+      REAL(wp) ::   zAngle_denom_gamma,  zAngle_denom_alpha 
+      REAL(wp) ::   zTany_1, zTany_2 
+      REAL(wp) ::   zx, zy, zdx, zdy, zda, zkxw, kyw, kaw
+      REAL(wp) ::   zinvdx, zinvdy, zinvda   
+      REAL(wp) ::   zdtemp1, zdtemp2, zatempprime, zinvsin
+
+      REAL(wp), PARAMETER ::   kfriction = 0.45_wp
+      !!---------------------------------------------------------------------
+      ! Factor to maintain the same stress as in EVP (see Section 3)
+      ! Can be set to 1 otherwise
+!      zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction)
+      zinvstressconviso = 1._wp
+ 
+      zinvsin = 1._wp/sin(2._wp*pphi) * zinvstressconviso 
+      !now uses phi as set in higher code
+       
+      ! compute eigenvalues, eigenvectors and angles for structure tensor, strain
+      ! rates
+
+      ! 1) structure tensor
+      za22 = 1._wp-pa11
+      zQ11Q11 = 1._wp
+      zQ12Q12 = rsmall
+      zQ11Q12 = rsmall
+ 
+      ! gamma: angle between general coordiantes and principal axis of A 
+      ! here Tan2gamma = 2 a12 / (a11 - a22) 
+      IF((ABS(pa11 - za22) > rsmall).OR.(ABS(pa12) > rsmall)) THEN      
+         zAngle_denom_gamma = 1._wp/sqrt( ( pa11 - za22 )*( pa11 - za22) + &
+                              4._wp*pa12*pa12 )
+   
+         zQ11Q11 = 0.5_wp + ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma   !Cos^2 
+         zQ12Q12 = 0.5_wp - ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma   !Sin^2
+         zQ11Q12 = pa12*zAngle_denom_gamma                     !CosSin 
+      ENDIF
+ 
+      ! rotation Q*atemp*Q^T
+      zatempprime = zQ11Q11*pa11 + 2.0_wp*zQ11Q12*pa12 + zQ12Q12*za22 
+      
+      ! make first principal value the largest
+      zatempprime = max(zatempprime, 1.0_wp - zatempprime)
+ 
+      ! 2) strain rate
+      zdtemp11 = 0.5_wp*(pdivu + ptension)
+      zdtemp12 = pshear*0.5_wp
+      zdtemp22 = 0.5_wp*(pdivu - ptension)
+
+      ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22) 
+
+      zQd11Qd11 = 1.0_wp
+      zQd12Qd12 = rsmall
+      zQd11Qd12 = rsmall
+
+      IF((ABS( zdtemp11 - zdtemp22) > rsmall).OR. (ABS(zdtemp12) > rsmall)) THEN
+
+         zAngle_denom_alpha = 1.0_wp/sqrt( ( zdtemp11 - zdtemp22 )* &
+                             ( zdtemp11 - zdtemp22 ) + 4.0_wp*zdtemp12*zdtemp12)
+
+         zQd11Qd11 = 0.5_wp + ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Cos^2 
+         zQd12Qd12 = 0.5_wp - ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Sin^2
+         zQd11Qd12 = zdtemp12*zAngle_denom_alpha !CosSin
+      ENDIF
+
+      zdtemp1 = zQd11Qd11*zdtemp11 + 2.0_wp*zQd11Qd12*zdtemp12 + zQd12Qd12*zdtemp22
+      zdtemp2 = zQd12Qd12*zdtemp11 - 2.0_wp*zQd11Qd12*zdtemp12 + zQd11Qd11*zdtemp22
+      ! In cos and sin values
+      zx = 0._wp
+      IF ((ABS(zdtemp1) > rsmall).OR.(ABS(zdtemp2) > rsmall)) THEN
+         zx = atan2(zdtemp2,zdtemp1)
+      ENDIF      
+               
+      ! to ensure the angle lies between pi/4 and 9 pi/4 
+      IF (zx < rpi*0.25_wp) zx = zx + rpi*2.0_wp
+      
+      ! y: angle between major principal axis of strain rate and structure
+      ! tensor
+      ! y = gamma - alpha 
+      ! Expressesed componently with
+      ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha)
+      
+      zTany_1 = zQ11Q12 - zQd11Qd12
+      zTany_2 = zQ11Q11 - zQd12Qd12
+      
+      zy = 0._wp
+      
+      IF ((ABS(zTany_1) > rsmall).OR.(ABS(zTany_2) > rsmall)) THEN
+         zy = atan2(zTany_1,zTany_2)
+      ENDIF
+      
+      ! to make sure y is between 0 and pi
+      IF (zy > rpi) zy = zy - rpi
+      IF (zy < 0)  zy = zy + rpi
+
+      ! 3) update anisotropic stress tensor 
+      zdx   = rpi/real(nx_yield-1,kind=wp)
+      zdy   = rpi/real(ny_yield-1,kind=wp)
+      zda   = 0.5_wp/real(na_yield-1,kind=wp)
+      zinvdx = 1._wp/zdx
+      zinvdy = 1._wp/zdy
+      zinvda = 1._wp/zda
+
+      ! % need 8 coords and 8 weights
+      ! % range in kx
+      kx  = int((zx-rpi*0.25_wp-rpi)*zinvdx) + 1
+      !!clem kx  = MAX( 1, MIN( nx_yield-1, INT((zx-rpi*0.25_wp-rpi)*zinvdx) + 1  ) )
+      zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx  
+      
+      ky  = int(zy*zinvdy) + 1
+      !!clem ky  = MAX( 1, MIN( ny_yield-1, INT(zy*zinvdy) + 1 ) )
+      kyw = ky - zy*zinvdy 
+      
+      ka  = int((zatempprime-0.5_wp)*zinvda) + 1
+      !!clem ka  = MAX( 1, MIN( na_yield-1, INT((zatempprime-0.5_wp)*zinvda) + 1 ) )
+      kaw = ka - (zatempprime-0.5_wp)*zinvda
+      
+      ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013)
+!!$      zstemp11r =  zkxw  * kyw         * kaw         * s11r(kx  ,ky  ,ka  ) &
+!!$        & + (1._wp-zkxw) * kyw         * kaw         * s11r(kx+1,ky  ,ka  ) &
+!!$        & + zkxw         * (1._wp-kyw) * kaw         * s11r(kx  ,ky+1,ka  ) &
+!!$        & + zkxw         * kyw         * (1._wp-kaw) * s11r(kx  ,ky  ,ka+1) &
+!!$        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s11r(kx+1,ky+1,ka  ) &
+!!$        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s11r(kx+1,ky  ,ka+1) &
+!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s11r(kx  ,ky+1,ka+1) &
+!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1)
+!!$      zstemp12r =  zkxw  * kyw         * kaw         * s12r(kx  ,ky  ,ka  ) &
+!!$        & + (1._wp-zkxw) * kyw         * kaw         * s12r(kx+1,ky  ,ka  ) &
+!!$        & + zkxw         * (1._wp-kyw) * kaw         * s12r(kx  ,ky+1,ka  ) &
+!!$        & + zkxw         * kyw         * (1._wp-kaw) * s12r(kx  ,ky  ,ka+1) &
+!!$        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s12r(kx+1,ky+1,ka  ) &
+!!$        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s12r(kx+1,ky  ,ka+1) &
+!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s12r(kx  ,ky+1,ka+1) &
+!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1)
+!!$      zstemp22r =  zkxw  * kyw         * kaw         * s22r(kx  ,ky  ,ka  ) &
+!!$        & + (1._wp-zkxw) * kyw         * kaw         * s22r(kx+1,ky  ,ka  ) &
+!!$        & + zkxw         * (1._wp-kyw) * kaw         * s22r(kx  ,ky+1,ka  ) &
+!!$        & + zkxw         * kyw         * (1._wp-kaw) * s22r(kx  ,ky  ,ka+1) &
+!!$        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s22r(kx+1,ky+1,ka  ) &
+!!$        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s22r(kx+1,ky  ,ka+1) &
+!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s22r(kx  ,ky+1,ka+1) &
+!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1)
+!!$      
+!!$      zstemp11s =  zkxw  * kyw         * kaw         * s11s(kx  ,ky  ,ka  ) &
+!!$        & + (1._wp-zkxw) * kyw         * kaw         * s11s(kx+1,ky  ,ka  ) &
+!!$        & + zkxw         * (1._wp-kyw) * kaw         * s11s(kx  ,ky+1,ka  ) &
+!!$        & + zkxw         * kyw         * (1._wp-kaw) * s11s(kx  ,ky  ,ka+1) &
+!!$        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s11s(kx+1,ky+1,ka  ) &
+!!$        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s11s(kx+1,ky  ,ka+1) &
+!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s11s(kx  ,ky+1,ka+1) &
+!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1)
+!!$      zstemp12s =  zkxw  * kyw         * kaw         * s12s(kx  ,ky  ,ka  ) &
+!!$        & + (1._wp-zkxw) * kyw         * kaw         * s12s(kx+1,ky  ,ka  ) &
+!!$        & + zkxw         * (1._wp-kyw) * kaw         * s12s(kx  ,ky+1,ka  ) &
+!!$        & + zkxw         * kyw         * (1._wp-kaw) * s12s(kx  ,ky  ,ka+1) &
+!!$        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s12s(kx+1,ky+1,ka  ) &
+!!$        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s12s(kx+1,ky  ,ka+1) &
+!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s12s(kx  ,ky+1,ka+1) &
+!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1)
+!!$      zstemp22s =  zkxw  * kyw         * kaw         * s22s(kx  ,ky  ,ka  ) &
+!!$        & + (1._wp-zkxw) * kyw         * kaw         * s22s(kx+1,ky  ,ka  ) &
+!!$        & + zkxw         * (1._wp-kyw) * kaw         * s22s(kx  ,ky+1,ka  ) &
+!!$        & + zkxw         * kyw         * (1._wp-kaw) * s22s(kx  ,ky  ,ka+1) &
+!!$        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s22s(kx+1,ky+1,ka  ) &
+!!$        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s22s(kx+1,ky  ,ka+1) &
+!!$        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s22s(kx  ,ky+1,ka+1) &
+!!$        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1)
+
+      zstemp11r = s11r(kx,ky,ka)
+      zstemp12r = s12r(kx,ky,ka)
+      zstemp22r = s22r(kx,ky,ka)
+      zstemp11s = s11s(kx,ky,ka)
+      zstemp12s = s12s(kx,ky,ka)
+      zstemp22s = s22s(kx,ky,ka)
+      
+      
+      ! Calculate mean ice stress over a collection of floes (Equation 3 in
+      ! Tsamados 2013)
+
+      zsig11  = pstrength*(zstemp11r + kfriction*zstemp11s) * zinvsin
+      zsig12  = pstrength*(zstemp12r + kfriction*zstemp12s) * zinvsin
+      zsig22  = pstrength*(zstemp22r + kfriction*zstemp22s) * zinvsin
+
+      !WRITE(numout,*) 'principal axis stress inside loop', sig11, sig12, sig22
+
+      ! Back - rotation of the stress from principal axes into general coordinates
+
+      ! Update stress
+      zsgprm11 = zQ11Q11*zsig11 + zQ12Q12*zsig22 -      2._wp*zQ11Q12 *zsig12
+      zsgprm12 = zQ11Q12*zsig11 - zQ11Q12*zsig22 + (zQ11Q11 - zQ12Q12)*zsig12
+      zsgprm22 = zQ12Q12*zsig11 + zQ11Q11*zsig22 +      2._wp*zQ11Q12 *zsig12
+
+      pstressp  = zsgprm11 + zsgprm22
+      pstress12 = zsgprm12
+      pstressm  = zsgprm11 - zsgprm22
+
+      ! Compute ridging and sliding functions in general coordinates 
+      ! (Equation 11 in Tsamados 2013)
+      IF (ksub == kndte) THEN
+         zrotstemp11r = zQ11Q11*zstemp11r - 2._wp*zQ11Q12* zstemp12r &
+                      + zQ12Q12*zstemp22r
+         zrotstemp12r = zQ11Q11*zstemp12r +       zQ11Q12*(zstemp11r-zstemp22r) &
+                      - zQ12Q12*zstemp12r
+         zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r & 
+                      + zQ11Q11*zstemp22r
+
+         zrotstemp11s = zQ11Q11*zstemp11s - 2._wp*zQ11Q12* zstemp12s &
+                      + zQ12Q12*zstemp22s
+         zrotstemp12s = zQ11Q11*zstemp12s +       zQ11Q12*(zstemp11s-zstemp22s) &
+                      - zQ12Q12*zstemp12s
+         zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s & 
+                      + zQ11Q11*zstemp22s
+
+         palphar =  zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 &
+                  + zrotstemp22r*zdtemp22
+         palphas =  zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 &
+                  + zrotstemp22s*zdtemp22
+      ENDIF
+   END SUBROUTINE update_stress_rdg 
+
+!=======================================================================
+
+
+   SUBROUTINE calc_ffrac( pstressp, pstressm, pstress12, pa11, pa12, &
+      &                   pmresult11, pmresult12 )
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE calc_ffrac  ***
+      !!                     
+      !! ** Purpose :   Computes term in evolution equation for structure tensor
+      !!                which determines the ice floe re-orientation due to fracture
+      !! ** Method :    Eq. 7: Ffrac = -kf(A-S) or = 0 depending on sigma_1 and sigma_2
+      !!---------------------------------------------------------------------
+      REAL (wp), INTENT(in)  ::   pstressp, pstressm, pstress12, pa11, pa12
+      REAL (wp), INTENT(out) ::   pmresult11, pmresult12
+
+      ! local variables
+
+      REAL (wp) ::   zsigma11, zsigma12, zsigma22  ! stress tensor elements
+      REAL (wp) ::   zAngle_denom        ! angle with principal component axis
+      REAL (wp) ::   zsigma_1, zsigma_2   ! principal components of stress
+      REAL (wp) ::   zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12
+
+!!$      REAL (wp), PARAMETER ::   kfrac = 0.0001_wp   ! rate of fracture formation 
+      REAL (wp), PARAMETER ::   kfrac = 1.e-3_wp   ! rate of fracture formation 
+      REAL (wp), PARAMETER ::   threshold = 0.3_wp  ! critical confinement ratio 
+      !!---------------------------------------------------------------
+      !
+      zsigma11 = 0.5_wp*(pstressp+pstressm) 
+      zsigma12 = pstress12 
+      zsigma22 = 0.5_wp*(pstressp-pstressm) 
+
+      ! Here's the change - no longer calculate gamma,
+      ! use trig formulation, angle signs are kept correct, don't worry
+   
+      ! rotate tensor to get into sigma principal axis
+   
+      ! here Tan2gamma = 2 sig12 / (sig11 - sig12) 
+      ! add rsmall to the denominator to stop 1/0 errors, makes very little 
+      ! error to the calculated angles < rsmall
+
+      zQ11Q11 = 0.1_wp
+      zQ12Q12 = rsmall
+      zQ11Q12 = rsmall
+
+      IF((ABS( zsigma11 - zsigma22) > rsmall).OR.(ABS(zsigma12) > rsmall)) THEN
+
+         zAngle_denom = 1.0_wp/sqrt( ( zsigma11 - zsigma22 )*( zsigma11 - &
+                       zsigma22 ) + 4.0_wp*zsigma12*zsigma12)
+
+         zQ11Q11 = 0.5_wp + ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Cos^2 
+         zQ12Q12 = 0.5_wp - ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Sin^2
+         zQ11Q12 = zsigma12*zAngle_denom                      !CosSin 
+      ENDIF
+
+      zsigma_1 = zQ11Q11*zsigma11 + 2.0_wp*zQ11Q12*zsigma12 + zQ12Q12*zsigma22 ! S(1,1)
+      zsigma_2 = zQ12Q12*zsigma11 - 2.0_wp*zQ11Q12*zsigma12 + zQ11Q11*zsigma22 ! S(2,2)
+
+      ! Pure divergence
+      IF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 >= 0.0_wp))  THEN
+         pmresult11 = 0.0_wp
+         pmresult12 = 0.0_wp
+
+      ! Unconfined compression: cracking of blocks not along the axial splitting
+      ! direction
+      ! which leads to the loss of their shape, so we again model it through diffusion
+      ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp))  THEN
+         pmresult11 = - kfrac * (pa11 - zQ12Q12)
+         pmresult12 = - kfrac * (pa12 + zQ11Q12)
+
+      ! Shear faulting
+      ELSEIF (zsigma_2 == 0.0_wp) THEN
+         pmresult11 = 0.0_wp
+         pmresult12 = 0.0_wp
+      ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN
+         pmresult11 = - kfrac * (pa11 - zQ12Q12)
+         pmresult12 = - kfrac * (pa12 + zQ11Q12)
+
+      ! Horizontal spalling
+      ELSE 
+         pmresult11 = 0.0_wp
+         pmresult12 = 0.0_wp
+      ENDIF
+
+   END SUBROUTINE calc_ffrac
+
+
+   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, id4, id5   ! local integers
+      INTEGER  ::   ix, iy, ip, iz, n, ia     ! local integers
+    
+      INTEGER, PARAMETER            ::    nz = 100
+      
+      REAL(wp) ::   ainit, xinit, yinit, pinit, zinit
+      REAL(wp) ::   da, dx, dy, dp, dz, a1
+
+      !!clem
+      REAL(wp) ::   zw1, zw2, zfac, ztemp
+      REAL(wp) ::   idx, idy, idz
+
+      REAL(wp), PARAMETER           ::   eps6 = 1.0e-6_wp
+      !!----------------------------------------------------------------------
+      !
+      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. )
+            id4 = iom_varid( numrir, 'aniso_11'  , ldstop = .FALSE. )
+            id5 = iom_varid( numrir, 'aniso_12'  , ldstop = .FALSE. )
+            !
+            IF( MIN( id1, id2, id3, id4, id5 ) > 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 )
+               CALL iom_get( numrir, jpdom_autoglo, 'aniso_11'  , aniso_11   )
+               CALL iom_get( numrir, jpdom_autoglo, 'aniso_12'  , aniso_12   )
+            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
+               aniso_11  (:,:) = 0.5_wp
+               aniso_12  (:,:) = 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
+            aniso_11  (:,:) = 0.5_wp
+            aniso_12  (:,:) = 0._wp
+         ENDIF
+         !
+
+         da = 0.5_wp/real(na_yield-1,kind=wp)
+         ainit = 0.5_wp - da
+         dx = rpi/real(nx_yield-1,kind=wp)
+         xinit = rpi + 0.25_wp*rpi - dx
+         dz = rpi/real(nz,kind=wp)
+         zinit = -rpi*0.5_wp
+         dy = rpi/real(ny_yield-1,kind=wp)
+         yinit = -dy
+
+         s11r(:,:,:) = 0._wp
+         s12r(:,:,:) = 0._wp
+         s22r(:,:,:) = 0._wp
+         s11s(:,:,:) = 0._wp
+         s12s(:,:,:) = 0._wp
+         s22s(:,:,:) = 0._wp
+
+!!$         DO ia=1,na_yield
+!!$            DO ix=1,nx_yield
+!!$               DO iy=1,ny_yield
+!!$                  s11r(ix,iy,ia) = 0._wp
+!!$                  s12r(ix,iy,ia) = 0._wp
+!!$                  s22r(ix,iy,ia) = 0._wp
+!!$                  s11s(ix,iy,ia) = 0._wp
+!!$                  s12s(ix,iy,ia) = 0._wp
+!!$                  s22s(ix,iy,ia) = 0._wp
+!!$                  IF (ia <= na_yield-1) THEN
+!!$                     DO iz=1,nz
+!!$                        s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+!!$                           s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)
+!!$                        s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+!!$                           s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)
+!!$                        s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+!!$                           s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 
+!!$                        s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+!!$                           s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)
+!!$                        s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+!!$                           s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)
+!!$                        s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+!!$                           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+!!$                           s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi)
+!!$                     ENDDO
+!!$                     IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp
+!!$                     IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp
+!!$                     IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp
+!!$                     IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp
+!!$                     IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp
+!!$                     IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp
+!!$                  ELSE
+!!$                     s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi)
+!!$                     s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi)
+!!$                     s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi)
+!!$                     s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi)
+!!$                     s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi)
+!!$                     s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi)
+!!$                     IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp
+!!$                     IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp
+!!$                     IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp
+!!$                     IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp
+!!$                     IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp
+!!$                     IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp
+!!$                  ENDIF
+!!$               ENDDO
+!!$            ENDDO
+!!$         ENDDO
+
+         !! faster but still very slow => to be improved         
+         zfac = dz/sin(2._wp*pphi)
+         DO ia = 1, na_yield-1
+            zw1 = w1(ainit+ia*da)
+            zw2 = w2(ainit+ia*da)
+            DO iz = 1, nz
+               idz = zinit+iz*dz
+               ztemp = zw1 * EXP(-zw2*(zinit+iz*dz)*(zinit+iz*dz))
+               DO iy = 1, ny_yield
+                  idy = yinit+iy*dy
+                  DO ix = 1, nx_yield
+                     idx = xinit+ix*dx
+                     s11r(ix,iy,ia) = s11r(ix,iy,ia) + ztemp * s11kr(idx,idy,idz)*zfac
+                     s12r(ix,iy,ia) = s12r(ix,iy,ia) + ztemp * s12kr(idx,idy,idz)*zfac
+                     s22r(ix,iy,ia) = s22r(ix,iy,ia) + ztemp * s22kr(idx,idy,idz)*zfac 
+                     s11s(ix,iy,ia) = s11s(ix,iy,ia) + ztemp * s11ks(idx,idy,idz)*zfac
+                     s12s(ix,iy,ia) = s12s(ix,iy,ia) + ztemp * s12ks(idx,idy,idz)*zfac
+                     s22s(ix,iy,ia) = s22s(ix,iy,ia) + ztemp * s22ks(idx,idy,idz)*zfac
+                  END DO
+               END DO
+            END DO
+         END DO
+
+         zfac = 1._wp/sin(2._wp*pphi)
+         ia = na_yield
+         DO iy = 1, ny_yield
+            idy = yinit+iy*dy
+            DO ix = 1, nx_yield
+               idx = xinit+ix*dx                  
+               s11r(ix,iy,ia) = 0.5_wp*s11kr(idx,idy,0._wp)*zfac
+               s12r(ix,iy,ia) = 0.5_wp*s12kr(idx,idy,0._wp)*zfac
+               s22r(ix,iy,ia) = 0.5_wp*s22kr(idx,idy,0._wp)*zfac
+               s11s(ix,iy,ia) = 0.5_wp*s11ks(idx,idy,0._wp)*zfac
+               s12s(ix,iy,ia) = 0.5_wp*s12ks(idx,idy,0._wp)*zfac
+               s22s(ix,iy,ia) = 0.5_wp*s22ks(idx,idy,0._wp)*zfac
+            ENDDO
+         ENDDO
+         WHERE (ABS(s11r(:,:,:)) < eps6) s11r(:,:,:) = 0._wp
+         WHERE (ABS(s12r(:,:,:)) < eps6) s12r(:,:,:) = 0._wp
+         WHERE (ABS(s22r(:,:,:)) < eps6) s22r(:,:,:) = 0._wp
+         WHERE (ABS(s11s(:,:,:)) < eps6) s11s(:,:,:) = 0._wp
+         WHERE (ABS(s12s(:,:,:)) < eps6) s12s(:,:,:) = 0._wp
+         WHERE (ABS(s22s(:,:,:)) < eps6) s22s(:,:,:) = 0._wp
+
+
+      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 )
+         CALL iom_rstput( iter, nitrst, numriw, 'aniso_11'  , aniso_11   )
+         CALL iom_rstput( iter, nitrst, numriw, 'aniso_12'  , aniso_12   )
+         !
+      ENDIF
+      !
+   END SUBROUTINE rhg_eap_rst
+
+
+   FUNCTION w1(pa)
+      !!-------------------------------------------------------------------
+      !! Function : w1 (see Gaussian function psi in Tsamados et al 2013)
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   pa
+      REAL(wp) ::   w1
+      !!-------------------------------------------------------------------
+   
+      w1 = -   223.87569446_wp &
+       &   +  2361.21986630_wp*pa &
+       &   - 10606.56079975_wp*pa*pa &
+       &   + 26315.50025642_wp*pa*pa*pa &
+       &   - 38948.30444297_wp*pa*pa*pa*pa &
+       &   + 34397.72407466_wp*pa*pa*pa*pa*pa &
+       &   - 16789.98003081_wp*pa*pa*pa*pa*pa*pa &
+       &   +  3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa
+   
+   END FUNCTION w1
+
+
+   FUNCTION w2(pa)
+      !!-------------------------------------------------------------------
+      !! Function : w2 (see Gaussian function psi in Tsamados et al 2013)
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   pa
+      REAL(wp) ::   w2
+      !!-------------------------------------------------------------------
+   
+      w2 = -    6670.68911883_wp &
+       &   +   70222.33061536_wp*pa &
+       &   -  314871.71525448_wp*pa*pa &
+       &   +  779570.02793492_wp*pa*pa*pa &
+       &   - 1151098.82436864_wp*pa*pa*pa*pa &
+       &   + 1013896.59464498_wp*pa*pa*pa*pa*pa &
+       &   -  493379.44906738_wp*pa*pa*pa*pa*pa*pa &
+       &   +  102356.55151800_wp*pa*pa*pa*pa*pa*pa*pa
+   
+   END FUNCTION w2
+
+   FUNCTION s11kr(px,py,pz) 
+      !!-------------------------------------------------------------------
+      !! Function : s11kr
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   px,py,pz
+   
+      REAL(wp) ::   s11kr, zpih
+   
+      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22
+      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22
+      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22
+      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22
+      REAL(wp) ::   zd11, zd12, zd22
+      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2
+      REAL(wp) ::   zHen1t2, zHen2t1
+      !!-------------------------------------------------------------------
+   
+      zpih = 0.5_wp*rpi
+   
+      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi)
+      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi)
+      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi)
+      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi)
+      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi)
+      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi)
+      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi)
+      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi)
+      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi)
+      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi)
+      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi)
+      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi)
+      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi)
+      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi)
+      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi)
+      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi)
+   ! In expression of tensor d, with this formulatin d(x)=-d(x+pi)
+   ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else
+   ! x=x-pi/2
+      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py))
+      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px))
+      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py))
+      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22
+      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22
+      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22
+   
+      IF (-zIIn1t2>=rsmall) THEN
+      zHen1t2 = 1._wp
+      ELSE
+      zHen1t2 = 0._wp
+      ENDIF
+   
+      IF (-zIIn2t1>=rsmall) THEN
+      zHen2t1 = 1._wp
+      ELSE
+      zHen2t1 = 0._wp
+      ENDIF
+   
+      s11kr = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11)
+
+   END FUNCTION s11kr
+
+   FUNCTION s12kr(px,py,pz)
+      !!-------------------------------------------------------------------
+      !! Function : s12kr
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   px,py,pz
+
+      REAL(wp) ::   s12kr, zs12r0, zs21r0, zpih
+
+      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22
+      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22
+      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22
+      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22
+      REAL(wp) ::   zd11, zd12, zd22
+      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2
+      REAL(wp) ::   zHen1t2, zHen2t1
+      !!-------------------------------------------------------------------
+      zpih = 0.5_wp*rpi
+   
+      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi)
+      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi)
+      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi)
+      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi)
+      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi)
+      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi)
+      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi)
+      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi)
+      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi)
+      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi)
+      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi)
+      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi)
+      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi)
+      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi)
+      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi)
+      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi)
+      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py))
+      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px))
+      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py))
+      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22
+      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22
+      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22
+   
+      IF (-zIIn1t2>=rsmall) THEN
+      zHen1t2 = 1._wp
+      ELSE
+      zHen1t2 = 0._wp
+      ENDIF
+   
+      IF (-zIIn2t1>=rsmall) THEN
+      zHen2t1 = 1._wp
+      ELSE
+      zHen2t1 = 0._wp
+      ENDIF
+   
+      zs12r0 = (- zHen1t2 * zn1t2i12 - zHen2t1 * zn2t1i12)
+      zs21r0 = (- zHen1t2 * zn1t2i21 - zHen2t1 * zn2t1i21)
+      s12kr=0.5_wp*(zs12r0+zs21r0)
+   
+   END FUNCTION s12kr
+
+   FUNCTION s22kr(px,py,pz)
+      !!-------------------------------------------------------------------
+      !! Function : s22kr
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   px,py,pz
+
+      REAL(wp) ::   s22kr, zpih
+
+      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22
+      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22
+      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22
+      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22
+      REAL(wp) ::   zd11, zd12, zd22
+      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2
+      REAL(wp) ::   zHen1t2, zHen2t1
+      !!-------------------------------------------------------------------
+      zpih = 0.5_wp*rpi
+
+      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi)
+      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi)
+      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi)
+      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi)
+      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi)
+      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi)
+      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi)
+      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi)
+      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi)
+      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi)
+      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi)
+      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi)
+      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi)
+      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi)
+      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi)
+      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi)
+      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py))
+      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px))
+      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py))
+      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22
+      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22
+      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22
+
+      IF (-zIIn1t2>=rsmall) THEN
+      zHen1t2 = 1._wp
+      ELSE
+      zHen1t2 = 0._wp
+      ENDIF
+
+      IF (-zIIn2t1>=rsmall) THEN
+      zHen2t1 = 1._wp
+      ELSE
+      zHen2t1 = 0._wp
+      ENDIF
+
+      s22kr = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22)
+
+   END FUNCTION s22kr
+
+   FUNCTION s11ks(px,py,pz)
+      !!-------------------------------------------------------------------
+      !! Function : s11ks
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   px,py,pz
+
+      REAL(wp) ::   s11ks, zpih
+
+      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22
+      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22
+      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22
+      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22
+      REAL(wp) ::   zd11, zd12, zd22
+      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2
+      REAL(wp) ::   zHen1t2, zHen2t1
+      !!-------------------------------------------------------------------
+      zpih = 0.5_wp*rpi
+
+      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi)
+      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi)
+      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi)
+      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi)
+      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi)
+      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi)
+      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi)
+      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi)
+      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi)
+      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi)
+      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi)
+      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi)
+      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi)
+      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi)
+      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi)
+      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi)
+      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py))
+      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px))
+      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py))
+      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22
+      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22
+      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22
+
+      IF (-zIIn1t2>=rsmall) THEN
+      zHen1t2 = 1._wp
+      ELSE
+      zHen1t2 = 0._wp
+      ENDIF
+
+      IF (-zIIn2t1>=rsmall) THEN
+      zHen2t1 = 1._wp
+      ELSE
+      zHen2t1 = 0._wp
+      ENDIF
+
+      s11ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11)
+
+   END FUNCTION s11ks
+
+   FUNCTION s12ks(px,py,pz)
+      !!-------------------------------------------------------------------
+      !! Function : s12ks
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   px,py,pz
+
+      REAL(wp) ::   s12ks, zs12s0, zs21s0, zpih
+
+      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22
+      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22
+      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22
+      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22
+      REAL(wp) ::   zd11, zd12, zd22
+      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2
+      REAL(wp) ::   zHen1t2, zHen2t1
+      !!-------------------------------------------------------------------
+      zpih = 0.5_wp*rpi
+
+      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi)
+      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi)
+      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi)
+      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi)
+      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi)
+      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi)
+      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi)
+      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi)
+      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi)
+      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi)
+      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi)
+      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi)
+      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi)
+      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi)
+      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi)
+      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi)
+      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py))
+      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px))
+      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py))
+      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22
+      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22
+      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22
+
+      IF (-zIIn1t2>=rsmall) THEN
+      zHen1t2 = 1._wp
+      ELSE
+      zHen1t2 = 0._wp
+      ENDIF
+
+      IF (-zIIn2t1>=rsmall) THEN
+      zHen2t1 = 1._wp
+      ELSE
+      zHen2t1 = 0._wp
+      ENDIF
+
+      zs12s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i12 + zHen2t1 * zt2t1i12)
+      zs21s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i21 + zHen2t1 * zt2t1i21)
+      s12ks=0.5_wp*(zs12s0+zs21s0)
+
+   END FUNCTION s12ks
+
+   FUNCTION s22ks(px,py,pz) 
+      !!-------------------------------------------------------------------
+      !! Function : s22ks
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   px,py,pz
+
+      REAL(wp) ::   s22ks, zpih
+
+      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22
+      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22
+      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22
+      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22
+      REAL(wp) ::   zd11, zd12, zd22
+      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2
+      REAL(wp) ::   zHen1t2, zHen2t1
+      !!-------------------------------------------------------------------
+      zpih = 0.5_wp*rpi
+
+      zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi)
+      zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi)
+      zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi)
+      zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi)
+      zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi)
+      zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi)
+      zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi)
+      zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi)
+      zt1t2i11 = cos(pz-pphi) * cos(pz+pphi)
+      zt1t2i12 = cos(pz-pphi) * sin(pz+pphi)
+      zt1t2i21 = sin(pz-pphi) * cos(pz+pphi)
+      zt1t2i22 = sin(pz-pphi) * sin(pz+pphi)
+      zt2t1i11 = cos(pz+pphi) * cos(pz-pphi)
+      zt2t1i12 = cos(pz+pphi) * sin(pz-pphi)
+      zt2t1i21 = sin(pz+pphi) * cos(pz-pphi)
+      zt2t1i22 = sin(pz+pphi) * sin(pz-pphi)
+      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py))
+      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px))
+      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py))
+      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22
+      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22
+      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22
+
+      IF (-zIIn1t2>=rsmall) THEN
+      zHen1t2 = 1._wp
+      ELSE
+      zHen1t2 = 0._wp
+      ENDIF
+
+      IF (-zIIn2t1>=rsmall) THEN
+      zHen2t1 = 1._wp
+      ELSE
+      zHen2t1 = 0._wp
+      ENDIF
+
+      s22ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22)
+
+   END FUNCTION s22ks
+
+#else
+   !!----------------------------------------------------------------------
+   !!   Default option         Empty module           NO SI3 sea-ice model
+   !!----------------------------------------------------------------------
+#endif
+
+   !!==============================================================================
+END MODULE icedyn_rhg_eap
Index: /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg_evp.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg_evp.F90	(revision 15569)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icedyn_rhg_evp.F90	(revision 15570)
@@ -736,5 +736,5 @@
 
          ! convergence test
-         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
+         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg_evp( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice )
          !
          !                                                ! ==================== !
@@ -962,7 +962,7 @@
 
 
-   SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb )
+   SUBROUTINE rhg_cvg_evp( kt, kiter, kitermax, pu, pv, pub, pvb )
       !!----------------------------------------------------------------------
-      !!                    ***  ROUTINE rhg_cvg  ***
+      !!                    ***  ROUTINE rhg_cvg_evp  ***
       !!                     
       !! ** Purpose :   check convergence of oce rheology
@@ -990,5 +990,5 @@
          IF( lwp ) THEN
             WRITE(numout,*)
-            WRITE(numout,*) 'rhg_cvg : ice rheology convergence control'
+            WRITE(numout,*) 'rhg_cvg_evp : ice rheology convergence control'
             WRITE(numout,*) '~~~~~~~'
          ENDIF
@@ -1029,5 +1029,5 @@
       ENDIF
       
-   END SUBROUTINE rhg_cvg
+   END SUBROUTINE rhg_cvg_evp
 
 
Index: /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icewri.F90
===================================================================
--- /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icewri.F90	(revision 15569)
+++ /NEMO/branches/UKMO/NEMO_4.0.4_EAP_rheology_fix/src/ICE/icewri.F90	(revision 15570)
@@ -254,4 +254,5 @@
       CALL iom_rstput( 0, 0, kid, 'sivolu', vt_i         )   ! Ice volume
       CALL iom_rstput( 0, 0, kid, 'sidive', divu_i*1.0e8 )   ! Ice divergence
+      CALL iom_rstput( 0, 0, kid, 'sishea', shear_i*1.0e8 )   ! Ice shear
       CALL iom_rstput( 0, 0, kid, 'si_amp', at_ip        )   ! Melt pond fraction
       CALL iom_rstput( 0, 0, kid, 'si_vmp', vt_ip        )   ! Melt pond volume
