Index: /NEMO/trunk/src/ICE/icedyn_rhg_eap.F90
===================================================================
--- /NEMO/trunk/src/ICE/icedyn_rhg_eap.F90	(revision 14116)
+++ /NEMO/trunk/src/ICE/icedyn_rhg_eap.F90	(revision 14117)
@@ -16,5 +16,5 @@
    !!                                           CICE code (Tsamados, Heorton)
    !!----------------------------------------------------------------------
-#if defined key_si3 && ! defined key_agrif
+#if defined key_si3 
    !!----------------------------------------------------------------------
    !!   'key_si3'                                       SI3 sea-ice model
@@ -66,5 +66,6 @@
    INTEGER ::   ncvgid   ! netcdf file id
    INTEGER ::   nvarid   ! netcdf variable id
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   aimsk00
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   eap_res  , aimsk15
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -202,10 +203,20 @@
       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) )
+      IF( kt == nit000 )  THEN 
+         !
+         ! for diagnostics 
+         ALLOCATE( aimsk00(jpi,jpj) )
+         ! for convergence tests
+         IF( nn_rhg_chkcvg > 0 ) ALLOCATE( eap_res(jpi,jpj), aimsk15(jpi,jpj) )
+      ENDIF
+      !
       DO_2D( 1, 1, 1, 1 )
-         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
+         aimsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
       END_2D
+      IF( nn_rhg_chkcvg > 0 ) THEN
+         DO_2D( 1, 1, 1, 1 )
+            aimsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
+         END_2D
+      ENDIF
       !
 !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
@@ -349,11 +360,11 @@
             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(ji,jj,Kmm)
+            zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
             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(ji,jj,Kmm)
+            zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
             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(ji,jj)
+            zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) )    ! if grounded icebergs are read: ocean depth = 0
             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_2D
@@ -819,17 +830,17 @@
             &                                  ztauy_ai, 'V', -1.0_wp, ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp )
          !
-         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 )
+         CALL iom_put( 'utau_oi' , ztaux_oi * aimsk00 )
+         CALL iom_put( 'vtau_oi' , ztauy_oi * aimsk00 )
+         CALL iom_put( 'utau_ai' , ztaux_ai * aimsk00 )
+         CALL iom_put( 'vtau_ai' , ztauy_ai * aimsk00 )
+         CALL iom_put( 'utau_bi' , ztaux_bi * aimsk00 )
+         CALL iom_put( 'vtau_bi' , ztauy_bi * aimsk00 )
       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
+      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * aimsk00 )   ! divergence
+      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * aimsk00 )   ! shear
+      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * aimsk00 )   ! delta
+      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * aimsk00 )   ! strength
 
       ! --- Stress tensor invariants (SIMIP diags) --- !
@@ -856,6 +867,6 @@
          !
          ! 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
+         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * aimsk00(:,:) ) ! Normal stress
+         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * aimsk00(:,:) ) ! Maximum shear stress
 
          DEALLOCATE ( zsig_I, zsig_II )
@@ -903,7 +914,7 @@
          CALL lbc_lnk_multi( 'icedyn_rhg_eap', zyield11, 'T', 1.0_wp, zyield22, 'T', 1.0_wp, zyield12, 'T', 1.0_wp )
 
-         CALL iom_put( 'yield11', zyield11 * zmsk00 )
-         CALL iom_put( 'yield22', zyield22 * zmsk00 )
-         CALL iom_put( 'yield12', zyield12 * zmsk00 )
+         CALL iom_put( 'yield11', zyield11 * aimsk00 )
+         CALL iom_put( 'yield22', zyield22 * aimsk00 )
+         CALL iom_put( 'yield12', zyield12 * aimsk00 )
       ENDIF
 
@@ -911,5 +922,5 @@
       IF( iom_use('aniso') ) THEN
          CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp )
-         CALL iom_put( 'aniso' , paniso_11 * zmsk00 )
+         CALL iom_put( 'aniso' , paniso_11 * aimsk00 )
       ENDIF
 
@@ -922,10 +933,10 @@
             &                                    zfU, 'U', -1.0_wp,   zfV, 'V', -1.0_wp )
 
-         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)
+         CALL iom_put( 'dssh_dx' , zspgU * aimsk00 )   ! Sea-surface tilt term in force balance (x)
+         CALL iom_put( 'dssh_dy' , zspgV * aimsk00 )   ! Sea-surface tilt term in force balance (y)
+         CALL iom_put( 'corstrx' , zCorU * aimsk00 )   ! Coriolis force term in force balance (x)
+         CALL iom_put( 'corstry' , zCorV * aimsk00 )   ! Coriolis force term in force balance (y)
+         CALL iom_put( 'intstrx' , zfU   * aimsk00 )   ! Internal force term in force balance (x)
+         CALL iom_put( 'intstry' , zfV   * aimsk00 )   ! Internal force term in force balance (y)
       ENDIF
 
@@ -938,6 +949,6 @@
          DO_2D( 0, 0, 0, 0 )
             ! 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)
+            zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * aimsk00(ji,jj)
+            zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * aimsk00(ji,jj)
 
             zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
@@ -973,13 +984,11 @@
             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(:,:) )
+                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * aimsk15(:,:) )
             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(:,:) )
