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 13127)
+++ /NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90	(revision 13128)
@@ -142,9 +142,9 @@
       REAL(wp) ::   zfac_x, zfac_y
       REAL(wp) ::   zshear, zdum1, zdum2
-      REAL(wp) ::   stressptmp, stressmtmp, stress12tmpF                ! anisotropic stress tensor components
-      REAL(wp) ::   alphar, alphas                                      ! for mechanical redistribution
+      REAL(wp) ::   zstressptmp, zstressmtmp, zstress12tmpF                ! anisotropic stress tensor components
+      REAL(wp) ::   zalphar, zalphas                                      ! for mechanical redistribution
       REAL(wp) ::   mresult11, mresult12, z1dtevpkth, p5kth, z1_dtevp_A        ! for structure tensor evolution
       !
-      REAL(wp), DIMENSION(jpi,jpj) ::   stress12tmp                     ! anisotropic stress tensor component for regridding
+      REAL(wp), DIMENSION(jpi,jpj) ::   zstress12tmp                     ! anisotropic stress tensor component for regridding
       REAL(wp), DIMENSION(jpi,jpj) ::   yield11, yield22, yield12       ! yield surface tensor for history
       REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
@@ -415,11 +415,11 @@
                CALL update_stress_rdg (jter, nn_nevp, phi, zdiv, zdt, &
                                        zdsT, paniso_11(ji,jj), paniso_12(ji,jj), &
-                                       stressptmp, stressmtmp, &
-                                       stress12tmp(ji,jj), strength(ji,jj), &
-                                       alphar, alphas)
+                                       zstressptmp, zstressmtmp, &
+                                       zstress12tmp(ji,jj), strength(ji,jj), &
+                                       zalphar, zalphas)
 
                ! structure tensor update
                IF (mod(jter,10) == 0) THEN
-                  CALL calc_ffrac(stressptmp, stressmtmp, stress12tmp(ji,jj), &
+                  CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), &
                                   paniso_11(ji,jj), paniso_12(ji,jj),           &
                                   mresult11,  mresult12)
@@ -431,8 +431,8 @@
 
                IF (jter == nn_nevp) THEN
-                  yield11(ji,jj) = 0.5_wp * (stressptmp + stressmtmp)
-                  yield22(ji,jj) = 0.5_wp * (stressptmp - stressmtmp)
-                  yield12(ji,jj) = stress12tmp(ji,jj)
-                  prdg_conv(ji,jj) = -min( alphar, 0._wp)    
+                  yield11(ji,jj) = 0.5_wp * (zstressptmp + zstressmtmp)
+                  yield22(ji,jj) = 0.5_wp * (zstressptmp - zstressmtmp)
+                  yield12(ji,jj) = zstress12tmp(ji,jj)
+                  prdg_conv(ji,jj) = -min( zalphar, 0._wp)    
                ENDIF
 
@@ -454,6 +454,6 @@
                !zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1
                !zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2
-               zs1(ji,jj) = ( zs1(ji,jj) + stressptmp * zalph1  ) * z1_alph1
-               zs2(ji,jj) = ( zs2(ji,jj) + stressmtmp * zalph1  ) * z1_alph1
+               zs1(ji,jj) = ( zs1(ji,jj) + zstressptmp * zalph1  ) * z1_alph1
+               zs2(ji,jj) = ( zs2(ji,jj) + zstressmtmp * zalph1  ) * z1_alph1
                !zs1(ji,jj) = ( stressptmp * zs1(ji,jj) + zalph1  ) * z1_alph1
                !zs2(ji,jj) = ( stressmtmp * zs2(ji,jj) + zalph1  ) * z1_alph1
@@ -462,11 +462,11 @@
          END DO
          !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. )
-         CALL lbc_lnk_multi( 'icedyn_rhg_eap', stress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.)
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', zstress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.)
 
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
                ! stress12tmp at F points 
