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 13105)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90	(revision 13113)
@@ -144,4 +144,5 @@
       REAL(wp) ::   stressptmp, stressmtmp, stress12tmpF                ! anisotropic stress tensor components
       REAL(wp) ::   alphar, alphas                                      ! 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
@@ -253,4 +254,9 @@
       zs2 (:,:) = pstress2_i (:,:)
       zs12(:,:) = pstress12_i(:,:)
+
+      ! constants for structure tensor
+      z1_dtevp_A = z1_dtevp/10.0_wp
+      z1dtevpkth = 1._wp / (z1_dtevp_A + 0.00002_wp)
+      p5kth = 0.5_wp * 0.00002_wp
 
       ! Ice strength
@@ -412,4 +418,16 @@
                                        stress12tmp(ji,jj), strength(ji,jj), &
                                        alphar, alphas)
+
+               ! structure tensor update
+               IF (mod(jter,10) == 0) THEN
+                  CALL calc_ffrac(stressptmp, stressmtmp, stress12tmp(ji,jj), &
+                                  paniso_11(ji,jj), paniso_12(ji,jj),           &
+                                  mresult11,  mresult12)
+
+                  paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A  + p5kth - mresult11) * z1dtevpkth ! implicit
+                  paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A  - mresult12) * z1dtevpkth ! implicit
+               ENDIF
+
+
                IF (jter == nn_nevp) THEN
                   yield11(ji,jj) = 0.5_wp * (stressptmp + stressmtmp)
@@ -428,8 +446,8 @@
                !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
                !   alpha = beta = sqrt(4*gamma)
-               IF( ln_aEVP ) THEN
-                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
-                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
-               ENDIF
+               !IF( ln_aEVP ) THEN
+               !   zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
+               !   z1_alph1 = 1._wp / ( zalph1 + 1._wp )
+               !ENDIF
                
                ! stress at T points (zkt/=0 if landfast)
@@ -444,5 +462,5 @@
          END DO
          !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. )
-         CALL lbc_lnk( 'icedyn_rhg_eap', stress12tmp, 'T', 1. )
+         CALL lbc_lnk_multi( 'icedyn_rhg_eap', stress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.)
 
          DO jj = 1, jpjm1
@@ -714,4 +732,5 @@
 !!$         ENDIF
          !
+
          !                                                ! ==================== !
       END DO                                              !  end loop over jter  !
@@ -851,6 +870,6 @@
       ! --- anisotropy tensor --- !
       IF( iom_use('aniso') ) THEN       
-         CALL lbc_lnk( 'icedyn_rhg_eap', aniso_11, 'T', 1. ) 
-         CALL iom_put( 'aniso' , aniso_11 * zmsk00 )
+         CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1. ) 
+         CALL iom_put( 'aniso' , paniso_11 * zmsk00 )
       ENDIF
       
@@ -919,6 +938,7 @@
       !!---------------------------------------------------------------------
       !!                   ***  SUBROUTINE rhg_eap_rst  ***
-      !! Updates the stress depending on values of strain rate and structure
-      !! tensor and for ksub=ndte it computes closing and sliding rate
+      !!                     
+      !! ** Purpose :   Updates the stress depending on values of strain rate and structure
+      !!                tensor and for ksub=ndte it computes closing and sliding rate
       !!---------------------------------------------------------------------
       INTEGER,  INTENT(in   ) ::   ksub, ndte
@@ -931,17 +951,17 @@
       INTEGER  ::   kx ,ky, ka
 
-      REAL(wp) ::   stemp11r, stemp12r, stemp22r   
-      REAL(wp) ::   stemp11s, stemp12s, stemp22s 
-      REAL(wp) ::   a22, Qd11Qd11, Qd11Qd12, Qd12Qd12 
-      REAL(wp) ::   Q11Q11, Q11Q12, Q12Q12 
-      REAL(wp) ::   dtemp11, dtemp12, dtemp22 
-      REAL(wp) ::   rotstemp11r, rotstemp12r, rotstemp22r   
-      REAL(wp) ::   rotstemp11s, rotstemp12s, rotstemp22s
-      REAL(wp) ::   sig11, sig12, sig22 
-      REAL(wp) ::   sgprm11, sgprm12, sgprm22 
-      REAL(wp) ::   invstressconviso 
-      REAL(wp) ::   Angle_denom_gamma,  Angle_denom_alpha 
-      REAL(wp) ::   Tany_1, Tany_2 
-      REAL(wp) ::   x, y, dx, dy, da, kxw, kyw, kaw
+      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) ::   invdx, invdy, invda   
       REAL(wp) ::   dtemp1, dtemp2, atempprime, a, invsin