+                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * aimsk15(:,:) )
             ENDIF
          ENDIF
       ENDIF
-      !
-      DEALLOCATE( zmsk00, zmsk15 )
       !
    END SUBROUTINE ice_dyn_rhg_eap
@@ -1006,5 +1015,4 @@
       REAL(wp)          ::   zresm           ! local real
       CHARACTER(len=20) ::   clname
-      REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence
       !!----------------------------------------------------------------------
 
@@ -1037,8 +1045,9 @@
       ELSE
          DO_2D( 1, 1, 1, 1 )
-            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)
+            eap_res(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) ) * aimsk15(ji,jj)
          END_2D
-         zresm = MAXVAL( zres )
+
+         zresm = MAXVAL( eap_res )
          CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
       ENDIF
@@ -1080,20 +1089,18 @@
       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) ::   zx, zy, zkxw, kyw, kaw
       REAL(wp) ::   zinvdx, zinvdy, zinvda
-      REAL(wp) ::   zdtemp1, zdtemp2, zatempprime, zinvsin
-
-      REAL(wp), PARAMETER ::   kfriction = 0.45_wp
-      !!---------------------------------------------------------------------
+      REAL(wp) ::   zdtemp1, zdtemp2, zatempprime
+
+      REAL(wp), PARAMETER ::   ppkfriction = 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
+!     REAL(wp), PARAMETER ::   ppinvstressconviso = 1._wp/(1._wp+ppkfriction*ppkfriction)
+      REAL(wp), PARAMETER ::   ppinvstressconviso = 1._wp
+
+      ! next statement uses pphi set in main module (icedyn_rhg_eap)
+      REAL(wp), PARAMETER ::   ppinvsin = 1._wp/sin(2._wp*pphi) * ppinvstressconviso
 
       ! compute eigenvalues, eigenvectors and angles for structure tensor, strain
@@ -1175,10 +1182,7 @@
 
       ! 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
+      zinvdx = real(nx_yield-1,kind=wp)/rpi
+      zinvdy = real(ny_yield-1,kind=wp)/rpi
+      zinvda = 2._wp*real(na_yield-1,kind=wp)
 
       ! % need 8 coords and 8 weights
@@ -1258,7 +1262,7 @@
       ! Tsamados 2013)
 
-      zsig11  = pstrength*(zstemp11r + kfriction*zstemp11s) * zinvsin
-      zsig12  = pstrength*(zstemp12r + kfriction*zstemp12s) * zinvsin
-      zsig22  = pstrength*(zstemp22r + kfriction*zstemp22s) * zinvsin
+      zsig11  = pstrength*(zstemp11r + ppkfriction*zstemp11s) * ppinvsin
+      zsig12  = pstrength*(zstemp12r + ppkfriction*zstemp12s) * ppinvsin
+      zsig22  = pstrength*(zstemp22r + ppkfriction*zstemp22s) * ppinvsin
 
       ! Back - rotation of the stress from principal axes into general coordinates
@@ -1319,7 +1323,7 @@
       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
+!!$   REAL (wp), PARAMETER ::   ppkfrac = 0.0001_wp   ! rate of fracture formation
+      REAL (wp), PARAMETER ::   ppkfrac = 1.e-3_wp    ! rate of fracture formation
+      REAL (wp), PARAMETER ::   ppthreshold = 0.3_wp  ! critical confinement ratio
       !!---------------------------------------------------------------
       !
@@ -1363,6 +1367,6 @@
       ! 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)
+         pmresult11 = - ppkfrac * (pa11 - zQ12Q12)
+         pmresult12 = - ppkfrac * (pa12 + zQ11Q12)
 
       ! Shear faulting
@@ -1370,7 +1374,7 @@
          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)
+      ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= ppthreshold)) THEN
+         pmresult11 = - ppkfrac * (pa11 - zQ12Q12)
+         pmresult12 = - ppkfrac * (pa12 + zQ11Q12)
 
       ! Horizontal spalling
@@ -1405,5 +1409,6 @@
       !!clem
       REAL(wp) ::   zw1, zw2, zfac, ztemp
-      REAL(wp) ::   idx, idy, idz
+      REAL(wp) ::   zidx, zidy, zidz
+      REAL(wp) ::   zsaak(6)                  ! temporary array
 
       REAL(wp), PARAMETER           ::   eps6 = 1.0e-6_wp
@@ -1522,33 +1527,36 @@
             zw2 = w2(ainit+ia*da)
             DO iz = 1, nz
-               idz = zinit+iz*dz
+               zidz = zinit+iz*dz
                ztemp = zw1 * EXP(-zw2*(zinit+iz*dz)*(zinit+iz*dz))
                DO iy = 1, ny_yield
-                  idy = yinit+iy*dy
+                  zidy = 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
+                     zidx = xinit+ix*dx
+                     call all_skr_sks(zidx,zidy,zidz,zsaak)
+                     zsaak = ztemp*zsaak*zfac
+                     s11r(ix,iy,ia) = s11r(ix,iy,ia) + zsaak(1)
+                     s12r(ix,iy,ia) = s12r(ix,iy,ia) + zsaak(2)
+                     s22r(ix,iy,ia) = s22r(ix,iy,ia) + zsaak(3)
+                     s11s(ix,iy,ia) = s11s(ix,iy,ia) + zsaak(4)
+                     s12s(ix,iy,ia) = s12s(ix,iy,ia) + zsaak(5)
+                     s22s(ix,iy,ia) = s22s(ix,iy,ia) + zsaak(6)
                   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