-               stress12tmpF = ( stress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + stress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  &
-                  &   + stress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + stress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  &
+               zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  &
+                  &   + zstress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + zstress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  &
                   &   ) * 0.25_wp * r1_e1e2f(ji,jj)
 
@@ -479,5 +479,5 @@
                ! stress at F points (zkt/=0 if landfast)
                !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
-               zs12(ji,jj) = ( zs12(ji,jj) + stress12tmpF * zalph1  ) * z1_alph1
+               zs12(ji,jj) = ( zs12(ji,jj) + zstress12tmpF * zalph1  ) * z1_alph1
                !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1  ) * z1_alph1
 
@@ -527,5 +527,5 @@
                   zvel  = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) )
                   zTauB = ztauy_base(ji,jj) / zvel
-                  !                 !--- OceanBottom-to-Ice stress
+                  !                 !--- OceanBottom-to-Ice stress 
                   ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj)
                   !
@@ -572,5 +572,5 @@
                   zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
                      &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
-                  !                 !--- Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress 
                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
                   !
@@ -625,5 +625,5 @@
                   zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  &
                      &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) )
-                  !                 !--- Ocean-to-Ice stress
+                  !                 !--- Ocean-to-Ice stress 
                   ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) )
                   !
@@ -841,5 +841,5 @@
 !!                                                                                                               ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11))
                zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) )          ! compressive stress, see Bouillon et al. 2015
-               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress
+               zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear )                     ! shear stress 
                zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 )
             END DO
@@ -853,5 +853,5 @@
          ! Stress tensor invariants (normal and shear stress N/m)
          IF( iom_use('normstr') )   CALL iom_put( 'normstr' ,       ( zs1(:,:) + zs2(:,:) )                       * zmsk00(:,:) ) ! Normal stress
-         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress
+         IF( iom_use('sheastr') )   CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 
 
          DEALLOCATE( zsig1 , zsig2 , zsig3 )
@@ -935,9 +935,9 @@
                                  stressp,  stressm, &
                                  stress12, strength, &
-                                 alphar, alphas)
+                                 palphar, palphas)
       !!---------------------------------------------------------------------
       !!                   ***  SUBROUTINE rhg_eap_rst  ***
       !!                     
-      !! ** Purpose :   Updates the stress depending on values of strain rate and structure
+      !! ** Purpose :   Updates the  depending on values of strain rate and structure
       !!                tensor and for ksub=ndte it computes closing and sliding rate
       !!---------------------------------------------------------------------
@@ -947,5 +947,5 @@
       REAL(wp), INTENT(in   ) ::   a11, a12
       REAL(wp), INTENT(  out) ::   stressp, stressm, stress12
-      REAL(wp), INTENT(  out) ::   alphar, alphas
+      REAL(wp), INTENT(  out) ::   palphar, palphas
 
       INTEGER  ::   kx ,ky, ka
@@ -964,6 +964,6 @@
       REAL(wp) ::   zTany_1, zTany_2 
       REAL(wp) ::   zx, zy, zdx, zdy, zda, zkxw, kyw, kaw
-      REAL(wp) ::   invdx, invdy, invda   
-      REAL(wp) ::   dtemp1, dtemp2, atempprime, a, invsin
+      REAL(wp) ::   zinvdx, zinvdy, zinvda   
+      REAL(wp) ::   zdtemp1, zdtemp2, zatempprime, zinvsin
 
       REAL(wp), PARAMETER ::   kfriction = 0.45_wp
@@ -973,5 +973,5 @@
       zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction)
  
-      invsin = 1._wp/sin(2._wp*phi) * zinvstressconviso 
+      zinvsin = 1._wp/sin(2._wp*phi) * zinvstressconviso 
       !now uses phi as set in higher code
        
@@ -997,8 +997,8 @@
  
       ! rotation Q*atemp*Q^T
