Changeset 13574 for NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90
- Timestamp:
- 2020-10-07T18:02:40+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2019/dev_r11842_SI3-10_EAP/tests/ICE_RHEO/MY_SRC/icedyn_rhg_eap.F90
r12263 r13574 49 49 PUBLIC rhg_eap_rst ! called by icedyn_rhg.F90 50 50 51 REAL(wp), PARAMETER :: p hi = 3.141592653589793_wp/12._wp ! diamond shaped floe smaller angle (default phi = 30 deg)51 REAL(wp), PARAMETER :: pphi = 3.141592653589793_wp/12._wp ! diamond shaped floe smaller angle (default phi = 30 deg) 52 52 53 53 ! look-up table for calculating structure tensor … … 67 67 CONTAINS 68 68 69 SUBROUTINE ice_dyn_rhg_eap( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, paniso_11, paniso_12 )69 SUBROUTINE ice_dyn_rhg_eap( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i, paniso_11, paniso_12, prdg_conv ) 70 70 !!------------------------------------------------------------------- 71 71 !! *** SUBROUTINE ice_dyn_rhg_eap *** … … 123 123 REAL(wp), DIMENSION(:,:), INTENT( out) :: pshear_i , pdivu_i , pdelta_i ! 124 124 REAL(wp), DIMENSION(:,:), INTENT(inout) :: paniso_11 , paniso_12 ! structure tensor components 125 REAL(wp), DIMENSION(:,:), INTENT(inout) :: prdg_conv ! for ridging 125 126 !! 126 INTEGER :: ji, jj , ij0, ij1, ii0, ii1! dummy loop indices127 INTEGER :: ji, jj ! dummy loop indices 127 128 INTEGER :: jter ! local integers 128 129 ! … … 141 142 REAL(wp) :: zfac_x, zfac_y 142 143 REAL(wp) :: zshear, zdum1, zdum2 143 REAL(wp) :: stressptmp, stressmtmp, stress12tmpF ! anisotropic stress tensor components 144 REAL(wp) :: alphar, alphas ! for mechanical redistribution 144 REAL(wp) :: zstressptmp, zstressmtmp, zstress12tmpF ! anisotropic stress tensor components 145 REAL(wp) :: zalphar, zalphas ! for mechanical redistribution 146 REAL(wp) :: zmresult11, zmresult12, z1dtevpkth, zp5kth, z1_dtevp_A ! for structure tensor evolution 145 147 REAL(wp) :: invw ! for test case 146 148 ! 147 REAL(wp), DIMENSION(jpi,jpj) :: stress12tmp ! anisotropic stress tensor component for regridding 149 REAL(wp), DIMENSION(jpi,jpj) :: zstress12tmp ! anisotropic stress tensor component for regridding 150 REAL(wp), DIMENSION(jpi,jpj) :: zyield11, zyield22, zyield12 ! yield surface tensor for history 148 151 REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points 149 152 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 … … 254 257 zs2 (:,:) = pstress2_i (:,:) 255 258 zs12(:,:) = pstress12_i(:,:) 259 260 ! constants for structure tensor 261 z1_dtevp_A = z1_dtevp/10.0_wp 262 z1dtevpkth = 1._wp / (z1_dtevp_A + 0.00002_wp) 263 zp5kth = 0.5_wp * 0.00002_wp 256 264 257 265 ! Ice strength … … 408 416 409 417 ! --- anisotropic stress calculation --- ! 410 CALL update_stress_rdg (jter, nn_nevp, phi,zdiv, zdt, &418 CALL update_stress_rdg (jter, nn_nevp, zdiv, zdt, & 411 419 zdsT, paniso_11(ji,jj), paniso_12(ji,jj), & 412 stressptmp, stressmtmp, & 413 stress12tmp(ji,jj), strength(ji,jj), & 414 alphar, alphas) 420 zstressptmp, zstressmtmp, & 421 zstress12tmp(ji,jj), strength(ji,jj), & 422 zalphar, zalphas) 423 424 ! structure tensor update 425 IF (mod(jter,10) == 0) THEN 426 CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), & 427 paniso_11(ji,jj), paniso_12(ji,jj), & 428 zmresult11, zmresult12) 429 430 paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A + zp5kth - zmresult11) * z1dtevpkth ! implicit 431 paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A - zmresult12) * z1dtevpkth ! implicit 432 ENDIF 433 434 415 435 IF (jter == nn_nevp) THEN 416 rdg_conv(ji,jj) = -min( alphar, 0._wp) 436 zyield11(ji,jj) = 0.5_wp * (zstressptmp + zstressmtmp) 437 zyield22(ji,jj) = 0.5_wp * (zstressptmp - zstressmtmp) 438 zyield12(ji,jj) = zstress12tmp(ji,jj) 439 prdg_conv(ji,jj) = -min( zalphar, 0._wp) 417 440 ENDIF 418 441 … … 426 449 ! gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 427 450 ! alpha = beta = sqrt(4*gamma) 428 IF( ln_aEVP ) THEN429 zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) )430 z1_alph1 = 1._wp / ( zalph1 + 1._wp )431 ENDIF451 !IF( ln_aEVP ) THEN 452 ! zalph1 = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 453 ! z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 454 !ENDIF 432 455 433 456 ! stress at T points (zkt/=0 if landfast) 434 457 !zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 435 458 !zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 436 zs1(ji,jj) = ( zs1(ji,jj) + stressptmp * zalph1 ) * z1_alph1437 zs2(ji,jj) = ( zs2(ji,jj) + stressmtmp * zalph1 ) * z1_alph1459 zs1(ji,jj) = ( zs1(ji,jj) + zstressptmp * zalph1 ) * z1_alph1 460 zs2(ji,jj) = ( zs2(ji,jj) + zstressmtmp * zalph1 ) * z1_alph1 438 461 !zs1(ji,jj) = ( stressptmp * zs1(ji,jj) + zalph1 ) * z1_alph1 439 462 !zs2(ji,jj) = ( stressmtmp * zs2(ji,jj) + zalph1 ) * z1_alph1 … … 442 465 END DO 443 466 !CALL lbc_lnk( 'icedyn_rhg_eap', zp_delt, 'T', 1. ) 444 CALL lbc_lnk ( 'icedyn_rhg_eap', stress12tmp, 'T', 1.)467 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zstress12tmp, 'T', 1. , paniso_11, 'T', 1. , paniso_12, 'T', 1.) 445 468 446 469 DO jj = 1, jpjm1 447 470 DO ji = 1, jpim1 448 471 ! stress12tmp at F points 449 stress12tmpF = ( stress12tmp(ji,jj+1) * e1e2t(ji,jj+1) +stress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1) &450 & + stress12tmp(ji,jj ) * e1e2t(ji,jj ) +stress12tmp(ji+1,jj ) * e1e2t(ji+1,jj ) &472 zstress12tmpF = ( zstress12tmp(ji,jj+1) * e1e2t(ji,jj+1) + zstress12tmp(ji+1,jj+1) * e1e2t(ji+1,jj+1) & 473 & + zstress12tmp(ji,jj ) * e1e2t(ji,jj ) + zstress12tmp(ji+1,jj ) * e1e2t(ji+1,jj ) & 451 474 & ) * 0.25_wp * r1_e1e2f(ji,jj) 452 475 … … 459 482 ! stress at F points (zkt/=0 if landfast) 460 483 !zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 * (1._wp + zkt) ) * 0.5_wp ) * z1_alph2 461 zs12(ji,jj) = ( zs12(ji,jj) + stress12tmpF * zalph1 ) * z1_alph1484 zs12(ji,jj) = ( zs12(ji,jj) + zstress12tmpF * zalph1 ) * z1_alph1 462 485 !zs12(ji,jj) = ( stress12tmpF * zs12(ji,jj) + zalph1 ) * z1_alph1 463 486 464 487 END DO 465 488 END DO 489 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 466 490 467 491 ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! … … 507 531 zvel = 5.e-05_wp + SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) 508 532 zTauB = ztauy_base(ji,jj) / zvel 509 ! !--- OceanBottom-to-Ice stress 533 ! !--- OceanBottom-to-Ice stress 510 534 ztauy_bi(ji,jj) = zTauB * v_ice(ji,jj) 511 535 ! … … 539 563 !extra code for test case boundary conditions 540 564 IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 541 v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) )565 v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 542 566 END IF 543 544 567 END DO 545 568 END DO … … 557 580 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 558 581 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 559 ! !--- Ocean-to-Ice stress 582 ! !--- Ocean-to-Ice stress 560 583 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 561 584 ! … … 595 618 !extra code for test case boundary conditions 596 619 IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 597 u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) )620 u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 598 621 END IF 599 600 622 END DO 601 623 END DO … … 615 637 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 616 638 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 617 ! !--- Ocean-to-Ice stress 639 ! !--- Ocean-to-Ice stress 618 640 ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 619 641 ! … … 653 675 !extra code for test case boundary conditions 654 676 IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 655 u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) )677 u_ice(ji,jj) = invw*(ztaux_ai(ji,jj) + zCorU(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj)) 656 678 END IF 657 658 679 END DO 659 680 END DO … … 709 730 !extra code for test case boundary conditions 710 731 IF (mjg(jj)<25 .or. mjg(jj)>975 .or. mig(ji)<25 .or. mig(ji)>975) THEN 711 v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) )732 v_ice(ji,jj) = invw*(ztauy_ai(ji,jj) + zCorV(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj)) 712 733 END IF 713 714 734 END DO 715 735 END DO … … 732 752 !!$ ENDIF 733 753 ! 754 734 755 ! ! ==================== ! 735 756 END DO ! end loop over jter ! 736 757 ! ! ==================== ! 758 ! 759 CALL lbc_lnk( 'icedyn_rhg_eap', prdg_conv, 'T', 1. ) ! only need this in ridging module after rheology completed 737 760 ! 738 761 !------------------------------------------------------------------------------! … … 838 861 !! ! (scheme converges if this value is ~1, see Bouillon et al 2009 (eq. 11)) 839 862 zsig1(ji,jj) = 0.5_wp * zdum2 * ( pstress1_i(ji,jj) ) ! compressive stress, see Bouillon et al. 2015 840 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 863 zsig2(ji,jj) = 0.5_wp * zdum2 * ( zshear ) ! shear stress 841 864 zsig3(ji,jj) = zdum2**2 * ( ( pstress1_i(ji,jj) + strength(ji,jj) )**2 + ( rn_ecc * zshear )**2 ) 842 865 END DO … … 850 873 ! Stress tensor invariants (normal and shear stress N/m) 851 874 IF( iom_use('normstr') ) CALL iom_put( 'normstr' , ( zs1(:,:) + zs2(:,:) ) * zmsk00(:,:) ) ! Normal stress 852 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 875 IF( iom_use('sheastr') ) CALL iom_put( 'sheastr' , SQRT( ( zs1(:,:) - zs2(:,:) )**2 + 4*zs12(:,:)**2 ) * zmsk00(:,:) ) ! Shear stress 853 876 854 877 DEALLOCATE( zsig1 , zsig2 , zsig3 ) 855 878 ENDIF 856 879 880 ! --- yieldcurve --- ! 881 IF( iom_use('yield11') .OR. iom_use('yield12') .OR. iom_use('yield22')) THEN 882 883 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zyield11, 'T', 1., zyield22, 'T', 1., zyield12, 'T', 1. ) 884 885 CALL iom_put( 'yield11', zyield11 * zmsk00 ) 886 CALL iom_put( 'yield22', zyield22 * zmsk00 ) 887 CALL iom_put( 'yield12', zyield12 * zmsk00 ) 888 ENDIF 889 857 890 ! --- anisotropy tensor --- ! 858 IF( iom_use('aniso') ) THEN 859 CALL iom_put( 'aniso' , aniso_11 * zmsk00 ) 891 IF( iom_use('aniso') ) THEN 892 CALL lbc_lnk( 'icedyn_rhg_eap', paniso_11, 'T', 1. ) 893 CALL iom_put( 'aniso' , paniso_11 * zmsk00 ) 860 894 ENDIF 861 895 … … 917 951 END SUBROUTINE ice_dyn_rhg_eap 918 952 919 SUBROUTINE update_stress_rdg (ksub, ndte, phi, divu,tension, &920 shear, a11,a12, &921 stressp,stressm, &922 stress12, strength, &923 alphar,alphas)953 SUBROUTINE update_stress_rdg (ksub, kndte, pdivu, ptension, & 954 pshear, pa11, pa12, & 955 pstressp, pstressm, & 956 pstress12, strength, & 957 palphar, palphas) 924 958 !!--------------------------------------------------------------------- 925 !! *** SUBROUTINE rhg_eap_rst *** 926 !! Updates the stress depending on values of strain rate and structure 927 !! tensor and for ksub=ndte it computes closing and sliding rate 959 !! *** SUBROUTINE update_stress_rdg *** 960 !! 961 !! ** Purpose : Updates the stress depending on values of strain rate and structure 962 !! tensor and for the last subcycle step it computes closing and sliding rate 928 963 !!--------------------------------------------------------------------- 929 INTEGER, INTENT(in ) :: ksub, ndte930 REAL(wp), INTENT(in ) :: phi,strength931 REAL(wp), INTENT(in ) :: divu, tension,shear932 REAL(wp), INTENT(in ) :: a11,a12933 REAL(wp), INTENT( out) :: stressp, stressm,stress12934 REAL(wp), INTENT( out) :: alphar,alphas964 INTEGER, INTENT(in ) :: ksub, kndte 965 REAL(wp), INTENT(in ) :: strength 966 REAL(wp), INTENT(in ) :: pdivu, ptension, pshear 967 REAL(wp), INTENT(in ) :: pa11, pa12 968 REAL(wp), INTENT( out) :: pstressp, pstressm, pstress12 969 REAL(wp), INTENT( out) :: palphar, palphas 935 970 936 971 INTEGER :: kx ,ky, ka 937 972 938 REAL(wp) :: stemp11r, stemp12r,stemp22r939 REAL(wp) :: stemp11s, stemp12s,stemp22s940 REAL(wp) :: a22, Qd11Qd11, Qd11Qd12,Qd12Qd12941 REAL(wp) :: Q11Q11, Q11Q12,Q12Q12942 REAL(wp) :: dtemp11, dtemp12,dtemp22943 REAL(wp) :: rotstemp11r, rotstemp12r,rotstemp22r944 REAL(wp) :: rotstemp11s, rotstemp12s,rotstemp22s945 REAL(wp) :: sig11, sig12,sig22946 REAL(wp) :: sgprm11, sgprm12,sgprm22947 REAL(wp) :: invstressconviso948 REAL(wp) :: Angle_denom_gamma,Angle_denom_alpha949 REAL(wp) :: Tany_1,Tany_2950 REAL(wp) :: x, y, dx, dy, da,kxw, kyw, kaw951 REAL(wp) :: invdx, invdy,invda952 REAL(wp) :: dtemp1, dtemp2, atempprime, a,invsin973 REAL(wp) :: zstemp11r, zstemp12r, zstemp22r 974 REAL(wp) :: zstemp11s, zstemp12s, zstemp22s 975 REAL(wp) :: za22, zQd11Qd11, zQd11Qd12, zQd12Qd12 976 REAL(wp) :: zQ11Q11, zQ11Q12, zQ12Q12 977 REAL(wp) :: zdtemp11, zdtemp12, zdtemp22 978 REAL(wp) :: zrotstemp11r, zrotstemp12r, zrotstemp22r 979 REAL(wp) :: zrotstemp11s, zrotstemp12s, zrotstemp22s 980 REAL(wp) :: zsig11, zsig12, zsig22 981 REAL(wp) :: zsgprm11, zsgprm12, zsgprm22 982 REAL(wp) :: zinvstressconviso 983 REAL(wp) :: zAngle_denom_gamma, zAngle_denom_alpha 984 REAL(wp) :: zTany_1, zTany_2 985 REAL(wp) :: zx, zy, zdx, zdy, zda, zkxw, kyw, kaw 986 REAL(wp) :: zinvdx, zinvdy, zinvda 987 REAL(wp) :: zdtemp1, zdtemp2, zatempprime, zinvsin 953 988 954 989 REAL(wp), PARAMETER :: kfriction = 0.45_wp … … 956 991 ! Factor to maintain the same stress as in EVP (see Section 3) 957 992 ! Can be set to 1 otherwise 958 invstressconviso = 1._wp/(1._wp+kfriction*kfriction)993 zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 959 994 960 invsin = 1._wp/sin(2._wp*phi) *invstressconviso995 zinvsin = 1._wp/sin(2._wp*pphi) * zinvstressconviso 961 996 !now uses phi as set in higher code 962 997 … … 965 1000 966 1001 ! 1) structure tensor 967 a22 = 1._wp-a11968 Q11Q11 = 1._wp969 Q12Q12 = rsmall970 Q11Q12 = rsmall1002 za22 = 1._wp-pa11 1003 zQ11Q11 = 1._wp 1004 zQ12Q12 = rsmall 1005 zQ11Q12 = rsmall 971 1006 972 1007 ! gamma: angle between general coordiantes and principal axis of A 973 1008 ! here Tan2gamma = 2 a12 / (a11 - a22) 974 if((ABS(a11 - a22) > rsmall).or.(ABS(a12) > rsmall)) then975 Angle_denom_gamma = 1._wp/sqrt( ( a11 - a22 )*( a11 -a22) + &976 4._wp* a12*a12 )1009 IF((ABS(pa11 - za22) > rsmall).OR.(ABS(pa12) > rsmall)) THEN 1010 zAngle_denom_gamma = 1._wp/sqrt( ( pa11 - za22 )*( pa11 - za22) + & 1011 4._wp*pa12*pa12 ) 977 1012 978 Q11Q11 = 0.5_wp + ( a11 - a22 )*0.5_wp*Angle_denom_gamma !Cos^2979 Q12Q12 = 0.5_wp - ( a11 - a22 )*0.5_wp*Angle_denom_gamma !Sin^2980 Q11Q12 = a12*Angle_denom_gamma !CosSin981 endif1013 zQ11Q11 = 0.5_wp + ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma !Cos^2 1014 zQ12Q12 = 0.5_wp - ( pa11 - za22 )*0.5_wp*zAngle_denom_gamma !Sin^2 1015 zQ11Q12 = pa12*zAngle_denom_gamma !CosSin 1016 ENDIF 982 1017 983 1018 ! rotation Q*atemp*Q^T 984 atempprime = Q11Q11*a11 + 2.0_wp*Q11Q12*a12 + Q12Q12*a221019 zatempprime = zQ11Q11*pa11 + 2.0_wp*zQ11Q12*pa12 + zQ12Q12*za22 985 1020 986 1021 ! make first principal value the largest 987 atempprime = max(atempprime, 1.0_wp -atempprime)1022 zatempprime = max(zatempprime, 1.0_wp - zatempprime) 988 1023 989 1024 ! 2) strain rate 990 dtemp11 = 0.5_wp*(divu + tension) 991 dtemp12 = shear*0.5_wp 992 dtemp22 = 0.5_wp*(divu - tension) 993 994 !WRITE(numout,*) 'div, tension, shear inside loop', ksub, divu, tension, shear 1025 zdtemp11 = 0.5_wp*(pdivu + ptension) 1026 zdtemp12 = pshear*0.5_wp 1027 zdtemp22 = 0.5_wp*(pdivu - ptension) 995 1028 996 1029 ! here Tan2alpha = 2 dtemp12 / (dtemp11 - dtemp22) 997 1030 998 Qd11Qd11 = 1.0_wp999 Qd12Qd12 = rsmall1000 Qd11Qd12 = rsmall1001 1002 if((ABS( dtemp11 - dtemp22) > rsmall).or. (ABS(dtemp12) > rsmall)) then1003 1004 Angle_denom_alpha = 1.0_wp/sqrt( ( dtemp11 -dtemp22 )* &1005 ( dtemp11 - dtemp22 ) + 4.0_wp*dtemp12*dtemp12)1006 1007 Qd11Qd11 = 0.5_wp + ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Cos^21008 Qd12Qd12 = 0.5_wp - ( dtemp11 - dtemp22 )*0.5_wp*Angle_denom_alpha !Sin^21009 Qd11Qd12 = dtemp12*Angle_denom_alpha !CosSin1010 endif1011 1012 dtemp1 = Qd11Qd11*dtemp11 + 2.0_wp*Qd11Qd12*dtemp12 + Qd12Qd12*dtemp221013 dtemp2 = Qd12Qd12*dtemp11 - 2.0_wp*Qd11Qd12*dtemp12 + Qd11Qd11*dtemp221031 zQd11Qd11 = 1.0_wp 1032 zQd12Qd12 = rsmall 1033 zQd11Qd12 = rsmall 1034 1035 IF((ABS( zdtemp11 - zdtemp22) > rsmall).OR. (ABS(zdtemp12) > rsmall)) THEN 1036 1037 zAngle_denom_alpha = 1.0_wp/sqrt( ( zdtemp11 - zdtemp22 )* & 1038 ( zdtemp11 - zdtemp22 ) + 4.0_wp*zdtemp12*zdtemp12) 1039 1040 zQd11Qd11 = 0.5_wp + ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Cos^2 1041 zQd12Qd12 = 0.5_wp - ( zdtemp11 - zdtemp22 )*0.5_wp*zAngle_denom_alpha !Sin^2 1042 zQd11Qd12 = zdtemp12*zAngle_denom_alpha !CosSin 1043 ENDIF 1044 1045 zdtemp1 = zQd11Qd11*zdtemp11 + 2.0_wp*zQd11Qd12*zdtemp12 + zQd12Qd12*zdtemp22 1046 zdtemp2 = zQd12Qd12*zdtemp11 - 2.0_wp*zQd11Qd12*zdtemp12 + zQd11Qd11*zdtemp22 1014 1047 ! In cos and sin values 1015 x = 0._wp1016 if ((ABS(dtemp1) > rsmall).or.(ABS(dtemp2) > rsmall)) then1017 x = atan2(dtemp2,dtemp1)1018 endif1048 zx = 0._wp 1049 IF ((ABS(zdtemp1) > rsmall).OR.(ABS(zdtemp2) > rsmall)) THEN 1050 zx = atan2(zdtemp2,zdtemp1) 1051 ENDIF 1019 1052 1020 1053 ! to ensure the angle lies between pi/4 and 9 pi/4 1021 if (x < rpi*0.25_wp) x =x + rpi*2.0_wp1054 IF (zx < rpi*0.25_wp) zx = zx + rpi*2.0_wp 1022 1055 1023 1056 ! y: angle between major principal axis of strain rate and structure … … 1027 1060 ! Tany = (Singamma*Cosgamma - Sinalpha*Cosgamma)/(Cos^2gamma - Sin^alpha) 1028 1061 1029 Tany_1 = Q11Q12 -Qd11Qd121030 Tany_2 = Q11Q11 -Qd12Qd121062 zTany_1 = zQ11Q12 - zQd11Qd12 1063 zTany_2 = zQ11Q11 - zQd12Qd12 1031 1064 1032 y = 0._wp1065 zy = 0._wp 1033 1066 1034 if ((ABS(Tany_1) > rsmall).or.(ABS(Tany_2) > rsmall)) then1035 y = atan2(Tany_1,Tany_2)1036 endif1067 IF ((ABS(zTany_1) > rsmall).OR.(ABS(zTany_2) > rsmall)) THEN 1068 zy = atan2(zTany_1,zTany_2) 1069 ENDIF 1037 1070 1038 1071 ! to make sure y is between 0 and pi 1039 if (y > rpi) y =y - rpi1040 if (y < 0) y =y + rpi1072 IF (zy > rpi) zy = zy - rpi 1073 IF (zy < 0) zy = zy + rpi 1041 1074 1042 1075 ! 3) update anisotropic stress tensor 1043 dx = rpi/real(nx_yield-1,kind=wp)1044 dy = rpi/real(ny_yield-1,kind=wp)1045 da = 0.5_wp/real(na_yield-1,kind=wp)1046 invdx = 1._wp/dx1047 invdy = 1._wp/dy1048 invda = 1._wp/da1076 zdx = rpi/real(nx_yield-1,kind=wp) 1077 zdy = rpi/real(ny_yield-1,kind=wp) 1078 zda = 0.5_wp/real(na_yield-1,kind=wp) 1079 zinvdx = 1._wp/zdx 1080 zinvdy = 1._wp/zdy 1081 zinvda = 1._wp/zda 1049 1082 1050 1083 ! % need 8 coords and 8 weights 1051 1084 ! % range in kx 1052 kx = int(( x-rpi*0.25_wp-rpi)*invdx) + 11053 kxw = kx - (x-rpi*0.25_wp-rpi)*invdx1085 kx = int((zx-rpi*0.25_wp-rpi)*zinvdx) + 1 1086 zkxw = kx - (zx-rpi*0.25_wp-rpi)*zinvdx 1054 1087 1055 ky = int( y*invdy) + 11056 kyw = ky - y*invdy1088 ky = int(zy*zinvdy) + 1 1089 kyw = ky - zy*zinvdy 1057 1090 1058 ka = int(( atempprime-0.5_wp)*invda) + 11059 kaw = ka - ( atempprime-0.5_wp)*invda1091 ka = int((zatempprime-0.5_wp)*zinvda) + 1 1092 kaw = ka - (zatempprime-0.5_wp)*zinvda 1060 1093 1061 1094 ! % Determine sigma_r(A1,Zeta,y) and sigma_s (see Section A1 of Tsamados 2013) 1062 stemp11r =kxw * kyw * kaw * s11r(kx ,ky ,ka ) &1063 & + (1._wp- kxw) * kyw * kaw * s11r(kx+1,ky ,ka ) &1064 & + kxw * (1._wp-kyw) * kaw * s11r(kx ,ky+1,ka ) &1065 & + kxw * kyw * (1._wp-kaw) * s11r(kx ,ky ,ka+1) &1066 & + (1._wp- kxw) * (1._wp-kyw) * kaw * s11r(kx+1,ky+1,ka ) &1067 & + (1._wp- kxw) * kyw * (1._wp-kaw) * s11r(kx+1,ky ,ka+1) &1068 & + kxw * (1._wp-kyw) * (1._wp-kaw) * s11r(kx ,ky+1,ka+1) &1069 & + (1._wp- kxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1)1070 stemp12r =kxw * kyw * kaw * s12r(kx ,ky ,ka ) &1071 & + (1._wp- kxw) * kyw * kaw * s12r(kx+1,ky ,ka ) &1072 & + kxw * (1._wp-kyw) * kaw * s12r(kx ,ky+1,ka ) &1073 & + kxw * kyw * (1._wp-kaw) * s12r(kx ,ky ,ka+1) &1074 & + (1._wp- kxw) * (1._wp-kyw) * kaw * s12r(kx+1,ky+1,ka ) &1075 & + (1._wp- kxw) * kyw * (1._wp-kaw) * s12r(kx+1,ky ,ka+1) &1076 & + kxw * (1._wp-kyw) * (1._wp-kaw) * s12r(kx ,ky+1,ka+1) &1077 & + (1._wp- kxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1)1078 stemp22r =kxw * kyw * kaw * s22r(kx ,ky ,ka ) &1079 & + (1._wp- kxw) * kyw * kaw * s22r(kx+1,ky ,ka ) &1080 & + kxw * (1._wp-kyw) * kaw * s22r(kx ,ky+1,ka ) &1081 & + kxw * kyw * (1._wp-kaw) * s22r(kx ,ky ,ka+1) &1082 & + (1._wp- kxw) * (1._wp-kyw) * kaw * s22r(kx+1,ky+1,ka ) &1083 & + (1._wp- kxw) * kyw * (1._wp-kaw) * s22r(kx+1,ky ,ka+1) &1084 & + kxw * (1._wp-kyw) * (1._wp-kaw) * s22r(kx ,ky+1,ka+1) &1085 & + (1._wp- kxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1)1095 zstemp11r = zkxw * kyw * kaw * s11r(kx ,ky ,ka ) & 1096 & + (1._wp-zkxw) * kyw * kaw * s11r(kx+1,ky ,ka ) & 1097 & + zkxw * (1._wp-kyw) * kaw * s11r(kx ,ky+1,ka ) & 1098 & + zkxw * kyw * (1._wp-kaw) * s11r(kx ,ky ,ka+1) & 1099 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s11r(kx+1,ky+1,ka ) & 1100 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s11r(kx+1,ky ,ka+1) & 1101 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s11r(kx ,ky+1,ka+1) & 1102 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11r(kx+1,ky+1,ka+1) 1103 zstemp12r = zkxw * kyw * kaw * s12r(kx ,ky ,ka ) & 1104 & + (1._wp-zkxw) * kyw * kaw * s12r(kx+1,ky ,ka ) & 1105 & + zkxw * (1._wp-kyw) * kaw * s12r(kx ,ky+1,ka ) & 1106 & + zkxw * kyw * (1._wp-kaw) * s12r(kx ,ky ,ka+1) & 1107 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s12r(kx+1,ky+1,ka ) & 1108 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s12r(kx+1,ky ,ka+1) & 1109 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s12r(kx ,ky+1,ka+1) & 1110 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12r(kx+1,ky+1,ka+1) 1111 zstemp22r = zkxw * kyw * kaw * s22r(kx ,ky ,ka ) & 1112 & + (1._wp-zkxw) * kyw * kaw * s22r(kx+1,ky ,ka ) & 1113 & + zkxw * (1._wp-kyw) * kaw * s22r(kx ,ky+1,ka ) & 1114 & + zkxw * kyw * (1._wp-kaw) * s22r(kx ,ky ,ka+1) & 1115 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s22r(kx+1,ky+1,ka ) & 1116 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s22r(kx+1,ky ,ka+1) & 1117 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s22r(kx ,ky+1,ka+1) & 1118 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22r(kx+1,ky+1,ka+1) 1086 1119 1087 stemp11s =kxw * kyw * kaw * s11s(kx ,ky ,ka ) &1088 & + (1._wp- kxw) * kyw * kaw * s11s(kx+1,ky ,ka ) &1089 & + kxw * (1._wp-kyw) * kaw * s11s(kx ,ky+1,ka ) &1090 & + kxw * kyw * (1._wp-kaw) * s11s(kx ,ky ,ka+1) &1091 & + (1._wp- kxw) * (1._wp-kyw) * kaw * s11s(kx+1,ky+1,ka ) &1092 & + (1._wp- kxw) * kyw * (1._wp-kaw) * s11s(kx+1,ky ,ka+1) &1093 & + kxw * (1._wp-kyw) * (1._wp-kaw) * s11s(kx ,ky+1,ka+1) &1094 & + (1._wp- kxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1)1095 stemp12s =kxw * kyw * kaw * s12s(kx ,ky ,ka ) &1096 & + (1._wp- kxw) * kyw * kaw * s12s(kx+1,ky ,ka ) &1097 & + kxw * (1._wp-kyw) * kaw * s12s(kx ,ky+1,ka ) &1098 & + kxw * kyw * (1._wp-kaw) * s12s(kx ,ky ,ka+1) &1099 & + (1._wp- kxw) * (1._wp-kyw) * kaw * s12s(kx+1,ky+1,ka ) &1100 & + (1._wp- kxw) * kyw * (1._wp-kaw) * s12s(kx+1,ky ,ka+1) &1101 & + kxw * (1._wp-kyw) * (1._wp-kaw) * s12s(kx ,ky+1,ka+1) &1102 & + (1._wp- kxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1)1103 stemp22s =kxw * kyw * kaw * s22s(kx ,ky ,ka ) &1104 & + (1._wp- kxw) * kyw * kaw * s22s(kx+1,ky ,ka ) &1105 & + kxw * (1._wp-kyw) * kaw * s22s(kx ,ky+1,ka ) &1106 & + kxw * kyw * (1._wp-kaw) * s22s(kx ,ky ,ka+1) &1107 & + (1._wp- kxw) * (1._wp-kyw) * kaw * s22s(kx+1,ky+1,ka ) &1108 & + (1._wp- kxw) * kyw * (1._wp-kaw) * s22s(kx+1,ky ,ka+1) &1109 & + kxw * (1._wp-kyw) * (1._wp-kaw) * s22s(kx ,ky+1,ka+1) &1110 & + (1._wp- kxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1)1120 zstemp11s = zkxw * kyw * kaw * s11s(kx ,ky ,ka ) & 1121 & + (1._wp-zkxw) * kyw * kaw * s11s(kx+1,ky ,ka ) & 1122 & + zkxw * (1._wp-kyw) * kaw * s11s(kx ,ky+1,ka ) & 1123 & + zkxw * kyw * (1._wp-kaw) * s11s(kx ,ky ,ka+1) & 1124 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s11s(kx+1,ky+1,ka ) & 1125 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s11s(kx+1,ky ,ka+1) & 1126 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s11s(kx ,ky+1,ka+1) & 1127 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s11s(kx+1,ky+1,ka+1) 1128 zstemp12s = zkxw * kyw * kaw * s12s(kx ,ky ,ka ) & 1129 & + (1._wp-zkxw) * kyw * kaw * s12s(kx+1,ky ,ka ) & 1130 & + zkxw * (1._wp-kyw) * kaw * s12s(kx ,ky+1,ka ) & 1131 & + zkxw * kyw * (1._wp-kaw) * s12s(kx ,ky ,ka+1) & 1132 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s12s(kx+1,ky+1,ka ) & 1133 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s12s(kx+1,ky ,ka+1) & 1134 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s12s(kx ,ky+1,ka+1) & 1135 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s12s(kx+1,ky+1,ka+1) 1136 zstemp22s = zkxw * kyw * kaw * s22s(kx ,ky ,ka ) & 1137 & + (1._wp-zkxw) * kyw * kaw * s22s(kx+1,ky ,ka ) & 1138 & + zkxw * (1._wp-kyw) * kaw * s22s(kx ,ky+1,ka ) & 1139 & + zkxw * kyw * (1._wp-kaw) * s22s(kx ,ky ,ka+1) & 1140 & + (1._wp-zkxw) * (1._wp-kyw) * kaw * s22s(kx+1,ky+1,ka ) & 1141 & + (1._wp-zkxw) * kyw * (1._wp-kaw) * s22s(kx+1,ky ,ka+1) & 1142 & + zkxw * (1._wp-kyw) * (1._wp-kaw) * s22s(kx ,ky+1,ka+1) & 1143 & + (1._wp-zkxw) * (1._wp-kyw) * (1._wp-kaw) * s22s(kx+1,ky+1,ka+1) 1111 1144 1112 1145 ! Calculate mean ice stress over a collection of floes (Equation 3 in 1113 1146 ! Tsamados 2013) 1114 1147 1115 sig11 = strength*(stemp11r + kfriction*stemp11s) * invsin 1116 sig12 = strength*(stemp12r + kfriction*stemp12s) * invsin 1117 sig22 = strength*(stemp22r + kfriction*stemp22s) * invsin 1118 1119 !WRITE(numout,*) 'strength inside loop', strength 1148 zsig11 = strength*(zstemp11r + kfriction*zstemp11s) * zinvsin 1149 zsig12 = strength*(zstemp12r + kfriction*zstemp12s) * zinvsin 1150 zsig22 = strength*(zstemp22r + kfriction*zstemp22s) * zinvsin 1151 1120 1152 !WRITE(numout,*) 'principal axis stress inside loop', sig11, sig12, sig22 1121 1153 … … 1123 1155 1124 1156 ! Update stress 1125 sgprm11 = Q11Q11*sig11 + Q12Q12*sig22 - 2._wp*Q11Q12 *sig12 1126 sgprm12 = Q11Q12*sig11 - Q11Q12*sig22 + (Q11Q11 - Q12Q12)*sig12 1127 sgprm22 = Q12Q12*sig11 + Q11Q11*sig22 + 2._wp*Q11Q12 *sig12 1128 1129 stressp = sgprm11 + sgprm22 1130 stress12 = sgprm12 1131 stressm = sgprm11 - sgprm22 1132 1133 !WRITE(numout,*) 'stress output inside loop', ksub, stressp 1157 zsgprm11 = zQ11Q11*zsig11 + zQ12Q12*zsig22 - 2._wp*zQ11Q12 *zsig12 1158 zsgprm12 = zQ11Q12*zsig11 - zQ11Q12*zsig22 + (zQ11Q11 - zQ12Q12)*zsig12 1159 zsgprm22 = zQ12Q12*zsig11 + zQ11Q11*zsig22 + 2._wp*zQ11Q12 *zsig12 1160 1161 pstressp = zsgprm11 + zsgprm22 1162 pstress12 = zsgprm12 1163 pstressm = zsgprm11 - zsgprm22 1134 1164 1135 1165 ! Compute ridging and sliding functions in general coordinates 1136 1166 ! (Equation 11 in Tsamados 2013) 1137 if (ksub == ndte) then1138 rotstemp11r = Q11Q11*stemp11r - 2._wp*Q11Q12*stemp12r &1139 + Q12Q12*stemp22r1140 rotstemp12r = Q11Q11*stemp12r + Q11Q12*(stemp11r-stemp22r) &1141 - Q12Q12*stemp12r1142 rotstemp22r = Q12Q12*stemp11r + 2._wp*Q11Q12*stemp12r &1143 + Q11Q11*stemp22r1144 1145 rotstemp11s = Q11Q11*stemp11s - 2._wp*Q11Q12*stemp12s &1146 + Q12Q12*stemp22s1147 rotstemp12s = Q11Q11*stemp12s + Q11Q12*(stemp11s-stemp22s) &1148 - Q12Q12*stemp12s1149 rotstemp22s = Q12Q12*stemp11s + 2._wp*Q11Q12*stemp12s &1150 + Q11Q11*stemp22s1151 1152 alphar = rotstemp11r*dtemp11 + 2._wp*rotstemp12r*dtemp12 &1153 + rotstemp22r*dtemp221154 alphas = rotstemp11s*dtemp11 + 2._wp*rotstemp12s*dtemp12 &1155 + rotstemp22s*dtemp221156 endif1167 IF (ksub == kndte) THEN 1168 zrotstemp11r = zQ11Q11*zstemp11r - 2._wp*zQ11Q12* zstemp12r & 1169 + zQ12Q12*zstemp22r 1170 zrotstemp12r = zQ11Q11*zstemp12r + zQ11Q12*(zstemp11r-zstemp22r) & 1171 - zQ12Q12*zstemp12r 1172 zrotstemp22r = zQ12Q12*zstemp11r + 2._wp*zQ11Q12* zstemp12r & 1173 + zQ11Q11*zstemp22r 1174 1175 zrotstemp11s = zQ11Q11*zstemp11s - 2._wp*zQ11Q12* zstemp12s & 1176 + zQ12Q12*zstemp22s 1177 zrotstemp12s = zQ11Q11*zstemp12s + zQ11Q12*(zstemp11s-zstemp22s) & 1178 - zQ12Q12*zstemp12s 1179 zrotstemp22s = zQ12Q12*zstemp11s + 2._wp*zQ11Q12* zstemp12s & 1180 + zQ11Q11*zstemp22s 1181 1182 palphar = zrotstemp11r*zdtemp11 + 2._wp*zrotstemp12r*zdtemp12 & 1183 + zrotstemp22r*zdtemp22 1184 palphas = zrotstemp11s*zdtemp11 + 2._wp*zrotstemp12s*zdtemp12 & 1185 + zrotstemp22s*zdtemp22 1186 ENDIF 1157 1187 END SUBROUTINE update_stress_rdg 1188 1189 !======================================================================= 1190 1191 1192 SUBROUTINE calc_ffrac (pstressp, pstressm, pstress12, & 1193 pa11, pa12, & 1194 pmresult11, pmresult12) 1195 !!--------------------------------------------------------------------- 1196 !! *** ROUTINE calc_ffrac *** 1197 !! 1198 !! ** Purpose : Computes term in evolution equation for structure tensor 1199 !! which determines the ice floe re-orientation due to fracture 1200 !! ** Method : Eq. 7: Ffrac = -kf(A-S) or = 0 depending on sigma_1 and sigma_2 1201 !!--------------------------------------------------------------------- 1202 REAL (wp), INTENT(in) :: pstressp, pstressm, pstress12, pa11, pa12 1203 REAL (wp), INTENT(out) :: pmresult11, pmresult12 1204 1205 ! local variables 1206 1207 REAL (wp) :: zsigma11, zsigma12, zsigma22 ! stress tensor elements 1208 REAL (wp) :: zAngle_denom ! angle with principal component axis 1209 REAL (wp) :: zsigma_1, zsigma_2 ! principal components of stress 1210 REAL (wp) :: zQ11, zQ12, zQ11Q11, zQ11Q12, zQ12Q12 1211 1212 REAL (wp), PARAMETER :: kfrac = 0.0001_wp ! rate of fracture formation 1213 REAL (wp), PARAMETER :: threshold = 0.3_wp ! critical confinement ratio 1214 !!--------------------------------------------------------------- 1215 ! 1216 zsigma11 = 0.5_wp*(pstressp+pstressm) 1217 zsigma12 = pstress12 1218 zsigma22 = 0.5_wp*(pstressp-pstressm) 1219 1220 ! Here's the change - no longer calculate gamma, 1221 ! use trig formulation, angle signs are kept correct, don't worry 1222 1223 ! rotate tensor to get into sigma principal axis 1224 1225 ! here Tan2gamma = 2 sig12 / (sig11 - sig12) 1226 ! add rsmall to the denominator to stop 1/0 errors, makes very little 1227 ! error to the calculated angles < rsmall 1228 1229 zQ11Q11 = 0.1_wp 1230 zQ12Q12 = rsmall 1231 zQ11Q12 = rsmall 1232 1233 IF((ABS( zsigma11 - zsigma22) > rsmall).OR.(ABS(zsigma12) > rsmall)) THEN 1234 1235 zAngle_denom = 1.0_wp/sqrt( ( zsigma11 - zsigma22 )*( zsigma11 - & 1236 zsigma22 ) + 4.0_wp*zsigma12*zsigma12) 1237 1238 zQ11Q11 = 0.5_wp + ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom !Cos^2 1239 zQ12Q12 = 0.5_wp - ( zsigma11 - zsigma22 )*0.5_wp*zAngle_denom !Sin^2 1240 zQ11Q12 = zsigma12*zAngle_denom !CosSin 1241 ENDIF 1242 1243 zsigma_1 = zQ11Q11*zsigma11 + 2.0_wp*zQ11Q12*zsigma12 & 1244 + zQ12Q12*zsigma22 ! S(1,1) 1245 zsigma_2 = zQ12Q12*zsigma11 - 2.0_wp*zQ11Q12*zsigma12 & 1246 + zQ11Q11*zsigma22 ! S(2,2) 1247 1248 ! Pure divergence 1249 IF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 >= 0.0_wp)) THEN 1250 pmresult11 = 0.0_wp 1251 pmresult12 = 0.0_wp 1252 1253 ! Unconfined compression: cracking of blocks not along the axial splitting 1254 ! direction 1255 ! which leads to the loss of their shape, so we again model it through diffusion 1256 ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp)) THEN 1257 pmresult11 = kfrac * (pa11 - zQ12Q12) 1258 pmresult12 = kfrac * (pa12 + zQ11Q12) 1259 1260 ! Shear faulting 1261 ELSEIF (zsigma_2 == 0.0_wp) THEN 1262 pmresult11 = 0.0_wp 1263 pmresult12 = 0.0_wp 1264 ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN 1265 pmresult11 = kfrac * (pa11 - zQ12Q12) 1266 pmresult12 = kfrac * (pa12 + zQ11Q12) 1267 1268 ! Horizontal spalling 1269 ELSE 1270 pmresult11 = 0.0_wp 1271 pmresult12 = 0.0_wp 1272 ENDIF 1273 1274 END SUBROUTINE calc_ffrac 1275 1158 1276 1159 1277 SUBROUTINE rhg_eap_rst( cdrw, kt ) … … 1234 1352 s12s(ix,iy,ia) = 0._wp 1235 1353 s22s(ix,iy,ia) = 0._wp 1236 IF (ia <= na_yield-1) then1354 IF (ia <= na_yield-1) THEN 1237 1355 DO iz=1,nz 1238 1356 s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1239 1357 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1240 s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz ,phi)*dz/sin(2._wp*phi)1358 s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1241 1359 s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1242 1360 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1243 s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz ,phi)*dz/sin(2._wp*phi)1361 s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1244 1362 s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1245 1363 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1246 s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz ,phi)*dz/sin(2._wp*phi)1364 s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1247 1365 s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1248 1366 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1249 s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz ,phi)*dz/sin(2._wp*phi)1367 s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1250 1368 s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1251 1369 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1252 s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz ,phi)*dz/sin(2._wp*phi)1370 s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1253 1371 s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1254 1372 exp(-w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1255 s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz ,phi)*dz/sin(2._wp*phi)1373 s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1256 1374 ENDDO 1257 1375 IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp … … 1262 1380 IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 1263 1381 ELSE 1264 s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp ,phi)/sin(2._wp*phi)1265 s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp ,phi)/sin(2._wp*phi)1266 s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp ,phi)/sin(2._wp*phi)1267 s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp ,phi)/sin(2._wp*phi)1268 s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp ,phi)/sin(2._wp*phi)1269 s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp ,phi)/sin(2._wp*phi)1382 s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1383 s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1384 s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1385 s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1386 s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1387 s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1270 1388 IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 1271 1389 IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp … … 1295 1413 1296 1414 1297 FUNCTION w1( a)1415 FUNCTION w1(pa) 1298 1416 !!------------------------------------------------------------------- 1299 1417 !! Function : w1 (see Gaussian function psi in Tsamados et al 2013) 1300 1418 !!------------------------------------------------------------------- 1301 REAL(wp), INTENT(in ) :: a1419 REAL(wp), INTENT(in ) :: pa 1302 1420 REAL(wp) :: w1 1303 1421 !!------------------------------------------------------------------- 1304 1422 1305 1423 w1 = - 223.87569446_wp & 1306 + 2361.2198663_wp*a &1307 - 10606.56079975_wp*a*a &1308 + 26315.50025642_wp*a*a*a &1309 - 38948.30444297_wp*a*a*a*a &1310 + 34397.72407466_wp*a*a*a*a*a &1311 - 16789.98003081_wp*a*a*a*a*a*a &1312 + 3495.82839237_wp*a*a*a*a*a*a*a1424 & + 2361.2198663_wp*pa & 1425 & - 10606.56079975_wp*pa*pa & 1426 & + 26315.50025642_wp*pa*pa*pa & 1427 & - 38948.30444297_wp*pa*pa*pa*pa & 1428 & + 34397.72407466_wp*pa*pa*pa*pa*pa & 1429 & - 16789.98003081_wp*pa*pa*pa*pa*pa*pa & 1430 & + 3495.82839237_wp*pa*pa*pa*pa*pa*pa*pa 1313 1431 1314 1432 END FUNCTION w1 1315 1433 1316 1434 1317 FUNCTION w2( a)1435 FUNCTION w2(pa) 1318 1436 !!------------------------------------------------------------------- 1319 1437 !! Function : w2 (see Gaussian function psi in Tsamados et al 2013) 1320 1438 !!------------------------------------------------------------------- 1321 REAL(wp), INTENT(in ) :: a1439 REAL(wp), INTENT(in ) :: pa 1322 1440 REAL(wp) :: w2 1323 1441 !!------------------------------------------------------------------- 1324 1442 1325 1443 w2 = - 6670.68911883_wp & 1326 + 70222.33061536_wp*a &1327 - 314871.71525448_wp*a*a &1328 + 779570.02793492_wp*a*a*a &1329 - 1151098.82436864_wp*a*a*a*a &1330 + 1013896.59464498_wp*a*a*a*a*a &1331 - 493379.44906738_wp*a*a*a*a*a*a &1332 + 102356.551518_wp*a*a*a*a*a*a*a1444 & + 70222.33061536_wp*pa & 1445 & - 314871.71525448_wp*pa*pa & 1446 & + 779570.02793492_wp*pa*pa*pa & 1447 & - 1151098.82436864_wp*pa*pa*pa*pa & 1448 & + 1013896.59464498_wp*pa*pa*pa*pa*pa & 1449 & - 493379.44906738_wp*pa*pa*pa*pa*pa*pa & 1450 & + 102356.551518_wp*pa*pa*pa*pa*pa*pa*pa 1333 1451 1334 1452 END FUNCTION w2 1335 1453 1336 FUNCTION s11kr( x,y,z,phi)1454 FUNCTION s11kr(px,py,pz) 1337 1455 !!------------------------------------------------------------------- 1338 1456 !! Function : s11kr 1339 1457 !!------------------------------------------------------------------- 1340 REAL(wp), INTENT(in ) :: x,y,z,phi1458 REAL(wp), INTENT(in ) :: px,py,pz 1341 1459 1342 REAL(wp) :: s11kr, p,pih1460 REAL(wp) :: s11kr, zpih 1343 1461 1344 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21,n1t2i221345 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21,n2t1i221346 REAL(wp) :: t1t2i11, t1t2i12, t1t2i21,t1t2i221347 REAL(wp) :: t2t1i11, t2t1i12, t2t1i21,t2t1i221348 REAL(wp) :: d11, d12,d221349 REAL(wp) :: IIn1t2, IIn2t1,IIt1t21350 REAL(wp) :: Hen1t2,Hen2t11462 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1463 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1464 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1465 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1466 REAL(wp) :: zd11, zd12, zd22 1467 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1468 REAL(wp) :: zHen1t2, zHen2t1 1351 1469 !!------------------------------------------------------------------- 1352 1470 1353 pih = 0.5_wp*rpi 1354 p = phi 1471 zpih = 0.5_wp*rpi 1355 1472 1356 n1t2i11 = cos(z+pih-p) * cos(z+p)1357 n1t2i12 = cos(z+pih-p) * sin(z+p)1358 n1t2i21 = sin(z+pih-p) * cos(z+p)1359 n1t2i22 = sin(z+pih-p) * sin(z+p)1360 n2t1i11 = cos(z-pih+p) * cos(z-p)1361 n2t1i12 = cos(z-pih+p) * sin(z-p)1362 n2t1i21 = sin(z-pih+p) * cos(z-p)1363 n2t1i22 = sin(z-pih+p) * sin(z-p)1364 t1t2i11 = cos(z-p) * cos(z+p)1365 t1t2i12 = cos(z-p) * sin(z+p)1366 t1t2i21 = sin(z-p) * cos(z+p)1367 t1t2i22 = sin(z-p) * sin(z+p)1368 t2t1i11 = cos(z+p) * cos(z-p)1369 t2t1i12 = cos(z+p) * sin(z-p)1370 t2t1i21 = sin(z+p) * cos(z-p)1371 t2t1i22 = sin(z+p) * sin(z-p)1473 zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 1474 zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 1475 zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 1476 zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 1477 zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 1478 zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 1479 zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 1480 zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 1481 zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 1482 zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 1483 zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 1484 zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 1485 zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 1486 zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 1487 zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 1488 zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 1372 1489 ! In expression of tensor d, with this formulatin d(x)=-d(x+pi) 1373 1490 ! Solution, when diagonalizing always check sgn(a11-a22) if > then keep x else 1374 1491 ! x=x-pi/2 1375 d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))1376 d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))1377 d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))1378 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 *d221379 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 *d221380 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 *d221492 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1493 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1494 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1495 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1496 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1497 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1381 1498 1382 IF (- IIn1t2>=rsmall) then1383 Hen1t2 = 1._wp1499 IF (-zIIn1t2>=rsmall) THEN 1500 zHen1t2 = 1._wp 1384 1501 ELSE 1385 Hen1t2 = 0._wp1502 zHen1t2 = 0._wp 1386 1503 ENDIF 1387 1504 1388 IF (- IIn2t1>=rsmall) then1389 Hen2t1 = 1._wp1505 IF (-zIIn2t1>=rsmall) THEN 1506 zHen2t1 = 1._wp 1390 1507 ELSE 1391 Hen2t1 = 0._wp1508 zHen2t1 = 0._wp 1392 1509 ENDIF 1393 1510 1394 s11kr = (- Hen1t2 * n1t2i11 - Hen2t1 *n2t1i11)1511 s11kr = (- zHen1t2 * zn1t2i11 - zHen2t1 * zn2t1i11) 1395 1512 1396 1513 END FUNCTION s11kr 1397 1514 1398 FUNCTION s12kr( x,y,z,phi)1515 FUNCTION s12kr(px,py,pz) 1399 1516 !!------------------------------------------------------------------- 1400 1517 !! Function : s12kr 1401 1518 !!------------------------------------------------------------------- 1402 REAL(wp), INTENT(in ) :: x,y,z,phi 1403 1404 REAL(wp) :: s12kr, s12r0, s21r0, p, pih 1405 1406 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21, n1t2i22 1407 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21, n2t1i22 1408 REAL(wp) :: t1t2i11, t1t2i12, t1t2i21, t1t2i22 1409 REAL(wp) :: t2t1i11, t2t1i12, t2t1i21, t2t1i22 1410 REAL(wp) :: d11, d12, d22 1411 REAL(wp) :: IIn1t2, IIn2t1, IIt1t2 1412 REAL(wp) :: Hen1t2, Hen2t1 1413 !!------------------------------------------------------------------- 1414 pih = 0.5_wp*rpi 1415 p = phi 1519 REAL(wp), INTENT(in ) :: px,py,pz 1520 1521 REAL(wp) :: s12kr, zs12r0, zs21r0, zpih 1522 1523 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1524 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1525 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1526 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1527 REAL(wp) :: zd11, zd12, zd22 1528 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1529 REAL(wp) :: zHen1t2, zHen2t1 1530 !!------------------------------------------------------------------- 1531 zpih = 0.5_wp*rpi 1416 1532 1417 n1t2i11 = cos(z+pih-p) * cos(z+p)1418 n1t2i12 = cos(z+pih-p) * sin(z+p)1419 n1t2i21 = sin(z+pih-p) * cos(z+p)1420 n1t2i22 = sin(z+pih-p) * sin(z+p)1421 n2t1i11 = cos(z-pih+p) * cos(z-p)1422 n2t1i12 = cos(z-pih+p) * sin(z-p)1423 n2t1i21 = sin(z-pih+p) * cos(z-p)1424 n2t1i22 = sin(z-pih+p) * sin(z-p)1425 t1t2i11 = cos(z-p) * cos(z+p)1426 t1t2i12 = cos(z-p) * sin(z+p)1427 t1t2i21 = sin(z-p) * cos(z+p)1428 t1t2i22 = sin(z-p) * sin(z+p)1429 t2t1i11 = cos(z+p) * cos(z-p)1430 t2t1i12 = cos(z+p) * sin(z-p)1431 t2t1i21 = sin(z+p) * cos(z-p)1432 t2t1i22 = sin(z+p) * sin(z-p)1433 d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y))1434 d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x))1435 d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y))1436 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 *d221437 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 *d221438 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 *d221533 zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 1534 zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 1535 zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 1536 zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 1537 zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 1538 zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 1539 zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 1540 zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 1541 zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 1542 zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 1543 zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 1544 zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 1545 zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 1546 zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 1547 zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 1548 zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 1549 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1550 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1551 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1552 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1553 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1554 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1439 1555 1440 IF (- IIn1t2>=rsmall) then1441 Hen1t2 = 1._wp1556 IF (-zIIn1t2>=rsmall) THEN 1557 zHen1t2 = 1._wp 1442 1558 ELSE 1443 Hen1t2 = 0._wp1559 zHen1t2 = 0._wp 1444 1560 ENDIF 1445 1561 1446 IF (- IIn2t1>=rsmall) then1447 Hen2t1 = 1._wp1562 IF (-zIIn2t1>=rsmall) THEN 1563 zHen2t1 = 1._wp 1448 1564 ELSE 1449 Hen2t1 = 0._wp1565 zHen2t1 = 0._wp 1450 1566 ENDIF 1451 1567 1452 s12r0 = (- Hen1t2 * n1t2i12 - Hen2t1 *n2t1i12)1453 s21r0 = (- Hen1t2 * n1t2i21 - Hen2t1 *n2t1i21)1454 s12kr=0.5_wp*( s12r0+s21r0)1568 zs12r0 = (- zHen1t2 * zn1t2i12 - zHen2t1 * zn2t1i12) 1569 zs21r0 = (- zHen1t2 * zn1t2i21 - zHen2t1 * zn2t1i21) 1570 s12kr=0.5_wp*(zs12r0+zs21r0) 1455 1571 1456 1572 END FUNCTION s12kr 1457 1573 1458 FUNCTION s22kr( x,y,z,phi)1574 FUNCTION s22kr(px,py,pz) 1459 1575 !!------------------------------------------------------------------- 1460 1576 !! Function : s22kr 1461 1577 !!------------------------------------------------------------------- 1462 REAL(wp), INTENT(in ) :: x,y,z,phi 1463 1464 REAL(wp) :: s22kr, p, pih 1465 1466 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21, n1t2i22 1467 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21, n2t1i22 1468 REAL(wp) :: t1t2i11, t1t2i12, t1t2i21, t1t2i22 1469 REAL(wp) :: t2t1i11, t2t1i12, t2t1i21, t2t1i22 1470 REAL(wp) :: d11, d12, d22 1471 REAL(wp) :: IIn1t2, IIn2t1, IIt1t2 1472 REAL(wp) :: Hen1t2, Hen2t1 1473 !!------------------------------------------------------------------- 1474 pih = 0.5_wp*rpi 1475 p = phi 1476 1477 n1t2i11 = cos(z+pih-p) * cos(z+p) 1478 n1t2i12 = cos(z+pih-p) * sin(z+p) 1479 n1t2i21 = sin(z+pih-p) * cos(z+p) 1480 n1t2i22 = sin(z+pih-p) * sin(z+p) 1481 n2t1i11 = cos(z-pih+p) * cos(z-p) 1482 n2t1i12 = cos(z-pih+p) * sin(z-p) 1483 n2t1i21 = sin(z-pih+p) * cos(z-p) 1484 n2t1i22 = sin(z-pih+p) * sin(z-p) 1485 t1t2i11 = cos(z-p) * cos(z+p) 1486 t1t2i12 = cos(z-p) * sin(z+p) 1487 t1t2i21 = sin(z-p) * cos(z+p) 1488 t1t2i22 = sin(z-p) * sin(z+p) 1489 t2t1i11 = cos(z+p) * cos(z-p) 1490 t2t1i12 = cos(z+p) * sin(z-p) 1491 t2t1i21 = sin(z+p) * cos(z-p) 1492 t2t1i22 = sin(z+p) * sin(z-p) 1493 d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 1494 d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 1495 d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 1496 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 1497 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 1498 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 1499 1500 IF (-IIn1t2>=rsmall) then 1501 Hen1t2 = 1._wp 1578 REAL(wp), INTENT(in ) :: px,py,pz 1579 1580 REAL(wp) :: s22kr, zpih 1581 1582 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1583 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1584 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1585 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1586 REAL(wp) :: zd11, zd12, zd22 1587 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1588 REAL(wp) :: zHen1t2, zHen2t1 1589 !!------------------------------------------------------------------- 1590 zpih = 0.5_wp*rpi 1591 1592 zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 1593 zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 1594 zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 1595 zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 1596 zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 1597 zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 1598 zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 1599 zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 1600 zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 1601 zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 1602 zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 1603 zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 1604 zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 1605 zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 1606 zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 1607 zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 1608 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1609 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1610 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1611 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1612 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1613 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1614 1615 IF (-zIIn1t2>=rsmall) THEN 1616 zHen1t2 = 1._wp 1502 1617 ELSE 1503 Hen1t2 = 0._wp1504 ENDIF 1505 1506 IF (- IIn2t1>=rsmall) then1507 Hen2t1 = 1._wp1618 zHen1t2 = 0._wp 1619 ENDIF 1620 1621 IF (-zIIn2t1>=rsmall) THEN 1622 zHen2t1 = 1._wp 1508 1623 ELSE 1509 Hen2t1 = 0._wp1510 ENDIF 1511 1512 s22kr = (- Hen1t2 * n1t2i22 - Hen2t1 *n2t1i22)1624 zHen2t1 = 0._wp 1625 ENDIF 1626 1627 s22kr = (- zHen1t2 * zn1t2i22 - zHen2t1 * zn2t1i22) 1513 1628 1514 1629 END FUNCTION s22kr 1515 1630 1516 FUNCTION s11ks( x,y,z,phi)1631 FUNCTION s11ks(px,py,pz) 1517 1632 !!------------------------------------------------------------------- 1518 1633 !! Function : s11ks 1519 1634 !!------------------------------------------------------------------- 1520 REAL(wp), INTENT(in ) :: x,y,z,phi 1521 1522 REAL(wp) :: s11ks, p, pih 1523 1524 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21, n1t2i22 1525 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21, n2t1i22 1526 REAL(wp) :: t1t2i11, t1t2i12, t1t2i21, t1t2i22 1527 REAL(wp) :: t2t1i11, t2t1i12, t2t1i21, t2t1i22 1528 REAL(wp) :: d11, d12, d22 1529 REAL(wp) :: IIn1t2, IIn2t1, IIt1t2 1530 REAL(wp) :: Hen1t2, Hen2t1 1531 !!------------------------------------------------------------------- 1532 pih = 0.5_wp*rpi 1533 p = phi 1534 1535 n1t2i11 = cos(z+pih-p) * cos(z+p) 1536 n1t2i12 = cos(z+pih-p) * sin(z+p) 1537 n1t2i21 = sin(z+pih-p) * cos(z+p) 1538 n1t2i22 = sin(z+pih-p) * sin(z+p) 1539 n2t1i11 = cos(z-pih+p) * cos(z-p) 1540 n2t1i12 = cos(z-pih+p) * sin(z-p) 1541 n2t1i21 = sin(z-pih+p) * cos(z-p) 1542 n2t1i22 = sin(z-pih+p) * sin(z-p) 1543 t1t2i11 = cos(z-p) * cos(z+p) 1544 t1t2i12 = cos(z-p) * sin(z+p) 1545 t1t2i21 = sin(z-p) * cos(z+p) 1546 t1t2i22 = sin(z-p) * sin(z+p) 1547 t2t1i11 = cos(z+p) * cos(z-p) 1548 t2t1i12 = cos(z+p) * sin(z-p) 1549 t2t1i21 = sin(z+p) * cos(z-p) 1550 t2t1i22 = sin(z+p) * sin(z-p) 1551 d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 1552 d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 1553 d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 1554 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 1555 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 1556 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 1557 1558 IF (-IIn1t2>=rsmall) then 1559 Hen1t2 = 1._wp 1635 REAL(wp), INTENT(in ) :: px,py,pz 1636 1637 REAL(wp) :: s11ks, zpih 1638 1639 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1640 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1641 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1642 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1643 REAL(wp) :: zd11, zd12, zd22 1644 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1645 REAL(wp) :: zHen1t2, zHen2t1 1646 !!------------------------------------------------------------------- 1647 zpih = 0.5_wp*rpi 1648 1649 zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 1650 zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 1651 zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 1652 zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 1653 zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 1654 zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 1655 zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 1656 zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 1657 zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 1658 zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 1659 zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 1660 zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 1661 zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 1662 zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 1663 zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 1664 zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 1665 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1666 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1667 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1668 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1669 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1670 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1671 1672 IF (-zIIn1t2>=rsmall) THEN 1673 zHen1t2 = 1._wp 1560 1674 ELSE 1561 Hen1t2 = 0._wp1562 ENDIF 1563 1564 IF (- IIn2t1>=rsmall) then1565 Hen2t1 = 1._wp1675 zHen1t2 = 0._wp 1676 ENDIF 1677 1678 IF (-zIIn2t1>=rsmall) THEN 1679 zHen2t1 = 1._wp 1566 1680 ELSE 1567 Hen2t1 = 0._wp1568 ENDIF 1569 1570 s11ks = sign(1._wp, IIt1t2+rsmall)*(Hen1t2 * t1t2i11 + Hen2t1 *t2t1i11)1681 zHen2t1 = 0._wp 1682 ENDIF 1683 1684 s11ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i11 + zHen2t1 * zt2t1i11) 1571 1685 1572 1686 END FUNCTION s11ks 1573 1687 1574 FUNCTION s12ks( x,y,z,phi)1688 FUNCTION s12ks(px,py,pz) 1575 1689 !!------------------------------------------------------------------- 1576 1690 !! Function : s12ks 1577 1691 !!------------------------------------------------------------------- 1578 REAL(wp), INTENT(in ) :: x,y,z,phi 1579 1580 REAL(wp) :: s12ks, s12s0, s21s0, p, pih 1581 1582 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21, n1t2i22 1583 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21, n2t1i22 1584 REAL(wp) :: t1t2i11, t1t2i12, t1t2i21, t1t2i22 1585 REAL(wp) :: t2t1i11, t2t1i12, t2t1i21, t2t1i22 1586 REAL(wp) :: d11, d12, d22 1587 REAL(wp) :: IIn1t2, IIn2t1, IIt1t2 1588 REAL(wp) :: Hen1t2, Hen2t1 1589 !!------------------------------------------------------------------- 1590 pih = 0.5_wp*rpi 1591 p =phi 1592 1593 n1t2i11 = cos(z+pih-p) * cos(z+p) 1594 n1t2i12 = cos(z+pih-p) * sin(z+p) 1595 n1t2i21 = sin(z+pih-p) * cos(z+p) 1596 n1t2i22 = sin(z+pih-p) * sin(z+p) 1597 n2t1i11 = cos(z-pih+p) * cos(z-p) 1598 n2t1i12 = cos(z-pih+p) * sin(z-p) 1599 n2t1i21 = sin(z-pih+p) * cos(z-p) 1600 n2t1i22 = sin(z-pih+p) * sin(z-p) 1601 t1t2i11 = cos(z-p) * cos(z+p) 1602 t1t2i12 = cos(z-p) * sin(z+p) 1603 t1t2i21 = sin(z-p) * cos(z+p) 1604 t1t2i22 = sin(z-p) * sin(z+p) 1605 t2t1i11 = cos(z+p) * cos(z-p) 1606 t2t1i12 = cos(z+p) * sin(z-p) 1607 t2t1i21 = sin(z+p) * cos(z-p) 1608 t2t1i22 = sin(z+p) * sin(z-p) 1609 d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 1610 d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 1611 d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 1612 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 1613 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 1614 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 1615 1616 IF (-IIn1t2>=rsmall) then 1617 Hen1t2 = 1._wp 1692 REAL(wp), INTENT(in ) :: px,py,pz 1693 1694 REAL(wp) :: s12ks, zs12s0, zs21s0, zpih 1695 1696 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1697 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1698 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1699 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1700 REAL(wp) :: zd11, zd12, zd22 1701 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1702 REAL(wp) :: zHen1t2, zHen2t1 1703 !!------------------------------------------------------------------- 1704 zpih = 0.5_wp*rpi 1705 1706 zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 1707 zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 1708 zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 1709 zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 1710 zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 1711 zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 1712 zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 1713 zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 1714 zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 1715 zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 1716 zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 1717 zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 1718 zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 1719 zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 1720 zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 1721 zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 1722 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1723 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1724 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1725 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1726 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1727 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1728 1729 IF (-zIIn1t2>=rsmall) THEN 1730 zHen1t2 = 1._wp 1618 1731 ELSE 1619 Hen1t2 = 0._wp1620 ENDIF 1621 1622 IF (- IIn2t1>=rsmall) then1623 Hen2t1 = 1._wp1732 zHen1t2 = 0._wp 1733 ENDIF 1734 1735 IF (-zIIn2t1>=rsmall) THEN 1736 zHen2t1 = 1._wp 1624 1737 ELSE 1625 Hen2t1 = 0._wp1626 ENDIF 1627 1628 s12s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i12 + Hen2t1 *t2t1i12)1629 s21s0 = sign(1._wp,IIt1t2+rsmall)*(Hen1t2 * t1t2i21 + Hen2t1 *t2t1i21)1630 s12ks=0.5_wp*( s12s0+s21s0)1738 zHen2t1 = 0._wp 1739 ENDIF 1740 1741 zs12s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i12 + zHen2t1 * zt2t1i12) 1742 zs21s0 = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i21 + zHen2t1 * zt2t1i21) 1743 s12ks=0.5_wp*(zs12s0+zs21s0) 1631 1744 1632 1745 END FUNCTION s12ks 1633 1746 1634 FUNCTION s22ks( x,y,z,phi)1747 FUNCTION s22ks(px,py,pz) 1635 1748 !!------------------------------------------------------------------- 1636 1749 !! Function : s22ks 1637 1750 !!------------------------------------------------------------------- 1638 REAL(wp), INTENT(in ) :: x,y,z,phi 1639 1640 REAL(wp) :: s22ks, p, pih 1641 1642 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21, n1t2i22 1643 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21, n2t1i22 1644 REAL(wp) :: t1t2i11, t1t2i12, t1t2i21, t1t2i22 1645 REAL(wp) :: t2t1i11, t2t1i12, t2t1i21, t2t1i22 1646 REAL(wp) :: d11, d12, d22 1647 REAL(wp) :: IIn1t2, IIn2t1, IIt1t2 1648 REAL(wp) :: Hen1t2, Hen2t1 1649 !!------------------------------------------------------------------- 1650 pih = 0.5_wp*rpi 1651 p =phi 1652 1653 n1t2i11 = cos(z+pih-p) * cos(z+p) 1654 n1t2i12 = cos(z+pih-p) * sin(z+p) 1655 n1t2i21 = sin(z+pih-p) * cos(z+p) 1656 n1t2i22 = sin(z+pih-p) * sin(z+p) 1657 n2t1i11 = cos(z-pih+p) * cos(z-p) 1658 n2t1i12 = cos(z-pih+p) * sin(z-p) 1659 n2t1i21 = sin(z-pih+p) * cos(z-p) 1660 n2t1i22 = sin(z-pih+p) * sin(z-p) 1661 t1t2i11 = cos(z-p) * cos(z+p) 1662 t1t2i12 = cos(z-p) * sin(z+p) 1663 t1t2i21 = sin(z-p) * cos(z+p) 1664 t1t2i22 = sin(z-p) * sin(z+p) 1665 t2t1i11 = cos(z+p) * cos(z-p) 1666 t2t1i12 = cos(z+p) * sin(z-p) 1667 t2t1i21 = sin(z+p) * cos(z-p) 1668 t2t1i22 = sin(z+p) * sin(z-p) 1669 d11 = cos(y)*cos(y)*(cos(x)+sin(x)*tan(y)*tan(y)) 1670 d12 = cos(y)*cos(y)*tan(y)*(-cos(x)+sin(x)) 1671 d22 = cos(y)*cos(y)*(sin(x)+cos(x)*tan(y)*tan(y)) 1672 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 * d22 1673 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 * d22 1674 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 * d22 1675 1676 IF (-IIn1t2>=rsmall) THEN 1677 Hen1t2 = 1._wp 1751 REAL(wp), INTENT(in ) :: px,py,pz 1752 1753 REAL(wp) :: s22ks, zpih 1754 1755 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1756 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1757 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1758 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1759 REAL(wp) :: zd11, zd12, zd22 1760 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1761 REAL(wp) :: zHen1t2, zHen2t1 1762 !!------------------------------------------------------------------- 1763 zpih = 0.5_wp*rpi 1764 1765 zn1t2i11 = cos(pz+zpih-pphi) * cos(pz+pphi) 1766 zn1t2i12 = cos(pz+zpih-pphi) * sin(pz+pphi) 1767 zn1t2i21 = sin(pz+zpih-pphi) * cos(pz+pphi) 1768 zn1t2i22 = sin(pz+zpih-pphi) * sin(pz+pphi) 1769 zn2t1i11 = cos(pz-zpih+pphi) * cos(pz-pphi) 1770 zn2t1i12 = cos(pz-zpih+pphi) * sin(pz-pphi) 1771 zn2t1i21 = sin(pz-zpih+pphi) * cos(pz-pphi) 1772 zn2t1i22 = sin(pz-zpih+pphi) * sin(pz-pphi) 1773 zt1t2i11 = cos(pz-pphi) * cos(pz+pphi) 1774 zt1t2i12 = cos(pz-pphi) * sin(pz+pphi) 1775 zt1t2i21 = sin(pz-pphi) * cos(pz+pphi) 1776 zt1t2i22 = sin(pz-pphi) * sin(pz+pphi) 1777 zt2t1i11 = cos(pz+pphi) * cos(pz-pphi) 1778 zt2t1i12 = cos(pz+pphi) * sin(pz-pphi) 1779 zt2t1i21 = sin(pz+pphi) * cos(pz-pphi) 1780 zt2t1i22 = sin(pz+pphi) * sin(pz-pphi) 1781 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1782 zd12 = cos(py)*cos(py)*tan(py)*(-cos(px)+sin(px)) 1783 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1784 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1785 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1786 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1787 1788 IF (-zIIn1t2>=rsmall) THEN 1789 zHen1t2 = 1._wp 1678 1790 ELSE 1679 Hen1t2 = 0._wp1680 ENDIF 1681 1682 IF (- IIn2t1>=rsmall) THEN1683 Hen2t1 = 1._wp1791 zHen1t2 = 0._wp 1792 ENDIF 1793 1794 IF (-zIIn2t1>=rsmall) THEN 1795 zHen2t1 = 1._wp 1684 1796 ELSE 1685 Hen2t1 = 0._wp1686 ENDIF 1687 1688 s22ks = sign(1._wp, IIt1t2+rsmall)*(Hen1t2 * t1t2i22 + Hen2t1 *t2t1i22)1797 zHen2t1 = 0._wp 1798 ENDIF 1799 1800 s22ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22) 1689 1801 1690 1802 END FUNCTION s22ks
Note: See TracChangeset
for help on using the changeset viewer.