+            zidy = 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
+               zidx = xinit+ix*dx
+               call all_skr_sks(zidx,zidy,0._wp,zsaak)
+               zsaak = 0.5_wp*zsaak*zfac
+               s11r(ix,iy,ia) = zsaak(1)
+               s12r(ix,iy,ia) = zsaak(2)
+               s22r(ix,iy,ia) = zsaak(3)
+               s11s(ix,iy,ia) = zsaak(4)
+               s12s(ix,iy,ia) = zsaak(5)
+               s22s(ix,iy,ia) = zsaak(6)
             ENDDO
          ENDDO
@@ -1616,12 +1624,12 @@
    END FUNCTION w2
 
-   FUNCTION s11kr(px,py,pz)
-      !!-------------------------------------------------------------------
-      !! Function : s11kr
-      !!-------------------------------------------------------------------
+   SUBROUTINE all_skr_sks( px, py, pz, allsk )
       REAL(wp), INTENT(in   ) ::   px,py,pz
-
-      REAL(wp) ::   s11kr, zpih
-
+      REAL(wp), INTENT(out  ) ::   allsk(6)
+
+      REAL(wp) ::   zs12r0, zs21r0
+      REAL(wp) ::   zs12s0, zs21s0
+
+      REAL(wp) ::   zpih
       REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22
       REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22
@@ -1673,296 +1681,34 @@
       ENDIF
 
-      s11kr = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11)
-
-   END FUNCTION s11kr
-
-   FUNCTION s12kr(px,py,pz)
+      !!-------------------------------------------------------------------
+      !! Function : s11kr
+      !!-------------------------------------------------------------------
+      allsk(1) = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11)
       !!-------------------------------------------------------------------
       !! 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)
+      allsk(2)=0.5_wp*(zs12r0+zs21r0)
       !!-------------------------------------------------------------------
       !! 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)
+      allsk(3) = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22)
       !!-------------------------------------------------------------------
       !! 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)
+
+      allsk(4) = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11)
       !!-------------------------------------------------------------------
       !! 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)
+      allsk(5)=0.5_wp*(zs12s0+zs21s0)
       !!-------------------------------------------------------------------
       !! 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
+      allsk(6) = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22)
+   END SUBROUTINE all_skr_sks
 
 #else
Index: /NEMO/trunk/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90
===================================================================
--- /NEMO/trunk/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90	(revision 14116)
+++ /NEMO/trunk/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90	(revision 14117)
@@ -6,7 +6,7 @@
    !! 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 
+   !!             -   !  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.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)
@@ -14,5 +14,5 @@
    !!            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) 
+   !!                                           CICE code (Tsamados, Heorton)
    !!----------------------------------------------------------------------
 #if defined key_si3
@@ -30,6 +30,6 @@
    USE icevar         ! ice_var_sshdyn
    USE icedyn_rdgrft  ! sea-ice: ice strength
-   USE bdy_oce , ONLY : ln_bdy 
-   USE bdyice 
+   USE bdy_oce , ONLY : ln_bdy
+   USE bdyice
 #if defined key_agrif
    USE agrif_ice_interp
@@ -66,5 +66,6 @@
    INTEGER ::   ncvgid   ! netcdf file id
    INTEGER ::   nvarid   ! netcdf variable id
-   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   aimsk00
+   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   eap_res  , aimsk15
    !!----------------------------------------------------------------------
    !! NEMO/ICE 4.0 , NEMO Consortium (2018)
@@ -80,7 +81,7 @@
       !!
       !! ** 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 .	
+      !!  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
@@ -90,9 +91,9 @@
       !!                                 |
       !!                      (ji-1,jj)  |  (ji,jj)
-      !!                             ---------   
+      !!                             ---------
       !!                            |         |
       !!                            | (ji,jj) |------(ji,jj)
       !!                            |         |
-      !!                             ---------   
+      !!                             ---------
       !!                     (ji-1,jj-1)     (ji,jj-1)
       !!
@@ -101,5 +102,5 @@
       !!                snow total volume (vt_s) per unit area
       !!
-      !! ** Action  : - compute u_ice, v_ice : the components of the 
+      !! ** Action  : - compute u_ice, v_ice : the components of the
       !!                sea-ice velocity vector
       !!              - compute delta_i, shear_i, divu_i, which are inputs
@@ -107,5 +108,5 @@
       !!
       !! ** Steps   : 0) compute mask at F point
-      !!              1) Compute ice snow mass, ice strength 
+      !!              1) Compute ice snow mass, ice strength
       !!              2) Compute wind, oceanic stresses, mass terms and
       !!                 coriolis terms of the momentum equation
@@ -166,5 +167,5 @@
       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) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points
       !
       REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear
@@ -172,5 +173,5 @@
       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   
+      !                                                                 !    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
@@ -192,5 +193,5 @@
       !! --- diags
       REAL(wp) ::   zsig1, zsig2, zsig12, zfac, z1_strength
