C------------------------------------------------------------
SUBROUTINE XSTPC(III,SCHGPC)
IMPLICIT REAL(A-H,O-Y)
IMPLICIT COMPLEX(Z)
INCLUDE 'PACVPP'
C STRESS FOR PARTIAL CORE CORRECTION
C-----ARRAYS FOR MFFT--------------------------------------------
DIMENSION IWL(8*IFX2+28),IWM(8*IFY2+28),IWN(8*IFZ2+28)
& ,IWORK(2*IFX2)
DIMENSION ZRB(IFX2,IFY2,IFZ2),ZRX(KFFT)
& ,ZRY(KFFT),ZRZ(KFFT)
C============================================================
EQUIVALENCE (ZV1D(1),ZRB(1,1,1))
C============================================================
ICOUNT=0
KFT1 = IFX2
KFT2 = IFY2
KFT3 = IFZ2
KSUM=KFT1*KFT2*KFT3
KVOL=(KFT1-1)*KFT2*KFT3
THIRD=1.D0/3.D0
ALFF=0.7D0
CO1=1.5D0*ALFF*(3.D0/PAI)**THIRD
CO2=(0.75D0/PAI)**THIRD
CO3=4.D0/3.D0
VCEL=ABS(UNIVOL/FLOAT(KVOL))
C*****---- SET UP C3FFT (FAST FOURIER TRANSFORMATION) -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,0,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(SET UP)] IERR = ',IERR
STOP
END IF
C----- X ---------------------------
DO 111 I=1,KSUM
ZV1D(I)=DCMPLX(0.0D0,0.0D0)
111 CONTINUE
DO 110 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZRB(L1,L2,L3)= ZRRPC(I)*GX(I)*GX(I)
110 CONTINUE
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,1,1,IWORK,IERRI)
IF (IERRI.NE.0) THEN
WRITE (6,*) ' C3FFT(INVERSE)] IERR = ',IERRI
STOP
END IF
DO 1005 I=1,KSUM
ZRX(I) = ZV1D(I)
1005 CONTINUE
C----- Y ---------------------------
DO 113 I=1,KSUM
ZV1D(I)=DCMPLX(0.0D0,0.0D0)
113 CONTINUE
DO 134 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZRB(L1,L2,L3)= ZRRPC(I)*GY(I)*GY(I)
134 CONTINUE
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,1,1,IWORK,IERRI)
IF (IERRI.NE.0) THEN
WRITE (6,*) ' C3FFT(INVERSE)] IERR = ',IERRI
STOP
END IF
DO 1006 I=1,KSUM
ZRY(I) = ZV1D(I)
1006 CONTINUE
C----- Z ---------------------------
DO 114 I=1,KSUM
ZV1D(I)=DCMPLX(0.0D0,0.0D0)
114 CONTINUE
DO 124 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZRB(L1,L2,L3)= ZRRPC(I)*GZ(I)*GZ(I)
124 CONTINUE
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,1,1,IWORK,IERRI)
IF (IERRI.NE.0) THEN
WRITE (6,*) ' C3FFT(INVERSE)] IERR = ',IERRI
STOP
END IF
DO 1007 I=1,KSUM
ZRZ(I) = ZV1D(I)
1007 CONTINUE
C-------------------------------
DO 112 I=1,KSUM
ZV1D(I)=DCMPLX(0.0D0,0.0D0)
112 CONTINUE
DO 210 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZRB(L1,L2,L3)=ZCHG(I)+ZRHPC(I)
210 CONTINUE
C*****---- INVERSE FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,1,1,IWORK,IERRI)
IF (IERRI.NE.0) THEN
WRITE (6,*) ' C3FFT(INVERSE)] IERR = ',IERRI
STOP
END IF
DO 1000 I=1,KSUM
CHGB1(I) = DREAL(ZV1D(I))
1000 CONTINUE
DO 20 I=1,KSUM
IF( CHGB1(I).LE.0.0D0 ) THEN
IF ( MOD(I,IFX2).NE.0 ) THEN
IF ( CHGB1(I).LE.-1.0D-5 ) THEN
ICOUNT=ICOUNT+1
C WRITE(6,610) I,CHGB1(I)
C 610 FORMAT(1H ,'**WARNING CHG.DEN<0.0 AT ',
C & I5,2D15.7,'***')
END IF
CHGB1(I) = 1.D-40
ELSE
CHGB1(I) = 1.D-40
END IF
END IF
20 CONTINUE
S = 0.D0
DO 1010 I=1,KSUM
S = S+CHGB1(I)
1010 CONTINUE
S=S*VCEL - SCHGPC
WRITE (6,*) 'REAL TOTAL CHARGE = ',S,' IN XSTPC'
C -- X ------------
DO 1020 I=1,KSUM
RS=CO2*( 1.0D0/CHGB1(I) )**THIRD
RH1 = -0.458D0*CO3/RS-(0.44D0*CO3*RS + 3.432D0)/
& (RS+7.8D0)**2
RH2 = -0.458D0/RS-0.44D0/(RS+7.8D0)
C
RHV(I)=RH1 - RH2
ZV1D(I)=ZRX(I)*(RHV(I))/CHGB1(I)
1020 CONTINUE
C*****---- FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,-1,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(FFT)] IERR = ',IERR
STOP
END IF
DENOM=1.0D0/DBLE(KVOL)
DO 30 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZXXX(I)=ZRB(L1,L2,L3)*DENOM
30 CONTINUE
C -- Y -------------
DO 1021 I=1,KSUM
ZV1D(I)=ZRY(I)*(RHV(I))/CHGB1(I)
1021 CONTINUE
C*****---- FAST FOURIER TRANSFORMATION -----
ZYYY(I)=ZRB(L1,L2,L3)*DENOM
31 CONTINUE
C -- Z -------------
DO 1022 I=1,KSUM
ZV1D(I)=ZRZ(I)*(RHV(I))/CHGB1(I)
1022 CONTINUE
C*****---- FAST FOURIER TRANSFORMATION -----
CALL C3FFT(ZRB,KFT1,KFT1-1,KFT2,KFT3,IWL,IWM,IWN
& ,0,-1,1,IWORK,IERR)
IF (IERR.NE.0) THEN
WRITE (6,*) ' C3FFT(FFT)] IERR = ',IERR
STOP
END IF
DO 32 I=1,KG
L1=IGF1(I)
L2=IGF2(I)
L3=IGF3(I)
ZZZZ(I)=ZRB(L1,L2,L3)*DENOM
32 CONTINUE
WRITE (6,*) 'ICOUNT = ',ICOUNT,' XSTPC END'
RETURN
END