@@ -951,7 +971,7 @@
       ! Factor to maintain the same stress as in EVP (see Section 3)
       ! Can be set to 1 otherwise
-      invstressconviso = 1._wp/(1._wp+kfriction*kfriction)
+      zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction)
  
-      invsin = 1._wp/sin(2._wp*phi) * invstressconviso 
+      invsin = 1._wp/sin(2._wp*phi) * zinvstressconviso 
       !now uses phi as set in higher code
        
@@ -960,22 +980,22 @@
 
       ! 1) structure tensor
-      a22 = 1._wp-a11
-      Q11Q11 = 1._wp
-      Q12Q12 = rsmall
-      Q11Q12 = rsmall
+      za22 = 1._wp-a11
+      zQ11Q11 = 1._wp
+      zQ12Q12 = rsmall
+      zQ11Q12 = rsmall
  
       ! gamma: angle between general coordiantes and principal axis of A 
       ! here Tan2gamma = 2 a12 / (a11 - a22) 
-      if((ABS(a11 - a22) > rsmall).or.(ABS(a12) > rsmall)) then      
-         Angle_denom_gamma = 1._wp/sqrt( ( a11 - a22 )*( a11 - a22) + &
+      IF((ABS(a11 - za22) > rsmall).OR.(ABS(a12) > rsmall)) THEN      
+         zAngle_denom_gamma = 1._wp/sqrt( ( a11 - za22 )*( a11 - za22) + &
                              4._wp*a12*a12 )
    
-         Q11Q11 = 0.5_wp + ( a11 - a22 )*0.5_wp*Angle_denom_gamma   !Cos^2 
-         Q12Q12 = 0.5_wp - ( a11 - a22 )*0.5_wp*Angle_denom_gamma   !Sin^2
-         Q11Q12 = a12*Angle_denom_gamma                     !CosSin 
-      endif
+         zQ11Q11 = 0.5_wp + ( a11 - za22 )*0.5_wp*zAngle_denom_gamma   !Cos^2 
+         zQ12Q12 = 0.5_wp - ( a11 - za22 )*0.5_wp*zAngle_denom_gamma   !Sin^2
+         zQ11Q12 = a12*zAngle_denom_gamma                     !CosSin 
+      ENDIF
  
       ! rotation Q*atemp*Q^T
-      atempprime = Q11Q11*a11 + 2.0_wp*Q11Q12*a12 + Q12Q12*a22 
+      atempprime = zQ11Q11*a11 + 2.0_wp*zQ11Q12*a12 + zQ12Q12*za22 
       
       ! make first principal value the largest
@@ -983,7 +1003,7 @@
  
       ! 2) strain rate
-      dtemp11 = 0.5_wp*(divu + tension)
-      dtemp12 = shear*0.5_wp
-      dtemp22 = 0.5_wp*(divu - tension)
+      zdtemp11 = 0.5_wp*(divu + tension)
+      zdtemp12 = shear*0.5_wp
+      zdtemp22 = 0.5_wp*(divu - tension)
 
       !WRITE(numout,*) 'div, tension, shear inside loop', ksub, divu, tension, shear
@@ -991,28 +1011,28 @@
       ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22) 
 