-      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig_I, zsig_II, zsig1_p, zsig2_p         
+      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)
@@ -199,15 +200,25 @@
       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)      
+      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) )
+      IF( kt == nit000 )  THEN 
+         !
+         ! for diagnostics 
+         ALLOCATE( aimsk00(jpi,jpj) )
+         ! for convergence tests
+         IF( nn_rhg_chkcvg > 0 ) ALLOCATE( eap_res(jpi,jpj), aimsk15(jpi,jpj) )
+      ENDIF
+      !
       DO_2D( 1, 1, 1, 1 )
-         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
+         aimsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice
       END_2D
+      IF( nn_rhg_chkcvg > 0 ) THEN
+         DO_2D( 1, 1, 1, 1 )
+            aimsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less
+         END_2D
+      ENDIF
       !
 !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization....
@@ -249,5 +260,5 @@
       ! 1) define some variables and initialize arrays
       !------------------------------------------------------------------------------!
-      zrhoco = rho0 * rn_cio 
+      zrhoco = rho0 * rn_cio
 !extra code for test case boundary conditions
       zinvw=1._wp/(zrhoco*0.5_wp)
@@ -270,7 +281,7 @@
       ENDIF
       z1_dtevp = 1._wp / zdtevp
-         
-      ! Initialise stress tensor 
-      zs1 (:,:) = pstress1_i (:,:) 
+
+      ! Initialise stress tensor
+      zs1 (:,:) = pstress1_i (:,:)
       zs2 (:,:) = pstress2_i (:,:)
       zs12(:,:) = pstress12_i(:,:)
@@ -319,9 +330,9 @@
          ! 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)
@@ -353,11 +364,11 @@
             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(ji,jj,Kmm)
+            zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
             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(ji,jj,Kmm)
+            zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0
             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(ji,jj)
+            zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) )    ! if grounded icebergs are read: ocean depth = 0
             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_2D
@@ -377,5 +388,5 @@
       !                                               ! ==================== !
       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
          !
@@ -404,5 +415,5 @@
                &   + 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)   &
@@ -410,5 +421,5 @@
                &    ) * 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)   &
@@ -418,9 +429,9 @@
 
             ! delta at T points
-            zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
+            zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )
 
          END_2D
          CALL lbc_lnk( 'icedyn_rhg_eap', zdelta, 'T', 1.0_wp )
-               
+
          ! P/delta at T points
          DO_2D( 1, 1, 1, 1 )
@@ -430,14 +441,14 @@
          DO_2D( 0, 1, 0, 1 )   ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2
 
-             ! shear at T points 
+             ! 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)
-            
+
            ! 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)   &
@@ -459,5 +470,5 @@
                zyield22(ji,jj) = 0.5_wp * (zstressptmp - zstressmtmp)
                zyield12(ji,jj) = zstress12tmp(ji,jj)
-               prdg_conv(ji,jj) = -min( zalphar, 0._wp)    
+               prdg_conv(ji,jj) = -min( zalphar, 0._wp)
             ENDIF
 
@@ -491,5 +502,5 @@
 
          DO_2D( 1, 0, 1, 0 )
-            ! stress12tmp at F points 
+            ! 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  )  &
@@ -504,5 +515,5 @@
                ! zalph2 = zalph2 - 1._wp
             ENDIF
-            
+
             ! stress at F points (zkt/=0 if landfast)
             zs12(ji,jj) = ( zs12(ji,jj) * zalph1 + zstress12tmpF ) * z1_alph1
@@ -570,5 +581,5 @@
                      &                                    + 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)                                                   & 
+                     &            + ( 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 )                                              &
@@ -630,5 +641,5 @@
                      &                                     + 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 
+                     &             ) * 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)
@@ -637,5 +648,5 @@
                      &                                    ) / 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 
+                     &              ) * 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
@@ -689,5 +700,5 @@
                      &                                     + 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 
+                     &             ) * 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)
@@ -744,5 +755,5 @@
                      &            + ( 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 )                                              & 
+                     &                                    ) / ( 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)
@@ -781,5 +792,5 @@
       !
       !------------------------------------------------------------------------------!
-      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 
+      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)
       !------------------------------------------------------------------------------!
       DO_2D( 1, 0, 1, 0 )
@@ -791,7 +802,7 @@
 
       END_2D
-      
+
       DO_2D( 0, 0, 0, 0 )
-         
+
          ! 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)   &
@@ -799,5 +810,5 @@
             &   ) * r1_e1e2t(ji,jj)
          zdt2 = zdt * zdt
-         
+
          zten_i(ji,jj) = zdt
 
@@ -806,5 +817,5 @@
             &   + 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 )
@@ -814,7 +825,7 @@
             &             + 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  
+         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
@@ -824,5 +835,5 @@
          &                                    zten_i, 'T', 1.0_wp, zs1    , 'T', 1.0_wp, zs2     , 'T', 1.0_wp, &
          &                                      zs12, 'F', 1.0_wp )
-      
+
       ! --- Store the stress tensor for the next time step --- !
       pstress1_i (:,:) = zs1 (:,:)
@@ -841,17 +852,17 @@
             &                                  ztauy_ai, 'V', -1.0_wp, ztaux_bi, 'U', -1.0_wp, ztauy_bi, 'V', -1.0_wp )
          !
