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 11920)
+++ NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90	(revision 11933)
@@ -13,4 +13,6 @@
    !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017)
    !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube]
+   !!                 !  2019     (S. Rynders)  change into eap rheology from
+   !!                                           CICE code (Tsamados, Heorton) 
    !!----------------------------------------------------------------------
 #if defined key_si3
@@ -46,4 +48,6 @@
    PUBLIC   ice_dyn_rhg_eap   ! called by icedyn_rhg.F90
    PUBLIC   rhg_eap_rst       ! called by icedyn_rhg.F90
+
+   REAL(wp), PARAMETER           ::   phi = 3.141592653589793_wp/12._wp    ! diamond shaped floe smaller angle (default phi = 30 deg)
 
    ! look-up table for calculating structure tensor
@@ -126,7 +130,7 @@
       REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling
       REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity
-      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017
+      REAL(wp) ::   zalph1, z1_alph1                                    ! alpha coef from Bouillon 2009 or Kimmritz 2017
       REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume
-      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars
+      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2, zdsT ! temporary scalars
       REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars
       REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice
@@ -137,5 +141,8 @@
       REAL(wp) ::   zfac_x, zfac_y
       REAL(wp) ::   zshear, zdum1, zdum2
+      REAL(wp) ::   stressptmp, stressmtmp                              ! anisotropic stress tensor components
+      REAL(wp) ::   alphar, alphas                                      ! for mechanical redistribution
       !
+      REAL(wp), DIMENSION(jpi,jpj) ::   stress12tmp                     ! anisotropic stress tensor component for regridding
       REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points
       REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017
@@ -237,8 +244,5 @@
       IF( .NOT. ln_aEVP ) THEN
          zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp
-         zalph2 = zalph1 * z1_ecc2
-
          z1_alph1 = 1._wp / ( zalph1 + 1._wp )
-         z1_alph2 = 1._wp / ( zalph2 + 1._wp )
       ENDIF
          
@@ -376,4 +380,9 @@
             DO ji = 2, jpi ! no vector loop
 
+               ! shear at T points 
+               zdsT = ( zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
+                  &   + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  &
+                  &   ) * 0.25_wp * r1_e1e2t(ji,jj)
+
                ! shear**2 at T points (doc eq. A16)
                zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  &
@@ -393,43 +402,52 @@
                zdt2 = zdt * zdt
                
+               ! --- anisotropic stress calculation --- !
+               CALL update_stress_rdg (jter, nn_nevp, phi, zdiv, zdt, &
+                                       zdsT, aniso_11(ji,jj), aniso_12(ji,jj), &
+                                       stressptmp, stressmtmp, &
+                                       stress12tmp(ji,jj), strength(ji,jj), &
+                                       alphar, alphas)
+
                ! delta at T points
-               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
+               !zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )  
 
                ! P/delta at T points
-               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
+               !zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl )
 
                ! alpha & beta for aEVP
                !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m
                !   alpha = beta = sqrt(4*gamma)
-               IF( ln_aEVP ) THEN
-                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
-                  z1_alph1 = 1._wp / ( zalph1 + 1._wp )
-                  zalph2   = zalph1
-                  z1_alph2 = z1_alph1
-               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)
-               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) * 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
              
             END DO
          END DO
-         CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. )
+         !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. )
+         CALL lbc_lnk( 'icedyn_rhg_eap', stress12tmp, 'T', 1. )
 
          DO jj = 1, jpjm1
             DO ji = 1, jpim1
+               ! stress12tmp at F points 
+               zdsT = ( 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  )  &
+                  &   ) * 0.25_wp * r1_e1e2f(ji,jj)
 
                ! alpha & beta for aEVP