-      atempprime = zQ11Q11*a11 + 2.0_wp*zQ11Q12*a12 + zQ12Q12*za22 
+      zatempprime = zQ11Q11*a11 + 2.0_wp*zQ11Q12*a12 + zQ12Q12*za22 
       
       ! make first principal value the largest
-      atempprime = max(atempprime, 1.0_wp - atempprime)
+      zatempprime = max(zatempprime, 1.0_wp - zatempprime)
  
       ! 2) strain rate
@@ -1025,10 +1025,10 @@
       ENDIF
 
-      dtemp1 = zQd11Qd11*zdtemp11 + 2.0_wp*zQd11Qd12*zdtemp12 + zQd12Qd12*zdtemp22
-      dtemp2 = zQd12Qd12*zdtemp11 - 2.0_wp*zQd11Qd12*zdtemp12 + zQd11Qd11*zdtemp22
+      zdtemp1 = zQd11Qd11*zdtemp11 + 2.0_wp*zQd11Qd12*zdtemp12 + zQd12Qd12*zdtemp22
+      zdtemp2 = zQd12Qd12*zdtemp11 - 2.0_wp*zQd11Qd12*zdtemp12 + zQd11Qd11*zdtemp22
       ! In cos and sin values
       zx = 0._wp
-      IF ((ABS(dtemp1) > rsmall).OR.(ABS(dtemp2) > rsmall)) THEN
-         zx = atan2(dtemp2,dtemp1)
+      IF ((ABS(zdtemp1) > rsmall).OR.(ABS(zdtemp2) > rsmall)) THEN
+         zx = atan2(zdtemp2,zdtemp1)
       ENDIF      
                
@@ -1059,18 +1059,18 @@
       zdy   = rpi/real(ny_yield-1,kind=wp)
       zda   = 0.5_wp/real(na_yield-1,kind=wp)
-      invdx = 1._wp/zdx
-      invdy = 1._wp/zdy
-      invda = 1._wp/zda
+      zinvdx = 1._wp/zdx
+      zinvdy = 1._wp/zdy
+      zinvda = 1._wp/zda
 
       ! % need 8 coords and 8 weights
       ! % range in kx
-      kx  = int((zx-rpi*0.25_wp-rpi)*invdx) + 1
-      zkxw = kx - (zx-rpi*0.25_wp-rpi)*invdx 
+      kx  = int((zx-rpi*0.25_wp-rpi)*zinvdx) + 1
+      zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx 
       
-      ky  = int(zy*invdy) + 1
-      kyw = ky - zy*invdy 
+      ky  = int(zy*zinvdy) + 1
+      kyw = ky - zy*zinvdy 
       
-      ka  = int((atempprime-0.5_wp)*invda) + 1
-      kaw = ka - (atempprime-0.5_wp)*invda
+      ka  = 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)
@@ -1128,7 +1128,7 @@
       ! Tsamados 2013)
 
-      zsig11  = strength*(zstemp11r + kfriction*zstemp11s) * invsin
-      zsig12  = strength*(zstemp12r + kfriction*zstemp12s) * invsin
-      zsig22  = strength*(zstemp22r + kfriction*zstemp22s) * invsin
+      zsig11  = strength*(zstemp11r + kfriction*zstemp11s) * zinvsin
+      zsig12  = strength*(zstemp12r + kfriction*zstemp12s) * zinvsin
+      zsig22  = strength*(zstemp22r + kfriction*zstemp22s) * zinvsin
 
       !WRITE(numout,*) 'strength inside loop', strength
@@ -1146,5 +1146,5 @@
       stressm  = zsgprm11 - zsgprm22
 
-      !WRITE(numout,*) 'stress output inside loop', ksub, stressp
+      !WRITE(numout,*) ' output inside loop', ksub, stressp
 
       ! Compute ridging and sliding functions in general coordinates 
@@ -1165,7 +1165,7 @@
                      + zQ11Q11*zstemp22s
 
-         alphar =  zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 &
+         palphar =  zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 &
                  + zrotstemp22r*zdtemp22
