Index: NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90
===================================================================
--- NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90	(revision 11911)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90	(revision 11920)
@@ -46,4 +46,11 @@
    PUBLIC   ice_dyn_rhg_eap   ! called by icedyn_rhg.F90
    PUBLIC   rhg_eap_rst       ! called by icedyn_rhg.F90
+
+   ! look-up table for calculating structure tensor
+   INTEGER, PARAMETER ::   nx_yield = 41
+   INTEGER, PARAMETER ::   ny_yield = 41
+   INTEGER, PARAMETER ::   na_yield = 21
+
+   REAL(wp), DIMENSION(nx_yield, ny_yield, na_yield) ::   s11r, s12r, s22r, s11s, s12s, s22s
 
    !! * Substitutions
@@ -872,6 +879,15 @@
       INTEGER, OPTIONAL, INTENT(in) ::   kt     ! ice time-step
       !
-      INTEGER  ::   iter            ! local integer
+      INTEGER  ::   iter                      ! local integer
       INTEGER  ::   id1, id2, id3, id4, id5   ! local integers
+      INTEGER  ::   ix, iy, ip, iz, n, ia     ! local integers
+    
+      INTEGER, PARAMETER            ::    nz = 100
+      
+      REAL(wp) ::   ainit, xinit, yinit, pinit, zinit
+      REAL(wp) ::   da, dx, dy, dp, dz, a1
+
+      REAL(wp), PARAMETER           ::   eps6 = 1.0e-6_wp
+      REAL(wp), PARAMETER           ::   phi = 3.141592653589793_wp/12._wp    ! diamond shaped floe smaller angle (default phi = 30 deg)
       !!----------------------------------------------------------------------
       !
@@ -911,4 +927,68 @@
          ENDIF
          !
+
+      da = 0.5_wp/real(na_yield-1,kind=wp)
+      ainit = 0.5_wp - da
+      dx = rpi/real(nx_yield-1,kind=wp)
+      xinit = rpi + 0.25_wp*rpi - dx
+      dz = rpi/real(nz,kind=wp)
+      zinit = -rpi*0.5_wp
+      dy = rpi/real(ny_yield-1,kind=wp)
+      yinit = -dy
+
+      DO ia=1,na_yield
+       DO ix=1,nx_yield
+        DO iy=1,ny_yield
+         s11r(ix,iy,ia) = 0._wp
+         s12r(ix,iy,ia) = 0._wp
+         s22r(ix,iy,ia) = 0._wp
+         s11s(ix,iy,ia) = 0._wp
+         s12s(ix,iy,ia) = 0._wp
+         s22s(ix,iy,ia) = 0._wp
+         IF (ia <= na_yield-1) then
+          DO iz=1,nz
+          s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+           s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)
+          s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+           s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)
+          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,phi)*dz/sin(2._wp*phi) 
+          s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+           s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)
+          s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+           s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)
+          s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* &
+           exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* &
+           s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz,phi)*dz/sin(2._wp*phi)
+          ENDDO
+          IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp
+          IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp
+          IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp
+          IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp
+          IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp
+          IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp
+         ELSE
+          s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi)
+          s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi)
+          s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi)
+          s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi)
+          s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi)
+          s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp,phi)/sin(2._wp*phi)
+          IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp
+          IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp
+          IF (abs(s22r(ix,iy,ia)) < eps6) s22r(ix,iy,ia) = 0._wp
+          IF (abs(s11s(ix,iy,ia)) < eps6) s11s(ix,iy,ia) = 0._wp
+          IF (abs(s12s(ix,iy,ia)) < eps6) s12s(ix,iy,ia) = 0._wp
+          IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp
+         ENDIF
+        ENDDO
+       ENDDO
+      ENDDO
+
       ELSEIF( TRIM(cdrw) == 'WRITE' ) THEN   ! Create restart file
          !                                   ! -------------------