-      Qd11Qd11 = 1.0_wp
-      Qd12Qd12 = rsmall
-      Qd11Qd12 = rsmall
-
-      if((ABS( dtemp11 - dtemp22) > rsmall).or. (ABS(dtemp12) > rsmall)) then
-
-         Angle_denom_alpha = 1.0_wp/sqrt( ( dtemp11 - dtemp22 )* &
-                             ( dtemp11 - dtemp22 ) + 4.0_wp*dtemp12*dtemp12)
-
-         Qd11Qd11 = 0.5_wp + ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Cos^2 
-         Qd12Qd12 = 0.5_wp - ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Sin^2
-         Qd11Qd12 = dtemp12*Angle_denom_alpha !CosSin
-      endif
-
-      dtemp1 = Qd11Qd11*dtemp11 + 2.0_wp*Qd11Qd12*dtemp12 + Qd12Qd12*dtemp22
-      dtemp2 = Qd12Qd12*dtemp11 - 2.0_wp*Qd11Qd12*dtemp12 + Qd11Qd11*dtemp22
+      zQd11Qd11 = 1.0_wp
+      zQd12Qd12 = rsmall
+      zQd11Qd12 = rsmall
+
+      IF((ABS( zdtemp11 - zdtemp22) > rsmall).OR. (ABS(zdtemp12) > rsmall)) THEN
+
+         zAngle_denom_alpha = 1.0_wp/sqrt( ( zdtemp11 - zdtemp22 )* &
+                             ( zdtemp11 - zdtemp22 ) + 4.0_wp*zdtemp12*zdtemp12)
+
+         zQd11Qd11 = 0.5_wp + ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Cos^2 
+         zQd12Qd12 = 0.5_wp - ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Sin^2
+         zQd11Qd12 = zdtemp12*zAngle_denom_alpha !CosSin
+      ENDIF
+
+      dtemp1 = zQd11Qd11*zdtemp11 + 2.0_wp*zQd11Qd12*zdtemp12 + zQd12Qd12*zdtemp22
+      dtemp2 = zQd12Qd12*zdtemp11 - 2.0_wp*zQd11Qd12*zdtemp12 + zQd11Qd11*zdtemp22
       ! In cos and sin values
-      x = 0._wp
-      if ((ABS(dtemp1) > rsmall).or.(ABS(dtemp2) > rsmall)) then
-         x = atan2(dtemp2,dtemp1)
-      endif      
+      zx = 0._wp
+      IF ((ABS(dtemp1) > rsmall).OR.(ABS(dtemp2) > rsmall)) THEN
+         zx = atan2(dtemp2,dtemp1)
+      ENDIF      
                
       ! to ensure the angle lies between pi/4 and 9 pi/4 
-      if (x < rpi*0.25_wp) x = x + rpi*2.0_wp
+      IF (zx < rpi*0.25_wp) zx = zx + rpi*2.0_wp
       
       ! y: angle between major principal axis of strain rate and structure
@@ -1022,32 +1042,32 @@
       ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha)
       
-      Tany_1 = Q11Q12 - Qd11Qd12
-      Tany_2 = Q11Q11 - Qd12Qd12
+      zTany_1 = zQ11Q12 - zQd11Qd12
+      zTany_2 = zQ11Q11 - zQd12Qd12
       
-      y = 0._wp
+      zy = 0._wp
       
-      if ((ABS(Tany_1) > rsmall).or.(ABS(Tany_2) > rsmall)) then
-         y = atan2(Tany_1,Tany_2)
-      endif
+      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 (y > rpi) y = y - rpi
-      if (y < 0)  y = y + rpi
+      IF (zy > rpi) zy = zy - rpi
+      IF (zy < 0)  zy = zy + rpi
 
       ! 3) update anisotropic stress tensor 
-      dx   = rpi/real(nx_yield-1,kind=wp)
-      dy   = rpi/real(ny_yield-1,kind=wp)
-      da   = 0.5_wp/real(na_yield-1,kind=wp)
-      invdx = 1._wp/dx
-      invdy = 1._wp/dy
-      invda = 1._wp/da
+      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)
+      invdx = 1._wp/zdx
+      invdy = 1._wp/zdy
+      invda = 1._wp/zda
 
       ! % need 8 coords and 8 weights
       ! % range in kx
-      kx  = int((x-rpi*0.25_wp-rpi)*invdx) + 1
-      kxw = kx - (x-rpi*0.25_wp-rpi)*invdx 
+      kx  = int((zx-rpi*0.25_wp-rpi)*invdx) + 1
+      zkxw = kx - (zx-rpi*0.25_wp-rpi)*invdx 
       
-      ky  = int(y*invdy) + 1
-      kyw = ky - y*invdy 
+      ky  = int(zy*invdy) + 1
+      kyw = ky - zy*invdy 
       
       ka  = int((atempprime-0.5_wp)*invda) + 1
@@ -1055,60 +1075,60 @@
       
       ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013)
