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