-         alphas =  zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 &
+         palphas =  zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 &
                  + zrotstemp22s*zdtemp22
       ENDIF
@@ -1190,8 +1190,8 @@
       ! local variables
 
-      REAL (wp) ::   sigma11, sigma12, sigma22  ! stress tensor elements
-      REAL (wp) ::   Angle_denom        ! angle with principal component axis
-      REAL (wp) ::   sigma_1, sigma_2   ! principal components of stress
-      REAL (wp) ::   Q11, Q12, Q11Q11, Q11Q12, Q12Q12
+      REAL (wp) ::   zsigma11, zsigma12, zsigma22  ! stress tensor elements
+      REAL (wp) ::   zAngle_denom        ! angle with principal component axis
+      REAL (wp) ::   zsigma_1, zsigma_2   ! principal components of stress
+      REAL (wp) ::   zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12
 
       REAL (wp), PARAMETER ::   kfrac = 0.0001_wp   ! rate of fracture formation 
@@ -1199,7 +1199,7 @@
       !!---------------------------------------------------------------
       !
-      sigma11 = 0.5_wp*(stressp+stressm) 
-      sigma12 = stress12 
-      sigma22 = 0.5_wp*(stressp-stressm) 
+      zsigma11 = 0.5_wp*(stressp+stressm) 
+      zsigma12 = stress12 
+      zsigma22 = 0.5_wp*(stressp-stressm) 
 
       ! Here's the change - no longer calculate gamma,
@@ -1212,25 +1212,25 @@
       ! error to the calculated angles < rsmall
 
-      Q11Q11 = 0.1_wp
-      Q12Q12 = rsmall
-      Q11Q12 = rsmall
-
-      IF((ABS( sigma11 - sigma22) > rsmall).OR.(ABS(sigma12) > rsmall)) THEN
-
-         Angle_denom = 1.0_wp/sqrt( ( sigma11 - sigma22 )*( sigma11 - &
-                       sigma22 ) + 4.0_wp*sigma12*sigma12)
-
-         Q11Q11 = 0.5_wp + ( sigma11 - sigma22 )*0.5_wp*Angle_denom   !Cos^2 
-         Q12Q12 = 0.5_wp - ( sigma11 - sigma22 )*0.5_wp*Angle_denom   !Sin^2
-         Q11Q12 = sigma12*Angle_denom                      !CosSin 
-      ENDIF
-
-      sigma_1 = Q11Q11*sigma11 + 2.0_wp*Q11Q12*sigma12  &
-              + Q12Q12*sigma22 ! S(1,1)
-      sigma_2 = Q12Q12*sigma11 - 2.0_wp*Q11Q12*sigma12  &
-              + Q11Q11*sigma22 ! S(2,2)
+      zQ11Q11 = 0.1_wp
+      zQ12Q12 = rsmall
+      zQ11Q12 = rsmall
+
+      IF((ABS( zsigma11 - zsigma22) > rsmall).OR.(ABS(zsigma12) > rsmall)) THEN
+
+         zAngle_denom = 1.0_wp/sqrt( ( zsigma11 - zsigma22 )*( zsigma11 - &
+                       zsigma22 ) + 4.0_wp*zsigma12*zsigma12)
+
+         zQ11Q11 = 0.5_wp + ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Cos^2 
+         zQ12Q12 = 0.5_wp - ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom   !Sin^2
+         zQ11Q12 = zsigma12*zAngle_denom                      !CosSin 
+      ENDIF
+
+      zsigma_1 = zQ11Q11*zsigma11 + 2.0_wp*zQ11Q12*zsigma12  &
+              + zQ12Q12*zsigma22 ! S(1,1)
+      zsigma_2 = zQ12Q12*zsigma11 - 2.0_wp*zQ11Q12*zsigma12  &
+              + zQ11Q11*zsigma22 ! S(2,2)
 
       ! Pure divergence