-         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
-       
+         CALL iom_put( 'utau_oi' , ztaux_oi * aimsk00 )
+         CALL iom_put( 'vtau_oi' , ztauy_oi * aimsk00 )
+         CALL iom_put( 'utau_ai' , ztaux_ai * aimsk00 )
+         CALL iom_put( 'vtau_ai' , ztauy_ai * aimsk00 )
+         CALL iom_put( 'utau_bi' , ztaux_bi * aimsk00 )
+         CALL iom_put( 'vtau_bi' , ztauy_bi * aimsk00 )
+      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
+      IF( iom_use('icediv') )   CALL iom_put( 'icediv' , pdivu_i  * aimsk00 )   ! divergence
+      IF( iom_use('iceshe') )   CALL iom_put( 'iceshe' , pshear_i * aimsk00 )   ! shear
+      IF( iom_use('icedlt') )   CALL iom_put( 'icedlt' , pdelta_i * aimsk00 )   ! delta
+      IF( iom_use('icestr') )   CALL iom_put( 'icestr' , strength * aimsk00 )   ! strength
 
       ! --- Stress tensor invariants (SIMIP diags) --- !
@@ -859,28 +870,28 @@
          !
          ALLOCATE( zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
-         !         
+         !
          DO_2D( 1, 1, 1, 1 )
-         
+
             ! 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 ) 
+            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_2D
          !
          ! 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
-         
+         IF( iom_use('normstr') )   CALL iom_put( 'normstr', zsig_I (:,:) * aimsk00(:,:) ) ! Normal stress
+         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr', zsig_II(:,:) * aimsk00(:,:) ) ! Maximum shear stress
+
          DEALLOCATE ( zsig_I, zsig_II )
-         
+
       ENDIF
 
@@ -891,9 +902,9 @@
       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) )         
-         !         
+         ALLOCATE( zsig1_p(jpi,jpj) , zsig2_p(jpi,jpj) , zsig_I(jpi,jpj) , zsig_II(jpi,jpj) )
+         !
          DO_2D( 1, 1, 1, 1 )
-         
-            ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates 
+
+            ! Ice stresses computed with **viscosities** (delta, p/delta) at **previous** iterates
             !                        and **deformations** at current iterates
             !                        following Lemieux & Dupont (2020)
@@ -902,9 +913,9 @@
             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) )
@@ -913,9 +924,9 @@
          END_2D
          !
-         CALL iom_put( 'sig1_pnorm' , zsig1_p ) 
-         CALL iom_put( 'sig2_pnorm' , zsig2_p ) 
-      
+         CALL iom_put( 'sig1_pnorm' , zsig1_p )
+         CALL iom_put( 'sig2_pnorm' , zsig2_p )
+
          DEALLOCATE( zsig1_p , zsig2_p , zsig_I, zsig_II )
-         
+
       ENDIF
 
@@ -925,15 +936,15 @@
          CALL lbc_lnk_multi( 'icedyn_rhg_eap', zyield11, 'T', 1.0_wp, zyield22, 'T', 1.0_wp, zyield12, 'T', 1.0_wp )
 
-         CALL iom_put( 'yield11', zyield11 * zmsk00 )
-         CALL iom_put( 'yield22', zyield22 * zmsk00 )
-         CALL iom_put( 'yield12', zyield12 * zmsk00 )
+         CALL iom_put( 'yield11', zyield11 * aimsk00 )
+         CALL iom_put( 'yield22', zyield22 * aimsk00 )
+         CALL iom_put( 'yield12', zyield12 * aimsk00 )
       ENDIF
 
       ! --- anisotropy tensor --- !
-      IF( iom_use('aniso') ) THEN       
-         CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp ) 
-         CALL iom_put( 'aniso' , paniso_11 * zmsk00 )
-      ENDIF
-      
+      IF( iom_use('aniso') ) THEN
+         CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1.0_wp )
+         CALL iom_put( 'aniso' , paniso_11 * aimsk00 )
+      ENDIF
+
       ! --- SIMIP --- !
       IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. &
@@ -944,10 +955,10 @@
             &                                    zfU, 'U', -1.0_wp,   zfV, 'V', -1.0_wp )
 
-         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)
+         CALL iom_put( 'dssh_dx' , zspgU * aimsk00 )   ! Sea-surface tilt term in force balance (x)
+         CALL iom_put( 'dssh_dy' , zspgV * aimsk00 )   ! Sea-surface tilt term in force balance (y)
+         CALL iom_put( 'corstrx' , zCorU * aimsk00 )   ! Coriolis force term in force balance (x)
+         CALL iom_put( 'corstry' , zCorV * aimsk00 )   ! Coriolis force term in force balance (y)
+         CALL iom_put( 'intstrx' , zfU   * aimsk00 )   ! Internal force term in force balance (x)
+         CALL iom_put( 'intstry' , zfV   * aimsk00 )   ! Internal force term in force balance (y)
       ENDIF
 
@@ -960,6 +971,6 @@
          DO_2D( 0, 0, 0, 0 )
             ! 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)