-      stemp11r =  kxw   * kyw         * kaw         * s11r(kx  ,ky  ,ka  ) &
-        & + (1._wp-kxw) * kyw         * kaw         * s11r(kx+1,ky  ,ka  ) &
-        & + kxw         * (1._wp-kyw) * kaw         * s11r(kx  ,ky+1,ka  ) &
-        & + kxw         * kyw         * (1._wp-kaw) * s11r(kx  ,ky  ,ka+1) &
-        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s11r(kx+1,ky+1,ka  ) &
-        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s11r(kx+1,ky  ,ka+1) &
-        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s11r(kx  ,ky+1,ka+1) &
-        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1)
-      stemp12r =  kxw   * kyw         * kaw         * s12r(kx  ,ky  ,ka  ) &
-        & + (1._wp-kxw) * kyw         * kaw         * s12r(kx+1,ky  ,ka  ) &
-        & + kxw         * (1._wp-kyw) * kaw         * s12r(kx  ,ky+1,ka  ) &
-        & + kxw         * kyw         * (1._wp-kaw) * s12r(kx  ,ky  ,ka+1) &
-        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s12r(kx+1,ky+1,ka  ) &
-        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s12r(kx+1,ky  ,ka+1) &
-        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s12r(kx  ,ky+1,ka+1) &
-        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1)
-      stemp22r = kxw    * kyw         * kaw         * s22r(kx  ,ky  ,ka  ) &
-        & + (1._wp-kxw) * kyw         * kaw         * s22r(kx+1,ky  ,ka  ) &
-        & + kxw         * (1._wp-kyw) * kaw         * s22r(kx  ,ky+1,ka  ) &
-        & + kxw         * kyw         * (1._wp-kaw) * s22r(kx  ,ky  ,ka+1) &
-        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s22r(kx+1,ky+1,ka  ) &
-        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s22r(kx+1,ky  ,ka+1) &
-        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s22r(kx  ,ky+1,ka+1) &
-        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1)
+      zstemp11r =  zkxw   * kyw         * kaw         * s11r(kx  ,ky  ,ka  ) &
+        & + (1._wp-zkxw) * kyw         * kaw         * s11r(kx+1,ky  ,ka  ) &
+        & + zkxw         * (1._wp-kyw) * kaw         * s11r(kx  ,ky+1,ka  ) &
+        & + zkxw         * kyw         * (1._wp-kaw) * s11r(kx  ,ky  ,ka+1) &
+        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s11r(kx+1,ky+1,ka  ) &
+        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s11r(kx+1,ky  ,ka+1) &
+        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s11r(kx  ,ky+1,ka+1) &
+        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1)
+      zstemp12r =  zkxw   * kyw         * kaw         * s12r(kx  ,ky  ,ka  ) &
+        & + (1._wp-zkxw) * kyw         * kaw         * s12r(kx+1,ky  ,ka  ) &
+        & + zkxw         * (1._wp-kyw) * kaw         * s12r(kx  ,ky+1,ka  ) &
+        & + zkxw         * kyw         * (1._wp-kaw) * s12r(kx  ,ky  ,ka+1) &
+        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s12r(kx+1,ky+1,ka  ) &
+        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s12r(kx+1,ky  ,ka+1) &
+        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s12r(kx  ,ky+1,ka+1) &
+        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1)
+      zstemp22r = zkxw    * kyw         * kaw         * s22r(kx  ,ky  ,ka  ) &
+        & + (1._wp-zkxw) * kyw         * kaw         * s22r(kx+1,ky  ,ka  ) &
+        & + zkxw         * (1._wp-kyw) * kaw         * s22r(kx  ,ky+1,ka  ) &
+        & + zkxw         * kyw         * (1._wp-kaw) * s22r(kx  ,ky  ,ka+1) &
+        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s22r(kx+1,ky+1,ka  ) &
+        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s22r(kx+1,ky  ,ka+1) &
+        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s22r(kx  ,ky+1,ka+1) &
+        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1)
       
