Changeset 11933


Ignore:
Timestamp:
2019-11-19T20:10:30+01:00 (6 months ago)
Author:
stefryn
Message:

add eap stresses, constant anisotropy for now

File:
1 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90

    r11920 r11933  
    1313   !!            3.7  !  2017     (C. Rousset)  add aEVP (Kimmritz 2016-2017) 
    1414   !!            4.0  !  2018     (many people) SI3 [aka Sea Ice cube] 
     15   !!                 !  2019     (S. Rynders)  change into eap rheology from 
     16   !!                                           CICE code (Tsamados, Heorton)  
    1517   !!---------------------------------------------------------------------- 
    1618#if defined key_si3 
     
    4648   PUBLIC   ice_dyn_rhg_eap   ! called by icedyn_rhg.F90 
    4749   PUBLIC   rhg_eap_rst       ! called by icedyn_rhg.F90 
     50 
     51   REAL(wp), PARAMETER           ::   phi = 3.141592653589793_wp/12._wp    ! diamond shaped floe smaller angle (default phi = 30 deg) 
    4852 
    4953   ! look-up table for calculating structure tensor 
     
    126130      REAL(wp) ::   zdtevp, z1_dtevp                                    ! time step for subcycling 
    127131      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity 
    128       REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017 
     132      REAL(wp) ::   zalph1, z1_alph1                                    ! alpha coef from Bouillon 2009 or Kimmritz 2017 
    129133      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume 
    130       REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2      ! temporary scalars 
     134      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2, zdsT ! temporary scalars 
    131135      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars 
    132136      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
     
    137141      REAL(wp) ::   zfac_x, zfac_y 
    138142      REAL(wp) ::   zshear, zdum1, zdum2 
     143      REAL(wp) ::   stressptmp, stressmtmp                              ! anisotropic stress tensor components 
     144      REAL(wp) ::   alphar, alphas                                      ! for mechanical redistribution 
    139145      ! 
     146      REAL(wp), DIMENSION(jpi,jpj) ::   stress12tmp                     ! anisotropic stress tensor component for regridding 
    140147      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
    141148      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
     
    237244      IF( .NOT. ln_aEVP ) THEN 
    238245         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 
    239          zalph2 = zalph1 * z1_ecc2 
    240  
    241246         z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
    242          z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
    243247      ENDIF 
    244248          
     
    376380            DO ji = 2, jpi ! no vector loop 
    377381 
     382               ! shear at T points  
     383               zdsT = ( zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
     384                  &   + zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * e1e2f(ji-1,jj-1)  & 
     385                  &   ) * 0.25_wp * r1_e1e2t(ji,jj) 
     386 
    378387               ! shear**2 at T points (doc eq. A16) 
    379388               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e1e2f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e1e2f(ji-1,jj  )  & 
     
    393402               zdt2 = zdt * zdt 
    394403                
     404               ! --- anisotropic stress calculation --- ! 
     405               CALL update_stress_rdg (jter, nn_nevp, phi, zdiv, zdt, & 
     406                                       zdsT, aniso_11(ji,jj), aniso_12(ji,jj), & 
     407                                       stressptmp, stressmtmp, & 
     408                                       stress12tmp(ji,jj), strength(ji,jj), & 
     409                                       alphar, alphas) 
     410 
    395411               ! delta at T points 
    396                zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
     412               !zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
    397413 
    398414               ! P/delta at T points 
    399                zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
     415               !zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
    400416 
    401417               ! alpha & beta for aEVP 
    402418               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 
    403419               !   alpha = beta = sqrt(4*gamma) 
    404                IF( ln_aEVP ) THEN 
    405                   zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
    406                   z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
    407                   zalph2   = zalph1 
    408                   z1_alph2 = z1_alph1 
    409                ENDIF 
     420               !IF( ln_aEVP ) THEN 
     421               !   zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     422               !   z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     423               !ENDIF 
    410424                
    411425               ! stress at T points (zkt/=0 if landfast) 
    412                zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 
    413                zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
     426               !zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 
     427               !zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
     428               zs1(ji,jj) = ( zs1(ji,jj) + stressptmp * zalph1  ) * z1_alph1 
     429               zs2(ji,jj) = ( zs2(ji,jj) + stressmtmp * zalph1  ) * z1_alph1 
    414430              
    415431            END DO 
    416432         END DO 
    417          CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. ) 
     433         !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. ) 
     434         CALL lbc_lnk( 'icedyn_rhg_eap', stress12tmp, 'T', 1. ) 
    418435 
    419436         DO jj = 1, jpjm1 
    420437            DO ji = 1, jpim1 
     438               ! stress12tmp at F points  
     439               zdsT = ( stress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + stress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1)  & 
     440                  &   + stress12tmp(ji,jj  ) * e1e2t(ji,jj  ) + stress12tmp(ji+1,jj  ) * e1e2t(ji+1,jj  )  & 
     441                  &   ) * 0.25_wp * r1_e1e2f(ji,jj) 
    421442 
    422443               ! alpha & beta for aEVP 
    423                IF( ln_aEVP ) THEN 
    424                   zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
    425                   z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
    426                   zbeta(ji,jj) = zalph2 
    427                ENDIF 
    428                 
    429                ! P/delta at F points 
    430                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) ) 
     444               !IF( ln_aEVP ) THEN 
     445               !   zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     446               !   z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     447               !ENDIF 
    431448                
    432449               ! stress at F points (zkt/=0 if landfast) 
    433                zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 
     450               !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 
     451               zs12(ji,jj) = ( zs12(ji,jj) + stress12tmp(ji,jj) * zalph1  ) * z1_alph1 
    434452 
    435453            END DO 
     
    867885   END SUBROUTINE ice_dyn_rhg_eap 
    868886 
     887   SUBROUTINE update_stress_rdg (ksub, ndte, phi, divu, tension, & 
     888                                 shear, a11, a12, & 
     889                                 stressp,  stressm, & 
     890                                 stress12, strength, & 
     891                                 alphar, alphas) 
     892      !!--------------------------------------------------------------------- 
     893      !!                   ***  SUBROUTINE rhg_eap_rst  *** 
     894      !! Updates the stress depending on values of strain rate and structure 
     895      !! tensor and for ksub=ndte it computes closing and sliding rate 
     896      !!--------------------------------------------------------------------- 
     897      INTEGER,  INTENT(in   ) ::   ksub, ndte 
     898      REAL(wp), INTENT(in   ) ::   phi, strength 
     899      REAL(wp), INTENT(in   ) ::   divu, tension, shear 
     900      REAL(wp), INTENT(in   ) ::   a11, a12 
     901      REAL(wp), INTENT(  out) ::   stressp, stressm, stress12 
     902      REAL(wp), INTENT(  out) ::   alphar, alphas 
     903 
     904      INTEGER  ::   kx ,ky, ka 
     905 
     906      REAL(wp) ::   stemp11r, stemp12r, stemp22r    
     907      REAL(wp) ::   stemp11s, stemp12s, stemp22s  
     908      REAL(wp) ::   a22, Qd11Qd11, Qd11Qd12, Qd12Qd12  
     909      REAL(wp) ::   Q11Q11, Q11Q12, Q12Q12  
     910      REAL(wp) ::   dtemp11, dtemp12, dtemp22  
     911      REAL(wp) ::   rotstemp11r, rotstemp12r, rotstemp22r    
     912      REAL(wp) ::   rotstemp11s, rotstemp12s, rotstemp22s 
     913      REAL(wp) ::   sig11, sig12, sig22  
     914      REAL(wp) ::   sgprm11, sgprm12, sgprm22  
     915      REAL(wp) ::   invstressconviso  
     916      REAL(wp) ::   Angle_denom_gamma,  Angle_denom_alpha  
     917      REAL(wp) ::   Tany_1, Tany_2  
     918      REAL(wp) ::   x, y, dx, dy, da, kxw, kyw, kaw 
     919      REAL(wp) ::   invdx, invdy, invda    
     920      REAL(wp) ::   dtemp1, dtemp2, atempprime, a, invsin 
     921 
     922      REAL(wp), PARAMETER ::   kfriction = 0.45_wp 
     923      !!--------------------------------------------------------------------- 
     924      ! 1) structure tensor 
     925      a22 = 1._wp-a11 
     926      Q11Q11 = 1._wp 
     927      Q12Q12 = rsmall 
     928      Q11Q12 = rsmall 
     929  
     930      ! gamma: angle between general coordiantes and principal axis of A  
     931      ! here Tan2gamma = 2 a12 / (a11 - a22)  
     932      if((ABS(a11 - a22) > rsmall).or.(ABS(a12) > rsmall)) then       
     933         Angle_denom_gamma = 1._wp/sqrt( ( a11 - a22 )*( a11 - a22) + & 
     934                             4._wp*a12*a12 ) 
     935    
     936         Q11Q11 = 0.5_wp + ( a11 - a22 )*0.5_wp*Angle_denom_gamma   !Cos^2  
     937         Q12Q12 = 0.5_wp - ( a11 - a22 )*0.5_wp*Angle_denom_gamma   !Sin^2 
     938         Q11Q12 = a12*Angle_denom_gamma                     !CosSin  
     939      endif 
     940  
     941      ! rotation Q*atemp*Q^T 
     942      atempprime = Q11Q11*a11 + 2.0_wp*Q11Q12*a12 + Q12Q12*a22  
     943       
     944      ! make first principal value the largest 
     945      atempprime = max(atempprime, 1.0_wp - atempprime) 
     946  
     947      ! 2) strain rate 
     948      dtemp11 = 0.5_wp*(divu + tension) 
     949      dtemp12 = shear*0.5_wp 
     950      dtemp22 = 0.5_wp*(divu - tension) 
     951 
     952      ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22)  
     953 
     954      Qd11Qd11 = 1.0_wp 
     955      Qd12Qd12 = rsmall 
     956      Qd11Qd12 = rsmall 
     957 
     958      if((ABS( dtemp11 - dtemp22) > rsmall).or. (ABS(dtemp12) > rsmall)) then 
     959 
     960         Angle_denom_alpha = 1.0_wp/sqrt( ( dtemp11 - dtemp22 )* & 
     961                             ( dtemp11 - dtemp22 ) + 4.0_wp*dtemp12*dtemp12) 
     962 
     963         Qd11Qd11 = 0.5_wp + ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Cos^2  
     964         Qd12Qd12 = 0.5_wp - ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Sin^2 
     965         Qd11Qd12 = dtemp12*Angle_denom_alpha !CosSin 
     966      endif 
     967 
     968      dtemp1 = Qd11Qd11*dtemp11 + 2.0_wp*Qd11Qd12*dtemp12 + Qd12Qd12*dtemp22 
     969      dtemp2 = Qd12Qd12*dtemp11 - 2.0_wp*Qd11Qd12*dtemp12 + Qd11Qd11*dtemp22 
     970      ! In cos and sin values 
     971      x = 0._wp 
     972      if ((ABS(dtemp1) > rsmall).or.(ABS(dtemp2) > rsmall)) then 
     973         x = atan2(dtemp2,dtemp1) 
     974      endif       
     975                
     976      ! to ensure the angle lies between pi/4 and 9 pi/4  
     977      if (x < rpi*0.25_wp) x = x + rpi*2.0_wp 
     978       
     979      ! y: angle between major principal axis of strain rate and structure 
     980      ! tensor 
     981      ! y = gamma - alpha  
     982      ! Expressesed componently with 
     983      ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha) 
     984       
     985      Tany_1 = Q11Q12 - Qd11Qd12 
     986      Tany_2 = Q11Q11 - Qd12Qd12 
     987       
     988      y = 0._wp 
     989       
     990      if ((ABS(Tany_1) > rsmall).or.(ABS(Tany_2) > rsmall)) then 
     991         y = atan2(Tany_1,Tany_2) 
     992      endif 
     993       
     994      ! to make sure y is between 0 and pi 
     995      if (y > rpi) y = y - rpi 
     996      if (y < 0)  y = y + rpi 
     997 
     998      ! 3) update anisotropic stress tensor  
     999      dx   = rpi/real(nx_yield-1,kind=wp) 
     1000      dy   = rpi/real(ny_yield-1,kind=wp) 
     1001      da   = 0.5_wp/real(na_yield-1,kind=wp) 
     1002      invdx = 1._wp/dx 
     1003      invdy = 1._wp/dy 
     1004      invda = 1._wp/da 
     1005 
     1006      ! % need 8 coords and 8 weights 
     1007      ! % range in kx 
     1008      kx  = int((x-rpi*0.25_wp-rpi)*invdx) + 1 
     1009      kxw = kx - (x-rpi*0.25_wp-rpi)*invdx  
     1010       
     1011      ky  = int(y*invdy) + 1 
     1012      kyw = ky - y*invdy  
     1013       
     1014      ka  = int((atempprime-0.5_wp)*invda) + 1 
     1015      kaw = ka - (atempprime-0.5_wp)*invda 
     1016       
     1017      ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013) 
     1018      stemp11r =  kxw   * kyw         * kaw         * s11r(kx  ,ky  ,ka  ) & 
     1019        & + (1._wp-kxw) * kyw         * kaw         * s11r(kx+1,ky  ,ka  ) & 
     1020        & + kxw         * (1._wp-kyw) * kaw         * s11r(kx  ,ky+1,ka  ) & 
     1021        & + kxw         * kyw         * (1._wp-kaw) * s11r(kx  ,ky  ,ka+1) & 
     1022        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s11r(kx+1,ky+1,ka  ) & 
     1023        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s11r(kx+1,ky  ,ka+1) & 
     1024        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s11r(kx  ,ky+1,ka+1) & 
     1025        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1) 
     1026      stemp12r =  kxw   * kyw         * kaw         * s12r(kx  ,ky  ,ka  ) & 
     1027        & + (1._wp-kxw) * kyw         * kaw         * s12r(kx+1,ky  ,ka  ) & 
     1028        & + kxw         * (1._wp-kyw) * kaw         * s12r(kx  ,ky+1,ka  ) & 
     1029        & + kxw         * kyw         * (1._wp-kaw) * s12r(kx  ,ky  ,ka+1) & 
     1030        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s12r(kx+1,ky+1,ka  ) & 
     1031        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s12r(kx+1,ky  ,ka+1) & 
     1032        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s12r(kx  ,ky+1,ka+1) & 
     1033        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1) 
     1034      stemp22r = kxw    * kyw         * kaw         * s22r(kx  ,ky  ,ka  ) & 
     1035        & + (1._wp-kxw) * kyw         * kaw         * s22r(kx+1,ky  ,ka  ) & 
     1036        & + kxw         * (1._wp-kyw) * kaw         * s22r(kx  ,ky+1,ka  ) & 
     1037        & + kxw         * kyw         * (1._wp-kaw) * s22r(kx  ,ky  ,ka+1) & 
     1038        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s22r(kx+1,ky+1,ka  ) & 
     1039        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s22r(kx+1,ky  ,ka+1) & 
     1040        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s22r(kx  ,ky+1,ka+1) & 
     1041        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 
     1042       
     1043      stemp11s =  kxw   * kyw         * kaw         * s11s(kx  ,ky  ,ka  ) & 
     1044        & + (1._wp-kxw) * kyw         * kaw         * s11s(kx+1,ky  ,ka  ) & 
     1045        & + kxw         * (1._wp-kyw) * kaw         * s11s(kx  ,ky+1,ka  ) & 
     1046        & + kxw         * kyw         * (1._wp-kaw) * s11s(kx  ,ky  ,ka+1) & 
     1047        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s11s(kx+1,ky+1,ka  ) & 
     1048        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s11s(kx+1,ky  ,ka+1) & 
     1049        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s11s(kx  ,ky+1,ka+1) & 
     1050        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1) 
     1051      stemp12s =  kxw   * kyw         * kaw         * s12s(kx  ,ky  ,ka  ) & 
     1052        & + (1._wp-kxw) * kyw         * kaw         * s12s(kx+1,ky  ,ka  ) & 
     1053        & + kxw         * (1._wp-kyw) * kaw         * s12s(kx  ,ky+1,ka  ) & 
     1054        & + kxw         * kyw         * (1._wp-kaw) * s12s(kx  ,ky  ,ka+1) & 
     1055        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s12s(kx+1,ky+1,ka  ) & 
     1056        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s12s(kx+1,ky  ,ka+1) & 
     1057        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s12s(kx  ,ky+1,ka+1) & 
     1058        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1) 
     1059      stemp22s =  kxw   * kyw         * kaw         * s22s(kx  ,ky  ,ka  ) & 
     1060        & + (1._wp-kxw) * kyw         * kaw         * s22s(kx+1,ky  ,ka  ) & 
     1061        & + kxw         * (1._wp-kyw) * kaw         * s22s(kx  ,ky+1,ka  ) & 
     1062        & + kxw         * kyw         * (1._wp-kaw) * s22s(kx  ,ky  ,ka+1) & 
     1063        & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s22s(kx+1,ky+1,ka  ) & 
     1064        & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s22s(kx+1,ky  ,ka+1) & 
     1065        & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s22s(kx  ,ky+1,ka+1) & 
     1066        & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1) 
     1067                    
     1068      ! Calculate mean ice stress over a collection of floes (Equation 3 in 
     1069      ! Tsamados 2013) 
     1070 
     1071      sig11  = strength*(stemp11r + kfriction*stemp11s) * invsin 
     1072      sig12  = strength*(stemp12r + kfriction*stemp12s) * invsin 
     1073      sig22  = strength*(stemp22r + kfriction*stemp22s) * invsin 
     1074 
     1075      ! Back - rotation of the stress from principal axes into general coordinates 
     1076 
     1077      ! Update stress 
     1078      sgprm11 = Q11Q11*sig11 + Q12Q12*sig22 -        2._wp*Q11Q12 *sig12 
     1079      sgprm12 = Q11Q12*sig11 - Q11Q12*sig22 + (Q11Q11 - Q12Q12)*sig12 
     1080      sgprm22 = Q12Q12*sig11 + Q11Q11*sig22 +        2._wp*Q11Q12 *sig12 
     1081 
     1082      stressp  = sgprm11 + sgprm22 
     1083      stress12 = sgprm12 
     1084      stressm  = sgprm11 - sgprm22 
     1085 
     1086      ! Compute ridging and sliding functions in general coordinates  
     1087      ! (Equation 11 in Tsamados 2013) 
     1088      if (ksub == ndte) then 
     1089         rotstemp11r = Q11Q11*stemp11r - 2._wp*Q11Q12* stemp12r & 
     1090                     + Q12Q12*stemp22r 
     1091         rotstemp12r = Q11Q11*stemp12r +    Q11Q12*(stemp11r-stemp22r) & 
     1092                     - Q12Q12*stemp12r 
     1093         rotstemp22r = Q12Q12*stemp11r + 2._wp*Q11Q12* stemp12r &  
     1094                     + Q11Q11*stemp22r 
     1095 
     1096         rotstemp11s = Q11Q11*stemp11s - 2._wp*Q11Q12* stemp12s & 
     1097                     + Q12Q12*stemp22s 
     1098         rotstemp12s = Q11Q11*stemp12s +    Q11Q12*(stemp11s-stemp22s) & 
     1099                     - Q12Q12*stemp12s 
     1100         rotstemp22s = Q12Q12*stemp11s + 2._wp*Q11Q12* stemp12s &  
     1101                     + Q11Q11*stemp22s 
     1102 
     1103         alphar =  rotstemp11r*dtemp11 + 2._wp*rotstemp12r*dtemp12 & 
     1104                 + rotstemp22r*dtemp22 
     1105         alphas =  rotstemp11s*dtemp11 + 2._wp*rotstemp12s*dtemp12 & 
     1106                 + rotstemp22s*dtemp22 
     1107      endif 
     1108   END SUBROUTINE update_stress_rdg  
    8691109 
    8701110   SUBROUTINE rhg_eap_rst( cdrw, kt ) 
     
    8731113      !!                      
    8741114      !! ** Purpose :   Read or write RHG file in restart file 
    875       !! 
     1115      !!       
    8761116      !! ** Method  :   use of IOM library 
    8771117      !!---------------------------------------------------------------------- 
     
    8891129 
    8901130      REAL(wp), PARAMETER           ::   eps6 = 1.0e-6_wp 
    891       REAL(wp), PARAMETER           ::   phi = 3.141592653589793_wp/12._wp    ! diamond shaped floe smaller angle (default phi = 30 deg) 
    8921131      !!---------------------------------------------------------------------- 
    8931132      ! 
Note: See TracChangeset for help on using the changeset viewer.