+            zfac_x = 0.5 * u_ice(ji,jj) * e2u(ji,jj) * aimsk00(ji,jj)
+            zfac_y = 0.5 * v_ice(ji,jj) * e1v(ji,jj) * aimsk00(ji,jj)
 
             zdiag_xmtrp_ice(ji,jj) = rhoi * zfac_x * ( vt_i(ji+1,jj) + vt_i(ji,jj) ) ! ice mass transport, X-component
@@ -979,5 +990,5 @@
 
          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( '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
@@ -995,21 +1006,19 @@
             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(:,:) )
+                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * aimsk15(:,:) )
             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(:,:) )
+                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * aimsk15(:,:) )
             ENDIF
          ENDIF
-      ENDIF      
-      !
-      DEALLOCATE( zmsk00, zmsk15 )
+      ENDIF
       !
    END SUBROUTINE ice_dyn_rhg_eap
 
-   
+
    SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb )
       !!----------------------------------------------------------------------
       !!                    ***  ROUTINE rhg_cvg  ***
-      !!                     
+      !!
       !! ** Purpose :   check convergence of oce rheology
       !!
@@ -1019,5 +1028,5 @@
       !!                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)   
+      !! ** 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
@@ -1026,7 +1035,6 @@
       INTEGER           ::   it, idtime, istatus
       INTEGER           ::   ji, jj          ! dummy loop indices
-      REAL(wp)          ::   zresm           ! local real 
+      REAL(wp)          ::   zresm           ! local real
       CHARACTER(len=20) ::   clname
-      REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence
       !!----------------------------------------------------------------------
 
@@ -1053,5 +1061,5 @@
       ! time
       it = ( kt - 1 ) * kitermax + kiter
-      
+
       ! convergence
       IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large)
@@ -1059,13 +1067,13 @@
       ELSE
          DO_2D( 1, 1, 1, 1 )
-            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)
+            eap_res(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) ) * aimsk15(ji,jj)
+            ! cut of the boundary of the box (forced velocities)
+            IF (mjg(jj)<=30 .or. mjg(jj)>970 .or. mig(ji)<=30 .or. mig(ji)>970) THEN
+               eap_res(ji,jj) = 0._wp
+            END IF
          END_2D
 
-         ! cut of the boundary of the box (forced velocities)
-         IF (mjg(jj)<=30 .or. mjg(jj)>970 .or. mig(ji)<=30 .or. mig(ji)>970) THEN
-            zres(ji,jj) = 0._wp
-         END IF
-         zresm = MAXVAL( zres )
+         zresm = MAXVAL( eap_res )
          CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain
       ENDIF
@@ -1077,5 +1085,5 @@
          IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid)
       ENDIF
-      
+
    END SUBROUTINE rhg_cvg
 
@@ -1085,5 +1093,5 @@
       !!---------------------------------------------------------------------
       !!                   ***  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
@@ -1098,30 +1106,28 @@
       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) ::   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
-      !!---------------------------------------------------------------------
+      REAL(wp) ::   zsig11, zsig12, zsig22
+      REAL(wp) ::   zsgprm11, zsgprm12, zsgprm22
+      REAL(wp) ::   zAngle_denom_gamma,  zAngle_denom_alpha
+      REAL(wp) ::   zTany_1, zTany_2
+      REAL(wp) ::   zx, zy, zkxw, kyw, kaw
+      REAL(wp) ::   zinvdx, zinvdy, zinvda
+      REAL(wp) ::   zdtemp1, zdtemp2, zatempprime
+
+      REAL(wp), PARAMETER ::   ppkfriction = 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
-       
+!     REAL(wp), PARAMETER ::   ppinvstressconviso = 1._wp/(1._wp+ppkfriction*ppkfriction)
+      REAL(wp), PARAMETER ::   ppinvstressconviso = 1._wp
+
+      ! next statement uses pphi set in main module (icedyn_rhg_eap)
+      REAL(wp), PARAMETER ::   ppinvsin = 1._wp/sin(2._wp*pphi) * ppinvstressconviso
+
       ! compute eigenvalues, eigenvectors and angles for structure tensor, strain
       ! rates
@@ -1132,22 +1138,22 @@
       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      
+
+      ! 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 
+
+         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
- 
+         zQ11Q12 = pa12*zAngle_denom_gamma                     !CosSin
+      ENDIF
+
       ! rotation Q*atemp*Q^T
-      zatempprime = zQ11Q11*pa11 + 2.0_wp*zQ11Q12*pa12 + zQ12Q12*za22 
-      
+      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)
@@ -1155,5 +1161,5 @@
       zdtemp22 = 0.5_wp*(pdivu - ptension)
 
-      ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22) 
+      ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22)
 
       zQd11Qd11 = 1.0_wp
@@ -1166,5 +1172,5 @@
                              ( zdtemp11 - zdtemp22 ) + 4.0_wp*zdtemp12*zdtemp12)
 
-         zQd11Qd11 = 0.5_wp + ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Cos^2 
+         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
@@ -1177,35 +1183,32 @@
       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 
+      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 
+      ! 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
+      ! 3) update anisotropic stress tensor
+      zinvdx = real(nx_yield-1,kind=wp)/rpi
+      zinvdy = real(ny_yield-1,kind=wp)/rpi
+      zinvda = 2._wp*real(na_yield-1,kind=wp)
 
       ! % need 8 coords and 8 weights
@@ -1213,14 +1216,14 @@
       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  