-      stemp11s =  kxw   * kyw         * kaw         * s11s(kx  ,ky  ,ka  ) &
-        & + (1._wp-kxw) * kyw         * kaw         * s11s(kx+1,ky  ,ka  ) &
-        & + kxw         * (1._wp-kyw) * kaw         * s11s(kx  ,ky+1,ka  ) &
-        & + kxw         * kyw         * (1._wp-kaw) * s11s(kx  ,ky  ,ka+1) &
-        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s11s(kx+1,ky+1,ka  ) &
-        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s11s(kx+1,ky  ,ka+1) &
-        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s11s(kx  ,ky+1,ka+1) &
-        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1)
-      stemp12s =  kxw   * kyw         * kaw         * s12s(kx  ,ky  ,ka  ) &
-        & + (1._wp-kxw) * kyw         * kaw         * s12s(kx+1,ky  ,ka  ) &
-        & + kxw         * (1._wp-kyw) * kaw         * s12s(kx  ,ky+1,ka  ) &
-        & + kxw         * kyw         * (1._wp-kaw) * s12s(kx  ,ky  ,ka+1) &
-        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s12s(kx+1,ky+1,ka  ) &
-        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s12s(kx+1,ky  ,ka+1) &
-        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s12s(kx  ,ky+1,ka+1) &
-        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1)
-      stemp22s =  kxw   * kyw         * kaw         * s22s(kx  ,ky  ,ka  ) &
-        & + (1._wp-kxw) * kyw         * kaw         * s22s(kx+1,ky  ,ka  ) &
-        & + kxw         * (1._wp-kyw) * kaw         * s22s(kx  ,ky+1,ka  ) &
-        & + kxw         * kyw         * (1._wp-kaw) * s22s(kx  ,ky  ,ka+1) &
-        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s22s(kx+1,ky+1,ka  ) &
-        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s22s(kx+1,ky  ,ka+1) &
-        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s22s(kx  ,ky+1,ka+1) &
-        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1)
+      zstemp11s =  zkxw   * kyw         * kaw         * s11s(kx  ,ky  ,ka  ) &
+        & + (1._wp-zkxw) * kyw         * kaw         * s11s(kx+1,ky  ,ka  ) &
+        & + zkxw         * (1._wp-kyw) * kaw         * s11s(kx  ,ky+1,ka  ) &
+        & + zkxw         * kyw         * (1._wp-kaw) * s11s(kx  ,ky  ,ka+1) &
+        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s11s(kx+1,ky+1,ka  ) &
+        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s11s(kx+1,ky  ,ka+1) &
+        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s11s(kx  ,ky+1,ka+1) &
+        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1)
+      zstemp12s =  zkxw   * kyw         * kaw         * s12s(kx  ,ky  ,ka  ) &
+        & + (1._wp-zkxw) * kyw         * kaw         * s12s(kx+1,ky  ,ka  ) &
+        & + zkxw         * (1._wp-kyw) * kaw         * s12s(kx  ,ky+1,ka  ) &
+        & + zkxw         * kyw         * (1._wp-kaw) * s12s(kx  ,ky  ,ka+1) &
+        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s12s(kx+1,ky+1,ka  ) &
+        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s12s(kx+1,ky  ,ka+1) &
+        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s12s(kx  ,ky+1,ka+1) &
+        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1)
+      zstemp22s =  zkxw   * kyw         * kaw         * s22s(kx  ,ky  ,ka  ) &
+        & + (1._wp-zkxw) * kyw         * kaw         * s22s(kx+1,ky  ,ka  ) &
+        & + zkxw         * (1._wp-kyw) * kaw         * s22s(kx  ,ky+1,ka  ) &
+        & + zkxw         * kyw         * (1._wp-kaw) * s22s(kx  ,ky  ,ka+1) &
+        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s22s(kx+1,ky+1,ka  ) &
+        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s22s(kx+1,ky  ,ka+1) &
+        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s22s(kx  ,ky+1,ka+1) &
+        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1)
                    
       ! Calculate mean ice stress over a collection of floes (Equation 3 in
       ! Tsamados 2013)
 
-      sig11  = strength*(stemp11r + kfriction*stemp11s) * invsin
-      sig12  = strength*(stemp12r + kfriction*stemp12s) * invsin
-      sig22  = strength*(stemp22r + kfriction*stemp22s) * invsin
+      zsig11  = strength*(zstemp11r + kfriction*zstemp11s) * invsin
+      zsig12  = strength*(zstemp12r + kfriction*zstemp12s) * invsin
+      zsig22  = strength*(zstemp22r + kfriction*zstemp22s) * invsin
 
       !WRITE(numout,*) 'strength inside loop', strength
@@ -1118,11 +1138,11 @@
 
       ! Update stress