@@ -926,10 +1006,406 @@
    END SUBROUTINE rhg_eap_rst
 
+
+   FUNCTION w1(a)
+      !!-------------------------------------------------------------------
+      !! Function : w1 (see Gaussian function psi in Tsamados et al 2013)
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   a
+      REAL(wp) ::   w1
+      !!-------------------------------------------------------------------
+   
+      w1 = - 223.87569446_wp &
+           + 2361.2198663_wp*a &
+           - 10606.56079975_wp*a*a &
+           + 26315.50025642_wp*a*a*a &
+           - 38948.30444297_wp*a*a*a*a &
+           + 34397.72407466_wp*a*a*a*a*a &
+           - 16789.98003081_wp*a*a*a*a*a*a &
+           + 3495.82839237_wp*a*a*a*a*a*a*a
+   
+   END FUNCTION w1
+
+
+   FUNCTION w2(a)
+      !!-------------------------------------------------------------------
+      !! Function : w2 (see Gaussian function psi in Tsamados et al 2013)
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   a
+      REAL(wp) ::   w2
+      !!-------------------------------------------------------------------
+   
+      w2 = - 6670.68911883_wp &
+           + 70222.33061536_wp*a &
+           - 314871.71525448_wp*a*a &
+           + 779570.02793492_wp*a*a*a &
+           - 1151098.82436864_wp*a*a*a*a &
+           + 1013896.59464498_wp*a*a*a*a*a &
+           - 493379.44906738_wp*a*a*a*a*a*a &
+           + 102356.551518_wp*a*a*a*a*a*a*a
+   
+   END FUNCTION w2
+
+   FUNCTION s11kr(x,y,z,phi) 
+      !!-------------------------------------------------------------------
+      !! Function : s11kr
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   x,y,z,phi
+   
+      REAL(wp) ::   s11kr, p, pih
+   
+      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
+      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
+      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22
+      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22
+      REAL(wp) ::   d11, d12, d22
+      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2
+      REAL(wp) ::   Hen1t2, Hen2t1
+      !!-------------------------------------------------------------------
+   
+      pih = 0.5_wp*rpi
+      p = phi
+   
+      n1t2i11 = cos(z+pih-p) * cos(z+p)
+      n1t2i12 = cos(z+pih-p) * sin(z+p)
+      n1t2i21 = sin(z+pih-p) * cos(z+p)
+      n1t2i22 = sin(z+pih-p) * sin(z+p)
+      n2t1i11 = cos(z-pih+p) * cos(z-p)
+      n2t1i12 = cos(z-pih+p) * sin(z-p)
+      n2t1i21 = sin(z-pih+p) * cos(z-p)
+      n2t1i22 = sin(z-pih+p) * sin(z-p)
+      t1t2i11 = cos(z-p) * cos(z+p)
+      t1t2i12 = cos(z-p) * sin(z+p)
+      t1t2i21 = sin(z-p) * cos(z+p)
+      t1t2i22 = sin(z-p) * sin(z+p)
+      t2t1i11 = cos(z+p) * cos(z-p)
+      t2t1i12 = cos(z+p) * sin(z-p)
+      t2t1i21 = sin(z+p) * cos(z-p)
+      t2t1i22 = sin(z+p) * sin(z-p)
+   ! In expression of tensor d, with this formulatin d(x)=-d(x+pi)
+   ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else
+   ! x=x-pi/2
+      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
+      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
+      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
+      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
+      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
+      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
+   
+      IF (-IIn1t2>=rsmall) then
+      Hen1t2 = 1._wp
+      ELSE
+      Hen1t2 = 0._wp
+      ENDIF
+   
+      IF (-IIn2t1>=rsmall) then
+      Hen2t1 = 1._wp
+      ELSE
+      Hen2t1 = 0._wp
+      ENDIF
+   
+      s11kr = (- Hen1t2 * n1t2i11 - Hen2t1 * n2t1i11)
+
+   END FUNCTION s11kr
+
+   FUNCTION s12kr(x,y,z,phi)
+      !!-------------------------------------------------------------------
+      !! Function : s12kr
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   x,y,z,phi
+
+      REAL(wp) ::   s12kr, s12r0, s21r0, p, pih
+
+      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
+      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
+      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22
+      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22
+      REAL(wp) ::   d11, d12, d22
+      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2
+      REAL(wp) ::   Hen1t2, Hen2t1
+      !!-------------------------------------------------------------------
+      pih = 0.5_wp*rpi
+      p = phi
+   
+      n1t2i11 = cos(z+pih-p) * cos(z+p)
+      n1t2i12 = cos(z+pih-p) * sin(z+p)
+      n1t2i21 = sin(z+pih-p) * cos(z+p)
+      n1t2i22 = sin(z+pih-p) * sin(z+p)
+      n2t1i11 = cos(z-pih+p) * cos(z-p)
+      n2t1i12 = cos(z-pih+p) * sin(z-p)
+      n2t1i21 = sin(z-pih+p) * cos(z-p)
+      n2t1i22 = sin(z-pih+p) * sin(z-p)
+      t1t2i11 = cos(z-p) * cos(z+p)
+      t1t2i12 = cos(z-p) * sin(z+p)
+      t1t2i21 = sin(z-p) * cos(z+p)
+      t1t2i22 = sin(z-p) * sin(z+p)
+      t2t1i11 = cos(z+p) * cos(z-p)
+      t2t1i12 = cos(z+p) * sin(z-p)
+      t2t1i21 = sin(z+p) * cos(z-p)
+      t2t1i22 = sin(z+p) * sin(z-p)
+      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
+      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
+      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
+      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
+      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
+      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
+   
+      IF (-IIn1t2>=rsmall) then
+      Hen1t2 = 1._wp
+      ELSE
+      Hen1t2 = 0._wp
+      ENDIF
+   
+      IF (-IIn2t1>=rsmall) then
+      Hen2t1 = 1._wp
+      ELSE
+      Hen2t1 = 0._wp
+      ENDIF
+   
+      s12r0 = (- Hen1t2 * n1t2i12 - Hen2t1 * n2t1i12)
+      s21r0 = (- Hen1t2 * n1t2i21 - Hen2t1 * n2t1i21)
+      s12kr=0.5_wp*(s12r0+s21r0)
+   
+   END FUNCTION s12kr
+
+   FUNCTION s22kr(x,y,z,phi)
+      !!-------------------------------------------------------------------
+      !! Function : s22kr
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   x,y,z,phi
+
+      REAL(wp) ::   s22kr, p, pih
+
+      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
+      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
+      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22
+      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22
+      REAL(wp) ::   d11, d12, d22
+      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2
+      REAL(wp) ::   Hen1t2, Hen2t1
+      !!-------------------------------------------------------------------
+      pih = 0.5_wp*rpi
+      p = phi
+
+      n1t2i11 = cos(z+pih-p) * cos(z+p)
+      n1t2i12 = cos(z+pih-p) * sin(z+p)
+      n1t2i21 = sin(z+pih-p) * cos(z+p)
+      n1t2i22 = sin(z+pih-p) * sin(z+p)
+      n2t1i11 = cos(z-pih+p) * cos(z-p)
+      n2t1i12 = cos(z-pih+p) * sin(z-p)
+      n2t1i21 = sin(z-pih+p) * cos(z-p)
+      n2t1i22 = sin(z-pih+p) * sin(z-p)
+      t1t2i11 = cos(z-p) * cos(z+p)
+      t1t2i12 = cos(z-p) * sin(z+p)
+      t1t2i21 = sin(z-p) * cos(z+p)
+      t1t2i22 = sin(z-p) * sin(z+p)
+      t2t1i11 = cos(z+p) * cos(z-p)
+      t2t1i12 = cos(z+p) * sin(z-p)
+      t2t1i21 = sin(z+p) * cos(z-p)
+      t2t1i22 = sin(z+p) * sin(z-p)
+      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
+      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
+      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
+      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
+      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
+      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
+
+      IF (-IIn1t2>=rsmall) then
+      Hen1t2 = 1._wp
+      ELSE
+      Hen1t2 = 0._wp
+      ENDIF
+
+      IF (-IIn2t1>=rsmall) then
+      Hen2t1 = 1._wp
+      ELSE
+      Hen2t1 = 0._wp
+      ENDIF
+
+      s22kr = (- Hen1t2 * n1t2i22 - Hen2t1 * n2t1i22)
+
+   END FUNCTION s22kr
+
+   FUNCTION s11ks(x,y,z,phi)
+      !!-------------------------------------------------------------------
+      !! Function : s11ks
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   x,y,z,phi
+
+      REAL(wp) ::   s11ks, p, pih
+
+      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
+      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
+      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22
+      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22
+      REAL(wp) ::   d11, d12, d22
+      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2
+      REAL(wp) ::   Hen1t2, Hen2t1
+      !!-------------------------------------------------------------------
+      pih = 0.5_wp*rpi
+      p = phi
+
+      n1t2i11 = cos(z+pih-p) * cos(z+p)
+      n1t2i12 = cos(z+pih-p) * sin(z+p)
+      n1t2i21 = sin(z+pih-p) * cos(z+p)
+      n1t2i22 = sin(z+pih-p) * sin(z+p)
+      n2t1i11 = cos(z-pih+p) * cos(z-p)
+      n2t1i12 = cos(z-pih+p) * sin(z-p)
+      n2t1i21 = sin(z-pih+p) * cos(z-p)
+      n2t1i22 = sin(z-pih+p) * sin(z-p)
+      t1t2i11 = cos(z-p) * cos(z+p)
+      t1t2i12 = cos(z-p) * sin(z+p)
+      t1t2i21 = sin(z-p) * cos(z+p)
+      t1t2i22 = sin(z-p) * sin(z+p)
+      t2t1i11 = cos(z+p) * cos(z-p)
+      t2t1i12 = cos(z+p) * sin(z-p)
+      t2t1i21 = sin(z+p) * cos(z-p)
+      t2t1i22 = sin(z+p) * sin(z-p)
+      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
+      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
+      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
+      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
+      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
+      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
+
+      IF (-IIn1t2>=rsmall) then
+      Hen1t2 = 1._wp
+      ELSE
+      Hen1t2 = 0._wp
+      ENDIF
+
+      IF (-IIn2t1>=rsmall) then
+      Hen2t1 = 1._wp
+      ELSE
+      Hen2t1 = 0._wp
+      ENDIF
+
+      s11ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i11 + Hen2t1 * t2t1i11)
+
+   END FUNCTION s11ks
+
+   FUNCTION s12ks(x,y,z,phi)
+      !!-------------------------------------------------------------------
+      !! Function : s12ks
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   x,y,z,phi
+
+      REAL(wp) ::   s12ks, s12s0, s21s0, p, pih
+
+      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
+      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
+      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22
+      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22
+      REAL(wp) ::   d11, d12, d22
+      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2
+      REAL(wp) ::   Hen1t2, Hen2t1
+      !!-------------------------------------------------------------------
+      pih = 0.5_wp*rpi
+      p =phi
+
+      n1t2i11 = cos(z+pih-p) * cos(z+p)
+      n1t2i12 = cos(z+pih-p) * sin(z+p)
+      n1t2i21 = sin(z+pih-p) * cos(z+p)
+      n1t2i22 = sin(z+pih-p) * sin(z+p)
+      n2t1i11 = cos(z-pih+p) * cos(z-p)
+      n2t1i12 = cos(z-pih+p) * sin(z-p)
+      n2t1i21 = sin(z-pih+p) * cos(z-p)
+      n2t1i22 = sin(z-pih+p) * sin(z-p)
+      t1t2i11 = cos(z-p) * cos(z+p)
+      t1t2i12 = cos(z-p) * sin(z+p)
+      t1t2i21 = sin(z-p) * cos(z+p)
+      t1t2i22 = sin(z-p) * sin(z+p)
+      t2t1i11 = cos(z+p) * cos(z-p)
+      t2t1i12 = cos(z+p) * sin(z-p)
+      t2t1i21 = sin(z+p) * cos(z-p)
+      t2t1i22 = sin(z+p) * sin(z-p)
+      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
+      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
+      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
+      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
+      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
+      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
+
+      IF (-IIn1t2>=rsmall) then
+      Hen1t2 = 1._wp
+      ELSE
+      Hen1t2 = 0._wp
+      ENDIF
+
+      IF (-IIn2t1>=rsmall) then
+      Hen2t1 = 1._wp
+      ELSE
+      Hen2t1 = 0._wp
+      ENDIF
+
+      s12s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i12 + Hen2t1 * t2t1i12)
+      s21s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i21 + Hen2t1 * t2t1i21)
+      s12ks=0.5_wp*(s12s0+s21s0)
+
+   END FUNCTION s12ks
+
+   FUNCTION s22ks(x,y,z,phi) 
+      !!-------------------------------------------------------------------
+      !! Function : s22ks
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   x,y,z,phi
+
+      REAL(wp) ::   s22ks, p, pih
+
+      REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
+      REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
+      REAL(wp) ::   t1t2i11, t1t2i12, t1t2i21, t1t2i22
+      REAL(wp) ::   t2t1i11, t2t1i12, t2t1i21, t2t1i22
+      REAL(wp) ::   d11, d12, d22
+      REAL(wp) ::   IIn1t2, IIn2t1, IIt1t2
+      REAL(wp) ::   Hen1t2, Hen2t1
+      !!-------------------------------------------------------------------
+      pih = 0.5_wp*rpi
+      p =phi
+
+      n1t2i11 = cos(z+pih-p) * cos(z+p)
+      n1t2i12 = cos(z+pih-p) * sin(z+p)
+      n1t2i21 = sin(z+pih-p) * cos(z+p)
+      n1t2i22 = sin(z+pih-p) * sin(z+p)
+      n2t1i11 = cos(z-pih+p) * cos(z-p)
+      n2t1i12 = cos(z-pih+p) * sin(z-p)
+      n2t1i21 = sin(z-pih+p) * cos(z-p)
+      n2t1i22 = sin(z-pih+p) * sin(z-p)
+      t1t2i11 = cos(z-p) * cos(z+p)
+      t1t2i12 = cos(z-p) * sin(z+p)
+      t1t2i21 = sin(z-p) * cos(z+p)
+      t1t2i22 = sin(z-p) * sin(z+p)
+      t2t1i11 = cos(z+p) * cos(z-p)
+      t2t1i12 = cos(z+p) * sin(z-p)
+      t2t1i21 = sin(z+p) * cos(z-p)
+      t2t1i22 = sin(z+p) * sin(z-p)
+      d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))
+      d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))
+      d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))
+      IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22
+      IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22
+      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
+
+      IF (-IIn1t2>=rsmall) THEN
+      Hen1t2 = 1._wp
+      ELSE
+      Hen1t2 = 0._wp
+      ENDIF
+
+      IF (-IIn2t1>=rsmall) THEN
+      Hen2t1 = 1._wp
+      ELSE
+      Hen2t1 = 0._wp
+      ENDIF
+
+      s22ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i22 + Hen2t1 * t2t1i22)
+
+   END FUNCTION s22ks
+
 #else
    !!----------------------------------------------------------------------
    !!   Default option         Empty module           NO SI3 sea-ice model
    !!----------------------------------------------------------------------
+
 #endif
-
    !!==============================================================================
 END MODULE icedyn_rhg_eap