-      
+      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 
-      
+      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  ) &
@@ -1248,5 +1251,5 @@
 !!$        & + 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  ) &
@@ -1280,12 +1283,12 @@
       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
+      zsig11  = pstrength*(zstemp11r + ppkfriction*zstemp11s) * ppinvsin
+      zsig12  = pstrength*(zstemp12r + ppkfriction*zstemp12s) * ppinvsin
+      zsig22  = pstrength*(zstemp22r + ppkfriction*zstemp22s) * ppinvsin
 
       ! Back - rotation of the stress from principal axes into general coordinates
@@ -1300,5 +1303,5 @@
       pstressm  = zsgprm11 - zsgprm22
 
-      ! Compute ridging and sliding functions in general coordinates 
+      ! Compute ridging and sliding functions in general coordinates
       ! (Equation 11 in Tsamados 2013)
       IF (ksub == kndte) THEN
@@ -1307,5 +1310,5 @@
          zrotstemp12r = zQ11Q11*zstemp12r +       zQ11Q12*(zstemp11r-zstemp22r) &
                       - zQ12Q12*zstemp12r
-         zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r & 
+         zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r &
                       + zQ11Q11*zstemp22r
 
@@ -1314,5 +1317,5 @@
          zrotstemp12s = zQ11Q11*zstemp12s +       zQ11Q12*(zstemp11s-zstemp22s) &
                       - zQ12Q12*zstemp12s
-         zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s & 
+         zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s &
                       + zQ11Q11*zstemp22s
 
@@ -1322,5 +1325,5 @@
                   + zrotstemp22s*zdtemp22
       ENDIF
-   END SUBROUTINE update_stress_rdg 
+   END SUBROUTINE update_stress_rdg
 
 !=======================================================================
@@ -1331,5 +1334,5 @@
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE calc_ffrac  ***
-      !!                     
+      !!
       !! ** Purpose :   Computes term in evolution equation for structure tensor
       !!                which determines the ice floe re-orientation due to fracture
@@ -1346,20 +1349,20 @@
       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 
+!!$   REAL (wp), PARAMETER ::   ppkfrac = 0.0001_wp   ! rate of fracture formation
+      REAL (wp), PARAMETER ::   ppkfrac = 1.e-3_wp    ! rate of fracture formation
+      REAL (wp), PARAMETER ::   ppthreshold = 0.3_wp  ! critical confinement ratio
       !!---------------------------------------------------------------
       !
-      zsigma11 = 0.5_wp*(pstressp+pstressm) 
-      zsigma12 = pstress12 
-      zsigma22 = 0.5_wp*(pstressp-pstressm) 
+      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 
+
+      ! 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
 
@@ -1373,7 +1376,7 @@
                        zsigma22 ) + 4.0_wp*zsigma12*zsigma12)
 
-         zQ11Q11 = 0.5_wp + ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Cos^2 
+         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 
+         zQ11Q12 = zsigma12*zAngle_denom                      !CosSin
       ENDIF
 
@@ -1390,6 +1393,6 @@
       ! 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)
+         pmresult11 = - ppkfrac * (pa11 - zQ12Q12)
+         pmresult12 = - ppkfrac * (pa12 + zQ11Q12)
 
       ! Shear faulting
@@ -1397,10 +1400,10 @@
          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)
+      ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= ppthreshold)) THEN
+         pmresult11 = - ppkfrac * (pa11 - zQ12Q12)
+         pmresult12 = - ppkfrac * (pa12 + zQ11Q12)
 
       ! Horizontal spalling
-      ELSE 
+      ELSE
          pmresult11 = 0.0_wp
          pmresult12 = 0.0_wp
@@ -1413,7 +1416,7 @@
       !!---------------------------------------------------------------------
       !!                   ***  ROUTINE rhg_eap_rst  ***
-      !!                     
+      !!
       !! ** Purpose :   Read or write RHG file in restart file
-      !!      
+      !!
       !! ** Method  :   use of IOM library
       !!----------------------------------------------------------------------
@@ -1424,7 +1427,7 @@
       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
@@ -1432,5 +1435,6 @@
       !!clem
       REAL(wp) ::   zw1, zw2, zfac, ztemp
-      REAL(wp) ::   idx, idy, idz
+      REAL(wp) ::   zidx, zidy, zidz
+      REAL(wp) ::   zsaak(6)                  ! temporary array
 
       REAL(wp), PARAMETER           ::   eps6 = 1.0e-6_wp
@@ -1508,5 +1512,5 @@
 !!$                        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) 
+!!$                           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))* &
@@ -1543,5 +1547,5 @@
 !!$         ENDDO
 
-         !! faster but still very slow => to be improved         
+         !! faster but still very slow => to be improved
          zfac = dz/sin(2._wp*pphi)
          DO ia = 1, na_yield-1
@@ -1549,33 +1553,36 @@
             zw2 = w2(ainit+ia*da)
             DO iz = 1, nz
-               idz = zinit+iz*dz
+               zidz = zinit+iz*dz
                ztemp = zw1 * EXP(-zw2*(zinit+iz*dz)*(zinit+iz*dz))
                DO iy = 1, ny_yield