-      sgprm11 = Q11Q11*sig11 + Q12Q12*sig22 -        2._wp*Q11Q12 *sig12
-      sgprm12 = Q11Q12*sig11 - Q11Q12*sig22 + (Q11Q11 - Q12Q12)*sig12
-      sgprm22 = Q12Q12*sig11 + Q11Q11*sig22 +        2._wp*Q11Q12 *sig12
-
-      stressp  = sgprm11 + sgprm22
-      stress12 = sgprm12
-      stressm  = sgprm11 - sgprm22
+      zsgprm11 = zQ11Q11*zsig11 + zQ12Q12*zsig22 -        2._wp*zQ11Q12 *zsig12
+      zsgprm12 = zQ11Q12*zsig11 - zQ11Q12*zsig22 + (zQ11Q11 - zQ12Q12)*zsig12
+      zsgprm22 = zQ12Q12*zsig11 + zQ11Q11*zsig22 +        2._wp*zQ11Q12 *zsig12
+
+      stressp  = zsgprm11 + zsgprm22
+      stress12 = zsgprm12
+      stressm  = zsgprm11 - zsgprm22
 
       !WRITE(numout,*) 'stress output inside loop', ksub, stressp
@@ -1130,25 +1150,113 @@
       ! Compute ridging and sliding functions in general coordinates 
       ! (Equation 11 in Tsamados 2013)
-      if (ksub == ndte) then
-         rotstemp11r = Q11Q11*stemp11r - 2._wp*Q11Q12* stemp12r &
-                     + Q12Q12*stemp22r
-         rotstemp12r = Q11Q11*stemp12r +    Q11Q12*(stemp11r-stemp22r) &
-                     - Q12Q12*stemp12r
-         rotstemp22r = Q12Q12*stemp11r + 2._wp*Q11Q12* stemp12r & 
-                     + Q11Q11*stemp22r
-
-         rotstemp11s = Q11Q11*stemp11s - 2._wp*Q11Q12* stemp12s &
-                     + Q12Q12*stemp22s
-         rotstemp12s = Q11Q11*stemp12s +    Q11Q12*(stemp11s-stemp22s) &
-                     - Q12Q12*stemp12s
-         rotstemp22s = Q12Q12*stemp11s + 2._wp*Q11Q12* stemp12s & 
-                     + Q11Q11*stemp22s
-
-         alphar =  rotstemp11r*dtemp11 + 2._wp*rotstemp12r*dtemp12 &
-                 + rotstemp22r*dtemp22
-         alphas =  rotstemp11s*dtemp11 + 2._wp*rotstemp12s*dtemp12 &
-                 + rotstemp22s*dtemp22
-      endif
+      IF (ksub == ndte) THEN
+         zrotstemp11r = zQ11Q11*zstemp11r - 2._wp*zQ11Q12* zstemp12r &
+                     + zQ12Q12*zstemp22r
+         zrotstemp12r = zQ11Q11*zstemp12r +    zQ11Q12*(zstemp11r-zstemp22r) &
+                     - zQ12Q12*zstemp12r
+         zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r & 
+                     + zQ11Q11*zstemp22r
+
+         zrotstemp11s = zQ11Q11*zstemp11s - 2._wp*zQ11Q12* zstemp12s &
+                     + zQ12Q12*zstemp22s
+         zrotstemp12s = zQ11Q11*zstemp12s +    zQ11Q12*(zstemp11s-zstemp22s) &
+                     - zQ12Q12*zstemp12s
+         zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s & 
+                     + zQ11Q11*zstemp22s
+
+         alphar =  zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 &
+                 + zrotstemp22r*zdtemp22
+         alphas =  zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 &
+                 + zrotstemp22s*zdtemp22
+      ENDIF
    END SUBROUTINE update_stress_rdg 