-      IF ((sigma_1 >= 0.0_wp).AND.(sigma_2 >= 0.0_wp))  THEN
+      IF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 >= 0.0_wp))  THEN
          mresult11 = 0.0_wp
          mresult12 = 0.0_wp
@@ -1239,15 +1239,15 @@
       ! direction
       ! which leads to the loss of their shape, so we again model it through diffusion
-      ELSEIF ((sigma_1 >= 0.0_wp).AND.(sigma_2 < 0.0_wp))  THEN
-         mresult11 = kfrac * (a11 - Q12Q12)
-         mresult12 = kfrac * (a12 + Q11Q12)
+      ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp))  THEN
+         mresult11 = kfrac * (a11 - zQ12Q12)
+         mresult12 = kfrac * (a12 + zQ11Q12)
 
       ! Shear faulting
-      ELSEIF (sigma_2 == 0.0_wp) THEN
+      ELSEIF (zsigma_2 == 0.0_wp) THEN
          mresult11 = 0.0_wp
          mresult12 = 0.0_wp
-      ELSEIF ((sigma_1 <= 0.0_wp).AND.(sigma_1/sigma_2 <= threshold)) THEN
-         mresult11 = kfrac * (a11 - Q12Q12)
-         mresult12 = kfrac * (a12 + Q11Q12)
+      ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN
+         mresult11 = kfrac * (a11 - zQ12Q12)
+         mresult12 = kfrac * (a12 + zQ11Q12)
 
       ! Horizontal spalling
@@ -1398,51 +1398,344 @@
 
 
-   FUNCTION w1(a)
+   FUNCTION w1(pa)
       !!-------------------------------------------------------------------
       !! Function : w1 (see Gaussian function psi in Tsamados et al 2013)
       !!-------------------------------------------------------------------
-      REAL(wp), INTENT(in   ) ::   a
+      REAL(wp), INTENT(in   ) ::   pa
       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
+       &   + 2361.2198663_wp*pa &
+       &   - 10606.56079975_wp*pa*pa &
+       &   + 26315.50025642_wp*pa*pa*pa &
+       &   - 38948.30444297_wp*pa*pa*pa*pa &
+       &   + 34397.72407466_wp*pa*pa*pa*pa*pa &
+       &   - 16789.98003081_wp*pa*pa*pa*pa*pa*pa &
+       &   + 3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa
    
    END FUNCTION w1
 
 
-   FUNCTION w2(a)
+   FUNCTION w2(pa)
       !!-------------------------------------------------------------------
       !! Function : w2 (see Gaussian function psi in Tsamados et al 2013)
       !!-------------------------------------------------------------------
-      REAL(wp), INTENT(in   ) ::   a
+      REAL(wp), INTENT(in   ) ::   pa
       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
+       &   + 70222.33061536_wp*pa &
+       &   - 314871.71525448_wp*pa*pa &
+       &   + 779570.02793492_wp*pa*pa*pa &
+       &   - 1151098.82436864_wp*pa*pa*pa*pa &
+       &   + 1013896.59464498_wp*pa*pa*pa*pa*pa &
+       &   - 493379.44906738_wp*pa*pa*pa*pa*pa*pa &
+       &   + 102356.551518_wp*pa*pa*pa*pa*pa*pa*pa
    
    END FUNCTION w2
 
-   FUNCTION s11kr(x,y,z,phi) 
+   FUNCTION s11kr(px,py,pz,phi) 
       !!-------------------------------------------------------------------
       !! Function : s11kr
       !!-------------------------------------------------------------------
