C Copy of INTEGRAL FORTRAN WRITTEN BY Jay Rasaiah
C MODIFIED INPUT TO RUN ON IMAC by E.LOVEJOY June 18 2024, Rev 8.01
C
C Compiled with:
C gfortran -g -fcheck=all -std=legacy -o integral integral.f
C -fcheck=all enables array bound checks at run time
C **********************************************************************INT00020
C THIS IS VERSION 1.57 WHICH USES LESS MEMORY AND IS FASTER THAN *INT00030
C VERSION 1.55 AND CALC THE SHORT-RANGED DIR CORR FUNCTION(DCS) INT00040
C PROGRAM RESTRICTED TO SYMMETRICALLY CHARGED ELECTROLYTES. *INT00050
C SOLVES HNC EQ BY ITERATIONUSING ISEINGS FAST FOURIER TRANSFORM FFT3S*INT00060
C IF (CONC.EQ.0) READS A NEW SET OF ION AND WELL SIZE PARAMETERS. *INT00070
C IF REM=0.0, NO REMAINDER ADDED TO G(I,J), T(2,4), G4(I,J) *INT00080
C IF REM.NE.0.0, DEBYE-HUCKEL REMAINDER ADDED *INT00090
C PRIMITIVE MODEL PLUS TWO SQUARE WELLS *INT00100
C SPAC IS NEAREST FRACTION OF A(1,2) THAT IS CLOSEST TO *INT00110
C MAX(COEFF1*DEBYE/FNN,COEFF2*A(1,2)/FNN) *INT00120
C DELA = DLN(ACTCO)/DLNC, CALC FROM COMPRESSIBILITY EQUATION. *INT00130
C DELB = LIMIT(DELB) AS BETA GOES TO ZERO *INT00140
C BETA MEAS DEVIATION FROM ELECTRONEUTRALITY *INT00150
C OSM IS CALCULATED FROM THE PRESSURE EQUATION *INT00160
C XENERGY AND DFX/DP FROM INDEP EQNS ASSUMING NO TEMP AND PRESS DEP *INT00170
C OF THE SHORT RANGE POTENTIAL XENERGY IS IN CALS PER MOLE OF SOLUTE *INT00180
C XVOLUME IN CC PER MOLE OF SOLUTE *INT00190
C DRF IS THE DIELECTRIC RESPONSE FUNCTION FOR THE HNC EQN *INT00200
C DRDF IS THE DIELECTRIC RESPONSE FUNCTION FOR THE DHLL *INT00210
C CONVGA = (C(1)*G(1,1)-C(2)*G(2,2))/GM12-C(1)+C(2) *INT00220
C CONVGB = 100*(G(1,2)-G(2,1))/GM12 *INT00230
C DCS(=ZETA) IS THE SHORT-RANGED DIRECT CORRELATION FUNCTION *INT00240
C CONVERGENCE CRITERION USEFUL ONLY FOR UNSYMMETRICAL SYSTEMS *INT00250
C IF EX.EQ.2.0, PRINTS XR,DCS,H *INT00260
C IF EX.EQ.3.0, PRINTS XR,DCS,H AND CHARGE DENSITY *INT00270
C IF EX.EQ.0.0, QUITS *INT00280
C IF CHARGE = 0.0, REPEATS CALC WITH SAME SPACING AS CHARGED SYSTEM *INT00290
C BUT WITH Z1 = Z2 = 0.0.FICTITIOUS CHARGE DENSITY ALSO CALC *INT00300
C IF ICC.GT.0, ITERATION BEGINS WITH TAU OF PRECEDING CONC *INT00310
C IF ICC.LT.0, READ TAU GENERATED IN A PREVIOUS EXECUTION *INT00320
C IF ICC.EQ.0, TAU CAL IN LAMBDA APPROXIMATION *INT00330
C WRITES TAU ON NDSK2 IN ALL CASES *INT00340
C SIX SUBRUTNIES ADJ,TRAP,SIMPSON,FFT3S,SETFT,RSA
C NIONS is number positive and negative ions 2 for NaCl and ZnSO4
C **************************DIMENSIONS**********************************INT00350
INTEGER MNM
PARAMETER (MNM = 2048)
REAL*8 DATE,TIME INT00360
DIMENSION A(2,2)
DIMENSION A1(2,2)
DIMENSION B(80)
DIMENSION BA(2,2)
DIMENSION B1(2,2)
DIMENSION C(2)
DIMENSION CB(2,2)
DIMENSION CONV(MNM,2,2)
DIMENSION C1(2,2)
DIMENSION DC(MNM,2,2)
DIMENSION DRF(MNM)
DIMENSION DRDF(MNM)
DIMENSION E(2,2)
DIMENSION FNU(2)
DIMENSION F2(2)
DIMENSION F4(2)
DIMENSION G(2,2)
DIMENSION GA(2,2)
DIMENSION GB(2,2)
DIMENSION GC(2,2)
DIMENSION GM(2,2)
DIMENSION GW(2,2)
DIMENSION G4(2,2)
C DIMENSION H(MNM,2,2)
REAL H(MNM,2,2)
DIMENSION LINT(2,2)
DIMENSION MINT(2,2)
DIMENSION NINT(2,2)
DIMENSION PHI(MNM,2,2)
DIMENSION P1(MNM,2,2)
DIMENSION Q(MNM,2,2)
DIMENSION Q1(MNM,2,2)
DIMENSION RHO(2,MNM)
DIMENSION RQO(2,MNM)
DIMENSION R2(2)
DIMENSION R4(2)
DIMENSION T(2,2)
DIMENSION TAU(MNM,2,2)
DIMENSION TCONV(MNM,2,2)
DIMENSION TH(MNM,2,2)
DIMENSION TM(2,2)
DIMENSION TPHI(MNM,2,2)
DIMENSION TQ(MNM,2,2)
DIMENSION TTAU(MNM,2,2)
DIMENSION TW(2,2)
DIMENSION TZETA(MNM,2,2)
DIMENSION W(2,2)
C DIMENSION XR(MNM)
REAL XR(MNM)
DIMENSION Y(2)
DIMENSION Y1(MNM)
DIMENSION Z(2)
DIMENSION ZA(2)
DIMENSION ZETA(MNM,2,2)
EQUIVALENCE (H,TH,PHI,TPHI,TTAU)
EQUIVALENCE (ZETA,TZETA,RHO,DRF,DC)
EQUIVALENCE (TAU,P1,CONV,TCONV)
EQUIVALENCE (RQO,DRDF)
EQUIVALENCE (XR,Y1)
EQUIVALENCE (B(6),Z(1)),(B(8),D),(B(9),TEMP),(B(10),A1(1,1)),
*(B(14),COEFF1),(B(15),DEBYE),(B(16),SPAC),(B(17),FNN),(B(18),CONC)
*,(B(19),DELA),(B(20),OSM),(B(21),CONTAC),(B(22),REMAIN),
*(B(23),FITER),(B(24),G(1,1)),(B(28),T(1,1)),(B(32),FNGARB),
*(B(33),DELB),(B(34),DIFF),(B(35),CONVG1),(B(36),CONVG2),
*(B(37),GA(1,1)),(B(5),DT),(B(4),DP),(B(41),GB(1,1)),
*(B(45),G4(1,1)),(B(49),F2(1)),(B(51),F4(1)),(B(52),COEFF2)
COMMON /JAY/ H,XR
C******************************* Hunting for issues********************
C Gfortran has some other nifty flags for ID of various issues:
C -ffpe-trap=invalid,zero,overflow
C gfortran -g -ffpe-trap=invalid,zero,overflow -std=legacy -fcheck=all
C******************************* DATE TIME FIX ************************
CHARACTER*25 date_string
CHARACTER*100 command
command = 'date "+%Y-%m-%d %H:%M:%S" > date.txt'
CALL SYSTEM(command)
OPEN(UNIT=1, FILE='date.txt', STATUS='OLD')
READ(1, '(A)') date_string
CLOSE(1)
C ******************************* CONSTANTS ****************************INT00600
C PROGRAM CONSTANTS *INT00610
DATA PI/3.141593/ INT00620
PRINT *, date_string
C CALL TIMER(DATE,TIME,IVIRT,ITOT) INT00630
NRDR = 5 INT00640
NPRT = 7 INT00650
NPUN = 7 INT00660
NDSK = 8 INT00670
NDSK2 = 9 INT00680
NIONS = 2 INT00690
PI2 = 2.*PI INT00700
PI4 = 4.*PI INT00710
PI43 = PI4/3. INT00720
PIPI = PI*PI INT00730
C **********************************************************************INT00740
C CONSTANTS FOR SUBROUTINE FFT3S *INT01430
F0 = 0. INT01440
RL = SPAC1*FNN INT01450
C **********************************************************************INT01460
C*****used to stop crashing of H MINT PHI TPHI TCONV CONV ZETA*******
I = 1
J = 1
IF (J .LT. 1) THEN
J = J + 1
END IF
IF (I .LT. 1) THEN
I = I + 1
END IF
DO 10 I = 1,80 INT00750
10 B(I) = 0.0 INT00760
B(3) = 1.57 INT00770
C20 READ(NRDR,9001) M1 INT00780
C20 READ(NRDR,*) M1 INT00780
C READ(NRDR,*) (Z(I),FNU(I),I=1,NIONS)
M1 = 11
C Z(1) = 1.0
Z(1) = 2.0
Z(2) = -1.0 INT00790
ZE = Z(1) INT00800
ZA(1) = Z(1) INT00810
ZA(2) = Z(2) INT00820
FNU(1) = 1.0
FNU(2) = 2.0
C25 READ(NRDR,*) D,TEMP,DT,DP
D = 78.358
T = 298.15
TEMP = 298.15
DT = 0.0
DP = 0.0
C30 READ(NRDR,*) A(1,1),BA(1,1),CB(1,1),A(1,2),BA(1,2),CB(1,2),A(2, INT00840
C *2),BA(2,2),CB(2,2) INT00850
A(1,1) = 4.2
A(1,2) = 4.2
A(2,1) = A(1,2) INT00860
A(2,2) = A (1,1)
BA(1,1) = 0
BA(2,2) = BA(1,1)
BA(1,2) = 0.0
BA(2,1) = BA(1,2)
CB(1,1) = 0.0
CB(2,2) = CB(1,1)
CB(1,2) = 0.0 INT00870
CB(2,1) = CB(1,2) INT00880
C35 READ(NRDR,*) E(1,1),W(1,1),E(1,2),W(1,2),E(2,2),W(2,2) INT00890
E(1,1)= 0.0
E(1,2)= 0.0
E(2,1) = E(1,2) INT00900
E(2,2) = E(1,1) INT00900
W(1,1) = 0.0
W(1,2) = 0.0
W(2,1) = W(1,2)
W(2,2) = W(1,1)
C40 READ(NRDR,*)CONC,COEFF1,COEFF2,EX,REM,CHARGE,MAX,I2,NCD,IC, INT00920
C *NR,IR,ICC
CONC = 0.1 INT00930
COEFF1=1.0
COEFF2=1.0
EX = 3.0
REM = 1.0
CHARGE = 1.0
MAX = 6
I2 =2
ICC= 1
C NCD =
C IC =
IR = 1
NR = 5
C REWIND NDSK2 INT00940
C IF(ICC.LT.0)READ(NDSK2,9003)(((TAU(K,I,J),I=1,NIONS),J=1,NIONS),K=INT00950
C *1,N) INT00960
C REWIND NDSK2 INT00970
IF (EX.EQ.0.0) GO TO 5050 INT00980
C IF(CONC.EQ.0.0) GO TO 30 INT00990
9001 FORMAT(I3) INT01000
9002 FORMAT(8F10.5) INT01010
9003 FORMAT(4F10.5) INT01020
9004 FORMAT(9F5.2) INT01030
9005 FORMAT(6F10.5) INT01040
9006 FORMAT(1F10.7,5F5.2,7I5) INT01050
WRITE(NPRT,9100)B(3), date_string
9100 FORMAT(1H1, 38X, 'HNC VERSION ', F4.2, ' - IBM 370/148 - ', A25)
C WRITE(NPRT,9100)B(3),DATE,TIME INT01060
C9100 FORMAT(1H1,38X,'HNC VERSION ',F4.2,' - IBM 370/148 - ',A8,2X,A8) INT01070
C**** N=2**M1=2048 for M1=11, ***************************************
N=2**M1 INT01090
N2=N/2 INT01100
FNN = N INT01120
FC = 0.0 INT01130
FNO = 0.0 INT01140
CC = 0.0 INT01150
DO 50 I = 1,NIONS INT01160
FNO = FNU(I)+FNO INT01170
C(I) = FNU(I)*CONC*6.0233*1.E-04 INT01180
CC = C(I)+CC INT01190
50 FC = C(I)*Z(I)**2+FC
C********************************************************************** INT01200
C Electron charge e=4.8032x10E10 cm**3/2 g**1/2 s**-2 (esu)
C Boltzconst k=1.3803x10E-16 cm**2 s**2/K INT01200
C Bjerrum length |z+z+|e**2/DkT=16.714x1.E+04/DT
C***********************************************************************
FL1 = 16.7094*1.E+04/(D*TEMP) INT01210
FL = PI4*FL1 INT01220
FK = (FC*FL)**0.5 INT01230
FKFK = FK*FK INT01240
DEBYE = 1./FK INT01250
FKA = A(1,2)/DEBYE INT01260
C REWIND NDSK INT01270
C WRITE(NDSK,9101)A(1,1),A(1,2),A(2,2),D,TEMP,CONC,FKA,DEBYE,Z(1),Z(INT01280
C *2) INT01290
9101 FORMAT(3F8.5,2F8.4,3F10.6,2F5.1) INT01300
SPAC1 = COEFF1*DEBYE/FNN INT01310
SPAC2 = COEFF2*A(1,2)/FNN INT01320
IF(SPAC2.GT.SPAC1) SPAC1=SPAC2 INT01330
S1 = A(1,2)/SPAC1 INT01340
L = S1+1. INT01350
S2 = L INT01360
IF (S2-S1.GT.0.5) GO TO 60 INT01370
SPAC1 = A(1,2)/S2 INT01380
GO TO 70 INT01390
60 SPAC1 = A(1,2)/(S2-1.) INT01400
70 CONTINUE INT01410
C **********************************************************************INT01420
C CONSTANTS FOR SUBROUTINE FFT3S *INT01430
F0 = 0. INT01440
RL = SPAC1*FNN INT01450
C FX,FY and FZ are Debye Huckel remainders for ?
C **********************************************************************INT01460
FI = 1.0 INT01470
SPAC = SPAC1 INT01480
DO 80 K = 1,N INT01490
XR(K) = SPAC1*FI INT01500
80 FI = FI+1. INT01510
XN = XR(N) INT01520
FKN = FK*XN INT01530
IF (REM.EQ.0.0) GO TO 90 INT01540
FX = FL*(1.+FKN)*EXP(-FKN)/FKFK INT01550
FY = FL*EXP(-FKN)/FK INT01560
FZ = FL*EXP(-FKN)*(FKN**3+3.*FKN**2+6.*FKN+6.)/(FKFK*FKFK) INT01570
GO TO 100 INT01580
90 FX = 0.0 INT01590
FY = 0.0 INT01600
FZ = 0.0 INT01610
100 CONTINUE INT01620
DO 110 I = 1,NIONS INT01630
DO 110 J = 1,NIONS INT01640
CALL ADJ (A(I,J),A1(I,J),MINT(I,J),SPAC1) INT01650
CALL ADJ (BA(I,J),B1(I,J),NINT(I,J),SPAC1) INT01660
CALL ADJ (CB(I,J),C1(I,J),LINT(I,J),SPAC1) INT01670
110 CONTINUE INT01680
GO TO 130 INT01690
120 PRINT *, date_string
C 120 CALL TIMER(DATE,TIME,IVIRT,ITOT) INT01700
WRITE(NPRT,9100)B(3),date_string INT01710
C WRITE(NPRT,9100)B(3),DATE,TIME INT01710
Z(1) = 0.0 INT01720
Z(2) = 0.0 INT01730
ZE = 0.0 INT01740
CHARGE = 20.0 INT01750
REM = 0.0 INT01760
DEBYE = 0.0 INT01770
BETA = 0.0 INT01780
CONVGA = 0.0 INT01790
130 WRITE(NPRT,9103) Z(1),Z(2),D,TEMP,DT,DP,A1(1,1),A1(2,2),A1(1,2), INT01800
*B1(1,1),B1(2,2),B1(1,2),C1(1,1),C1(2,2),C1(1,2),E(1,1),E(2,2), INT01810
*E(1,2),W(1,1),W(2,2),W(1,2),COEFF1,COEFF2,SPAC1,N,REM,CONC,DEBYE INT01820
9103 FORMAT('Z1=',F5.1,3X,'Z2=',F5.1,/,'DIEL CONST=',F7.3,3X,'TEMP=', INT01830
*F7.3,'DLND/DLNT=',F10.5,3X,'DLND/DP=',F7.3,/,'A11=',F6.3,3X,'A22='INT01840
*,F6.3,3X,'A12=',F6.3,/,'B11=',F6.3,3X,'B22=',F6.3,3X,'B12=',F6.3,/INT01850
*'C11=',F6.3,3X,'C22=',F6.3,3X,'C12=',F6.3,/, INT01860
*'E11=',F6.3,3X,'E22=',F6.3,3X,'E12=',F6.3,/, INT01870
*'W11=',F6.3,3X,'W22=',F6.3,3X,'W12=',F6.3,/, INT01880
*'COEFF1=',F7.4,3X,'COEFF2=',F7.4,3X,'SPAC1=',F7.4,3X,'N=',I4,3X, INT01890
*'REM=',F7.4,/,'CONC=',F10.6,3X,'DEBYE=',F10.5,/) INT01900
C *** CALC OF Q,Q1,AND TQ INT01910
C WRITE(6,5558) INT01920
5558 FORMAT(1X,'HO') INT01930
DO 180 I = 1,NIONS INT01940
DO 180 J = 1,NIONS INT01950
AB = A1(I,J)+B1(I,J) INT01960
BC = AB+C1(I,J) INT01970
DO 180 K = 1,N INT01980
Q(K,I,J) = -Z(I)*Z(J)*FL*EXP(-FK*XR(K))/(PI4*XR(K)) INT01990
IF (XR(K) .LT. A1(I,J)) GO TO 170 INT02000
IF(ABS(XR(K)-AB).LT.1.E-04) GO TO 160 INT02010
IF ((XR(K)-AB).LT.1.E-07) GO TO 150 INT02020
IF(ABS(XR(K)-BC).LT.1.E-04) GO TO 140 INT02030
IF((XR(K)-BC).LT.1.E-07) GO TO 160 INT02040
140 Q1(K,I,J) = EXP(Q(K,I,J)) INT02050
GO TO 180 INT02060
150 Q1(K,I,J) = EXP(Q(K,I,J)-E(I,J)) INT02070
GO TO 180 INT02080
160 Q1(K,I,J) = EXP(Q(K,I,J)-W(I,J)) INT02090
GO TO 180 INT02100
170 Q1(K,I,J) = 0.0 INT02110
180 CONTINUE INT02120
C WRITE(6,5557) INT02130
5557 FORMAT(1X,'HA') INT02140
DO 190 I = 1,NIONS INT02150
DO 190 J = 1,NIONS INT02160
190 CALL FFT3S(F0,Q(1,I,J),M1,RL,TQ(1,I,J),1) INT02170
C***************** Distributution ************************************
C CALC OF APPROX DISTRIBUTION FUNCTIONS after counter I1 =0 *INT02190
1010 I1 = 0 INT02210
IF (ICC.NE.0) GO TO 1070 INT02220
DO 1020 I = 1,NIONS INT02230
DO 1020 J = 1,NIONS INT02240
DO 1020 K = 1,N INT02250
1020 PHI(K,I,J) = Q1(K,I,J)-1. INT02260
DO 1030 I = 1,NIONS INT02270
DO 1030 J = 1,NIONS INT02280
PHI(MINT(I,J),I,J) = (PHI(MINT(I,J),I,J)-1.)*.5 INT02290
1030 CALL FFT3S(F0,PHI(1,I,J),M1,RL,TPHI(1,I,J),1) INT02300
DO 1040 I = 1,NIONS INT02310
DO 1040 J = 1,NIONS INT02320
DO 1040 K = 1,N INT02330
TCONV(K,I,J) = 0.0 INT02340
DO 1040 L = 1,NIONS INT02350
1040 TCONV(K,I,J) = TPHI(K,I,L)*TPHI(K,L,J)*C(L)+TCONV(K,I,J) INT02360
DO 1050 I = 1,NIONS INT02370
DO 1050 J = 1,NIONS INT02380
1050 CALL FFT3S(F0,TCONV(1,I,J),M1,RL,CONV(1,I,J),-1) INT02390
DO 1060 I = 1,NIONS INT02400
DO 1060 J = 1,NIONS INT02410
DO 1060 K = 1,N INT02420
1060 TAU(K,I,J) = CONV(K,I,J)+Q(K,I,J)*XR(K)*FK*.5 INT02430
1070 DO 1080 I = 1,NIONS INT02440
DO 1080 J = 1,NIONS INT02450
DO 1075 K = 1,N INT02460
1075 H(K,I,J) = Q1(K,I,J)*(1.+TAU(K,I,J))-1. INT02470
IF (MINT(I,J) < 1) THEN
MINT(I,J) = 1
END IF
H(MINT(I,J),I,J) = (H(MINT(I,J),I,J)-1.)*.5 INT02480
DO 1080 K = 1,N INT02490
1080 ZETA(K,I,J) = H(K,I,J)-Q(K,I,J)-TAU(K,I,J) INT02500
IF (ICC.NE.0) GO TO 1090 INT02510
WRITE(NPRT,9104) INT02520
9104 FORMAT(38H0RESULTS OF TRIANGLE TERM IN DIST FUNC) INT02530
GO TO 1100 INT02540
1090 IF(ICC.GT.0) GO TO 1112 INT02550
WRITE(NPRT,9105) INT02560
9105 FORMAT(1H0,'RESULTS USING PRECEDING TAU AT PRESENT CONC') INT02570
GO TO 1100 INT02580
1112 WRITE(NPRT,9102) INT02590
9102 FORMAT(1H0,'RESULTS USING PREVIOUS TAU STORED ON DISK') INT02600
1100 WRITE(NPRT,9106) INT02610
9106 FORMAT(6X,4HDELA,9X,4HDELB,9X,4HDIFF,10X,3HOSM,7X,7HXENERGY,5X, INT02620
*7HXVOLUME,4X,4HITER,6X,4HBETA,7X,6HCONVGA,7X,6HCONVGB) INT02630
GO TO 3010 INT02640
C***************** HNC ***********************************************
C SOLVE HNC INTEGRAL EQUATION BY ITERATION *INT02660
2010 WRITE(NPRT,9107) CONTAC,SQ1,SQ2,REMAIN,FG INT02680
WRITE(NPRT,5555)ZE INT02690
IF(ZE.EQ.0.0) GO TO 2020 INT02700
5555 FORMAT(1X,F10.5) INT02710
WRITE(NPRT,9108) INT02720
WRITE(NPRT,9109) (Z(I),F2(I),R2(I),Y(I),F4(I),R4(I),I=1,NIONS) INT02730
WRITE(NPRT,9110) S2,CS2,E2,S4,CS4,E4 INT02740
2020 WRITE(NPRT,9111) INT02750
9111 FORMAT(19H0RESULTS OF HNC EQN) INT02760
WRITE(NPRT,9106) INT02770
2030 CONTINUE INT02780
DO 2040 I = 1,NIONS INT02790
DO 2040 J = 1,NIONS INT02800
2040 CALL FFT3S(F0,H(1,I,J),M1,RL,TH(1,I,J),1) INT02810
DO 2050 I = 1,NIONS INT02820
DO 2050 J = 1,NIONS INT02830
2050 CALL FFT3S(F0,ZETA(1,I,J),M1,RL,TZETA(1,I,J),1) INT02840
DO 2070 I = 1,NIONS INT02850
DO 2070 J = 1,NIONS INT02860
DO 2070 K = 1,N INT02870
SUM1 = 0.0 INT02880
DO 2060 L = 1,NIONS INT02890
2060 SUM1 = TZETA(K,I,L)*TH(K,L,J)*C(L)+SUM1 INT02900
2070 P1(K,I,J) = SUM1 INT02910
DO 2090 I = 1,NIONS INT02920
DO 2090 J = 1,NIONS INT02930
DO 2090 K = 1,N INT02940
P2 = 0.0 INT02950
P3 = 0.0 INT02960
DO 2080 L = 1,NIONS INT02970
P2 = TQ(K,I,L)*TZETA(K,L,J)*C(L)+P2 INT02980
2080 P3 = TQ(K,I,L)*P1(K,L,J)*C(L)+P3 INT02990
2090 TTAU(K,I,J) = P1(K,I,J)+P2+P3 INT03000
DO 2100 I = 1,NIONS INT03010
DO 2100 J = 1,NIONS INT03020
2100 CALL FFT3S(F0,TTAU(1,I,J),M1,RL,TAU(1,I,J),-1) INT03030
I1 = I1+1 INT03040
2110 DO 2120 I = 1,NIONS INT03050
DO 2120 J = 1,NIONS INT03060
C WRITE(7,5556) INT03070
5556 FORMAT(1X,'HI') INT03080
DO 2115 K = 1,N INT03090
C somethign about this raises the floating point exception
C looking deeper, the math library doesnt like that something here has
C an inf value. if is a valid float though. line 39 of
C math_errf.c library.
C it probably would preffer IMPLICIT NONE
C
2115 H(K,I,J) = Q1(K,I,J)*EXP(TAU(K,I,J))-1. INT03100
H(MINT(I,J),I,J) = (H(MINT(I,J),I,J)-1.)*.5 INT03110
DO 2120 K = 1,N INT03120
2120 ZETA(K,I,J) = H(K,I,J)-Q(K,I,J)-TAU(K,I,J) INT03130
IF (I1.EQ.MAX) GO TO 3010 INT03140
IF ((I1/I2)*I2.EQ.I1) GO TO 3010 INT03150
GO TO 2030
C**************** XEnergy*********************************************
C DELA,DELB,OSM,XENERGY,DFX/DP *INT03180
3010 CONTAC = 0.0 INT03200
SQ1 = 0.0 INT03210
SQ2 = 0.0 INT03220
REMAIN = 0.0 INT03230
XE = 0.0 INT03240
DO 3020 I = 1,NIONS INT03250
DO 3020 J = 1,NIONS INT03260
NA = MINT(I,J) INT03270
HNA = H(NA,I,J) INT03280
H(NA,I,J) = HNA*2.+1. INT03290
AA = A(I,J)**3 INT03300
NB = NA +NINT(I,J) INT03310
BB = (A(I,J)+BA(I,J))**3 INT03320
NC = NB+LINT(I,J) INT03330
AC = (A(I,J)+BA(I,J)+CB(I,J))**3 INT03340
AD = PI4*(A1(I,J)-A(I,J)) INT03350
CALL TRAP(NA,NB,I,J,SUM1,SUM2,SUM3) INT03360
T(I,J) = PI4*SUM1*SPAC1 INT03370
G(I,J) = PI4*SUM2*SPAC1 INT03380
G4(I,J) = PI4*SUM3*SPAC1 INT03390
GM(I,J) = G(I,J) INT03400
TM(I,J) = T(I,J) INT03410
XE = XE + FNU(I)*FNU(J)*E(I,J)*(PI43*(BB-AA)+G(I,J)) INT03420
CALL TRAP(NB,NC,I,J,SUM1,SUM2,SUM3) INT03430
T(I,J) = PI4*SUM1*SPAC1 INT03440
G(I,J) = PI4*SUM2*SPAC1 INT03450
G4(I,J) = PI4*SUM3*SPAC1+G4(I,J) INT03460
GW(I,J) = G(I,J) INT03470
TW(I,J) = T(I,J) INT03480
XE = XE+FNU(I)*FNU(J)*W(I,J)*(PI43*(AC-BB)+G(I,J)) INT03490
ERROR1 = -Z(I)*Z(J)*FY+AD*H(NA,I,J)*XR(NA) INT03500
ERROR2 = -Z(I)*Z(J)*FX+AD*H(NA,I,J)*XR(NA)**2 INT03510
ERROR3 = -Z(I)*Z(J)*FZ+AD*H(NA,I,J)*XR(NA)**4 INT03520
C PRINT *, "*****From xEnergy*******"
C PRINT *, "NC",NC,"N",N,"I",I,"J",J
C PRINT *,"SUM1",SUM1,"SUM2",SUM2,"SUM3",SUM3
CALL SIMP(NC,N,I,J,SUM1,SUM2,SUM3) INT03530
T(I,J) = PI43*SPAC1*SUM1-PI2*A(I,J)**2+TW(I,J)+TM(I,J)+ERROR1 INT03540
G(I,J) = PI43*(SPAC1*SUM2-A(I,J)**3)+GM(I,J)+GW(I,J)+ERROR2 INT03560
G4(I,J) = PI43*SPAC1*SUM3-2.513234*A(I,J)**5+G4(I,J)+ERROR3 INT03570
GA(I,J) = H(NA,I,J)+1. INT03580
GB(I,J) = H(NB,I,J)+1. INT03590
GC(I,J) = H(NC,I,J)+1. INT03600
CONTAC = CONTAC + FNU(I)*FNU(J)*GA(I,J)*AA INT03610
SQ1 = SQ1+FNU(I)*FNU(J)*(GB(I,J)*BB*(1.-EXP(-E(I,J)+W(I,J)))) INT03620
SQ2 = SQ2+FNU(I)*FNU(J)*(GC(I,J)*AC*(1.-EXP(-W(I,J)))) INT03630
REMAIN = REMAIN + Z(I)*Z(J)*FNU(I)*FNU(J)*T(I,J) INT03640
IF(I1.EQ.MAX) GO TO 3011 INT03650
H(NA,I,J) = HNA INT03660
GO TO 3020 INT03670
3011 DC(NA,I,J) = DC(NA,I,J)+HNA+1. INT03680
3020 CONTINUE INT03690
CONTAC = CONTAC*2.09439*CC/FNO**2 INT03700
SQ1 = SQ1*2.09439*CC/FNO**2 INT03710
SQ2 = SQ2*2.09439*CC/FNO**2 INT03720
REMAIN = REMAIN*FL1*CC/(6.*FNO**2) INT03730
OSM = 1.+CONTAC+SQ1+SQ2+REMAIN INT03740
XE = 1.97827*TEMP*(3.*(1.+DT)*REMAIN*FNO+XE*CC*0.5/FNO) INT03750
XFP = -8.31467*1.E-05*TEMP*DP*3.*REMAIN*FNO INT03760
GM12 = (G(1,2)+G(2,1))*.5 INT03770
CONVGB = 100*(G(1,2)-G(2,1))/GM12 INT03780
IF (ZE.EQ.0.0) GO TO 3050 INT03790
S2 = 0.0 INT03800
S4 = 0.0 INT03810
DO 3040 I = 1,NIONS INT03820
SUM1 = 0.0 INT03830
SUM2 = 0.0 INT03840
DO 3030 J = 1,NIONS INT03850
SUM1 = SUM1 + C(J)*Z(J)*G(I,J) INT03860
3030 SUM2 = SUM2 + C(J)*Z(J)*G4(I,J) INT03870
F2(I) = SUM1 INT03880
F4(I) = SUM2 INT03890
Y(I) = 6.*Z(I)/FKFK INT03900
R2(I) = 1.+F2(I)/Z(I) INT03910
R4(I) = 1. + F4(I)/Y(I) INT03920
S2 = C(I)*Z(I)*F2(I)+S2 INT03930
3040 S4 = C(I)*Z(I)*F4(I)+S4 INT03940
CONVGA = (C(1)*G(1,1)-C(2)*G(2,2))/GM12 - C(1)+C(2) INT03950
CS2 = FKFK/FL INT03960
CS4 = 6./FL INT03970
E2 = 1.+S2/CS2 INT03980
E4 = 1.+S4/CS4 INT03990
BETA = FNO*R2(1)*R2(2)/(FNU(1)*R2(2)+FNU(2)*R2(1)) INT04000
3050 CONTINUE INT04010
FG = 1.0/(CC*GM12) INT04020
DELA = (FG-1.0-FG*BETA)/(1.0+FG*BETA) INT04030
DELB = FG-1.0 INT04040
DIFF = DELA-DELB INT04050
ITER = I1 INT04060
FITER = ITER INT04070
WRITE(NPRT,9112) DELA,DELB,DIFF,OSM,XE,XFP,ITER,BETA,CONVGA,CONVGBINT04080
9112 FORMAT( 6E13.4,I5,3E13.4) INT04090
IF (I1.EQ.MAX) GO TO 3060 INT04100
IF(I1.EQ.0) GO TO 2010 INT04110
IF ((I1/I2)*I2.EQ.I1) GO TO 2030 INT04120
3060 WRITE(NPRT,9113) INT04130
9113 FORMAT(15H0MATRIX G(I,J),21X,15HMATRIX GM(I,J),20X, INT04140
*15HMATRIX GW(I,J),20X,15HMATRIX G4(I,J)) INT04150
WRITE(NPRT,9114) ((G(I,J),J=1,NIONS),(GM(I,L),L=1,NIONS),(GW(I,M),INT04160
*M=1,NIONS),(G4(I,K),K=1,NIONS),I=1,NIONS) INT04170
9114 FORMAT(1X,2E14.6,6X,2E14.6,6X,2E14.6,6X,2E14.6) INT04180
WRITE(NPRT,9115) INT04190
9115 FORMAT(15H0MATRIX T(I,J),21X,15HMATRIX TM(I,J),20X, INT04200
*15HMATRIX TW(I,J)) INT04210
WRITE(NPRT,9116)((T(I,J),J=1,NIONS),(TM(I,L),L=1,NIONS),(TW(I,M),MINT04220
*=1,NIONS),I=1,NIONS) INT04230
9116 FORMAT(1X,2E14.6,6X,2E14.6,6X,2E14.6) INT04240
WRITE(NPRT,9117) INT04250
9117 FORMAT(16H0MATRIX GA(I,J),20X,15HMATRIX GB(I,J),20X, INT04260
*15HMATRIX GC(I,J)) INT04270
WRITE(NPRT,9118) ((GA(I,J),J=1,NIONS),(GB(I,L),L=1,NIONS),(GC(I,M)INT04280
*,M=1,NIONS),I=1,NIONS) INT04290
9118 FORMAT(1X,2E14.6,6X,2E14.6,6X,2E14.6) INT04300
WRITE(NPRT,9119) INT04310
9119 FORMAT(17H0MATRIX MINT(I,J),3X,16HMATRIX NINT(I,J),3X, INT04320
*16HMATRIX LINT(I,J)) INT04330
WRITE(NPRT,9120) ((MINT(I,J),J=1,NIONS),(NINT(I,L),L=1,NIONS),(LININT04340
*T(I,M),M=1,NIONS),I=1,NIONS) INT04350
9120 FORMAT(1X,2I4,13X,2I4,13X,2I4) INT04360
WRITE(NPRT,9107) CONTAC,SQ1,SQ2,REMAIN,FG INT04370
9107 FORMAT(8H0CONTAC=,1P,E12.4,3X,'SQUARE1=',1P,E12.4,3X,'SQUARE2=', INT04380
*1P,E12.4,3X,'REMAIN=',1P,E12.4,3X,'GAMMA=',1P,E12.4) INT04390
IF (ZE.EQ.0.0) GO TO 3070 INT04400
WRITE(NPRT,9108) INT04410
9108 FORMAT(6X,1HZ,13X,2HF2,12X,2HR2,12X,1HY,13X,2HF4,12X,2HR4) INT04420
WRITE(NPRT,9109) (Z(I),F2(I),R2(I),Y(I),F4(I),R4(I),I=1,NIONS) INT04430
9109 FORMAT(1X,6E14.6) INT04440
WRITE(NPRT,9110) S2,CS2,E2,S4,CS4,E4 INT04450
9110 FORMAT('S2=',E14.6,2X,'CS2=',1P,E14.6,2X,'E2=',1P,E14.6,2X,'S4=' INT04460
*,1P,E14.6,2X,'CS4=',1P,E14.6,2X,'E4=',1P,E14.6) INT04470
3070 CONTINUE INT04480
WRITE(NPRT,9121) INT04490
9121 FORMAT(7X,5HXR(N),7X,9H Q(1,1,N),5X,8HH(1,1,N),5X,9H Q(1,2,N), INT04500
*6X,8HH(1,2,N),5X,9H Q(2,2,N),6X,8HH(2,2,N)) INT04510
WRITE(NPRT,9122) XR(N),Q(N,1,1),H(N,1,1),Q(N,1,2),H(N,1,2), INT04520
*Q(N,2,2),H(N,2,2) INT04530
9122 FORMAT(7E14.4//) INT04540
IF(EX.NE.2.0.AND.EX.NE.3.0) GO TO 4070 INT04550
3080 WRITE(NPRT,9123) INT04560
9123 FORMAT(7X,5HXR(K),6X,10H DC(1,1,K),5X,8HH(1,1,K),5X,10H DC(1,2,K),INT04570
*5X,8HH(1,2,K),5X,10H DC(2,1,K),5X,8HH(2,1,K),5X,10H DC(2,2,K), INT04580
*5X,8HH(2,2,K)) INT04590
WRITE(NPRT,9124)(XR(K),((DC(K,I,J),H(K,I,J),J=1,NIONS),I=1,NIONS) INT04600
*,K=1,N) INT04610
WRITE(NDSK,9125)C1,C2 INT04620
WRITE(NDSK,9133)FO,N,M1,RL,((MINT(I,J),I=1,2),J=1,2) INT04630
9133 FORMAT(F8.4,2I4,F8.4,4I15.6) INT04640
WRITE(NDSK,9126)(XR(K),H(K,1,1),H(K,1,2),K=1,N) INT04650
WRITE(NDSK2,9003)(((TAU(K,I,J),I=1,NIONS),J=1,NIONS),K=1,N) INT04660
REWIND NDSK2 INT04670
9124 FORMAT(1X,1P,9E14.4) INT04680
C ERIC some issue with this line format
9125 FORMAT(1X,A8,2E16.6) INT04690
9126 FORMAT(2E16.7) INT04700
IF (ZE.EQ.0.0) GO TO 4030 INT04710
IF (EX.NE.3.0) GO TO 4070 INT04720
C **********************************************************************INT04730
C CALC OF CHARGE DENSITY *INT04740
C **********************************************************************INT04750
DO 4020 K = 1,N INT04760
DO 4020 I = 1,NIONS INT04770
SUM1 = 0.0 INT04780
SUM2 = 0.0 INT04790
DO 4010 J = 1,NIONS INT04800
SUM1 = SUM1 + FNU(J)*Z(J)*H(K,I,J) INT04810
4010 SUM2 = SUM2 + FNU(J)*Z(J)*Q(K,I,J) INT04820
RHO(I,K) = SUM1*XR(K)**2 INT04830
4020 RQO(I,K) = SUM2*XR(K)**2 INT04840
GO TO 4060 INT04850
4030 Z(1) = 1.0 INT04860
Z(2) = -1.0 INT04870
DO 4050 K = 1,N INT04880
DO 4050 I = 1,NIONS INT04890
SUM1 = 0.0 INT04900
DO 4040 J = 1,NIONS INT04910
4040 SUM1 = SUM1+FNU(J)*Z(J)*H(K,I,J) INT04920
RHO(I,K) = SUM1*XR(K)**2 INT04930
4050 RQO(I,K) = 0.0 INT04940
4060 WRITE(NPRT,9127) INT04950
9127 FORMAT('1CHARGE DENSITY'///1X,9X,'XR',11X,'RHO1',10X,'RQO1',11X, INT04960
*'RHO2',10X,'RQO2') INT04970
C PRINT *, 'NCD:' , NCD
C PRINT *, 'K:' , K
C PRINT *, 'IC:' , IC
WRITE(NPRT,9128)(XR(K),RHO(1,K),RQO(1,K),RHO(2,K),RQO(2,K),K=1,NCDINT04980
*,IC) INT04990
PRINT *, 'NR:' , NR
9128 FORMAT(1X,1P,5E14.4) INT05000
WRITE(NDSK,9126)(RHO(1,K),RHO(2,K),K=1,N2) INT05010
9129 FORMAT(2E16.7) INT05020
REWIND NDSK INT05030
4070 CONTINUE INT05040
B(1) = 0.0 INT05050
B(2) = 1.0 INT05060
Z(1) = ZA(1) INT05070
Z(2) = ZA(2) INT05080
ZE = Z(1) INT05090
IF (ZE.EQ.0.0) GO TO 5040 INT05100
C **********************************************************************INT05110
C CALC OF DIELECTRIC RESPONSE FUNCTION INT05120
C **********************************************************************INT05130
DO 5010 I = 1,NIONS INT05140
DO 5010 J = 1,NIONS INT05150
H(MINT(I,J),I,J) = (H(MINT(I,J),I,J)-1.)*.5 INT05160
5010 CALL FFT3S(F0,H(1,I,J),M1,RL,TH(1,I,J),1) INT05170
SPACT = PI2/(2.*FNN-1.)/SPAC1 INT05180
TK = 1.0 INT05190
DO 5030 K = 1,N INT05200
Y1(K) = TK*SPACT INT05210
SUM = 0.0 INT05220
DO 5020 I = 1,NIONS INT05230
DO 5020 J = 1,NIONS INT05240
5020 SUM = SUM + C(I)*C(J)*Z(I)*Z(J)*TH(K,I,J) INT05250
Y1K2 = Y1(K)**2 INT05260
DRF(K) = 1.-(FKFK/Y1K2)*(1.+SUM/FC) INT05270
DRDF(K) = 1.-FKFK/(FKFK+Y1K2) INT05280
5030 TK = TK+1. INT05290
WRITE(NPRT,9130) INT05300
9130 FORMAT('1DIELECTRIC RESPONSE FUNCTION'///1X,3(9X,'YK',12X,'DRF', INT05310
*11X,'DRDF')) INT05320
WRITE(NPRT,9131)(Y1(K),DRF(K),DRDF(K),K=1,NR,IR) INT05330
9131 FORMAT(1X,1P,9E14.4) INT05340
5040 CONTINUE
C5040 CALL TIMER(DATE,TIME,IVIRT,ITOT) INT05350
C RVIRT=FLOAT(IVIRT)*1.E-6 INT05360
C WRITE(NPRT,9132)TIME,DATE,RVIRT INT05370
WRITE(NPRT,9132)date_string INT05370
C9132 FORMAT(//////////'END OF EXECUTION AT',A8,' ON ',A8,5X,F8.3, INT05380
C *' SECONDS USED') INT05390
9132 FORMAT(//////////'END EXECUTION ON ',A20) INT05390
IF (ABS(CHARGE-0.0).LT.1.E-04) GO TO 120 INT05400
C GOTO 40 INT05410
C Program loop
5050 CONTINUE INT05420
STOP INT05430
END INT05440
C ****Program stops and ends ****************************************
C************************ SUBROUTINES *******************************
C************************ ADJ ***************************************
SUBROUTINE ADJ(A,A1,M,SPAC1) INT05450
S1 = A/SPAC1 INT05460
L = S1+1 INT05470
S2 = L INT05480
IF(S2-S1.GE.0.5) GO TO 10 INT05490
A1 = S2*SPAC1 INT05500
M = L INT05510
GO TO 20 INT05520
10 A1 = (S2-1.)*SPAC1 INT05530
M = L-1 INT05540
20 CONTINUE INT05550
RETURN INT05560
END
C************************ TRAP ***************************************
SUBROUTINE TRAP(N1,N11,I,J,SUM1,SUM2,SUM3) INT05580
DIMENSION H(2048,2,2),XR(2048) INT05590
COMMON/JAY/ H,XR INT05600
IF (N11.EQ.N1) GO TO 20 INT05610
SUM1 = (H(N1,I,J)*XR(N1)+H(N11,I,J)*XR(N11))*.5 INT05620
SUM2 = (H(N1,I,J)*XR(N1)**2+H(N11,I,J)*XR(N11)**2)*.5 INT05630
SUM3 = (H(N1,I,J)*XR(N1)**4+H(N11,I,J)*XR(N11)**4)*.5 INT05640
L1 = N1+1 INT05650
L2 = N11-1 INT05660
IF (N11.EQ.L1) GO TO 30 INT05670
DO 10 K = L1,L2 INT05680
SUM1 = SUM1 + H(K,I,J)*XR(K) INT05690
SUM3 = SUM3 + H(K,I,J)*XR(K)**4 INT05700
10 SUM2 = SUM2 + H(K,I,J)*XR(K)**2 INT05710
GO TO 30 INT05720
20 SUM1 = 0.0 INT05730
SUM2 = 0.0 INT05740
SUM3 = 0.0 INT05750
30 RETURN INT05760
END
C********************** SIMPSON **************************************
SUBROUTINE SIMP(NC,N,I,J,SUM1,SUM2,SUM3) INT05780
DIMENSION H(2048,2,2),XR(2048) INT05790
COMMON /JAY/ H,XR INT05800
SUM1 = H(NC,I,J)*XR(NC) INT05810
SUM2 = SUM1*XR(NC) INT05820
SUM3 = SUM2*XR(NC)**2 INT05830
C PRINT *, 'NC:',NC
C Setting NC to 1 for testing??
NC = 1
N2 = (N-NC-2)/2
C PRINT *, '**N2:', N2
DO 10 K = 1,N2 INT05850
K4 = NC+2*K-1 INT05860
K5 = K4+1 INT05870
K6 = NC+2*N2+1 INT05880
K7 = K6+1 INT05890
SUM1 = SUM1 + 4.*H(K4,I,J)*XR(K4)+2.*H(K5,I,J)*XR(K5) INT05900
SUM3 = SUM3 + 4.*H(K4,I,J)*XR(K4)**4+2.*H(K5,I,J)*XR(K5)**4 INT05910
10 SUM2 = SUM2 + 4.*H(K4,I,J)*XR(K4)**2+2.*H(K5,I,J)*XR(K5)**2 INT05920
SUM1 = SUM1 + 4.*H(K6,I,J)*XR(K6)+2.*H(K7,I,J)*XR(K7) INT05930
SUM2 = SUM2 + 4.*H(K6,I,J)*XR(K6)**2+2.*H(K7,I,J)*XR(K7)**2 INT05940
SUM3 = SUM3 + 4.*H(K6,I,J)*XR(K6)**4+2.*H(K7,I,J)*XR(K7)**4 INT05950
RETURN INT05960
END
C**************** Fast Fourier Transform *************************
SUBROUTINE FFT3S(F0,F,MM,RL,T,MODE) INT05990
C FAST FOURIER TRANSFORM OF 3-DIMENSIONAL CENTROSYMMETRIC FUNCTION INT06000
DIMENSION F(2048),T(1) INT06010
COMMON /FDATA/ N,N2,FF,RTTWO,NC,NS,ND,NR INT06020
COMMON /FWORK/ W(12288) INT06030
C DIMENSION OF W SHOULD BE .GE. 3*N INT06040
C FIXED DATA PREPARATION INT06050
DATA PI/3.141593/,M/-1/ INT06060
MD = MM+1 INT06070
CALL SETFT(MD,M) INT06080
N = 2**MM INT06090
RN = FLOAT(N) INT06100
DR = RL/RN INT06110
RL3 = RL*RL*RL INT06120
C = 4.*RL3*FF/RN INT06130
CC = 0.0 INT06140
A = 1.0 INT06150
IF(IABS(MODE)-2)202,200,201 INT06160
200 CC = 0.5*PI*DR*DR*DR*F(1) INT06170
A = .125*RN*RN*RN/(RL3*RL3) INT06180
GOTO 202 INT06190
201 CC = F0 INT06200
202 IF(MODE.GT.0)GOTO 203 INT06210
C = 1.0/C INT06220
CC = A*CC INT06230
203 CONTINUE INT06240
C COPY INPUT VECTOR INTO WORKING STORAGE INT06250
DO 204 I=1,N INT06260
204 W(ND+I+1) = C*FF*F(I)*I INT06270
C Eric Starts debugging here
C RETURN
C END
CALL RSA(M) INT06280
C COPY INTO OUTPUT VECTOR
C N eq 2048
C DO 205 I=1,N
DO 205 I=1,1
205 T(I) = W(NR+I+1)/FLOAT(I)+CC INT06310
RETURN INT06320
END INT06330
C**************************** SETFT **********************************
SUBROUTINE SETFT(MM,M) INT06340
C PREPARE FIXED DATA FOR FAST FOURIER TRANSFORM INT06350
COMMON /FDATA/ N,N2,F,RTTWO,NC,NS,NDATA,NRES INT06360
COMMON /FWORK/ W(12288) INT06370
IF(MM.LT.0)RETURN INT06380
IF(MM.EQ.M)RETURN INT06390
M = MM INT06400
N = 2**M INT06410
N2 = N/2 INT06420
F = 1.0/FLOAT(N) INT06430
DA = 6.2831853071796D0*F INT06440
RTTWO = SQRT(2.0) INT06450
F = SQRT(F) INT06460
N1 = (N2-1)/2 INT06470
NC = 0 INT06480
NS = N1 INT06490
NDATA = N1+N1 INT06500
IF(N1.LE.0)RETURN INT06510
DO 10 I=1,N1 INT06520
A = FLOAT(I)*DA INT06530
W(I) = COS(A) INT06540
10 W(N1+I) = SIN(A) INT06550
RETURN INT06560
END
C************************** RSA *************************************
SUBROUTINE RSA(M) INT06580
C REAL ODD ANALYSIS OR SYNTHESIS INT06590
COMMON /FDATA/ N,N2,F,RTTWO,NC,NS,NDATA,NRES INT06600
COMMON /FWORK/ W(12288) INT06610
NQ = N2 INT06620
NIB = NDATA +1 INT06630
NOB = NDATA+N2+2 INT06640
DO 200 L = 1,M INT06650
NIM = NIB+NQ INT06660
NIE = NIM+NQ INT06670
NQH = NQ/2 INT06680
NQL = (NQ-1)/2 INT06690
NOM = NOB+NQL+1 INT06700
C TPRIME = 0 INT06710
IF(NQL.EQ.0)GOTO 20 INT06720
DO 10 IT = 1,NQL INT06730
W(NOB+IT)=W(NIB+IT)-W(NIM-IT) INT06740
10 W(NOM+IT) = W(NIB+IT)+W(NIM-IT) INT06750
20 IF(NQH.NE.NQL) W(NOM) = 2.0*W(NIB+NQH)
C ?? 0.LT.TPRIME.LT.NP/2. INT06770
NO1 = NQ INT06780
NO2 = N2-NQ INT06790
IF(NO1-NO2)40,80,120 INT06800
40 NI=NO1+NO1 INT06810
W(NOB+NO1) = W(NIB+NI)+W(NIM+NI) INT06820
W(NOB+NO2) = -W(NIB+NI)+W(NIM+NI) INT06830
IF(NQL.EQ.0)GOTO 60 INT06840
CC = W(NC+NO1) INT06850
SS = W(NS+NO1) INT06860
DO 50 IT = 1,NQL INT06870
RE = CC*W(NIE+NI-IT)-SS*W(NIM+NI-IT) INT06880
AI = SS*W(NIE+NI-IT)+CC*W(NIM+NI-IT) INT06890
W(NOM+NO1+IT) = W(NIM+NI+IT)-RE INT06900
W(NOM+NO2+IT) = W(NIM+NI+IT)+RE INT06910
W(NOB+NO1+IT) = W(NIB+NI+IT)+AI INT06920
50 W(NOB+NO2+IT) = -W(NIB+NI+IT)+AI INT06930
60 IF(NQH.EQ.NQL)GOTO 70 INT06940
NR = NO1/2 INT06950
W(NOM+NO1) = 2.0*(W(NS+NR)*W(NIM+NI+NQH)+W(NC+NR)*W(NIB+NI+NQH)) INT06960
W(NOM+NO2) = 2.0*(W(NC+NR)*W(NIM+NI+NQH)-W(NS+NR)*W(NIB+NI+NQH)) INT06970
70 NO1 = NO1+NQ INT06980
NO2 = NO2-NQ INT06990
IF(NO1-NO2)40,80,120 INT07000
80 W(NOB+NO1)=W(NIM) INT07010
IF(NQL.EQ.0)GOTO 100 INT07020
DO 90 IT = 1,NQL INT07030
W(NOM+NO1+IT) = W(NIM+IT) INT07040
90 W(NOB+NO1+IT) = W(NIE-IT) INT07050
100 IF(NQH.NE.NQL)W(NOM+NO1) = RTTWO*W(NIM+NQH) INT07060
120 NT = NIB INT07070
NIB = NOB INT07080
NOB = NT INT07090
NQ = NQH INT07100
200 CONTINUE INT07110
NRES = NIB-1 INT07120
RETURN INT07130
END INT07140