+
+!=======================================================================
+
+
+   SUBROUTINE calc_ffrac (stressp, stressm, stress12, &
+                          a11, a12,                   &
+                          mresult11, mresult12)
+      !!---------------------------------------------------------------------
+      !!                   ***  ROUTINE calc_ffrac  ***
+      !!                     
+      !! ** Purpose :   Computes term in evolution equation for structure tensor
+      !!                which determines the ice floe re-orientation due to fracture
+      !! ** Method :    Eq. 7: Ffrac = -kf(A-S) or = 0 depending on sigma_1 and sigma_2
+      !!---------------------------------------------------------------------
+      REAL (wp), INTENT(in)  ::   stressp, stressm, stress12, a11, a12
+      REAL (wp), INTENT(out) ::   mresult11, mresult12
+
+      ! 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), PARAMETER ::   kfrac = 0.0001_wp   ! rate of fracture formation 
+      REAL (wp), PARAMETER ::   threshold = 0.3_wp  ! critical confinement ratio 
+      !!---------------------------------------------------------------
+      !
+      sigma11 = 0.5_wp*(stressp+stressm) 
+      sigma12 = stress12 
+      sigma22 = 0.5_wp*(stressp-stressm) 
+
+      ! Here's the change - no longer calculate gamma,
+      ! use trig formulation, angle signs are kept correct, don't worry
+   
+      ! rotate tensor to get into sigma principal axis
+   
+      ! here Tan2gamma = 2 sig12 / (sig11 - sig12) 
+      ! add rsmall to the denominator to stop 1/0 errors, makes very little 
+      ! error to the calculated angles < rsmall
+
+      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)
+
+      ! Pure divergence
+      IF ((sigma_1 >= 0.0_wp).AND.(sigma_2 >= 0.0_wp))  THEN
+         mresult11 = 0.0_wp
+         mresult12 = 0.0_wp
+
+      ! Unconfined compression: cracking of blocks not along the axial splitting
+      ! direction
+      ! which leads to the loss of their shape, so we again model it through diffusion
+      ELSEIF ((sigma_1 >= 0.0_wp).AND.(sigma_2 < 0.0_wp))  THEN
+         mresult11 = kfrac * (a11 - Q12Q12)
+         mresult12 = kfrac * (a12 + Q11Q12)
+
+      ! Shear faulting
+      ELSEIF (sigma_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)
+
+      ! Horizontal spalling
+      ELSE 
+         mresult11 = 0.0_wp
+         mresult12 = 0.0_wp
+      ENDIF
+
+   END SUBROUTINE calc_ffrac
+
 
    SUBROUTINE rhg_eap_rst( cdrw, kt )
@@ -1229,5 +1337,5 @@
          s12s(ix,iy,ia) = 0._wp
          s22s(ix,iy,ia) = 0._wp
-         IF (ia <= na_yield-1) then
+         IF (ia <= na_yield-1) THEN
           DO iz=1,nz
           s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* &
@@ -1299,11 +1407,11 @@
    
       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*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
@@ -1319,11 +1427,11 @@
    
       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*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
@@ -1375,5 +1483,5 @@
       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
    
-      IF (-IIn1t2>=rsmall) then
+      IF (-IIn1t2>=rsmall) THEN
       Hen1t2 = 1._wp
       ELSE
@@ -1381,5 +1489,5 @@
       ENDIF
    
-      IF (-IIn2t1>=rsmall) then
+      IF (-IIn2t1>=rsmall) THEN
       Hen2t1 = 1._wp
       ELSE
@@ -1433,5 +1541,5 @@
       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
    
-      IF (-IIn1t2>=rsmall) then
+      IF (-IIn1t2>=rsmall) THEN
       Hen1t2 = 1._wp
       ELSE
@@ -1439,5 +1547,5 @@
       ENDIF
    
-      IF (-IIn2t1>=rsmall) then
+      IF (-IIn2t1>=rsmall) THEN
       Hen2t1 = 1._wp
       ELSE
@@ -1493,5 +1601,5 @@
       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
 
-      IF (-IIn1t2>=rsmall) then
+      IF (-IIn1t2>=rsmall) THEN
       Hen1t2 = 1._wp
       ELSE
@@ -1499,5 +1607,5 @@
       ENDIF
 
-      IF (-IIn2t1>=rsmall) then
+      IF (-IIn2t1>=rsmall) THEN
       Hen2t1 = 1._wp
       ELSE
@@ -1551,5 +1659,5 @@
       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
 
-      IF (-IIn1t2>=rsmall) then
+      IF (-IIn1t2>=rsmall) THEN
       Hen1t2 = 1._wp
       ELSE
@@ -1557,5 +1665,5 @@
       ENDIF
 
-      IF (-IIn2t1>=rsmall) then
+      IF (-IIn2t1>=rsmall) THEN
       Hen2t1 = 1._wp
       ELSE
@@ -1609,5 +1717,5 @@
       IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22
 
-      IF (-IIn1t2>=rsmall) then
+      IF (-IIn1t2>=rsmall) THEN
       Hen1t2 = 1._wp
       ELSE
@@ -1615,5 +1723,5 @@
       ENDIF
 
-      IF (-IIn2t1>=rsmall) then
+      IF (-IIn2t1>=rsmall) THEN
       Hen2t1 = 1._wp
       ELSE