-      REAL(wp), INTENT(in   ) ::   x,y,z,phi
-   
-      REAL(wp) ::   s11kr, p, pih
-   
+      REAL(wp), INTENT(in   ) ::   px,py,pz,phi
+   
+      REAL(wp) ::   s11kr, zpih
+   
+      REAL(wp) ::   zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22
+      REAL(wp) ::   zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22
+      REAL(wp) ::   zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22
+      REAL(wp) ::   zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22
+      REAL(wp) ::   zd11, zd12, zd22
+      REAL(wp) ::   zIIn1t2, zIIn2t1, zIIt1t2
+      REAL(wp) ::   zHen1t2, zHen2t1
+      !!-------------------------------------------------------------------
+   
+      zpih = 0.5_wp*rpi
+   
+      zn1t2i11 = cos(pz+zpih-phi) * cos(pz+phi)
+      zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi)
+      zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi)
+      zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi)
+      zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi)
+      zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi)
+      zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi)
+      zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi)
+      zt1t2i11 = cos(pz-phi) * cos(pz+phi)
+      zt1t2i12 = cos(pz-phi) * sin(pz+phi)
+      zt1t2i21 = sin(pz-phi) * cos(pz+phi)
+      zt1t2i22 = sin(pz-phi) * sin(pz+phi)
+      zt2t1i11 = cos(pz+phi) * cos(pz-phi)
+      zt2t1i12 = cos(pz+phi) * sin(pz-phi)
+      zt2t1i21 = sin(pz+phi) * cos(pz-phi)
+      zt2t1i22 = sin(pz+phi) * sin(pz-phi)
+   ! In expression of tensor d, with this formulatin d(x)=-d(x+pi)
+   ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else
+   ! x=x-pi/2
+      zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py))
+      zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px))
+      zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py))
+      zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22
+      zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22
+      zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22
+   
+      IF (-zIIn1t2>=rsmall) THEN
+      zHen1t2 = 1._wp
+      ELSE
+      zHen1t2 = 0._wp
+      ENDIF
+   
+      IF (-zIIn2t1>=rsmall) THEN
+      zHen2t1 = 1._wp
+      ELSE
+      zHen2t1 = 0._wp
+      ENDIF
+   
+      s11kr = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11)
+
+   END FUNCTION s11kr
+
+   FUNCTION s12kr(px,py,pz,phi)
+      !!-------------------------------------------------------------------
+      !! Function : s12kr
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   px,py,pz,phi
+
+      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-phi) * cos(pz+phi)
+      zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi)
+      zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi)
+      zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi)
+      zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi)
+      zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi)
+      zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi)
+      zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi)
+      zt1t2i11 = cos(pz-phi) * cos(pz+phi)
+      zt1t2i12 = cos(pz-phi) * sin(pz+phi)
+      zt1t2i21 = sin(pz-phi) * cos(pz+phi)
+      zt1t2i22 = sin(pz-phi) * sin(pz+phi)
+      zt2t1i11 = cos(pz+phi) * cos(pz-phi)
+      zt2t1i12 = cos(pz+phi) * sin(pz-phi)
+      zt2t1i21 = sin(pz+phi) * cos(pz-phi)
+      zt2t1i22 = sin(pz+phi) * sin(pz-phi)
+      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,phi)
+      !!-------------------------------------------------------------------
+      !! Function : s22kr
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   px,py,pz,phi
+
+      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-phi) * cos(pz+phi)
+      zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi)
+      zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi)
+      zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi)
+      zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi)
+      zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi)
+      zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi)
+      zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi)
+      zt1t2i11 = cos(pz-phi) * cos(pz+phi)
+      zt1t2i12 = cos(pz-phi) * sin(pz+phi)
+      zt1t2i21 = sin(pz-phi) * cos(pz+phi)
+      zt1t2i22 = sin(pz-phi) * sin(pz+phi)
+      zt2t1i11 = cos(pz+phi) * cos(pz-phi)
+      zt2t1i12 = cos(pz+phi) * sin(pz-phi)
+      zt2t1i21 = sin(pz+phi) * cos(pz-phi)
+      zt2t1i22 = sin(pz+phi) * sin(pz-phi)
+      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,phi)
+      !!-------------------------------------------------------------------
+      !! Function : s11ks
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   px,py,pz,phi
+
+      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-phi) * cos(pz+phi)
+      zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi)
+      zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi)
+      zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi)
+      zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi)
+      zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi)
+      zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi)
+      zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi)
+      zt1t2i11 = cos(pz-phi) * cos(pz+phi)
+      zt1t2i12 = cos(pz-phi) * sin(pz+phi)
+      zt1t2i21 = sin(pz-phi) * cos(pz+phi)
+      zt1t2i22 = sin(pz-phi) * sin(pz+phi)
+      zt2t1i11 = cos(pz+phi) * cos(pz-phi)
+      zt2t1i12 = cos(pz+phi) * sin(pz-phi)
+      zt2t1i21 = sin(pz+phi) * cos(pz-phi)
+      zt2t1i22 = sin(pz+phi) * sin(pz-phi)
+      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,phi)
+      !!-------------------------------------------------------------------
+      !! Function : s12ks
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   px,py,pz,phi
+
+      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-phi) * cos(pz+phi)
+      zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi)
+      zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi)
+      zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi)
+      zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi)
+      zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi)
+      zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi)
+      zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi)
+      zt1t2i11 = cos(pz-phi) * cos(pz+phi)
+      zt1t2i12 = cos(pz-phi) * sin(pz+phi)
+      zt1t2i21 = sin(pz-phi) * cos(pz+phi)
+      zt1t2i22 = sin(pz-phi) * sin(pz+phi)
+      zt2t1i11 = cos(pz+phi) * cos(pz-phi)
+      zt2t1i12 = cos(pz+phi) * sin(pz-phi)
+      zt2t1i21 = sin(pz+phi) * cos(pz-phi)
+      zt2t1i22 = sin(pz+phi) * sin(pz-phi)
+      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,phi) 
+      !!-------------------------------------------------------------------
+      !! Function : s22ks
+      !!-------------------------------------------------------------------
+      REAL(wp), INTENT(in   ) ::   px,py,pz,phi
+
+      REAL(wp) ::   s22ks, zpih
+
       REAL(wp) ::   n1t2i11, n1t2i12, n1t2i21, n1t2i22
       REAL(wp) ::   n2t1i11, n2t1i12, n2t1i21, n2t1i22