-                  idy = yinit+iy*dy
+                  zidy = 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
+                     zidx = xinit+ix*dx
+                     call all_skr_sks(zidx,zidy,zidz,zsaak)
+                     zsaak = ztemp*zsaak*zfac
+                     s11r(ix,iy,ia) = s11r(ix,iy,ia) + zsaak(1)
+                     s12r(ix,iy,ia) = s12r(ix,iy,ia) + zsaak(2)
+                     s22r(ix,iy,ia) = s22r(ix,iy,ia) + zsaak(3)
+                     s11s(ix,iy,ia) = s11s(ix,iy,ia) + zsaak(4)
+                     s12s(ix,iy,ia) = s12s(ix,iy,ia) + zsaak(5)
+                     s22s(ix,iy,ia) = s22s(ix,iy,ia) + zsaak(6)
                   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
+            zidy = 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
+               zidx = xinit+ix*dx
+               call all_skr_sks(zidx,zidy,0._wp,zsaak)
+               zsaak = 0.5_wp*zsaak*zfac
+               s11r(ix,iy,ia) = zsaak(1)
+               s12r(ix,iy,ia) = zsaak(2)
+               s22r(ix,iy,ia) = zsaak(3)
+               s11s(ix,iy,ia) = zsaak(4)
+               s12s(ix,iy,ia) = zsaak(5)
+               s22s(ix,iy,ia) = zsaak(6)
             ENDDO
          ENDDO
@@ -1611,5 +1618,5 @@
       REAL(wp) ::   w1
       !!-------------------------------------------------------------------
-   
+
       w1 = -   223.87569446_wp &
        &   +  2361.21986630_wp*pa &
@@ -1620,5 +1627,5 @@
        &   - 16789.98003081_wp*pa*pa*pa*pa*pa*pa &
        &   +  3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa
-   
+
    END FUNCTION w1
 
@@ -1631,5 +1638,5 @@
       REAL(wp) ::   w2
       !!-------------------------------------------------------------------
-   
+
       w2 = -    6670.68911883_wp &
        &   +   70222.33061536_wp*pa &
@@ -1640,15 +1647,15 @@
        &   -  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
-      !!-------------------------------------------------------------------
+   SUBROUTINE all_skr_sks( px, py, pz, allsk )
       REAL(wp), INTENT(in   ) ::   px,py,pz
-   
-      REAL(wp) ::   s11kr, zpih
-   
+      REAL(wp), INTENT(out  ) ::   allsk(6)
+
+      REAL(wp) ::   zs12r0, zs21r0
+      REAL(wp) ::   zs12s0, zs21s0
+
+      REAL(wp) ::   zpih
       REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22
       REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22
@@ -1659,7 +1666,7 @@
       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)
@@ -1687,5 +1694,5 @@
       zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22
       zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22
-   
+
       IF (-zIIn1t2>=rsmall) THEN
       zHen1t2 = 1._wp
@@ -1693,5 +1700,5 @@
       zHen1t2 = 0._wp
       ENDIF
-   
+
       IF (-zIIn2t1>=rsmall) THEN
       zHen2t1 = 1._wp
@@ -1699,297 +1706,35 @@
       zHen2t1 = 0._wp
       ENDIF
-   
-      s11kr = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11)
-
-   END FUNCTION s11kr
-
-   FUNCTION s12kr(px,py,pz)
+
+      !!-------------------------------------------------------------------
+      !! Function : s11kr
+      !!-------------------------------------------------------------------
+      allsk(1) = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11)
       !!-------------------------------------------------------------------
       !! 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)
+      allsk(2)=0.5_wp*(zs12r0+zs21r0)
       !!-------------------------------------------------------------------
       !! 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)
+      allsk(3) = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22)
       !!-------------------------------------------------------------------
       !! 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)
+
+      allsk(4) = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11)
       !!-------------------------------------------------------------------
       !! 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) 
+      allsk(5)=0.5_wp*(zs12s0+zs21s0)
       !!-------------------------------------------------------------------
       !! 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
+      allsk(6) = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22)
+   END SUBROUTINE all_skr_sks
 
 #else
@@ -1997,4 +1742,20 @@
    !!   Default option         Empty module           NO SI3 sea-ice model
    !!----------------------------------------------------------------------
+   USE par_kind
+   USE lib_mpp
+   CONTAINS
+   SUBROUTINE ice_dyn_rhg_eap( kt, Kmm, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, paniso_11, paniso_12, prdg_conv )
+      INTEGER                 , INTENT(in   ) ::   kt                                    ! time step
+      INTEGER                 , INTENT(in   ) ::   Kmm                                   ! ocean time level index
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pstress1_i, pstress2_i, pstress12_i   !
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   pshear_i  , pdivu_i   , pdelta_i      !
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   paniso_11 , paniso_12                 ! structure tensor components
+      REAL(wp), DIMENSION(:,:), INTENT(in   ) ::   prdg_conv                             ! for ridging
+      CALL ctl_stop('EAP rheology is currently dsabled due to issues with AGRIF preprocessing')
+   END SUBROUTINE ice_dyn_rhg_eap
+   SUBROUTINE rhg_eap_rst( cdrw, kt )
+      CHARACTER(len=*) , INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag
+      INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
+   END SUBROUTINE rhg_eap_rst
 #endif
 