-               IF( ln_aEVP ) THEN
-                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )
-                  z1_alph2 = 1._wp / ( zalph2 + 1._wp )
-                  zbeta(ji,jj) = zalph2
-               ENDIF
-               
-               ! P/delta at F points
-               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) )
+               !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 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) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2
+               zs12(ji,jj) = ( zs12(ji,jj) + stress12tmp(ji,jj) * zalph1  ) * z1_alph1
 
             END DO
@@ -867,4 +885,226 @@
    END SUBROUTINE ice_dyn_rhg_eap
 
+   SUBROUTINE update_stress_rdg (ksub, ndte, phi, divu, tension, &
+                                 shear, a11, a12, &
+                                 stressp,  stressm, &
+                                 stress12, strength, &
+                                 alphar, alphas)
+      !!---------------------------------------------------------------------
+      !!                   ***  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
+      !!---------------------------------------------------------------------
+      INTEGER,  INTENT(in   ) ::   ksub, ndte
+      REAL(wp), INTENT(in   ) ::   phi, strength
+      REAL(wp), INTENT(in   ) ::   divu, tension, shear
+      REAL(wp), INTENT(in   ) ::   a11, a12
+      REAL(wp), INTENT(  out) ::   stressp, stressm, stress12
+      REAL(wp), INTENT(  out) ::   alphar, alphas
+
+      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) ::   invdx, invdy, invda   
+      REAL(wp) ::   dtemp1, dtemp2, atempprime, a, invsin
+
+      REAL(wp), PARAMETER ::   kfriction = 0.45_wp
+      !!---------------------------------------------------------------------
+      ! 1) structure tensor
+      a22 = 1._wp-a11
+      Q11Q11 = 1._wp
+      Q12Q12 = rsmall
+      Q11Q12 = 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) + &
+                             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
+ 
+      ! rotation Q*atemp*Q^T
+      atempprime = Q11Q11*a11 + 2.0_wp*Q11Q12*a12 + Q12Q12*a22 
+      
+      ! make first principal value the largest
+      atempprime = max(atempprime, 1.0_wp - atempprime)
+ 
+      ! 2) strain rate
+      dtemp11 = 0.5_wp*(divu + tension)
+      dtemp12 = shear*0.5_wp
+      dtemp22 = 0.5_wp*(divu - tension)
+
+      ! 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
+      ! In cos and sin values
+      x = 0._wp
+      if ((ABS(dtemp1) > rsmall).or.(ABS(dtemp2) > rsmall)) then
+         x = 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
+      
+      ! y: angle between major principal axis of strain rate and structure
+      ! tensor
+      ! y = gamma - alpha 
+      ! Expressesed componently with
+      ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha)
+      
+      Tany_1 = Q11Q12 - Qd11Qd12
+      Tany_2 = Q11Q11 - Qd12Qd12
+      
+      y = 0._wp
+      
+      if ((ABS(Tany_1) > rsmall).or.(ABS(Tany_2) > rsmall)) then
+         y = atan2(Tany_1,Tany_2)
+      endif
+      
+      ! to make sure y is between 0 and pi
+      if (y > rpi) y = y - rpi
+      if (y < 0)  y = y + 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
+
+      ! % 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 
+      
+      ky  = int(y*invdy) + 1
+      kyw = ky - y*invdy 
+      
+      ka  = int((atempprime-0.5_wp)*invda) + 1
+      kaw = ka - (atempprime-0.5_wp)*invda
+      
+      ! % 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)
+      
+      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)
+                   
+      ! 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
+
+      ! Back - rotation of the stress from principal axes into general coordinates
+
+      ! 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
+
+      ! 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
+   END SUBROUTINE update_stress_rdg 
 
    SUBROUTINE rhg_eap_rst( cdrw, kt )
@@ -873,5 +1113,5 @@
       !!                     
       !! ** Purpose :   Read or write RHG file in restart file
-      !!
+      !!      
       !! ** Method  :   use of IOM library
       !!----------------------------------------------------------------------
@@ -889,5 +1129,4 @@
 
       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)
       !!----------------------------------------------------------------------
       !