@@ -1453,34 +1746,29 @@
       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))
+      zpih = 0.5_wp*rpi
+
+      n1t2i11 = cos(pz+zpih-phi) * cos(pz+phi)
+      n1t2i12 = cos(pz+zpih-phi) * sin(pz+phi)
+      n1t2i21 = sin(pz+zpih-phi) * cos(pz+phi)
+      n1t2i22 = sin(pz+zpih-phi) * sin(pz+phi)
+      n2t1i11 = cos(pz-zpih+phi) * cos(pz-phi)
+      n2t1i12 = cos(pz-zpih+phi) * sin(pz-phi)
+      n2t1i21 = sin(pz-zpih+phi) * cos(pz-phi)
+      n2t1i22 = sin(pz-zpih+phi) * sin(pz-phi)
+      t1t2i11 = cos(pz-phi) * cos(pz+phi)
+      t1t2i12 = cos(pz-phi) * sin(pz+phi)
+      t1t2i21 = sin(pz-phi) * cos(pz+phi)
+      t1t2i22 = sin(pz-phi) * sin(pz+phi)
+      t2t1i11 = cos(pz+phi) * cos(pz-phi)
+      t2t1i12 = cos(pz+phi) * sin(pz-phi)
+      t2t1i21 = sin(pz+phi) * cos(pz-phi)
+      t2t1i22 = sin(pz+phi) * sin(pz-phi)
+      d11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py))
+      d12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px))
+      d22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py))
       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
@@ -1488,5 +1776,5 @@
       Hen1t2 = 0._wp
       ENDIF
-   
+
       IF (-IIn2t1>=rsmall) THEN
       Hen2t1 = 1._wp
@@ -1494,298 +1782,4 @@
       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)
