Changeset 13128
- Timestamp:
- 2020-06-18T16:55:01+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11842_SI3-10_EAP/src/ICE/icedyn_rhg_eap.F90
r13113 r13128 142 142 REAL(wp) :: zfac_x, zfac_y 143 143 REAL(wp) :: zshear, zdum1, zdum2 144 REAL(wp) :: stressptmp, stressmtmp,stress12tmpF ! anisotropic stress tensor components145 REAL(wp) :: alphar,alphas ! for mechanical redistribution144 REAL(wp) :: zstressptmp, zstressmtmp, zstress12tmpF ! anisotropic stress tensor components 145 REAL(wp) :: zalphar, zalphas ! for mechanical redistribution 146 146 REAL(wp) :: mresult11, mresult12, z1dtevpkth, p5kth, z1_dtevp_A ! for structure tensor evolution 147 147 ! 148 REAL(wp), DIMENSION(jpi,jpj) :: stress12tmp ! anisotropic stress tensor component for regridding148 REAL(wp), DIMENSION(jpi,jpj) :: zstress12tmp ! anisotropic stress tensor component for regridding 149 149 REAL(wp), DIMENSION(jpi,jpj) :: yield11, yield22, yield12 ! yield surface tensor for history 150 150 REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points … … 415 415 CALL update_stress_rdg (jter, nn_nevp, phi, zdiv, zdt, & 416 416 zdsT, paniso_11(ji,jj), paniso_12(ji,jj), & 417 stressptmp,stressmtmp, &418 stress12tmp(ji,jj), strength(ji,jj), &419 alphar,alphas)417 zstressptmp, zstressmtmp, & 418 zstress12tmp(ji,jj), strength(ji,jj), & 419 zalphar, zalphas) 420 420 421 421 ! structure tensor update 422 422 IF (mod(jter,10) == 0) THEN 423 CALL calc_ffrac( stressptmp, stressmtmp,stress12tmp(ji,jj), &423 CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), & 424 424 paniso_11(ji,jj), paniso_12(ji,jj), & 425 425 mresult11, mresult12) … … 431 431 432 432 IF (jter == nn_nevp) THEN 433 yield11(ji,jj) = 0.5_wp * ( stressptmp +stressmtmp)434 yield22(ji,jj) = 0.5_wp * ( stressptmp -stressmtmp)435 yield12(ji,jj) = stress12tmp(ji,jj)436 prdg_conv(ji,jj) = -min( alphar, 0._wp)433 yield11(ji,jj) = 0.5_wp * (zstressptmp + zstressmtmp) 434 yield22(ji,jj) = 0.5_wp * (zstressptmp - zstressmtmp) 435 yield12(ji,jj) = zstress12tmp(ji,jj) 436 prdg_conv(ji,jj) = -min( zalphar, 0._wp) 437 437 ENDIF 438 438 … … 454 454 !zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 455 455 !zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 456 zs1(ji,jj) = ( zs1(ji,jj) + stressptmp * zalph1 ) * z1_alph1457 zs2(ji,jj) = ( zs2(ji,jj) + stressmtmp * zalph1 ) * z1_alph1456 zs1(ji,jj) = ( zs1(ji,jj) + zstressptmp * zalph1 ) * z1_alph1 457 zs2(ji,jj) = ( zs2(ji,jj) + zstressmtmp * zalph1 ) * z1_alph1 458 458 !zs1(ji,jj) = ( stressptmp * zs1(ji,jj) + zalph1 ) * z1_alph1 459 459 !zs2(ji,jj) = ( stressmtmp * zs2(ji,jj) + zalph1 ) * z1_alph1 … … 462 462 END DO 463 463 !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. ) 464 CALL lbc_lnk_multi( 'icedyn_rhg_eap', stress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.)464 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zstress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.) 465 465 466 466 DO jj = 1, jpjm1 467 467 DO ji = 1, jpim1 468 468 ! stress12tmp at F points 469 stress12tmpF = ( stress12tmp(ji,jj+1) * e1e2t(ji,jj+1) +stress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1) &470 & + stress12tmp(ji,jj ) * e1e2t(ji,jj ) +stress12tmp(ji+1,jj ) * e1e2t(ji+1,jj ) &469 zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1) & 470 & + zstress12tmp(ji,jj ) * e1e2t(ji,jj ) + zstress12tmp(ji+1,jj ) * e1e2t(ji+1,jj ) & 471 471 & ) * 0.25_wp * r1_e1e2f(ji,jj) 472 472 … … 479 479 ! stress at F points (zkt/=0 if landfast) 480 480 !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 481 zs12(ji,jj) = ( zs12(ji,jj) + stress12tmpF * zalph1 ) * z1_alph1481 zs12(ji,jj) = ( zs12(ji,jj) + zstress12tmpF * zalph1 ) * z1_alph1 482 482 !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1 ) * z1_alph1 483 483 … … 527 527 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 528 528 zTauB = ztauy_base(ji,jj) / zvel 529 ! !--- OceanBottom-to-Ice stress 529 ! !--- OceanBottom-to-Ice stress 530 530 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 531 531 ! … … 572 572 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 573 573 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 574 ! !--- Ocean-to-Ice stress 574 ! !--- Ocean-to-Ice stress 575 575 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 576 576 ! … … 625 625 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 626 626 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 627 ! !--- Ocean-to-Ice stress 627 ! !--- Ocean-to-Ice stress 628 628 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 629 629 ! … … 841 841 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 842 842 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 843 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 843 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 844 844 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 845 845 END DO … … 853 853 ! Stress tensor invariants (normal and shear stress N/m) 854 854 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 855 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 855 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 856 856 857 857 DEALLOCATE( zsig1 , zsig2 , zsig3 ) … … 935 935 stressp, stressm, & 936 936 stress12, strength, & 937 alphar,alphas)937 palphar, palphas) 938 938 !!--------------------------------------------------------------------- 939 939 !! *** SUBROUTINE rhg_eap_rst *** 940 940 !! 941 !! ** Purpose : Updates the stressdepending on values of strain rate and structure941 !! ** Purpose : Updates the depending on values of strain rate and structure 942 942 !! tensor and for ksub=ndte it computes closing and sliding rate 943 943 !!--------------------------------------------------------------------- … … 947 947 REAL(wp), INTENT(in ) :: a11, a12 948 948 REAL(wp), INTENT( out) :: stressp, stressm, stress12 949 REAL(wp), INTENT( out) :: alphar,alphas949 REAL(wp), INTENT( out) :: palphar, palphas 950 950 951 951 INTEGER :: kx ,ky, ka … … 964 964 REAL(wp) :: zTany_1, zTany_2 965 965 REAL(wp) :: zx, zy, zdx, zdy, zda, zkxw, kyw, kaw 966 REAL(wp) :: invdx, invdy,invda967 REAL(wp) :: dtemp1, dtemp2, atempprime, a,invsin966 REAL(wp) :: zinvdx, zinvdy, zinvda 967 REAL(wp) :: zdtemp1, zdtemp2, zatempprime, zinvsin 968 968 969 969 REAL(wp), PARAMETER :: kfriction = 0.45_wp … … 973 973 zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 974 974 975 invsin = 1._wp/sin(2._wp*phi) * zinvstressconviso975 zinvsin = 1._wp/sin(2._wp*phi) * zinvstressconviso 976 976 !now uses phi as set in higher code 977 977 … … 997 997 998 998 ! rotation Q*atemp*Q^T 999 atempprime = zQ11Q11*a11 + 2.0_wp*zQ11Q12*a12 + zQ12Q12*za22999 zatempprime = zQ11Q11*a11 + 2.0_wp*zQ11Q12*a12 + zQ12Q12*za22 1000 1000 1001 1001 ! make first principal value the largest 1002 atempprime = max(atempprime, 1.0_wp -atempprime)1002 zatempprime = max(zatempprime, 1.0_wp - zatempprime) 1003 1003 1004 1004 ! 2) strain rate … … 1025 1025 ENDIF 1026 1026 1027 dtemp1 = zQd11Qd11*zdtemp11 + 2.0_wp*zQd11Qd12*zdtemp12 + zQd12Qd12*zdtemp221028 dtemp2 = zQd12Qd12*zdtemp11 - 2.0_wp*zQd11Qd12*zdtemp12 + zQd11Qd11*zdtemp221027 zdtemp1 = zQd11Qd11*zdtemp11 + 2.0_wp*zQd11Qd12*zdtemp12 + zQd12Qd12*zdtemp22 1028 zdtemp2 = zQd12Qd12*zdtemp11 - 2.0_wp*zQd11Qd12*zdtemp12 + zQd11Qd11*zdtemp22 1029 1029 ! In cos and sin values 1030 1030 zx = 0._wp 1031 IF ((ABS( dtemp1) > rsmall).OR.(ABS(dtemp2) > rsmall)) THEN1032 zx = atan2( dtemp2,dtemp1)1031 IF ((ABS(zdtemp1) > rsmall).OR.(ABS(zdtemp2) > rsmall)) THEN 1032 zx = atan2(zdtemp2,zdtemp1) 1033 1033 ENDIF 1034 1034 … … 1059 1059 zdy = rpi/real(ny_yield-1,kind=wp) 1060 1060 zda = 0.5_wp/real(na_yield-1,kind=wp) 1061 invdx = 1._wp/zdx1062 invdy = 1._wp/zdy1063 invda = 1._wp/zda1061 zinvdx = 1._wp/zdx 1062 zinvdy = 1._wp/zdy 1063 zinvda = 1._wp/zda 1064 1064 1065 1065 ! % need 8 coords and 8 weights 1066 1066 ! % range in kx 1067 kx = int((zx-rpi*0.25_wp-rpi)* invdx) + 11068 zkxw = kx - (zx-rpi*0.25_wp-rpi)* invdx1067 kx = int((zx-rpi*0.25_wp-rpi)*zinvdx) + 1 1068 zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx 1069 1069 1070 ky = int(zy* invdy) + 11071 kyw = ky - zy* invdy1070 ky = int(zy*zinvdy) + 1 1071 kyw = ky - zy*zinvdy 1072 1072 1073 ka = int(( atempprime-0.5_wp)*invda) + 11074 kaw = ka - ( atempprime-0.5_wp)*invda1073 ka = int((zatempprime-0.5_wp)*zinvda) + 1 1074 kaw = ka - (zatempprime-0.5_wp)*zinvda 1075 1075 1076 1076 ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013) … … 1128 1128 ! Tsamados 2013) 1129 1129 1130 zsig11 = strength*(zstemp11r + kfriction*zstemp11s) * invsin1131 zsig12 = strength*(zstemp12r + kfriction*zstemp12s) * invsin1132 zsig22 = strength*(zstemp22r + kfriction*zstemp22s) * invsin1130 zsig11 = strength*(zstemp11r + kfriction*zstemp11s) * zinvsin 1131 zsig12 = strength*(zstemp12r + kfriction*zstemp12s) * zinvsin 1132 zsig22 = strength*(zstemp22r + kfriction*zstemp22s) * zinvsin 1133 1133 1134 1134 !WRITE(numout,*) 'strength inside loop', strength … … 1146 1146 stressm = zsgprm11 - zsgprm22 1147 1147 1148 !WRITE(numout,*) ' stressoutput inside loop', ksub, stressp1148 !WRITE(numout,*) ' output inside loop', ksub, stressp 1149 1149 1150 1150 ! Compute ridging and sliding functions in general coordinates … … 1165 1165 + zQ11Q11*zstemp22s 1166 1166 1167 alphar = zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 &1167 palphar = zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 & 1168 1168 + zrotstemp22r*zdtemp22 1169 alphas = zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 &1169 palphas = zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 & 1170 1170 + zrotstemp22s*zdtemp22 1171 1171 ENDIF … … 1190 1190 ! local variables 1191 1191 1192 REAL (wp) :: sigma11, sigma12,sigma22 ! stress tensor elements1193 REAL (wp) :: Angle_denom ! angle with principal component axis1194 REAL (wp) :: sigma_1,sigma_2 ! principal components of stress1195 REAL (wp) :: Q11, Q12, Q11Q11, Q11Q12,Q12Q121192 REAL (wp) :: zsigma11, zsigma12, zsigma22 ! stress tensor elements 1193 REAL (wp) :: zAngle_denom ! angle with principal component axis 1194 REAL (wp) :: zsigma_1, zsigma_2 ! principal components of stress 1195 REAL (wp) :: zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12 1196 1196 1197 1197 REAL (wp), PARAMETER :: kfrac = 0.0001_wp ! rate of fracture formation … … 1199 1199 !!--------------------------------------------------------------- 1200 1200 ! 1201 sigma11 = 0.5_wp*(stressp+stressm)1202 sigma12 = stress121203 sigma22 = 0.5_wp*(stressp-stressm)1201 zsigma11 = 0.5_wp*(stressp+stressm) 1202 zsigma12 = stress12 1203 zsigma22 = 0.5_wp*(stressp-stressm) 1204 1204 1205 1205 ! Here's the change - no longer calculate gamma, … … 1212 1212 ! error to the calculated angles < rsmall 1213 1213 1214 Q11Q11 = 0.1_wp1215 Q12Q12 = rsmall1216 Q11Q12 = rsmall1217 1218 IF((ABS( sigma11 - sigma22) > rsmall).OR.(ABS(sigma12) > rsmall)) THEN1219 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^21224 Q12Q12 = 0.5_wp - ( sigma11 - sigma22 )*0.5_wp*Angle_denom !Sin^21225 Q11Q12 = sigma12*Angle_denom !CosSin1226 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)1214 zQ11Q11 = 0.1_wp 1215 zQ12Q12 = rsmall 1216 zQ11Q12 = rsmall 1217 1218 IF((ABS( zsigma11 - zsigma22) > rsmall).OR.(ABS(zsigma12) > rsmall)) THEN 1219 1220 zAngle_denom = 1.0_wp/sqrt( ( zsigma11 - zsigma22 )*( zsigma11 - & 1221 zsigma22 ) + 4.0_wp*zsigma12*zsigma12) 1222 1223 zQ11Q11 = 0.5_wp + ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom !Cos^2 1224 zQ12Q12 = 0.5_wp - ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom !Sin^2 1225 zQ11Q12 = zsigma12*zAngle_denom !CosSin 1226 ENDIF 1227 1228 zsigma_1 = zQ11Q11*zsigma11 + 2.0_wp*zQ11Q12*zsigma12 & 1229 + zQ12Q12*zsigma22 ! S(1,1) 1230 zsigma_2 = zQ12Q12*zsigma11 - 2.0_wp*zQ11Q12*zsigma12 & 1231 + zQ11Q11*zsigma22 ! S(2,2) 1232 1232 1233 1233 ! Pure divergence 1234 IF (( sigma_1 >= 0.0_wp).AND.(sigma_2 >= 0.0_wp)) THEN1234 IF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 >= 0.0_wp)) THEN 1235 1235 mresult11 = 0.0_wp 1236 1236 mresult12 = 0.0_wp … … 1239 1239 ! direction 1240 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)) THEN1242 mresult11 = kfrac * (a11 - Q12Q12)1243 mresult12 = kfrac * (a12 + Q11Q12)1241 ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp)) THEN 1242 mresult11 = kfrac * (a11 - zQ12Q12) 1243 mresult12 = kfrac * (a12 + zQ11Q12) 1244 1244 1245 1245 ! Shear faulting 1246 ELSEIF ( sigma_2 == 0.0_wp) THEN1246 ELSEIF (zsigma_2 == 0.0_wp) THEN 1247 1247 mresult11 = 0.0_wp 1248 1248 mresult12 = 0.0_wp 1249 ELSEIF (( sigma_1 <= 0.0_wp).AND.(sigma_1/sigma_2 <= threshold)) THEN1250 mresult11 = kfrac * (a11 - Q12Q12)1251 mresult12 = kfrac * (a12 + Q11Q12)1249 ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN 1250 mresult11 = kfrac * (a11 - zQ12Q12) 1251 mresult12 = kfrac * (a12 + zQ11Q12) 1252 1252 1253 1253 ! Horizontal spalling … … 1398 1398 1399 1399 1400 FUNCTION w1( a)1400 FUNCTION w1(pa) 1401 1401 !!------------------------------------------------------------------- 1402 1402 !! Function : w1 (see Gaussian function psi in Tsamados et al 2013) 1403 1403 !!------------------------------------------------------------------- 1404 REAL(wp), INTENT(in ) :: a1404 REAL(wp), INTENT(in ) :: pa 1405 1405 REAL(wp) :: w1 1406 1406 !!------------------------------------------------------------------- 1407 1407 1408 1408 w1 = - 223.87569446_wp & 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*a1409 & + 2361.2198663_wp*pa & 1410 & - 10606.56079975_wp*pa*pa & 1411 & + 26315.50025642_wp*pa*pa*pa & 1412 & - 38948.30444297_wp*pa*pa*pa*pa & 1413 & + 34397.72407466_wp*pa*pa*pa*pa*pa & 1414 & - 16789.98003081_wp*pa*pa*pa*pa*pa*pa & 1415 & + 3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa 1416 1416 1417 1417 END FUNCTION w1 1418 1418 1419 1419 1420 FUNCTION w2( a)1420 FUNCTION w2(pa) 1421 1421 !!------------------------------------------------------------------- 1422 1422 !! Function : w2 (see Gaussian function psi in Tsamados et al 2013) 1423 1423 !!------------------------------------------------------------------- 1424 REAL(wp), INTENT(in ) :: a1424 REAL(wp), INTENT(in ) :: pa 1425 1425 REAL(wp) :: w2 1426 1426 !!------------------------------------------------------------------- 1427 1427 1428 1428 w2 = - 6670.68911883_wp & 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*a1429 & + 70222.33061536_wp*pa & 1430 & - 314871.71525448_wp*pa*pa & 1431 & + 779570.02793492_wp*pa*pa*pa & 1432 & - 1151098.82436864_wp*pa*pa*pa*pa & 1433 & + 1013896.59464498_wp*pa*pa*pa*pa*pa & 1434 & - 493379.44906738_wp*pa*pa*pa*pa*pa*pa & 1435 & + 102356.551518_wp*pa*pa*pa*pa*pa*pa*pa 1436 1436 1437 1437 END FUNCTION w2 1438 1438 1439 FUNCTION s11kr( x,y,z,phi)1439 FUNCTION s11kr(px,py,pz,phi) 1440 1440 !!------------------------------------------------------------------- 1441 1441 !! Function : s11kr 1442 1442 !!------------------------------------------------------------------- 1443 REAL(wp), INTENT(in ) :: x,y,z,phi 1444 1445 REAL(wp) :: s11kr, p, pih 1446 1443 REAL(wp), INTENT(in ) :: px,py,pz,phi 1444 1445 REAL(wp) :: s11kr, zpih 1446 1447 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1448 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1449 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1450 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1451 REAL(wp) :: zd11, zd12, zd22 1452 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1453 REAL(wp) :: zHen1t2, zHen2t1 1454 !!------------------------------------------------------------------- 1455 1456 zpih = 0.5_wp*rpi 1457 1458 zn1t2i11 = cos(pz+zpih-phi) * cos(pz+phi) 1459 zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi) 1460 zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi) 1461 zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi) 1462 zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi) 1463 zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi) 1464 zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi) 1465 zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi) 1466 zt1t2i11 = cos(pz-phi) * cos(pz+phi) 1467 zt1t2i12 = cos(pz-phi) * sin(pz+phi) 1468 zt1t2i21 = sin(pz-phi) * cos(pz+phi) 1469 zt1t2i22 = sin(pz-phi) * sin(pz+phi) 1470 zt2t1i11 = cos(pz+phi) * cos(pz-phi) 1471 zt2t1i12 = cos(pz+phi) * sin(pz-phi) 1472 zt2t1i21 = sin(pz+phi) * cos(pz-phi) 1473 zt2t1i22 = sin(pz+phi) * sin(pz-phi) 1474 ! In expression of tensor d, with this formulatin d(x)=-d(x+pi) 1475 ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else 1476 ! x=x-pi/2 1477 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1478 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1479 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1480 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1481 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1482 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1483 1484 IF (-zIIn1t2>=rsmall) THEN 1485 zHen1t2 = 1._wp 1486 ELSE 1487 zHen1t2 = 0._wp 1488 ENDIF 1489 1490 IF (-zIIn2t1>=rsmall) THEN 1491 zHen2t1 = 1._wp 1492 ELSE 1493 zHen2t1 = 0._wp 1494 ENDIF 1495 1496 s11kr = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11) 1497 1498 END FUNCTION s11kr 1499 1500 FUNCTION s12kr(px,py,pz,phi) 1501 !!------------------------------------------------------------------- 1502 !! Function : s12kr 1503 !!------------------------------------------------------------------- 1504 REAL(wp), INTENT(in ) :: px,py,pz,phi 1505 1506 REAL(wp) :: s12kr, zs12r0, zs21r0, zpih 1507 1508 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1509 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1510 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1511 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1512 REAL(wp) :: zd11, zd12, zd22 1513 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1514 REAL(wp) :: zHen1t2, zHen2t1 1515 !!------------------------------------------------------------------- 1516 zpih = 0.5_wp*rpi 1517 1518 zn1t2i11 = cos(pz+zpih-phi) * cos(pz+phi) 1519 zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi) 1520 zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi) 1521 zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi) 1522 zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi) 1523 zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi) 1524 zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi) 1525 zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi) 1526 zt1t2i11 = cos(pz-phi) * cos(pz+phi) 1527 zt1t2i12 = cos(pz-phi) * sin(pz+phi) 1528 zt1t2i21 = sin(pz-phi) * cos(pz+phi) 1529 zt1t2i22 = sin(pz-phi) * sin(pz+phi) 1530 zt2t1i11 = cos(pz+phi) * cos(pz-phi) 1531 zt2t1i12 = cos(pz+phi) * sin(pz-phi) 1532 zt2t1i21 = sin(pz+phi) * cos(pz-phi) 1533 zt2t1i22 = sin(pz+phi) * sin(pz-phi) 1534 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1535 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1536 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1537 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1538 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1539 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1540 1541 IF (-zIIn1t2>=rsmall) THEN 1542 zHen1t2 = 1._wp 1543 ELSE 1544 zHen1t2 = 0._wp 1545 ENDIF 1546 1547 IF (-zIIn2t1>=rsmall) THEN 1548 zHen2t1 = 1._wp 1549 ELSE 1550 zHen2t1 = 0._wp 1551 ENDIF 1552 1553 zs12r0 = (- zHen1t2 * zn1t2i12 - zHen2t1 * zn2t1i12) 1554 zs21r0 = (- zHen1t2 * zn1t2i21 - zHen2t1 * zn2t1i21) 1555 s12kr=0.5_wp*(zs12r0+zs21r0) 1556 1557 END FUNCTION s12kr 1558 1559 FUNCTION s22kr(px,py,pz,phi) 1560 !!------------------------------------------------------------------- 1561 !! Function : s22kr 1562 !!------------------------------------------------------------------- 1563 REAL(wp), INTENT(in ) :: px,py,pz,phi 1564 1565 REAL(wp) :: s22kr, zpih 1566 1567 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1568 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1569 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1570 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1571 REAL(wp) :: zd11, zd12, zd22 1572 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1573 REAL(wp) :: zHen1t2, zHen2t1 1574 !!------------------------------------------------------------------- 1575 zpih = 0.5_wp*rpi 1576 1577 zn1t2i11 = cos(pz+zpih-phi) * cos(pz+phi) 1578 zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi) 1579 zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi) 1580 zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi) 1581 zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi) 1582 zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi) 1583 zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi) 1584 zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi) 1585 zt1t2i11 = cos(pz-phi) * cos(pz+phi) 1586 zt1t2i12 = cos(pz-phi) * sin(pz+phi) 1587 zt1t2i21 = sin(pz-phi) * cos(pz+phi) 1588 zt1t2i22 = sin(pz-phi) * sin(pz+phi) 1589 zt2t1i11 = cos(pz+phi) * cos(pz-phi) 1590 zt2t1i12 = cos(pz+phi) * sin(pz-phi) 1591 zt2t1i21 = sin(pz+phi) * cos(pz-phi) 1592 zt2t1i22 = sin(pz+phi) * sin(pz-phi) 1593 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1594 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1595 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1596 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1597 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1598 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1599 1600 IF (-zIIn1t2>=rsmall) THEN 1601 zHen1t2 = 1._wp 1602 ELSE 1603 zHen1t2 = 0._wp 1604 ENDIF 1605 1606 IF (-zIIn2t1>=rsmall) THEN 1607 zHen2t1 = 1._wp 1608 ELSE 1609 zHen2t1 = 0._wp 1610 ENDIF 1611 1612 s22kr = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22) 1613 1614 END FUNCTION s22kr 1615 1616 FUNCTION s11ks(px,py,pz,phi) 1617 !!------------------------------------------------------------------- 1618 !! Function : s11ks 1619 !!------------------------------------------------------------------- 1620 REAL(wp), INTENT(in ) :: px,py,pz,phi 1621 1622 REAL(wp) :: s11ks, zpih 1623 1624 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1625 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1626 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1627 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1628 REAL(wp) :: zd11, zd12, zd22 1629 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1630 REAL(wp) :: zHen1t2, zHen2t1 1631 !!------------------------------------------------------------------- 1632 zpih = 0.5_wp*rpi 1633 1634 zn1t2i11 = cos(pz+zpih-phi) * cos(pz+phi) 1635 zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi) 1636 zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi) 1637 zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi) 1638 zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi) 1639 zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi) 1640 zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi) 1641 zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi) 1642 zt1t2i11 = cos(pz-phi) * cos(pz+phi) 1643 zt1t2i12 = cos(pz-phi) * sin(pz+phi) 1644 zt1t2i21 = sin(pz-phi) * cos(pz+phi) 1645 zt1t2i22 = sin(pz-phi) * sin(pz+phi) 1646 zt2t1i11 = cos(pz+phi) * cos(pz-phi) 1647 zt2t1i12 = cos(pz+phi) * sin(pz-phi) 1648 zt2t1i21 = sin(pz+phi) * cos(pz-phi) 1649 zt2t1i22 = sin(pz+phi) * sin(pz-phi) 1650 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1651 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1652 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1653 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1654 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1655 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1656 1657 IF (-zIIn1t2>=rsmall) THEN 1658 zHen1t2 = 1._wp 1659 ELSE 1660 zHen1t2 = 0._wp 1661 ENDIF 1662 1663 IF (-zIIn2t1>=rsmall) THEN 1664 zHen2t1 = 1._wp 1665 ELSE 1666 zHen2t1 = 0._wp 1667 ENDIF 1668 1669 s11ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11) 1670 1671 END FUNCTION s11ks 1672 1673 FUNCTION s12ks(px,py,pz,phi) 1674 !!------------------------------------------------------------------- 1675 !! Function : s12ks 1676 !!------------------------------------------------------------------- 1677 REAL(wp), INTENT(in ) :: px,py,pz,phi 1678 1679 REAL(wp) :: s12ks, zs12s0, zs21s0, zpih 1680 1681 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1682 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1683 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1684 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1685 REAL(wp) :: zd11, zd12, zd22 1686 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1687 REAL(wp) :: zHen1t2, zHen2t1 1688 !!------------------------------------------------------------------- 1689 zpih = 0.5_wp*rpi 1690 1691 zn1t2i11 = cos(pz+zpih-phi) * cos(pz+phi) 1692 zn1t2i12 = cos(pz+zpih-phi) * sin(pz+phi) 1693 zn1t2i21 = sin(pz+zpih-phi) * cos(pz+phi) 1694 zn1t2i22 = sin(pz+zpih-phi) * sin(pz+phi) 1695 zn2t1i11 = cos(pz-zpih+phi) * cos(pz-phi) 1696 zn2t1i12 = cos(pz-zpih+phi) * sin(pz-phi) 1697 zn2t1i21 = sin(pz-zpih+phi) * cos(pz-phi) 1698 zn2t1i22 = sin(pz-zpih+phi) * sin(pz-phi) 1699 zt1t2i11 = cos(pz-phi) * cos(pz+phi) 1700 zt1t2i12 = cos(pz-phi) * sin(pz+phi) 1701 zt1t2i21 = sin(pz-phi) * cos(pz+phi) 1702 zt1t2i22 = sin(pz-phi) * sin(pz+phi) 1703 zt2t1i11 = cos(pz+phi) * cos(pz-phi) 1704 zt2t1i12 = cos(pz+phi) * sin(pz-phi) 1705 zt2t1i21 = sin(pz+phi) * cos(pz-phi) 1706 zt2t1i22 = sin(pz+phi) * sin(pz-phi) 1707 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1708 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1709 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1710 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1711 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1712 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1713 1714 IF (-zIIn1t2>=rsmall) THEN 1715 zHen1t2 = 1._wp 1716 ELSE 1717 zHen1t2 = 0._wp 1718 ENDIF 1719 1720 IF (-zIIn2t1>=rsmall) THEN 1721 zHen2t1 = 1._wp 1722 ELSE 1723 zHen2t1 = 0._wp 1724 ENDIF 1725 1726 zs12s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i12 + zHen2t1 * zt2t1i12) 1727 zs21s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i21 + zHen2t1 * zt2t1i21) 1728 s12ks=0.5_wp*(zs12s0+zs21s0) 1729 1730 END FUNCTION s12ks 1731 1732 FUNCTION s22ks(px,py,pz,phi) 1733 !!------------------------------------------------------------------- 1734 !! Function : s22ks 1735 !!------------------------------------------------------------------- 1736 REAL(wp), INTENT(in ) :: px,py,pz,phi 1737 1738 REAL(wp) :: s22ks, zpih 1739 1447 1740 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21, n1t2i22 1448 1741 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21, n2t1i22 … … 1453 1746 REAL(wp) :: Hen1t2, Hen2t1 1454 1747 !!------------------------------------------------------------------- 1455 1456 pih = 0.5_wp*rpi 1457 p = phi 1458 1459 n1t2i11 = cos(z+pih-p) * cos(z+p) 1460 n1t2i12 = cos(z+pih-p) * sin(z+p) 1461 n1t2i21 = sin(z+pih-p) * cos(z+p) 1462 n1t2i22 = sin(z+pih-p) * sin(z+p) 1463 n2t1i11 = cos(z-pih+p) * cos(z-p) 1464 n2t1i12 = cos(z-pih+p) * sin(z-p) 1465 n2t1i21 = sin(z-pih+p) * cos(z-p) 1466 n2t1i22 = sin(z-pih+p) * sin(z-p) 1467 t1t2i11 = cos(z-p) * cos(z+p) 1468 t1t2i12 = cos(z-p) * sin(z+p) 1469 t1t2i21 = sin(z-p) * cos(z+p) 1470 t1t2i22 = sin(z-p) * sin(z+p) 1471 t2t1i11 = cos(z+p) * cos(z-p) 1472 t2t1i12 = cos(z+p) * sin(z-p) 1473 t2t1i21 = sin(z+p) * cos(z-p) 1474 t2t1i22 = sin(z+p) * sin(z-p) 1475 ! In expression of tensor d, with this formulatin d(x)=-d(x+pi) 1476 ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else 1477 ! x=x-pi/2 1478 d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 1479 d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 1480 d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 1748 zpih = 0.5_wp*rpi 1749 1750 n1t2i11 = cos(pz+zpih-phi) * cos(pz+phi) 1751 n1t2i12 = cos(pz+zpih-phi) * sin(pz+phi) 1752 n1t2i21 = sin(pz+zpih-phi) * cos(pz+phi) 1753 n1t2i22 = sin(pz+zpih-phi) * sin(pz+phi) 1754 n2t1i11 = cos(pz-zpih+phi) * cos(pz-phi) 1755 n2t1i12 = cos(pz-zpih+phi) * sin(pz-phi) 1756 n2t1i21 = sin(pz-zpih+phi) * cos(pz-phi) 1757 n2t1i22 = sin(pz-zpih+phi) * sin(pz-phi) 1758 t1t2i11 = cos(pz-phi) * cos(pz+phi) 1759 t1t2i12 = cos(pz-phi) * sin(pz+phi) 1760 t1t2i21 = sin(pz-phi) * cos(pz+phi) 1761 t1t2i22 = sin(pz-phi) * sin(pz+phi) 1762 t2t1i11 = cos(pz+phi) * cos(pz-phi) 1763 t2t1i12 = cos(pz+phi) * sin(pz-phi) 1764 t2t1i21 = sin(pz+phi) * cos(pz-phi) 1765 t2t1i22 = sin(pz+phi) * sin(pz-phi) 1766 d11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1767 d12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1768 d22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1481 1769 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 1482 1770 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 1483 1771 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 1484 1772 1485 1773 IF (-IIn1t2>=rsmall) THEN 1486 1774 Hen1t2 = 1._wp … … 1488 1776 Hen1t2 = 0._wp 1489 1777 ENDIF 1490 1778 1491 1779 IF (-IIn2t1>=rsmall) THEN 1492 1780 Hen2t1 = 1._wp … … 1494 1782 Hen2t1 = 0._wp 1495 1783 ENDIF 1496 1497 s11kr = (- Hen1t2 * n1t2i11 - Hen2t1 * n2t1i11)1498 1499 END FUNCTION s11kr1500 1501 FUNCTION s12kr(x,y,z,phi)1502 !!-------------------------------------------------------------------1503 !! Function : s12kr1504 !!-------------------------------------------------------------------1505 REAL(wp), INTENT(in ) :: x,y,z,phi1506 1507 REAL(wp) :: s12kr, s12r0, s21r0, p, pih1508 1509 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21, n1t2i221510 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21, n2t1i221511 REAL(wp) :: t1t2i11, t1t2i12, t1t2i21, t1t2i221512 REAL(wp) :: t2t1i11, t2t1i12, t2t1i21, t2t1i221513 REAL(wp) :: d11, d12, d221514 REAL(wp) :: IIn1t2, IIn2t1, IIt1t21515 REAL(wp) :: Hen1t2, Hen2t11516 !!-------------------------------------------------------------------1517 pih = 0.5_wp*rpi1518 p = phi1519 1520 n1t2i11 = cos(z+pih-p) * cos(z+p)1521 n1t2i12 = cos(z+pih-p) * sin(z+p)1522 n1t2i21 = sin(z+pih-p) * cos(z+p)1523 n1t2i22 = sin(z+pih-p) * sin(z+p)1524 n2t1i11 = cos(z-pih+p) * cos(z-p)1525 n2t1i12 = cos(z-pih+p) * sin(z-p)1526 n2t1i21 = sin(z-pih+p) * cos(z-p)1527 n2t1i22 = sin(z-pih+p) * sin(z-p)1528 t1t2i11 = cos(z-p) * cos(z+p)1529 t1t2i12 = cos(z-p) * sin(z+p)1530 t1t2i21 = sin(z-p) * cos(z+p)1531 t1t2i22 = sin(z-p) * sin(z+p)1532 t2t1i11 = cos(z+p) * cos(z-p)1533 t2t1i12 = cos(z+p) * sin(z-p)1534 t2t1i21 = sin(z+p) * cos(z-p)1535 t2t1i22 = sin(z+p) * sin(z-p)1536 d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))1537 d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))1538 d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))1539 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d221540 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d221541 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d221542 1543 IF (-IIn1t2>=rsmall) THEN1544 Hen1t2 = 1._wp1545 ELSE1546 Hen1t2 = 0._wp1547 ENDIF1548 1549 IF (-IIn2t1>=rsmall) THEN1550 Hen2t1 = 1._wp1551 ELSE1552 Hen2t1 = 0._wp1553 ENDIF1554 1555 s12r0 = (- Hen1t2 * n1t2i12 - Hen2t1 * n2t1i12)1556 s21r0 = (- Hen1t2 * n1t2i21 - Hen2t1 * n2t1i21)1557 s12kr=0.5_wp*(s12r0+s21r0)1558 1559 END FUNCTION s12kr1560 1561 FUNCTION s22kr(x,y,z,phi)1562 !!-------------------------------------------------------------------1563 !! Function : s22kr1564 !!-------------------------------------------------------------------1565 REAL(wp), INTENT(in ) :: x,y,z,phi1566 1567 REAL(wp) :: s22kr, p, pih1568 1569 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21, n1t2i221570 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21, n2t1i221571 REAL(wp) :: t1t2i11, t1t2i12, t1t2i21, t1t2i221572 REAL(wp) :: t2t1i11, t2t1i12, t2t1i21, t2t1i221573 REAL(wp) :: d11, d12, d221574 REAL(wp) :: IIn1t2, IIn2t1, IIt1t21575 REAL(wp) :: Hen1t2, Hen2t11576 !!-------------------------------------------------------------------1577 pih = 0.5_wp*rpi1578 p = phi1579 1580 n1t2i11 = cos(z+pih-p) * cos(z+p)1581 n1t2i12 = cos(z+pih-p) * sin(z+p)1582 n1t2i21 = sin(z+pih-p) * cos(z+p)1583 n1t2i22 = sin(z+pih-p) * sin(z+p)1584 n2t1i11 = cos(z-pih+p) * cos(z-p)1585 n2t1i12 = cos(z-pih+p) * sin(z-p)1586 n2t1i21 = sin(z-pih+p) * cos(z-p)1587 n2t1i22 = sin(z-pih+p) * sin(z-p)1588 t1t2i11 = cos(z-p) * cos(z+p)1589 t1t2i12 = cos(z-p) * sin(z+p)1590 t1t2i21 = sin(z-p) * cos(z+p)1591 t1t2i22 = sin(z-p) * sin(z+p)1592 t2t1i11 = cos(z+p) * cos(z-p)1593 t2t1i12 = cos(z+p) * sin(z-p)1594 t2t1i21 = sin(z+p) * cos(z-p)1595 t2t1i22 = sin(z+p) * sin(z-p)1596 d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))1597 d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))1598 d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))1599 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d221600 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d221601 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d221602 1603 IF (-IIn1t2>=rsmall) THEN1604 Hen1t2 = 1._wp1605 ELSE1606 Hen1t2 = 0._wp1607 ENDIF1608 1609 IF (-IIn2t1>=rsmall) THEN1610 Hen2t1 = 1._wp1611 ELSE1612 Hen2t1 = 0._wp1613 ENDIF1614 1615 s22kr = (- Hen1t2 * n1t2i22 - Hen2t1 * n2t1i22)1616 1617 END FUNCTION s22kr1618 1619 FUNCTION s11ks(x,y,z,phi)1620 !!-------------------------------------------------------------------1621 !! Function : s11ks1622 !!-------------------------------------------------------------------1623 REAL(wp), INTENT(in ) :: x,y,z,phi1624 1625 REAL(wp) :: s11ks, p, pih1626 1627 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21, n1t2i221628 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21, n2t1i221629 REAL(wp) :: t1t2i11, t1t2i12, t1t2i21, t1t2i221630 REAL(wp) :: t2t1i11, t2t1i12, t2t1i21, t2t1i221631 REAL(wp) :: d11, d12, d221632 REAL(wp) :: IIn1t2, IIn2t1, IIt1t21633 REAL(wp) :: Hen1t2, Hen2t11634 !!-------------------------------------------------------------------1635 pih = 0.5_wp*rpi1636 p = phi1637 1638 n1t2i11 = cos(z+pih-p) * cos(z+p)1639 n1t2i12 = cos(z+pih-p) * sin(z+p)1640 n1t2i21 = sin(z+pih-p) * cos(z+p)1641 n1t2i22 = sin(z+pih-p) * sin(z+p)1642 n2t1i11 = cos(z-pih+p) * cos(z-p)1643 n2t1i12 = cos(z-pih+p) * sin(z-p)1644 n2t1i21 = sin(z-pih+p) * cos(z-p)1645 n2t1i22 = sin(z-pih+p) * sin(z-p)1646 t1t2i11 = cos(z-p) * cos(z+p)1647 t1t2i12 = cos(z-p) * sin(z+p)1648 t1t2i21 = sin(z-p) * cos(z+p)1649 t1t2i22 = sin(z-p) * sin(z+p)1650 t2t1i11 = cos(z+p) * cos(z-p)1651 t2t1i12 = cos(z+p) * sin(z-p)1652 t2t1i21 = sin(z+p) * cos(z-p)1653 t2t1i22 = sin(z+p) * sin(z-p)1654 d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))1655 d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))1656 d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))1657 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d221658 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d221659 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d221660 1661 IF (-IIn1t2>=rsmall) THEN1662 Hen1t2 = 1._wp1663 ELSE1664 Hen1t2 = 0._wp1665 ENDIF1666 1667 IF (-IIn2t1>=rsmall) THEN1668 Hen2t1 = 1._wp1669 ELSE1670 Hen2t1 = 0._wp1671 ENDIF1672 1673 s11ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i11 + Hen2t1 * t2t1i11)1674 1675 END FUNCTION s11ks1676 1677 FUNCTION s12ks(x,y,z,phi)1678 !!-------------------------------------------------------------------1679 !! Function : s12ks1680 !!-------------------------------------------------------------------1681 REAL(wp), INTENT(in ) :: x,y,z,phi1682 1683 REAL(wp) :: s12ks, s12s0, s21s0, p, pih1684 1685 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21, n1t2i221686 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21, n2t1i221687 REAL(wp) :: t1t2i11, t1t2i12, t1t2i21, t1t2i221688 REAL(wp) :: t2t1i11, t2t1i12, t2t1i21, t2t1i221689 REAL(wp) :: d11, d12, d221690 REAL(wp) :: IIn1t2, IIn2t1, IIt1t21691 REAL(wp) :: Hen1t2, Hen2t11692 !!-------------------------------------------------------------------1693 pih = 0.5_wp*rpi1694 p =phi1695 1696 n1t2i11 = cos(z+pih-p) * cos(z+p)1697 n1t2i12 = cos(z+pih-p) * sin(z+p)1698 n1t2i21 = sin(z+pih-p) * cos(z+p)1699 n1t2i22 = sin(z+pih-p) * sin(z+p)1700 n2t1i11 = cos(z-pih+p) * cos(z-p)1701 n2t1i12 = cos(z-pih+p) * sin(z-p)1702 n2t1i21 = sin(z-pih+p) * cos(z-p)1703 n2t1i22 = sin(z-pih+p) * sin(z-p)1704 t1t2i11 = cos(z-p) * cos(z+p)1705 t1t2i12 = cos(z-p) * sin(z+p)1706 t1t2i21 = sin(z-p) * cos(z+p)1707 t1t2i22 = sin(z-p) * sin(z+p)1708 t2t1i11 = cos(z+p) * cos(z-p)1709 t2t1i12 = cos(z+p) * sin(z-p)1710 t2t1i21 = sin(z+p) * cos(z-p)1711 t2t1i22 = sin(z+p) * sin(z-p)1712 d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))1713 d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))1714 d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))1715 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d221716 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d221717 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d221718 1719 IF (-IIn1t2>=rsmall) THEN1720 Hen1t2 = 1._wp1721 ELSE1722 Hen1t2 = 0._wp1723 ENDIF1724 1725 IF (-IIn2t1>=rsmall) THEN1726 Hen2t1 = 1._wp1727 ELSE1728 Hen2t1 = 0._wp1729 ENDIF1730 1731 s12s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i12 + Hen2t1 * t2t1i12)1732 s21s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i21 + Hen2t1 * t2t1i21)1733 s12ks=0.5_wp*(s12s0+s21s0)1734 1735 END FUNCTION s12ks1736 1737 FUNCTION s22ks(x,y,z,phi)1738 !!-------------------------------------------------------------------1739 !! Function : s22ks1740 !!-------------------------------------------------------------------1741 REAL(wp), INTENT(in ) :: x,y,z,phi1742 1743 REAL(wp) :: s22ks, p, pih1744 1745 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21, n1t2i221746 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21, n2t1i221747 REAL(wp) :: t1t2i11, t1t2i12, t1t2i21, t1t2i221748 REAL(wp) :: t2t1i11, t2t1i12, t2t1i21, t2t1i221749 REAL(wp) :: d11, d12, d221750 REAL(wp) :: IIn1t2, IIn2t1, IIt1t21751 REAL(wp) :: Hen1t2, Hen2t11752 !!-------------------------------------------------------------------1753 pih = 0.5_wp*rpi1754 p =phi1755 1756 n1t2i11 = cos(z+pih-p) * cos(z+p)1757 n1t2i12 = cos(z+pih-p) * sin(z+p)1758 n1t2i21 = sin(z+pih-p) * cos(z+p)1759 n1t2i22 = sin(z+pih-p) * sin(z+p)1760 n2t1i11 = cos(z-pih+p) * cos(z-p)1761 n2t1i12 = cos(z-pih+p) * sin(z-p)1762 n2t1i21 = sin(z-pih+p) * cos(z-p)1763 n2t1i22 = sin(z-pih+p) * sin(z-p)1764 t1t2i11 = cos(z-p) * cos(z+p)1765 t1t2i12 = cos(z-p) * sin(z+p)1766 t1t2i21 = sin(z-p) * cos(z+p)1767 t1t2i22 = sin(z-p) * sin(z+p)1768 t2t1i11 = cos(z+p) * cos(z-p)1769 t2t1i12 = cos(z+p) * sin(z-p)1770 t2t1i21 = sin(z+p) * cos(z-p)1771 t2t1i22 = sin(z+p) * sin(z-p)1772 d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))1773 d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))1774 d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))1775 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d221776 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d221777 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d221778 1779 IF (-IIn1t2>=rsmall) THEN1780 Hen1t2 = 1._wp1781 ELSE1782 Hen1t2 = 0._wp1783 ENDIF1784 1785 IF (-IIn2t1>=rsmall) THEN1786 Hen2t1 = 1._wp1787 ELSE1788 Hen2t1 = 0._wp1789 ENDIF1790 1784 1791 1785 s22ks = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i22 + Hen2t1 * t2t1i22)
Note: See TracChangeset
for help on using the changeset viewer.
