C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ 5/1, 2009
SUBROUTINE KBINT
C^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
CCC!$OMP THREADPRIVATE( /PSSNL/, /DIV/, /WORK8/, /XY123/ )
C!$OMP THREADPRIVATE( /PSSNL/ )
C
DIMENSION XL(KNG1),YL1(KNG1),YL2(KNG1),YL3(KNG1)
& , YL4(KNG1),YLD(KNG1)
& ,QLX(KNG1),QLY(KNG1),QLZ(KNG1)
& ,QLA1(KNG1),QLA2(KNG1),QLA3(KNG1)
& ,QLA4(KNG1),QLA5(KNG1)
& ,ALK2(KNG1)
DIMENSION SNP(KNG1,KNV3,KTYP,10)
& ,SNP2(KNG1,KNV3,KTYP,9)
& ,SNP3(KNG1,KNV3,KTYP,4)
CX,Y1,Y2,Y3,Y4,YD,QX,QY,QZ, QA1,QA2,QA3,QA4,QA5, AK2
C
C K-B TYPE ADDITION
C INTEGRATION OF K-B TYPE POTENTIAL
C MODIFIED TO STRESS CALCULATION FOR KB-SEPARABLE FORM
C
NSEK=2
CO0=DSQRT(2.0D0)
CO1=DSQRT(3.0D0)
SC=DSQRT(PAI4)
PC=DSQRT(3.0D0*PAI4)
DC=DSQRT(5.0D0*PAI4)
C
C WRITE (6,*) 'CHECK POINT 0 IN KBINT'
C
C DO 1003 IK=1,KV3
C WRITE (6,*) 'CHECK POINT 01 IN KBINT, IK = ',IK
C AKX = VX(IK)
C AKY = VY(IK)
C AKZ = VZ(IK)
C
C DO 1004 I=1,KNG1
C RAK(I,IK) = 0.0D0
C 1004 CONTINUE
C IIBA = IBA(IK)
C DO 1005 I=1,IBA(IK)
C RK=SQRT((AKX+GX(NBASE(I,IK)))**2
C & +(AKY+GY(NBASE(I,IK)))**2 + (AKZ+GZ(NBASE(I,IK)))**2 )
C IF (RK.NE.0.0D0) THEN
C RAK(I,IK)=1.0D0/RK
C ELSE
C RAK(I,IK)=0.0D0
C END IF
C 1005 CONTINUE
C 1003 CONTINUE
C
!$OMP PARALLEL DEFAULT(NONE)
C!$OMP& COPYIN( SNL,SNL2,SNL3 )
!$OMP& FIRSTPRIVATE( NSEK, CO1, SC,PC,DC )
!$OMP& PRIVATE(XL,YL1,YL2,YL3,YL4,YLD,QLX,QLY,QLZ,
!$OMP& QLA1,QLA2,QLA3,QLA4,QLA5, ALK2,
!$OMP& AKX,AKY,AKZ, R, RK, IIBA,
!$OMP& FACS,FACP,FACD, FAC1,FAC2,FAC3 )
!$OMP& SHARED( KV3, RAK, SNP,SNP2,SNP3,
!$OMP& VX,VY,VZ,IBA,NBASE,
!$OMP& MESHR,WVS,WVP,WVD, DX,RAD, OMO,
!$OMP& GX,GY,GZ, NLSPD )
!$OMP DO
DO 1006 NNN=1,KV3
C WRITE (6,*) 'CHECK POINT 1 IN KBINT, NNN,DX = ',NNN,DX
AKX = VX(NNN)
AKY = VY(NNN)
AKZ = VZ(NNN)
C DO 1225 I = 1,KNG1
C XL(I) = 0.0D0
C YL1(I) = 0.0D0
C YL2(I) = 0.0D0
C YL3(I) = 0.0D0
C YL4(I) = 0.0D0
C YLD(I) = 0.0D0
C QLX(I) = 0.0D0
C QLY(I) = 0.0D0
C QLZ(I) = 0.0D0
C QLA1(I) = 0.0D0
C QLA2(I) = 0.0D0
C QLA3(I) = 0.0D0
C QLA4(I) = 0.0D0
C QLA5(I) = 0.0D0
C ALK2(I) = 0.0D0
C 1225 CONTINUE
C IIBA = IBA(NNN)
DO 1004 I=1,KNG1
RAK(I,NNN) = 0.0D0
1004 CONTINUE
DO 1005 I=1,IBA(NNN)
RK=SQRT((AKX+GX(NBASE(I,NNN)))**2
& +(AKY+GY(NBASE(I,NNN)))**2 + (AKZ+GZ(NBASE(I,NNN)))**2 )
IF (RK.NE.0.0D0) THEN
RAK(I,NNN)=1.0D0/RK
ELSE
RAK(I,NNN)=0.0D0
END IF
1005 CONTINUE
C WRITE (6,*) 'CHECK POINT 2 IN KBINT, NNN = ',NNN
IIBA = IBA(NNN)
C
DO 7500 J =1,IIBA
QLX(J)=AKX+GX(NBASE(J,NNN))
QLY(J)=AKY+GY(NBASE(J,NNN))
QLZ(J)=AKZ+GZ(NBASE(J,NNN))
QLA1(J)=0.5D0*(2.0D0*QLZ(J)*QLZ(J)-QLX(J)*QLX(J)-QLY(J)*QLY(J))
QLA2(J)=0.5D0*CO1*(QLX(J)*QLX(J)-QLY(J)*QLY(J))
QLA3(J)=CO1*QLX(J)*QLY(J)
QLA4(J)=CO1*QLX(J)*QLZ(J)
QLA5(J)=CO1*QLY(J)*QLZ(J)
ALK2(J) = SQRT( QLX(J)**2 + QLY(J)**2 + QLZ(J)**2 )
7500 CONTINUE
C
DO 7300 ITY=1,KTYP
C ATTN] IL=1,4 (S,P ONLY)
DO 7311 KN=1,KNG1
SNP(KN,NNN,ITY,1)=0.0D0
SNP(KN,NNN,ITY,2)=0.0D0
SNP(KN,NNN,ITY,3)=0.0D0
SNP(KN,NNN,ITY,4)=0.0D0
SNP2(KN,NNN,ITY,1)=0.0D0
SNP2(KN,NNN,ITY,2)=0.0D0
SNP2(KN,NNN,ITY,3)=0.0D0
SNP2(KN,NNN,ITY,4)=0.0D0
SNP3(KN,NNN,ITY,1)=0.0D0
7311 CONTINUE
C
DO 7320 N=1,MESHR
R=RAD(N)
FACS= OMO(N)*WVS(N,ITY)*R**NSEK*SC*DX
FACP= OMO(N)*WVP(N,ITY)*R**NSEK*PC*DX
FAC1=-OMO(N)*WVS(N,ITY)*R**(NSEK+1)*SC*DX
FAC2= OMO(N)*WVP(N,ITY)*R**(NSEK+1)*PC*DX
DO 7777 KN=1,IIBA
XL(KN)=ALK2(KN)*R
7777 CONTINUE
CALL DSJNV(0,IIBA,XL,YL1)
CALL DSJNV(1,IIBA,XL,YL2)
CALL DSJNV(2,IIBA,XL,YL3)
DO 7878 KN=1,IIBA
YL3(KN)=(YL1(KN)-2.0D0*YL3(KN))/3.0D0
7878 CONTINUE
DO 7330 KN=1,IIBA
C RKS(N)=$DMSJM(0,AK2(KN)*RAD(N))
SNP(KN,NNN,ITY,1) =SNP(KN,NNN,ITY,1) +FACS*YL1(KN)
SNP(KN,NNN,ITY,2) =SNP(KN,NNN,ITY,2)
& +FACP*YL2(KN)*QLX(KN)*RAK(KN,NNN)
SNP(KN,NNN,ITY,3) =SNP(KN,NNN,ITY,3)
& +FACP*YL2(KN)*QLY(KN)*RAK(KN,NNN)
SNP(KN,NNN,ITY,4) =SNP(KN,NNN,ITY,4)
& +FACP*YL2(KN)*QLZ(KN)*RAK(KN,NNN)
SNP2(KN,NNN,ITY,1) =SNP2(KN,NNN,ITY,1)+FAC1*YL2(KN)
SNP2(KN,NNN,ITY,2) =SNP2(KN,NNN,ITY,2)
& +FAC2*YL3(KN)*QLX(KN)*RAK(KN,NNN)
SNP2(KN,NNN,ITY,3) =SNP2(KN,NNN,ITY,3)
& +FAC2*YL3(KN)*QLY(KN)*RAK(KN,NNN)
SNP2(KN,NNN,ITY,4) =SNP2(KN,NNN,ITY,4)
& +FAC2*YL3(KN)*QLZ(KN)*RAK(KN,NNN)
SNP3(KN,NNN,ITY,1) =SNP3(KN,NNN,ITY,1)+FACP*YL2(KN)
7330 CONTINUE
7320 CONTINUE
C
IF (NLSPD(ITY).EQ.2) THEN
C
DO 7314 KN=1,KNG1
SNP(KN,NNN,ITY,5)=0.0D0
SNP(KN,NNN,ITY,6)=0.0D0
SNP(KN,NNN,ITY,7)=0.0D0
SNP(KN,NNN,ITY,8)=0.0D0
SNP(KN,NNN,ITY,9)=0.0D0
SNP(KN,NNN,ITY,10)=0.0D0
SNP2(KN,NNN,ITY,5)=0.0D0
SNP2(KN,NNN,ITY,6)=0.0D0
SNP2(KN,NNN,ITY,7)=0.0D0
SNP2(KN,NNN,ITY,8)=0.0D0
SNP2(KN,NNN,ITY,9)=0.0D0
SNP3(KN,NNN,ITY,2)=0.0D0
SNP3(KN,NNN,ITY,3)=0.0D0
SNP3(KN,NNN,ITY,4)=0.0D0
7314 CONTINUE
C
DO 7321 N=1,MESHR
R=RAD(N)
FACS= OMO(N)*WVS(N,ITY)*R**NSEK*SC*DX
FACP= OMO(N)*WVP(N,ITY)*R**NSEK*PC*DX
FACD= OMO(N)*WVD(N,ITY)*R**NSEK*DC*DX
FAC1=-OMO(N)*WVS(N,ITY)*R**(NSEK+1)*SC*DX
FAC2= OMO(N)*WVP(N,ITY)*R**(NSEK+1)*PC*DX
FAC3= OMO(N)*WVD(N,ITY)*R**(NSEK+1)*DC*DX
DO 8778 KN=1,IIBA
XL(KN)=ALK2(KN)*R
8778 CONTINUE
CALL DSJNV(0,IIBA,XL,YL1)
CALL DSJNV(1,IIBA,XL,YL2)
CALL DSJNV(2,IIBA,XL,YL3)
CALL DSJNV(3,IIBA,XL,YL4)
DO 7788 KN=1,IIBA
YLD(KN)= YL3(KN)
YL3(KN)=(YL1(KN)-2.0D0*YL3(KN))/3.0D0
YL4(KN)=(2.0D0*YL2(KN)-3.0D0*YL4(KN))/5.0D0
7788 CONTINUE
DO 7331 KN=1,IIBA
C ATTN] FOR SSL2 (IN ISSP)
C---*----1----*----2----*----3----*----4----*----5----*----6----*----7
SNP(KN,NNN,ITY,5) =SNP(KN,NNN,ITY,5)
& +FACD*YLD(KN)*QLA1(KN)*RAK(KN,NNN)*RAK(KN,NNN)
SNP(KN,NNN,ITY,6) =SNP(KN,NNN,ITY,6)
& +FACD*YLD(KN)*QLA2(KN)*RAK(KN,NNN)*RAK(KN,NNN)
SNP(KN,NNN,ITY,7) =SNP(KN,NNN,ITY,7)
& +FACD*YLD(KN)*QLA3(KN)*RAK(KN,NNN)*RAK(KN,NNN)
SNP(KN,NNN,ITY,8) =SNP(KN,NNN,ITY,8)
& +FACD*YLD(KN)*QLA4(KN)*RAK(KN,NNN)*RAK(KN,NNN)
SNP(KN,NNN,ITY,9) =SNP(KN,NNN,ITY,9)
& +FACD*YLD(KN)*QLA5(KN)*RAK(KN,NNN)*RAK(KN,NNN)
SNP(KN,NNN,ITY,10)=SNP(KN,NNN,ITY,10)+FACD*YLD(KN)
SNP2(KN,NNN,ITY,5)=SNP2(KN,NNN,ITY,5)
& +FAC3*YL4(KN)*QLA1(KN)*RAK(KN,NNN)*RAK(KN,NNN)
SNP2(KN,NNN,ITY,6)=SNP2(KN,NNN,ITY,6)
& +FAC3*YL4(KN)*QLA2(KN)*RAK(KN,NNN)*RAK(KN,NNN)
SNP2(KN,NNN,ITY,7)=SNP2(KN,NNN,ITY,7)
& +FAC3*YL4(KN)*QLA3(KN)*RAK(KN,NNN)*RAK(KN,NNN)
SNP2(KN,NNN,ITY,8)=SNP2(KN,NNN,ITY,8)
& +FAC3*YL4(KN)*QLA4(KN)*RAK(KN,NNN)*RAK(KN,NNN)
SNP2(KN,NNN,ITY,9)=SNP2(KN,NNN,ITY,9)
& +FAC3*YL4(KN)*QLA5(KN)*RAK(KN,NNN)*RAK(KN,NNN)
SNP3(KN,NNN,ITY,2)=SNP3(KN,NNN,ITY,2)+CO1*FACD*QLX(KN)*YLD(KN)
& *RAK(KN,NNN)
SNP3(KN,NNN,ITY,3)=SNP3(KN,NNN,ITY,3)+CO1*FACD*QLY(KN)*YLD(KN)
& *RAK(KN,NNN)
SNP3(KN,NNN,ITY,4)=SNP3(KN,NNN,ITY,4)+CO1*FACD*QLZ(KN)*YLD(KN)
& *RAK(KN,NNN)
7331 CONTINUE
7321 CONTINUE
END IF
7300 CONTINUE
1006 CONTINUE
!$OMP END DO
!$OMP END PARALLEL
C
C!$OMP DO
DO 4000 NNN=1,KV3
C WRITE (6,*) 'CHECK POINT 3 IN KBINT, NNN,DX = ',NNN
DO 4010 ITY=1,KTYP
DO 7351 KN=1,KNG1
SNL(KN,NNN,ITY,1)=SNP(KN,NNN,ITY,1)
SNL(KN,NNN,ITY,2)=SNP(KN,NNN,ITY,2)
SNL(KN,NNN,ITY,3)=SNP(KN,NNN,ITY,3)
SNL(KN,NNN,ITY,4)=SNP(KN,NNN,ITY,4)
SNL2(KN,NNN,ITY,1)=SNP2(KN,NNN,ITY,1)
SNL2(KN,NNN,ITY,2)=SNP2(KN,NNN,ITY,2)
SNL2(KN,NNN,ITY,3)=SNP2(KN,NNN,ITY,3)
SNL2(KN,NNN,ITY,4)=SNP2(KN,NNN,ITY,4)
SNL3(KN,NNN,ITY,1)=SNP3(KN,NNN,ITY,1)
7351 CONTINUE
C
IF (NLSPD(ITY).EQ.2) THEN
C WRITE (6,*) 'CHECK POINT 4 IN KBINT'
DO 7354 KN=1,KNG1
SNL(KN,NNN,ITY,5)=SNP(KN,NNN,ITY,5)
SNL(KN,NNN,ITY,6)=SNP(KN,NNN,ITY,6)
SNL(KN,NNN,ITY,7)=SNP(KN,NNN,ITY,7)
SNL(KN,NNN,ITY,8)=SNP(KN,NNN,ITY,8)
SNL(KN,NNN,ITY,9)=SNP(KN,NNN,ITY,9)
SNL(KN,NNN,ITY,10)=SNP(KN,NNN,ITY,10)
SNL2(KN,NNN,ITY,5)=SNP2(KN,NNN,ITY,5)
SNL2(KN,NNN,ITY,6)=SNP2(KN,NNN,ITY,6)
SNL2(KN,NNN,ITY,7)=SNP2(KN,NNN,ITY,7)
SNL2(KN,NNN,ITY,8)=SNP2(KN,NNN,ITY,8)
SNL2(KN,NNN,ITY,9)=SNP2(KN,NNN,ITY,9)
SNL3(KN,NNN,ITY,2)=SNP3(KN,NNN,ITY,2)
SNL3(KN,NNN,ITY,3)=SNP3(KN,NNN,ITY,3)
SNL3(KN,NNN,ITY,4)=SNP3(KN,NNN,ITY,4)
7354 CONTINUE
END IF
4010 CONTINUE
4000 CONTINUE
C!$OMP END DO
C!$OMP END PARALLEL
C
RETURN
END