Changeset 13113
- Timestamp:
- 2020-06-16T17:18:24+02:00 (5 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90
r13105 r13113 144 144 REAL(wp) :: stressptmp, stressmtmp, stress12tmpF ! anisotropic stress tensor components 145 145 REAL(wp) :: alphar, alphas ! for mechanical redistribution 146 REAL(wp) :: mresult11, mresult12, z1dtevpkth, p5kth, z1_dtevp_A ! for structure tensor evolution 146 147 ! 147 148 REAL(wp), DIMENSION(jpi,jpj) :: stress12tmp ! anisotropic stress tensor component for regridding … … 253 254 zs2 (:,:) = pstress2_i (:,:) 254 255 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 255 261 256 262 ! Ice strength … … 412 418 stress12tmp(ji,jj), strength(ji,jj), & 413 419 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 414 432 IF (jter == nn_nevp) THEN 415 433 yield11(ji,jj) = 0.5_wp * (stressptmp + stressmtmp) … … 428 446 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 429 447 ! alpha = beta = sqrt(4*gamma) 430 IF( ln_aEVP ) THEN431 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 ENDIF448 !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 434 452 435 453 ! stress at T points (zkt/=0 if landfast) … … 444 462 END DO 445 463 !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.) 447 465 448 466 DO jj = 1, jpjm1 … … 714 732 !!$ ENDIF 715 733 ! 734 716 735 ! ! ==================== ! 717 736 END DO ! end loop over jter ! … … 851 870 ! --- anisotropy tensor --- ! 852 871 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 ) 855 874 ENDIF 856 875 … … 919 938 !!--------------------------------------------------------------------- 920 939 !! *** 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 923 943 !!--------------------------------------------------------------------- 924 944 INTEGER, INTENT(in ) :: ksub, ndte … … 931 951 INTEGER :: kx ,ky, ka 932 952 933 REAL(wp) :: stemp11r, stemp12r,stemp22r934 REAL(wp) :: stemp11s, stemp12s,stemp22s935 REAL(wp) :: a22, Qd11Qd11, Qd11Qd12,Qd12Qd12936 REAL(wp) :: Q11Q11, Q11Q12,Q12Q12937 REAL(wp) :: dtemp11, dtemp12,dtemp22938 REAL(wp) :: rotstemp11r, rotstemp12r,rotstemp22r939 REAL(wp) :: rotstemp11s, rotstemp12s,rotstemp22s940 REAL(wp) :: sig11, sig12,sig22941 REAL(wp) :: sgprm11, sgprm12,sgprm22942 REAL(wp) :: invstressconviso943 REAL(wp) :: Angle_denom_gamma,Angle_denom_alpha944 REAL(wp) :: Tany_1,Tany_2945 REAL(wp) :: x, y, dx, dy, da,kxw, kyw, kaw953 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 946 966 REAL(wp) :: invdx, invdy, invda 947 967 REAL(wp) :: dtemp1, dtemp2, atempprime, a, invsin … … 951 971 ! Factor to maintain the same stress as in EVP (see Section 3) 952 972 ! Can be set to 1 otherwise 953 invstressconviso = 1._wp/(1._wp+kfriction*kfriction)973 zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 954 974 955 invsin = 1._wp/sin(2._wp*phi) * invstressconviso975 invsin = 1._wp/sin(2._wp*phi) * zinvstressconviso 956 976 !now uses phi as set in higher code 957 977 … … 960 980 961 981 ! 1) structure tensor 962 a22 = 1._wp-a11963 Q11Q11 = 1._wp964 Q12Q12 = rsmall965 Q11Q12 = rsmall982 za22 = 1._wp-a11 983 zQ11Q11 = 1._wp 984 zQ12Q12 = rsmall 985 zQ11Q12 = rsmall 966 986 967 987 ! gamma: angle between general coordiantes and principal axis of A 968 988 ! here Tan2gamma = 2 a12 / (a11 - a22) 969 if((ABS(a11 - a22) > rsmall).or.(ABS(a12) > rsmall)) then970 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) + & 971 991 4._wp*a12*a12 ) 972 992 973 Q11Q11 = 0.5_wp + ( a11 - a22 )*0.5_wp*Angle_denom_gamma !Cos^2974 Q12Q12 = 0.5_wp - ( a11 - a22 )*0.5_wp*Angle_denom_gamma !Sin^2975 Q11Q12 = a12*Angle_denom_gamma !CosSin976 endif993 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 977 997 978 998 ! rotation Q*atemp*Q^T 979 atempprime = Q11Q11*a11 + 2.0_wp*Q11Q12*a12 + Q12Q12*a22999 atempprime = zQ11Q11*a11 + 2.0_wp*zQ11Q12*a12 + zQ12Q12*za22 980 1000 981 1001 ! make first principal value the largest … … 983 1003 984 1004 ! 2) strain rate 985 dtemp11 = 0.5_wp*(divu + tension)986 dtemp12 = shear*0.5_wp987 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) 988 1008 989 1009 !WRITE(numout,*) 'div, tension, shear inside loop', ksub, divu, tension, shear … … 991 1011 ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22) 992 1012 993 Qd11Qd11 = 1.0_wp994 Qd12Qd12 = rsmall995 Qd11Qd12 = rsmall996 997 if((ABS( dtemp11 - dtemp22) > rsmall).or. (ABS(dtemp12) > rsmall)) then998 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^21003 Qd12Qd12 = 0.5_wp - ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Sin^21004 Qd11Qd12 = dtemp12*Angle_denom_alpha !CosSin1005 endif1006 1007 dtemp1 = Qd11Qd11*dtemp11 + 2.0_wp*Qd11Qd12*dtemp12 + Qd12Qd12*dtemp221008 dtemp2 = Qd12Qd12*dtemp11 - 2.0_wp*Qd11Qd12*dtemp12 + Qd11Qd11*dtemp221013 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 1009 1029 ! In cos and sin values 1010 x = 0._wp1011 if ((ABS(dtemp1) > rsmall).or.(ABS(dtemp2) > rsmall)) then1012 x = atan2(dtemp2,dtemp1)1013 endif1030 zx = 0._wp 1031 IF ((ABS(dtemp1) > rsmall).OR.(ABS(dtemp2) > rsmall)) THEN 1032 zx = atan2(dtemp2,dtemp1) 1033 ENDIF 1014 1034 1015 1035 ! to ensure the angle lies between pi/4 and 9 pi/4 1016 if (x < rpi*0.25_wp) x =x + rpi*2.0_wp1036 IF (zx < rpi*0.25_wp) zx = zx + rpi*2.0_wp 1017 1037 1018 1038 ! y: angle between major principal axis of strain rate and structure … … 1022 1042 ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha) 1023 1043 1024 Tany_1 = Q11Q12 -Qd11Qd121025 Tany_2 = Q11Q11 -Qd12Qd121044 zTany_1 = zQ11Q12 - zQd11Qd12 1045 zTany_2 = zQ11Q11 - zQd12Qd12 1026 1046 1027 y = 0._wp1047 zy = 0._wp 1028 1048 1029 if ((ABS(Tany_1) > rsmall).or.(ABS(Tany_2) > rsmall)) then1030 y = atan2(Tany_1,Tany_2)1031 endif1049 IF ((ABS(zTany_1) > rsmall).OR.(ABS(zTany_2) > rsmall)) THEN 1050 zy = atan2(zTany_1,zTany_2) 1051 ENDIF 1032 1052 1033 1053 ! to make sure y is between 0 and pi 1034 if (y > rpi) y =y - rpi1035 if (y < 0) y =y + rpi1054 IF (zy > rpi) zy = zy - rpi 1055 IF (zy < 0) zy = zy + rpi 1036 1056 1037 1057 ! 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/ dx1042 invdy = 1._wp/ dy1043 invda = 1._wp/ da1058 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 1044 1064 1045 1065 ! % need 8 coords and 8 weights 1046 1066 ! % range in kx 1047 kx = int(( x-rpi*0.25_wp-rpi)*invdx) + 11048 kxw = kx - (x-rpi*0.25_wp-rpi)*invdx1067 kx = int((zx-rpi*0.25_wp-rpi)*invdx) + 1 1068 zkxw = kx - (zx-rpi*0.25_wp-rpi)*invdx 1049 1069 1050 ky = int( y*invdy) + 11051 kyw = ky - y*invdy1070 ky = int(zy*invdy) + 1 1071 kyw = ky - zy*invdy 1052 1072 1053 1073 ka = int((atempprime-0.5_wp)*invda) + 1 … … 1055 1075 1056 1076 ! % 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) 1081 1101 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) 1106 1126 1107 1127 ! Calculate mean ice stress over a collection of floes (Equation 3 in 1108 1128 ! Tsamados 2013) 1109 1129 1110 sig11 = strength*(stemp11r + kfriction*stemp11s) * invsin1111 sig12 = strength*(stemp12r + kfriction*stemp12s) * invsin1112 sig22 = strength*(stemp22r + kfriction*stemp22s) * invsin1130 zsig11 = strength*(zstemp11r + kfriction*zstemp11s) * invsin 1131 zsig12 = strength*(zstemp12r + kfriction*zstemp12s) * invsin 1132 zsig22 = strength*(zstemp22r + kfriction*zstemp22s) * invsin 1113 1133 1114 1134 !WRITE(numout,*) 'strength inside loop', strength … … 1118 1138 1119 1139 ! Update stress 1120 sgprm11 = Q11Q11*sig11 + Q12Q12*sig22 - 2._wp*Q11Q12 *sig121121 sgprm12 = Q11Q12*sig11 - Q11Q12*sig22 + (Q11Q11 - Q12Q12)*sig121122 sgprm22 = Q12Q12*sig11 + Q11Q11*sig22 + 2._wp*Q11Q12 *sig121123 1124 stressp = sgprm11 +sgprm221125 stress12 = sgprm121126 stressm = sgprm11 -sgprm221140 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 1127 1147 1128 1148 !WRITE(numout,*) 'stress output inside loop', ksub, stressp … … 1130 1150 ! Compute ridging and sliding functions in general coordinates 1131 1151 ! (Equation 11 in Tsamados 2013) 1132 if (ksub == ndte) then1133 rotstemp11r = Q11Q11*stemp11r - 2._wp*Q11Q12*stemp12r &1134 + Q12Q12*stemp22r1135 rotstemp12r = Q11Q11*stemp12r + Q11Q12*(stemp11r-stemp22r) &1136 - Q12Q12*stemp12r1137 rotstemp22r = Q12Q12*stemp11r + 2._wp*Q11Q12*stemp12r &1138 + Q11Q11*stemp22r1139 1140 rotstemp11s = Q11Q11*stemp11s - 2._wp*Q11Q12*stemp12s &1141 + Q12Q12*stemp22s1142 rotstemp12s = Q11Q11*stemp12s + Q11Q12*(stemp11s-stemp22s) &1143 - Q12Q12*stemp12s1144 rotstemp22s = Q12Q12*stemp11s + 2._wp*Q11Q12*stemp12s &1145 + Q11Q11*stemp22s1146 1147 alphar = rotstemp11r*dtemp11 + 2._wp*rotstemp12r*dtemp12 &1148 + rotstemp22r*dtemp221149 alphas = rotstemp11s*dtemp11 + 2._wp*rotstemp12s*dtemp12 &1150 + rotstemp22s*dtemp221151 endif1152 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 1152 1172 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 1153 1261 1154 1262 SUBROUTINE rhg_eap_rst( cdrw, kt ) … … 1229 1337 s12s(ix,iy,ia) = 0._wp 1230 1338 s22s(ix,iy,ia) = 0._wp 1231 IF (ia <= na_yield-1) then1339 IF (ia <= na_yield-1) THEN 1232 1340 DO iz=1,nz 1233 1341 s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* & … … 1299 1407 1300 1408 w1 = - 223.87569446_wp & 1301 1302 1303 1304 1305 1306 1307 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 1308 1416 1309 1417 END FUNCTION w1 … … 1319 1427 1320 1428 w2 = - 6670.68911883_wp & 1321 1322 1323 1324 1325 1326 1327 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 1328 1436 1329 1437 END FUNCTION w2 … … 1375 1483 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 1376 1484 1377 IF (-IIn1t2>=rsmall) then1485 IF (-IIn1t2>=rsmall) THEN 1378 1486 Hen1t2 = 1._wp 1379 1487 ELSE … … 1381 1489 ENDIF 1382 1490 1383 IF (-IIn2t1>=rsmall) then1491 IF (-IIn2t1>=rsmall) THEN 1384 1492 Hen2t1 = 1._wp 1385 1493 ELSE … … 1433 1541 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 1434 1542 1435 IF (-IIn1t2>=rsmall) then1543 IF (-IIn1t2>=rsmall) THEN 1436 1544 Hen1t2 = 1._wp 1437 1545 ELSE … … 1439 1547 ENDIF 1440 1548 1441 IF (-IIn2t1>=rsmall) then1549 IF (-IIn2t1>=rsmall) THEN 1442 1550 Hen2t1 = 1._wp 1443 1551 ELSE … … 1493 1601 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 1494 1602 1495 IF (-IIn1t2>=rsmall) then1603 IF (-IIn1t2>=rsmall) THEN 1496 1604 Hen1t2 = 1._wp 1497 1605 ELSE … … 1499 1607 ENDIF 1500 1608 1501 IF (-IIn2t1>=rsmall) then1609 IF (-IIn2t1>=rsmall) THEN 1502 1610 Hen2t1 = 1._wp 1503 1611 ELSE … … 1551 1659 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 1552 1660 1553 IF (-IIn1t2>=rsmall) then1661 IF (-IIn1t2>=rsmall) THEN 1554 1662 Hen1t2 = 1._wp 1555 1663 ELSE … … 1557 1665 ENDIF 1558 1666 1559 IF (-IIn2t1>=rsmall) then1667 IF (-IIn2t1>=rsmall) THEN 1560 1668 Hen2t1 = 1._wp 1561 1669 ELSE … … 1609 1717 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 1610 1718 1611 IF (-IIn1t2>=rsmall) then1719 IF (-IIn1t2>=rsmall) THEN 1612 1720 Hen1t2 = 1._wp 1613 1721 ELSE … … 1615 1723 ENDIF 1616 1724 1617 IF (-IIn2t1>=rsmall) then1725 IF (-IIn2t1>=rsmall) THEN 1618 1726 Hen2t1 = 1._wp 1619 1727 ELSE
Note: See TracChangeset
for help on using the changeset viewer.
