Changeset 13155
 Timestamp:
 20200624T17:18:26+02:00 (3 years ago)
 Location:
 NEMO/branches/2019/dev_r11842_SI310_EAP/src/ICE
 Files:

 2 edited
Legend:
 Unmodified
 Added
 Removed

NEMO/branches/2019/dev_r11842_SI310_EAP/src/ICE/icedyn_rdgrft.F90
r12261 r13155 139 139 INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not 140 140 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 141 REAL(wp), DIMENSION(jpij) :: conv! 1D rdg_conv (if EAP rheology)141 REAL(wp), DIMENSION(jpij) :: zconv ! 1D rdg_conv (if EAP rheology) 142 142 ! 143 143 INTEGER, PARAMETER :: jp_itermax = 20 … … 177 177 ! just needed here 178 178 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i ) 179 CALL tab_2d_1d( npti, nptidx(1:npti), conv(1:npti) , rdg_conv )179 CALL tab_2d_1d( npti, nptidx(1:npti), zconv (1:npti) , rdg_conv ) 180 180 ! needed here and in the iteration loop 181 181 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i) ! zdivu is used as a work array here (no change in divu_i) … … 189 189 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji)  ABS( zdivu(ji) ) )  MIN( zdivu(ji), 0._wp ) 190 190 IF( ln_rhg_EVP ) closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji)  ABS( zdivu(ji) ) )  MIN( zdivu(ji), 0._wp ) 191 IF( ln_rhg_EAP ) closing_net(ji) = conv(ji)191 IF( ln_rhg_EAP ) closing_net(ji) = zconv(ji) 192 192 ! 193 193 IF( zdivu(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), zdivu(ji) ) ! make sure the closing rate is large enough 
NEMO/branches/2019/dev_r11842_SI310_EAP/src/ICE/icedyn_rhg_eap.F90
r13128 r13155 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 ! lookup table for calculating structure tensor … … 144 144 REAL(wp) :: zstressptmp, zstressmtmp, zstress12tmpF ! anisotropic stress tensor components 145 145 REAL(wp) :: zalphar, zalphas ! for mechanical redistribution 146 REAL(wp) :: mresult11, mresult12, z1dtevpkth,p5kth, z1_dtevp_A ! for structure tensor evolution146 REAL(wp) :: zmresult11, zmresult12, z1dtevpkth, zp5kth, z1_dtevp_A ! for structure tensor evolution 147 147 ! 148 148 REAL(wp), DIMENSION(jpi,jpj) :: zstress12tmp ! anisotropic stress tensor component for regridding 149 REAL(wp), DIMENSION(jpi,jpj) :: yield11, yield22,yield12 ! yield surface tensor for history149 REAL(wp), DIMENSION(jpi,jpj) :: zyield11, zyield22, zyield12 ! yield surface tensor for history 150 150 REAL(wp), DIMENSION(jpi,jpj) :: zp_delt ! P/delta at T points 151 151 REAL(wp), DIMENSION(jpi,jpj) :: zbeta ! beta coef from Kimmritz 2017 … … 258 258 z1_dtevp_A = z1_dtevp/10.0_wp 259 259 z1dtevpkth = 1._wp / (z1_dtevp_A + 0.00002_wp) 260 p5kth = 0.5_wp * 0.00002_wp260 zp5kth = 0.5_wp * 0.00002_wp 261 261 262 262 ! Ice strength … … 413 413 414 414 !  anisotropic stress calculation  ! 415 CALL update_stress_rdg (jter, nn_nevp, phi,zdiv, zdt, &415 CALL update_stress_rdg (jter, nn_nevp, zdiv, zdt, & 416 416 zdsT, paniso_11(ji,jj), paniso_12(ji,jj), & 417 417 zstressptmp, zstressmtmp, & … … 423 423 CALL calc_ffrac(zstressptmp, zstressmtmp, zstress12tmp(ji,jj), & 424 424 paniso_11(ji,jj), paniso_12(ji,jj), & 425 mresult11,mresult12)426 427 paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A + p5kth mresult11) * z1dtevpkth ! implicit428 paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A  mresult12) * z1dtevpkth ! implicit425 zmresult11, zmresult12) 426 427 paniso_11(ji,jj) = (paniso_11(ji,jj)*z1_dtevp_A + zp5kth  zmresult11) * z1dtevpkth ! implicit 428 paniso_12(ji,jj) = (paniso_12(ji,jj)*z1_dtevp_A  zmresult12) * z1dtevpkth ! implicit 429 429 ENDIF 430 430 431 431 432 432 IF (jter == nn_nevp) THEN 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)433 zyield11(ji,jj) = 0.5_wp * (zstressptmp + zstressmtmp) 434 zyield22(ji,jj) = 0.5_wp * (zstressptmp  zstressmtmp) 435 zyield12(ji,jj) = zstress12tmp(ji,jj) 436 436 prdg_conv(ji,jj) = min( zalphar, 0._wp) 437 437 ENDIF … … 484 484 END DO 485 485 END DO 486 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 486 487 487 488 !  Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)  ! … … 861 862 IF( iom_use('yield11') .OR. iom_use('yield12') .OR. iom_use('yield22')) THEN 862 863 863 CALL lbc_lnk_multi( 'icedyn_rhg_eap', yield11, 'T', 1., yield22, 'T', 1.,yield12, 'T', 1. )864 865 CALL iom_put( 'yield11', yield11 * zmsk00 )866 CALL iom_put( 'yield22', yield22 * zmsk00 )867 CALL iom_put( 'yield12', yield12 * zmsk00 )864 CALL lbc_lnk_multi( 'icedyn_rhg_eap', zyield11, 'T', 1., zyield22, 'T', 1., zyield12, 'T', 1. ) 865 866 CALL iom_put( 'yield11', zyield11 * zmsk00 ) 867 CALL iom_put( 'yield22', zyield22 * zmsk00 ) 868 CALL iom_put( 'yield12', zyield12 * zmsk00 ) 868 869 ENDIF 869 870 … … 931 932 END SUBROUTINE ice_dyn_rhg_eap 932 933 933 SUBROUTINE update_stress_rdg (ksub, ndte, phi, divu,tension, &934 shear, a11,a12, &935 stressp,stressm, &936 stress12, strength, &934 SUBROUTINE update_stress_rdg (ksub, kndte, pdivu, ptension, & 935 pshear, pa11, pa12, & 936 pstressp, pstressm, & 937 pstress12, strength, & 937 938 palphar, palphas) 938 939 !! 939 !! *** SUBROUTINE rhg_eap_rst***940 !! *** SUBROUTINE update_stress_rdg *** 940 941 !! 941 !! ** Purpose : Updates the depending on values of strain rate and structure942 !! tensor and for ksub=ndteit computes closing and sliding rate942 !! ** Purpose : Updates the stress depending on values of strain rate and structure 943 !! tensor and for the last subcycle step it computes closing and sliding rate 943 944 !! 944 INTEGER, INTENT(in ) :: ksub, ndte945 REAL(wp), INTENT(in ) :: phi,strength946 REAL(wp), INTENT(in ) :: divu, tension,shear947 REAL(wp), INTENT(in ) :: a11,a12948 REAL(wp), INTENT( out) :: stressp, stressm,stress12945 INTEGER, INTENT(in ) :: ksub, kndte 946 REAL(wp), INTENT(in ) :: strength 947 REAL(wp), INTENT(in ) :: pdivu, ptension, pshear 948 REAL(wp), INTENT(in ) :: pa11, pa12 949 REAL(wp), INTENT( out) :: pstressp, pstressm, pstress12 949 950 REAL(wp), INTENT( out) :: palphar, palphas 950 951 … … 973 974 zinvstressconviso = 1._wp/(1._wp+kfriction*kfriction) 974 975 975 zinvsin = 1._wp/sin(2._wp*p hi) * zinvstressconviso976 zinvsin = 1._wp/sin(2._wp*pphi) * zinvstressconviso 976 977 !now uses phi as set in higher code 977 978 … … 980 981 981 982 ! 1) structure tensor 982 za22 = 1._wp a11983 za22 = 1._wppa11 983 984 zQ11Q11 = 1._wp 984 985 zQ12Q12 = rsmall … … 987 988 ! gamma: angle between general coordiantes and principal axis of A 988 989 ! here Tan2gamma = 2 a12 / (a11  a22) 989 IF((ABS( a11  za22) > rsmall).OR.(ABS(a12) > rsmall)) THEN990 zAngle_denom_gamma = 1._wp/sqrt( ( a11  za22 )*(a11  za22) + &991 4._wp* a12*a12 )992 993 zQ11Q11 = 0.5_wp + ( a11  za22 )*0.5_wp*zAngle_denom_gamma !Cos^2994 zQ12Q12 = 0.5_wp  ( a11  za22 )*0.5_wp*zAngle_denom_gamma !Sin^2995 zQ11Q12 = a12*zAngle_denom_gamma !CosSin990 IF((ABS(pa11  za22) > rsmall).OR.(ABS(pa12) > rsmall)) THEN 991 zAngle_denom_gamma = 1._wp/sqrt( ( pa11  za22 )*( pa11  za22) + & 992 4._wp*pa12*pa12 ) 993 994 zQ11Q11 = 0.5_wp + ( pa11  za22 )*0.5_wp*zAngle_denom_gamma !Cos^2 995 zQ12Q12 = 0.5_wp  ( pa11  za22 )*0.5_wp*zAngle_denom_gamma !Sin^2 996 zQ11Q12 = pa12*zAngle_denom_gamma !CosSin 996 997 ENDIF 997 998 998 999 ! rotation Q*atemp*Q^T 999 zatempprime = zQ11Q11* a11 + 2.0_wp*zQ11Q12*a12 + zQ12Q12*za221000 zatempprime = zQ11Q11*pa11 + 2.0_wp*zQ11Q12*pa12 + zQ12Q12*za22 1000 1001 1001 1002 ! make first principal value the largest … … 1003 1004 1004 1005 ! 2) strain rate 1005 zdtemp11 = 0.5_wp*(divu + tension) 1006 zdtemp12 = shear*0.5_wp 1007 zdtemp22 = 0.5_wp*(divu  tension) 1008 1009 !WRITE(numout,*) 'div, tension, shear inside loop', ksub, divu, tension, shear 1006 zdtemp11 = 0.5_wp*(pdivu + ptension) 1007 zdtemp12 = pshear*0.5_wp 1008 zdtemp22 = 0.5_wp*(pdivu  ptension) 1010 1009 1011 1010 ! here Tan2alpha = 2 dtemp12 / (dtemp11  dtemp22) … … 1132 1131 zsig22 = strength*(zstemp22r + kfriction*zstemp22s) * zinvsin 1133 1132 1134 !WRITE(numout,*) 'strength inside loop', strength1135 1133 !WRITE(numout,*) 'principal axis stress inside loop', sig11, sig12, sig22 1136 1134 … … 1142 1140 zsgprm22 = zQ12Q12*zsig11 + zQ11Q11*zsig22 + 2._wp*zQ11Q12 *zsig12 1143 1141 1144 stressp = zsgprm11 + zsgprm22 1145 stress12 = zsgprm12 1146 stressm = zsgprm11  zsgprm22 1147 1148 !WRITE(numout,*) ' output inside loop', ksub, stressp 1142 pstressp = zsgprm11 + zsgprm22 1143 pstress12 = zsgprm12 1144 pstressm = zsgprm11  zsgprm22 1149 1145 1150 1146 ! Compute ridging and sliding functions in general coordinates 1151 1147 ! (Equation 11 in Tsamados 2013) 1152 IF (ksub == ndte) THEN1148 IF (ksub == kndte) THEN 1153 1149 zrotstemp11r = zQ11Q11*zstemp11r  2._wp*zQ11Q12* zstemp12r & 1154 1150 + zQ12Q12*zstemp22r … … 1175 1171 1176 1172 1177 SUBROUTINE calc_ffrac ( stressp, stressm,stress12, &1178 a11,a12, &1179 mresult11,mresult12)1173 SUBROUTINE calc_ffrac (pstressp, pstressm, pstress12, & 1174 pa11, pa12, & 1175 pmresult11, pmresult12) 1180 1176 !! 1181 1177 !! *** ROUTINE calc_ffrac *** … … 1185 1181 !! ** Method : Eq. 7: Ffrac = kf(AS) or = 0 depending on sigma_1 and sigma_2 1186 1182 !! 1187 REAL (wp), INTENT(in) :: stressp, stressm, stress12, a11,a121188 REAL (wp), INTENT(out) :: mresult11,mresult121183 REAL (wp), INTENT(in) :: pstressp, pstressm, pstress12, pa11, pa12 1184 REAL (wp), INTENT(out) :: pmresult11, pmresult12 1189 1185 1190 1186 ! local variables … … 1199 1195 !! 1200 1196 ! 1201 zsigma11 = 0.5_wp*( stressp+stressm)1202 zsigma12 = stress121203 zsigma22 = 0.5_wp*( stresspstressm)1197 zsigma11 = 0.5_wp*(pstressp+pstressm) 1198 zsigma12 = pstress12 1199 zsigma22 = 0.5_wp*(pstressppstressm) 1204 1200 1205 1201 ! Here's the change  no longer calculate gamma, … … 1233 1229 ! Pure divergence 1234 1230 IF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 >= 0.0_wp)) THEN 1235 mresult11 = 0.0_wp1236 mresult12 = 0.0_wp1231 pmresult11 = 0.0_wp 1232 pmresult12 = 0.0_wp 1237 1233 1238 1234 ! Unconfined compression: cracking of blocks not along the axial splitting … … 1240 1236 ! which leads to the loss of their shape, so we again model it through diffusion 1241 1237 ELSEIF ((zsigma_1 >= 0.0_wp).AND.(zsigma_2 < 0.0_wp)) THEN 1242 mresult11 = kfrac * (a11  zQ12Q12)1243 mresult12 = kfrac * (a12 + zQ11Q12)1238 pmresult11 = kfrac * (pa11  zQ12Q12) 1239 pmresult12 = kfrac * (pa12 + zQ11Q12) 1244 1240 1245 1241 ! Shear faulting 1246 1242 ELSEIF (zsigma_2 == 0.0_wp) THEN 1247 mresult11 = 0.0_wp1248 mresult12 = 0.0_wp1243 pmresult11 = 0.0_wp 1244 pmresult12 = 0.0_wp 1249 1245 ELSEIF ((zsigma_1 <= 0.0_wp).AND.(zsigma_1/zsigma_2 <= threshold)) THEN 1250 mresult11 = kfrac * (a11  zQ12Q12)1251 mresult12 = kfrac * (a12 + zQ11Q12)1246 pmresult11 = kfrac * (pa11  zQ12Q12) 1247 pmresult12 = kfrac * (pa12 + zQ11Q12) 1252 1248 1253 1249 ! Horizontal spalling 1254 1250 ELSE 1255 mresult11 = 0.0_wp1256 mresult12 = 0.0_wp1251 pmresult11 = 0.0_wp 1252 pmresult12 = 0.0_wp 1257 1253 ENDIF 1258 1254 … … 1341 1337 s11r(ix,iy,ia) = s11r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1342 1338 exp(w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1343 s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz ,phi)*dz/sin(2._wp*phi)1339 s11kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1344 1340 s12r(ix,iy,ia) = s12r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1345 1341 exp(w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1346 s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz ,phi)*dz/sin(2._wp*phi)1342 s12kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1347 1343 s22r(ix,iy,ia) = s22r(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1348 1344 exp(w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1349 s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz ,phi)*dz/sin(2._wp*phi)1345 s22kr(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1350 1346 s11s(ix,iy,ia) = s11s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1351 1347 exp(w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1352 s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz ,phi)*dz/sin(2._wp*phi)1348 s11ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1353 1349 s12s(ix,iy,ia) = s12s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1354 1350 exp(w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1355 s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz ,phi)*dz/sin(2._wp*phi)1351 s12ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1356 1352 s22s(ix,iy,ia) = s22s(ix,iy,ia) + 1*w1(ainit+ia*da)* & 1357 1353 exp(w2(ainit+ia*da)*(zinit+iz*dz)*(zinit+iz*dz))* & 1358 s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz ,phi)*dz/sin(2._wp*phi)1354 s22ks(xinit+ix*dx,yinit+iy*dy,zinit+iz*dz)*dz/sin(2._wp*pphi) 1359 1355 ENDDO 1360 1356 IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp … … 1365 1361 IF (abs(s22s(ix,iy,ia)) < eps6) s22s(ix,iy,ia) = 0._wp 1366 1362 ELSE 1367 s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp ,phi)/sin(2._wp*phi)1368 s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp ,phi)/sin(2._wp*phi)1369 s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp ,phi)/sin(2._wp*phi)1370 s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp ,phi)/sin(2._wp*phi)1371 s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp ,phi)/sin(2._wp*phi)1372 s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp ,phi)/sin(2._wp*phi)1363 s11r(ix,iy,ia) = 0.5_wp*s11kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1364 s12r(ix,iy,ia) = 0.5_wp*s12kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1365 s22r(ix,iy,ia) = 0.5_wp*s22kr(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1366 s11s(ix,iy,ia) = 0.5_wp*s11ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1367 s12s(ix,iy,ia) = 0.5_wp*s12ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1368 s22s(ix,iy,ia) = 0.5_wp*s22ks(xinit+ix*dx,yinit+iy*dy,0._wp)/sin(2._wp*pphi) 1373 1369 IF (abs(s11r(ix,iy,ia)) < eps6) s11r(ix,iy,ia) = 0._wp 1374 1370 IF (abs(s12r(ix,iy,ia)) < eps6) s12r(ix,iy,ia) = 0._wp … … 1437 1433 END FUNCTION w2 1438 1434 1439 FUNCTION s11kr(px,py,pz ,phi)1435 FUNCTION s11kr(px,py,pz) 1440 1436 !! 1441 1437 !! Function : s11kr 1442 1438 !! 1443 REAL(wp), INTENT(in ) :: px,py,pz ,phi1439 REAL(wp), INTENT(in ) :: px,py,pz 1444 1440 1445 1441 REAL(wp) :: s11kr, zpih … … 1456 1452 zpih = 0.5_wp*rpi 1457 1453 1458 zn1t2i11 = cos(pz+zpihp hi) * cos(pz+phi)1459 zn1t2i12 = cos(pz+zpihp hi) * sin(pz+phi)1460 zn1t2i21 = sin(pz+zpihp hi) * cos(pz+phi)1461 zn1t2i22 = sin(pz+zpihp hi) * sin(pz+phi)1462 zn2t1i11 = cos(pzzpih+p hi) * cos(pzphi)1463 zn2t1i12 = cos(pzzpih+p hi) * sin(pzphi)1464 zn2t1i21 = sin(pzzpih+p hi) * cos(pzphi)1465 zn2t1i22 = sin(pzzpih+p hi) * sin(pzphi)1466 zt1t2i11 = cos(pzp hi) * cos(pz+phi)1467 zt1t2i12 = cos(pzp hi) * sin(pz+phi)1468 zt1t2i21 = sin(pzp hi) * cos(pz+phi)1469 zt1t2i22 = sin(pzp hi) * sin(pz+phi)1470 zt2t1i11 = cos(pz+p hi) * cos(pzphi)1471 zt2t1i12 = cos(pz+p hi) * sin(pzphi)1472 zt2t1i21 = sin(pz+p hi) * cos(pzphi)1473 zt2t1i22 = sin(pz+p hi) * sin(pzphi)1454 zn1t2i11 = cos(pz+zpihpphi) * cos(pz+pphi) 1455 zn1t2i12 = cos(pz+zpihpphi) * sin(pz+pphi) 1456 zn1t2i21 = sin(pz+zpihpphi) * cos(pz+pphi) 1457 zn1t2i22 = sin(pz+zpihpphi) * sin(pz+pphi) 1458 zn2t1i11 = cos(pzzpih+pphi) * cos(pzpphi) 1459 zn2t1i12 = cos(pzzpih+pphi) * sin(pzpphi) 1460 zn2t1i21 = sin(pzzpih+pphi) * cos(pzpphi) 1461 zn2t1i22 = sin(pzzpih+pphi) * sin(pzpphi) 1462 zt1t2i11 = cos(pzpphi) * cos(pz+pphi) 1463 zt1t2i12 = cos(pzpphi) * sin(pz+pphi) 1464 zt1t2i21 = sin(pzpphi) * cos(pz+pphi) 1465 zt1t2i22 = sin(pzpphi) * sin(pz+pphi) 1466 zt2t1i11 = cos(pz+pphi) * cos(pzpphi) 1467 zt2t1i12 = cos(pz+pphi) * sin(pzpphi) 1468 zt2t1i21 = sin(pz+pphi) * cos(pzpphi) 1469 zt2t1i22 = sin(pz+pphi) * sin(pzpphi) 1474 1470 ! In expression of tensor d, with this formulatin d(x)=d(x+pi) 1475 1471 ! Solution, when diagonalizing always check sgn(a11a22) if > then keep x else … … 1498 1494 END FUNCTION s11kr 1499 1495 1500 FUNCTION s12kr(px,py,pz ,phi)1496 FUNCTION s12kr(px,py,pz) 1501 1497 !! 1502 1498 !! Function : s12kr 1503 1499 !! 1504 REAL(wp), INTENT(in ) :: px,py,pz ,phi1500 REAL(wp), INTENT(in ) :: px,py,pz 1505 1501 1506 1502 REAL(wp) :: s12kr, zs12r0, zs21r0, zpih … … 1516 1512 zpih = 0.5_wp*rpi 1517 1513 1518 zn1t2i11 = cos(pz+zpihp hi) * cos(pz+phi)1519 zn1t2i12 = cos(pz+zpihp hi) * sin(pz+phi)1520 zn1t2i21 = sin(pz+zpihp hi) * cos(pz+phi)1521 zn1t2i22 = sin(pz+zpihp hi) * sin(pz+phi)1522 zn2t1i11 = cos(pzzpih+p hi) * cos(pzphi)1523 zn2t1i12 = cos(pzzpih+p hi) * sin(pzphi)1524 zn2t1i21 = sin(pzzpih+p hi) * cos(pzphi)1525 zn2t1i22 = sin(pzzpih+p hi) * sin(pzphi)1526 zt1t2i11 = cos(pzp hi) * cos(pz+phi)1527 zt1t2i12 = cos(pzp hi) * sin(pz+phi)1528 zt1t2i21 = sin(pzp hi) * cos(pz+phi)1529 zt1t2i22 = sin(pzp hi) * sin(pz+phi)1530 zt2t1i11 = cos(pz+p hi) * cos(pzphi)1531 zt2t1i12 = cos(pz+p hi) * sin(pzphi)1532 zt2t1i21 = sin(pz+p hi) * cos(pzphi)1533 zt2t1i22 = sin(pz+p hi) * sin(pzphi)1514 zn1t2i11 = cos(pz+zpihpphi) * cos(pz+pphi) 1515 zn1t2i12 = cos(pz+zpihpphi) * sin(pz+pphi) 1516 zn1t2i21 = sin(pz+zpihpphi) * cos(pz+pphi) 1517 zn1t2i22 = sin(pz+zpihpphi) * sin(pz+pphi) 1518 zn2t1i11 = cos(pzzpih+pphi) * cos(pzpphi) 1519 zn2t1i12 = cos(pzzpih+pphi) * sin(pzpphi) 1520 zn2t1i21 = sin(pzzpih+pphi) * cos(pzpphi) 1521 zn2t1i22 = sin(pzzpih+pphi) * sin(pzpphi) 1522 zt1t2i11 = cos(pzpphi) * cos(pz+pphi) 1523 zt1t2i12 = cos(pzpphi) * sin(pz+pphi) 1524 zt1t2i21 = sin(pzpphi) * cos(pz+pphi) 1525 zt1t2i22 = sin(pzpphi) * sin(pz+pphi) 1526 zt2t1i11 = cos(pz+pphi) * cos(pzpphi) 1527 zt2t1i12 = cos(pz+pphi) * sin(pzpphi) 1528 zt2t1i21 = sin(pz+pphi) * cos(pzpphi) 1529 zt2t1i22 = sin(pz+pphi) * sin(pzpphi) 1534 1530 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1535 1531 zd12 = cos(py)*cos(py)*tan(py)*(cos(px)+sin(px)) … … 1557 1553 END FUNCTION s12kr 1558 1554 1559 FUNCTION s22kr(px,py,pz ,phi)1555 FUNCTION s22kr(px,py,pz) 1560 1556 !! 1561 1557 !! Function : s22kr 1562 1558 !! 1563 REAL(wp), INTENT(in ) :: px,py,pz ,phi1559 REAL(wp), INTENT(in ) :: px,py,pz 1564 1560 1565 1561 REAL(wp) :: s22kr, zpih … … 1575 1571 zpih = 0.5_wp*rpi 1576 1572 1577 zn1t2i11 = cos(pz+zpihp hi) * cos(pz+phi)1578 zn1t2i12 = cos(pz+zpihp hi) * sin(pz+phi)1579 zn1t2i21 = sin(pz+zpihp hi) * cos(pz+phi)1580 zn1t2i22 = sin(pz+zpihp hi) * sin(pz+phi)1581 zn2t1i11 = cos(pzzpih+p hi) * cos(pzphi)1582 zn2t1i12 = cos(pzzpih+p hi) * sin(pzphi)1583 zn2t1i21 = sin(pzzpih+p hi) * cos(pzphi)1584 zn2t1i22 = sin(pzzpih+p hi) * sin(pzphi)1585 zt1t2i11 = cos(pzp hi) * cos(pz+phi)1586 zt1t2i12 = cos(pzp hi) * sin(pz+phi)1587 zt1t2i21 = sin(pzp hi) * cos(pz+phi)1588 zt1t2i22 = sin(pzp hi) * sin(pz+phi)1589 zt2t1i11 = cos(pz+p hi) * cos(pzphi)1590 zt2t1i12 = cos(pz+p hi) * sin(pzphi)1591 zt2t1i21 = sin(pz+p hi) * cos(pzphi)1592 zt2t1i22 = sin(pz+p hi) * sin(pzphi)1573 zn1t2i11 = cos(pz+zpihpphi) * cos(pz+pphi) 1574 zn1t2i12 = cos(pz+zpihpphi) * sin(pz+pphi) 1575 zn1t2i21 = sin(pz+zpihpphi) * cos(pz+pphi) 1576 zn1t2i22 = sin(pz+zpihpphi) * sin(pz+pphi) 1577 zn2t1i11 = cos(pzzpih+pphi) * cos(pzpphi) 1578 zn2t1i12 = cos(pzzpih+pphi) * sin(pzpphi) 1579 zn2t1i21 = sin(pzzpih+pphi) * cos(pzpphi) 1580 zn2t1i22 = sin(pzzpih+pphi) * sin(pzpphi) 1581 zt1t2i11 = cos(pzpphi) * cos(pz+pphi) 1582 zt1t2i12 = cos(pzpphi) * sin(pz+pphi) 1583 zt1t2i21 = sin(pzpphi) * cos(pz+pphi) 1584 zt1t2i22 = sin(pzpphi) * sin(pz+pphi) 1585 zt2t1i11 = cos(pz+pphi) * cos(pzpphi) 1586 zt2t1i12 = cos(pz+pphi) * sin(pzpphi) 1587 zt2t1i21 = sin(pz+pphi) * cos(pzpphi) 1588 zt2t1i22 = sin(pz+pphi) * sin(pzpphi) 1593 1589 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1594 1590 zd12 = cos(py)*cos(py)*tan(py)*(cos(px)+sin(px)) … … 1614 1610 END FUNCTION s22kr 1615 1611 1616 FUNCTION s11ks(px,py,pz ,phi)1612 FUNCTION s11ks(px,py,pz) 1617 1613 !! 1618 1614 !! Function : s11ks 1619 1615 !! 1620 REAL(wp), INTENT(in ) :: px,py,pz ,phi1616 REAL(wp), INTENT(in ) :: px,py,pz 1621 1617 1622 1618 REAL(wp) :: s11ks, zpih … … 1632 1628 zpih = 0.5_wp*rpi 1633 1629 1634 zn1t2i11 = cos(pz+zpihp hi) * cos(pz+phi)1635 zn1t2i12 = cos(pz+zpihp hi) * sin(pz+phi)1636 zn1t2i21 = sin(pz+zpihp hi) * cos(pz+phi)1637 zn1t2i22 = sin(pz+zpihp hi) * sin(pz+phi)1638 zn2t1i11 = cos(pzzpih+p hi) * cos(pzphi)1639 zn2t1i12 = cos(pzzpih+p hi) * sin(pzphi)1640 zn2t1i21 = sin(pzzpih+p hi) * cos(pzphi)1641 zn2t1i22 = sin(pzzpih+p hi) * sin(pzphi)1642 zt1t2i11 = cos(pzp hi) * cos(pz+phi)1643 zt1t2i12 = cos(pzp hi) * sin(pz+phi)1644 zt1t2i21 = sin(pzp hi) * cos(pz+phi)1645 zt1t2i22 = sin(pzp hi) * sin(pz+phi)1646 zt2t1i11 = cos(pz+p hi) * cos(pzphi)1647 zt2t1i12 = cos(pz+p hi) * sin(pzphi)1648 zt2t1i21 = sin(pz+p hi) * cos(pzphi)1649 zt2t1i22 = sin(pz+p hi) * sin(pzphi)1630 zn1t2i11 = cos(pz+zpihpphi) * cos(pz+pphi) 1631 zn1t2i12 = cos(pz+zpihpphi) * sin(pz+pphi) 1632 zn1t2i21 = sin(pz+zpihpphi) * cos(pz+pphi) 1633 zn1t2i22 = sin(pz+zpihpphi) * sin(pz+pphi) 1634 zn2t1i11 = cos(pzzpih+pphi) * cos(pzpphi) 1635 zn2t1i12 = cos(pzzpih+pphi) * sin(pzpphi) 1636 zn2t1i21 = sin(pzzpih+pphi) * cos(pzpphi) 1637 zn2t1i22 = sin(pzzpih+pphi) * sin(pzpphi) 1638 zt1t2i11 = cos(pzpphi) * cos(pz+pphi) 1639 zt1t2i12 = cos(pzpphi) * sin(pz+pphi) 1640 zt1t2i21 = sin(pzpphi) * cos(pz+pphi) 1641 zt1t2i22 = sin(pzpphi) * sin(pz+pphi) 1642 zt2t1i11 = cos(pz+pphi) * cos(pzpphi) 1643 zt2t1i12 = cos(pz+pphi) * sin(pzpphi) 1644 zt2t1i21 = sin(pz+pphi) * cos(pzpphi) 1645 zt2t1i22 = sin(pz+pphi) * sin(pzpphi) 1650 1646 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1651 1647 zd12 = cos(py)*cos(py)*tan(py)*(cos(px)+sin(px)) … … 1671 1667 END FUNCTION s11ks 1672 1668 1673 FUNCTION s12ks(px,py,pz ,phi)1669 FUNCTION s12ks(px,py,pz) 1674 1670 !! 1675 1671 !! Function : s12ks 1676 1672 !! 1677 REAL(wp), INTENT(in ) :: px,py,pz ,phi1673 REAL(wp), INTENT(in ) :: px,py,pz 1678 1674 1679 1675 REAL(wp) :: s12ks, zs12s0, zs21s0, zpih … … 1689 1685 zpih = 0.5_wp*rpi 1690 1686 1691 zn1t2i11 = cos(pz+zpihp hi) * cos(pz+phi)1692 zn1t2i12 = cos(pz+zpihp hi) * sin(pz+phi)1693 zn1t2i21 = sin(pz+zpihp hi) * cos(pz+phi)1694 zn1t2i22 = sin(pz+zpihp hi) * sin(pz+phi)1695 zn2t1i11 = cos(pzzpih+p hi) * cos(pzphi)1696 zn2t1i12 = cos(pzzpih+p hi) * sin(pzphi)1697 zn2t1i21 = sin(pzzpih+p hi) * cos(pzphi)1698 zn2t1i22 = sin(pzzpih+p hi) * sin(pzphi)1699 zt1t2i11 = cos(pzp hi) * cos(pz+phi)1700 zt1t2i12 = cos(pzp hi) * sin(pz+phi)1701 zt1t2i21 = sin(pzp hi) * cos(pz+phi)1702 zt1t2i22 = sin(pzp hi) * sin(pz+phi)1703 zt2t1i11 = cos(pz+p hi) * cos(pzphi)1704 zt2t1i12 = cos(pz+p hi) * sin(pzphi)1705 zt2t1i21 = sin(pz+p hi) * cos(pzphi)1706 zt2t1i22 = sin(pz+p hi) * sin(pzphi)1687 zn1t2i11 = cos(pz+zpihpphi) * cos(pz+pphi) 1688 zn1t2i12 = cos(pz+zpihpphi) * sin(pz+pphi) 1689 zn1t2i21 = sin(pz+zpihpphi) * cos(pz+pphi) 1690 zn1t2i22 = sin(pz+zpihpphi) * sin(pz+pphi) 1691 zn2t1i11 = cos(pzzpih+pphi) * cos(pzpphi) 1692 zn2t1i12 = cos(pzzpih+pphi) * sin(pzpphi) 1693 zn2t1i21 = sin(pzzpih+pphi) * cos(pzpphi) 1694 zn2t1i22 = sin(pzzpih+pphi) * sin(pzpphi) 1695 zt1t2i11 = cos(pzpphi) * cos(pz+pphi) 1696 zt1t2i12 = cos(pzpphi) * sin(pz+pphi) 1697 zt1t2i21 = sin(pzpphi) * cos(pz+pphi) 1698 zt1t2i22 = sin(pzpphi) * sin(pz+pphi) 1699 zt2t1i11 = cos(pz+pphi) * cos(pzpphi) 1700 zt2t1i12 = cos(pz+pphi) * sin(pzpphi) 1701 zt2t1i21 = sin(pz+pphi) * cos(pzpphi) 1702 zt2t1i22 = sin(pz+pphi) * sin(pzpphi) 1707 1703 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1708 1704 zd12 = cos(py)*cos(py)*tan(py)*(cos(px)+sin(px)) … … 1730 1726 END FUNCTION s12ks 1731 1727 1732 FUNCTION s22ks(px,py,pz ,phi)1728 FUNCTION s22ks(px,py,pz) 1733 1729 !! 1734 1730 !! Function : s22ks 1735 1731 !! 1736 REAL(wp), INTENT(in ) :: px,py,pz ,phi1732 REAL(wp), INTENT(in ) :: px,py,pz 1737 1733 1738 1734 REAL(wp) :: s22ks, zpih 1739 1735 1740 REAL(wp) :: n1t2i11, n1t2i12, n1t2i21,n1t2i221741 REAL(wp) :: n2t1i11, n2t1i12, n2t1i21,n2t1i221742 REAL(wp) :: t1t2i11, t1t2i12, t1t2i21,t1t2i221743 REAL(wp) :: t2t1i11, t2t1i12, t2t1i21,t2t1i221744 REAL(wp) :: d11, d12,d221745 REAL(wp) :: IIn1t2, IIn2t1,IIt1t21746 REAL(wp) :: Hen1t2,Hen2t11736 REAL(wp) :: zn1t2i11, zn1t2i12, zn1t2i21, zn1t2i22 1737 REAL(wp) :: zn2t1i11, zn2t1i12, zn2t1i21, zn2t1i22 1738 REAL(wp) :: zt1t2i11, zt1t2i12, zt1t2i21, zt1t2i22 1739 REAL(wp) :: zt2t1i11, zt2t1i12, zt2t1i21, zt2t1i22 1740 REAL(wp) :: zd11, zd12, zd22 1741 REAL(wp) :: zIIn1t2, zIIn2t1, zIIt1t2 1742 REAL(wp) :: zHen1t2, zHen2t1 1747 1743 !! 1748 1744 zpih = 0.5_wp*rpi 1749 1745 1750 n1t2i11 = cos(pz+zpihphi) * cos(pz+phi)1751 n1t2i12 = cos(pz+zpihphi) * sin(pz+phi)1752 n1t2i21 = sin(pz+zpihphi) * cos(pz+phi)1753 n1t2i22 = sin(pz+zpihphi) * sin(pz+phi)1754 n2t1i11 = cos(pzzpih+phi) * cos(pzphi)1755 n2t1i12 = cos(pzzpih+phi) * sin(pzphi)1756 n2t1i21 = sin(pzzpih+phi) * cos(pzphi)1757 n2t1i22 = sin(pzzpih+phi) * sin(pzphi)1758 t1t2i11 = cos(pzphi) * cos(pz+phi)1759 t1t2i12 = cos(pzphi) * sin(pz+phi)1760 t1t2i21 = sin(pzphi) * cos(pz+phi)1761 t1t2i22 = sin(pzphi) * sin(pz+phi)1762 t2t1i11 = cos(pz+phi) * cos(pzphi)1763 t2t1i12 = cos(pz+phi) * sin(pzphi)1764 t2t1i21 = sin(pz+phi) * cos(pzphi)1765 t2t1i22 = sin(pz+phi) * sin(pzphi)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))1769 IIn1t2 = n1t2i11 * d11 + (n1t2i12 + n1t2i21) * d12 + n1t2i22 *d221770 IIn2t1 = n2t1i11 * d11 + (n2t1i12 + n2t1i21) * d12 + n2t1i22 *d221771 IIt1t2 = t1t2i11 * d11 + (t1t2i12 + t1t2i21) * d12 + t1t2i22 *d221772 1773 IF ( IIn1t2>=rsmall) THEN1774 Hen1t2 = 1._wp1746 zn1t2i11 = cos(pz+zpihpphi) * cos(pz+pphi) 1747 zn1t2i12 = cos(pz+zpihpphi) * sin(pz+pphi) 1748 zn1t2i21 = sin(pz+zpihpphi) * cos(pz+pphi) 1749 zn1t2i22 = sin(pz+zpihpphi) * sin(pz+pphi) 1750 zn2t1i11 = cos(pzzpih+pphi) * cos(pzpphi) 1751 zn2t1i12 = cos(pzzpih+pphi) * sin(pzpphi) 1752 zn2t1i21 = sin(pzzpih+pphi) * cos(pzpphi) 1753 zn2t1i22 = sin(pzzpih+pphi) * sin(pzpphi) 1754 zt1t2i11 = cos(pzpphi) * cos(pz+pphi) 1755 zt1t2i12 = cos(pzpphi) * sin(pz+pphi) 1756 zt1t2i21 = sin(pzpphi) * cos(pz+pphi) 1757 zt1t2i22 = sin(pzpphi) * sin(pz+pphi) 1758 zt2t1i11 = cos(pz+pphi) * cos(pzpphi) 1759 zt2t1i12 = cos(pz+pphi) * sin(pzpphi) 1760 zt2t1i21 = sin(pz+pphi) * cos(pzpphi) 1761 zt2t1i22 = sin(pz+pphi) * sin(pzpphi) 1762 zd11 = cos(py)*cos(py)*(cos(px)+sin(px)*tan(py)*tan(py)) 1763 zd12 = cos(py)*cos(py)*tan(py)*(cos(px)+sin(px)) 1764 zd22 = cos(py)*cos(py)*(sin(px)+cos(px)*tan(py)*tan(py)) 1765 zIIn1t2 = zn1t2i11 * zd11 + (zn1t2i12 + zn1t2i21) * zd12 + zn1t2i22 * zd22 1766 zIIn2t1 = zn2t1i11 * zd11 + (zn2t1i12 + zn2t1i21) * zd12 + zn2t1i22 * zd22 1767 zIIt1t2 = zt1t2i11 * zd11 + (zt1t2i12 + zt1t2i21) * zd12 + zt1t2i22 * zd22 1768 1769 IF (zIIn1t2>=rsmall) THEN 1770 zHen1t2 = 1._wp 1775 1771 ELSE 1776 Hen1t2 = 0._wp1777 ENDIF 1778 1779 IF ( IIn2t1>=rsmall) THEN1780 Hen2t1 = 1._wp1772 zHen1t2 = 0._wp 1773 ENDIF 1774 1775 IF (zIIn2t1>=rsmall) THEN 1776 zHen2t1 = 1._wp 1781 1777 ELSE 1782 Hen2t1 = 0._wp1783 ENDIF 1784 1785 s22ks = sign(1._wp, IIt1t2+rsmall)*(Hen1t2 * t1t2i22 + Hen2t1 *t2t1i22)1778 zHen2t1 = 0._wp 1779 ENDIF 1780 1781 s22ks = sign(1._wp,zIIt1t2+rsmall)*(zHen1t2 * zt1t2i22 + zHen2t1 * zt2t1i22) 1786 1782 1787 1783 END FUNCTION s22ks
Note: See TracChangeset
for help on using the changeset viewer.