Changeset 13113


Ignore:
Timestamp:
2020-06-16T17:18:24+02:00 (7 weeks ago)
Author:
stefryn
Message:

added structure tensor evolution + some changes for NEMO coding rules, to be continued

File:
1 edited

Legend:

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

    r13105 r13113  
    144144      REAL(wp) ::   stressptmp, stressmtmp, stress12tmpF                ! anisotropic stress tensor components 
    145145      REAL(wp) ::   alphar, alphas                                      ! for mechanical redistribution 
     146      REAL(wp) ::   mresult11, mresult12, z1dtevpkth, p5kth, z1_dtevp_A        ! for structure tensor evolution 
    146147      ! 
    147148      REAL(wp), DIMENSION(jpi,jpj) ::   stress12tmp                     ! anisotropic stress tensor component for regridding 
     
    253254      zs2 (:,:) = pstress2_i (:,:) 
    254255      zs12(:,:) = pstress12_i(:,:) 
     256 
     257      ! constants for structure tensor 
     258      z1_dtevp_A = z1_dtevp/10.0_wp 
     259      z1dtevpkth = 1._wp / (z1_dtevp_A + 0.00002_wp) 
     260      p5kth = 0.5_wp * 0.00002_wp 
    255261 
    256262      ! Ice strength 
     
    412418                                       stress12tmp(ji,jj), strength(ji,jj), & 
    413419                                       alphar, alphas) 
     420 
     421               ! structure tensor update 
     422               IF (mod(jter,10) == 0) THEN 
     423                  CALL calc_ffrac(stressptmp, stressmtmp, stress12tmp(ji,jj), & 
     424                                  paniso_11(ji,jj), paniso_12(ji,jj),           & 
     425                                  mresult11,  mresult12) 
     426 
     427                  paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A  + p5kth - mresult11) * z1dtevpkth ! implicit 
     428                  paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A  - mresult12) * z1dtevpkth ! implicit 
     429               ENDIF 
     430 
     431 
    414432               IF (jter == nn_nevp) THEN 
    415433                  yield11(ji,jj) = 0.5_wp * (stressptmp + stressmtmp) 
     
    428446               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 
    429447               !   alpha = beta = sqrt(4*gamma) 
    430                IF( ln_aEVP ) THEN 
    431                   zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
    432                   z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
    433                ENDIF 
     448               !IF( ln_aEVP ) THEN 
     449               !   zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     450               !   z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     451               !ENDIF 
    434452                
    435453               ! stress at T points (zkt/=0 if landfast) 
     
    444462         END DO 
    445463         !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. ) 
    446          CALL lbc_lnk( 'icedyn_rhg_eap', stress12tmp, 'T', 1. ) 
     464         CALL lbc_lnk_multi( 'icedyn_rhg_eap', stress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.) 
    447465 
    448466         DO jj = 1, jpjm1 
     
    714732!!$         ENDIF 
    715733         ! 
     734 
    716735         !                                                ! ==================== ! 
    717736      END DO                                              !  end loop over jter  ! 
     
    851870      ! --- anisotropy tensor --- ! 
    852871      IF( iom_use('aniso') ) THEN        
    853          CALL lbc_lnk( 'icedyn_rhg_eap', aniso_11, 'T', 1. )  
    854          CALL iom_put( 'aniso' , aniso_11 * zmsk00 ) 
     872         CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1. )  
     873         CALL iom_put( 'aniso' , paniso_11 * zmsk00 ) 
    855874      ENDIF 
    856875       
     
    919938      !!--------------------------------------------------------------------- 
    920939      !!                   ***  SUBROUTINE rhg_eap_rst  *** 
    921       !! Updates the stress depending on values of strain rate and structure 
    922       !! tensor and for ksub=ndte it computes closing and sliding rate 
     940      !!                      
     941      !! ** Purpose :   Updates the stress depending on values of strain rate and structure 
     942      !!                tensor and for ksub=ndte it computes closing and sliding rate 
    923943      !!--------------------------------------------------------------------- 
    924944      INTEGER,  INTENT(in   ) ::   ksub, ndte 
     
    931951      INTEGER  ::   kx ,ky, ka 
    932952 
    933       REAL(wp) ::   stemp11r, stemp12r, stemp22r    
    934       REAL(wp) ::   stemp11s, stemp12s, stemp22s  
    935       REAL(wp) ::   a22, Qd11Qd11, Qd11Qd12, Qd12Qd12  
    936       REAL(wp) ::   Q11Q11, Q11Q12, Q12Q12  
    937       REAL(wp) ::   dtemp11, dtemp12, dtemp22  
    938       REAL(wp) ::   rotstemp11r, rotstemp12r, rotstemp22r    
    939       REAL(wp) ::   rotstemp11s, rotstemp12s, rotstemp22s 
    940       REAL(wp) ::   sig11, sig12, sig22  
    941       REAL(wp) ::   sgprm11, sgprm12, sgprm22  
    942       REAL(wp) ::   invstressconviso  
    943       REAL(wp) ::   Angle_denom_gamma,  Angle_denom_alpha  
    944       REAL(wp) ::   Tany_1, Tany_2  
    945       REAL(wp) ::   x, y, dx, dy, da, kxw, kyw, kaw 
     953      REAL(wp) ::   zstemp11r, zstemp12r, zstemp22r    
     954      REAL(wp) ::   zstemp11s, zstemp12s, zstemp22s  
     955      REAL(wp) ::   za22, zQd11Qd11, zQd11Qd12, zQd12Qd12  
     956      REAL(wp) ::   zQ11Q11, zQ11Q12, zQ12Q12  
     957      REAL(wp) ::   zdtemp11, zdtemp12, zdtemp22  
     958      REAL(wp) ::   zrotstemp11r, zrotstemp12r, zrotstemp22r    
     959      REAL(wp) ::   zrotstemp11s, zrotstemp12s, zrotstemp22s 
     960      REAL(wp) ::   zsig11, zsig12, zsig22  
     961      REAL(wp) ::   zsgprm11, zsgprm12, zsgprm22  
     962      REAL(wp) ::   zinvstressconviso  
     963      REAL(wp) ::   zAngle_denom_gamma,  zAngle_denom_alpha  
     964      REAL(wp) ::   zTany_1, zTany_2  
     965      REAL(wp) ::   zx, zy, zdx, zdy, zda, zkxw, kyw, kaw 
    946966      REAL(wp) ::   invdx, invdy, invda    
    947967      REAL(wp) ::   dtemp1, dtemp2, atempprime, a, invsin 
     
    951971      ! Factor to maintain the same stress as in EVP (see Section 3) 
    952972      ! Can be set to 1 otherwise 
    953       invstressconviso = 1._wp/(1._wp+kfriction*kfriction) 
     973      zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 
    954974  
    955       invsin = 1._wp/sin(2._wp*phi) * invstressconviso  
     975      invsin = 1._wp/sin(2._wp*phi) * zinvstressconviso  
    956976      !now uses phi as set in higher code 
    957977        
     
    960980 
    961981      ! 1) structure tensor 
    962       a22 = 1._wp-a11 
    963       Q11Q11 = 1._wp 
    964       Q12Q12 = rsmall 
    965       Q11Q12 = rsmall 
     982      za22 = 1._wp-a11 
     983      zQ11Q11 = 1._wp 
     984      zQ12Q12 = rsmall 
     985      zQ11Q12 = rsmall 
    966986  
    967987      ! gamma: angle between general coordiantes and principal axis of A  
    968988      ! here Tan2gamma = 2 a12 / (a11 - a22)  
    969       if((ABS(a11 - a22) > rsmall).or.(ABS(a12) > rsmall)) then       
    970          Angle_denom_gamma = 1._wp/sqrt( ( a11 - a22 )*( a11 - a22) + & 
     989      IF((ABS(a11 - za22) > rsmall).OR.(ABS(a12) > rsmall)) THEN       
     990         zAngle_denom_gamma = 1._wp/sqrt( ( a11 - za22 )*( a11 - za22) + & 
    971991                             4._wp*a12*a12 ) 
    972992    
    973          Q11Q11 = 0.5_wp + ( a11 - a22 )*0.5_wp*Angle_denom_gamma   !Cos^2  
    974          Q12Q12 = 0.5_wp - ( a11 - a22 )*0.5_wp*Angle_denom_gamma   !Sin^2 
    975          Q11Q12 = a12*Angle_denom_gamma                     !CosSin  
    976       endif 
     993         zQ11Q11 = 0.5_wp + ( a11 - za22 )*0.5_wp*zAngle_denom_gamma   !Cos^2  
     994         zQ12Q12 = 0.5_wp - ( a11 - za22 )*0.5_wp*zAngle_denom_gamma   !Sin^2 
     995         zQ11Q12 = a12*zAngle_denom_gamma                     !CosSin  
     996      ENDIF 
    977997  
    978998      ! rotation Q*atemp*Q^T 
    979       atempprime = Q11Q11*a11 + 2.0_wp*Q11Q12*a12 + Q12Q12*a22  
     999      atempprime = zQ11Q11*a11 + 2.0_wp*zQ11Q12*a12 + zQ12Q12*za22  
    9801000       
    9811001      ! make first principal value the largest 
     
    9831003  
    9841004      ! 2) strain rate 
    985       dtemp11 = 0.5_wp*(divu + tension) 
    986       dtemp12 = shear*0.5_wp 
    987       dtemp22 = 0.5_wp*(divu - tension) 
     1005      zdtemp11 = 0.5_wp*(divu + tension) 
     1006      zdtemp12 = shear*0.5_wp 
     1007      zdtemp22 = 0.5_wp*(divu - tension) 
    9881008 
    9891009      !WRITE(numout,*) 'div, tension, shear inside loop', ksub, divu, tension, shear 
     
    9911011      ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22)  
    9921012 
    993       Qd11Qd11 = 1.0_wp 
    994       Qd12Qd12 = rsmall 
    995       Qd11Qd12 = rsmall 
    996  
    997       if((ABS( dtemp11 - dtemp22) > rsmall).or. (ABS(dtemp12) > rsmall)) then 
    998  
    999          Angle_denom_alpha = 1.0_wp/sqrt( ( dtemp11 - dtemp22 )* & 
    1000                              ( dtemp11 - dtemp22 ) + 4.0_wp*dtemp12*dtemp12) 
    1001  
    1002          Qd11Qd11 = 0.5_wp + ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Cos^2  
    1003          Qd12Qd12 = 0.5_wp - ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Sin^2 
    1004          Qd11Qd12 = dtemp12*Angle_denom_alpha !CosSin 
    1005       endif 
    1006  
    1007       dtemp1 = Qd11Qd11*dtemp11 + 2.0_wp*Qd11Qd12*dtemp12 + Qd12Qd12*dtemp22 
    1008       dtemp2 = Qd12Qd12*dtemp11 - 2.0_wp*Qd11Qd12*dtemp12 + Qd11Qd11*dtemp22 
     1013      zQd11Qd11 = 1.0_wp 
     1014      zQd12Qd12 = rsmall 
     1015      zQd11Qd12 = rsmall 
     1016 
     1017      IF((ABS( zdtemp11 - zdtemp22) > rsmall).OR. (ABS(zdtemp12) > rsmall)) THEN 
     1018 
     1019         zAngle_denom_alpha = 1.0_wp/sqrt( ( zdtemp11 - zdtemp22 )* & 
     1020                             ( zdtemp11 - zdtemp22 ) + 4.0_wp*zdtemp12*zdtemp12) 
     1021 
     1022         zQd11Qd11 = 0.5_wp + ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Cos^2  
     1023         zQd12Qd12 = 0.5_wp - ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Sin^2 
     1024         zQd11Qd12 = zdtemp12*zAngle_denom_alpha !CosSin 
     1025      ENDIF 
     1026 
     1027      dtemp1 = zQd11Qd11*zdtemp11 + 2.0_wp*zQd11Qd12*zdtemp12 + zQd12Qd12*zdtemp22 
     1028      dtemp2 = zQd12Qd12*zdtemp11 - 2.0_wp*zQd11Qd12*zdtemp12 + zQd11Qd11*zdtemp22 
    10091029      ! In cos and sin values 
    1010       x = 0._wp 
    1011       if ((ABS(dtemp1) > rsmall).or.(ABS(dtemp2) > rsmall)) then 
    1012          x = atan2(dtemp2,dtemp1) 
    1013       endif       
     1030      zx = 0._wp 
     1031      IF ((ABS(dtemp1) > rsmall).OR.(ABS(dtemp2) > rsmall)) THEN 
     1032         zx = atan2(dtemp2,dtemp1) 
     1033      ENDIF       
    10141034                
    10151035      ! to ensure the angle lies between pi/4 and 9 pi/4  
    1016       if (x < rpi*0.25_wp) x = x + rpi*2.0_wp 
     1036      IF (zx < rpi*0.25_wp) zx = zx + rpi*2.0_wp 
    10171037       
    10181038      ! y: angle between major principal axis of strain rate and structure 
     
    10221042      ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha) 
    10231043       
    1024       Tany_1 = Q11Q12 - Qd11Qd12 
    1025       Tany_2 = Q11Q11 - Qd12Qd12 
     1044      zTany_1 = zQ11Q12 - zQd11Qd12 
     1045      zTany_2 = zQ11Q11 - zQd12Qd12 
    10261046       
    1027       y = 0._wp 
     1047      zy = 0._wp 
    10281048       
    1029       if ((ABS(Tany_1) > rsmall).or.(ABS(Tany_2) > rsmall)) then 
    1030          y = atan2(Tany_1,Tany_2) 
    1031       endif 
     1049      IF ((ABS(zTany_1) > rsmall).OR.(ABS(zTany_2) > rsmall)) THEN 
     1050         zy = atan2(zTany_1,zTany_2) 
     1051      ENDIF 
    10321052       
    10331053      ! to make sure y is between 0 and pi 
    1034       if (y > rpi) y = y - rpi 
    1035       if (y < 0)  y = y + rpi 
     1054      IF (zy > rpi) zy = zy - rpi 
     1055      IF (zy < 0)  zy = zy + rpi 
    10361056 
    10371057      ! 3) update anisotropic stress tensor  
    1038       dx   = rpi/real(nx_yield-1,kind=wp) 
    1039       dy   = rpi/real(ny_yield-1,kind=wp) 
    1040       da   = 0.5_wp/real(na_yield-1,kind=wp) 
    1041       invdx = 1._wp/dx 
    1042       invdy = 1._wp/dy 
    1043       invda = 1._wp/da 
     1058      zdx   = rpi/real(nx_yield-1,kind=wp) 
     1059      zdy   = rpi/real(ny_yield-1,kind=wp) 
     1060      zda   = 0.5_wp/real(na_yield-1,kind=wp) 
     1061      invdx = 1._wp/zdx 
     1062      invdy = 1._wp/zdy 
     1063      invda = 1._wp/zda 
    10441064 
    10451065      ! % need 8 coords and 8 weights 
    10461066      ! % range in kx 
    1047       kx  = int((x-rpi*0.25_wp-rpi)*invdx) + 1 
    1048       kxw = kx - (x-rpi*0.25_wp-rpi)*invdx  
     1067      kx  = int((zx-rpi*0.25_wp-rpi)*invdx) + 1 
     1068      zkxw = kx - (zx-rpi*0.25_wp-rpi)*invdx  
    10491069       
    1050       ky  = int(y*invdy) + 1 
    1051       kyw = ky - y*invdy  
     1070      ky  = int(zy*invdy) + 1 
     1071      kyw = ky - zy*invdy  
    10521072       
    10531073      ka  = int((atempprime-0.5_wp)*invda) + 1 
     
    10551075       
    10561076      ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013) 
    1057       stemp11r =  kxw   * kyw         * kaw         * s11r(kx  ,ky  ,ka  ) & 
    1058         & + (1._wp-kxw) * kyw         * kaw         * s11r(kx+1,ky  ,ka  ) & 
    1059         & + kxw         * (1._wp-kyw) * kaw         * s11r(kx  ,ky+1,ka  ) & 
    1060         & + kxw         * kyw         * (1._wp-kaw) * s11r(kx  ,ky  ,ka+1) & 
    1061         & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s11r(kx+1,ky+1,ka  ) & 
    1062         & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s11r(kx+1,ky  ,ka+1) & 
    1063         & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s11r(kx  ,ky+1,ka+1) & 
    1064         & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1) 
    1065       stemp12r =  kxw   * kyw         * kaw         * s12r(kx  ,ky  ,ka  ) & 
    1066         & + (1._wp-kxw) * kyw         * kaw         * s12r(kx+1,ky  ,ka  ) & 
    1067         & + kxw         * (1._wp-kyw) * kaw         * s12r(kx  ,ky+1,ka  ) & 
    1068         & + kxw         * kyw         * (1._wp-kaw) * s12r(kx  ,ky  ,ka+1) & 
    1069         & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s12r(kx+1,ky+1,ka  ) & 
    1070         & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s12r(kx+1,ky  ,ka+1) & 
    1071         & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s12r(kx  ,ky+1,ka+1) & 
    1072         & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1) 
    1073       stemp22r = kxw    * kyw         * kaw         * s22r(kx  ,ky  ,ka  ) & 
    1074         & + (1._wp-kxw) * kyw         * kaw         * s22r(kx+1,ky  ,ka  ) & 
    1075         & + kxw         * (1._wp-kyw) * kaw         * s22r(kx  ,ky+1,ka  ) & 
    1076         & + kxw         * kyw         * (1._wp-kaw) * s22r(kx  ,ky  ,ka+1) & 
    1077         & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s22r(kx+1,ky+1,ka  ) & 
    1078         & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s22r(kx+1,ky  ,ka+1) & 
    1079         & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s22r(kx  ,ky+1,ka+1) & 
    1080         & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 
     1077      zstemp11r =  zkxw   * kyw         * kaw         * s11r(kx  ,ky  ,ka  ) & 
     1078        & + (1._wp-zkxw) * kyw         * kaw         * s11r(kx+1,ky  ,ka  ) & 
     1079        & + zkxw         * (1._wp-kyw) * kaw         * s11r(kx  ,ky+1,ka  ) & 
     1080        & + zkxw         * kyw         * (1._wp-kaw) * s11r(kx  ,ky  ,ka+1) & 
     1081        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s11r(kx+1,ky+1,ka  ) & 
     1082        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s11r(kx+1,ky  ,ka+1) & 
     1083        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s11r(kx  ,ky+1,ka+1) & 
     1084        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1) 
     1085      zstemp12r =  zkxw   * kyw         * kaw         * s12r(kx  ,ky  ,ka  ) & 
     1086        & + (1._wp-zkxw) * kyw         * kaw         * s12r(kx+1,ky  ,ka  ) & 
     1087        & + zkxw         * (1._wp-kyw) * kaw         * s12r(kx  ,ky+1,ka  ) & 
     1088        & + zkxw         * kyw         * (1._wp-kaw) * s12r(kx  ,ky  ,ka+1) & 
     1089        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s12r(kx+1,ky+1,ka  ) & 
     1090        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s12r(kx+1,ky  ,ka+1) & 
     1091        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s12r(kx  ,ky+1,ka+1) & 
     1092        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1) 
     1093      zstemp22r = zkxw    * kyw         * kaw         * s22r(kx  ,ky  ,ka  ) & 
     1094        & + (1._wp-zkxw) * kyw         * kaw         * s22r(kx+1,ky  ,ka  ) & 
     1095        & + zkxw         * (1._wp-kyw) * kaw         * s22r(kx  ,ky+1,ka  ) & 
     1096        & + zkxw         * kyw         * (1._wp-kaw) * s22r(kx  ,ky  ,ka+1) & 
     1097        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s22r(kx+1,ky+1,ka  ) & 
     1098        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s22r(kx+1,ky  ,ka+1) & 
     1099        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s22r(kx  ,ky+1,ka+1) & 
     1100        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 
    10811101       
    1082       stemp11s =  kxw   * kyw         * kaw         * s11s(kx  ,ky  ,ka  ) & 
    1083         & + (1._wp-kxw) * kyw         * kaw         * s11s(kx+1,ky  ,ka  ) & 
    1084         & + kxw         * (1._wp-kyw) * kaw         * s11s(kx  ,ky+1,ka  ) & 
    1085         & + kxw         * kyw         * (1._wp-kaw) * s11s(kx  ,ky  ,ka+1) & 
    1086         & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s11s(kx+1,ky+1,ka  ) & 
    1087         & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s11s(kx+1,ky  ,ka+1) & 
    1088         & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s11s(kx  ,ky+1,ka+1) & 
    1089         & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1) 
    1090       stemp12s =  kxw   * kyw         * kaw         * s12s(kx  ,ky  ,ka  ) & 
    1091         & + (1._wp-kxw) * kyw         * kaw         * s12s(kx+1,ky  ,ka  ) & 
    1092         & + kxw         * (1._wp-kyw) * kaw         * s12s(kx  ,ky+1,ka  ) & 
    1093         & + kxw         * kyw         * (1._wp-kaw) * s12s(kx  ,ky  ,ka+1) & 
    1094         & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s12s(kx+1,ky+1,ka  ) & 
    1095         & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s12s(kx+1,ky  ,ka+1) & 
    1096         & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s12s(kx  ,ky+1,ka+1) & 
    1097         & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1) 
    1098       stemp22s =  kxw   * kyw         * kaw         * s22s(kx  ,ky  ,ka  ) & 
    1099         & + (1._wp-kxw) * kyw         * kaw         * s22s(kx+1,ky  ,ka  ) & 
    1100         & + kxw         * (1._wp-kyw) * kaw         * s22s(kx  ,ky+1,ka  ) & 
    1101         & + kxw         * kyw         * (1._wp-kaw) * s22s(kx  ,ky  ,ka+1) & 
    1102         & + (1._wp-kxw) * (1._wp-kyw) * kaw         * s22s(kx+1,ky+1,ka  ) & 
    1103         & + (1._wp-kxw) * kyw         * (1._wp-kaw) * s22s(kx+1,ky  ,ka+1) & 
    1104         & + kxw         * (1._wp-kyw) * (1._wp-kaw) * s22s(kx  ,ky+1,ka+1) & 
    1105         & + (1._wp-kxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1) 
     1102      zstemp11s =  zkxw   * kyw         * kaw         * s11s(kx  ,ky  ,ka  ) & 
     1103        & + (1._wp-zkxw) * kyw         * kaw         * s11s(kx+1,ky  ,ka  ) & 
     1104        & + zkxw         * (1._wp-kyw) * kaw         * s11s(kx  ,ky+1,ka  ) & 
     1105        & + zkxw         * kyw         * (1._wp-kaw) * s11s(kx  ,ky  ,ka+1) & 
     1106        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s11s(kx+1,ky+1,ka  ) & 
     1107        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s11s(kx+1,ky  ,ka+1) & 
     1108        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s11s(kx  ,ky+1,ka+1) & 
     1109        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1) 
     1110      zstemp12s =  zkxw   * kyw         * kaw         * s12s(kx  ,ky  ,ka  ) & 
     1111        & + (1._wp-zkxw) * kyw         * kaw         * s12s(kx+1,ky  ,ka  ) & 
     1112        & + zkxw         * (1._wp-kyw) * kaw         * s12s(kx  ,ky+1,ka  ) & 
     1113        & + zkxw         * kyw         * (1._wp-kaw) * s12s(kx  ,ky  ,ka+1) & 
     1114        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s12s(kx+1,ky+1,ka  ) & 
     1115        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s12s(kx+1,ky  ,ka+1) & 
     1116        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s12s(kx  ,ky+1,ka+1) & 
     1117        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1) 
     1118      zstemp22s =  zkxw   * kyw         * kaw         * s22s(kx  ,ky  ,ka  ) & 
     1119        & + (1._wp-zkxw) * kyw         * kaw         * s22s(kx+1,ky  ,ka  ) & 
     1120        & + zkxw         * (1._wp-kyw) * kaw         * s22s(kx  ,ky+1,ka  ) & 
     1121        & + zkxw         * kyw         * (1._wp-kaw) * s22s(kx  ,ky  ,ka+1) & 
     1122        & + (1._wp-zkxw) * (1._wp-kyw) * kaw         * s22s(kx+1,ky+1,ka  ) & 
     1123        & + (1._wp-zkxw) * kyw         * (1._wp-kaw) * s22s(kx+1,ky  ,ka+1) & 
     1124        & + zkxw         * (1._wp-kyw) * (1._wp-kaw) * s22s(kx  ,ky+1,ka+1) & 
     1125        & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1) 
    11061126                    
    11071127      ! Calculate mean ice stress over a collection of floes (Equation 3 in 
    11081128      ! Tsamados 2013) 
    11091129 
    1110       sig11  = strength*(stemp11r + kfriction*stemp11s) * invsin 
    1111       sig12  = strength*(stemp12r + kfriction*stemp12s) * invsin 
    1112       sig22  = strength*(stemp22r + kfriction*stemp22s) * invsin 
     1130      zsig11  = strength*(zstemp11r + kfriction*zstemp11s) * invsin 
     1131      zsig12  = strength*(zstemp12r + kfriction*zstemp12s) * invsin 
     1132      zsig22  = strength*(zstemp22r + kfriction*zstemp22s) * invsin 
    11131133 
    11141134      !WRITE(numout,*) 'strength inside loop', strength 
     
    11181138 
    11191139      ! Update stress 
    1120       sgprm11 = Q11Q11*sig11 + Q12Q12*sig22 -        2._wp*Q11Q12 *sig12 
    1121       sgprm12 = Q11Q12*sig11 - Q11Q12*sig22 + (Q11Q11 - Q12Q12)*sig12 
    1122       sgprm22 = Q12Q12*sig11 + Q11Q11*sig22 +        2._wp*Q11Q12 *sig12 
    1123  
    1124       stressp  = sgprm11 + sgprm22 
    1125       stress12 = sgprm12 
    1126       stressm  = sgprm11 - sgprm22 
     1140      zsgprm11 = zQ11Q11*zsig11 + zQ12Q12*zsig22 -        2._wp*zQ11Q12 *zsig12 
     1141      zsgprm12 = zQ11Q12*zsig11 - zQ11Q12*zsig22 + (zQ11Q11 - zQ12Q12)*zsig12 
     1142      zsgprm22 = zQ12Q12*zsig11 + zQ11Q11*zsig22 +        2._wp*zQ11Q12 *zsig12 
     1143 
     1144      stressp  = zsgprm11 + zsgprm22 
     1145      stress12 = zsgprm12 
     1146      stressm  = zsgprm11 - zsgprm22 
    11271147 
    11281148      !WRITE(numout,*) 'stress output inside loop', ksub, stressp 
     
    11301150      ! Compute ridging and sliding functions in general coordinates  
    11311151      ! (Equation 11 in Tsamados 2013) 
    1132       if (ksub == ndte) then 
    1133          rotstemp11r = Q11Q11*stemp11r - 2._wp*Q11Q12* stemp12r & 
    1134                      + Q12Q12*stemp22r 
    1135          rotstemp12r = Q11Q11*stemp12r +    Q11Q12*(stemp11r-stemp22r) & 
    1136                      - Q12Q12*stemp12r 
    1137          rotstemp22r = Q12Q12*stemp11r + 2._wp*Q11Q12* stemp12r &  
    1138                      + Q11Q11*stemp22r 
    1139  
    1140          rotstemp11s = Q11Q11*stemp11s - 2._wp*Q11Q12* stemp12s & 
    1141                      + Q12Q12*stemp22s 
    1142          rotstemp12s = Q11Q11*stemp12s +    Q11Q12*(stemp11s-stemp22s) & 
    1143                      - Q12Q12*stemp12s 
    1144          rotstemp22s = Q12Q12*stemp11s + 2._wp*Q11Q12* stemp12s &  
    1145                      + Q11Q11*stemp22s 
    1146  
    1147          alphar =  rotstemp11r*dtemp11 + 2._wp*rotstemp12r*dtemp12 & 
    1148                  + rotstemp22r*dtemp22 
    1149          alphas =  rotstemp11s*dtemp11 + 2._wp*rotstemp12s*dtemp12 & 
    1150                  + rotstemp22s*dtemp22 
    1151       endif 
     1152      IF (ksub == ndte) THEN 
     1153         zrotstemp11r = zQ11Q11*zstemp11r - 2._wp*zQ11Q12* zstemp12r & 
     1154                     + zQ12Q12*zstemp22r 
     1155         zrotstemp12r = zQ11Q11*zstemp12r +    zQ11Q12*(zstemp11r-zstemp22r) & 
     1156                     - zQ12Q12*zstemp12r 
     1157         zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r &  
     1158                     + zQ11Q11*zstemp22r 
     1159 
     1160         zrotstemp11s = zQ11Q11*zstemp11s - 2._wp*zQ11Q12* zstemp12s & 
     1161                     + zQ12Q12*zstemp22s 
     1162         zrotstemp12s = zQ11Q11*zstemp12s +    zQ11Q12*(zstemp11s-zstemp22s) & 
     1163                     - zQ12Q12*zstemp12s 
     1164         zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s &  
     1165                     + zQ11Q11*zstemp22s 
     1166 
     1167         alphar =  zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 & 
     1168                 + zrotstemp22r*zdtemp22 
     1169         alphas =  zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 & 
     1170                 + zrotstemp22s*zdtemp22 
     1171      ENDIF 
    11521172   END SUBROUTINE update_stress_rdg  
     1173 
     1174!======================================================================= 
     1175 
     1176 
     1177   SUBROUTINE calc_ffrac (stressp, stressm, stress12, & 
     1178                          a11, a12,                   & 
     1179                          mresult11, mresult12) 
     1180      !!--------------------------------------------------------------------- 
     1181      !!                   ***  ROUTINE calc_ffrac  *** 
     1182      !!                      
     1183      !! ** Purpose :   Computes term in evolution equation for structure tensor 
     1184      !!                which determines the ice floe re-orientation due to fracture 
     1185      !! ** Method :    Eq. 7: Ffrac = -kf(A-S) or = 0 depending on sigma_1 and sigma_2 
     1186      !!--------------------------------------------------------------------- 
     1187      REAL (wp), INTENT(in)  ::   stressp, stressm, stress12, a11, a12 
     1188      REAL (wp), INTENT(out) ::   mresult11, mresult12 
     1189 
     1190      ! local variables 
     1191 
     1192      REAL (wp) ::   sigma11, sigma12, sigma22  ! stress tensor elements 
     1193      REAL (wp) ::   Angle_denom        ! angle with principal component axis 
     1194      REAL (wp) ::   sigma_1, sigma_2   ! principal components of stress 
     1195      REAL (wp) ::   Q11, Q12, Q11Q11, Q11Q12, Q12Q12 
     1196 
     1197      REAL (wp), PARAMETER ::   kfrac = 0.0001_wp   ! rate of fracture formation  
     1198      REAL (wp), PARAMETER ::   threshold = 0.3_wp  ! critical confinement ratio  
     1199      !!--------------------------------------------------------------- 
     1200      ! 
     1201      sigma11 = 0.5_wp*(stressp+stressm)  
     1202      sigma12 = stress12  
     1203      sigma22 = 0.5_wp*(stressp-stressm)  
     1204 
     1205      ! Here's the change - no longer calculate gamma, 
     1206      ! use trig formulation, angle signs are kept correct, don't worry 
     1207    
     1208      ! rotate tensor to get into sigma principal axis 
     1209    
     1210      ! here Tan2gamma = 2 sig12 / (sig11 - sig12)  
     1211      ! add rsmall to the denominator to stop 1/0 errors, makes very little  
     1212      ! error to the calculated angles < rsmall 
     1213 
     1214      Q11Q11 = 0.1_wp 
     1215      Q12Q12 = rsmall 
     1216      Q11Q12 = rsmall 
     1217 
     1218      IF((ABS( sigma11 - sigma22) > rsmall).OR.(ABS(sigma12) > rsmall)) THEN 
     1219 
     1220         Angle_denom = 1.0_wp/sqrt( ( sigma11 - sigma22 )*( sigma11 - & 
     1221                       sigma22 ) + 4.0_wp*sigma12*sigma12) 
     1222 
     1223         Q11Q11 = 0.5_wp + ( sigma11 - sigma22 )*0.5_wp*Angle_denom   !Cos^2  
     1224         Q12Q12 = 0.5_wp - ( sigma11 - sigma22 )*0.5_wp*Angle_denom   !Sin^2 
     1225         Q11Q12 = sigma12*Angle_denom                      !CosSin  
     1226      ENDIF 
     1227 
     1228      sigma_1 = Q11Q11*sigma11 + 2.0_wp*Q11Q12*sigma12  & 
     1229              + Q12Q12*sigma22 ! S(1,1) 
     1230      sigma_2 = Q12Q12*sigma11 - 2.0_wp*Q11Q12*sigma12  & 
     1231              + Q11Q11*sigma22 ! S(2,2) 
     1232 
     1233      ! Pure divergence 
     1234      IF ((sigma_1 >= 0.0_wp).AND.(sigma_2 >= 0.0_wp))  THEN 
     1235         mresult11 = 0.0_wp 
     1236         mresult12 = 0.0_wp 
     1237 
     1238      ! Unconfined compression: cracking of blocks not along the axial splitting 
     1239      ! direction 
     1240      ! which leads to the loss of their shape, so we again model it through diffusion 
     1241      ELSEIF ((sigma_1 >= 0.0_wp).AND.(sigma_2 < 0.0_wp))  THEN 
     1242         mresult11 = kfrac * (a11 - Q12Q12) 
     1243         mresult12 = kfrac * (a12 + Q11Q12) 
     1244 
     1245      ! Shear faulting 
     1246      ELSEIF (sigma_2 == 0.0_wp) THEN 
     1247         mresult11 = 0.0_wp 
     1248         mresult12 = 0.0_wp 
     1249      ELSEIF ((sigma_1 <= 0.0_wp).AND.(sigma_1/sigma_2 <= threshold)) THEN 
     1250         mresult11 = kfrac * (a11 - Q12Q12) 
     1251         mresult12 = kfrac * (a12 + Q11Q12) 
     1252 
     1253      ! Horizontal spalling 
     1254      ELSE  
     1255         mresult11 = 0.0_wp 
     1256         mresult12 = 0.0_wp 
     1257      ENDIF 
     1258 
     1259   END SUBROUTINE calc_ffrac 
     1260 
    11531261 
    11541262   SUBROUTINE rhg_eap_rst( cdrw, kt ) 
     
    12291337         s12s(ix,iy,ia) = 0._wp 
    12301338         s22s(ix,iy,ia) = 0._wp 
    1231          IF (ia <= na_yield-1) then 
     1339         IF (ia <= na_yield-1) THEN 
    12321340          DO iz=1,nz 
    12331341          s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 
     
    12991407    
    13001408      w1 = - 223.87569446_wp & 
    1301            + 2361.2198663_wp*a & 
    1302            - 10606.56079975_wp*a*a & 
    1303            + 26315.50025642_wp*a*a*a & 
    1304            - 38948.30444297_wp*a*a*a*a & 
    1305            + 34397.72407466_wp*a*a*a*a*a & 
    1306            - 16789.98003081_wp*a*a*a*a*a*a & 
    1307            + 3495.82839237_wp*a*a*a*a*a*a*a 
     1409       &   + 2361.2198663_wp*a & 
     1410       &   - 10606.56079975_wp*a*a & 
     1411       &   + 26315.50025642_wp*a*a*a & 
     1412       &   - 38948.30444297_wp*a*a*a*a & 
     1413       &   + 34397.72407466_wp*a*a*a*a*a & 
     1414       &   - 16789.98003081_wp*a*a*a*a*a*a & 
     1415       &   + 3495.82839237_wp*a*a*a*a*a*a*a 
    13081416    
    13091417   END FUNCTION w1 
     
    13191427    
    13201428      w2 = - 6670.68911883_wp & 
    1321            + 70222.33061536_wp*a & 
    1322            - 314871.71525448_wp*a*a & 
    1323            + 779570.02793492_wp*a*a*a & 
    1324            - 1151098.82436864_wp*a*a*a*a & 
    1325            + 1013896.59464498_wp*a*a*a*a*a & 
    1326            - 493379.44906738_wp*a*a*a*a*a*a & 
    1327            + 102356.551518_wp*a*a*a*a*a*a*a 
     1429       &   + 70222.33061536_wp*a & 
     1430       &   - 314871.71525448_wp*a*a & 
     1431       &   + 779570.02793492_wp*a*a*a & 
     1432       &   - 1151098.82436864_wp*a*a*a*a & 
     1433       &   + 1013896.59464498_wp*a*a*a*a*a & 
     1434       &   - 493379.44906738_wp*a*a*a*a*a*a & 
     1435       &   + 102356.551518_wp*a*a*a*a*a*a*a 
    13281436    
    13291437   END FUNCTION w2 
     
    13751483      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    13761484    
    1377       IF (-IIn1t2>=rsmall) then 
     1485      IF (-IIn1t2>=rsmall) THEN 
    13781486      Hen1t2 = 1._wp 
    13791487      ELSE 
     
    13811489      ENDIF 
    13821490    
    1383       IF (-IIn2t1>=rsmall) then 
     1491      IF (-IIn2t1>=rsmall) THEN 
    13841492      Hen2t1 = 1._wp 
    13851493      ELSE 
     
    14331541      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    14341542    
    1435       IF (-IIn1t2>=rsmall) then 
     1543      IF (-IIn1t2>=rsmall) THEN 
    14361544      Hen1t2 = 1._wp 
    14371545      ELSE 
     
    14391547      ENDIF 
    14401548    
    1441       IF (-IIn2t1>=rsmall) then 
     1549      IF (-IIn2t1>=rsmall) THEN 
    14421550      Hen2t1 = 1._wp 
    14431551      ELSE 
     
    14931601      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    14941602 
    1495       IF (-IIn1t2>=rsmall) then 
     1603      IF (-IIn1t2>=rsmall) THEN 
    14961604      Hen1t2 = 1._wp 
    14971605      ELSE 
     
    14991607      ENDIF 
    15001608 
    1501       IF (-IIn2t1>=rsmall) then 
     1609      IF (-IIn2t1>=rsmall) THEN 
    15021610      Hen2t1 = 1._wp 
    15031611      ELSE 
     
    15511659      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    15521660 
    1553       IF (-IIn1t2>=rsmall) then 
     1661      IF (-IIn1t2>=rsmall) THEN 
    15541662      Hen1t2 = 1._wp 
    15551663      ELSE 
     
    15571665      ENDIF 
    15581666 
    1559       IF (-IIn2t1>=rsmall) then 
     1667      IF (-IIn2t1>=rsmall) THEN 
    15601668      Hen2t1 = 1._wp 
    15611669      ELSE 
     
    16091717      IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 
    16101718 
    1611       IF (-IIn1t2>=rsmall) then 
     1719      IF (-IIn1t2>=rsmall) THEN 
    16121720      Hen1t2 = 1._wp 
    16131721      ELSE 
     
    16151723      ENDIF 
    16161724 
    1617       IF (-IIn2t1>=rsmall) then 
     1725      IF (-IIn2t1>=rsmall) THEN 
    16181726      Hen2t1 = 1._wp 
    16191727      ELSE 
Note: See TracChangeset for help on using the changeset viewer.